ADMM Reconstruction

Algorithms
,
1

I successfully ran the FirstCudaReconstruction.py file on my project. But sometimes the scanning angle may only be 90 °, and the reconstruction effect using FDK is not good. I checked the rtk and found that ADMM reconstruction may be used for non full week scan reconstruction.

I referred to rtkadmmtotalvariationtest.cxx, but encountered an error.

ImageType = itk.CudaImage[itk.F, 3]
GradientOutputImageType = itk.CudaImage[itk.CovariantVector[itk.F, 3], 3]
ADMMTotalVariationType = itk.RTK.ADMMTotalVariationConeBeamReconstructionFilter[ImageType, GradientOutputImageType]

And the error prompt looks like this.

itk.support.extras.TemplateTypeError: itk.ADMMTotalVariationConeBeamReconstructionFilter is not wrapped for input type `itk.CudaImage[itk.F,3], itk.CudaImage[itk.CovariantVector[itk.F,3],3]`.

How can I create an correct object about ADMMTotalVariationConeBeamReconstructionFilter? Any help would be appreciated.

2

The wrapping deduces some of the templates from default values:

It should work with itk.RTK.ADMMTotalVariationConeBeamReconstructionFilter[ImageType].
Don’t you have the following full error message to indicate the reason for the problem?

Traceback (most recent call last):
  File "/home/srit/src/itk/lin64-PythonWrapping-bindings/itk/support/template_class.py", line 525, in __getitem__
    this_item = self.__template__[key]
KeyError: (<class 'itk.itkCudaImagePython.itkCudaImageF3'>, <class 'itk.itkCudaImagePython.itkCudaImageCVF33'>)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/srit/src/itk/lin64-PythonWrapping-bindings/itk/support/template_class.py", line 529, in __getitem__
    raise itk.TemplateTypeError(self, key)
itk.support.extras.TemplateTypeError: itk.ADMMTotalVariationConeBeamReconstructionFilter is not wrapped for input type `itk.CudaImage[itk.F,3], itk.CudaImage[itk.CovariantVector[itk.F,3],3]`.

To limit the size of the package, only a limited number of
types are available in ITK Python. To print the supported
types, run the following command in your python environment:

    itk.ADMMTotalVariationConeBeamReconstructionFilter.GetTypes()

Possible solutions:
* If you are an application user:
** Convert your input image into a supported format (see below).
** Contact developer to report the issue.
* If you are an application developer, force input images to be
loaded in a supported pixel type.

    e.g.: instance = itk.ADMMTotalVariationConeBeamReconstructionFilter[itk.CudaImage[itk.F,3]].New(my_input)

* (Advanced) If you are an application developer, build ITK Python yourself and
turned to `ON` the corresponding CMake option to wrap the pixel type or image
dimension you need. When configuring ITK with CMake, you can set
`ITK_WRAP_${type}` (replace ${type} with appropriate pixel type such as
`double`). If you need to support images with 4 or 5 dimensions, you can add
these dimensions to the list of dimensions in the CMake variable
`ITK_WRAP_IMAGE_DIMS`.

Supported input types:

itk.CudaImage[itk.F,3]
itk.Image[itk.F,3]
itk.Image[itk.D,3]
1 Like
3

Thank you for your reply. I have tried itk.RTK.ADMMTotalVariationConeBeamReconstructionFilter[ImageType]. But I’m not sure is there need pip other wrapping like use itk.CudaImage, so I asked this question there for some suggestions. It’s really help for me.
By the way, does RTK support reconstructing partial data first, so that scanned data can be reconstructed simultaneously with CT scans, achieving the goal of accelerating reconstruction speed.

4

On the rtkadmmtotalvariationtest.cxx, it’s seems that the calculation method of origins is not formula origin = spacing*(size-1)*(-0.5) like FDK reconstruction. How can I set the value of origin?
On the other word, I want run ADMM on GPU, so I use SetPixelContainer to transfer my data to GPU. There is an error happened, and that error didn’t happen on CPU.

    admmtotalvariation.Update()
RuntimeError: C:\P\IPP\ITK-source\ITK\Modules\Core\Common\include\itkImageConstIterator.h:211:
ITK ERROR: Region ImageRegion (000000CDBE4B2330)
  Dimension: 3
  Index: [0, 468, 111]
  Size: [752, 156, 1]
 is outside of buffered region ImageRegion (000002B7B213EEA0)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [0, 0, 0]
5

Yes, you can start the reconstruction before the end of the acquisition but it is a bit tricky and I removed the example showing it with this commit (to reduce the maintenance burden).

The formula is used to center the detector around the origin (0,0) but it does not have to be perfectly centered as long as the simulation is consistent with the reconstruction.

There is not enough context to understand what is going here. Did you update the reader beforehand?

1 Like
6

Yes, I have update the reader. There is my code about this section:

def loaImage(rawImages, width, height, spacing):
    itk_images = itk.GetImageFromArray(rawImages)
    itk_images.SetSpacing([spacing, spacing, spacing])
    itk_images.SetOrigin([spacing*(width-1)*(-0.5), spacing*(height-1)*(-0.5), 0])
    itk_images.Update()
    print(itk_images)
    return itk_images
reader = loaImage(images, width, height, spacing_origin)
newReader = itk.CastImageFilter[itk.Image[itk.US, 3], itk.Image[itk.F, 3]].New()
newReader.SetInput(reader)
newReader.Update()
projections = CudaImage.New()
projections.SetPixelContainer(newReader.GetOutput().GetPixelContainer())
projections.CopyInformation(newReader.GetOutput())
projections.SetBufferedRegion(newReader.GetOutput().GetBufferedRegion())
projections.SetRequestedRegion(newReader.GetOutput().GetRequestedRegion())

ConstantImageSourceType = itk.RTK.ConstantImageSource[CudaImage]
# Create reconstructed image
constantImageSource2 = ConstantImageSourceType.New()
sizeOutput = [outX, outY, outZ]  # [128, 128, 128]
scale = [images.shape[2] / outX, images.shape[2] / outY, images.shape[1] / outZ]
spacing = [spacing_origin / (sdd / sid) * scale[0], spacing_origin / (sdd / sid) * scale[1],
               spacing_origin / (sdd / sid) * scale[2]]
origin = [-spacing[0] * (outX - 1) * 0.5, -spacing[1] * (outY - 1) * 0.5,
              -spacing[2] * (outZ - 1) * 0.5]  # VolspaingX*(VolX-1)*-0.5
constantImageSource2.SetOrigin(origin)
constantImageSource2.SetSpacing(spacing)
constantImageSource2.SetSize(sizeOutput)
constantImageSource2.SetConstant(0.)

ADMMTotalVariationType = itk.RTK.ADMMTotalVariationConeBeamReconstructionFilter[CudaImage]
admmtotalvariation = ADMMTotalVariationType.New()
admmtotalvariation.SetInput(0, constantImageSource2.GetOutput())
admmtotalvariation.SetInput(1, projections)
admmtotalvariation.SetGeometry(geometry)
admmtotalvariation.SetAlpha(100)
admmtotalvariation.SetBeta(1000)
admmtotalvariation.SetAL_iterations(3)
admmtotalvariation.SetCG_iterations(2)
admmtotalvariation.SetForwardProjectionFilter(ADMMTotalVariationType.ForwardProjectionType_FP_CUDARAYCAST)
admmtotalvariation.SetBackProjectionFilter(ADMMTotalVariationType.BackProjectionType_BP_CUDAVOXELBASED)
admmtotalvariation.Update()
7

I don’t really know… First comment, in

I don’t think this Update() is required. Second, can you try if disconnecting the output of newReader from the pipeline helps? Something like:

newReader.Update()
img = newReader.GetOutput()
img.DisconnectPipeline()
projections = CudaImage.New()
projections.SetPixelContainer(img.GetPixelContainer())
projections.CopyInformation(img)
projections.SetBufferedRegion(img.GetBufferedRegion())
projections.SetRequestedRegion(img.GetRequestedRegion())
8

I tried this method but it didn’t solve the problem.
I use this pipeline on ADMM CPU reconstruction, there’s no problem here either.

ADMMTotalVariationType = itk.RTK.ADMMTotalVariationConeBeamReconstructionFilter[ImageType]
admmtotalvariation = ADMMTotalVariationType.New()
admmtotalvariation.SetInput(constantImageSource2.GetOutput())
admmtotalvariation.SetInput(1, newReader.GetOutput())
admmtotalvariation.SetGeometry(geometry)
admmtotalvariation.SetAlpha(100)
admmtotalvariation.SetBeta(1000)
admmtotalvariation.SetAL_iterations(3)
admmtotalvariation.SetCG_iterations(2)
admmtotalvariation.SetForwardProjectionFilter(ADMMTotalVariationType.ForwardProjectionType_FP_JOSEPH)
admmtotalvariation.SetBackProjectionFilter(ADMMTotalVariationType.BackProjectionType_BP_VOXELBASED)
admmtotalvariation.Update()

So is there a problem with the way I transfer data to the GPU? I don’t know where the problem is because this pipepline was useful on FDK GPU reconstruction.
By monitoring the system memory, there was a momentary high GPU usage, and it is unclear if this is the problem. My data size is about 920KB*180, GPU size is 15GB.
20240306152409

9

May I ask if it is convenient to send me a copy of this file? My email is 1533168486@qq.com, very appriciate for you.

10

You can download the original file in the commit :slight_smile: https://raw.githubusercontent.com/SimonRit/RTK/38a9f44e2c8a3311735f627e1150b3ff7a94d169/applications/rtkinlinefdk/rtkinlinefdk.cxx

11

I’m sorry that I can’t open the link as you provided. It will be helpful if you can send the file by email.

12

The link can be opened by this method

13

I just checked your data and code that you sent privately. The volume you’re trying to reconstruct has a sizeOutput of [752, 752, 624], i.e., 10 GB. Make it [128, 128, 128] and your code runs fine. Don’t forget that advanced reconstruction algorithms require multiple copies of the volume.

1 Like
14

Thank you very much for your help.