Hi,
Thanks for sharing your experience. I tried to use your code to reproduce
your problem but I couldn't. Enclosed is the code that I run with
import daye
simulator = daye.simulator()
simulator.runCBCT('proj.mha', 0, 0)
simulator.buildVolumeFromFile('proj.mha','recon.mha',simulator)
Does it work for you? If yes, can you send your full example?
Simon
On Mon, Feb 18, 2019 at 7:51 AM Pierre Daye <[email protected]>
wrote:
>
>
> Hello,
>
>
>
> First of all congrats for the update of RTK and the new Python
> integration. This is really great for Python users who don’t want the whole
>
> C++ armada to test a small reconstruction principle!
>
> I could easily run the first reconstruction and it worked perfectly!
> However, I tried then to first save the images using these method within
>
> a simulator python class:
>
>
>
> def buildGeometry(self, offsetZ, offsetY):
>
> geometry = rtk.ThreeDCircularProjectionGeometry.New()
>
> for x in range(0, self.param['numberOfProjections']):
>
> angle = self.param['firstAngle'] + x *
> self.param['angularArc'] / self.param['numberOfProjections']
>
> geometry.AddProjection(self.param['sid'] + offsetZ,
>
> self.param['sdd'], angle, self.param['isox'],
>
> self.param['isoy'], self.param['outOfPlaneAngle'],
>
> self.param['inPlaneAngle'], self.param['sourceOffsetX'],
> offsetY)
>
>
>
> return geometry
>
>
>
> def runCBCT(self, filename, offsetZ, offsetY):
>
> geometry = self.buildGeometry(offsetZ, offsetY)
>
>
>
> constantImageSource =
> rtk.ConstantImageSource[self.param['TImageType']].New()
>
> constantImageSource.SetOrigin( self.param['origin'] )
>
> constantImageSource.SetSpacing( self.param['spacing'] )
>
> constantImageSource.SetSize( self.param['sizeOutput'] )
>
> constantImageSource.SetConstant(0.0)
>
> source = constantImageSource.GetOutput()
>
>
>
> self.rei.SetGeometry(geometry)
>
> self.rei.SetInput(source)
>
>
>
> projections = self.rei.GetOutput()
>
>
>
> writer = itk.ImageFileWriter[self.param['TImageType']].New()
>
> writer.SetFileName(filename)
>
> writer.SetInput(projections)
>
> writer.Update()
>
>
>
> This parts works perfectly and the generated images are OK (checked with
> IMageJ).
>
> The problems start when I want to reconstruct the volume.
>
>
>
> def buildVolumeFromFile(self, file, fileOUT, CBCT):
>
> geometry = CBCT.buildGeometry(0, 0)
>
> print("Performing reconstruction")
>
> TImageType = self.param['TImageType']
>
> reader = itk.ImageFileReader[TImageType].New()
>
> reader.SetFileName(file)
>
> reiImage = reader.GetOutput()
>
>
>
> # Create reconstructed image
>
> constantImageSource2 = rtk.ConstantImageSource[TImageType].New()
>
> origin = [ -63.5, -63.5, -63.5 ]
>
> sizeOutput = [ 128, 128, 128 ]
>
> constantImageSource2.SetOrigin( origin )
>
> constantImageSource2.SetSpacing( [1, 1, 1] )
>
> constantImageSource2.SetSize( sizeOutput )
>
> constantImageSource2.SetConstant(0.0)
>
> source2 = constantImageSource2.GetOutput()
>
>
>
> print("Performing reconstruction")
>
> feldkamp = rtk.FDKConeBeamReconstructionFilter[TImageType].New()
>
> feldkamp.SetGeometry( geometry )
>
> feldkamp.SetInput(0, source2)
>
> feldkamp.SetInput(1, reiImage)
>
> image = feldkamp.GetOutput()
>
>
>
> print("Masking field-of-view")
>
> fov = rtk.FieldOfViewImageFilter[TImageType, TImageType].New()
>
> fov.SetGeometry(geometry)
>
> fov.SetProjectionsStack(reiImage)
>
> fov.SetInput(image)
>
> image = fov.GetOutput()
>
>
>
> writer = itk.ImageFileWriter[TImageType].New()
>
> writer.SetFileName ( fileOUT )
>
> writer.SetInput(image)
>
> writer.Update()
>
>
>
> I cannot run this script. It generates a criptic segfault… What am I doing
> wrong here?
>
> If one of you RTK gurus could help me, I would greatly appreciate it!
>
>
>
> Thanks,
>
>
>
> Pierre
>
> Disclaimer | Use of IBA e-communication
> <https://iba-worldwide.com/disclaimer>
>
> The contents of this e-mail message and any attachments are intended
> solely for the recipient (s) named above. This communication is intended to
> be and to remain confidential and may be protected by intellectual property
> rights. Any use of the information contained herein (including but not
> limited to, total or partial reproduction, communication or distribution of
> any form) by persons other than the designated recipient(s) is prohibited.
> Please notify the sender immediately by e-mail if you have received this
> e-mail by mistake and delete this e-mail from your system. E-mail
> transmission cannot be guaranteed to be secure or error-free. Ion Beam
> Applications does not accept liability for any such errors. Thank you for
> your cooperation.
> _______________________________________________
> Rtk-users mailing list
> [email protected]
> https://public.kitware.com/mailman/listinfo/rtk-users
>
import itk
from itk import RTK as rtk
class simulator:
def __init__(self):
self.param = dict()
self.param['firstAngle'] = 0
self.param['angularArc'] = 360
self.param['numberOfProjections'] = 32
self.param['sid'] = 1000
self.param['sdd'] = 2000
self.param['isox'] = 0
self.param['isoy'] = 0
self.param['inPlaneAngle'] = 0
self.param['outOfPlaneAngle'] = 0
self.param['sourceOffsetX'] = 0
self.param['origin'] = -255
self.param['spacing'] = 2
self.param['sizeOutput'] = [256,256,self.param['numberOfProjections']]
self.param['TImageType'] = itk.Image[itk.F,3]
self.rei = rtk.RayEllipsoidIntersectionImageFilter[itk.Image[itk.F,3],itk.Image[itk.F,3]].New()
self.rei.SetAxis(20)
self.rei.SetCenter([20,0,0])
def buildGeometry(self, offsetZ, offsetY):
geometry = rtk.ThreeDCircularProjectionGeometry.New()
for x in range(0, self.param['numberOfProjections']):
angle = self.param['firstAngle'] + x * self.param['angularArc'] / self.param['numberOfProjections']
geometry.AddProjection(self.param['sid'] + offsetZ,
self.param['sdd'], angle, self.param['isox'],
self.param['isoy'], self.param['outOfPlaneAngle'],
self.param['inPlaneAngle'], self.param['sourceOffsetX'], offsetY)
return geometry
def runCBCT(self, filename, offsetZ, offsetY):
geometry = self.buildGeometry(offsetZ, offsetY)
constantImageSource = rtk.ConstantImageSource[self.param['TImageType']].New()
constantImageSource.SetOrigin( self.param['origin'] )
constantImageSource.SetSpacing( self.param['spacing'] )
constantImageSource.SetSize( self.param['sizeOutput'] )
constantImageSource.SetConstant(0.0)
constantImageSource.Update()
source = constantImageSource.GetOutput()
self.rei.SetGeometry(geometry)
self.rei.SetInput(source)
projections = self.rei.GetOutput()
writer = itk.ImageFileWriter[self.param['TImageType']].New()
writer.SetFileName(filename)
writer.SetInput(projections)
writer.Update()
def buildVolumeFromFile(self, file, fileOUT, CBCT):
geometry = CBCT.buildGeometry(0, 0)
print("Performing reconstruction")
TImageType = self.param['TImageType']
reader = itk.ImageFileReader[TImageType].New()
reader.SetFileName(file)
reiImage = reader.GetOutput()
# Create reconstructed image
constantImageSource2 = rtk.ConstantImageSource[TImageType].New()
origin = [ -63.5, -63.5, -63.5 ]
sizeOutput = [ 128, 128, 128 ]
constantImageSource2.SetOrigin( origin )
constantImageSource2.SetSpacing( [1, 1, 1] )
constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetConstant(0.0)
source2 = constantImageSource2.GetOutput()
print("Performing reconstruction")
feldkamp = rtk.FDKConeBeamReconstructionFilter[TImageType].New()
feldkamp.SetGeometry( geometry )
feldkamp.SetInput(0, source2)
feldkamp.SetInput(1, reiImage)
image = feldkamp.GetOutput()
print("Masking field-of-view")
fov = rtk.FieldOfViewImageFilter[TImageType, TImageType].New()
fov.SetGeometry(geometry)
fov.SetProjectionsStack(reiImage)
fov.SetInput(image)
image = fov.GetOutput()
writer = itk.ImageFileWriter[TImageType].New()
writer.SetFileName ( fileOUT )
writer.SetInput(image)
writer.Update()
_______________________________________________
Rtk-users mailing list
[email protected]
https://public.kitware.com/mailman/listinfo/rtk-users
