Coverage for crip\preprocess.py: 84%

1''' 

2 Preprocess module of crip. 

3 

4 https://github.com/SEU-CT-Recon/crip 

5''' 

6 

7import numpy as np 

8import warnings 

9from scipy.interpolate import interpn 

10import cv2 

11 

12from .shared import * 

13from ._typing import * 

14from .utils import * 

15 

16warnings.simplefilter("ignore", DeprecationWarning) 

17 

18 

19@ConvertListNDArray 

20def averageProjections(projections: ThreeD) -> TwoD: 

21 ''' Average projections along the first axis (views) to get a single projection. 

22 ''' 

23 cripAssert(is3D(projections), f'`projections` should be 3D, but got {len(projections.shape)}-D.') 

24 

25 projections = asFloat(projections) 

26 res = projections.sum(axis=0) / projections.shape[0] 

27 

28 return res 

29 

30 

31@ConvertListNDArray 

32def correctFlatDarkField(projections: TwoOrThreeD, 

33 flat: Or[TwoD, ThreeD, float, None] = None, 

34 dark: Or[TwoD, ThreeD, float] = 0, 

35 fillNaN: Or[float, None] = 0) -> TwoOrThreeD: 

36 ''' Perform flat field (air) and dark field correction to get post-log projections. 

37 I.e., `- log [(X - D) / (F - D)]`. 

38 Usually `flat` and `dark` are 2D. 

39 If `flat` and `projections` are both 3D, perform view by view. 

40 If `flat` is None, estimate it by brightest pixel each view. 

41 ''' 

42 if flat is None: 

43 cripWarning(False, '`flat` is None. Use the maximum value of each view instead.') 

44 flat = np.max(projections.reshape((projections.shape[0], -1)), axis=1) 

45 flat = np.ones_like(projections) * flat[:, np.newaxis, np.newaxis] 

46 

47 if not isType(flat, np.ndarray): 

48 flat = np.ones_like(projections) * flat 

49 if not isType(dark, np.ndarray): 

50 dark = np.ones_like(projections) * dark 

51 

52 sampleProjection = projections if is2D(projections) else projections[0] 

53 

54 def checkShape(haystack, needle, needleName): 

55 cripAssert( 

56 isOfSameShape(haystack, needle), 

57 f'`projections` and `{needleName}` should have same shape, but got {haystack.shape} and {needle.shape}.') 

58 

59 checkShape(sampleProjection if is2D(flat) else projections, flat, 'flat') 

60 checkShape(sampleProjection if is2D(dark) else projections, dark, 'dark') 

61 

62 numerator = projections - dark 

63 denominator = flat - dark 

64 cripWarning(np.min(numerator > 0), 'Some `projections` values are not greater than zero after canceling `dark`.') 

65 cripWarning(np.min(denominator > 0), 'Some `flat` values are not greater than zero after canceling `dark`.') 

66 

67 res = -np.log(numerator / denominator) 

68 

69 if fillNaN is not None: 

70 res = np.nan_to_num(res, nan=fillNaN) 

71 

72 return res 

73 

74 

75@ConvertListNDArray 

76def projectionsToSinograms(projections: ThreeD): 

77 ''' Permute projections to sinograms by axes swapping `(views, h, w) -> (h, views, w)`. 

78 The width direction is along detector channels of a row. 

79 ''' 

80 cripAssert(is3D(projections), 'projections should be 3D.') 

81 

82 (views, h, w) = projections.shape 

83 sinograms = np.zeros((h, views, w), dtype=projections.dtype) 

84 for i in range(views): 

85 sinograms[:, i, :] = projections[i, :, :] 

86 

87 return sinograms 

88 

89 

90@ConvertListNDArray 

91def sinogramsToProjections(sinograms: ThreeD): 

92 ''' Permute sinograms back to projections by axes swapping `(h, views, w) -> (views, h, w)`. 

93 The width direction is along detector channels of a row. 

94 ''' 

95 cripAssert(is3D(sinograms), 'projections should be 3D.') 

96 

97 (h, views, w) = sinograms.shape 

98 projections = np.zeros((views, h, w), dtype=sinograms.dtype) 

99 for i in range(views): 

100 projections[i, :, :] = sinograms[:, i, :] 

101 

102 return projections 

103 

104 

105@ConvertListNDArray 

106def padImage(img: TwoOrThreeD, 

107 padding: Tuple[int, int, int, int], 

108 mode: str = 'reflect', 

109 decay: bool = False, 

110 cval: float = 0): 

111 ''' 

112 Pad each 2D image on four directions using symmetric `padding` (Up, Right, Down, Left). 

113 `mode` determines the border value, can be `edge`, `constant` (with `cval`), `reflect`. 

114 `decay` can be True to smooth padded border. 

115 ''' 

116 cripAssert(mode in ['edge', 'constant', 'reflect'], f'Invalid mode: {mode}.') 

117 

118 cosineDecay = lambda ascend, dot: np.cos( 

119 np.linspace(-np.pi / 2, 0, dot) if ascend else np.linspace(0, np.pi / 2, dot)), 

120 

121 h, w = getHnW(img) 

122 nPadU, nPadR, nPadD, nPadL = padding 

123 padH = h + nPadU + nPadD 

124 padW = w + nPadL + nPadR 

125 

126 def decayLR(xPad, w, nPadL, nPadR): 

127 xPad[:, 0:nPadL] *= cosineDecay(True, nPadL)[:] 

128 xPad[:, w - nPadR:w] *= cosineDecay(False, nPadR)[:] 

129 return xPad 

130 

131 def procOne(img): 

132 kwargs = {'constant_values': cval} if mode == 'constant' else {} 

133 xPad = np.pad(img, ((nPadU, nPadD), (nPadL, nPadR)), mode=mode, **kwargs) 

134 if decay: 

135 xPad = decayLR(xPad, padW, nPadL, nPadR) 

136 xPad = decayLR(xPad.T, padH, nPadU, nPadD) 

137 xPad = xPad.T 

138 return xPad 

139 

140 if is3D(img): 

141 res = [procOne(img[i, ...]) for i in range(img.shape[0])] 

142 return np.array(res) 

143 else: 

144 return procOne(img) 

145 

146 

147@ConvertListNDArray 

148def correctRingArtifactProjLi(sgm: TwoOrThreeD, sigma: float, ksize: Or[int, None] = None): 

149 ''' 

150 Apply the ring artifact correction method in projection domain (input postlog sinogram), 

151 using gaussian filter in sinogram detector direction [1]. 

152 [1] 李俊江,胡少兴,李保磊等.CT图像环状伪影校正方法[J].北京航空航天大学学报,2007,(11):1378-1382.DOI:10.13700/j.bh.1001-5965.2007.11.020 

153 https://kns.cnki.net/kcms2/article/abstract?v=xBNwvqFr00JMj5WzBbZMcj9N-kBm9Pi08enuNi8ymtGWVZuwGHdLWwttkgzSivkJSBVk0CVQZxo7DgBmujqhLCaVvvBMoig5RV0B4fDLUnjGQcCPo3O4KfAo4iX4vCwZ&uniplatform=NZKPT&flag=copy 

154 ''' 

155 ksize = ksize or int(2 * np.ceil(2 * sigma) + 1) 

156 kernel = np.squeeze(cv2.getGaussianKernel(ksize, sigma)) 

157 

158 def procOne(sgm: TwoD): 

159 Pc = np.mean(sgm, axis=0) 

160 Rc = np.convolve(Pc, kernel, mode='same') 

161 Ec = Pc - Rc 

162 return sgm - Ec[np.newaxis, :] 

163 

164 if is3D(sgm): 

165 res = np.array(list(map(procOne, sgm))) 

166 else: 

167 res = procOne(sgm) 

168 

169 return res 

170 

171 

172def fanToPara(sgm: TwoD, gammas: NDArray, betas: NDArray, sid: float, oThetas: Tuple[float], 

173 oLines: Tuple[float]) -> TwoD: 

174 ''' 

175 Re-order a Fan-Beam sinogram to Parallel-Beam's. 

176 `gammas` is fan angles from min to max [RAD], computed by `arctan(elementOffcenter / SDD)` for each element. 

177 `betas` is system rotation angles from min to max [RAD]. 

178 `sid` is Source-Isocenter-Distance [mm]. 

179 `oThetas` is output rotation angle range as (min, delta, max) tuple [RAD] (excluding max) 

180 `oLines` is output detector element physical locations range as (min, delta, max) tuple [mm] (excluding max) 

181 ``` 

182 /| <- gamma for detector element X 

183 / | <- SID 

184 / | 

185 ====X============ <- detector 

186 ^<--^ <- offcenter 

187 ''' 

188 nThetas = np.round((oThetas[2] - oThetas[0]) / oThetas[1]).astype(np.int32) 

189 nLines = np.round((oLines[2] - oLines[0]) / oLines[1]).astype(np.int32) 

190 

191 thetas1 = np.linspace(oThetas[0], oThetas[2], nThetas).reshape((1, -1)) 

192 R1 = np.linspace(oLines[0], oLines[2], nLines).reshape((1, -1)) 

193 

194 thetas = thetas1.T @ np.ones((1, nLines)) 

195 Rs = np.ones((nThetas, 1)) @ R1 

196 

197 oGammas = np.arcsin(Rs / sid) 

198 oBetas = thetas + oGammas - np.pi / 2 

199 

200 minBetas, maxBetas = np.min(betas), np.max(betas) 

201 inds = oBetas > maxBetas 

202 if len(inds) > 0: 

203 oBetas[inds] = oBetas[inds] - 2 * np.pi 

204 inds = oBetas < minBetas 

205 if len(inds) > 0: 

206 oBetas[inds] = oBetas[inds] + 2 * np.pi 

207 oBetas = np.minimum(oBetas, maxBetas) 

208 

209 interpolator = interpn((betas, gammas), sgm, (oBetas, oGammas), method='linear', bounds_error=False, fill_value=0) 

210 out = interpolator.reshape(oBetas.shape) 

211 

212 return out