Discrepancy (in Cone-beam geometry) between applying filterProjections then backproject and directly using FBP

[hhuang91]

Hi,

First of all, thank you for the amazing work!

I'm encountering some discrepancy issue for conebeam geometry. Basically, I tried to filter the projection first and then backproject the filtered projection. However, the results are different from directly calling FBP function (the latter seems to give me correct attenuation coefficient value in reconstruction, but the former was about ~0.1x the correct value). I tried to find what FBP does uniquely, and thought that maybe weightbackprojection was the cause, but even applying weightedbackprojection to filtered projection, the issue was still there. I would really appreciate some help on this.

Here is some part of the code that I use for filter and then backproject:

  projector = leaptorch.Projector(forward_project=False,use_static=False,use_gpu=True,gpu_device=torch.device(device),batch_size=1)
  projector.set_conebeam(numAngles, 
                              numRows, 
                              numCols, 
                              pixelHeight, 
                              pixelWidth, 
                              centerRow, 
                              centerCol, 
                              phis, 
                              sod, 
                              sdd)
  projector.set_volume(numX, numY, numZ, voxelWidth, voxelHeight)
  projector.allocate_batch_data()
  projF = projector.leapct.filterProjections(proj.clone())
  vol = projector(projF)

when trying to use the weightedbackprojector, I replaced the last line of code with this vol = projector.leapct.weightedBackproject(projF[0], projector.vol_data[0]) But it did not help

[hhuang91]

oh, the reason why I need to do filteration and backprojection separately is that my current projects work with partial angle reconstruction.

[kylechampley]

Hello. Yes, I put in the filterProjections command in there so that people could run the filtering and backprojection steps separately.

You are correct in that you should be using filterProjections followed by weightedBackprojection in general, but for axial cone-beam data there is no difference between backproject and weightedBackproject. Anyway, in my test case I do get the correct reconstruction doing these separately. Would you mind running the following command and sending me the output: projector.leapct.print_parameters()

[kylechampley]

I found and fixed a bug associated with the weightedBackproject function. This may resolve your problem. The change is in the champley_dev branch here: https://github.com/LLNL/LEAP/tree/champley_dev

This will get merged into the main branch for the next release.

[hhuang91]

Hi,

My apologies for lack of communication. I have been occupied by recent deadlines. I tried to use the projector.leapct.print_parameters(), but nothing comes up (I'm using python 3.9 with Spyder 5.1.5). However, I was able to use save_parameters() function and here is what's in the saved file

# CT volume parameters
numX = 256
numY = 256
numZ = 256
voxelWidth = 1.000000e+00
voxelHeight = 1.000000e+00
offsetX = 0.000000e+00
offsetY = 0.000000e+00
offsetZ = 0.000000e+00

# CT geometry parameters
geometry = cone
numAngles = 304
numRows = 660
numCols = 864
pixelHeight = 6.160000e-01
pixelWidth = 6.160000e-01
centerRow = 3.295000e+02
centerCol = 4.315000e+02
angularRange = 2.598548e+02
sod = 7.850000e+02
sdd = 1.200000e+03
tau = 0.000000e+00
helicalPitch = 0.000000e+00

I will try out the dev branch, see it gets resolved, and report back.

Thank you so much for the help!

[kylechampley]

It seems to work out for me when using the latest branch. I got a perfect match (difference image was exactly zero) when I run the text below. FYI, you will get better computational performance if you use smaller voxels. Here is my test script:

import sys import os import time import numpy as np from leapctype import * leapct = tomographicModels()

leapct.set_conebeam(304, 660, 864, 6.160000e-01, 6.160000e-01, 3.295000e+02, 4.315000e+02, leapct.setAngleArray(304, 2.598548e+02), 7.850000e+02, 1.200000e+03) leapct.set_volume(256, 256, 256, 1.0, 1.0)

g = leapct.allocateProjections() f = leapct.allocateVolume() leapct.set_FORBILD(f,True)

device_name = "cuda:" + str(leapct.get_gpu()) device = torch.device(device_name) g = torch.from_numpy(g).to(device) f = torch.from_numpy(f).to(device) f_2 = f.clone() f_2[:] = 0.0

"Simulate" projection data

leapct.project(g,f)

leapct.FBP(g,f)

q=leapct.filterProjections(g) leapct.weightedBackproject(q,f_2)

diff = f-f_2 print(torch.max(torch.abs(diff)))

[kylechampley]

Did this issue get resolved?

[hhuang91]

Hi,

Yes. The dev branch has resolved the issue. Thank you!

Apologies for delayed update.