dyang37
commented
2 years ago Ok I think I found out the cause of the nan problem. It is because of the sinogram weight. The way demo worked is that we added transmission noise to the synthetic sinogram, and perform recon with transmission weight type calculated from noisy sinogram:
sino = mbircone.cone3D.project(phantom, angles,
num_det_rows, num_det_channels,
dist_source_detector, magnification)
sino_weights = mbircone.cone3D.calc_weights(sino, weight_type='transmission')
sino_noise_sigma = 0.01 # transmission noise level
noise = sino_noise_sigma * 1. / np.sqrt(sino_weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))
sino_noisy = sino + noise
recon_qGGMRF = mbircone.cone3D.recon(sino_noisy, angles, dist_source_detector, magnification,
sharpness=sharpness, weight_type='unweighted',
verbose=1)It appears that the transmission weight calculated from noisy sinogram has some extreme values. The maximum values of the calculated weight is 7.845908754559244e+96, and that caused the nan problem we saw above. (Before we changed the detector pixel pitch, the max weights value was 0.9602327231700871)
As a sanity check, I simply clipped the maximum value of weights to be 10.0, and the nan problem goes away:
************************** Iteration 0 (max. 20) **************************
* Cost = 2.4964202496e+11
* Rel. Update = 0.0000000000e+00 % (threshold = 1.9999999553e-02 %)
* RWFE = ||e||_W/||y||_W = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* RUFE = ||e|| / ||y|| = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
* 1/M ||e||^2_W = 3.5201702118e+01 = 1/0.0284077171
* weightScaler_value = 3.5534198105e-05 = 1/28141.9042968750
* ----------------------------------------------------------------------------
* voxelsPerSecond = 0.0000000000e+00
* time icd update = -3.8534309715e-02 s
* ratioUpdated = 0.0000000000e+00 %
* totalEquits = 0.0000000000e+00
******************************************************************************
************************** Iteration 1 (max. 20) **************************
* Cost = 2.4892886221e+11
* Rel. Update = 1.0000000000e+02 % (threshold = 1.9999999553e-02 %)
* RWFE = ||e||_W/||y||_W = 9.9853851318e+01 % (threshold = 0.0000000000e+00 %)
* RUFE = ||e|| / ||y|| = 7.2232589722e+01 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
* 1/M ||e||^2_W = 3.5098880768e+01 = 1/0.0284909364
* weightScaler_value = 3.5534198105e-05 = 1/28141.9042968750
* ----------------------------------------------------------------------------
* voxelsPerSecond = 6.9297718750e+04
* time icd update = 9.8886108398e+00 s
* ratioUpdated = 1.0000000000e+02 %
* totalEquits = 1.0000000000e+00
*****************************************************************************Potential approach to fix the issue:
Previously the transmission noise level is hard coded to be sino_noise_sigma = 0.01, and transmission noise is calculated as
noise = sino_noise_sigma * 1. / np.sqrt(sino_weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))This noise should be scaled w.r.t. the detector pixel pitch as well, since by changing detector pixel pitch from 0.25 to 1, the sinogram value is also scaled up:
Below is the clean synthetic sinogram using pixel pitch = 0.25
Below is the clean synthetic sinogram u
sing pixel pitch = 1.0
As a result, the sino_noise_sigma value should also be scaled down by 4x, so that the sinogram noise level noise = sino_noise_sigma * 1. / np.sqrt(sino_weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels)) stays the same.
The demo script
demo_mace3D.pydoes not run correctly. The QGGMRF part appears to generate inf values. This needs to be tracked down and fixed.