Hi,
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.
Simon
On Fri, Apr 7, 2023 at 11:10 AM BALADASTAGIRI ROOPANAGUDI <
dastagiri051...@gmail.com> wrote:
> Hi all,
>
> 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.
>
> Below is the code used for anisotropic voxel for slice thickness of 1mm
>
> size[0] = 1536; // size along X
> size[1] = 1536; // size along Y
> size[2] = angleLength; // size along Z
> ImportFilterType::IndexType start;
> start.Fill(0);
> ImportFilterType::RegionType region;
> region.SetIndex( start );
> region.SetSize( size );
>
> importFilter->SetRegion( region );
>
> double spacing[ dimension ];
> spacing[0] = 0.278; // along X direction
> spacing[1] = 0.278; // along Y direction
> spacing[2] = 1.0;
> double origin[ dimension ];
> origin[0] =-0.5* spacing[0]*(size[0]-1);// X coordinate
> origin[1] = -0.5* spacing[1]*(size[1]-1); // Y coordinate
> origin[2] =-0.5*spacing[2]*size[2]; // Z coordinate
>
> importFilter->SetOrigin( origin );// along Z direction
>
> importFilter->SetSpacing( spacing );
>
> const unsigned int numberOfPixels = size[0] * size[1] * size[2];
> const bool importImageFilterWillOwnTheBuffer = false;
> importFilter->SetImportPointer( ExtBuffer, numberOfPixels,
> importImageFilterWillOwnTheBuffer );
> try{
> importFilter->Update();
> }
> catch(itk::ExceptionObject e){
>
> throw(e);
> return 0;
> }
>
> ProjectionsType::Pointer rval = importFilter->GetOutput();
>
> itk::OrientImageFilter<ProjectionsType, ProjectionsType>::Pointer
> orienter =
> itk::OrientImageFilter<ProjectionsType,
> ProjectionsType>::New();
> orienter->SetInput(rval);
>
> try
> {
>
> orienter->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
> orienter->GenerateOutputInformation();
> rval = orienter->GetOutput();
> typedef itk::ChangeInformationImageFilter< ProjectionsType>
> FilterType;
> FilterType::Pointer filter = FilterType::New();
> ProjectionsType::SpacingType pSpacing;
> pSpacing[0] = 0.278;
> pSpacing[1] = 0.278;
> pSpacing[2] =1;
> ProjectionsType::PointType pOrigin;
> pOrigin[0] = -0.5* pSpacing[0]*(1536-1);
> pOrigin[1] = -0.5* pSpacing[0]*(1536-1);
> pOrigin[2] = -0.5* pSpacing[2]*angleLength;
>
> filter->SetOutputOrigin(pOrigin);
> ProjectionsType::DirectionType direction
> =importFilter->GetOutput()->GetDirection();
> filter->SetOutputDirection(direction);
> filter->SetOutputSpacing(pSpacing);
> filter->ChangeAll();
> filter->SetInput(rval);
> try
> {
> filter->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
>
> typedef rtk::ThreeDCircularProjectionGeometry Geometry;
> Geometry::Pointer geo=Geometry::New();
> angleCount=mDetailsD2.at(8);
> float angleValue = angles[0];
> 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));
>
> try
> {
> geo->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
>
>
> typedef rtk::ThreeDCircularProjectionGeometryXMLFileWriter
> GeometryWriterType;
> GeometryWriterType::Pointer geometryWriter = GeometryWriterType::New();
> geometryWriter->SetFilename("D:\\geo_cirs.xml");
> geometryWriter->SetObject(geo);
> geometryWriter->WriteFile();
> try
> {
> geo->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
>
> ConstantImageSourceType::Pointer outVol=ConstantImageSourceType::New();
>
> //generate a constant image for reconstruction.
> ConstantImageSourceType::PointType origin_out;
> ConstantImageSourceType::SizeType size_out;
> ConstantImageSourceType::SpacingType spacing_out;
>
> size_out[0] = ctDIM1;//512
> size_out[1] = ctDIM3;//120
> size_out[2] = ctDIM2;//512
>
>
>
> * spacing_out[0] = 256.0/(float)ctDIM1; spacing_out[1] =
> 256.0/(float)ctDIM3;//If spacing_out is made 0.5 in all direction it
> works fine spacing_out[2] = 256.0/(float)ctDIM2;*
>
> origin_out[0] = -0.5*(size_out[0]-1)*(spacing_out[0]);
> origin_out[1] = -0.5*(size_out[1]-1)*(spacing_out[1]);
> origin_out[2] = -0.5*(size_out[2]-1)*(spacing_out[2]);
> outVol->SetOrigin(origin_out);
> outVol->SetSpacing(spacing_out);
> outVol->SetSize(size_out);
> outVol->SetConstant(0);
> auto start5 = high_resolution_clock::now();
>
> typedef
> rtk::BoellaardScatterCorrectionImageFilter<ProjectionsType,ProjectionsType>
> LagCorrectionImageFilterType1;
> LagCorrectionImageFilterType1::Pointer LagCorrection =
> LagCorrectionImageFilterType1::New();
> LagCorrection->SetInput(filter->GetOutput());
>
> LagCorrection->SetScatterToPrimaryRatio(1.2);//0.8
>
> try
> {
> LagCorrection->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
>
> typedef
> rtk::I0EstimationProjectionFilter<ProjectionsType,ProjectionsType,3>IoFilterType;
> IoFilterType::Pointer ioFilter=IoFilterType::New();
> ioFilter->SetInput(LagCorrection->GetOutput());
> try
> {
> ioFilter->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
> typedef
> rtk::LUTbasedVariableI0RawToAttenuationImageFilter<ProjectionsType,
> VolumeType> ConvertFilterType;
> ConvertFilterType::Pointer convert = ConvertFilterType::New();
> convert->SetInput(ioFilter->GetOutput());
>
> try
> {
> convert->Update();
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
>
> using DDFOFFFOVType =
> rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<VolumeType>;
> DDFOFFFOVType::Pointer ddf;
> ddf = DDFOFFFOVType::New();
> ddf->SetInput(convert->GetOutput());
> ddf->SetGeometry(geo);
> PSSFType::Pointer pssf = PSSFType::New();
> pssf->SetInput(ddf->GetOutput());
>
> pssf->SetGeometry(geo);
> pssf->InPlaceOff();
> try
> {
>
> pssf->Update();
>
> }
> catch(itk::ExceptionObject e)
> {
> throw(e);
> return 0;
> }
>
> FDKType::Pointer feldkamp=FDKType::New();
> feldkamp->SetInput(0, outVol->GetOutput());
> feldkamp->SetInput(1, pssf->GetOutput());
> feldkamp->SetGeometry(geo);
> feldkamp->GetRampFilter()->SetTruncationCorrection(1.0);
> feldkamp->GetRampFilter()->SetHannCutFrequency(0.8);//high pass
> filter.
> feldkamp->SetGPUEnabled(1);
> feldkamp->Update();
>
> I am not sure where i am going wrong,
>
> I have attached two images with same and variable pixelspacing for your
> references,
>
>
> With Regards,
> Dastagiri R
> _______________________________________________
> Rtk-users mailing list
> rtk-us...@openrtk.org
> https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
>
_______________________________________________
Rtk-users mailing list
rtk-us...@openrtk.org
https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
