Coverage for crip\postprocess.py: 68%

1''' 

2 Postprocess module of crip. 

3 

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

5''' 

6 

7import numpy as np 

8 

9from .shared import * 

10from ._typing import * 

11from .utils import * 

12 

13 

14def fovCropRadius(sid: float, sdd: float, detectorWidth: float, pixelSize: float, roundOff=True) -> Or[int, float]: 

15 ''' Compute the radius [pixel] of the valid FOV of the reconstruction. 

16 `sid` and `sdd` are Source-Isocenter-Distance and Source-Detector-Distance, respectively. 

17 `detectorWidth` is the width of the detector, i.e., `#elements * elementWidth`. 

18 `pixelSize` is the pixel size of the reconstructed image. 

19 Recommended length unit is [mm]. 

20 ''' 

21 halfDW = detectorWidth / 2 

22 L = np.sqrt(halfDW**2 + sdd**2) 

23 

24 r1 = sid * np.arcsin(halfDW / L) # table as arc 

25 r2 = sid / (sdd / detectorWidth) / 2 # table as plane 

26 r3 = sid / (L / halfDW) 

27 

28 r = min(r1, r2, r3) / pixelSize 

29 

30 return round(r) if roundOff else r 

31 

32 

33@ConvertListNDArray 

34def fovCrop(img: TwoOrThreeD, radius: int, fill: Or[int, float] = 0) -> TwoOrThreeD: 

35 ''' Crop a circle FOV on reconstructed image `img` with `radius` [pixel] 

36 and `fill` outside FOV. 

37 ''' 

38 cripAssert(radius >= 1 and isInt(radius), 'Invalid radius. Radius should be positive int.') 

39 cripAssert(is2or3D(img), 'img should be 2D or 3D.') 

40 

41 N, M = img.shape[-2:] 

42 x = np.array(range(N), dtype=DefaultFloatDType) - N / 2 - 0.5 

43 y = np.array(range(M), dtype=DefaultFloatDType) - M / 2 - 0.5 

44 xx, yy = np.meshgrid(x, y) 

45 outside = xx**2 + yy**2 > radius**2 

46 

47 cropped = img.copy() 

48 if is3D(img): 

49 cropped[:, outside] = fill 

50 else: 

51 cropped[outside] = fill 

52 

53 return cropped 

54 

55 

56@ConvertListNDArray 

57def muToHU(image: TwoOrThreeD, muWater: float, b=1000) -> TwoOrThreeD: 

58 ''' Convert mu to HU using `HU = (mu - muWater) / muWater * b` 

59 ''' 

60 return (image - muWater) / muWater * b 

61 

62 

63@ConvertListNDArray 

64def huToMu(image: TwoOrThreeD, muWater: float, b=1000) -> TwoOrThreeD: 

65 ''' Convert HU to mu (invert `muToHU`). 

66 ''' 

67 return image / b * muWater + muWater 

68 

69 

70@ConvertListNDArray 

71def huNoRescale(image: TwoOrThreeD, b: float = -1000, k: float = 1) -> TwoOrThreeD: 

72 ''' Invert the rescale-slope (y = kx + b) of HU value to get linear relationship between HU and mu. 

73 ''' 

74 return (image - b) / k 

75 

76 

77@ConvertListNDArray 

78def postlogsToRaws(postlogs: TwoOrThreeD, flat: Or[TwoD, float]) -> TwoOrThreeD: 

79 ''' Invert `postlog` images to the original raw projections according to `flat` field. 

80 ''' 

81 return np.exp(-postlogs) * flat