CERN / TIGRE

TIGRE: Tomographic Iterative GPU-based Reconstruction Toolbox
BSD 3-Clause "New" or "Revised" License
563 stars 186 forks source link

Difference in values before projection and after reconstruction #181

Closed Elvira-Zainulina closed 4 years ago

Elvira-Zainulina commented 4 years ago

Hello,

I wanted to apply TIGRE (python) for the reconstruction of CT projections but the values in the reconstructed scans differed from the expected values significantly. So, I decided to conduct the experiment. The code of the experiment is below.

mins = []
maxs = []
for i in tqdm(range(100, 10000, 100)):
  a = np.zeros((120, 512, 512), dtype=np.float32)
  a[10:110, 15:35, 15:35] = i
  a[20:100, 22:28, 22:28] = 0.0
  prjs = Ax(a, geo, angles)  # geo.mode = "cone", the type of projection is spiral (helical)
  b = algs.fdk(prjs, geo, angles)
  mins.append(-b.min())
  maxs.append(a.max() - b.max())

Then I have plotted the mins and maxs and obtained the following plot, where the orange line is maxs and the blue line is mins

diff.png

Is this is expected behavior? If it is then can you advise something that can preserve a range of values of the image before projection and after reconstruction?

Specifications

AnderBiguri commented 4 years ago

Can you share your geometry too? You say it's spiral, but FDK is not a valid algorithm for helical CT

Elvira-Zainulina commented 4 years ago

Here is my geometry.

class Geometry(tigre.utilities.geometry.Geometry):

    def __init__(self):
        # VARIABLE                                          DESCRIPTION                    UNITS
        # -------------------------------------------------------------------------------------
        self.DSD = np.array(DSD)                                   # Distance Source Detector      (mm)
        self.DSO = np.array(DSO)                                  # Distance Source Origin        (mm)
        # Detector parameters
        self.nDetector = np.array((nRows, nCols))               # number of pixels              (px)
        self.dDetector = np.array((mmRows, mmCols))               # size of each pixel            (mm)
        self.sDetector = self.nDetector * self.dDetector    # total size of the detector    (mm)
        # self.rotDetector = np.array((0, 0, 0)) 

        # Image parameters
        self.nVoxel = np.array((120, 512, 512))             # number of voxels              (vx)
        self.sVoxel = np.array((120, 512 * 0.5859375, 512 * 0.5859375))             # total size of the image       (mm)
        self.dVoxel = self.sVoxel/self.nVoxel               # size of each voxel            (mm)
        # Offsets
        self.offDetector = np.array((-offDetRow, -offDetCol))                 # Offset of Detector            (mm)
        self.offOrigin = np.stack([offsets_zero, offsets_zero, offOriginZ], axis = 0).T

        # Auxiliary
        self.accuracy = 0.5                                 # Accuracy of FWD proj          (vx/sample)
        # Mode
        self.mode = 'cone'                                  # parallel, cone                ...
        self.filter = None

I have taken values for geometry from real CT projection geometry. offOriginZ is array of the offset along Z axis.

AnderBiguri commented 4 years ago

@Elvira-Zainulina can you see the effect you describe if you remove the origin offset in the experiment above?

I think this is because FDK is not a valid algorithm for Z offsets in the image, thus you may be seeing an unwanted effect there. Can you test that please?

Elvira-Zainulina commented 4 years ago

I have removed the origin offset but the difference is still rather big. Here is the result, where the orange line is -maxs and the blue line is mins.

diff0.png

Actually, I have made a typo in the first comment. The orange line on the plot is also -maxs, not maxs.

AnderBiguri commented 4 years ago

Oh well, I was wrong then!

Let me investigate...

AnderBiguri commented 4 years ago

@Elvira-Zainulina can you do another test for me? I have the feeling that you may be comparing noise. FDK does bring artifacts, and it makes sense that the value of those artifacts change with the value of the pixel values that create them. Bigger image value, bigger absolute error, of course.

So, can you compute the relative error instead? And can you confirm that voxel values are properly reconstructed? i.e. that b[10:110, 15:35, 15:35] is approximately i each time.

Elvira-Zainulina commented 4 years ago

Yes, you are right. The max values corresponded to the noise. Well, I have calculated relative error for max values. (In my setting it approximately equals to 0.3309.) Also, I decided to calculate MAE and MSE between the regions of images containing the square: RMSE / i ≈ 0.35, MAE / i ≈ 0.15 for all considered i.

Then I have checked the values of the reconstructed central slice. Here are the results (normalization by min and max values of the original image):

slice.png

It preserves for all considered i. 30% difference between the original image and the reconstructed is rather big for my problem. So, is it possible to avoid this difference or mitigate it somehow?

AnderBiguri commented 4 years ago

@Elvira-Zainulina Ah, welcome to the complicated world of CT reconstruction then! A big question you ask!

Few comments on that: FDK is not the best algorithm.... FDK is a good algorithm if you have infinite noiseless projections. Of course, this is never true. It particularly fails on very sharp edges as the one you show there, and it shows artifacts very early with quite decent data. Of course, as said before, FDK also only works for circular (not helical) CT. If you can obtain approximately ~3000 projections, its likely that FDK will result in a very noiseless image, but otherwise, you risk the effect you are seeing. If you want to try this, do the same experiment with different amount of angles.

One of the main purposes of releasing TIGRE is indeed showcasing that FDK is not a great algorithm, and providing alternative algorithms. You can find ways of calling different algorithms on the demo files, so you can compare their result with FDK.

I suggest giving OS_SART a try (demo here). Its a simple algorithm, which means its not amazing, but should be able to reconstruct images with much less error.

Let me know if you have more questions.

AnderBiguri commented 4 years ago

Hi @Elvira-Zainulina !

Do you have more questions? Otherwise I will close the issue. But feel free to send me an email or open new issues in the future.

Elvira-Zainulina commented 4 years ago

No, I don't have any questions now. Thank you for your help!

AnderBiguri commented 4 years ago

@Elvira-Zainulina no worries! Feel free to ask questions anytime!