[ITK-users] Multiresolution Registration error while trying to improve a journal paper
vishal
itkhelpacc at gmail.com
Mon Feb 8 03:03:23 EST 2016
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.
More information about the Insight-users
mailing list