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

1# Library of SPAM filters. 

2# Copyright (C) 2020 SPAM Contributors 

3# 

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

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

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

7# any later version. 

8# 

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

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

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

12# more details. 

13# 

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

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

16 

17import numpy 

18import spam.mesh 

19 

20# Default structural element 

21# 0 0 0 

22# 0 1 0 

23# 0 0 0 

24# 0 1 0 

25# 1 2 1 

26# 0 1 0 

27# 0 0 0 

28# 0 1 0 

29# 0 0 0 

30structEl = spam.mesh.structuringElement(radius=1, order=1).astype("<f4") 

31structEl[1, 1, 1] = 2.0 

32 

33 

34def average(im, structEl=structEl): 

35 """ 

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

37 It works for 3D and 2D images 

38 

39 Parameters 

40 ---------- 

41 im : 3D or 2D numpy array 

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

43 

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

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

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

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

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

49 

50 Returns 

51 ------- 

52 imFiltered : 3D or 2D numpy array 

53 The averaged image 

54 """ 

55 import spambind.filters.filtersToolkit as mft 

56 

57 # Detect 2D image: 

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

59 # pad it 

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

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

62 

63 imFiltered = numpy.zeros_like(im).astype("<f4") 

64 mft.average(im, imFiltered, structEl) 

65 

66 # Return back 2D image: 

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

68 imFiltered = imFiltered[0, :, :] 

69 

70 return imFiltered 

71 

72 

73def variance(im, structEl=structEl): 

74 """ " 

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

76 It works for 3D and 2D images 

77 

78 Parameters 

79 ---------- 

80 im : 3D or 2D numpy array 

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

82 

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

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

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

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

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

88 

89 Returns 

90 ------- 

91 imFiltered : 3D or 2D numpy array 

92 The variance image 

93 """ 

94 import spambind.filters.filtersToolkit as mft 

95 

96 # Detect 2D image: 

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

98 # pad it 

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

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

101 

102 imFiltered = numpy.zeros_like(im).astype("<f4") 

103 mft.variance(im, imFiltered, structEl) 

104 

105 # Return back 2D image: 

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

107 imFiltered = imFiltered[0, :, :] 

108 

109 return imFiltered 

110 

111 

112def hessian(im): 

113 """ 

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

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

116 I hope you have a lot of memory! 

117 

118 Parameters 

119 ----------- 

120 im: 3D numpy array 

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

122 

123 Returns 

124 -------- 

125 list containing two lists: 

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

127 Maximum, Intermediate, Minimum eigenvalues 

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

129 Components Z, Y, X of Maximum 

130 Components Z, Y, X of Intermediate 

131 Components Z, Y, X of Minimum eigenvalues 

132 """ 

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

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

135 import spambind.filters.filtersToolkit as ftk 

136 

137 im = im.astype("<f4") 

138 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

162 

163 # gaussian 

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

165 

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

167 gradient = numpy.gradient(im) 

168 

169 # second derivate 

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

171 hzz = tmp[0] 

172 hzy = tmp[1] 

173 hzx = tmp[2] 

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

175 hyz = tmp[0] 

176 hyy = tmp[1] 

177 hyx = tmp[2] 

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

179 hxz = tmp[0] 

180 hxy = tmp[1] 

181 hxx = tmp[2] 

182 

183 del tmp, gradient 

184 

185 # run eigen solver for each pixel!!! 

186 ftk.hessian(hzz, hzy, hzx, hyz, hyy, hyx, hxz, hxy, hxx, eigValA, eigValB, eigValC, eigVecAz, eigVecAy, eigVecAx, eigVecBz, eigVecBy, eigVecBx, eigVecCz, eigVecCy, eigVecCx) 

187 

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

189 

190 

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

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

193# """ 

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

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

196# PARAMETERS: 

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

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

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

200# RETURNS: 

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

202# HISTORY: 

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

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

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

206# """ 

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

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

209# # Step 2: structutral element 

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

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

212# if prlv>5: 

213# import progressbar 

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

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

216# bar.start() 

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

218# # convert structural element coordinates into shift to apply 

219# shift_to_apply = c-structural_element_size 

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

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

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

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

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

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

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

227# return o_im 

228# 

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

230# """ 

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

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

233# PARAMETERS: 

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

235# RETURNS: 

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

237# HISTORY: 

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

239# """ 

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

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