RTK/Examples/FirstReconstruction

From Openrtk
Revision as of 10:17, 29 January 2014 by Julien.jomier (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
Error creating thumbnail: File missing

FirstReconstruction.cxx

// RTK includes
#include "rtkConfiguration.h"
#include <rtkFDKBackProjectionImageFilter.h>
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometry.h>
#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkDisplacedDetectorImageFilter.h>
#include <rtkParkerShortScanImageFilter.h>
#include <rtkFDKConeBeamReconstructionFilter.h>

// ITK includes
#include <itkImageFileWriter.h>
#include <itkStreamingImageFilter.h>

int main(int argc, char **argv)
{
  // Defines the image type
  typedef itk::Image< float, 3 > ImageType;

  // Defines the RTK geometry object
  typedef rtk::ThreeDCircularProjectionGeometry GeometryType;
  GeometryType::Pointer geometry = GeometryType::New();

  // Projection matrices
  unsigned int numberOfProjections = 360;
  unsigned int firstAngle = 0;
  unsigned int angularArc = 360;
  unsigned int sid = 600; // source to isocenter distance in mm
  unsigned int sdd = 1200; // source to detector distance in mm
  int isox = 0; // X coordinate on the projection image of isocenter
  int isoy = 0; // Y coordinate on the projection image of isocenter

  for(unsigned int noProj=0; noProj<numberOfProjections; noProj++)
    {
    double angle = (float)firstAngle + (float)noProj * angularArc / (float)numberOfProjections;
    geometry->AddProjection(sid,
                            sdd,
                            angle,
                            isox,
                            isoy);
    }

  // Create a stack of empty projection images
  typedef rtk::ConstantImageSource< ImageType > ConstantImageSourceType;
  ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
  ConstantImageSourceType::PointType origin;
  ConstantImageSourceType::SpacingType spacing;
  ConstantImageSourceType::SizeType sizeOutput;

  origin[0] = -127.;
  origin[1] = -127.;
  origin[2] = -127.;

  // Adjust size according to geometry
  sizeOutput[0] = 256;
  sizeOutput[1] = 256;
  sizeOutput[2] = numberOfProjections;

  spacing[0] = 1.;
  spacing[1] = 1.;
  spacing[2] = 1.;
  constantImageSource->SetOrigin( origin );
  constantImageSource->SetSpacing( spacing );
  constantImageSource->SetSize( sizeOutput );
  constantImageSource->SetConstant( 0. );
  
  // Create the projector
  typedef rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType> REIType;
  REIType::Pointer rei = REIType::New();
  REIType::VectorType semiprincipalaxis, center;
  semiprincipalaxis.Fill(50.);
  center.Fill(0.);
  //Set GrayScale value, axes, center...
  rei->SetDensity(2.);
  rei->SetAngle(0.);
  rei->SetCenter(center);
  rei->SetAxis(semiprincipalaxis);
  rei->SetGeometry( geometry );
  rei->SetInput( constantImageSource->GetOutput() );
  rei->Update();

  // Create reconstructed image
  ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();

  // Adjust size according to geometry
  sizeOutput[0] = 256;
  sizeOutput[1] = 256;
  sizeOutput[2] = 256;

  constantImageSource2->SetOrigin( origin );
  constantImageSource2->SetSpacing( spacing );
  constantImageSource2->SetSize( sizeOutput );
  constantImageSource2->SetConstant( 0. );

  // FDK reconstruction filtering
  typedef rtk::FDKConeBeamReconstructionFilter< ImageType > FDKCPUType;
  FDKCPUType::Pointer feldkamp = FDKCPUType::New();
  feldkamp->SetInput( 0, constantImageSource2->GetOutput() );
  feldkamp->SetInput( 1, rei->GetOutput() );
  feldkamp->SetGeometry( geometry );
  feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
  feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
  feldkamp->Update();

  // Writer
  typedef itk::ImageFileWriter<ImageType> WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( "output.mha" );
  writer->SetUseCompression(true);
  writer->SetInput( feldkamp->GetOutput() );
  writer->Update();
  
  return 0;
}