[Rtk-users] data loss in volume when reconstructing a anisotropic voxel

[BALADASTAGIRI ROOPANAGUDI]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230407/80f6d02e/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: samepixel.PNG
Type: image/png
Size: 146750 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230407/80f6d02e/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vaiablepixelspac.PNG
Type: image/png
Size: 201165 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230407/80f6d02e/attachment-0003.png>