General debugging tips

[ajensen1234]

Hi! I'm wondering if there are any general guidelines for debugging that might help to steer toward a good reconstruction! Tangentially, this would be related to the sensitivity of the reconstruction, and how close the geometry needs to be defined to get even a “bad” reconstruction.

For context, I have a collection of cone-beam fluoro images of a hip phantom captured from an OEC-9800:

In order to determine the geo parameters, we used motion capture rigidly affixed to both the phantom and the OEC C-ARM to track motion.

Here's a rough view I've generated of the imaging geometry during a “sweep” of the c-arm:

And so, back to the original question. How sensitive would the reconstruction (say, CGLS) be of minor perturbations in the geometry definitions? If I knew the source to detector (DSD) and source to origin (DSO), would just using those be sufficient for generating a “bad” reconstruction? Using all combinations of geometry offsets, center of rotation offsets, etc, the only results I've achieved look like the images attached below:

Is there a “best” way to go about debugging these types of issues? Should I even expect that getting things “close enough” should get even a bad reconstruction? Thanks!

[AnderBiguri]

Reconstructions are extremely sensitive to geometric inaccuracies, but generally a DSD/DSO mismatch is the one that is less important.

How to debug? ah, I don't know. I have build enough heuristic from years at looking at images to start slowly thinking about what things could be, but I have not found a fixed way.

Few things that my brain things when seeing your images:

1- They've gone super bad. Negative values in the image are bad!
2- Above is quite possibly enhanced/caused by all those zeroes in your projection data. Zero is a valid number. FDK does not suffer from this, because it only takes the high frequency values of data, but any iterative algorithm will enforce detector values. If they see zeros, it can only be caused by negative values in the image (to compensate the other views where the detector has a value), and a negative value means a source of X-rays, which a patient should not be. Basically you need to figure out a way to either inpaint or crop all that area if you want to use iteratives.
3- maybe your angles are negative sign. You get concave artifacts with edges pointing outwards, this often suggest some sign change (or projection flip) required.
4- Your projections are not log transformed. TIGRE (and any recon) requires X-ray absorbtion (i.e. after Beer-law) to reconstruct, not X-ray transmission (i.e. raw measurement, light amount), otherwise its not a linear problem. The real data demo has some code to do so. You want "air" to be black, "bone" to be white.

I don't fully understand your geometry diagram, so not sure what else to suggest, but if you explain how it differs from a standard circular trajectory, I may have suggestions on changing the TIGRE geo.

[ajensen1234]

Hi Ander,
Thanks for the reply!
I think I was able to figure out some of the dimensional transforms and shifts going on (I'm guessing having the dimensions end up "unitless" makes getting things into CUDA kernels a lot easier!)

One thing I was a little bit confused about was defining the location of the source after a rotation.
In the order of the source getting interacted with in the fuction:

1. Create source object, and move it along the x-axis to DSO (lines 633-636)

    Point3D S;
    S.x=geo.DSO[i];
    S.y=0;
    S.z=0;

2. Rotate the object by $\theta$ (line 694). For now, we'll assume just one rotational degree of freedom, so it ends up just being a z rotation.

    eulerZYZ(geo,&S);

3. Negative source offset (same relative transform as positive image offset) (line 702)

S.x=S.x-geo.offOrigX[i]; 
S.y=S.y-geo.offOrigY[i];  
S.z=S.z-geo.offOrigZ[i];

4. Adjusting pixel locations to get the center of the image as the origin and scaling units (708, 714)

S.x       =   S.x+geo.sVoxelX/2-geo.dVoxelX/2;           
S.y       =   S.y+geo.sVoxelY/2-geo.dVoxelY/2;               
S.z       =   S.z  +geo.sVoxelZ/2-geo.dVoxelZ/2;
S.x       =    S.x/geo.dVoxelX;           
S.y       =    S.y/geo.dVoxelY;            
S.z       =    S.z/geo.dVoxelZ;

5. Adjusting for center of rotation (lines 720-726). You do this by rotating a y-vector about the same rotation $\theta$

    float CORx, CORy;
    CORx=-geo.COR[i]*sin(geo.alpha)/geo.dVoxelX;
    CORy= geo.COR[i]*cos(geo.alpha)/geo.dVoxelY;
    Pfinal.x+=CORx;   Pfinal.y+=CORy;
    Pfinalu0.x+=CORx;   Pfinalu0.y+=CORy;
    Pfinalv0.x+=CORx;   Pfinalv0.y+=CORy;
    S.x+=CORx; S.y+=CORy;

So, if we remove the dimensional scaling, and the relative transforms that are occuring to account for the image offset, the transforms that soley describe the location of the detector is broken down into the $\theta$ rotation and the center of rotation offset.

However, I think there's either (1) and error in the way that the final location of the source is being determined, or (2) I have a misunderstanding of what the center of rotation offset means.

Regarding #2, my understanding of the center of rotation offset is that there is a point that does not lay on the TIGRE x-axis about which the source-detector unit is spinning.

Here's a diagram:

For this understanding of the center of rotation, the equivalent transformation is:

1. Move the source by COR.
2. Apply the rotation
3. Move the source by -COR along the rotated axis.

Here's a diagram:

And so, if my understanding of both a rotation about COR is correct, and I've read the computeDeltas correctly, then I think there's an issue. Mainly, that the application of this rotation is just:

1. Rotate $(DSO, 0, 0)$ about $\theta_y$ (line 694)
2. Apply rotated COR vector (lines 720-726)

Here's a diagram:

It's distinctly possible that I'm missing something here, and perhaps some of this is getting handled in the dimensional adjustments that I didn't quite pick up. I do think, if your COR was sufficiently small, then these effects would probably not be noticable. But, unfortunately for us, we had a large COR, so my diagrams just assumed that it was big. The math would remain the same, though.

Please do let me know what you think; I'm sorry the message is so long, but I wanted to be thorough.

[AnderBiguri]

I'll try to check all this logic when I have a little bit of time. I tend to just lose track of what the exact math is here, its clearly non-trivial and not super easy to follow. However, its extremely unlikely there is a mistake, as I have verified my operators with many many other recosntruction tools over the years, and I can get the same results given the same inputs.

[ajensen1234]

That sounds good. Thanks ahead of time for taking the time to look it over!

[ajensen1234]

Hi Ander,
I just found the ASTRA toolbox (https://github.com/astra-toolbox/astra-toolbox), which seems to support some arbitrary geometries, at the cost of less flexibility in the algorithm construction. Having done so, I found the following issue #515 and commit 66cfc43.

I think this might be exactly the solution that I was looking for, as it was a bit more straightforward to directly calculate the transformation matrices of our system.

I will let you know if this works. If it does, I can close this discussion, and hopefully all the dialogue will be useful for anyone else in the future trying to do fun things with CT reconstruction!

[AnderBiguri]

Hi @ajensen1234 . ASTRA==TIGRE for flexibility of geometries, the only difference really is how they are defined :)

I had completely forgotten about that PR. Indeed, that seems to be the thing you need. Do let me know if it works, I never tested it myself!