Coverage for larch/math/tomography.py: 21%

101 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-10-16 21:04 +0000

1''' 

2This module defines a functions necessary for tomography calculations. 

3 

4Authors/Modifications: 

5---------------------- 

6* Margaret Koker, koker@cars.uchicago.edu 

7''' 

8 

9import logging 

10logger = logging.getLogger(__name__) 

11logger.level = logging.ERROR 

12logger = logging.getLogger('tomopy.recon') 

13logger.level = logging.ERROR 

14 

15import numpy as np 

16from scipy.optimize import leastsq, minimize 

17 

18HAS_tomopy = False 

19try: 

20 import tomopy 

21 HAS_tomopy = True 

22except ImportError: 

23 pass 

24 

25TAU = 2.0 * np.pi 

26 

27# GLOBAL VARIABLES 

28TOMOPY_ALG = ['gridrec', 'art', 'bart', 'mlem', 'osem', 'ospml_hybrid', 

29 'ospml_quad', 'pml_hybrid', 'pml_quad', 'sirt' ] 

30 

31TOMOPY_FILT = ['shepp', 'ramlak', 'butterworth','parzen', 'cosine', 'hann', 

32 'hamming', 'None'] 

33 

34def ensure_radians(a): 

35 """ensure angle data is in radians, not degrees, 

36 converts degrees to radians if a peak-to-peak > 32 or step size > 0.2 

37 """ 

38 if np.ptp(a) > 32 or np.diff(a).mean() > 0.20: 

39 a = np.radians(a) 

40 return a 

41 

42 

43def reshape_sinogram(A, x, omega): 

44 

45 ## == INPUTS == 

46 ## A : array from .get_roimap() 

47 ## x : x array for checking shape of A 

48 ## omega : omega array for checking shape of A 

49 ## 

50 ## == RETURNS == 

51 ## A : A in shape/format needed for tomopy reconstruction 

52 ## sinogram_order : flag/argument for tomopy reconstruction (shape dependent) 

53 

54 if len(x) == len(omega): 

55 print("Warning: guessing that 2nd axis is omega") 

56 # Cannot reorder sinogram based on length of positional 

57 # arrays when same length. Acceptable orders: 

58 # sinogram_order = False : sino = [ 2th , slice, X ] 

59 # sinogram_order = True : sino = [ slice , 2th , X ]''') 

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

61 A = A.reshape(1, A.shape[0], A.shape[1]) 

62 return A, True 

63 

64 if len(x) < 1 or len(omega) < 1: 

65 return A,False 

66 if len(A.shape) != 3: 

67 if len(A.shape) == 2: 

68 A = A.reshape(1,A.shape[0],A.shape[1]) 

69 

70 if len(A.shape) == 3: 

71 if len(x) == A.shape[0]: 

72 A = np.einsum('kij->ijk', A) 

73 if len(x) == A.shape[1]: 

74 A = np.einsum('ikj->ijk', A) 

75 sinogram_order = len(omega) == A.shape[1] 

76 

77 return A, sinogram_order 

78 

79def trim_sinogram(sino, x, omega, pixel_trim=10): 

80 if len(omega) == sino.shape[-1]: 

81 omega = omega[pixel_trim:-1*(pixel_trim+1)] 

82 elif len(x) == sino.shape[-1]: 

83 x = x[pixel_trim:-1*(pixel_trim+1)] 

84 

85 sino = sino[:,pixel_trim:-1*(pixel_trim+1)] 

86 

87 return sino, x, omega 

88 

89def find_tomo_center(sino, omega, center=None, tol=0.25, blur_weight=2.0, 

90 sinogram_order=True): 

91 

92 """find rotation axis center for a sinogram, 

93 mixing negative entropy (as tomopy uses) and other focusing scores 

94 

95 Arguments 

96 --------- 

97 sino : ndarray for sinogram 

98 omega: ndarray of angles in radians 

99 center: initial value for center [mid-point] 

100 tol: fit tolerance for center pixel [0.25] 

101 sinogram_order: bool for axis order of sinogram 

102 

103 Returns 

104 ------- 

105 pixel value for refined center 

106 

107 Notes 

108 ------ 

109 

110 The algormithm combines a few focusing scores from Y. Sun, S. Duthaler, and B. Nelson, 

111 MICROSCOPY RESEARCH AND TECHNIQUE 65:139–149 (2004) (doi: 10.1002/jemt.20118a) 

112 

113 For a reconstructed image `img` the variance is calculated as 

114 

115 blur = -((img - img.mean())**2).sum()/img.size 

116 

117 and is combined with negative-entropy is calculated as 

118 ioff = (img.max() - img.min())/25.0 

119 imin = img.min() - ioff 

120 imax = img.max() + ioff 

121 hist, _ = np.histogram(img, bins=512, range=[imin, imax]) 

122 hist = hist/(2*(imax-imin)) 

123 hist[np.where(hist==0)] = 1.e-20 

124 negent = -np.dot(hist, np.log(hist)) 

125 

126 

127 """ 

128 xmax = sino.shape[0] 

129 if sinogram_order: 

130 xmax = sino.shape[2] 

131 if center is None: 

132 center = xmax/2.0 

133 

134 omega = ensure_radians(omega) 

135 

136 img = tomopy.recon(sino, omega, center, 

137 sinogram_order=sinogram_order, 

138 algorithm='gridrec', filter_name='shepp') 

139 img = tomopy.circ_mask(img, axis=0) 

140 ioff = (img.max() - img.min())/25.0 

141 imin = img.min() - ioff 

142 imax = img.max() + ioff 

143 

144 out = minimize(center_score, center, method='Nelder-Mead', tol=tol, 

145 args=(sino, omega, blur_weight, sinogram_order, imin, imax)) 

146 return out.x[0] 

147 

148def center_score(center, sino, omega, blur_weight=2.0, sinogram_order=True, 

149 imin=None, imax=None, verbose=False): 

150 """Cost function used for the ``find_center`` routine: 

151 combines a few focusing scores from 

152 Y. Sun, S. Duthaler, and B. Nelson, 

153 MICROSCOPY RESEARCH AND TECHNIQUE 65:139–149 (2004) 

154 (doi: 10.1002/jemt.20118a) 

155 

156 name formula in paper python 

157 

158 blur (F 10) -((img - img.mean())**2).sum()/img.size 

159 

160 negent is calculated as 

161 hist, _ = np.histogram(img, bins=512, range=[imin, imax]) 

162 hist = hist/(2*(imax-imin)) 

163 hist[np.where(hist==0)] = 1.e-20 

164 negent = -np.dot(hist, np.log(hist)) 

165 

166 """ 

167 ns, nang, nx = sino.shape 

168 if isinstance(center, (list, tuple, np.ndarray)): 

169 center = center[0] 

170 

171 img = tomopy.recon(sino, ensure_radians(omega), center, 

172 sinogram_order=sinogram_order, 

173 algorithm='gridrec', filter_name='shepp') 

174 img = tomopy.circ_mask(img, axis=0) 

175 blur = -((img - img.mean())**2).sum()/img.size 

176 

177 if imin is None or imax is None: 

178 ioff = (img.max() - img.min())/25.0 

179 if imin is None: 

180 imin = img.min() - ioff 

181 if imax is None: 

182 imax = img.max() + ioff 

183 try: 

184 hist, _ = np.histogram(img, bins=512, range=[imin, imax]) 

185 hist = hist/(2*(imax-imin)) 

186 hist[np.where(hist==0)] = 1.e-20 

187 negent = -np.dot(hist, np.log(hist)) 

188 except: 

189 negent = blur 

190 if verbose: 

191 print("Center %.3f %13.5g, %13.5g" % (center, blur, negent)) 

192 return blur*blur_weight + negent 

193 

194 

195def tomo_reconstruction(sino, omega, algorithm='gridrec', 

196 filter_name='shepp', num_iter=1, center=None, 

197 refine_center=False, sinogram_order=True): 

198 ''' 

199 INPUT -> sino : slice, 2th, x OR 2th, slice, x (with flag sinogram_order=True/False) 

200 OUTPUT -> tomo : slice, x, y 

201 ''' 

202 if center is None: 

203 center = sino.shape[1]/2. 

204 

205 x = tomopy.init_tomo(sino, sinogram_order) 

206 nomega = len(omega) 

207 ns, nth, nx = sino.shape 

208 if nth > nomega: 

209 sino = sino[:, :nomega, :] 

210 romega = ensure_radians(omega) 

211 

212 if refine_center: 

213 center = find_tomo_center(sino, romega, center=center, 

214 sinogram_order=sinogram_order) 

215 print(">> Refine Center done>> ", center, sinogram_order) 

216 algorithm = algorithm.lower() 

217 recon_kws = {} 

218 if algorithm.startswith('gridr'): 

219 recon_kws['filter_name'] = filter_name 

220 else: 

221 recon_kws['num_iter'] = num_iter 

222 

223 tomo = tomopy.recon(sino, romega, algorithm=algorithm, 

224 center=center, sinogram_order=sinogram_order, **recon_kws) 

225 return center, tomo*(sino.mean()/tomo.mean())