Coverage for crip\spec.py: 27%

1''' 

2 Spectrum CT module of crip. 

3 

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

5''' 

6 

7import numpy as np 

8from scipy.ndimage import uniform_filter 

9 

10from .postprocess import gaussianSmooth 

11from .utils import ConvertListNDArray, cripAssert, cripWarning, is2D, isOfSameShape 

12from ._typing import * 

13from .physics import Atten, Spectrum, computeAttenedSpectrum 

14from .shared import applyPolyV2D2, fitPolyV2D2 

15 

16 

17def deDecompProjCoeff(lowSpec: Spectrum, highSpec: Spectrum, base1: Atten, len1: Or[NDArray, Iterable], base2: Atten, 

18 len2: Or[NDArray, Iterable]): 

19 ''' Compute the decomposing coefficient (Order 2 with bias term) of two spectra onto two material bases 

20 used in the projection domain. 

21 ''' 

22 

23 def computePostlog(spec, attens, Ls): 

24 flat = np.sum(spec.omega) 

25 attenedSpec = computeAttenedSpectrum(spec, attens, Ls) 

26 

27 return -np.log(attenedSpec.omega / flat) 

28 

29 lenCombo = [] 

30 postlogLow = [] 

31 postlogHigh = [] 

32 for i in len1: 

33 for j in len2: 

34 lenCombo.append([i, j]) 

35 postlogLow.append(computePostlog(lowSpec, [base1, base2], [i, j])) 

36 postlogHigh.append(computePostlog(highSpec, [base1, base2], [i, j])) 

37 

38 beta, gamma = fitPolyV2D2(np.array(postlogLow), np.array(postlogHigh), np.array(lenCombo)).T 

39 

40 return beta, gamma 

41 

42 

43@ConvertListNDArray 

44def deDecompProj(lowProj: TwoOrThreeD, highProj: TwoOrThreeD, coeff1: NDArray, 

45 coeff2: NDArray) -> Tuple[TwoOrThreeD, TwoOrThreeD]: 

46 ''' Perform dual-energy decompose in projection domain point-by-point using coeffs. 

47 Coefficients can be generated using @see `deDecompProjCoeff`. 

48 ''' 

49 cripAssert(isOfSameShape(lowProj, highProj), 'Two projection sets should have same shape.') 

50 cripAssert( 

51 len(coeff1) == 6 and len(coeff2) == 6, 

52 'Decomposing coefficients should have length 6 (2 variable, order 2 with bias).') 

53 

54 return applyPolyV2D2(coeff1, lowProj, highProj), applyPolyV2D2(coeff2, lowProj, highProj) 

55 

56 

57@ConvertListNDArray 

58def deDecompRecon(low: TwoOrThreeD, 

59 high: TwoOrThreeD, 

60 muB1: List[float], 

61 muB2: List[float], 

62 checkCond: bool = True) -> TwoOrThreeD: 

63 ''' Perform Dual-Energy Two-Material decomposition in image domain. 

64 The outputs are decomposing coefficients. Used values can be LAC (mu), or HU+1000. 

65 `muB*` stores the value of base* in (low-energy, high-energy) order. 

66 ''' 

67 cripAssert(isOfSameShape(low, high), 'Two volumes should have same shape.') 

68 COND_TOLERANCE = 1000 

69 

70 A = np.array([ 

71 [muB1[0], muB2[0]], 

72 [muB1[1], muB2[1]], 

73 ]) 

74 checkCond and cripWarning( 

75 np.linalg.cond(A) <= COND_TOLERANCE, 'Material decomposition matrix possesses high condition number.') 

76 M = np.linalg.inv(A) 

77 

78 def decomp1(low, high): 

79 c1 = M[0, 0] * low + M[0, 1] * high 

80 c2 = M[1, 0] * low + M[1, 1] * high 

81 return c1, c2 

82 

83 if is2D(low): 

84 return decomp1(low, high) 

85 else: 

86 return np.array(list(map(lambda args: decomp1(*args), zip(low, high)))).transpose((1, 0, 2, 3)) 

87 

88 

89@ConvertListNDArray 

90def teDecompRecon(low: TwoOrThreeD, mid: TwoOrThreeD, high: TwoOrThreeD, muB1: List[float], muB2: List[float], 

91 muB3: List[float]) -> TwoOrThreeD: 

92 ''' Perform Triple-Energy Three-Material decomposition in image domain. 

93 The outputs are decomposing coefficients. Used values can be LAC (mu), or HU+1000. 

94 `muB*` stores the value of base* in low-energy, mid-energy, high-energy order. 

95 ''' 

96 cripAssert(isOfSameShape(low, mid) and isOfSameShape(low, high), 'Volumes should have same shape.') 

97 

98 A = np.array([ 

99 [muB1[0], muB2[0], muB3[0]], 

100 [muB1[1], muB2[1], muB3[1]], 

101 [muB1[2], muB2[2], muB3[2]], 

102 ]) 

103 M = np.linalg.inv(A) 

104 

105 def decomp1(low, mid, high): 

106 c1 = M[0, 0] * low + M[0, 1] * mid + M[0, 2] * high 

107 c2 = M[1, 0] * low + M[1, 1] * mid + M[1, 2] * high 

108 c3 = M[2, 0] * low + M[2, 1] * mid + M[2, 2] * high 

109 

110 return c1, c2, c3 

111 

112 if is2D(low): 

113 return decomp1(low, mid, high) 

114 else: 

115 return np.array(list(map(lambda args: decomp1(*args), zip(low, mid, high)))).transpose((1, 0, 2, 3)) 

116 

117 

118@ConvertListNDArray 

119def deDecompReconVolCon(low: TwoOrThreeD, high: TwoOrThreeD, muB1: List[float], muB2: List[float], 

120 muB3: List[float]) -> TwoOrThreeD: 

121 ''' Perform Dual-Energy Three-Material decomposition with Volume Conservation constraint in image domain. 

122 `muB*` stores the value of base* in low-energy, high-energy order. 

123 ''' 

124 pseudoMid = np.zeros_like(low) 

125 

126 return teDecompRecon(low, pseudoMid, high, [*muB1, 1.0], [*muB2, 1.0], [*muB3, 1.0]) 

127 

128 

129def genMaterialPhantom(img, zsmooth=3, sigma=1, l=80, h=300, base=1000): 

130 ''' Generate the phantom of material bases (one is water) from SECT using soft-thresholding. 

131 zsmooth: smooth window in slice direction. 

132 sigma: Gaussian smooth sigma in single slice. 

133 [l, h] defines the fuzzy range of another material, e.g., bone. 

134 base: the reference HU of another material. 

135 ''' 

136 cripAssert(np.min(img)) < 0 # HU 

137 

138 def softThreshold(img: NDArray, l, h, mode='lower'): 

139 shape = img.shape 

140 img = img.flatten() 

141 

142 lower = img < l 

143 upper = img >= h 

144 

145 transitional = (img >= l) * (img < h) 

146 if mode == 'upper': 

147 transitional = transitional * (img - l) / (h - l) 

148 res = upper + transitional 

149 elif mode == 'lower': 

150 transitional = transitional * (h - img) / (h - l) 

151 res = lower + transitional 

152 

153 return res.reshape(shape) 

154 

155 kernel = (zsmooth, 1, 1) if zsmooth is not None else (1, 1) 

156 img = uniform_filter(img, kernel, mode='reflect') 

157 img = gaussianSmooth(img, sigma) 

158 

159 water = (img + 1000) / 1000 * softThreshold(img, l, h, 'lower') 

160 b2 = (img + 1000) / (base + 1000) * softThreshold(img, l, h, 'upper') 

161 

162 return water, b2 

163 

164 

165def compose2(b1: float, b2: float, v1: NDArray, v2: NDArray) -> NDArray: 

166 ''' Compose two vectors `(v1, v2)` by coeffs `(b1, b2)`. 

167 ''' 

168 return b1 * v1 + b2 * v2 

169 

170 

171def compose3(b1: float, b2: float, b3: float, v1: NDArray, v2: NDArray, v3: NDArray) -> NDArray: 

172 ''' Compose three vectors `(v1, v2, v3)` by coeffs `(b1, b2, b3)`. 

173 ''' 

174 return b1 * v1 + b2 * v2 + b3 * v3 

175 

176 

177def vmi2Mat(b1: TwoOrThreeD, b2, b1Mat: Atten, b2Mat: Atten, E: int): 

178 ''' Virtual Monoenergetic Imaging using two-material decomposition at energy `E` [keV]. 

179 ''' 

180 return b1 * b1Mat.mu[E] + b2 * b2Mat.mu[E] 

181 

182 

183def vmi3Mat(b1, b2, b3, b1Mat, b2Mat, b3Mat, E): 

184 ''' Virtual Monoenergetic Imaging using three-material decomposition. 

185 ''' 

186 return b1 * b1Mat.mu[E] + b2 * b2Mat.mu[E] + b3 * b3Mat.mu[E] 

187 

188 

189def deSubtration(low: TwoOrThreeD, high: TwoOrThreeD, kLow: float, kHigh: float = 1) -> TwoOrThreeD: 

190 ''' Dual-Energy subtraction for e.g. bone removal in CT slices or X-ray DR projections. 

191 ''' 

192 return kHigh * high - kLow * low 

193 

194 

195def vncBasis(contrastHULow: float, contrastHUHigh: float) -> NDArray: 

196 ''' Construct Virtual Non-Contrast basis on CT images with HU as value. 

197 ''' 

198 pass