Helical CT reconstruction does not work in Python

[resonance20]

Expected Behavior

Tried to implement Demo 13 from the MATLAB Code within Python

Actual Behavior

A ValueError is raised while computing V. The entire error message:

ValueError                                Traceback (most recent call last)
c:\Users\z003zv1a\Desktop\Code\denoise_dl\backproject.py in 
     64 #print(help(algs.fdk))
     65 
---> 66 imgfdk1=algs.sirt(projections,geo,angles,niter=20)
     67 
     68 # The look quite similar:

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pytigre-0.1.8-py3.7-win-amd64.egg\tigre\algorithms\iterative_recon_alg.py in iterativereconalg(proj, geo, angles, niter, **kwargs)
    406 
    407     def iterativereconalg(proj, geo, angles, niter, **kwargs):
--> 408         alg = IterativeReconAlg(proj, geo, angles, niter, **kwargs)
    409         if name is not None:
    410             alg.name = name

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pytigre-0.1.8-py3.7-win-amd64.egg\tigre\algorithms\art_family_algorithms.py in __init__(self, proj, geo, angles, niter, **kwargs)
     37     def __init__(self, proj, geo, angles, niter, **kwargs):
     38         kwargs.update(dict(blocksize=angles.shape[0]))
---> 39         IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs)
     40 
     41 

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pytigre-0.1.8-py3.7-win-amd64.egg\tigre\algorithms\iterative_recon_alg.py in __init__(self, proj, geo, angles, niter, **kwargs)
    169             self.set_w()
    170         if not hasattr(self, 'V'):
--> 171             self.set_v()
    172         if not hasattr(self, 'res'):
    173             self.set_res()

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pytigre-0.1.8-py3.7-win-amd64.egg\tigre\algorithms\iterative_recon_alg.py in set_v(self)
    208             step = -geo.dVoxel[1]
    209 
--> 210             xv = np.arange(start, stop + step, step)
    211 
    212             start = geo.sVoxel[2] / 2 - geo.dVoxel[2] / 2 + geo.offOrigin[2]

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Code to reproduce the problem (If applicable)

import tigre
import numpy as np

class Geometry(tigre.utilities.geometry.Geometry):

    def __init__(self):
        # VARIABLE                                          DESCRIPTION                    UNITS
        # -------------------------------------------------------------------------------------
        self.DSD = 1536                                     # Distance Source Detector      (mm)
        self.DSO = 1000                                     # Distance Source Origin        (mm)
        # Detector parameters
        self.nDetector = np.array((512, 512))               # number of pixels              (px)
        self.dDetector = np.array((0.8, 0.8))               # size of each pixel            (mm)
        self.sDetector = self.nDetector * self.dDetector    # total size of the detector    (mm)
        self.rotDetector = np.array((0, 0, 0)) 
        # Image parameters
        self.nVoxel = np.array((512, 512, 120))             # number of voxels              (vx)
        self.sVoxel = np.array((256, 256, 120))             # total size of the image       (mm)
        self.dVoxel = self.sVoxel/self.nVoxel               # size of each voxel            (mm)
        # Offsets
        self.offDetector = np.array((0, 0))                 # Offset of Detector            (mm)

        # Auxiliary
        self.accuracy = 0.5                                 # Accuracy of FWD proj          (vx/sample)
        # Mode
        self.mode = 'cone'                                  # parallel, cone                ...

geo=Geometry()
from tigre.utilities.Ax import Ax
from tigre.demos.Test_data import data_loader

# define angles
angles=np.linspace(0, 2*np.pi, num=180, dtype=np.float32)
angles = np.stack([angles, angles, angles], axis = 0).transpose()

#Define offsets
offsets_dist = np.linspace(-256, 256, num=180)
offsets_zero = np.zeros(offsets_dist.shape)
geo.offOrigin = np.stack([offsets_zero, offsets_zero, offsets_dist], axis = 0).transpose()

# load head phantom data
head=data_loader.load_head_phantom(number_of_voxels=geo.nVoxel)
# generate projections
projections=Ax(head,geo,angles,'interpolated')

import tigre.algorithms as algs
#print(help(algs.fdk))

imgfdk1=algs.sirt(projections,geo,angles,niter=20)

# The look quite similar:
tigre.plotimg(imgfdk1,slice=32,dim='x')

Specifications

• Python version: 3.7 (Anaconda)
• OS: Windows 10
• CUDA version: 10.1

[AnderBiguri]

The problem is that compute_V does not handle well arbitrary offOrigin, I think.

For helical, the easy fix (and as far as I can think, with correct values) is to change in the following code:

https://github.com/CERN/TIGRE/blob/master/Python/tigre/algorithms/iterative_recon_alg.py#L206-L214

geo.offOrigin[1] to geo.offOrigin[1,0] and geo.offOrigin[2] to geo.offOrigin[2,0]

The "hard" but proper fix is to compute all inside the for loop, slower, but more flexible to arbitrary offsets:

https://github.com/CERN/TIGRE/blob/master/Python/tigre/algorithms/iterative_recon_alg.py#L222

---

You may have this problem in other places though, as I am surprised this was not tested before... ah apologies, the python version is quite a mess and currently I dont have the time it requires to fix it.

[genusn]

@AnderBiguri The commit db8dfd2 seems to be incompatible with example.py. It stopped with the following error massage:

  File "....\Python\tigre\algorithms\iterative_recon_alg.py", line 206, in set_v
    start = geo.sVoxel[1] / 2 - geo.dVoxel[1] / 2 + geo.offOrigin[1,0]
IndexError: too many indices for array

[AnderBiguri]

Ah, my bad. I did not test enough. I think if eaxh algorithm would check_geometry to start with then it would be solved, though this was happening. I'll investigate more how to avoid this

[AnderBiguri]

@genusn fixed properly this time.

[genusn]

@AnderBiguri Thank you. I confirmed that example.py is executable without error.