The OSSART_TV algorithm in Python has problems

[hy010227]

The results of reconstruction using OSSART and OSSART_TV are indistinguishable, and adjusting the TV parameters also produces no changes.

Moreover, I am unable to view the minTV function here; it tells me "cannot find declaration to go to."

import tigre
import numpy as np
import tigre.algorithms as algs
import os

geo = tigre.geometry()

# Distances
geo.DSD = 1300 # Distance Source Detector      (mm)
geo.DSO = 1000  # Distance Source Origin        (mm)
# Detector parameters
geo.nDetector = np.array([822, 824])  # number of pixels              (px)
geo.dDetector = np.array([0.192, 0.192])  # size of each pixel            (mm)
geo.sDetector = geo.nDetector * geo.dDetector  # total size of the detector    (mm)
# Image parameters
geo.nVoxel = np.array([512, 512, 512])  # number of voxels              (vx)
geo.sVoxel = np.array([512, 512, 512]) * np.array([0.1,0.1,0.1])  # total size of the image       (mm)
geo.dVoxel = np.array([0.1,0.1,0.1])  # size of each voxel            (mm)
# Offsets
geo.offOrigin = np.array([0,0, 0])  # Offset of image from origin   (mm)
geo.offDetector = np.array([0, 0])  # Offset of Detector            (mm)
geo.mode = "cone"
geo.COR = 0
geo.accuracy = 0.5
geo.rotDetector = np.array([0, 0, 0])
print(geo)

#%% Load data and generate projections

#angles
clockwise = -np.pi / 2
angles = np.linspace(0, 2 * np.pi, 720, endpoint=False)+ clockwise
angles1 = angles[0:30]
angles2 = angles[150:210]
angles3 = angles[690:720]
angles = np.concatenate((angles1,angles2, angles3), axis=0)

#projections
proj = np.fromfile('./proj_simu.raw',dtype=np.float32).reshape(720,822,824)
proj1 = proj[0:30, :, :]
proj2 = proj[150:210, :, :]
proj3 = proj[690:720, :, :]
proj = np.concatenate((proj1, proj2, proj3),axis=0)

# TV
niter = 20
imgOSSART_TV = algs.ossart_tv(proj, geo, angles, niter,tvlambda=100,tviter=20)
imgOSSART_TV.astype(np.float32).tofile(os.path.join('./TV/','TV_simuiter20_-30-30_150-210.raw'))

I adjusted the TV parameters to a large value to observe whether TV has any effect

Specifications

• MATLAB/python version:
• OS:
• CUDA version:

[AnderBiguri]

How much have you adjusted the parameters?

minTV is a cuda code, so you can not open it in python. If you show a piece of code with the TIGRE data, I can help you debug your script, but I can't with your data.

[hy010227]

The following is the reconstruction code. I have made various changes to tvlambda and tviter, but I found no improvement in the images compared to OSSART. I have tried the parameters 0.2, 50, and 100 for tvlambda, but there was no difference.

import tigre

import numpy as np

import tigre.algorithms as algs

import os

geo = tigre.geometry()

# Distances

geo.DSD = 1300 # Distance Source Detector      (mm)

geo.DSO = 1000  # Distance Source Origin        (mm)

# Detector parameters

geo.nDetector = np.array([822, 824])  # number of pixels              (px)

geo.dDetector = np.array([0.192, 0.192])  # size of each pixel            (mm)

geo.sDetector = geo.nDetector * geo.dDetector  # total size of the detector    (mm)

# Image parameters

geo.nVoxel = np.array([512, 512, 512])  # number of voxels              (vx)

geo.sVoxel = np.array([512, 512, 512]) * np.array([0.1,0.1,0.1])  # total size of the image       (mm)

geo.dVoxel = np.array([0.1,0.1,0.1])  # size of each voxel            (mm)

# Offsets

geo.offOrigin = np.array([0,0, 0])  # Offset of image from origin   (mm)

geo.offDetector = np.array([0, 0])  # Offset of Detector            (mm)

geo.mode = "cone"

geo.COR = 0

geo.accuracy = 0.5

geo.rotDetector = np.array([0, 0, 0])

print(geo)

#%% Load data and generate projections

#angles

clockwise = -np.pi / 2

angles = np.linspace(0, 2 * np.pi, 720, endpoint=False)+ clockwise

angles1 = angles[0:30]

angles2 = angles[150:210]

angles3 = angles[690:720]

angles = np.concatenate((angles1,angles2, angles3), axis=0)

#projections

proj = np.fromfile('./proj_simu.raw',dtype=np.float32).reshape(720,822,824)

proj1 = proj[0:30, :, :]

proj2 = proj[150:210, :, :]

proj3 = proj[690:720, :, :]

proj = np.concatenate((proj1, proj2, proj3),axis=0)

# SART

niter = 20

imgOSSART_TV = algs.ossart_tv(proj, geo, angles, niter,tvlambda=100,tviter=20)

imgOSSART_TV.astype(np.float32).tofile(os.path.join('./TV/','TV_simuiter20_-30-30_150-210.raw'))

this is result of SART.

![Uploading image.png…]()

this is result of TV

[hy010227]

I have now conducted experiments using the built-in data from TIGRE, and I found that there is no difference between the reconstruction results of SART and TV. I subtracted the reconstructed images of the two algorithms and found no error. I'm not sure what could be causing this. The following is the code I used.

import tigre
import numpy as np
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise
import tigre.algorithms as algs

geo = tigre.geometry_default(high_resolution=True)

#%% Load data and generate projections
# define angles
angles = np.linspace(0, 2 * np.pi, 30)
# Load thorax phantom 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]))
#fdk
imgFDK = algs.fdk(noise_projections, geo, angles)
imgFDK = imgFDK.astype(np.float32).tofile('FDK_head.raw')

#SART
niter = 50
imgOSSART = algs.ossart(noise_projections, geo, angles, niter)
imgOSSART.astype(np.float32).tofile('sart_head.raw')

#TV
niter = 50
imgOSSART_TV = algs.ossart_tv(noise_projections, geo, angles, niter)
imgOSSART_TV.astype(np.float32).tofile('sart_tv_head.raw')

[hy010227]

hello, I think the issue may be with the OSSART_TV algorithm. I just experimented using the SART and SART_TV algorithms, and I found that TV is effective.

[AnderBiguri]

You are right, working on fixing it right now.

[hy010227]

Hello, I modified the OSSART_TV algorithm like this, and now the TV is effective.

class OSSART_TV(IterativeReconAlg):
doc = ( "OSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n" "Simultaneous Algebraic Reconstruction Technique with TV regularization algorithm\n" "OSSART_TV(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20,TVLAMBDA=50,TVITER=50) \n" "solves the reconstruction problem using the projection data PROJ taken\n" "over ALPHA angles, corresponding to the geometry described in GEO,\n" "using NITER iterations.\n" ) + IterativeReconAlg.doc

def __init__(self, proj, geo, angles, niter, **kwargs):

    self.blocksize = 20 if 'blocksize' not in kwargs else kwargs['blocksize']
    self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda']
    self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter']
    # these two settings work well for nVoxel=[254,254,254]

    IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs)

"""HXY_modified"""

# Override
def run_main_iter(self):
    """
    Goes through the main iteration for the given configuration.
    :return: None
    """
    Quameasopts = self.Quameasopts

    for i in range(self.niter):

        res_prev = None
        if Quameasopts is not None:
            res_prev = copy.deepcopy(self.res)
        if self.verbose:
            self._estimate_time_until_completion(i)

        getattr(self, self.dataminimizing)()
        # print("run_main_iter: gpuids = {}", self.gpuids)
        self.res = im3ddenoise(self.res, self.tviter, self.tvlambda, self.gpuids)
        if Quameasopts is not None:
            self.error_measurement(res_prev, i)

ossart_tv = decorator(OSSART_TV, name="ossart_tv")