SimonRit
commented
1 year ago This does not seem to be an issue but a question for the mailing list: https://www.openrtk.org/RTK/project/contactus.html. Converting a numpy array to an image and back can be done with itk.image_from_array and itk.array_from_image.
!/usr/bin/env python
import sys import itk import numpy as np from itk import RTK as rtk from crip.io import imreadRaws import SimpleITK as sitk
############ 2、投影数据 ############ projections = imreadRaws(raw_dir, 640, 640, dtype=np.uint16) projections = projections[20:557] projections = np.transpose(projections, (1,2,0)) projections = np.float64(projections)
projections = projections[:, :, ::-1]
####################################
Defines the image type
ImageType = itk.Image[itk.F,3]
Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New() numberOfProjections = 537 firstAngle = 0. angularArc = 360. sid = 433. # source to isocenter distance sdd = 660. # source to detector distance
isox = -34.3 # X coordinate on the projection image of isocenter
isoy = -40.64 # Y coordinate on the projection image of isocenter
for x in range(0,numberOfProjections): angle = firstAngle + x * angularArc / numberOfProjections geometry.AddProjection(sid,sdd,angle)
geometry.AddProjection(sid,sdd,angle,isox,isoy)
Create a stack of empty projection images
ConstantImageSourceType = rtk.ConstantImageSource[ImageType]
constantImageSource = ConstantImageSourceType.New()
origin = [ -34.3, -40.64, 0. ]
sizeOutput = [ 640, 640, numberOfProjections ]
spacing = [ 0.254, 0.254, 0.254 ]
constantImageSource.SetOrigin( origin )
constantImageSource.SetSpacing( spacing )
constantImageSource.SetSize( sizeOutput )
constantImageSource.SetConstant(0.)
output = constantImageSource.GetOutput()
projections_sitk = sitk.GetImageFromArray(projections) origin = [ -34.3, -40.64, 0. ] spacing = [ 0.254, 0.254, 0.254 ] projections_sitk.SetOrigin(origin) projections_sitk.SetSpacing(spacing)
Create reconstructed image
ConstantImageSourceType = rtk.ConstantImageSource[ImageType] constantImageSource2 = ConstantImageSourceType.New() sizeOutput = [ 672, 672, 448 ] origin = [ 0, 0, 0 ] spacing = [ 0.223, 0.223, 0.201 ] constantImageSource2.SetOrigin( origin ) constantImageSource2.SetSpacing( spacing ) constantImageSource2.SetSize( sizeOutput ) constantImageSource2.SetConstant(0.)
FDK reconstruction
print("Reconstructing...") FDKCPUType = rtk.FDKConeBeamReconstructionFilter[ImageType] feldkamp = FDKCPUType.New() feldkamp.SetInput(0, constantImageSource2.GetOutput()) feldkamp.SetInput(1, projections_sitk) feldkamp.SetGeometry(geometry)
feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
feldkamp.GetRampFilter().SetHannCutFrequency(0.0)
feldkamp.Update()
output = feldkamp.GetOutput() print(output.shape)