search

LLNL / LEAP

comprehensive library of 3D transmission Computed Tomography (CT) algorithms with Python and C++ APIs, a PyQt GUI, and fully integrated with PyTorch
https://leapct.readthedocs.io
MIT License
120 stars 12 forks source link

Question: When I use different param_id for postprocessing, image will be zero. #107

Closed Ellen-D216 closed 2 months ago

Ellen-D216 commented 2 months ago

I want to simulate different imaging conditions with different param_id. However, when I postprocessed the image under the first condition, the reconstruction using the second imaging condition resulted in zero images.

leapct = tomographicModels(param_id=1)

# geo1
geo = Geometry()
geo.angles = np.linspace(0, 360, 720).astype(np.float32) 
geo.voxel_numx, geo.voxel_numy, geo.voxel_numz = tuple(im.size)
geo.voxel_height, geo.voxel_width = im.spacing[2], im.spacing[0]
geo.det_rows, geo.det_cols = geo.voxel_numz, 768 
geo.det_center_row, geo.det_center_col = geo.voxel_numz , 768 / 2
geo.pixel_height, geo.pixel_width = im.spacing[2], im.spacing[0]

leapct.set_parallelBeam(
    len(geo.angles), 
    geo.det_rows, geo.det_cols,
    geo.pixel_height, geo.pixel_width,
    geo.det_center_row, geo.det_center_col,
    geo.angles, 
)
leapct.set_volume(
    geo.voxel_numx, geo.voxel_numy, geo.voxel_numz,
    geo.voxel_width, geo.voxel_height,
    geo.offsetx, geo.offsety, geo.offsetz
)
sino1 = np.zeros((len(geo.angles), geo.det_rows, geo.det_cols), dtype=np.float32)
leapct.project(sino1, im.array)
recon1 = leapct.FBP(sino1)
# post process
leapct.diffuse(recon1, 0.02/20, 4)
leapct.BilateralFilter(recon_post, 4, 0.02)

# geo2, just change pixel size
geo = Geometry()
geo.angles = np.linspace(0, 360, 720).astype(np.float32) 
geo.voxel_numx, geo.voxel_numy, geo.voxel_numz = tuple(im.size)
geo.voxel_height, geo.voxel_width = im.spacing[2], im.spacing[0]
geo.det_rows, geo.det_cols = geo.voxel_numz, 768 * 2
geo.det_center_row, geo.det_center_col = geo.voxel_numz , 768
geo.pixel_height, geo.pixel_width = im.spacing[2], im.spacing[0] / 2

# change param_id and repeat the process
leapct.param_id = 2
leapct.set_model()
leapct.set_parallelBeam(
    len(geo.angles), 
    geo.det_rows, geo.det_cols,
    geo.pixel_height, geo.pixel_width,
    geo.det_center_row, geo.det_center_col,
    geo.angles, 
)
leapct.set_volume(
    geo.voxel_numx, geo.voxel_numy, geo.voxel_numz,
    geo.voxel_width, geo.voxel_height,
    geo.offsetx, geo.offsety, geo.offsetz
)
sino2 = np.zeros((len(geo.angles), geo.det_rows, geo.det_cols), dtype=np.float32)
leapct.project(sino2, im.array)
recon2 = leapct.FBP(sino2) # recon2 is zero
kylechampley commented 2 months ago

Thats for the question. I feel like our documentation is not so clear as the use of param_id does not align with our intensions.

The param_id argument is for instances that already exist. You should not use this argument to make a new set of parameters. In addition, one should consider the set_model() as a private function, i.e., users should not use it.

I recommend that you specify two parameter sets as follows. leapct1 = tomographicModels() leapct2 = tomographicModels() These two objects will have independent parameters that can be different.

In addition, I think you have the order of the dimensions incorrect for the volume. The order is as follows (using notation like you have used: recon = np.zeros(voxel_numz, voxel_numy, voxel_numx, dtype=np.float32)

In your script it appears that you have these dimensions reversed.

Let me know if you have further questions.

hws203 commented 2 months ago

How about making the param_id in private member? So user can modify it by direct member access. I think that any class has a public member function of accessing a param, then that param should be private.

Ellen-D216 commented 2 months ago

Thats for the question. I feel like our documentation is not so clear as the use of param_id does not align with our intensions.

The param_id argument is for instances that already exist. You should not use this argument to make a new set of parameters. In addition, one should consider the set_model() as a private function, i.e., users should not use it.

I recommend that you specify two parameter sets as follows. leapct1 = tomographicModels() leapct2 = tomographicModels() These two objects will have independent parameters that can be different.

In addition, I think you have the order of the dimensions incorrect for the volume. The order is as follows (using notation like you have used: recon = np.zeros(voxel_numz, voxel_numy, voxel_numx, dtype=np.float32)

In your script it appears that you have these dimensions reversed.

Let me know if you have further questions.

Thank you for your advices, I've modified my code, but the result hasn't changed. After my attempts, I found that as long as the BilateralFilter function is called, the results of the tomographicModels created afterwards are all zero during the forward and reverse projection process.

leapct1 = tomographicModels()

# geo1
geo = Geometry()
geo.angles = np.linspace(0, 360, 720).astype(np.float32) 
geo.voxel_numz, geo.voxel_numy, geo.voxel_numx = tuple(im.shape)
geo.voxel_height, geo.voxel_width = im.spacing[2], im.spacing[0]
geo.det_rows, geo.det_cols = geo.voxel_numz, 768 
geo.det_center_row, geo.det_center_col = geo.voxel_numz , 768 / 2
geo.pixel_height, geo.pixel_width = im.spacing[2], im.spacing[0]

leapct1.set_parallelBeam(
    len(geo.angles), 
    geo.det_rows, geo.det_cols,
    geo.pixel_height, geo.pixel_width,
    geo.det_center_row, geo.det_center_col,
    geo.angles, 
)
leapct1.set_volume(
    geo.voxel_numx, geo.voxel_numy, geo.voxel_numz,
    geo.voxel_width, geo.voxel_height,
    geo.offsetx, geo.offsety, geo.offsetz
)
sino1 = np.zeros((len(geo.angles), geo.det_rows, geo.det_cols), dtype=np.float32)
leapct1.project(sino1, im.array.astype(np.float32))
recon1 = leapct1.FBP(sino1)
# post process
recon_post1 = recon1.copy().astype(np.float32)
leapct1.diffuse(recon_post1, 0.02/20, 4)
# If I call the function, sino2 and recon2 are all zero.
# leapct1.BilateralFilter(recon_post1, 4, 0.02)

leapct2 = tomographicModels()
# geo2, just change pixel size
geo = Geometry()
geo.angles = np.linspace(0, 360, 720).astype(np.float32) 
geo.voxel_numz, geo.voxel_numy, geo.voxel_numx = tuple(im.shape)
geo.voxel_height, geo.voxel_width = im.spacing[2], im.spacing[0]
geo.det_rows, geo.det_cols = geo.voxel_numz, 768 * 2
geo.det_center_row, geo.det_center_col = geo.voxel_numz , 768
geo.pixel_height, geo.pixel_width = im.spacing[2], im.spacing[0] / 2

leapct2.set_parallelBeam(
    len(geo.angles), 
    geo.det_rows, geo.det_cols,
    geo.pixel_height, geo.pixel_width,
    geo.det_center_row, geo.det_center_col,
    geo.angles, 
)
leapct2.set_volume(
    geo.voxel_numx, geo.voxel_numy, geo.voxel_numz,
    geo.voxel_width, geo.voxel_height,
    geo.offsetx, geo.offsety, geo.offsetz
)
sino2 = np.zeros((len(geo.angles), geo.det_rows, geo.det_cols), dtype=np.float32)
leapct2.project(sino2, im.array.astype(np.float32))
recon2 = leapct2.FBP(sino2)

recon_post2 = recon2.copy().astype(np.float32)
leapct2.diffuse(recon_post2, 0.02/20, 4)
# leapct2.BilateralFilter(recon_post2, 4, 0.02)
kylechampley commented 2 months ago

It's hard for me to reproduce your issue because I do not have access to this "Geometry" class or "im". I tried to make modifications to your script to avoid these missing features and things seemed to work for me with one exception. I still see that you have the ordering of things reversed. For example, this line:

geo.pixel_height, geo.pixel_width = im.spacing[2], im.spacing[0]

needs to be like this:

geo.pixel_width, geo.pixel_height = im.spacing[2], im.spacing[0]

Don't forget to also change this for leapct2. Please change these lines and run your code again.

kylechampley commented 2 months ago

How about making the param_id in private member? So user can modify it by direct member access. I think that any class has a public member function of accessing a param, then that param should be private.

@hws203, this is a good suggestion; I might do this.

kylechampley commented 2 months ago

@Ellen-D216, one more thing. Do you see any warning or error messages printed to the screen?

Ellen-D216 commented 2 months ago

@Ellen-D216, one more thing. Do you see any warning or error messages printed to the screen?

There are no messages. I've compiled the Release version. Maybe I should try compiling the Debug version to see if there are any errors or warnings.

Ellen-D216 commented 2 months ago

This is likely an issue with my own code. The results I get using the native API are normal.

kylechampley commented 2 months ago

I'm not sure what you mean by the "release" and the "debug" versions, but this is not the issue. Did you swap indices of im.spacing like I suggested? I think this is the most likely reason for your issue.

Ellen-D216 commented 2 months ago

Yes, I changed the indices, but the results didn't meet my expectations. I'll keep looking into it.

Ellen-D216 commented 2 months ago

I reproduced the issue without "Geometry", here is my code. recon1 is normal, recon2 is zero.

import numpy as np 
import SimpleITK as sitk 
import matplotlib.pyplot as plt 
from leap_tools.leapctype import tomographicModels

im = sitk.ReadImage('image.nii.gz')
shape = im.GetSize() # order is x,y,z
spacing = im.GetSpacing() # order is x,y,z

angles = np.linspace(0, 360, 720, dtype=np.float32)
image_array = sitk.GetArrayFromImage(im) # array order is z,y,x

leapct1 = tomographicModels()
leapct1.set_parallelBeam(
    numAngles=len(angles),
    numRows=shape[2], numCols=shape[0], # numRows is equal to the number of z (shape[2])
    pixelHeight=spacing[2], pixelWidth=spacing[0], # pixelHeight is also equal to the spacing of z (spacing[2])
    centerRow=shape[2]/2, centerCol=shape[0]/2,
    phis=angles
)
leapct1.set_volume(
    numX=shape[0], numY=shape[1], numZ=shape[2],
    voxelHeight=spacing[2], voxelWidth=spacing[0]
)
sino1 = np.zeros((len(angles), shape[2], shape[0]), dtype=np.float32)
leapct1.project(sino1, image_array.astype(np.float32))
recon1 = np.zeros_like(image_array, dtype=np.float32)
leapct1.FBP(sino1, recon1)
leapct1.diffuse(recon1, 0.02/20, 4)
leapct1.BilateralFilter(recon1, 4, 0.2) # True

leapct2 = tomographicModels()
leapct2.set_parallelBeam(
    numAngles=len(angles),
    numRows=shape[2], numCols=shape[0]*2,
    pixelHeight=spacing[2], pixelWidth=spacing[0]/2,
    centerRow=shape[2]/2, centerCol=shape[0],
    phis=angles
)
leapct2.set_volume(
    numX=shape[0], numY=shape[1], numZ=shape[2],
    voxelHeight=spacing[2], voxelWidth=spacing[0]
)
sino2 = np.zeros((len(angles), shape[2], shape[0]*2), dtype=np.float32)
leapct2.project(sino2, image_array.astype(np.float32))
recon2 = np.zeros_like(image_array, dtype=np.float32)
leapct2.FBP(sino2, recon2)
leapct2.diffuse(recon2, 0.02/20, 4)
leapct2.BilateralFilter(recon2, 4, 0.2) # False
Ellen-D216 commented 2 months ago

I have another question: if I have many simulation conditions (around a few hundred), can tomographicModels handle that many simulations at the same time? I'd like to use multi-threading for simulation.

kylechampley commented 2 months ago

It seems to me that you need to swap the axes of image_array like this:

image_array = np.ascontiguousarray(np.swapaxes(image_array, 0, 2), dtype=np.float32)

before running any LEAP algorithms.

Ellen-D216 commented 2 months ago

The exchange did solve the issue, but I've been simulating it this way before and it never caused any problems.