<div dir="ltr"><div>Hi,</div><div>Maybe you can first check whether the projections are read adequately by adding a line itk.imwrite(projectionsSource)? The ProjectionsReader tries to automatically convert the projections to line integral and it might do something wrong... I don't see what could be wrong in your code.<br></div><div>Simon<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Feb 28, 2023 at 11:53 PM Eldridge, Bryce D <<a href="mailto:bdeldri@sandia.gov">bdeldri@sandia.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg8033272720856042634">





<div style="overflow-wrap: break-word;" lang="EN-US">
<div class="m_8033272720856042634WordSection1">
<p class="MsoNormal">Hi,<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">I am trying to modify the first reconstruction example to read a set of 16-bit grayscale TIFF projections that I captured, run the 3D reconstruction, and save the output. The projects are taken at 5 degree increments over a 180 degree arc.
<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">It runs without errors, but the reconstructed volume looks very strange which makes me think I am doing something wrong.<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Here is the basic python code I am trying to use (modified from the example), if anyone can see what I might be doing wrong I would appreciate it.
<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"># Defines the image type<u></u><u></u></p>
<p class="MsoNormal">ImageType = itk.Image[itk.F,3]<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"># Defines the RTK geometry object<u></u><u></u></p>
<p class="MsoNormal">geometry = rtk.ThreeDCircularProjectionGeometry.New()<u></u><u></u></p>
<p class="MsoNormal">numberOfProjections = 37<u></u><u></u></p>
<p class="MsoNormal">firstAngle = 0.<u></u><u></u></p>
<p class="MsoNormal">angularArc = 185.<u></u><u></u></p>
<p class="MsoNormal">sid = 914 # source to isocenter distance (mm)<u></u><u></u></p>
<p class="MsoNormal">sdd = 1219 # source to detector distance (mm)<u></u><u></u></p>
<p class="MsoNormal">for x in range(0,numberOfProjections):<u></u><u></u></p>
<p class="MsoNormal">  angle = firstAngle + x * angularArc / numberOfProjections<u></u><u></u></p>
<p class="MsoNormal">  geometry.AddProjection(sid,sdd,angle)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Skipping some code … (file name list setup, geometry output writing) …<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">tiffio = itk.TIFFImageIO.New()<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">ProjectionsReaderType = rtk.ProjectionsReader[ImageType]<u></u><u></u></p>
<p class="MsoNormal">projectionsSource = ProjectionsReaderType.New()<u></u><u></u></p>
<p class="MsoNormal">projectionsSource.SetImageIO(tiffio)<u></u><u></u></p>
<p class="MsoNormal">projectionsSource.SetFileNames(fileNameList)<u></u><u></u></p>
<p class="MsoNormal">projOrigin = [ -0.14*(3072-1)/2, -0.14*(2560-1)/2, 0 ] #input images are 3072x2560 pixels with a 0.14mm pixel size<u></u><u></u></p>
<p class="MsoNormal">projSpacing = [ 0.14, 0.14, 1.0 ]<u></u><u></u></p>
<p class="MsoNormal">projectionsSource.SetOrigin( projOrigin )<u></u><u></u></p>
<p class="MsoNormal">projectionsSource.SetSpacing( projSpacing )<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">ConstantImageSourceType = rtk.ConstantImageSource[ImageType]<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"># Create reconstructed image<u></u><u></u></p>
<p class="MsoNormal">constantImageSource2 = ConstantImageSourceType.New()<u></u><u></u></p>
<p class="MsoNormal">sizeOutput = [ 512, 512, 512 ]<u></u><u></u></p>
<p class="MsoNormal">origin = [ -255.5, -255.5, -255.5 ]<u></u><u></u></p>
<p class="MsoNormal">spacing = [ 1.0, 1.0, 1.0 ]<u></u><u></u></p>
<p class="MsoNormal">constantImageSource2.SetOrigin( origin )<u></u><u></u></p>
<p class="MsoNormal">constantImageSource2.SetSpacing( spacing )<u></u><u></u></p>
<p class="MsoNormal">constantImageSource2.SetSize( sizeOutput )<u></u><u></u></p>
<p class="MsoNormal">constantImageSource2.SetConstant(0.)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"># FDK reconstruction<u></u><u></u></p>
<p class="MsoNormal">print("Reconstructing...")<u></u><u></u></p>
<p class="MsoNormal">FDKCPUType = rtk.FDKConeBeamReconstructionFilter[ImageType]<u></u><u></u></p>
<p class="MsoNormal">feldkamp = FDKCPUType.New()<u></u><u></u></p>
<p class="MsoNormal">feldkamp.SetInput(0, constantImageSource2.GetOutput()) # this is the template for the output image type<u></u><u></u></p>
<p class="MsoNormal">feldkamp.SetInput(1, projectionsSource.GetOutput()) # this is the projection stack from rtk.ProjectionsReader<u></u><u></u></p>
<p class="MsoNormal">feldkamp.SetGeometry(geometry)<u></u><u></u></p>
<p class="MsoNormal">feldkamp.GetRampFilter().SetTruncationCorrection(0.0)<u></u><u></u></p>
<p class="MsoNormal">feldkamp.GetRampFilter().SetHannCutFrequency(0.0)<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"># Writer<u></u><u></u></p>
<p class="MsoNormal">print("Writing output image...")<u></u><u></u></p>
<p class="MsoNormal">WriterType = rtk.ImageFileWriter[ImageType]<u></u><u></u></p>
<p class="MsoNormal">writer = WriterType.New()<u></u><u></u></p>
<p class="MsoNormal">writer.SetFileName(sys.argv[1])<u></u><u></u></p>
<p class="MsoNormal">writer.SetInput(feldkamp.GetOutput())<u></u><u></u></p>
<p class="MsoNormal">writer.Update()<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Thanks<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
</div>

_______________________________________________<br>
Rtk-users mailing list<br>
<a href="mailto:rtk-users@openrtk.org" target="_blank">rtk-users@openrtk.org</a><br>
<a href="https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users" rel="noreferrer" target="_blank">https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users</a><br>
</div></blockquote></div>