search

CERN / TIGRE

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

Unexpected half-beam reconstruction result #384

Closed Enorielle closed 8 months ago

Enorielle commented 2 years ago

Hello,

I'm having issues reconstructing CT data from a half-beam scan. My question may not be specific to TIGRE so i apologize in advance. I'm rather new to CT in general, and to TIGRE in particular...

Many thanks in advance for your help.

Enorielle

Expected Behavior

When I reconstruct the data with RTK, here is what I obtain (and what i expect):

image

Actual Behavior

Performing the same reconstruction with TIGRE, using the FDK algorithm, here is what I get:

image

After the initial surprise and a bit of digging, I realized that Wang weighting is not (yet?) implemented in the Python version of TIGRE. I then moved to using iterative algorithms (and in this particular case, CGLS 30 iterations) and here is the result:

image

As you can see, it is better. However, I still have the central ring artifact, characteristic of half-beam scans not properly accounted for. I'm very surprised by this result as I expected iterative algorithms to not need any particular weighting to account for half-beam scans.

I tried looking around in the documentation but could not find any explanation for this behavior. Any pointers would be much appreciated!

Code to reproduce the problem (If applicable)

Here is the code I use. CT_geo is a Python dictionary containing all relevant CT parameters. ` 

geo = tigre.geometry()
geo.DSD = CT_geo['FD']
geo.DSO = CT_geo['FO']

geo.nDetector = np.array([CT_geo['HEIGHT'], CT_geo['WIDTH']])
geo.dDetector = np.array([CT_geo['PIX'], CT_geo['PIX']])
geo.sDetector = geo.nDetector * geo.dDetector
angles = np.linspace(CT_geo['THETA_START'], math.radians(CT_geo['THETA_STOP']), CT_geo['NPROJS'], False)

geo.nVoxel = np.array([CT_geo['NVOX_Z'], CT_geo['NVOX_Y'], CT_geo['NVOX_X']])
geo.dVoxel = np.array([CT_geo['VOX'], CT_geo['VOX'], CT_geo['VOX']])
geo.sVoxel = geo.nVoxel * geo.dVoxel

geo.offOrigin = np.array([0,0,0])
geo.offDetector = np.array([0,0])
geo.accuracy = 0.5

geo.COR = -CT_geo['COR'] 
geo.mode = "cone"

im_cgls      = algs.cgls(projections, geo, angles,30)
recons_fdk = algs.fdk(projections, geo, angles)

`

Specifications

AnderBiguri commented 2 years ago

Hi!

Thanks for the issue!

So, you are right in most cases indeed!

TIGRE python is 95% complete and wang weigths is one of the missing pieces. Unfortunately I am super busy with other things, but if you feel brave enough, adding them should be quite a simple thing: Take the code from here (full functions lower in the same file) https://github.com/CERN/TIGRE/blob/master/MATLAB/Algorithms/FDK.m#L46 And add it (translate it) here: https://github.com/CERN/TIGRE/blob/master/Python/tigre/algorithms/single_pass_algorithms.py#L81

Its not a lot of code, just some weights in the horizontal direction. Translation should be practically a copy-paste with minor changes.

Now, indeed Iterative algorithms should not suffer from this. And in some way, they don't, if you compare both of your results. However, nuances: Sometimes corrections (e.g. from flat-field correction) are worse in the corners of the detector. So the recon may be good for your data, its the data that its a bit off. In fact, if you do the same test but with synthetic data, you will see that the rings are not there. Another comment: its been shown in few papers that some iterative algorithms (I don't think CGLS is one of those, due to how the update happens mathematically) can benefit from wang weights too. See #376 that we are trying to add and comments in #373 with papers.

But I guess the real question is: what should you do? My best advice is to write the Wang weigths for python. They are easy and you have code to copy, and its the best way you will ensure good quality results. You can also try other iteratives, such as OS-SART. CGLS tends to diverge when data is misaligned with the model, much easier than OS-SART does.

Enorielle commented 2 years ago

Hi Ander,

First of all: thanks a lot for the quick reply !

Regarding Wang weights, I already started implementing them for the python version of TIGRE. If you want, I could share when I'm done. Just keep in my mind that i'm no python professional !

Thanks for the explanation. Actually, i do not use the full detector size. I crop the borders to avoid corner issues... but it might no be enough. I launched a simulation with synthetic data just to check for this effect nevertheless.

To be thorough, I tested other algorithms but the rings remain, just or less pronounced depending on the algorithm used.

Cheers

AnderBiguri commented 2 years ago

@Enorielle yes you are right. #362 also shows that, you are correct. I am a bit confused on why this happens admittedly. It doesn't always happen, so I am not sure. For some algorithms, the V matrix can be the culprit, but CGLS does not have it, so unsure really.

This requires some deep investigation, but unfortunately I don't have the time now for it. worrying nevertheless.

In any case, applying wang weigths to the iterative algorithms too seem to mitigate much of the effect.

hyaihjq commented 10 months ago

Hello,

I'm having issues reconstructing CT data from a half-beam scan. My question may not be specific to TIGRE so i apologize in advance. I'm rather new to CT in general, and to TIGRE in particular...

Many thanks in advance for your help.

Enorielle

Expected Behavior

When I reconstruct the data with RTK, here is what I obtain (and what i expect):

image

Actual Behavior

Performing the same reconstruction with TIGRE, using the FDK algorithm, here is what I get:

image

After the initial surprise and a bit of digging, I realized that Wang weighting is not (yet?) implemented in the Python version of TIGRE. I then moved to using iterative algorithms (and in this particular case, CGLS 30 iterations) and here is the result:

image

As you can see, it is better. However, I still have the central ring artifact, characteristic of half-beam scans not properly accounted for. I'm very surprised by this result as I expected iterative algorithms to not need any particular weighting to account for half-beam scans.

I tried looking around in the documentation but could not find any explanation for this behavior. Any pointers would be much appreciated!

Code to reproduce the problem (If applicable)

Here is the code I use. CT_geo is a Python dictionary containing all relevant CT parameters. `

geo = tigre.geometry()
geo.DSD = CT_geo['FD']
geo.DSO = CT_geo['FO']

geo.nDetector = np.array([CT_geo['HEIGHT'], CT_geo['WIDTH']])
geo.dDetector = np.array([CT_geo['PIX'], CT_geo['PIX']])
geo.sDetector = geo.nDetector * geo.dDetector
angles = np.linspace(CT_geo['THETA_START'], math.radians(CT_geo['THETA_STOP']), CT_geo['NPROJS'], False)

geo.nVoxel = np.array([CT_geo['NVOX_Z'], CT_geo['NVOX_Y'], CT_geo['NVOX_X']])
geo.dVoxel = np.array([CT_geo['VOX'], CT_geo['VOX'], CT_geo['VOX']])
geo.sVoxel = geo.nVoxel * geo.dVoxel

geo.offOrigin = np.array([0,0,0])
geo.offDetector = np.array([0,0])
geo.accuracy = 0.5

geo.COR = -CT_geo['COR'] 
geo.mode = "cone"

im_cgls      = algs.cgls(projections, geo, angles,30)
recons_fdk = algs.fdk(projections, geo, angles)

`

Specifications

  • Python version: 3.9.7
  • OS: Debian GNU/Linux 10
  • CUDA version: 11.4

Can you share the RTK implementation code? I encountered the same problem. slice

AnderBiguri commented 10 months ago

@hyaihjq just to be clear, this is almost certainly that the geometry being given to Tigre is not correct.

Likely angles, shift of the detector, are either not given in the appropriate units, or with an opposite sign.

Perhaps is the Wang weights, but not sure.

hyaihjq commented 10 months ago

@hyaihjq just to be clear, this is almost certainly that the geometry being given to Tigre is not correct.

Likely angles, shift of the detector, are either not given in the appropriate units, or with an opposite sign.

Perhaps is the Wang weights, but not sure.

Thank you for your reply. I will try adjusting the parameters to see if it can be improved

Enorielle commented 10 months ago

@hyaihjq : should you still need the RTK source code, it is available on their website https://www.openrtk.org/Doxygen/index.html

hyaihjq commented 10 months ago

@hyaihjq : should you still need the RTK source code, it is available on their website https://www.openrtk.org/Doxygen/index.html

There are very few RTK related examples, and I have also referred to the website you sent, but I did not obtain the correct results. If it is convenient, can you provide me with the code for reference? Could you please send it to this email “hyhjq@mail.ustc.edu.cn”. Thank you!!!

Enorielle commented 10 months ago

I'm really sorry but as we embedded the code to RTK in some proprietary code, I'm not allowed to share it...

In our case, we used rtksimulatedgeometry to construct the xml geomtery file and rtkfdk for the reconstruction.

Best of luck!

hyaihjq commented 10 months ago

I'm really sorry but as we embedded the code to RTK in some proprietary code, I'm not allowed to share it...

In our case, we used rtksimulatedgeometry to construct the xml geomtery file and rtkfdk for the reconstruction.

Best of luck!

thanks for your reply，i encounter the same problem by use "https://github.com/astra-toolbox"，I suspect there may be a problem with the projection data preprocessing. Can you share data preprocessing methods if it's convenient?

AnderBiguri commented 10 months ago

@hyaihjq if you want help with how to do it in Tigre, please open a new issue with a description of your scan, code and result. For example, if you have shifted detector, you need to describe this in the geometry appropriately, in any toolbox that you may use.

Making the geometry you provide to Tigre (or Astra or anything) the exact real one is the most important thing.

hyaihjq commented 8 months ago

@AnderBiguri Hi AnderBiguri, I found that the order of the angle and projection data is reversed，and i apply the “wang weighting” to the projection data. slice1 How to calculate the radius of a circle to crop the surrounding area.

AnderBiguri commented 8 months ago

@hyaihjq, there is a function called cropCBCT, is that what you are looking for?

hyaihjq commented 8 months ago

@hyaihjq, there is a function called cropCBCT, is that what you are looking for?

def cropCBCT(img, geo):
cropR = (geo.sDetector[0] / 2 * geo.DSO) / geo.DSD
maxD = min(geo.nVoxel - 1) / 2
cropR = min([cropR / geo.dVoxel[0], maxD]) x, y = np.meshgrid(range(img.shape[1]), range(img.shape[2]))
inM = (x - img.shape[1] / 2) 2 + (y - img.shape[2] / 2) 2 < cropR ** 2
img = np.multiply(img, inM)
return img

image

i translate matlab to python，i got this，it's false

AnderBiguri commented 8 months ago

@hyaihjq try:

def cropCBCT(img, geo):
      cropR = img.shape[1]/2
      x, y = np.meshgrid(range(img.shape[1]), range(img.shape[2]))
      inM = (x - img.shape[1] / 2) ** 2 + (y - img.shape[2] / 2) ** 2 < cropR ** 2
      img = np.multiply(img, inM)
      return img

Should work for your case.

hyaihjq commented 8 months ago

@hyaihjq try:

def cropCBCT(img, geo):
      cropR = img.shape[1]/2
      x, y = np.meshgrid(range(img.shape[1]), range(img.shape[2]))
      inM = (x - img.shape[1] / 2) ** 2 + (y - img.shape[2] / 2) ** 2 < cropR ** 2
      img = np.multiply(img, inM)
      return img

Should work for your case.

image

AnderBiguri commented 8 months ago

Feel free do adjust cropR as desired :) Glad it worked

hyaihjq commented 8 months ago

Feel free do adjust cropR as desired :) Glad it worked

sure, thanks for your help