Fix memory init in pocs_tv and aw_pocs_tv

[tsadakane]

Summary

Fixed a bug that the pocs algorithms give different results for each run.

Detail

Running the program below...

#%%Initialize
import tigre
import numpy as np
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise
import tigre.algorithms as algs
from matplotlib import pyplot as plt
from tigre.utilities.im3Dnorm import im3DNORM
import os

#%% Geometry
geo = tigre.geometry_default(high_resolution=False)

#%% Load data and generate projections
# define angles
angles = np.linspace(0, 2 * np.pi, 100)
if os.path.isfile("noise_projections.npy"):
    print("Loading noise_projections.npy")
    noise_projections = np.load("noise_projections.npy")
else:
    # Load thorax phatom data
    head = sample_loader.load_head_phantom(geo.nVoxel)
    # generate projections
    projections = tigre.Ax(head, geo, angles)
    # add noise
    noise_projections = CTnoise.add(projections, Poisson=1e5, Gaussian=np.array([0, 10]))
    np.save("noise_projections.npy", noise_projections)

epsilon = (
    im3DNORM(tigre.Ax(algs.fdk(noise_projections, geo, angles), geo, angles) - noise_projections, 2)
    * 0.15
)
alpha = 0.002
niter = 3
ng = 3#5

lmbda = 1
lambdared = 0.9999  # you generally want 1

#   'alpha_red':   Defines the reduction rate of the TV hyperparameter
alpha_red = 0.95

#   'Ratio':       The maximum allowed image/TV update ration. If the TV
#                  update changes the image more than this, the parameter
#                  will be reduced.default is 0.95
ratio = 0.94

#   'Verbose'      1 or 0. Default is 1. Gives information about the
#                  progress of the algorithm.

verb = True

imgASDPOCS = algs.asd_pocs(
    noise_projections,
    geo,
    angles,
    niter,  # these are very important
    tviter=ng,
    maxl2err=epsilon,
    alpha=alpha,  # less important.
    lmbda=lmbda,
    lmbda_red=lambdared,
    rmax=ratio,
    verbose=verb,
)

imgOSASDPOCS = algs.os_asd_pocs(
    noise_projections,
    geo,
    angles,
    niter,  # these are very important
    tviter=ng,
    maxl2err=epsilon,
    alpha=alpha,  # less important.
    lmbda=lmbda,
    lmbda_red=lambdared,
    rmax=ratio,
    verbose=verb,
    # OSC params
    blocksize=10,
)

imgAWASDPOCS = algs.awasd_pocs(
    noise_projections,
    geo,
    angles,
    niter,  # these are very important
    tviter=ng,
    maxl2err=epsilon,
    alpha=alpha,  # less important.
    lmbda=lmbda,
    lmbda_red=lambdared,
    rmax=ratio,
    verbose=verb,  # AwASD_POCS params
    delta=np.array([-0.005]),
)

#%% plot results
# np.save("imgASDPOCS.npy"  , imgASDPOCS  *1000)
# np.save("imgOSASDPOCS.npy", imgOSASDPOCS*1000)
# np.save("imgAWASDPOCS.npy", imgAWASDPOCS*1000)

# plot images
tigre.plotimg(
    np.concatenate([imgAWASDPOCS, imgOSASDPOCS, imgASDPOCS], axis=1), dim="z", step=2,
    # savegif="result_base.gif"
    savegif="result_PR.gif"
)

The result of the code is like

POCS Algorithm in progress.
Estimated time until completion : 00:00:23

     Stop criteria met:
     c = -0.06136499
     beta = 0.9996000599960002
     iter = 4

POCS Algorithm in progress.
Estimated time until completion : 00:00:02

     Stop criteria met:
     c = -0.28363463
     beta = 0.9996000599960002
     iter = 4

POCS Algorithm in progress.
Estimated time until completion : 00:00:23

     Stop criteria met:
     c = -0.08796015
     beta = 0.9996000599960002
     iter = 4

It is expected that the values c are the same for each run, but different.

run	asd_pocs	os_asd_pocs	awasd_pocs
1	-0.06136499	-0.28363463	-0.08796015
2	-0.7518723	-0.20649311	-0.092016704
3	-0.24782406	-0.19568634	-0.038304713

It is easily confirmed that it stems from the fact that the value of dg in ASD_POCS.run_main_iter(self) is sometimes different. The value dg is calculated from pocs_tv function in pocs_tv.cu.

I found that some of the memories are not initialized after its allocation, and added some lines as this PR.

After this fix, the values c are the same for each run (for each algorithm) as expected:

run	asd_pocs	os_asd_pocs	awasd_pocs
1	-0.43480712	-0.33600685	-0.13818325
2	-0.43480712	-0.33600685	-0.13818325
3	-0.43480712	-0.33600685	-0.13818325

Environment

Windows 10 CUDA 11.1 GeForce GTX 1060 6GB x1 Python 3.8.10 (anaconda) Visual Studio 2019 16.11.0

Additionally...

Noise disappeared...? It seems to be plausible.

Base

PR