Index error

[justingm]

Hi!

I am trying to reconstruct my own data using TIGRE. I am sure this is just due to my ignorance but I am getting the following error

IndexError: index 2 is out of bounds for axis 0 with size 2 Exception ignored in: 'tigre.utilities.cuda_interface._types.convert_to_c_geometry' Traceback (most recent call last): File "/home/justin/.local/lib/python3.8/site-packages/pytigre-2.2.0-py3.8-linux-x86_64.egg/tigre/utilities/Atb.py", line 36, in Atb return _Atb_ext(projections, geox, geox.angles, backprojection_type, geox.mode, gpuids=gpuids) IndexError: index 2 is out of bounds for axis 0 with size 2

I am not quite sure where things went wrong since I was just following the examples. Thanks in advance!

Best, Justin

[justingm]

Please ignore my previous post. I managed to find out what was causing the error. I just mistyped one the parameters.

I do have another problem. I was able to reconstruct my CT but I am seeing weird artifacts. Here's a comparison of the reconstructed image from TIGRE and the vendor software. In TIGRE, there seems to be some sort of shadow around the rods.

Just to give you the full picture, I have the projections of a CT that I took and each has a dimension of 3200 x 3680. When reconstructed at full resolution, the resulting voxel size is ~9 um. Since this is of very high resolution, I downsampled them by a factor of 8 to fit the memory that I have. Looking at the reconstruction log file from the vendor software, I found the following details:

Source to detector array distance: 336.67 mm Source to center of rotation distance: 102.08 mm Size of detector array (Z x Y): 100.948 x 116.09 mm -> 3200 x 3680 pixels Size of image voxels: 9.564 x 9.564 x 9.564 um

Using these, I tried reconstructing the downsampled images in TIGRE

--- TIGRE Reconstruction ---

% Create geometry object

geo = tigre.geometry() geo.mode = "cone"

% Distance

geo. DSD = 336.67 geo.DSO = 102.08

% Detector parameters

% NOTE: upon downsampling, size is 400x460 but the last part is full of zeros so I just excluded that so I have 459

geo.nDetector = np.array([400, 459]) geo.sDetector = np.array([100.948, 116.09]) geo.dDetector = geo.sDetector/geo.nDetector

% Image parameters

% NOTE: I am not really sure what to put in the nVoxel. From the reconstruction log that I have, it was x=3680, y=3680, z=3200 wherein z is the number of slices. I just divided that by 8 since I downsampled the data and then, set nVoxel as 400x459x459

orig_dVoxel = 0.0095645 # [mm] geo.nVoxel = np.array([400,459,459]) geo.dVoxel = np.array([orig_dVoxel, orig_dVoxel, orig_dVoxel])8 geo.sVoxel = geo.nVoxelgeo.dVoxel

% Offsets

geo.offOrigin = np.array([0,0,0]) geo.offDetector = np.array([0,0])

% Angles

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

% Import projections

image = sitk.ReadImage(file) projection = sitk.GetArrayFromImage(image) # [360,400, 459]

% Run reconstruction

imageFDK = alg.fdk(projection, geo, angles)

I know that I set the offsets to zero and in the reconstruction log of the vendor, there's Yoffset = -41.5159 mm and Zoffset=0, which corresponds to the CT projection average center offset in mm. I tried plugging it in the geo.offDetector =[0,-41.5159] but the results become very weird (see image below). If I do it the other way around geo.offDetector =[-41.5159, 0], then there is no difference from before and still with the circular shadows. Now, I am a bit stucked.

[AnderBiguri]

Glad that you found a fix!

For your current issue, it seems that you are likely using a machine with a rotating stage, instead of a rotating source-detector. Your artifacts look exactly like geo.cor artifacts, caused by a shift on the center of rotation. Unfortunately as for now, TIGRE does not have an a utomatic COR detection capability, so you will need to manually input it. I suggest having a look at your vendor software metadata, as they tend to store their aproximated COR value in the recon metadata.

[justingm]

Hello Ander,

You are right. The geo.cor did the trick! I have another question though. Is TIGRE capable of reconstructing the image by parts? I guess I can also just do that by myself and just stitch it afterwards. Although I am not sure if it is that straightforward.

[AnderBiguri]

@justingm Well, depends what you mean by parts! In cone beam CT, the problem is only divisible by upper and lower half, the rest is all (mostly) coupled, i.e. you need the entire upper half of the projection data to reconstruct any slice of the upper half of the image.

In any case, yes you can do it. Just give the correct image/projection offset and sizes. For example, for COR estimations, I tend to reconstruct 10 different CORs and select the best, but I only need about 40 center slices of the projection and 2 slices of the volume to do so (as COR is z-independent), so I describe a geometry like that and can reconstruct really fast.

There are demos on how to define complex geometries, have a look!

[justingm]

I wanted to reconstruct it at the highest resolution but due to memory issues, I couldn't, so I was thinking of splitting it to smaller volumes and just stitching the reconstructed image afterwards. Anyway, I will look at the demos and see if I can work it out. Thanks!

[AnderBiguri]

@justingm ah yes. You can always reconstruct some center slices, or as I said, top half and bottom half. But you can not split the rest much more unless you have a very very high DSD/DSO values, such that the cone is closer to a parallel beam. Its all a bit of geometry I guess.