<div dir="ltr">Hi Simon,<br><i>I got a better reconstruction result with your advice last week,but there still are two mainly questions.on the one hand,
<span style="color:rgb(0,0,0)">as the picture1 showing,</span> the question <span style="color:rgb(0,0,0)">that the shape of the region corresponding to the reconstructed target changes from a rectangle to a diamond in the slice is still exist.on the other hand,there are much noise and artifact in the result. the projection images used were processed with dark image and flat image.The picture 2 is the reconstruction slice using astra toolbox with the same iamges.</span> </i> I want to get the same effect as Figure 2 by using the RTK , please help me to point out which code i need to fix in my code<i> .</i><br><i>.</i><img src="cid:ii_lxufyivz0" alt="picture1.png" width="190" height="190" style="margin-right: 0px;"> <img src="cid:ii_lxufyiwq1" alt="picture2.png" width="190" height="190" style="margin-right: 0px;"><div> picture1 picture2<br></div><div><br></div><div>The codes as follows:<br> int CT_ReconAndDetec(string projPath,string ReconPath,<br> int beginAng,int endAng,int rotateAngle,<br> int projNum, int DetecX,int DetecZ,<br> int nx, int ny, int nz,<br> float pixel_size){ <br>
using ImageType = itk::CudaImage<float, 3>;<br> using ProjImageType=itk::Image<unsigned short,3>;<br> using GeometryType = rtk::ThreeDCircularProjectionGeometry;<br> GeometryType::Pointer geometry = GeometryType::New();<br> unsigned int numberOfProjections = 1201; <br> double firstAngle = 0;<br> double angularArc = 360; <br><br> unsigned int sid = 66;<br> unsigned int sdd = 199;<br> for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)<br> {<br> double angle = firstAngle + noProj * angularArc / numberOfProjections;<br> double ProjOffsetX = 0.;<br> double ProjOffsetY = 0.;<br> double OutoffPlanAngle = 0;<br> double InPlanAngle = 0;<br> double SourceOffSetX = 0.;<br> double SourceOffSetY = 0.;<br> geometry->AddProjection(sid, sdd, angle, ProjOffsetX, ProjOffsetY,<br> OutoffPlanAngle, InPlanAngle, SourceOffSetX, SourceOffSetY);<br> } </div><div> <br> using NamesGeneratorType = itk::NumericSeriesFileNames;<br> NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();<br> nameGenerator->SetSeriesFormat(projPath+"/scan_%06d.tif");<br> nameGenerator->SetStartIndex(0);<br> nameGenerator->SetEndIndex(1200);<br> nameGenerator->SetIncrementIndex(1);<br> std::vector<std::string> fileNames = nameGenerator->GetFileNames();</div><div><br></div> using Reader3DType = rtk::ProjectionsReader<ProjImageType>;<br><div> Reader3DType::Pointer readerT = Reader3DType::New(); <br> using TIFFIOType = itk::TIFFImageIO;<br> TIFFIOType::Pointer tiffIO = TIFFIOType::New(); <br> <br> readerT->SetImageIO(tiffIO);<br> readerT->SetFileNames(fileNames);<br> try<br> {<br> readerT->Update();<br> }<br> catch (itk::ExceptionObject & error)<br> {<br> std::cerr << "ExceptionObject caught !" << std::endl;<br> std::cerr << error << std::endl;<br> return EXIT_FAILURE;<br> } <br> using CastFilterType = itk::CastImageFilter<ProjImageType, ImageType>;<br> CastFilterType::Pointer castFilter = CastFilterType::New();<br> castFilter->SetInput(readerT->GetOutput());<br> try<br> {<br> castFilter->Update();<br> }<br> catch (itk::ExceptionObject & error)<br> {<br> std::cerr << "Error casting the image: " << error << std::endl;<br> return EXIT_FAILURE;<br> }<br><br> ImageType::Pointer inputImage = castFilter->GetOutput(); <br> ImageType::PointType newOrigin;<br> newOrigin[0]= -DetecX / 2.0 * pixel_size; <br> newOrigin[1] = -DetecZ / 2.0 * pixel_size;<br> newOrigin[2] = 0.;<br> inputImage->SetOrigin(newOrigin);</div><div><br> ImageType::DirectionType newDirection;<br> newDirection.SetIdentity()<br> inputImage->SetDirection(newDirection);<br><br> ImageType::SpacingType newSpacing;<br> newSpacing[0] = pixel_size; <br> newSpacing[1] = pixel_size; <br> newSpacing[2] = pixel_size; <br> inputImage->SetSpacing(newSpacing); </div><div><br> using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;<br> ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();<br> ConstantImageSourceType::PointType origin;<br> ConstantImageSourceType::SpacingType spacing;<br> ConstantImageSourceType::SizeType sizeOutput; <br> sizeOutput[0] = nx;<br> sizeOutput[1] = ny;<br> sizeOutput[2] = nz;<br> spacing.Fill(pixel_size);<br> origin[0] = -(sizeOutput[0] / 2.0)*spacing[0]; <br> origin[1] = -(sizeOutput[1] / 2.0)*spacing[1]; <br> origin[2] = -(sizeOutput[2] / 2.0)*spacing[2]; <br> constantImageSource->SetOrigin(origin);<br> constantImageSource->SetSpacing(spacing);<br> constantImageSource->SetSize(sizeOutput);<br> constantImageSource->SetConstant(0.5); </div><div><br> using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter; <br> FDKGPUType::Pointer feldkamp = FDKGPUType::New();<br> feldkamp->SetInput(0, constantImageSource->GetOutput());<br> feldkamp->SetInput(1, inputImage); <br> feldkamp->SetGeometry(geometry);<br> feldkamp->GetRampFilter()->SetTruncationCorrection(0.5);<br> feldkamp->GetRampFilter()->SetHannCutFrequency(0.86); </div><div><br> using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;<br> FOVFilterType::Pointer fieldofview = FOVFilterType::New();<br> fieldofview->SetInput(0, feldkamp->GetOutput());<br> fieldofview->SetProjectionsStack(inputImage); // readerT->GetOutput()<br> fieldofview->SetGeometry(geometry); </div><div><br> using WriterType = itk::ImageFileWriter<ImageType>;<br> WriterType::Pointer writer = WriterType::New();<br> writer->SetFileName(ReconPath);<br> writer->SetInput(fieldofview->GetOutput());<br> try<br> {<br> writer->Update();<br> }<br> catch (itk::ExceptionObject & error)<br> {<br> std::cerr << error << std::endl;<br> }<br> std::cout << "Done!" << std::endl;<br> return 0; <br>}<br><br></div></div>