<div dir="ltr">Hi,<div><br></div><div>Thank you for your answer. I still had the issue with your updated code but I noticed I was still using itk 5.0.1 because of an <a href="https://tinyurl.com/y3dm3h86">issue</a> with itk 5.1.2 due to numpy 1.9.4 with Windows 10 recent update. I downgraded numpy to 1.9.3 so itk could install 5.1.2 and then installed the artifact you pointed for RTK. </div><div>Everything works fine now, time to play with the parameters!</div><div><br></div><div>Clément</div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le ven. 18 déc. 2020 à 08:00, Simon Rit <<a href="mailto:simon.rit@creatis.insa-lyon.fr" target="_blank">simon.rit@creatis.insa-lyon.fr</a>> a écrit :<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 dir="ltr"><div>Hi,</div><div>I slightly modified your script and it worked for me. What version are you using? I used the latest ITK v5.1.2 with, for RTK, an artifact published <a href="https://github.com/SimonRit/RTK/actions/runs/344449575" target="_blank">here</a> (visible only if you're logged in github).</div><div>Simon</div><div><br></div><div>PS: I don't think you need the displaced detector, it is centered on the central ray in your simulation.<br></div><div><br></div><div><span style="font-family:monospace">import itk<br>from itk import RTK as rtk<br><br># Defines the image type<br>ImageType = itk.Image[itk.F,3]<br><br># Defines the RTK geometry object<br>geometry = rtk.ThreeDCircularProjectionGeometry.New()<br>numberOfProjections = 180<br>firstAngle = 90<br>sid = 600 # source to isocenter distance<br>sdd = 1200 # source to detector distance<br>angularArc = 180<br>for x in range(0,numberOfProjections):<br>  angle = firstAngle + x * angularArc / numberOfProjections<br>  geometry.AddProjection(sid,sdd,angle)<br><br># Create a stack of empty projection images<br>ConstantImageSourceType = rtk.ConstantImageSource[ImageType]<br>constantImageSource = ConstantImageSourceType.New()<br>origin = [ -199, -199, 0. ]<br>sizeOutput = [ 100, 100, 180 ]<br>spacing = [ 4.0, 4.0, 2.0 ]<br>constantImageSource.SetOrigin( origin )<br>constantImageSource.SetSpacing( spacing )<br>constantImageSource.SetSize( sizeOutput )<br>constantImageSource.SetConstant(0.)<br><br># Forward projection of input volume<br>proj = rtk.JosephForwardProjectionImageFilter[ImageType, ImageType].New()<br>proj.SetInput(constantImageSource)<br>proj.SetInput(1, itk.imread('test.mha')) # input is ImageType from a reader<br>proj.SetGeometry(geometry)<br>proj.Update()<br><br># Create weights<br>weightsSource = ConstantImageSourceType.New()<br>weightsSource.SetInformationFromImage(proj.GetOutput())<br>weightsSource.SetConstant(1.0)<br><br># Create output like volume<br>outputImageSource = ConstantImageSourceType.New()<br>sizeOutput = [ 100, 100, 100]<br>origin = [ -200, -200, -200 ]<br>spacing = [ 4.0, 4.0, 4.0 ]<br>outputImageSource.SetOrigin(origin)<br>outputImageSource.SetSpacing(spacing)<br>outputImageSource.SetSize(sizeOutput)<br>outputImageSource.SetConstant(0.)<br><br># Create reconstructed image<br>RCGType = rtk.ConjugateGradientConeBeamReconstructionFilter[ImageType]<br>conjugateGradient = RCGType.New()<br>conjugateGradient.SetInput(outputImageSource.GetOutput())<br>conjugateGradient.SetInput(1, proj.GetOutput())<br>conjugateGradient.SetInput(2, weightsSource.GetOutput())<br>conjugateGradient.SetGeometry(geometry)<br>conjugateGradient.SetGamma(2.0)<br>conjugateGradient.SetNumberOfIterations(2)<br>conjugateGradient.SetDisableDisplacedDetectorFilter(False)<br>conjugateGradient.Update()</span></div><div><span style="font-family:monospace"><br></span></div><div><span style="font-family:monospace">itk.imwrite(conjugateGradient.GetOutput(), 'output.mha') </span></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Dec 16, 2020 at 10:36 AM Clément Mirabel <<a href="mailto:clement.mirabel@gmail.com" target="_blank">clement.mirabel@gmail.com</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 dir="ltr"><div>Hi,</div><div><br></div><div>I created another topic since the content is different from the previous one and people could land here with a different search in the archives.</div><div><br></div><div>While trying to reconstruct a set of images from 0 and 180° without having to add the fan angle, Simon recommended to use the conjugate gradient algorithm in RTK. So, I first tried to replicate in Python the conjugate gradient application written in C++ (<a href="https://github.com/SimonRit/RTK/tree/master/applications/rtkconjugategradient" target="_blank">https://github.com/SimonRit/RTK/tree/master/applications/rtkconjugategradient</a>). But whatever types of input I set, my python kernel dies when updating the filter without any error or debug message. I spent some time thinking about it but I can't find out what is wrong.</div><div>I have attached my code at the end of this message if you have suggestions about what would be wrong.</div><div><br></div><div>Thanks for your help,</div><div><br></div><div>Clément</div><div><br></div><div>------------------------------------------------------------------------------------------------------------</div><div><br></div><div># Defines the image type<br>ImageType = itk.Image[itk.F,3]<br></div><div><br></div><div># Defines the RTK geometry object<br>geometry = rtk.ThreeDCircularProjectionGeometry.New()<br>numberOfProjections = 180<br>firstAngle = 90<br>sid = 600 # source to isocenter distance<br>sdd = 1200 # source to detector distance<br>angularArc = 180<br>for x in range(0,numberOfProjections):<br>  angle = firstAngle + x * angularArc / numberOfProjections<br>  geometry.AddProjection(sid,sdd,angle)<br><br># Create a stack of empty projection images<br>ConstantImageSourceType = rtk.ConstantImageSource[ImageType]<br>constantImageSource = ConstantImageSourceType.New()<br>origin = [ -199, -199, 0. ]<br>sizeOutput = [ 200, 200, 180 ]<br>spacing = [ 2.0, 2.0, 2.0 ]<br>constantImageSource.SetOrigin( origin )<br>constantImageSource.SetSpacing( spacing )<br>constantImageSource.SetSize( sizeOutput )<br>constantImageSource.SetConstant(0.)<br><br># Forward projection of input volume<br>proj = rtk.JosephForwardProjectionImageFilter[ImageType, ImageType].New()<br>proj.SetInput(constantImageSource)<br>proj.SetInput(1, input) # input is ImageType from a reader<br>proj.SetGeometry(geometry)<br></div><div><br></div><div># Create weights<br>weightsSource = ConstantImageSourceType.New()<br>weightsSource.SetInformationFromImage(proj.GetOutput())<br>weightsSource.SetConstant(1.0)</div><div><br></div><div># Create output like volume<br>outputImageSource = ConstantImageSourceType.New()<br>sizeOutput = [ 400, 400, 400]<br>origin = [ -200, -200, -200 ]<br>spacing = [ 1.0, 1.0, 1.0 ]<br>outputImageSource.SetOrigin(origin)<br>outputImageSource.SetSpacing(spacing)<br>outputImageSource.SetSize(sizeOutput)<br>outputImageSource.SetConstant(0.)<br><br># Create reconstructed image</div>RCGType = rtk.ConjugateGradientConeBeamReconstructionFilter[ImageType]<br>conjugateGradient = RCGType.New()<br>conjugateGradient.SetInput(outputImageSource.GetOutput())<br>conjugateGradient.SetInput(1, proj.GetOutput())<br>conjugateGradient.SetInput(2, weightsSource.GetOutput())<br>conjugateGradient.SetGeometry(geometry)<br>conjugateGradient.SetGamma(2.0)<br>conjugateGradient.SetNumberOfIterations(30)<br>conjugateGradient.SetDisableDisplacedDetectorFilter(True)<br><div>conjugateGradient.Update()</div></div>
_______________________________________________<br>
Rtk-users mailing list<br>
<a href="mailto:Rtk-users@public.kitware.com" target="_blank">Rtk-users@public.kitware.com</a><br>
<a href="https://public.kitware.com/mailman/listinfo/rtk-users" rel="noreferrer" target="_blank">https://public.kitware.com/mailman/listinfo/rtk-users</a><br>
</blockquote></div>
</blockquote></div>