Coverage for /usr/local/lib/python3.8/site-packages/spam/filters/movingFilters.py: 100%

1""" 

2Library of SPAM filters. 

3Copyright (C) 2020 SPAM Contributors 

4 

5This program is free software: you can redistribute it and/or modify it 

6under the terms of the GNU General Public License as published by the Free 

7Software Foundation, either version 3 of the License, or (at your option) 

8any later version. 

9 

10This program is distributed in the hope that it will be useful, but WITHOUT 

11ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

13more details. 

14 

15You should have received a copy of the GNU General Public License along with 

16this program. If not, see <http://www.gnu.org/licenses/>. 

17""" 

18 

19import numpy 

20import spam.mesh 

21 

22# Default structural element 

23# 0 0 0 

24# 0 1 0 

25# 0 0 0 

26# 0 1 0 

27# 1 2 1 

28# 0 1 0 

29# 0 0 0 

30# 0 1 0 

31# 0 0 0 

32structEl = spam.mesh.structuringElement(radius=1, order=1).astype('<f4') 

33structEl[1, 1, 1] = 2.0 

34 

35 

36def average(im, structEl=structEl): 

37 """ 

38 This function calculates the average map of a grey scale image over a structuring element 

39 It works for 3D and 2D images 

40 

41 Parameters 

42 ---------- 

43 im : 3D or 2D numpy array 

44 The grey scale image for which the average map will be calculated 

45 

46 structEl : 3D or 2D numpy array, optional 

47 The structural element defining the local window-size of the operation 

48 Note that the value of each component is considered as a weight (ponderation) for the operation 

49 (see `spam.mesh.structured.structuringElement` for details about the structural element) 

50 Default = radius = 1 (3x3x3 array), order = 1 (`diamond` shape) 

51 

52 Returns 

53 ------- 

54 imFiltered : 3D or 2D numpy array 

55 The averaged image 

56 """ 

57 import spam.filters.filtersToolkit as mft 

58 

59 # Detect 2D image: 

60 if len(im.shape) == 2: 

61 # pad it 

62 im = im[numpy.newaxis, ...] 

63 structEl = structEl[numpy.newaxis, ...] 

64 

65 imFiltered = numpy.zeros_like(im).astype('<f4') 

66 mft.average(im, imFiltered, structEl) 

67 

68 # Return back 2D image: 

69 if im.shape[0] == 1: 

70 imFiltered = imFiltered[0, :, :] 

71 

72 return imFiltered 

73 

74 

75def variance(im, structEl=structEl): 

76 """" 

77 This function calculates the variance map of a grey scale image over a structuring element 

78 It works for 3D and 2D images 

79 

80 Parameters 

81 ---------- 

82 im : 3D or 2D numpy array 

83 The grey scale image for which the variance map will be calculated 

84 

85 structEl : 3D or 2D numpy array, optional 

86 The structural element defining the local window-size of the operation 

87 Note that the value of each component is considered as a weight (ponderation) for the operation 

88 (see `spam.mesh.structured.structuringElement` for details about the structural element) 

89 Default = radius = 1 (3x3x3 array), order = 1 (`diamond` shape) 

90 

91 Returns 

92 ------- 

93 imFiltered : 3D or 2D numpy array 

94 The variance image 

95 """ 

96 import spam.filters.filtersToolkit as mft 

97 

98 # Detect 2D image: 

99 if len(im.shape) == 2: 

100 # pad it 

101 im = im[numpy.newaxis, ...] 

102 structEl = structEl[numpy.newaxis, ...] 

103 

104 imFiltered = numpy.zeros_like(im).astype('<f4') 

105 mft.variance(im, imFiltered, structEl) 

106 

107 # Return back 2D image: 

108 if im.shape[0] == 1: 

109 imFiltered = imFiltered[0, :, :] 

110 

111 return imFiltered 

112 

113 

114def hessian(im): 

115 """ 

116 This function computes the hessian matrix of grey values (matrix of second derivatives) 

117 and returns eigenvalues and eigenvectors of the hessian matrix for each voxel... 

118 I hope you have a lot of memory! 

119 

120 Parameters 

121 ----------- 

122 im: 3D numpy array 

123 The grey scale image for which the hessian will be calculated 

124 

125 Returns 

126 -------- 

127 list containing two lists: 

128 List 1: contains 3 different 3D arrays (same size as im): 

129 Maximum, Intermediate, Minimum eigenvalues 

130 List 2: contains 9 different 3D arrays (same size as im): 

131 Components Z, Y, X of Maximum 

132 Components Z, Y, X of Intermediate 

133 Components Z, Y, X of Minimum eigenvalues 

134 """ 

135 # 2018-10-24 EA OS MCB "double" hessian fracture filter 

136 # There is already an imageJ implementation, but it does not output eigenvectors 

137 import spam.filters.filtersToolkit as ftk 

138 

139 im = im.astype('<f4') 

140 

141 # The hessian matrix in 3D is a 3x3 matrix of gradients embedded in each point... 

142 # hessian = numpy.zeros((3,3,im.shape[0],im.shape[1],im.shape[2]), dtype='<f4' ) 

143 hzz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

144 hzy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

145 hzx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

146 hyz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

147 hyy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

148 hyx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

149 hxz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

150 hxy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

151 hxx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

152 eigValA = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

153 eigValB = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

154 eigValC = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

155 eigVecAz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

156 eigVecAy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

157 eigVecAx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

158 eigVecBz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

159 eigVecBy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

160 eigVecBx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

161 eigVecCz = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

162 eigVecCy = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

163 eigVecCx = numpy.zeros((im.shape[0], im.shape[1], im.shape[2]), dtype='<f4') 

164 

165 # gaussian 

166 # gauss = ndi.gaussian_filter(image,sigma=0.5, mode='constant', cval=0) 

167 

168 # first greylevel derivative, return Z, Y, X gradients 

169 gradient = numpy.gradient(im) 

170 

171 # second derivate 

172 tmp = numpy.gradient(gradient[0]) 

173 hzz = tmp[0] 

174 hzy = tmp[1] 

175 hzx = tmp[2] 

176 tmp = numpy.gradient(gradient[1]) 

177 hyz = tmp[0] 

178 hyy = tmp[1] 

179 hyx = tmp[2] 

180 tmp = numpy.gradient(gradient[2]) 

181 hxz = tmp[0] 

182 hxy = tmp[1] 

183 hxx = tmp[2] 

184 

185 del tmp, gradient 

186 

187 # run eigen solver for each pixel!!! 

188 ftk.hessian(hzz, hzy, hzx, hyz, hyy, hyx, hxz, hxy, hxx, 

189 eigValA, eigValB, eigValC, 

190 eigVecAz, eigVecAy, eigVecAx, 

191 eigVecBz, eigVecBy, eigVecBx, 

192 eigVecCz, eigVecCy, eigVecCx) 

193 

194 return [[eigValA, eigValB, eigValC], [eigVecAz, eigVecAy, eigVecAx, eigVecBz, eigVecBy, eigVecBx, eigVecCz, eigVecCy, eigVecCx]] 

195 

196 

197# old equivalent 100% python stuff... way slower 

198# def _moving_average( im, struct=default_struct ): 

199# """ 

200# Calculate the average of a grayscale image over a 3x3x3 structuring element 

201# The output is float32 since results is sometimes outof the uint bounds during computation 

202# PARAMETERS: 

203# - im (numpy.array): The grayscale image to be treated 

204# - struct (array of int): the structural element. 

205# Note that the value of each component is considerred as a weight (ponderation) for the operation 

206# RETURNS: 

207# - o_im (numpy.array float32): The averaged image 

208# HISTORY: 

209# 2016-03-24 (JBC): First version of the function 

210# 2016-04-05 (ER): generalisation using structural elements 

211# 2016-05-03 (ER): add progress bar 

212# """ 

213# # Step 1: init output_image as float 32 bits 

214# o_im = numpy.zeros( im.shape, dtype='<f4' ) 

215# # Step 2: structutral element 

216# structural_element_size = int( len( struct )/2 ) 

217# structural_element_weight = float( struct.sum() ) 

218# if prlv>5: 

219# import progressbar 

220# max_values = len( numpy.argwhere( struct ) ) 

221# bar = progressbar.ProgressBar( maxval=max_values, widgets=['Average filter: ', progressbar.Percentage(), progressbar.Bar('=', '[', ']')] ) 

222# bar.start() 

223# for i, c in enumerate( numpy.argwhere( struct ) ): 

224# # convert structural element coordinates into shift to apply 

225# shift_to_apply = c-structural_element_size 

226# # get the local weight from the structural element value 

227# current_voxel_weight = float( struct[c[0], c[1], c[2]] ) 

228# # if prlv>5: print ' Shift {} of weight {}'.format( shift_to_apply, current_voxel_weight ) 

229# # output_image = output_image + ( voxel_weight * image / element_weight ) 

230# o_im += current_voxel_weight*sman.multiple_shifts( im, shifts=shift_to_apply )/structural_element_weight 

231# if prlv>5: bar.update( i+1 ) 

232# if prlv>5: bar.finish() 

233# return o_im 

234# 

235# def _moving_variance( im, struct=default_struct ): 

236# """ 

237# Calculate the variance of a grayscale image over a 3x3x3 structuring element 

238# The output is float32 since results is sometimes outof the uint bounds during computation 

239# PARAMETERS: 

240# - image (numpy.array): The grayscale image to be treat 

241# RETURNS: 

242# - o_im (numpy.array): The varianced image 

243# HISTORY: 

244# 2016-04-05 (ER): First version of the function 

245# """ 

246# # Step 1: return E(im**2) - E(im)**2 

247# return moving_average( numpy.square( im.astype('<f4') ), struct=struct ) - numpy.square( moving_average( im, struct=struct ) )