The OSSART_TV algorithm in Python has problems
Closed this issue · 6 comments
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:
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.
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.
this is result of TV
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')
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.
You are right, working on fixing it right now.
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")