RTK/Examples/FirstReconstruction: Difference between revisions

From Openrtk
Jump to navigation Jump to search
(Created page with "<div class="floatcenter">File:ITK_Examples_Baseline_ImageProcessing_TestConnectedComponentImageFilter.png</div> ==FirstReconstruction.cxx== <source lang="cpp"> // RTK inc...")
 
No edit summary
 
(2 intermediate revisions by the same user not shown)
Line 1: Line 1:
<div class="floatcenter">[[File:ITK_Examples_Baseline_ImageProcessing_TestConnectedComponentImageFilter.png]]</div>
<div class="floatcenter">[[File:FirstReconstructionTutorial.png|200px]]</div>


==FirstReconstruction.cxx==
==FirstReconstruction.cxx==
Line 69: Line 69:
   constantImageSource->SetConstant( 0. );
   constantImageSource->SetConstant( 0. );
    
    
  // Create the projector
  // Create the projector
   typedef rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType> REIType;
   typedef rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType> REIType;
   REIType::Pointer rei = REIType::New();
   REIType::Pointer rei = REIType::New();
   
  REIType::VectorType semiprincipalaxis, center;
  semiprincipalaxis.Fill(50.);
  center.Fill(0.);
   //Set GrayScale value, axes, center...
   //Set GrayScale value, axes, center...
   rei->SetMultiplicativeConstant(2.);
   rei->SetDensity(2.);
   rei->SetSemiPrincipalAxisX(50.);
   rei->SetAngle(0.);
  rei->SetSemiPrincipalAxisY(50.);
   rei->SetCenter(center);
  rei->SetSemiPrincipalAxisZ(50.);
   rei->SetAxis(semiprincipalaxis);
  rei->SetCenterX(0.);
  rei->SetCenterY(0.);
   rei->SetCenterZ(0.);
   rei->SetRotationAngle(0.);
   rei->SetGeometry( geometry );
   rei->SetGeometry( geometry );
   rei->SetInput( constantImageSource->GetOutput() );
   rei->SetInput( constantImageSource->GetOutput() );

Latest revision as of 10:17, 29 January 2014

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;
}