<div dir="ltr"><div>Hi,</div><div>If I'm following correctly, your reconstructed images both have 120 pixels in the Y direction. If you use .5 mm spacing, the length of your volume in the Y direction is 60 mm. If you use 256/120=2.1 mm, the length is 256 mm, as instructed. You need to keep the length (spacing_out[1]*size_out[1]) constant if you want to keep the same length in the Y direction.</div><div>Simon<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Apr 7, 2023 at 11:10 AM BALADASTAGIRI ROOPANAGUDI <<a href="mailto:dastagiri051428@gmail.com">dastagiri051428@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 all,</div><div><br></div><div>
When i tried to reconstruct a anisotropic voxel there is loss of data in
the top & bottom portion of volume and the same doesnt happen if i
try to reconstruct it with equal spacing out in all three direction.</div><div><br></div><div>Below is the code used for
anisotropic voxel for slice thickness of 1mm<br></div><div><br></div><div> size[0] = 1536; // size along X<br> size[1] = 1536; // size along Y<br> size[2] = angleLength; // size along Z<br> ImportFilterType::IndexType start;<br> start.Fill(0);<br> ImportFilterType::RegionType region;<br> region.SetIndex( start );<br> region.SetSize( size );<br><br> importFilter->SetRegion( region );<br><br> double spacing[ dimension ];<br> spacing[0] = 0.278; // along X direction<br> spacing[1] = 0.278; // along Y direction<br> spacing[2] = 1.0;<br> double origin[ dimension ];<br> origin[0] =-0.5* spacing[0]*(size[0]-1);// X coordinate<br> origin[1] = -0.5* spacing[1]*(size[1]-1); // Y coordinate<br> origin[2] =-0.5*spacing[2]*size[2]; // Z coordinate<br><br> importFilter->SetOrigin( origin );// along Z direction<br><br> importFilter->SetSpacing( spacing );<br><br> const unsigned int numberOfPixels = size[0] * size[1] * size[2];<br> const bool importImageFilterWillOwnTheBuffer = false; <br> importFilter->SetImportPointer( ExtBuffer, numberOfPixels,<br> importImageFilterWillOwnTheBuffer );<br> try{<br> importFilter->Update();<br> }<br> catch(itk::ExceptionObject e){<br> <br> throw(e);<br> return 0;<br> }<br><br> ProjectionsType::Pointer rval = importFilter->GetOutput();<br><br> itk::OrientImageFilter<ProjectionsType, ProjectionsType>::Pointer orienter =<br> itk::OrientImageFilter<ProjectionsType, ProjectionsType>::New();<br> orienter->SetInput(rval);<br><br> try<br> {<br><br> orienter->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br> orienter->GenerateOutputInformation();<br> rval = orienter->GetOutput();<br> typedef itk::ChangeInformationImageFilter< ProjectionsType> FilterType;<br> FilterType::Pointer filter = FilterType::New();<br> ProjectionsType::SpacingType pSpacing;<br> pSpacing[0] = 0.278;<br> pSpacing[1] = 0.278;<br> pSpacing[2] =1;<br> ProjectionsType::PointType pOrigin;<br> pOrigin[0] = -0.5* pSpacing[0]*(1536-1);<br> pOrigin[1] = -0.5* pSpacing[0]*(1536-1);<br> pOrigin[2] = -0.5* pSpacing[2]*angleLength;<br><br> filter->SetOutputOrigin(pOrigin);<br> ProjectionsType::DirectionType direction =importFilter->GetOutput()->GetDirection();<br> filter->SetOutputDirection(direction);<br> filter->SetOutputSpacing(pSpacing);<br> filter->ChangeAll();<br> filter->SetInput(rval);<br> try<br> {<br> filter->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br><br> typedef rtk::ThreeDCircularProjectionGeometry Geometry;<br> Geometry::Pointer geo=Geometry::New();<br> angleCount=mDetailsD2.at(8);<br> float angleValue = angles[0];<br> geo->AddProjection(mDetailsD2.at(0), mDetailsD2.at(1),angles[i],mDetailsD2.at(2),mDetailsD2.at(3),mDetailsD2.at(4),mDetailsD2.at(5),mDetailsD2.at(6),mDetailsD2.at(7));<br><br> try<br> {<br> geo->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br><br><br> typedef rtk::ThreeDCircularProjectionGeometryXMLFileWriter GeometryWriterType;<br> GeometryWriterType::Pointer geometryWriter = GeometryWriterType::New();<br> geometryWriter->SetFilename("D:\\geo_cirs.xml");<br> geometryWriter->SetObject(geo);<br> geometryWriter->WriteFile();<br> try<br> {<br> geo->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br><br> ConstantImageSourceType::Pointer outVol=ConstantImageSourceType::New();<br><br> //generate a constant image for reconstruction.<br> ConstantImageSourceType::PointType origin_out;<br> ConstantImageSourceType::SizeType size_out;<br> ConstantImageSourceType::SpacingType spacing_out;<br><br> size_out[0] = ctDIM1;//512<br> size_out[1] = ctDIM3;//120<br> size_out[2] = ctDIM2;//512<br><b><br> spacing_out[0] = 256.0/(float)ctDIM1;<br> spacing_out[1] = 256.0/(float)ctDIM3;//If spacing_out is made 0.5 in all direction it works fine<br> spacing_out[2] = 256.0/(float)ctDIM2;</b><br><br> origin_out[0] = -0.5*(size_out[0]-1)*(spacing_out[0]);<br> origin_out[1] = -0.5*(size_out[1]-1)*(spacing_out[1]);<br> origin_out[2] = -0.5*(size_out[2]-1)*(spacing_out[2]);<br> outVol->SetOrigin(origin_out);<br> outVol->SetSpacing(spacing_out);<br> outVol->SetSize(size_out);<br> outVol->SetConstant(0);<br> auto start5 = high_resolution_clock::now();<br><br> typedef rtk::BoellaardScatterCorrectionImageFilter<ProjectionsType,ProjectionsType> LagCorrectionImageFilterType1;<br> LagCorrectionImageFilterType1::Pointer LagCorrection = LagCorrectionImageFilterType1::New();<br> LagCorrection->SetInput(filter->GetOutput());<br> <br> LagCorrection->SetScatterToPrimaryRatio(1.2);//0.8<br><br> try<br> {<br> LagCorrection->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br><br> typedef rtk::I0EstimationProjectionFilter<ProjectionsType,ProjectionsType,3>IoFilterType;<br> IoFilterType::Pointer ioFilter=IoFilterType::New();<br> ioFilter->SetInput(LagCorrection->GetOutput());<br> try<br> {<br> ioFilter->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br> typedef rtk::LUTbasedVariableI0RawToAttenuationImageFilter<ProjectionsType, VolumeType> ConvertFilterType;<br> ConvertFilterType::Pointer convert = ConvertFilterType::New();<br> convert->SetInput(ioFilter->GetOutput());<br><br> try<br> {<br> convert->Update();<br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br> <br><br> using DDFOFFFOVType = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<VolumeType>;<br> DDFOFFFOVType::Pointer ddf;<br> ddf = DDFOFFFOVType::New();<br> ddf->SetInput(convert->GetOutput());<br> ddf->SetGeometry(geo);<br> PSSFType::Pointer pssf = PSSFType::New();<br> pssf->SetInput(ddf->GetOutput());<br><br> pssf->SetGeometry(geo);<br> pssf->InPlaceOff();<br> try<br> {<br><br> pssf->Update();<br><br> }<br> catch(itk::ExceptionObject e)<br> {<br> throw(e);<br> return 0;<br> }<br><br> FDKType::Pointer feldkamp=FDKType::New();<br> feldkamp->SetInput(0, outVol->GetOutput());<br> feldkamp->SetInput(1, pssf->GetOutput());<br> feldkamp->SetGeometry(geo);<br> feldkamp->GetRampFilter()->SetTruncationCorrection(1.0);<br> feldkamp->GetRampFilter()->SetHannCutFrequency(0.8);//high pass filter. <br> feldkamp->SetGPUEnabled(1); <br> feldkamp->Update();</div><div><br></div><div>
I am not sure where i am going wrong,<br><br><div>I have attached two images with same and variable pixelspacing for your references,</div>
</div><div><br></div><div><br></div><div>
<div>With Regards,</div><div>Dastagiri R</div>
</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>
</blockquote></div>