CERN / TIGRE

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

Use original geometry for FDK redundancy weighting #416

Closed twhitbread closed 1 year ago

twhitbread commented 1 year ago

Following https://github.com/CERN/TIGRE/pull/415 I noticed another bug in FDK offset weighting. It may have been introduced by me in https://github.com/CERN/TIGRE/pull/415, or in the original PR https://github.com/CERN/TIGRE/pull/376.

The redundancy_weighting function was using the new zgeo structure that comes out of zeropadding, which means the detector offset is zero (so the weighting is turned off automatically) and theta was being calculated incorrectly anyway. I fixed this by passing in the original geo structure as an additional input for FDK only; for all other algorithms there is no zero-padding as far as I'm aware, so the correct geo structure is used.

Also fixed a weight bug in MLEM.

AnderBiguri commented 1 year ago

Thanks @twhitbread ! Is there a chance you can give me a geometry to reproduce the bug? The solution is correct but I may try to find a cleaner way of obtaining this behavior than passing 2 geos.

twhitbread commented 1 year ago

Sure:

angles = linspace(0, 2*pi, 101); angles(end) = []; geo = defaultGeometry('nVoxel', [128; 128; 128]); geo.rotDetector = [0; 0; 0]; geo.offDetector(1) = 100;

twhitbread commented 1 year ago

I brought in recent changes from master and MLEM is failing. Will look for a solution and push to this branch if I can find one.

twhitbread commented 1 year ago

The bug is in the computation of W. I think this is a mis-named variable: it's usually a forward-projection weight, but in MLEM it was computed by doing a back-projection (so more like V weight?). In the recent PR the code was replaced with the new computeW function, which does a forward-projection, so W is not the correct size, and we get an error on line 65 of MLEM.

Do you have a fix suggestion? Revert back to the old way? Or use the computeV function?

AnderBiguri commented 1 year ago

@twhitbread I think that should be compteV, yes...

twhitbread commented 1 year ago

@AnderBiguri computeV requires variables like alphablocks which are not currently available in MLEM. I've just reverted to the old calculation for now, and changed the variable name to V.

AnderBiguri commented 1 year ago

@twhitbread sure!

For reference (for me or you to fix), the computeV for MLEM is the same trace as SIRT, the same input params.

twhitbread commented 1 year ago

@AnderBiguri using the line from SIRT changes the scale of the reconstruction quite a bit:

image

Is this a problem or are you happy with this?

AnderBiguri commented 1 year ago

@twhitbread hum, I must have screwed this in some other way.... Certainly should not happen. I will investigate