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
« 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.
4Authors/Modifications:
5----------------------
6* Margaret Koker, koker@cars.uchicago.edu
7'''
9import logging
10logger = logging.getLogger(__name__)
11logger.level = logging.ERROR
12logger = logging.getLogger('tomopy.recon')
13logger.level = logging.ERROR
15import numpy as np
16from scipy.optimize import leastsq, minimize
18HAS_tomopy = False
19try:
20 import tomopy
21 HAS_tomopy = True
22except ImportError:
23 pass
25TAU = 2.0 * np.pi
27# GLOBAL VARIABLES
28TOMOPY_ALG = ['gridrec', 'art', 'bart', 'mlem', 'osem', 'ospml_hybrid',
29 'ospml_quad', 'pml_hybrid', 'pml_quad', 'sirt' ]
31TOMOPY_FILT = ['shepp', 'ramlak', 'butterworth','parzen', 'cosine', 'hann',
32 'hamming', 'None']
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
43def reshape_sinogram(A, x, omega):
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)
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
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])
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]
77 return A, sinogram_order
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)]
85 sino = sino[:,pixel_trim:-1*(pixel_trim+1)]
87 return sino, x, omega
89def find_tomo_center(sino, omega, center=None, tol=0.25, blur_weight=2.0,
90 sinogram_order=True):
92 """find rotation axis center for a sinogram,
93 mixing negative entropy (as tomopy uses) and other focusing scores
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
103 Returns
104 -------
105 pixel value for refined center
107 Notes
108 ------
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)
113 For a reconstructed image `img` the variance is calculated as
115 blur = -((img - img.mean())**2).sum()/img.size
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))
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
134 omega = ensure_radians(omega)
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
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]
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)
156 name formula in paper python
158 blur (F 10) -((img - img.mean())**2).sum()/img.size
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))
166 """
167 ns, nang, nx = sino.shape
168 if isinstance(center, (list, tuple, np.ndarray)):
169 center = center[0]
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
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
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.
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)
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
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())