Inaccurate FDK reconstruction for low voxel size
simonCtBern opened this issue · 21 comments
I'm trying to use tigre for reconstruction cone-beam CT images with variable scan geometries.
However, I encountered unexpected results for a test data set with simulated projections of a simple circular scan trajectory.
Expected Behavior
The quality of the reconstructed volume should not depend on "linear scaling of the object / projection magnification".
Actual Behavior
Reconstruction of the simulated scan using FDK gives good results for the following set of geometry parameters:
- Distance from source to detector (DSD) = 556.0 mm
- Distance from source to origin (DSO)= 200.0 mm
- Size of each voxel (dVoxel) = [0.07194245, ...] mm
However, reconstruction of the same data provides poor results if a larger magfinication is provided as geometry parameters:
- Distance from source to detector (DSD) = 556.0 mm
- Distance from source to origin (DSO)= 6.0 mm
- Size of each voxel (dVoxel) = [0.00215827, ...] mm
All other parameters (number of pixels, pixel size) and the projections are identical for both reconstructions:
- Number of pixels (nDetector) = [1200 2000]
- Size of each pixel (dDetector) = [0.2 0.2] mm
- Number of voxels (nVoxel) = [1200 2000 2000]
I tried changing the accuracy parameter of the geometry over a wide range of several orders of magnitudes without visible effect.
It seems that the is some "absolute numerical accuracy threshold" which creates problems at low voxel size?
Is there a way around this?
(Typing a wrong magnification as a workaround may not matter for an ideal circular scan trajectory, but my goal is to test variations/deviations in the geometry for all projections individually, which would be influenced by the magnification.)
Scan setup
The test data for the reconstruction is a simulated circular cone-beam scan of a hexagonal object with a hole in the center (roughly representing a hexagonal nut).
result with large DSO (low magnification):
result with low DSO (high magnification), "correct" simulation setup with bad reconstruction:
Specifications
- python version: 3.10.14
- OS: Windows 10 64 bit
- CUDA version: 11.7
- GPU: 2x GeForce GTX TITAN X (12 GB), compute capability 5.2
Hi @simonCtBern
Thanks for using TIGRE!
I need to be honest here: I understood almost nothing of what you are trying to experiment with here, nor the goals or how you changed the code. I also don't understand what I should be looking at in the images and what you expect to see. Can you share a little bit of context here, and some code?
My current guess is that you seem to be playing with ratios of pixels sizes in the detector vs voxel sizes in the image perhaps?
Just in case my guess was right:
You generally want to maitain a ratio of geo.dDetector=geo.dVoxel*geo.DSD/geo.DSO. This is ideally the one you want to accurately have image space for all the info in the detector, without wasting compute time. The way the code is made, particularly the backprojection, will suffer a little bit numerically at very large deviations from this ratio. If you want to know why I can easily point you to the corresponding resources that explain it (mostly my thesis).
Hi @AnderBiguri
Thanks for the quick reply!
I'm afraid I made my question sound more complicated than it is.
In a standard circular cone-beam CT, situations 1 and 2 in the following figure produce the same projections and the reconstructed volume for both situations must look identical (just with different voxel size geo.dVoxel), given that the size of the blue object is scaled according to its distance from the source (DSO 1 and DSO 2) and the detector remains unchanged:
However, it seems that using TIGRE, the reconstructed volume does not look the same, but there are numerical artifacts in situation 2 (second image in my initial post). We normally use a commercial software for reconstruction which does not produce these artifacts in situation 2.
I didn't share any code because I'm only using algs.fdk() without modifications and there is no error.
Is it possible that there is some form of numerical threshold in the FDK backprojection which produces bad results for small voxel size geo.dVoxel? Without knowing the details, I believe the algorithm calculates an intersection between x-rays and voxels and interpolates their projection on the detector..?
@simonCtBern Thanks, this makes more sense!
Are the images real scans, or simulated with TIGRE? For real scans this would not necesarily be true, of course, due to focus, and other real X-ray measurements.
Indeed, as you well point it out, this should not happen. Would you be able to share code/data to reproduce this myself? Ideally a cropped out version (as close to 2D as possible without removing the effects) so I can run it in a local computer, as my direct access to big computers for recon is a little bit limited these days.
Another test to run is to try to do Atb() directly on the sinos, rather than fdk. More to see if this is a numerical inaccuracy in the projectors or in the filtering. In any case, definetly the code to reproduce it would help me a lot pinpoint the issue and see what is off.
@AnderBiguri The screenshots were from projections simulated with aRTist, but it happened for real measured data as well (using simulated data was my problem searching approach).
I haven't tried running Atb() directly on the projections, but it would surprise me if it made a difference, since the (filtered or unfiltered) projections in my original datasets (i.e. not from TIGRE) don't depend on geo.DSO..?
Below, you can find code using TIGRE only for simulating and reconstructing a cube in almost 2D. On my system, the artifacts appear or not depending on the choice of geo.DSO:
geo.DSO = 100
geo.DSO = 6
geo.DSO = 1
#% testing reconstruction accuracy for small voxel sizes
import tigre
import tigre.algorithms as algs
import numpy as np
import matplotlib.pyplot as plt
#%% define projection parameters
detector_width = 2000
detector_height = 4
nb_projections = 1500
pixel_size = 0.2
# define the projection angles
angles = np.linspace(0, -2*np.pi, nb_projections, False).astype('float32')
# create a geomtry object
geo = tigre.geometry()
#%% Distance Source Object
geo.DSO = 1.0
# geo.DSO = 6.0
# geo.DSO = 100.0
#%% Distance Source Detector
geo.DSD = 600.0
# Detector parameters
geo.nDetector = np.array([detector_height, detector_width]) # number of pixels (px)
geo.dDetector = np.array([pixel_size, pixel_size]) # size of each pixel (mm)
geo.sDetector = geo.nDetector * geo.dDetector # total size of the detector (mm)
# Image parameters
geo.nVoxel = np.array([detector_height, detector_width, detector_width]) # number of voxels (vx)
geo.sVoxel = geo.nVoxel * pixel_size * geo.DSO / geo.DSD # total size of the image (mm)
geo.dVoxel = geo.sVoxel / geo.nVoxel # size of each voxel (mm)
# Offsets
geo.offOrigin = np.array([0, 0, 0]) # Offset of image from origin (mm)
geo.offDetector = np.array([0, 0]) # Offset of Detector (mm)
# Auxiliary
geo.accuracy = 0.1 # Variable to define accuracy of 'interpolated' projection.
# It defines the amount of samples per voxel. Recommended <=0.5 (vx/sample)
geo.rotDetector = np.array([0, 0, 0]) # Rotation of the detector, by
# X,Y and Z axis respectively. (rad)
# This can also be defined per angle
geo.mode = "cone" # Or 'parallel'. Geometry type.
#%% object definition
# create a volume with a cube in its center, filling 1/2 of each dimension
cube = np.zeros(
(detector_height, detector_width, detector_width), dtype='float32')
cube[(detector_height // 4):(detector_height * 3 // 4),
(detector_width // 4):(detector_width * 3 // 4),
(detector_width // 4):(detector_width * 3 // 4)] = 1
#%% Simulate forward projection.
projections = tigre.Ax(cube, geo, angles)
#%% FDK
volFDK = algs.fdk(projections, geo, angles, filter="shepp_logan")
#%% plot a slice
cmap = plt.cm.gray
volSlice = volFDK[volFDK.shape[0] // 2, :, :].copy()
norm = plt.Normalize(vmin=volSlice.min(), vmax=volSlice.max())
image = cmap(norm(volSlice))
plt.imshow(image)
Looking into this...
just FYI: geo.accuracy is a red herring here, its only used for one type of projector (interpolated) but the default is ray-voxel (Siddon), so its not being used, and it would also not have an influence, I think.
Hum, this seems to be related to numerical precission in the backprojection indeed (not filtering).
I thought it was filtering, because indeed I see no issue when I replace algs.fkd by tigre.Atb.
However, if I reconstruct using an algorithm that essentially only requires forward/backprojection, e.g. SIRT, I do see the effect too.
It is plausible that it is related to #355, which has a relatively easy fix, so I will try it.
It could also be that TIGRE uses floats and not doubles, as you are pushing the magnification by a lot (in real life it would be imposible to have a source so good that 600 mag would not be just a blur). Let me try to see if #355 fixes it....
In particular, this seems to be the solution:
Change the corresponding TIGRE function in the same location by the following version:
https://github.com/IbexInnovations/TIGRE/blob/8d8c8eb4e7be922a69512baa08b4dad5bf516cc3/Common/CUDA/Siddon_projection.cu#L680
I don't' have access to PCs to recompile TIGRE until end of the week, if you want to give this a shot before, just replace the function by the linked one.
Yes, we can't do scans with 600 mag :) but with 100, and there are artifacts already, so it seems to have some relevance for real life scans.
I tried re-installing TIGRE with the changed function in the source files, but it didn't solve the problem.
(I hope it recompiled properly for re-installation, I did pip install . as in the installation instructions and the verbose output of pip said "Successfully built pytigre", "Successfully uninstalled pytigre-2.4.0", "Successfully installed pytigre-2.4.0", so it looked promising).
It would surprise me if it was a limitation of 32-bit float, since it's still floating-point numbers and i'd expect the magnification to be "taken care of by the exponent" and not limited by the precision of the fraction of the floating-point representation, but that's just a guess. The commercial software we use that doen't produce the artifacts also stores at least the reconstruction volume with 32-bit float (we can tell by the total memory usage), but then again the ray-voxel-intersection might be implemented differently...
@simonCtBern thanks for checking! Indeed, preccissionwise I was talking about the ray-tracing precision, not image values, those should be fine.
This will certainly need some digging in the CUDA code. I think it has some high priority for me to fix it, but its quite a hard problem to tackle, and I am quite busy in the next 4 weeks, just to let you know. If yo feel brave enough to dig on the projection codes (Siddon_projection.cu and voxel_backprojection.cu), feel free to tackle it, but these are quite non-trivial files to be playing with.
Will let you know as soon as I can fix this.