Get "ValueError: operands could not be broadcast together with shapes " bug when using fbp after ossart_tv

[GreameLee]

when I run:

limited_angle = np.linspace(0, np.pi/2, 180)
limited_projection = tigre.Ax(img, geo, limited_angle)
print(limited_projection.shape)
img_tv = ossart_tv(limited_projection,geo,limited_angle,5)
img_tv.shape
img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)
img_fbp.shape

it get the bug,
when I run:

limited_angle = np.linspace(0, np.pi/2, 180)
limited_projection = tigre.Ax(img, geo, limited_angle)
print(limited_projection.shape)
img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)
img_fbp.shape
img_tv = ossart_tv(limited_projection,geo,limited_angle,5)
img_tv.shape

it goes well

[AnderBiguri]

Which line is the bug in?

[GreameLee]

img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)

[AnderBiguri]

@GreameLee I am quite sure the error message is much longer, and that it has more information :) Can you paste the entire error trace?

[GreameLee]

when I run

limited_angle = np.linspace(0, np.pi/2, 180)
limited_projection = tigre.Ax(img, geo, limited_angle)
print(limited_projection.shape)
img_tv = ossart_tv(limited_projection,geo,limited_angle,5)
img_tv.shape
img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)
img_fbp.shape

I get this bug:
'''

ValueError Traceback (most recent call last)
Cell In[40], line 6
4 img_tv = ossart_tv(limited_projection,geo,limited_angle,5)
5 img_tv.shape
----> 6 img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)
7 img_fbp.shape

File c:\Users\Haodong_Li\AppData\Local\anaconda3\envs\pytorch\lib\sitepackages\tigre\algorithms\single_pass_algorithms.py:195, in fbp(proj, geo, angles, **kwargs)
192 gpuids = kwargs["gpuids"] if "gpuids" in kwargs else None
193 proj_filt = filtering(copy.deepcopy(proj), geox,
194 angles, parker=False, verbose=verbose)
--> 195 return Atb(proj_filt, geo, angles, gpuids=gpuids) * geo.DSO / geo.DSD

ValueError: operands could not be broadcast together with shapes (1,256,256) (180,)
'''

[AnderBiguri]

I can not reproduce the error.

Can you show a minimal example using the TIGRE head data that also trhows an error, so I can copy paste it?

[GreameLee]

import tigre
import numpy as np
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise
import tigre.algorithms as algs
from scipy.io import loadmat
import matplotlib.pyplot as plt
from tigre.algorithms.iterative_recon_alg import IterativeReconAlg
from tigre.algorithms.iterative_recon_alg import decorator
from tigre.utilities.im_3d_denoise import im3ddenoise

class OSSART_TV(IterativeReconAlg):  
    __doc__ = (
        "OSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n"
        "Simultaneous Algebraic Reconxtruction 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 descrived 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)


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

geo = tigre.geometry()
# VARIABLE                                   DESCRIPTION                    UNITS
# -------------------------------------------------------------------------------------
# Distances
geo.DSD = 1536  # Distance Source Detector      (mm)
geo.DSO = 1000  # Distance Source Origin        (mm)
# Image parameters
geo.nVoxel = np.array([1, 256, 256])  # number of voxels              (vx)
geo.sVoxel = np.array([1, 256, 256])  # total size of the image       (mm)
geo.dVoxel = geo.sVoxel / geo.nVoxel  # size of each voxel            (mm)
print(geo.dVoxel)
# Detector parameters
geo.nDetector = np.array([1, 512])  # number of pixels              (px)
geo.dDetector = np.array([geo.dVoxel[0], 0.8])  # size of each pixel            (mm)
geo.sDetector = geo.nDetector * geo.dDetector  # total size of the detector    (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)
# MAKE SURE THAT THE DETECTOR PIXELS SIZE IN V IS THE SAME AS THE IMAGE!

geo.mode = "parallel"

#%% Define angles of projection and load phatom image

angles = np.linspace(0, 2 * np.pi, 180)

head = sample_loader.load_head_phantom(geo.nVoxel)

limited_angle = np.linspace(0, np.pi/2, 180)
limited_projection = tigre.Ax(img, geo, limited_angle)
print(limited_projection.shape)
img_tv = ossart_tv(limited_projection,geo,limited_angle,5)
print(img_tv.shape)
img_fbp = tigre.algorithms.fbp(limited_projection,geo,limited_angle)
print(img_fbp.shape)

[AnderBiguri]

I see.

try modifying the line that errors to:

return Atb(proj_filt, geo, angles, gpuids=gpuids) * geo.DSO[0] / geo.DSD[0]

This may solve the issue.

[GreameLee]

But after I modified the code as your comment, when I run single fbp algorithm I can see this bug:
'''
Traceback (most recent call last):
File "tigre_ct.py", line 72, in
fbp = tigre.algorithms.fbp(sino, geo, angles)
File "/home/haodong/anaconda3/envs/mbir/lib/python3.8/site-packages/tigre/algorithms/single_pass_algorithms.py", line 195, in fbp
return Atb(proj_filt, geo, angles, gpuids=gpuids) * geo.DSO[0] / geo.DSD[0]
TypeError: 'int' object is not subscriptable
'''

I guess there is the problem that you need delete gpuids after TIGRE using each CT reconstruction algorithm

[AnderBiguri]

Hi @GreameLee, I am currently traveling but I'll fix this ASAP. I think it's just a mix of vectors/ints there, should be an easy fix, but can't do it until I have full access to a computer with a gpu.

[GreameLee]

Hi @GreameLee, I am currently traveling but I'll fix this ASAP. I think it's just a mix of vectors/ints there, should be an easy fix, but can't do it until I have full access to a computer with a gpu.

I just try to fix this by myself, it worked if I changed the

 return Atb(proj_filt, geo, angles, gpuids=gpuids) * geo.DSO[0] / geo.DSD[0]

to

if isinstance(geo.DSO, int) == False:
        a = Atb(proj_filt, geo, angles, gpuids=gpuids)* geo.DSO[0] / geo.DSD[0]
    else:
        a = Atb(proj_filt, geo, angles, gpuids=gpuids)* geo.DSO / geo.DSD
    return a