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