Hi Simon,
I am currently trying to create an alternate version of
rtkforwardprojection that takes an image pointer as an input and returns
the projection without using the GGO parser.
It is very similar to the original with some minor adjustments to determine
size and spacing. It currently works as intended when the forwardprojection
flag type is set to use the JosephForwardProjectionImageFilter, however it
returns a blank image with the correct spacing and dimensions when using
the CudaForwardProjectionImageFilter.
Cuda is installed on my machine and the original rtkforwardprojection with
cuda also works as intended. I have included a snippet of code and attached
the header file containing the function if anyone wants to help me out on
this one.
Thanks!
Solomon
....
rtk::ForwardProjectionImageFilter<OutputImageType,
OutputImageType>::Pointer forwardProjection;
switch(fp_flag)
{
case(1):
forwardProjection =
rtk::JosephForwardProjectionImageFilter<OutputImageType,
OutputImageType>::New();
break;
case(2):
forwardProjection =
rtk::RayCastInterpolatorForwardProjectionImageFilter<OutputImageType,
OutputImageType>::New();
break;
case(3): //This section does not work as intended, even though it almost
exactly the same as original (SetStepSize argument is different)
#ifdef RTK_USE_CUDA
forwardProjection = rtk::CudaForwardProjectionImageFilter<OutputImageType,
OutputImageType>::New();
dynamic_cast<rtk::CudaForwardProjectionImageFilter<OutputImageType,
OutputImageType>*>( forwardProjection.GetPointer()
)->SetStepSize(spacingOutput[2]);
#else
std::cerr << "The program has not been compiled with cuda option" <<
std::endl;
return EXIT_FAILURE;
#endif
break;
default:
throw std::exception("Unhandled --method value.");
/* std::cerr << "Unhandled --method value." << std::endl;
return EXIT_FAILURE;*/
}
forwardProjection->SetInput( constantImageSource->GetOutput() );
forwardProjection->SetInput(1, input_image );
forwardProjection->SetGeometry( geometryReader->GetOutputObject() );
projProbe.Start();
if(lowmem_flag)
{
TRY_AND_EXIT_ON_ITK_EXCEPTION( forwardProjection->Update() );
}
projProbe.Stop();
.....
return forwardProjection->GetOutput();
/*=========================================================================
*
* Copyright RTK Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
//Functionally very similar to rtkforwardprojections but uses image pointer as input/output.
//Input image should have origin shifted to -(dimension*spacing/2)-1. See ShiftOffset.h
#pragma once
#include <stdexcept>
//#include "rtkforwardprojections_ggo.h"
//#include "rtkGgoFunctions.h"
#include "rtkThreeDCircularProjectionGeometryXMLFile.h"
#include "rtkJosephForwardProjectionImageFilter.h"
#ifdef RTK_USE_CUDA
#include "rtkCudaForwardProjectionImageFilter.h"
#endif
#include "rtkRayCastInterpolatorForwardProjectionImageFilter.h"
//#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkTimeProbe.h>
#ifdef RTK_USE_CUDA
typedef itk::CudaImage< float, 3U > OutputImageType;
#else
typedef itk::Image< float, 3U > OutputImageType;
#endif
OutputImageType::Pointer ForwardProject3(
OutputImageType::Pointer input_image,
std::string fn_geometry,
double scale=1,
int fp_flag=1,
int verbose_flag=1,
int lowmem_flag=1)
//int main(int argc, char * argv[])
{
//GGO(rtkforwardprojections, args_info);
// Geometry
if(verbose_flag)
std::cout << "Reading geometry information from "
<< fn_geometry
<< "..."
<< std::flush;
rtk::ThreeDCircularProjectionGeometryXMLFileReader::Pointer geometryReader;
geometryReader = rtk::ThreeDCircularProjectionGeometryXMLFileReader::New();
geometryReader->SetFilename(fn_geometry);
TRY_AND_EXIT_ON_ITK_EXCEPTION( geometryReader->GenerateOutputInformation() )
if(verbose_flag)
std::cout << " done." << std::endl;
// Input reader
/* if(verbose_flag)
std::cout << "Reading input volume "
<< fn_mask_shifted
<< "..."
<< std::flush;
itk::TimeProbe readerProbe;
typedef itk::ImageFileReader< OutputImageType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( fn_mask_shifted );
readerProbe.Start();
TRY_AND_EXIT_ON_ITK_EXCEPTION( reader->Update() )
readerProbe.Stop();
if(1)
std::cout << " done in "
<< readerProbe.GetMean() << ' ' << readerProbe.GetUnit()
<< '.' << std::endl;*/
// Create a stack of empty projection images
typedef rtk::ConstantImageSource< OutputImageType > ConstantImageSourceType;
ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
constantImageSource ->SetInformationFromImage(input_image);
//rtk::SetConstantImageSourceFromGgo<ConstantImageSourceType, args_info_rtkforwardprojections>(constantImageSource, args_info);
// Adjust size according to input image and geometry
ConstantImageSourceType::SizeType sizeOutput;
ConstantImageSourceType::SpacingType spacingOutput;
ConstantImageSourceType::PointType originOutput;
sizeOutput[0] = constantImageSource->GetSize()[0];
sizeOutput[1] = constantImageSource->GetSize()[1];
sizeOutput[2] = constantImageSource->GetSize()[2];
spacingOutput[0] = constantImageSource->GetSpacing()[0];
spacingOutput[1] = constantImageSource->GetSpacing()[1];
spacingOutput[2] = constantImageSource->GetSpacing()[2];
originOutput[0] = -249;
originOutput[1] = -249;
originOutput[2] = -249;
int max_idx=0,diagonal;
for (int i = 1; i < 3; i++){
if (sizeOutput[max_idx]<sizeOutput[i]) {
max_idx = i;
}
}
diagonal = ceil(sqrt(pow(double(sizeOutput[max_idx]/scale),2)*3)*1.2*(spacingOutput[0]/spacingOutput[2]));
sizeOutput[0] = diagonal;
sizeOutput[1] = diagonal;
sizeOutput[2] = geometryReader->GetOutputObject()->GetGantryAngles().size();
constantImageSource->SetSize( sizeOutput );
spacingOutput[0] = spacingOutput[2];
spacingOutput[1] = spacingOutput[2];
spacingOutput = spacingOutput*scale;
constantImageSource->SetSpacing(spacingOutput);
constantImageSource->SetOrigin(originOutput);
// Create forward projection image filter
if(verbose_flag)
std::cout << "Projecting volume..." << std::flush;
itk::TimeProbe projProbe;
rtk::ForwardProjectionImageFilter<OutputImageType, OutputImageType>::Pointer forwardProjection;
switch(fp_flag)
{
case(1):
forwardProjection = rtk::JosephForwardProjectionImageFilter<OutputImageType, OutputImageType>::New();
break;
case(2):
forwardProjection = rtk::RayCastInterpolatorForwardProjectionImageFilter<OutputImageType, OutputImageType>::New();
break;
case(3):
#ifdef RTK_USE_CUDA
forwardProjection = rtk::CudaForwardProjectionImageFilter<OutputImageType, OutputImageType>::New();
dynamic_cast<rtk::CudaForwardProjectionImageFilter<OutputImageType, OutputImageType>*>( forwardProjection.GetPointer() )->SetStepSize(0.4);
#else
std::cerr << "The program has not been compiled with cuda option" << std::endl;
return EXIT_FAILURE;
#endif
break;
default:
throw std::exception("Unhandled --method value.");
/* std::cerr << "Unhandled --method value." << std::endl;
return EXIT_FAILURE;*/
}
forwardProjection->SetInput( constantImageSource->GetOutput() );
forwardProjection->SetInput(1, input_image );
forwardProjection->SetGeometry( geometryReader->GetOutputObject() );
projProbe.Start();
if(lowmem_flag)
{
TRY_AND_EXIT_ON_ITK_EXCEPTION( forwardProjection->Update() );
}
projProbe.Stop();
if(verbose_flag)
std::cout << " done in "
<< projProbe.GetMean() << ' ' << projProbe.GetUnit()
<< '.' << std::endl;
//Write
//if(verbose_flag)
// std::cout << "Writing... " << std::flush;
//itk::TimeProbe writeProbe;
//typedef itk::ImageFileWriter< OutputImageType > WriterType;
//WriterType::Pointer writer = WriterType::New();
//writer->SetFileName( "Y:/Outputs/projectiontest_3.mha" );
//writer->SetInput( forwardProjection->GetOutput() );
//if(lowmem_flag)
// {
// writer->SetNumberOfStreamDivisions(sizeOutput[2]);
// }
//writeProbe.Start();
//TRY_AND_EXIT_ON_ITK_EXCEPTION( writer->Update() );
//writeProbe.Stop();
//if(verbose_flag)
// std::cout << " done in "
// << writeProbe.GetMean() << ' ' << writeProbe.GetUnit()
// << '.' << std::endl;
return forwardProjection->GetOutput();
}
_______________________________________________
Rtk-users mailing list
[email protected]
http://public.kitware.com/mailman/listinfo/rtk-users
