[Insight-users] Image to Spatial object registration problems with Scene Spatial Object
fjlb at duke.edu
fjlb at duke.edu
Mon Jul 14 12:22:19 EDT 2008
Hi,
I am attempting to use the image to spatial object registration framework.
Basically I have a scene spatial object that contains two lines and a black and
white phantom with 2 tube-like lines in white on a black background. I want to
register these two so that the spatial objects go to the center of the tubes
(axis extraction). I smoothed the images so that the center of the tube is
roughly in the middle. I do this using a bspline transform with the lbfgs
optimizer and I wrote my own metric for this (metric basically tries to
minimize the external energy on the lines). However, the registration method is
not accepting the scene as an input (it does accept a regular line spatial
object though). Please let me know if my approach is in the right direction or
if there's a better way to do this. I attached the error I'm getting and the
code below if you're interested.
Thanks,
Federico
Here's the error I get.
****************************************************************************
Bspline2Obj2D.cxx
E:\Federico\CT\Source\Insight\Code\Algorithms\itkImageToSpatialObjectMetric.h(85)
: error C2039: 'ObjectDimension' : is not a member of
'itk::SceneSpatialObject<TSpaceDimension>'
with
[
TSpaceDimension=2
]
e:\Federico\CT\Source\Insight\Code\Algorithms\itkImageToSpatialObjectRegistrationMethod.h(111)
: see reference to class template instantiation
'itk::ImageToSpatialObjectMetric<TFixedImage,TMovingSpatialObject>' being
compiled
with
[
TFixedImage=itk::ImageToSpatialObjectRegistrationMethod<ImageType,SceneType>::FixedImageType,
TMovingSpatialObject=itk::ImageToSpatialObjectRegistrationMethod<ImageType,SceneType>::MovingSpatialObjectType
]
\Federico\CT\C\Bspline2Obj2D\Bspline2Obj2D.cxx(553) : see reference to
class template instantiation
'itk::ImageToSpatialObjectRegistrationMethod<TFixedImage,TMovingSpatialObject>'
being compiled
with
[
TFixedImage=ImageType,
TMovingSpatialObject=SceneType
]
E:\Federico\CT\Source\Insight\Code\Algorithms\itkImageToSpatialObjectMetric.h(85)
: error C2065: 'ObjectDimension' : undeclared identifier
E:\Federico\CT\Source\Insight\Code\Algorithms\itkImageToSpatialObjectMetric.h(91)
: error C2975: 'NInputDimensions' : invalid template argument for
'itk::Transform', compile-time evaluatable constant expression expected
e:\Federico\CT\Source\Insight\Code\Common\itkTransform.h(64) : see
declaration of 'NInputDimensions'
E:\Federico\CT\Source\Insight\Code\Algorithms\itkImageToSpatialObjectMetric.h(108)
: error C2975: 'n' : invalid template argument for 'vnl_vector_fixed',
compile-time evaluatable constant expression expected
E:\Federico\CT\Source\Insight\Utilities\vxl\core\vnl\vnl_vector_fixed.h(81) :
see declaration of 'n'
\Federico\CT\C\Bspline2Obj2D\Bspline2Obj2D.cxx(649) : error C2664:
'itk::ImageToSpatialObjectRegistrationMethod<TFixedImage,TMovingSpatialObject>::SetTransform'
: cannot convert parameter 1 from
'itk::BSplineDeformableTransform<TScalarType,NDimensions,VSplineOrder>::Pointer'
to
'itk::ImageToSpatialObjectRegistrationMethod<TFixedImage,TMovingSpatialObject>::TransformType
*'
with
[
TFixedImage=ImageType,
TMovingSpatialObject=SceneType
]
and
[
TScalarType=CoordinateRepType,
NDimensions=2,
VSplineOrder=3
]
and
[
TFixedImage=ImageType,
TMovingSpatialObject=SceneType
]
No user-defined-conversion operator available that can perform this
conversion, or the operator cannot be called
Build log was saved at
"file://e:\Federico\Ct\C\Bspline2Obj2D\bin\Bspline2Obj2D.dir\Debug\BuildLog.htm"
Bspline2Obj2D - 5 error(s), 0 warning(s)
****************************************************************************
Code:
****************************************************************************
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <new>
#include "itkSpatialObjectReader.h"
#include "itkSpatialObjectWriter.h"
#include "itkLineSpatialObject.h"
#include "itkLineSpatialObjectPoint.h"
#include "itkSceneSpatialObject.h"
#include "itkLineSpatialObjectPoint.h"
#include "itkImage.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkImageToSpatialObjectRegistrationMethod.h"
#include "itkImageToSpatialObjectMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkNormalVariateGenerator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkGradientDescentOptimizer.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkDerivativeImageFilter.h"
#include "itkCommand.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkLaplacianImageFilter.h"
#include "itkGradientVectorFlowImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkCovariantVector.h"
#include "itkExtractImageFilter.h"
#include "itkThresholdImageFilter.h"
#include "itkBSplineDeformableTransform.h"
#include "itkLBFGSOptimizer.h"
#include "itkCastImageFilter.h"
using namespace std;
template < class TOptimizer >
class IterationCallback : public itk::Command
{
public:
typedef IterationCallback Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
itkTypeMacro( IterationCallback, Superclass );
itkNewMacro( Self );
/** Type defining the optimizer. */
typedef TOptimizer OptimizerType;
/** Method to specify the optimizer. */
void SetOptimizer( OptimizerType * optimizer )
{
m_Optimizer = optimizer;
m_Optimizer->AddObserver( itk::IterationEvent(), this );
}
/** Execute method will print data at each iteration */
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object *, const itk::EventObject & event)
{
if( typeid( event ) == typeid( itk::StartEvent ) )
{
//std::cout << std::endl << "Position Value";
//std::cout << std::endl << std::endl;
}
else if( typeid( event ) == typeid( itk::IterationEvent ) )
{
//std::cout << m_Optimizer->GetCurrentIteration() << " ";
//std::cout << m_Optimizer->GetValue() << " ";
// std::cout << m_Optimizer->GetCurrentPosition() << std::endl;
}
else if( typeid( event ) == typeid( itk::EndEvent ) )
{
//std::cout << std::endl << std::endl;
//std::cout << "After " << m_Optimizer->GetCurrentIteration();
//std::cout << " iterations " << std::endl;
//std::cout << "Solution is = " << m_Optimizer->GetCurrentPosition();
//std::cout << std::endl;
}
}
// Software Guide : EndCodeSnippet
protected:
IterationCallback() {};
itk::WeakPointer<OptimizerType> m_Optimizer;
};
template <typename TFixedImage, typename TMovingSpatialObject>
class SnakeImageToSpatialObjectMetric :
public itk::ImageToSpatialObjectMetric<TFixedImage, TMovingSpatialObject>
{
// Software Guide : EndCodeSnippet
public:
/** Standard class typedefs. */
typedef SnakeImageToSpatialObjectMetric Self;
typedef itk::ImageToSpatialObjectMetric< TFixedImage, TMovingSpatialObject
> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef itk::Point< double, 2 > PointType;
typedef std::list<PointType> PointListType;
typedef TMovingSpatialObject MovingSpatialObjectType;
typedef typename Superclass::ParametersType ParametersType;
typedef typename Superclass::DerivativeType DerivativeType;
typedef typename Superclass::MeasureType MeasureType;
typedef float PixelType;
//typedef unsigned char PixelType;
typedef itk::Image< PixelType, 2 > ImageType2D;
typedef itk::LineSpatialObject< 2 > LineType;
typedef itk::LineSpatialObjectPoint< 2 > LinePointType;
typedef itk::SceneSpatialObject< 2 > SceneSpatialObjectType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(SnakeImageToSpatialObjectMetric, ImageToSpatialObjectMetric);
itkStaticConstMacro( ParametricSpaceDimension, unsigned int, 2 );
/** Specify the fixed image and its gradient */
void SetFixedImage( const FixedImageType * image )
{
this->m_FixedImage = image;
GetImageDerivative();
alpha = 0.000;
beta = 0.0;
gamma = 1;
}
/** Specify the moving spatial object. */
void SetMovingSpatialObject( const MovingSpatialObjectType * object)
{
if(!this->m_FixedImage)
{
std::cout << "Please set the image before the moving spatial object"
<< std::endl;
return;
}
this->m_MovingSpatialObject = object;
MovingSpatialObjectType* scene = const_cast < MovingSpatialObjectType * >
(object);
SceneType::ObjectType * sceneobject;
sceneobject = scene->GetObjectById(1);
LineType * Line = static_cast < LineType * > (sceneobject);
sceneobject = scene->GetObjectById(2);
LineType * Line1 = static_cast <LineType * > (sceneobject);
SceneType::ObjectListType * myObjectList = scene->GetObjects();
std::cout << "Number of objects in SceneSpatialObject = ";
std::cout << myObjectList->size() << std::endl;
int N, i=0, N1;
LineType::PointListType pointlist = Line->GetPoints();
LineType::PointListType pointlist1 = Line1->GetPoints();
std::cout << "Number of points representing the line: ";
std::cout << "Number of points representing the line: ";
N = pointlist.size();
N1 = pointlist1.size();
std::cout << N << " " << N1 << std::endl;
LineType::PointListType::const_iterator it = Line->GetPoints().begin();
m_PointList.clear();
while ( it != Line->GetPoints().end() )
{
LinePointType::PointType p = (*it).GetPosition();
m_PointList.push_back(p);
std::cout << "Initial point " << p << std::endl;
it++;
}
LineType::PointListType::const_iterator it = Line1->GetPoints().begin();
while ( it != Line1->GetPoints().end() )
{
LinePointType::PointType p = (*it).GetPosition();
m_PointList.push_back(p);
std::cout << "Initial point " << p << std::endl;
it++;
}
std::cout << "Point list size (2 objects) N = " << static_cast<unsigned int>(
m_PointList.size() ) << std::endl;
}
//unsigned int GetNumberOfParameters(void) const {return m_PointList.size() *
2;};
unsigned int GetNumberOfParameters(void) const {return
m_Transform->GetNumberOfParameters();};
/** Get the Derivatives of the Match Measure */
void GetDerivative( const ParametersType &, DerivativeType & Derivative) const
{
//unsigned int parameterDimension = GetNumberOfParameters();
unsigned int parameterDimension = m_PointList.size() * 2;
Derivative = DerivativeType( parameterDimension );
//Get PointList into x,y,z, arrays
int i = 0;
//Get number of parameters
unsigned int N = static_cast<unsigned int>( m_PointList.size() );
std::cout << "Point list size N = " << N << std::endl;
double *x, *y, *dExt;
int M = N*2;
x = new (nothrow) double[N];
y = new (nothrow) double[N];
dExt = new (nothrow) double[M];
typedef TransformType::ParametersType ParametersType;
const unsigned int numberOfParameters = m_Transform->GetNumberOfParameters();
std::cout << "Number of Parameters = " << numberOfParameters << std::endl;
ParametersType parameters( numberOfParameters );
//Get transformed points into x,y
PointListType::const_iterator it = m_PointList.begin();
int j = 0;
while( it != m_PointList.end() )
{
PointType transformedPoint = this->m_Transform->TransformPoint(*it);
x[j] = transformedPoint[0];
y[j] = transformedPoint[1];
//parameters = m_Transform->GetParameters();
//std::cout <<"Parameters = " << parameters << std::endl;
std::cout <<"Transformed Points = " << transformedPoint << " " << x[j] << " "
<< y[j] << std::endl;
j++;
it++;
}
//Get Image size to check if transformed point is inside image
ImageType2D::RegionType region = m_FixedImage->GetLargestPossibleRegion();
ImageType2D::SizeType size = region.GetSize();
//std::cout << "Size = " << size << std::endl;
//std::cout << "Calculating Metric Derivative"<< std::endl;
float eps = -1;
//Calculate dExt
for( i = 0; i < N; i++)
{
ImageType2D::IndexType pixelIndex;
pixelIndex[0] = x[i];
pixelIndex[1] = y[i];
//std::cout <<"Pixel Index = " << pixelIndex << " " << pixelIndex[0] << " "
<< pixelIndex[1] << " " << x[i] << " " << y[i] << std::endl;
//Check if pixel is inside image boundaries
//if ( (pixelIndex[0] >= 0) && (pixelIndex[0] <= (size[0]-1)) &&
(pixelIndex[1] >= 0) && (pixelIndex[1] <= (size[1]-1)) )
if ( (x[i] >= 0) && (x[i] <= (size[0]-1)) && (y[i] >= 0) && (y[i] <=
(size[1]-1)) )
{
ImageType2D::PixelType pixelValue1 = m_xderiv->GetPixel(
pixelIndex );
dExt[ i*2 ] = (eps)*pixelValue1;
Derivative[ i*2 ] = (eps)*pixelValue1;
ImageType2D::PixelType pixelValue2 = m_yderiv->GetPixel( pixelIndex );
dExt[ i*2+1 ] = (eps)*pixelValue2;
Derivative[ i*2+1 ] = (eps)*pixelValue2;
//std::cout<< "dExt = " << dExt[i*2] << " " << dExt[i*2+1] << std::endl;
}
else
{
dExt[ i*2 ] = 0;
Derivative[ i*2 ] = 0;
dExt[ i*2+1 ] = 0;
Derivative[ i*2+1 ] = 0;
//std::cout << "Transformed Point is out of bounds " << std::endl;
}
}
return;
}
MeasureType GetValue( const ParametersType & parameters ) const
{
double value = 0, Eext = 0;
double y1der = 0, y2der = 0, all1der = 0, all2der = 0;
double x1der = 0, x2der = 0, temp = 0;
int i = 0, N;
//Dynamically allocate space for points
N = static_cast<unsigned int>( m_PointList.size() );
double *x, *y;
x = new (nothrow) double[N];
y = new (nothrow) double[N];
this->m_Transform->SetParameters( parameters );
PointListType::const_iterator it = m_PointList.begin();
//Get the transformed points
while( it != m_PointList.end() )
{
PointType transformedPoint = this->m_Transform->TransformPoint(*it);
x[i] = transformedPoint[0];
y[i] = transformedPoint[1];
//std::cout << "Transformed Point: " << x[i] << " " << y[i] << std::endl;
i++;
it++;
}
//std::cout << "Average Distance = " << dist/N << std::endl;
ImageType2D::IndexType pixelIndex;
ImageType2D::PixelType pixelValue;
ImageType2D::RegionType region = m_FixedImage->GetLargestPossibleRegion();
ImageType2D::SizeType size = region.GetSize();
//Get Eext for snake
for(unsigned int j=0; j < N; j++)
{
pixelIndex[0] = x[j];
pixelIndex[1] = y[j];
//if ( (pixelIndex[0] >= 0) && (pixelIndex[0] <= (size[0]-1)) &&
(pixelIndex[1] >= 0) && (pixelIndex[1] <= (size[1]-1)) )
if ( (x[j] >= 0) && (x[j] <= (size[0]-1)) && (y[j] >= 0) && (y[j] <=
(size[1]-1)) )
{
pixelValue = m_FixedImage->GetPixel( pixelIndex );
Eext = Eext - pixelValue;
}
else
{
Eext = Eext - 0;
//std::cout << "Transformed point is out of bounds " << std::endl;
}
}
//std::cout << "Eext = " << Eext << std::endl;
//std::cout << "Points read succesfully inside GetValue " << std::endl;
value = Eext;
//std::cout << Eint << " " << Eext << " " << value << std::endl;
std::cout << "GetValue exit successful - Value = " << value << std::endl;
return value;
}
// Software Guide : EndCodeSnippet
/** Get Value and Derivatives for MultipleValuedOptimizers */
void GetValueAndDerivative( const ParametersType & parameters,
MeasureType & Value, DerivativeType & Derivative ) const
{
Value = this->GetValue(parameters);
this->GetDerivative(parameters,Derivative);
}
private:
PointListType m_PointList;
ImageType2D::Pointer m_xderiv;
ImageType2D::Pointer m_yderiv;
double alpha, beta, gamma;
void GetImageDerivative( )
{
std::cout << "Calculating Image Derivative"<< std::endl;
typedef double ComponentType;
const unsigned int VDimension = 2;
typedef itk::CovariantVector< ComponentType, VDimension
> VectorPixelType;
typedef itk::Image< VectorPixelType, VDimension > VectorImageType;
typedef itk::Image< double, VDimension > InternalImageType;
typedef itk::LaplacianImageFilter< InternalImageType, InternalImageType
> LaplacianFilterType;
typedef itk::GradientRecursiveGaussianImageFilter< ImageType2D,
VectorImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetSigma( 1.0 );
filter->SetInput( m_FixedImage );
filter->Update();
typedef itk::GradientVectorFlowImageFilter< VectorImageType, VectorImageType
> FilterType2;
FilterType2::Pointer filter2 = FilterType2::New();
LaplacianFilterType::Pointer lapfilter = LaplacianFilterType::New();
filter2->SetInput( filter->GetOutput() );
filter2->SetLaplacianFilter( lapfilter );
filter2->SetIterationNum( 5 );
//filter2->SetNoiseLevel( 500 );
filter2->Update();
typedef itk::VectorIndexSelectionCastImageFilter< VectorImageType,
ImageType2D > FilterType3;
FilterType3::Pointer componentExtractor1 = FilterType3::New();
unsigned int indexOfComponentToExtract = 0;
componentExtractor1->SetIndex( indexOfComponentToExtract );
componentExtractor1->SetInput( filter2->GetOutput() );
componentExtractor1->Update();
m_xderiv = componentExtractor1->GetOutput();
FilterType3::Pointer componentExtractor2 = FilterType3::New();
indexOfComponentToExtract = 1;
componentExtractor2->SetIndex( indexOfComponentToExtract );
componentExtractor2->SetInput( filter2->GetOutput() );
componentExtractor2->Update();
m_yderiv = componentExtractor2->GetOutput();
typedef itk::RescaleIntensityImageFilter< ImageType2D, ImageType2D
> FilterType4;
FilterType4::Pointer filter4 = FilterType4::New();
filter4->SetInput( m_xderiv );
filter4->SetOutputMinimum( 0 );
filter4->SetOutputMaximum( 255 );
filter4->Update();
typedef unsigned char OPixelType;
typedef itk::Image < OPixelType, ImageDimension > OImageType;
typedef itk::CastImageFilter< ImageType2D, OImageType > CastFilterType;
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput( filter4->GetOutput() );
caster->Update();
//typedef itk::Image< unsigned char, 2 > ImageType2;
/*typedef itk::ImageFileWriter < OImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( "2dphantom_xder5.png" );
writer->SetInput( caster->GetOutput() );
writer->Update();*/
/*filter4->SetInput( m_yderiv );
filter4->Update();
caster->SetInput( filter4->GetOutput() );
caster->Update();
writer->SetFileName( "2dphantom_yder5.png" );
writer->SetInput( caster->GetOutput() );
writer->Update();*/
}
};
int main( int argc, char *argv[] )
{
if( argc < 3 )
{
std::cerr << "Too many parameters " << std::endl;
std::cerr << "Usage: VesselnessImage RCA_SpatialObject Results(txt file)" <<
std::endl;
return EXIT_FAILURE;
}
//Define Image Parameters
//typedef unsigned short PixelType;
typedef float PixelType;
const unsigned int ImageDimension = 2;
typedef itk::Image< PixelType, ImageDimension > ImageType;
//Define Image Reader
typedef itk::ImageFileReader< ImageType > ImageReaderType;
ImageReaderType::Pointer imagereader = ImageReaderType::New();
//Define Rescale Intensity Filter
typedef itk::RescaleIntensityImageFilter< ImageType, ImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
//Define Threshold Filter
typedef itk::ThresholdImageFilter< ImageType > FilterType2;
FilterType2::Pointer filter2 = FilterType2::New();
//Read in Vesselness (Cost) Image
imagereader->SetFileName( argv[1] );
//imagereader->Update();
filter->SetInput( imagereader->GetOutput() );
filter->SetOutputMinimum( 0 );
filter->SetOutputMaximum( 1 );
filter->Update();
filter2->SetInput( filter->GetOutput() );
filter2->SetOutsideValue( 0 );
filter2->ThresholdBelow( 0.01 );
filter2->Update();
ImageType::Pointer image = filter2->GetOutput();
//typedef itk::Image< unsigned char, ImageDimension > ImageType2;
typedef unsigned char OPixelType;
typedef itk::Image < OPixelType, ImageDimension > OImageType;
typedef itk::CastImageFilter< ImageType, OImageType > CastFilterType;
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput( imagereader->GetOutput() );
caster->Update();
typedef itk::ImageFileWriter < OImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
/*writer->SetFileName( "2dphantom_copy.png" );
writer->SetInput( caster->GetOutput() );
writer->Update();*/
//ImageType::Pointer image = imagereader->GetOutput();
//ImageType::RegionType region = image->GetLargestPossibleRegion();
//ImageType::SizeType size = region.GetSize();
//std::cout <<"Inside Main " << size << std::endl;
std::cout << "Vesselness Image read in successfully " << std::endl;
//Define Spatial Object Parameters
const unsigned int SODimension = 2;
typedef itk::SceneSpatialObject< SODimension > SceneType;
SceneType::ObjectType * sceneobject;
typedef itk::LineSpatialObject< SODimension > LineType;
typedef LineType::Pointer LinePointer;
typedef itk::LineSpatialObjectPoint< SODimension > LinePointType;
//Define Spatial Object Reader
typedef itk::SpatialObjectReader< SODimension > SOReaderType;
SOReaderType::Pointer soreader = SOReaderType::New();
//Read Spatial Object in
soreader->SetFileName( argv[2] );
soreader->Update();
SOReaderType::SceneType * scene = soreader->GetScene();
sceneobject = scene->GetObjectById(1);
LineType * Line = static_cast < LineType * > (sceneobject);
sceneobject = scene->GetObjectById(2);
LineType * Line1 = static_cast < LineType * > (sceneobject);
std::cout << "Spatial Object read in successfully" << std::endl;
//Set Spatial Object parameters into array
double *x, *y, *x1, *y1;
int N, i = 0, N1;
LineType::PointListType pointlist = Line->GetPoints();
LineType::PointListType pointlist1 = Line1->GetPoints();
N = pointlist.size();
N1 = pointlist1.size();
x = new (nothrow) double[N];
y = new (nothrow) double[N];
x1 = new (nothrow) double[N1];
y1 = new (nothrow) double[N1];
std::cout << "Number of points representing the lines: ";
std::cout << N << " " << N1 <<std::endl;
LineType::PointListType::const_iterator it = Line->GetPoints().begin();
while ( it != Line->GetPoints().end() )
{
LinePointType::PointType p = (*it).GetPosition();
x[i] = p[0];
y[i] = p[1];
it++;
i++;
}
it = Line1->GetPoints().begin();
i = 0;
while ( it != Line1->GetPoints().end() )
{
LinePointType::PointType p = (*it).GetPosition();
x1[i] = p[0];
y1[i] = p[1];
it++;
i++;
}
std::cout << "Point List succesfully stored in arrays x, y, x1, y1" <<
std::endl;
std::cout << "x = " << x << " y = " << y << std::endl;
std::cout << "x1 = " << x1 << " y1 = " << y1 << std::endl;
//Define Registration Method
typedef itk::ImageToSpatialObjectRegistrationMethod< ImageType, SceneType
> RegistrationType;
RegistrationType::Pointer registration = RegistrationType::New();
//Define Interpolation Method
typedef itk::LinearInterpolateImageFunction< ImageType, double
> InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
//Define Transform
const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineDeformableTransform <CoordinateRepType, SpaceDimension,
SplineOrder > TransformType;
TransformType::Pointer transform = TransformType::New();
ImageType::RegionType fixedRegion = image->GetBufferedRegion();
typedef TransformType::RegionType RegionType;
RegionType bsplineRegion;
RegionType::SizeType gridSizeOnImage;
RegionType::SizeType gridBorderSize;
RegionType::SizeType totalGridSize;
typedef TransformType::SpacingType SpacingType;
SpacingType spacing = image->GetSpacing();
typedef TransformType::OriginType OriginType;
OriginType origin = image->GetOrigin();
std::cout << "Fixed image spacing " << spacing << std::endl;
std::cout << "Fixed image origin " << origin << std::endl;
std::cout << "Buffered region size " << fixedRegion.GetSize() << std::endl;
gridSizeOnImage.Fill( 5 );
gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper
)
totalGridSize = gridSizeOnImage + gridBorderSize;
std::cout << "Grid size on Image " << gridSizeOnImage << std::endl;
bsplineRegion.SetSize( totalGridSize );
ImageType::SizeType fixedImageSize = fixedRegion.GetSize();
//Set grid uniform spacing
for(unsigned int r=0; r<ImageDimension; r++)
{
spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1) /
static_cast<double>(gridSizeOnImage[r] - 1) );
origin[r] -= spacing[r];
}
std::cout << "Grid spacing " << spacing << std::endl;
std::cout << "Grid origin " << origin << std::endl;
transform->SetGridSpacing( spacing );
transform->SetGridOrigin( origin );
transform->SetGridRegion( bsplineRegion );
std::cout << "Bsplineregion: " << bsplineRegion << std::endl;
//Define Optimizer Type
typedef itk::LBFGSOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
//Set optimizer parameters together
optimizer->SetGradientConvergenceTolerance( 0.05 );
optimizer->SetLineSearchAccuracy( 0.9 );
optimizer->SetDefaultStepLength( 0.01 ) ;
optimizer->TraceOn();
optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
//Define Optimizer Type
/*typedef itk::GradientDescentOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
//Set optimizer parameters together
double LearningRate = 1;
unsigned long numberOfIterations = 50;
optimizer->SetNumberOfIterations( numberOfIterations );
optimizer->SetLearningRate( LearningRate );*/
//Define Metric
typedef SnakeImageToSpatialObjectMetric< ImageType, SceneType > MetricType;
MetricType::Pointer metric = MetricType::New();
//Define Iteration Callback
typedef IterationCallback< OptimizerType > IterationCallbackType;
IterationCallbackType::Pointer callback = IterationCallbackType::New();
callback->SetOptimizer( optimizer );
//Put registration together
registration->SetFixedImage( image );
//registration->SetMovingSpatialObject( Line );
registration->SetMovingSpatialObject( scene );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetOptimizer( optimizer );
registration->SetMetric( metric );
typedef TransformType::ParametersType ParametersType;
const unsigned int numberOfParameters = transform->GetNumberOfParameters();
std::cout << "# of Parameters = " << numberOfParameters << std::endl;
ParametersType parameters( numberOfParameters );
parameters.Fill( 0.0 );
transform->SetParameters( parameters );
registration->SetInitialTransformParameters( transform->GetParameters() );
// Software Guide : EndCodeSnippet
std::cout << "Intial Parameters = " << std::endl;
std::cout << transform->GetParameters() << std::endl;
optimizer->MinimizeOn();
/*TransformType::ParametersType initialParameters;
TransformType::ParametersType parametersScale;
initialParameters.set_size( 2*N );
parametersScale.set_size( 2*N );
for ( i = 0; i < N; i++ )
{
initialParameters[ i*2 ] = x[i];
initialParameters[ i*2+1 ] = y[i];
parametersScale[ i*2 ] = 1;
parametersScale[ i*2+1 ] = 1;
}
transform->SetParameters( initialParameters );
optimizer->SetScales( parametersScale );
optimizer->MinimizeOn();
registration->SetInitialTransformParameters( initialParameters );
std::cout << "Initial Parameters: " << transform->GetNumberOfParameters() <<
std::endl;
std::cout << "Registration Set" << std::endl;*/
//Set transform initial parameters
//Trigger registration execution
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & exp )
{
std::cerr << "Exception caught ! " << std::endl;
std::cerr << exp << std::endl;
}
//Get final Parameters
RegistrationType::ParametersType finalParameters =
registration->GetLastTransformParameters();
//std::cout << "\n Final solution is: " << finalParameters << std::endl;
//Write results to file
ofstream myfile;
myfile.open ( argv[3] );
for ( unsigned int j = 0; j<N; j++)
{
myfile << finalParameters[ 2*j ] << "\t" << finalParameters[ 2*j + 1] <<
std::endl;
}
myfile.close();
std::cout << "Results succesfully written to: " << argv[3] << std::endl;
return 0;
}
More information about the Insight-users
mailing list