cabouman / mbircone

BSD 3-Clause "New" or "Revised" License
11 stars 9 forks source link

Additional nanfixs #43

Closed dyang37 closed 3 years ago

dyang37 commented 3 years ago

This PR contains the following two updates:

  1. In cone3D.py add corner case protection of negative weights entry.
  2. In src/recon3DCone.c, Fix incorrect reported ||e||_W/||y||_W = -nan when ||y||_W is 0. This is the root cause of the -nan issue for the case of weights=0. In this case ||e||_W is reported instead.

Tested both 3D shepp Logan and mace3D demos, both work. For the second fix, tested 3D shepp Logan with 0 weights. This time the reported e2 over y2 is correct:

************************** Iteration 0  (max. 20) **************************
*  Cost                   = 0.0000000000e+00
*  Rel. Update            = 0.0000000000e+00 % (threshold = 2.0000000000e-02 %)
*  RWFE = ||e||_W/||y||_W = 0.0000000000e+00 % (threshold = 0.0000000000e+00 %)
*  RUFE = ||e|| / ||y||   = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
*  1/M ||e||^2_W          = 0.0000000000e+00 = 1/inf
*  weightScaler_value     = 1.0000000000e+00 = 1/1.0000000000
* ----------------------------------------------------------------------------
*  voxelsPerSecond        = 0.0000000000e+00
*  time icd update        = 1.9073486328e-06 s
*  ratioUpdated           = 0.0000000000e+00 %
*  totalEquits            = 0.0000000000e+00
******************************************************************************

************************** Iteration 1  (max. 20) **************************
*  Cost                   = 0.0000000000e+00
*  Rel. Update            = 0.0000000000e+00 % (threshold = 2.0000000000e-02 %)
*  RWFE = ||e||_W/||y||_W = 0.0000000000e+00 % (threshold = 0.0000000000e+00 %)
*  RUFE = ||e|| / ||y||   = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
*  1/M ||e||^2_W          = 0.0000000000e+00 = 1/inf
*  weightScaler_value     = 1.0000000000e+00 = 1/1.0000000000
* ----------------------------------------------------------------------------
*  voxelsPerSecond        = 8.3719536628e+04
*  time icd update        = 1.6813638210e+00 s
*  ratioUpdated           = 1.0000000000e+02 %
*  totalEquits            = 1.0000000000e+00
******************************************************************************

[ticToc] MBIR3DCone = 1.821220e+00 s
DamonLee5 commented 3 years ago

I ran the demo and did not have ||e||_W/||y||_W = -nan after this fix.

************************** Iteration 0  (max. 20) **************************
*  Cost                   = 0.0000000000e+00
*  Rel. Update            = 0.0000000000e+00 % (threshold = 2.0000000000e-02 %)
*  RWFE = ||e||_W/||y||_W = 0.0000000000e+00 % (threshold = 0.0000000000e+00 %)
*  RUFE = ||e|| / ||y||   = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
*  1/M ||e||^2_W          = 0.0000000000e+00 = 1/inf       
*  weightScaler_value     = 1.0000000000e+00 = 1/1.0000000000
* ----------------------------------------------------------------------------
*  voxelsPerSecond        = 0.0000000000e+00 
*  time icd update        = 2.1457672119e-06 s
*  ratioUpdated           = 0.0000000000e+00 %
*  totalEquits            = 0.0000000000e+00 
******************************************************************************

************************** Iteration 1  (max. 20) **************************
*  Cost                   = 0.0000000000e+00
*  Rel. Update            = 0.0000000000e+00 % (threshold = 2.0000000000e-02 %)
*  RWFE = ||e||_W/||y||_W = 0.0000000000e+00 % (threshold = 0.0000000000e+00 %)
*  RUFE = ||e|| / ||y||   = 1.0000000000e+02 % (threshold = 0.0000000000e+00 %)
* ----------------------------------------------------------------------------
*  1/M ||e||^2_W          = 0.0000000000e+00 = 1/inf       
*  weightScaler_value     = 1.0000000000e+00 = 1/1.0000000000
* ----------------------------------------------------------------------------
*  voxelsPerSecond        = 6.4630340511e+04 
*  time icd update        = 2.1779708862e+00 s
*  ratioUpdated           = 1.0000000000e+02 %
*  totalEquits            = 1.0000000000e+00 
******************************************************************************

[ticToc] MBIR3DCone = 2.299892e+00 s
recon shape =  (103, 59, 59)
press Enter