[ITK] [ITK-users] Multiresolution Registration error while trying to improve a journal paper

Dženan Zukić dzenanz at gmail.com
Mon Feb 8 09:45:07 EST 2016


Can you attache the source code as a file? When you pasted it into the
email, some long comments were wrapped and that made them non-comments.
There is over a 100 compiler error messages.

On Mon, Feb 8, 2016 at 3:03 AM, vishal <itkhelpacc at gmail.com> wrote:

> hey Dženan,
> this is the code below is a multiresolution version of the registration
> code
> in this
> https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration
> project... first when i tried to build and execute the code it shows the
> error "vector subscript is out of range".. then i tried to modified the
> github code keeping IntensityBased2D3DRegistration(the code below) as a
> reference it still shows another error called "access violation
> location...." ... kindly help me out...
> PLS NOTE:- i modified the code to take only one 2d image and one 3d volume
> inorder to perform registration. to run this code u will have to build the
> GetSiddonRayCastTracing.cxx code and obtain a projection. then build the
> code below and give the following arguments:
> Something.exe -v theProjectionFromtheStepAbove.dcm 0 3Dvolume.dcm
> where
> 0 is the projection angle.
> regards
> Vishal
>
> /*=========================================================================
>  *
>  *  Copyright Insight Software 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.
>  *
>
>
> *=========================================================================*/
>
> /*=========================================================================
>
>  This program implements an intensity based 2D-3D registration algorithm
> using the
>  SiddonJacobsRayCastInterpolateImageFunction and
> NormalizedCorrelationTwoImageToOneImageMetric
>  similarity measure
>
>  PowellOptimizer is used as the optimization method to avoid gradient
> calculation.
>  Euler3DTransform instead of CenteredEuler3DTransform is used to avoid the
> shift of the center.
>
>  When generating DRRs, the program attempts to simulate a 2D imaging system
> attached to a linac
>  for radiation therapy. The program registers two 2D portal images with
> their corresponding DRR
>  images generated from the 3D dataset. The portal images may be acquired at
> any arbitrary projection
>  angles. The similarity measure is the summation of the measures calculated
> each projection.
>  Multiresolution strategy was not implemented.
>
>  This program was modified from the ITK
> application--IntensityBased2D3DRegistration.cxx
>
> =========================================================================*/
> //#include "itkTwoProjectionImageRegistrationMethod.h"
>
> // The transformation used is a rigid 3D Euler transform with the
> // provision of a center of rotation which defaults to the center of
> // the 3D volume. In practice the center of the particular feature of
> // interest in the 3D volume should be used.
> #include <iostream>
>                                         //WRITTEN USING
> IntesityBased2D3dRegistration.cxx &
> MultiResImageRegistration3.cxx
> #include "itkEuler3DTransform.h"
>
> // We have chosen to implement the registration using the normalized
> coorelation
> // similarity measure.
>
> //#include "itkGradientDifferenceTwoImageToOneImageMetric.h"
> //#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
>
> // This is an intensity based registration algorithm so ray casting is
> // used to project the 3D volume onto pixels in the target 2D image.
> #include "itkSiddonJacobsRayCastInterpolateImageFunction.h"
> #include "itkMultiResolutionImageRegistrationMethod.h"
> // Finally the Powell optimizer is used to avoid requiring gradient
> information.
> //#include "itkPowellOptimizer.h"
> #include "itkRegularStepGradientDescentOptimizer.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkRescaleIntensityImageFilter.h"
> #include "itkFlipImageFilter.h"
>
> #include "itkCommand.h"
> #include "itkTimeProbesCollectorBase.h"
>
> //MODIFIED HEADERS
> //#include "itkImageRegistrationMethodv4.h"
> //#include <itkImageRegistrationMethod.h>
> #include <itkNormalizedCorrelationImageToImageMetric.h>
> //#include "itkMattesMutualInformationImageToImageMetric.h"
> //#include <itkObjectToObjectOptimizerBase.h>
> //#include "itkEuler3DTransform.h"
> //#include "itkPowellOptimizer.h"
> //
> //#include "itkImageFileReader.h"
> //#include "itkImageFileWriter.h"
> //
> //#include "itkResampleImageFilter.h"
> //#include "itkCastImageFilter.h"
> //#include "itkRescaleIntensityImageFilter.h"
> //#include "itkFlipImageFilter.h"
> //
> //#include "itkCommand.h"
> //#include "itkTimeProbesCollectorBase.h"
> //#include "itkSiddonJacobsRayCastInterpolateImageFunction.h"
> // First we define the command class to allow us to monitor the
> registration.
>
> class CommandIterationUpdate : public itk::Command
> {
> public:
>   typedef  CommandIterationUpdate Self;
>   typedef  itk::Command           Superclass;
>   typedef itk::SmartPointer<Self> Pointer;
>   itkNewMacro( Self );
>
> protected:
>   CommandIterationUpdate() {};
>
> public:
>   typedef itk::RegularStepGradientDescentOptimizer  OptimizerType;
>   typedef const OptimizerType * OptimizerPointer;
>
>   void Execute(itk::Object *caller, const itk::EventObject & event)
>   {
>     Execute( (const itk::Object *)caller, event);
>   }
>
>   void Execute(const itk::Object *object, const itk::EventObject & event)
>   {
>     OptimizerPointer optimizer =
>       dynamic_cast< OptimizerPointer >( object );
>     if( typeid( event ) != typeid( itk::IterationEvent ) )
>       {
>       return;
>       }
>     std::cout << "Iteration: " << optimizer->GetCurrentIteration() <<
> std::endl;
>     std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
>     std::cout << "Position: " << optimizer->GetCurrentPosition() <<
> std::endl;
>         //std::cout<<"GETCURRENT LINE VALUE
> "<<optimizer->GetCurrentLineIteration()<< std::endl;
>   }
> };
>
>
> // Software Guide : BeginLatex
> //
> // A second \emph{Command/Observer} is used to reduce the minimum
> // registration step length on each occasion that the resolution of
> // the multi-scale registration is increased (see
> // \doxygen{MultiResImageRegistration1} for more info).
> //
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> template <typename TRegistration>
> class RegistrationInterfaceCommand : public itk::Command
> {
> public:
>   typedef  RegistrationInterfaceCommand   Self;
>   typedef  itk::Command                   Superclass;
>   typedef  itk::SmartPointer<Self>        Pointer;
>   itkNewMacro( Self );
>
> protected:
>   RegistrationInterfaceCommand() {};
>
> public:
>   typedef   TRegistration                              RegistrationType;
>   typedef   RegistrationType *                         RegistrationPointer;
>   typedef   itk::RegularStepGradientDescentOptimizer   OptimizerType;
>   typedef   OptimizerType *                            OptimizerPointer;
>
>   void Execute(itk::Object * object, const itk::EventObject & event)
> ITK_OVERRIDE
>   {
>    /* if( typeid( event ) != typeid( itk::IterationEvent ) )
>       {
>       return;
>       }*/ //old
>            if( !(itk::IterationEvent().CheckEvent( &event )) )
>       {
>       return;
>       }
>
>    // RegistrationPointer registration =
>           //                  dynamic_cast<RegistrationPointer>( object );
>
>     RegistrationPointer registration = static_cast<RegistrationPointer>(
> object );
>
>     // If this is the first resolution level we assume that the
>     // maximum step length (representing the first step size) and the
>     // minimum step length (representing the convergence criterion)
>     // have already been set.  At each subsequent resolution level, we
>     // will reduce the minimum step length by a factor of four in order
>     // to allow the optimizer to focus on progressively smaller
>     // regions. The maximum step length is set to the current step
>     // length. In this way, when the optimizer is reinitialized at the
>     // beginning of the registration process for the next level, the
>     // step length will simply start with the last value used for the
>     // previous level. This will guarantee the continuity of the path
>     // taken by the optimizer through the parameter space.
>           if(registration == ITK_NULLPTR)
>       {
>       return;
>       }
>     /*OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >(
>                        registration->GetOptimizer() );*/
>         OptimizerPointer optimizer = static_cast< OptimizerPointer
> >(registration->GetModifiableOptimizer() );
>
>         std::cout << "-------------------------------------" << std::endl;
>     std::cout << "MultiResolution Level : "
>               << registration->GetCurrentLevel()  << std::endl;
>     std::cout << std::endl;
>
>         //if ( registration->GetCurrentLevel() != 0 ) //old
>  //     {
>  //     optimizer->SetMaximumStepLength( optimizer->GetCurrentStepLength()
> );
>  //     optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength()
> /
> 4.0 );
>  //     }
>          if ( registration->GetCurrentLevel() == 0 )
>       {
>       optimizer->SetMaximumStepLength( 16.00 );
>       optimizer->SetMinimumStepLength( 0.01 );
>       }
>     else
>       {
>       optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() /
> 4.0 );
>       optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() /
> 10.0 );
>       } //mulit3
>     optimizer->Print( std::cout );
>   }
>
>  /* void Execute(const itk::Object * , const itk::EventObject & )
>     { return; }*/ //old
>         void Execute(const itk::Object * , const itk::EventObject & )
> ITK_OVERRIDE
>     { return; } //multi3
> };
> // Software Guide : EndCodeSnippet
>
>
>
>
> void exe_usage()
> {
>   std::cerr << "\n";
>   std::cerr << "Usage: TwoProjection2D3DRegistration <options> Image2D1
> ProjAngle1 Image2D2 ProjAngle2 Volume3D\n";
>   std::cerr << "       Registers two 2D images to a 3D volume. \n\n";
>   std::cerr << "   where <options> is one or more of the following:\n\n";
>   std::cerr << "       <-h>                     Display (this) usage
> information\n";
>   std::cerr << "       <-v>                     Verbose output [default:
> no]\n";
>   std::cerr << "       <-dbg>                   Debugging output [default:
> no]\n";
>  // old
>   std::cerr << "       <-n int>                 The number of scales to
> apply [default: 2]\n";
>   std::cerr << "       <-maxScale int>          The scale factor
> corresponding to max resolution [default: 1]\n";
>   std::cerr << "       <-step float float>      Maximum and minimum step
> sizes [default: 4 and 0.01]\n";
>  //end old
>   std::cerr << "       <-scd float>             Source to isocenter
> distance
> [default: 1000mm]\n";
>   std::cerr << "       <-t float float float>   CT volume translation in x,
> y, and z direction in mm \n";
>   std::cerr << "       <-rx float>              CT volume rotation about x
> axis in degrees \n";
>   std::cerr << "       <-ry float>              CT volume rotation about y
> axis in degrees \n";
>   std::cerr << "       <-rz float>              CT volume rotation about z
> axis in degrees \n";
>   std::cerr << "       <-2dcx float float float float>    Central axis
> positions of the 2D images in continuous indices \n";
>   std::cerr << "       <-res float float float float>     Pixel spacing of
> projection images in the isocenter plane [default: 1x1 mm]  \n";
>   std::cerr << "       <-iso float float float> Isocenter location in voxel
> in indices (center of rotation and projection center)\n";
>   std::cerr << "       <-threshold float>       Intensity threshold below
> which are ignore [default: 0]\n";
>   std::cerr << "       <-o file>                Output image filename\n\n";
>   std::cerr << "                                by  Jian Wu\n";
>   std::cerr << "                                eewujian at hotmail.com\n";
>   std::cerr << "                                (Univeristy of
> Florida)\n\n";
>   exit(EXIT_FAILURE);
> }
>
>
> int main( int argc, char *argv[] )
> {
>   char *fileImage2D1 = NULL;
>   double projAngle1 = -999;
>  /* char *fileImage2D2 = NULL;
>   double projAngle2 = -999;*/
>   char *fileVolume3D = NULL;
>   // Default output file names
>   const char *fileOutput1 = "Image2D1_Registered.tif";
>   //const char *fileOutput2 = "Image2D2_Registered.tif";
>
>   bool ok;
>   bool verbose = false;
>   bool debug = false;
>   bool customized_iso = false;
>   bool customized_2DCX = false; // Flag for customized 2D image central
> axis
> positions
>   bool customized_2DRES = false; // Flag for customized 2D image pixel
> spacing
>
>   unsigned int nScales = 2;     //OLD PROG
>   int maxScale = 1;                     //OLD PROG
>
>   double rx = 0.;
>   double ry = 0.;
>   double rz = 0.;
>
>   double tx = 0.;
>   double ty = 0.;
>   double tz = 0.;
>
>   double cx = 0.;
>   double cy = 0.;
>   double cz = 0.;
>
>   double scd = 1000.; // Source to isocenter distance
>
>
>   double maxStepSize = 4.;      //OLD PROG
>   double minStepSize = 1.;  //OLD PROG
>
>   double image1centerX = 0.0;
>   double image1centerY = 0.0;
>   double image2centerX = 0.0;
>   double image2centerY = 0.0;
>
>   double image1resX = 1.0;
>   double image1resY = 1.0;
>   double image2resX = 1.0;
>   double image2resY = 1.0;
>
>   double threshold = 0.0;
>
>   // Parse command line parameters
>
>   //if (argc <= 5) //ORIGINAL
>     if (argc <= 4)
>                 exe_usage();
>
>   while (argc > 1)
>     {
>     ok = false;
>
>     if ((ok == false) && (strcmp(argv[1], "-h") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       exe_usage();
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-v") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       verbose = true;
>           std::cout<<"INSIDE -v"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-dbg") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       debug = true;
>             std::cout<<"INSIDE -debug"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-scd") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       scd = atof(argv[1]);
>       argc--; argv++;
>             std::cout<<"INSIDE -scd"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-t") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       tx=atof(argv[1]);
>       argc--; argv++;
>       ty=atof(argv[1]);
>       argc--; argv++;
>       tz=atof(argv[1]);
>       argc--; argv++;
>             std::cout<<"INSIDE -t"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-rx") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       rx=atof(argv[1]);
>       argc--; argv++;
>             std::cout<<"INSIDE -rx"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-ry") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       ry=atof(argv[1]);
>       argc--; argv++;
>             std::cout<<"INSIDE -ry"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-rz") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       rz=atof(argv[1]);
>           //std::cout<<"RZ
> INTIALIZALIZES"<<std::endl;
>       argc--; argv++;
>             std::cout<<"INSIDE -rz"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-2dcx") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       image1centerX = atof(argv[1]);
>       argc--; argv++;
>       image1centerY = atof(argv[1]);
>       argc--; argv++;
>       image2centerX = atof(argv[1]);
>       argc--; argv++;
>       image2centerY = atof(argv[1]);
>       argc--; argv++;
>       customized_2DCX = true;
>             std::cout<<"INSIDE -2dcx"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-res") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       image1resX = atof(argv[1]);
>       argc--; argv++;
>       image1resY = atof(argv[1]);
>       argc--; argv++;
>       image2resX = atof(argv[1]);
>       argc--; argv++;
>       image2resY = atof(argv[1]);
>       argc--; argv++;
>       customized_2DRES = true;
>             std::cout<<"INSIDE -res"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-iso") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       cx=atof(argv[1]);
>       argc--; argv++;
>       cy=atof(argv[1]);
>       argc--; argv++;
>       cz=atof(argv[1]);
>       argc--; argv++;
>       customized_iso = true;
>             std::cout<<"INSIDE -iso"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-threshold")
> ==
> 0))
>       {
>       argc--; argv++;
>       ok = true;
>       threshold=atof(argv[1]);
>       argc--; argv++;
>             std::cout<<"INSIDE
> -threshold"<<std::endl;
>       }
>
>     if ((ok == false) && (strcmp(argv[1], "-o") == 0))
>       {
>       argc--; argv++;
>       ok = true;
>       fileOutput1 = argv[1];
>       argc--; argv++;
>           std::cout<<"INSIDE -o"<<std::endl;
>       /*fileOutput2 = argv[1];
>       argc--; argv++;*/
>       }
>
>
>     if (ok == false)
>       {
>
>       if (fileImage2D1 == NULL)
>         {
>         fileImage2D1 = argv[1];
>         argc--;
>         argv++;
>                   std::cout<<"INSIDE INTIALIZING FILE NAME
> fileImage2D1"<<std::endl;
>         }
>
>       if (projAngle1 == -999)
>         {
>         projAngle1 = atof(argv[1]);
>         argc--;
>         argv++;
>                  std::cout<<"INSIDE INTIALIZING FILE NAME
> projAngle1"<<std::endl;
>         }
>
>       //if (fileImage2D2 == NULL)
>       //  {
>       //  fileImage2D2 = argv[1];
>       //  argc--;
>       //  argv++;
>       //  }
>
>       //if (projAngle2 == -999)
>       //  {
>       //  projAngle2 = atof(argv[1]);
>       //  argc--;
>       //  argv++;
>       //  }
>
>       else if (fileVolume3D == NULL)
>         {
>         fileVolume3D = argv[1];
>         argc--;
>         argv++;
>                 std::cout<<"INSIDE INTIALIZING FILE NAME
> FileVolume3D"<<std::endl;
>         }
>
>       else
>         {
>         std::cerr << "ERROR: Cannot parse argument "
> << argv[1] << std::endl;
>         exe_usage();
>         }
>       }
>     }
>
>   if (verbose)
>     {
>     if (fileImage2D1)  std::cout << "Input 2D image 1: "
> << fileImage2D1  << std::endl;
>     if (fileImage2D1)  std::cout << "Projection Angle 1: "
> << projAngle1  << std::endl;
>     //if (fileImage2D2)  std::cout << "Input 2D image 2: "
> << fileImage2D2  << std::endl;
>     //if (fileImage2D2)  std::cout << "Projection Angle 2:
> "
> << projAngle2  << std::endl;
>     if (fileVolume3D) std::cout << "Input 3D image: "
> << fileVolume3D << std::endl;
>     if (fileOutput1)   std::cout << "Output image 1: "
> << fileOutput1   << std::endl;
>     //if (fileOutput2)   std::cout << "Output image 2: "
> << fileOutput2   << std::endl;
>     }
>
>
>   // We begin the program proper by defining the 2D and 3D images. The
>   // {TwoProjectionImageRegistrationMethod} requires that both
>   // images have the same dimension so the 2D image is given
>   // dimension 3 and the size of the {z} dimension is set to unity.
>   std::cout<<"OUT OF THE INLIAZING LOOPS PROGRAM
> STARTS!!"<<std::endl;
>   const    unsigned int    Dimension = 3;
>   typedef  float           InternalPixelType;
>   typedef  short           PixelType3D;
>
>   typedef itk::Image< PixelType3D, Dimension > ImageType3D;
>
>   typedef unsigned char                            OutputPixelType;
>   typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
>
>   // The following lines define each of the components used in the
>   // registration: The transform, optimizer, metric, interpolator and
>   // the registration method itself.
>
>   typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
>
>   typedef itk::Euler3DTransform< double >  TransformType;
>
>   typedef itk::RegularStepGradientDescentOptimizer      OptimizerType;
>
>   //typedef itk::GradientDifferenceTwoImageToOneImageMetric<
>  /* typedef itk::NormalizedCorrelationTwoImageToOneImageMetric<
>     InternalImageType,
>     InternalImageType >    MetricType;*/
>    typedef itk::NormalizedCorrelationImageToImageMetric<
>     InternalImageType,
>     InternalImageType >    MetricType;
>
>   typedef itk::SiddonJacobsRayCastInterpolateImageFunction<
>     InternalImageType,
>     double > InterpolatorType;
>
>
>  /* typedef itk::TwoProjectionImageRegistrationMethod<
>     InternalImageType,
>     InternalImageType >   RegistrationType;*/
>   typedef itk::MultiResolutionImageRegistrationMethod<
>                                     InternalImageType,
>                                     InternalImageType >
> RegistrationType;
>
>   // Software Guide : BeginCodeSnippet  //THESE ARE THE CHANGES MADE
>   typedef itk::MultiResolutionPyramidImageFilter<
>                                     InternalImageType,
>                                     InternalImageType > ImagePyramidType2D;
>   typedef itk::MultiResolutionPyramidImageFilter<
>                                     InternalImageType,
>                                     InternalImageType > ImagePyramidType3D;
>   // Software Guide : EndCodeSnippet
>
>
>   // Each of the registration components are instantiated in the
>   // usual way...
>
>     ImagePyramidType2D::Pointer imagePyramid2D =
>       ImagePyramidType2D::New();
>   ImagePyramidType3D::Pointer imagePyramid3D =
>       ImagePyramidType3D::New();
>
>
>   std::cout<<"TYPES INTIALIZED"<<std::endl;
>   MetricType::Pointer         metric        = MetricType::New();
>   TransformType::Pointer      transform     = TransformType::New();
>   OptimizerType::Pointer      optimizer     = OptimizerType::New();
>   InterpolatorType::Pointer   interpolator1  = InterpolatorType::New();
>   //InterpolatorType::Pointer   interpolator2  = InterpolatorType::New();
>   RegistrationType::Pointer   registration  = RegistrationType::New();
>
>  //metric->ComputeGradientOff(); //ORIGINAL
>   metric->ComputeGradientOn();
>   metric->SetSubtractMean( true );
>
>   // and passed to the registration method:
>
>   registration->SetMetric(    metric        );
>   registration->SetOptimizer(     optimizer     );
>   registration->SetTransform(     transform     );
>   registration->SetInterpolator(  interpolator1  );
>  // registration->SetInterpolator2(  interpolator2  );
>    registration->SetFixedImagePyramid( imagePyramid2D );
>   registration->SetMovingImagePyramid( imagePyramid3D ); //changes made
> from
> 2d/3d
>   std::cout<<"REGISTRATION METHODS SET"<<std::endl;
>   if (debug)
>     {
>     metric->DebugOn();
>     //transform->DebugOn();
>     //optimizer->DebugOn();
>     interpolator1->DebugOn();
>    // interpolator2->DebugOn();
>     //registration->DebugOn();
>     }
>
>
>   //  The 2- and 3-D images are read from files,
>
>   typedef itk::ImageFileReader< InternalImageType > ImageReaderType2D;
>   typedef itk::ImageFileReader< ImageType3D >       ImageReaderType3D;
>
>   ImageReaderType2D::Pointer imageReader2D1 = ImageReaderType2D::New();
>  // ImageReaderType2D::Pointer imageReader2D2 = ImageReaderType2D::New();
>   ImageReaderType3D::Pointer imageReader3D = ImageReaderType3D::New();
>
>   imageReader2D1->SetFileName( fileImage2D1 );
>   //imageReader2D2->SetFileName( fileImage2D2 );
>   imageReader3D->SetFileName( fileVolume3D );
>   imageReader3D->Update();
>   std::cout<<"2d 3d volume filename SET"<<std::endl;
>   ImageType3D::Pointer image3DIn = imageReader3D->GetOutput();
>
>   // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the
> origin of the CT image
>   // to be (0,0,0). Because we align the CT isocenter with the central
> axis,
> the projection
>   // geometry is fully defined. The origin of the CT image becomes
> irrelavent.
>   ImageType3D::PointType image3DOrigin;
>   image3DOrigin[0] = 0.0;
>   image3DOrigin[1] = 0.0;
>   image3DOrigin[2] = 0.0;
>   image3DIn->SetOrigin(image3DOrigin);
>   std::cout<<"VOLUME ORIGIN SET"<<std::endl;
>   InternalImageType::Pointer image_tmp1 = imageReader2D1->GetOutput();
>   //InternalImageType::Pointer image_tmp2 = imageReader2D2->GetOutput();
>
>   imageReader2D1->Update();
>  // imageReader2D2->Update();
>
>   if (customized_2DRES)
>     {
>     InternalImageType::SpacingType spacing;
>     spacing[0] = image1resX;
>     spacing[1] = image1resY;
>     spacing[2] = 1.0;
>     image_tmp1->SetSpacing( spacing );
>
>   //  spacing[0] = image2resX;
>   //  spacing[1] = image2resY;
>   //  image_tmp2->SetSpacing( spacing );
>
>     }
>   // The input 2D images were loaded as 3D images. They were considered
>   // as a single slice from a 3D volume. By default, images stored on the
>   // disk are treated as if they have RAI orientation. After view point
>   // transformation, the order of 2D image pixel reading is equivalent to
>   // from inferior to superior. This is contradictory to the traditional
>   // 2D x-ray image storage, in which a typical 2D image reader reads and
>   // writes images from superior to inferior. Thus the loaded 2D DICOM
>   // images should be flipped in y-direction. This was done by using a.
>   // FilpImageFilter.
>   typedef itk::FlipImageFilter< InternalImageType > FlipFilterType;
>   FlipFilterType::Pointer flipFilter1 = FlipFilterType::New();
>  // FlipFilterType::Pointer flipFilter2 = FlipFilterType::New();
>
>   typedef FlipFilterType::FlipAxesArrayType FlipAxesArrayType;
>   FlipAxesArrayType flipArray;
>   flipArray[0] = 0;
>   flipArray[1] = 1;
>   flipArray[2] = 0;
>
>   flipFilter1->SetFlipAxes( flipArray );
>  // flipFilter2->SetFlipAxes( flipArray );
>
>   flipFilter1->SetInput( imageReader2D1->GetOutput() );
>   //flipFilter2->SetInput( imageReader2D2->GetOutput() );
>
>   // The input 2D images may have 16 bits. We rescale the pixel value to
> between 0-255.
>   typedef itk::RescaleIntensityImageFilter<
>     InternalImageType, InternalImageType > Input2DRescaleFilterType;
>
>   Input2DRescaleFilterType::Pointer rescaler2D1 =
> Input2DRescaleFilterType::New();
>   rescaler2D1->SetOutputMinimum(   0 );
>   rescaler2D1->SetOutputMaximum( 255 );
>   rescaler2D1->SetInput( flipFilter1->GetOutput() );
>   std::cout<<"flipFilter1->GetOutput()  SET"<<std::endl;
>   /*Input2DRescaleFilterType::Pointer rescaler2D2 =
> Input2DRescaleFilterType::New();
>   rescaler2D2->SetOutputMinimum(   0 );
>   rescaler2D2->SetOutputMaximum( 255 );
>   rescaler2D2->SetInput( flipFilter2->GetOutput() );*/
>
>
>   //  The 3D CT dataset is casted to the internal image type using
>   //  {CastImageFilters}.
>
>   typedef itk::CastImageFilter<
>     ImageType3D, InternalImageType > CastFilterType3D;
>
>   CastFilterType3D::Pointer caster3D = CastFilterType3D::New();
>   caster3D->SetInput( image3DIn );
>
>   rescaler2D1->Update();
>   //rescaler2D2->Update();
>   caster3D->Update();
>     std::cout<<"caster3D->GetOutput()  SET"<<std::endl;
>
>   registration->SetFixedImage(  rescaler2D1->GetOutput() );
>   //registration->SetFixedImage2(  rescaler2D2->GetOutput() );
>   registration->SetMovingImage( caster3D->GetOutput() );
>
>     registration->SetFixedImageRegion(
>                 rescaler2D1->GetOutput()->GetBufferedRegion() ); //changes
> made from 2d/3d
>   // Initialise the transform
>   // ~~~~~~~~~~~~~~~~~~~~~~~~
>
>   // Set the order of the computation. Default ZXY
>   transform->SetComputeZYX(true);
>
>
>   // The transform is initialised with the translation [tx,ty,tz] and
>   // rotation [rx,ry,rz] specified on the command line
>
>   TransformType::OutputVectorType translation;
>
>   translation[0] = tx;
>   translation[1] = ty;
>   translation[2] = tz;
>
>   transform->SetTranslation(translation);
>     std::cout<<"transform->SetTranslation(translation);  SET"<<std::endl;
>   // constant for converting degrees to radians
>   const double dtr = ( atan(1.0) * 4.0 ) / 180.0;
>   transform->SetRotation(dtr*rx, dtr*ry, dtr*rz);
>    std::cout<<"transform->SetRotation(dtr*rx, dtr*ry, dtr*rz)
> SET"<<std::endl;
>   // The centre of rotation is set by default to the centre of the 3D
>   // volume but can be offset from this position using a command
>   // line specified translation [cx,cy,cz]
>
>   ImageType3D::PointType origin3D = image3DIn->GetOrigin();
>   const itk::Vector<double, 3> resolution3D = image3DIn->GetSpacing();
>
>   typedef ImageType3D::RegionType      ImageRegionType3D;
>   typedef ImageRegionType3D::SizeType  SizeType3D;
>
>   ImageRegionType3D region3D = caster3D->GetOutput()->GetBufferedRegion();
>   SizeType3D        size3D   = region3D.GetSize();
>
>   TransformType::InputPointType isocenter;
>   if (customized_iso)
>     {
>     // Isocenter location given by the user.
>     isocenter[0] = origin3D[0] + resolution3D[0] * cx;
>     isocenter[1] = origin3D[1] + resolution3D[1] * cy;
>     isocenter[2] = origin3D[2] + resolution3D[2] * cz;
>     }
>   else
>     {
>     // Set the center of the image as the isocenter.
>     isocenter[0] = origin3D[0] + resolution3D[0] * static_cast<double>(
> size3D[0] ) / 2.0;
>     isocenter[1] = origin3D[1] + resolution3D[1] * static_cast<double>(
> size3D[1] ) / 2.0;
>     isocenter[2] = origin3D[2] + resolution3D[2] * static_cast<double>(
> size3D[2] ) / 2.0;
>     }
>
>   transform->SetCenter(isocenter);
>
>
>   if (verbose)
>     {
>     std::cout << "3D image size: "
>               << size3D[0] << ", " << size3D[1] << ", " << size3D[2] <<
> std::endl
>               << "   resolution: "
>               << resolution3D[0] << ", " << resolution3D[1] << ", " <<
> resolution3D[2] << std::endl
>               << "Transform: " << transform << std::endl;
>     }
>
>
>   // Set the origin of the 2D image
>   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>   // For correct (perspective) projection of the 3D volume, the 2D
>   // image needs to be placed at a certain distance (the source-to-
>   // isocenter distance {scd} ) from the focal point, and the normal
>   // from the imaging plane to the focal point needs to be specified.
>   //
>   // By default, the imaging plane normal is set by default to the
>   // center of the 2D image but may be modified from this using the
>   // command line parameters [image1centerX, image1centerY,
>   // image2centerX, image2centerY].
>
>   double origin2D1[ Dimension ];
>  // double origin2D2[ Dimension ];
>
>   // Note: Two 2D images may have different image sizes and pixel
> dimensions, although
>   // scd are the same.
>
>   const itk::Vector<double, 3> resolution2D1 =
> imageReader2D1->GetOutput()->GetSpacing();
>   //const itk::Vector<double, 3> resolution2D2 =
> imageReader2D2->GetOutput()->GetSpacing();
>
>   typedef InternalImageType::RegionType ImageRegionType2D;
>   typedef ImageRegionType2D::SizeType   SizeType2D;
>
> //  ImageRegionType2D region2D1 =
> rescaler2D1->GetOutput()->GetBufferedRegion(); //ORIGINAL 2d/3d
>   ImageRegionType2D region2D1 =
> rescaler2D1->GetOutput()->GetLargestPossibleRegion();
>   //  ImageRegionType2D region2D2 =
> rescaler2D2->GetOutput()->GetBufferedRegion();
>   SizeType2D        size2D1   = region2D1.GetSize();
>  // SizeType2D        size2D2   = region2D2.GetSize();
>
>   if (!customized_2DCX)
>     { // Central axis positions are not given by the user. Use the image
> centers
>     // as the central axis position.
>     image1centerX = ((double) size2D1[0] - 1.)/2.;
>     image1centerY = ((double) size2D1[1] - 1.)/2.;
>     //image2centerX = ((double) size2D2[0] - 1.)/2.;
>    // image2centerY = ((double) size2D2[1] - 1.)/2.;
>     }
>
>   // 2D Image 1
>   origin2D1[0] = - resolution2D1[0] * image1centerX;
>   origin2D1[1] = - resolution2D1[1] * image1centerY;
>   origin2D1[2] = - scd;
>
>   imageReader2D1->GetOutput()->SetOrigin( origin2D1 );
>   rescaler2D1->GetOutput()->SetOrigin( origin2D1 );
>
>   //// 2D Image 2
>   //origin2D2[0] = - resolution2D2[0] * image2centerX;
>   //origin2D2[1] = - resolution2D2[1] * image2centerY;
>   //origin2D2[2] = - scd;
>
>   //imageReader2D2->GetOutput()->SetOrigin( origin2D2 );
>   //rescaler2D2->GetOutput()->SetOrigin( origin2D2 );
>
>   registration->SetFixedImageRegion(
> rescaler2D1->GetOutput()->GetBufferedRegion() );
>   //registration->SetFixedImageRegion2(
> rescaler2D2->GetOutput()->GetBufferedRegion() );
>
>   if (verbose)
>     {
>     std::cout << "2D image 1 size: "
>               << size2D1[0] << ", " << size2D1[1] << ", " << size2D1[2] <<
> std::endl
>               << "   resolution: "
>               << resolution2D1[0] << ", " << resolution2D1[1] << ", " <<
> resolution2D1[2] << std::endl
>               << "   and position: "
>               << origin2D1[0] << ", " << origin2D1[1] << ", " <<
> origin2D1[2] << std::endl;
>              // << "2D image 2 size: "
>               //<< size2D2[0] << ", " << size2D2[1] << ", " << size2D2[2]
> <<
> std::endl
>               //<< "   resolution: "
>              // << resolution2D2[0] << ", " << resolution2D2[1] << ", " <<
> resolution2D2[2] << std::endl
>             //  << "   and position: "
>            //   << origin2D2[0] << ", " << origin2D2[1] << ", " <<
> origin2D2[2] << std::endl;
>     }
>
>
>
> //----------------------------------------------------------------------------
> // Set the moving and fixed images' schedules
>
> //----------------------------------------------------------------------------
>   const unsigned int ResolutionLevels = 3;
>   std::cout<<"\nSet the moving and fixed images' schedules"<<std::endl;
>   RegistrationType::ScheduleType fixedSchedule( ResolutionLevels,Dimension
> );
>   fixedSchedule[0][0] = 4;
>   fixedSchedule[0][1] = 4;
>   fixedSchedule[0][2] = 1;
>   fixedSchedule[1][0] = 2;
>   fixedSchedule[1][1] = 2;
>   fixedSchedule[1][2] = 1;
>   fixedSchedule[2][0] = 1;
>   fixedSchedule[2][1] = 1;
>   fixedSchedule[2][2] = 1;
>
>   RegistrationType::ScheduleType movingSchedule(
> ResolutionLevels,Dimension);
>   movingSchedule[0][0] = 4;
>   movingSchedule[0][1] = 4;
>   movingSchedule[0][2] = 4;
>   movingSchedule[1][0] = 2;
>   movingSchedule[1][1] = 2;
>   movingSchedule[1][2] = 2;
>   movingSchedule[2][0] = 1;
>   movingSchedule[2][1] = 1;
>   movingSchedule[2][2] = 1;
>
>   registration->SetSchedules( fixedSchedule, movingSchedule );
>   std::cout<<"\nSet the moving and fixed images' schedules to
> registration"<<std::endl;
>
>
>
>
>
>         //                              //The old code start from here
>  // // Set the pyramid schedule for the 2D image
>  // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>  // // Software Guide : BeginLatex
>  // //
>  // // The multi-resolution hierachy is specified with respect to the 2D
>  // // image. The number of scales or levels defaults to two, with each
> level
>  // // differing in resolution by a factor of two. The highest resolution
> (final)
>  // // level corresponds to the resolution of the input 2D image but
>  // // this can be altered by the user via the \emph{maxScale} command
>  // // line parameter. The number of scales \emph{nScales} can also be
>  // // specified by the user.
>  // //
>  // // Software Guide : EndLatex
>  // std::vector<double> sampledRes2D;
>  // std::vector<ImageRegionType2D> imageRegionPyramid2D;
>
>  // imagePyramid2D->SetNumberOfLevels( nScales );
>  // std::cout<<"imagepyramid 2d"<<std::endl;
>  //
>  // sampledRes2D.reserve( nScales-1 );
>  // std::cout<<"sampledRes2D
> "<<sampledRes2D[0]<<std::endl;
>  // imageRegionPyramid2D.reserve( nScales );
>
>  // typedef ImagePyramidType2D::ScheduleType   ScheduleType2D;
>  // ScheduleType2D schedule2D = imagePyramid2D->GetSchedule();
>
>  // for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++ )
>  //   {
>  //   schedule2D[ nScales-1 ][ dim ] = maxScale;
>
> //std::cout<<"schedule2D["<<nScales-1<<"]["<<dim
> <<"]  "<<schedule2D[ nScales-1 ][ dim
> ]<<std::endl;
>  //   }
>
>  // for ( int level=nScales-2; level >= 0; level-- )
>  //   {
>  //   for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++
> )
>  //     {
>  //     schedule2D[ level ][ dim ] = 2*schedule2D[ level+1 ][ dim ];
>         //
> std::cout<<"schedule2D["<<level<<"]["<<dim
> <<"]  "<< 2*schedule2D[ level+1 ][ dim
> ]<<std::endl;
>  //     }
>  //   }
>  // std::cout<<"out of the imagepyramid type 2d
> "<<std::endl;
>  // // Compute the 2D ImageRegion corresponding to each level of the
>  // // pyramid.
>
>  // typedef ImageRegionType2D::IndexType       IndexType2D;
>  // IndexType2D       inputStart2D = region2D1.GetIndex();
>
>  // for ( unsigned int level=0; level < nScales; level++ )
>  //   {
>  //   SizeType2D  size;
>  //   IndexType2D start;
>  //   sampledRes2D[ level ] = 0.; //this is the error ORIGINAL
>         ////sampledRes2D[ level ] = 1.; //this is the error
>  //   for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension;
> dim++ )
>  //     {
>  //     float scaleFactor = static_cast<float>( schedule2D[ level ][ dim
> ] );
>
>  //     size[ dim ] = static_cast<SizeType2D::SizeValueType>(
>  //       floor( static_cast<float>( size2D1[ dim ] ) / scaleFactor ) );
>
>  //     if( size[ dim ] < 1 )
>  //       {
>  //       size[ dim ] = 1;
>  //       schedule2D[ level ][ dim ] = 1;
>  //       }
>
>  //      std::cout << level << " " << dim << " "
>  //               << size[ dim ] << " " << schedule2D[ level ][ dim ]
>  //               << std::endl;
>
>  //      scaleFactor = static_cast<float>( schedule2D[ level ][ dim ] );
>  //      start[ dim ] = static_cast<IndexType2D::IndexValueType>(
>  //        ceil(  static_cast<float>( inputStart2D[ dim ] ) / scaleFactor )
> );
>
>  //     if ( dim < 2 )
>  //       {
>  //       sampledRes2D[ level ] +=  scaleFactor*resolution2D1[ dim ]
>  //                                *scaleFactor*resolution2D1[ dim ];
>  //       }
>  //     }
>
>  //   sampledRes2D[ level ] = sqrt( sampledRes2D[ level ] );
>
>  //   imageRegionPyramid2D[ level ].SetSize( size );
>  //   imageRegionPyramid2D[ level ].SetIndex( start );
>  //   }
>
>  // imagePyramid2D->SetSchedule( schedule2D );
>
>
>   //// Set the pyramid schedule for the 3D image
>   //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>   //// Software Guide : BeginLatex
>   ////
>   //// The 3D image pyramid naturally contains the same number of levels
>   //// as the multi-resolution schedule for the 2D image. In addition
>   //// the sub-sampling factors for each level are set such that
>   //// resolution of the 3D image in each case is reduced to no greater
>   //// than the corresponding 2D image resolution at that scale. This
>   //// ensures that the 3D volume is reduced in size as far as possible,
>   //// minimising ray-casting computation time, whilst retaining
>   //// sufficient information to ensure accurate production of the DRR
>   //// at the resolution of the 2D image.
>   ////
>   //// Software Guide : EndLatex
>
>   //std::vector<ImageRegionType3D> imageRegionPyramid3D;
>
>   //imagePyramid3D->SetNumberOfLevels( nScales );
>   //
>   //imageRegionPyramid3D.reserve( nScales );
>
>   //typedef ImagePyramidType3D::ScheduleType   ScheduleType3D;
>   //ScheduleType3D schedule3D = imagePyramid3D->GetSchedule();
>   //
>   //// Compute the 2D image pyramid schedule such that the 3D volume
>   //// resolution is no greater than the 2D image resolution.
>
>   //for ( unsigned int level=0; level < nScales; level++ )
>   //  {
>   //  for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
>   //    {
>   //    double res = resolution3D[ dim ];
>   //    schedule3D[ level ][ dim ] = 1;
>   //    while ( res*2. < sampledRes2D[ level ] )
>   //      {
>   //      schedule3D[ level ][ dim ] *= 2;
>   //      res *= 2.;
>   //      }
>   //    }
>   //  }
>   //
>   //// Compute the 3D ImageRegion corresponding to each level of the
>   //// pyramid.
>
>   //typedef ImageRegionType3D::IndexType       IndexType3D;
>   //IndexType3D       inputStart3D = region3D.GetIndex();
>
>   //for ( unsigned int level=0; level < nScales; level++ )
>   //  {
>   //  SizeType3D  size;
>   //  IndexType3D start;
>   //  for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
>   //    {
>   //    float scaleFactor = static_cast<float>( schedule3D[ level ][ dim ]
> );
>
>   //    size[ dim ] = static_cast<SizeType3D::SizeValueType>(
>   //      floor( static_cast<float>( size3D[ dim ] ) / scaleFactor ) );
>
>   //    if( size[ dim ] < 1 )
>   //      {
>   //      size[ dim ] = 1;
>   //      schedule3D[ level ][ dim ] = 1;
>   //      }
>
>   //    scaleFactor = static_cast<float>( schedule3D[ level ][ dim ] );
>   //    start[ dim ] = static_cast<IndexType3D::IndexValueType>(
>   //      ceil(  static_cast<float>( inputStart3D[ dim ] ) / scaleFactor )
> );
>   //    }
>   //  imageRegionPyramid3D[ level ].SetSize( size );
>   //  imageRegionPyramid3D[ level ].SetIndex( start );
>   //  }
>
>   //imagePyramid3D->SetSchedule( schedule3D );
>
>
>   // Initialize the ray cast interpolator
>   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>   // The ray cast interpolator is used to project the 3D volume. It
>   // does this by casting rays from the (transformed) focal point to
>   // each (transformed) pixel coordinate in the 2D image.
>   //
>   // In addition a threshold may be specified to ensure that only
>   // intensities greater than a given value contribute to the
>   // projected volume. This can be used, for instance, to remove soft
>   // tissue from projections of CT data and force the registration
>   // to find a match which aligns bony structures in the images.
>
>   // 2D Image 1
>   interpolator1->SetProjectionAngle( dtr*projAngle1 );
>   interpolator1->SetFocalPointToIsocenterDistance(scd);
>   interpolator1->SetThreshold(threshold);
>   interpolator1->SetTransform(transform);
>
>   interpolator1->Initialize();
>   std::cout<<"\nInterpolator"<<std::endl;
>   // 2D Image 2
>  /* interpolator2->SetProjectionAngle( dtr*projAngle2 );
>   interpolator2->SetFocalPointToIsocenterDistance(scd);
>   interpolator2->SetThreshold(threshold);
>   interpolator2->SetTransform(transform);
>
>   interpolator2->Initialize();
> */
>
>   // Set up the transform and start position
>   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>   // The registration start position is intialised using the
>   // transformation parameters.
>
>
>
>   registration->SetInitialTransformParameters( transform->GetParameters()
> );
>                                         //trying to set the FIXED parameter
> //#if defined(ITK_FIXED_PARAMETERS_ARE_DOUBLE) // After 4.8.1
> //   TransformType::FixedParametersType fixedParameters(3);
> //#else                                         //Pre 4.8.1
> //  TransformType::ParametersType fixedParameters(6);
> //#endif
> //  fixedParameters[0] = 0.0;
> //  fixedParameters[1] = 0.0;
> //  fixedParameters[3] = 0.0;
> //  fixedParameters[4] = 0.0;
> //  fixedParameters[5] = 0.0;
> //  transform->SetFixedParameters( fixedParameters );
>
>   // We wish to minimize the negative normalized correlation similarity
> measure.
>
>   // optimizer->SetMaximize( true );  // for
> GradientDifferenceTwoImageToOneImageMetric
>   optimizer->SetMaximize( false );  // for NCC ORIGINAL
>   //optimizer->SetMaximize( true );             //NCC CHANGED
>
>   optimizer->SetMaximumStepLength( maxStepSize );
>   optimizer->SetMinimumStepLength( minStepSize );
>   //optimizer->SetStepLength(2);  // ORGINAL 4
>   //optimizer->SetStepTolerance( 0.01); //ORIGINAL 0.02
>   //optimizer->SetValueTolerance( 0.00001 );//ORIGINAL 0.001
>   //optimizer->SetMetricWorstPossibleValue(-0.2);
>   optimizer->SetNumberOfIterations(200);
>   // The optimizer weightings are set such that one degree equates to
>   // one millimeter.
>
>   itk::Optimizer::ScalesType weightings( transform->GetNumberOfParameters()
> );
>
>
> std::cout<<"transform->GetNumberOfParameters()==>"<<transform->GetNumberOfParameters()<<std::endl;
>  /* weightings[0] = 1./dtr;
>   weightings[1] = 1./dtr;               //ORIGINAL
>   weightings[2] = 1./dtr;
>   weightings[3] = 1.;
>   weightings[4] = 1.;
>   weightings[5] = 1.;*/
>
>   weightings[0] = 1./dtr;
>   weightings[1] = 1./dtr;
>   weightings[2] = 1./dtr;
>   weightings[3] = 1.0;
>   weightings[4] = 1.0;
>   weightings[5] = 1.0;
>
>   std::cout<<"WEIGHTINGS
> "<<weightings[0]<<std::endl;
>   optimizer->SetScales( weightings );
>  // optimizer->SetMetricWorstPossibleValue(0);//to change the stop
> condition
>   if (verbose)
>     {
>     optimizer->Print( std::cout );
>     }
>
>
>   // Create the observers
>   // ~~~~~~~~~~~~~~~~~~~~
>
>   CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>
>   optimizer->AddObserver( itk::IterationEvent(), observer );
>
>   // Software Guide : BeginLatex
>   //
>   //  Once all the registration components are in place we can create
>   //  an instance of the interface command and connect it to the
>   //  registration object using the \code{AddObserver()} method.
>   //
>   // Software Guide : EndLatex
>   // Software Guide : BeginCodeSnippet
>   typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
>   CommandType::Pointer command = CommandType::New();
>   registration->AddObserver( itk::IterationEvent(), command );
>   // Software Guide : EndCodeSnippet
>
>   // Software Guide : BeginLatex
>   //
>   // Finally we set the number of multi-resolution levels:
>   //
>   // Software Guide : EndLatex
>   // Software Guide : BeginCodeSnippet
>   //registration->SetNumberOfLevels( nScales );
>   //// Software Guide : EndCodeSnippet
>
>   //imagePyramid3D->Print(std::cout);
>   //imagePyramid2D->Print(std::cout);
>
>   // Start the registration
>   // ~~~~~~~~~~~~~~~~~~~~~~
>
>   // Create a timer to record calculation time.
>   itk::TimeProbesCollectorBase timer;
>
>   if (verbose)
>     {
>     std::cout << "Starting the registration now..." << std::endl;
>     }
>
>   try
>     {
>     timer.Start("Registration");
>     // Start the registration.
>     registration->StartRegistration();
>     timer.Stop("Registration");
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cout << "ExceptionObject caught !" << std::endl;
>     std::cout << err << std::endl;
>     return -1;
>     }
>
>   typedef RegistrationType::ParametersType ParametersType;
>   ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
>   const double RotationAlongX = finalParameters[0]/dtr; // Convert radian
> to
> degree
>   const double RotationAlongY = finalParameters[1]/dtr;
>   const double RotationAlongZ = finalParameters[2]/dtr;
>   const double TranslationAlongX = finalParameters[3];
>   const double TranslationAlongY = finalParameters[4];
>   const double TranslationAlongZ = finalParameters[5];
>
>   const int numberOfIterations = optimizer->GetCurrentIteration();
>
>   const double bestValue = optimizer->GetValue();
>
>   std::cout << "Result = " << std::endl;
>   std::cout << " Rotation Along X = " << RotationAlongX  << " deg" <<
> std::endl;
>   std::cout << " Rotation Along Y = " << RotationAlongY  << " deg" <<
> std::endl;
>   std::cout << " Rotation Along Z = " << RotationAlongZ  << " deg" <<
> std::endl;
>   std::cout << " Translation X = " << TranslationAlongX  << " mm" <<
> std::endl;
>   std::cout << " Translation Y = " << TranslationAlongY  << " mm" <<
> std::endl;
>   std::cout << " Translation Z = " << TranslationAlongZ  << " mm" <<
> std::endl;
>   std::cout << " Number Of Iterations = " << numberOfIterations <<
> std::endl;
>   std::cout << " Metric value  = " << bestValue          << std::endl;
>   std::cout<<"GET STOP
> "<<optimizer->GetStopConditionDescription()<<std::endl;
>
>   // Write out the projection images at the registration position
>   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
>   TransformType::Pointer finalTransform = TransformType::New();
>   // The transform is determined by the parameters and the rotation center.
>   finalTransform->SetParameters( finalParameters );
>   finalTransform->SetCenter(isocenter);
>
>   typedef itk::ResampleImageFilter< InternalImageType, InternalImageType >
> ResampleFilterType;
>
>   // The ResampleImageFilter is the driving force for the projection image
> generation.
>   ResampleFilterType::Pointer resampleFilter1 = ResampleFilterType::New();
>
>   resampleFilter1->SetInput( caster3D->GetOutput() ); // Link the 3D
> volume.
>   resampleFilter1->SetDefaultPixelValue( 0 );
>
>   // The parameters of interpolator1, such as ProjectionAngle and
> FocalPointToIsocenterDistance
>   // have been set before registration. Here we only need to replace the
> initial
>   // transform with the final transform.
>   interpolator1->SetTransform( finalTransform );
>   interpolator1->Initialize();
>   resampleFilter1->SetInterpolator( interpolator1 );
>
>   // The output 2D projection image has the same image size, origin, and
> the
> pixel spacing as
>   // those of the input 2D image.
>   resampleFilter1->SetSize(
> rescaler2D1->GetOutput()->GetLargestPossibleRegion().GetSize() );
>   resampleFilter1->SetOutputOrigin(  rescaler2D1->GetOutput()->GetOrigin()
> );
>   resampleFilter1->SetOutputSpacing( rescaler2D1->GetOutput()->GetSpacing()
> );
>
>   // Do the same thing for the output image 2.
>   ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New();
>   resampleFilter2->SetInput( caster3D->GetOutput() );
>   resampleFilter2->SetDefaultPixelValue( 0 );
>
>   // The parameters of interpolator2, such as ProjectionAngle and
> FocalPointToIsocenterDistance
>   // have been set before registration. Here we only need to replace the
> initial
>   // transform with the final transform.
>  /* interpolator2->SetTransform( finalTransform );
>   interpolator2->Initialize();
>   resampleFilter2->SetInterpolator( interpolator2 );
>
>   resampleFilter2->SetSize(
> rescaler2D2->GetOutput()->GetLargestPossibleRegion().GetSize() );
>   resampleFilter2->SetOutputOrigin(  rescaler2D2->GetOutput()->GetOrigin()
> );
>   resampleFilter2->SetOutputSpacing( rescaler2D2->GetOutput()->GetSpacing()
> );
> */
>
>
> /////////////////////////////---DEGUG--START----////////////////////////////////////
>   if (debug)
>     {
>     InternalImageType::PointType outputorigin2D1 =
> rescaler2D1->GetOutput()->GetOrigin();
>     std::cout << "Output image 1 origin: " << outputorigin2D1[0] << ", " <<
> outputorigin2D1[1]
>               << ", " << outputorigin2D1[2] << std::endl;
>     /*InternalImageType::PointType outputorigin2D2 =
> rescaler2D2->GetOutput()->GetOrigin();
>     std::cout << "Output image 2 origin: " << outputorigin2D2[0] << ", " <<
> outputorigin2D2[1]
>               << ", " << outputorigin2D2[2] << std::endl;*/
>     }
>
>
> /////////////////////////////---DEGUG--END----//////////////////////////////////////
>
>
>   // As explained before, the computed projection is upsided-down.
>   // Here we use a FilpImageFilter to flip the images in y-direction.
>   flipFilter1->SetInput( resampleFilter1->GetOutput() );
>   //flipFilter2->SetInput( resampleFilter2->GetOutput() ); //ORIGINAL
>   // Rescale the intensity of the projection images to 0-255 for output.
>   typedef itk::RescaleIntensityImageFilter<
>     InternalImageType, OutputImageType > RescaleFilterType;
>
>   RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
>   rescaler1->SetOutputMinimum(   0 );
>   rescaler1->SetOutputMaximum( 255 );
>   rescaler1->SetInput( flipFilter1->GetOutput() );      //ORIGINAL
>  // rescaler1->SetInput( resampleFilter1->GetOutput() );
>
>   //RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
>   //rescaler2->SetOutputMinimum(   0 );
>   //rescaler2->SetOutputMaximum( 255 );
>   //rescaler2->SetInput( flipFilter2->GetOutput() );    //ORIGINAL
>  // rescaler2->SetInput( resampleFilter2->GetOutput() );
>
>
>   typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>   WriterType::Pointer writer1 = WriterType::New();
>   //WriterType::Pointer writer2 = WriterType::New();
>
>   writer1->SetFileName( fileOutput1 );
>   writer1->SetInput( rescaler1->GetOutput() );
>
>   try
>     {
>     std::cout << "Writing image: " << fileOutput1 << std::endl;
>     writer1->Update();
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     }
>
>   //writer2->SetFileName( fileOutput2 );
>   //writer2->SetInput( rescaler2->GetOutput() );
>
>   //try
>   //  {
>   //  std::cout << "Writing image: " << fileOutput2 << std::endl;
>   //  writer2->Update();
>   //  }
>   //catch( itk::ExceptionObject & err )
>   //  {
>   //  std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
>   //  std::cerr << err << std::endl;
>   //  }
>   timer.Report();
>
>   return EXIT_SUCCESS;
> }
>
>
>
>
> --
> View this message in context:
> http://itk-users.7.n7.nabble.com/Multiresolution-Registration-error-while-trying-to-improve-a-journal-paper-tp36739p36751.html
> Sent from the ITK - Users mailing list archive at Nabble.com.
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20160208/c3a9bd41/attachment-0001.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list