[Insight-users] LineSpatialObject: problem with applying a
transformation
Hoveijn, I
i.hoveijn at rt.umcg.nl
Mon Mar 19 11:10:09 EST 2007
Dear ITK-users,
I am using BSplineDeformableTransform to find a transformation form one
image to another. Then I would like to use the transformation to map a
line in the first image to a line in the second image. By analogy I cut
and pasted the code below. Everything works fine except transforming the
line. Please let me know if more details are needed. I would be grateful
if someone could help me.
Remark. The header of the procedure may look strange. This code is
compiled and linked into a dll file and then called form Matlab.
Igor Hoveijn
########################################################################
############################################
/////////////////////////////////////////////////////////////
// Includes /////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// gebaseerd op DeformableRegistration7.cxx
// Include own header
#include "deform2.h"
// ITK includes
#include "itkBSplineDeformableTransform.h"
#include "itkCastImageFilter.h"
#include "itkCommand.h"
#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkImageRegistrationMethod.h"
#include "itkImportImageFilter.h"
#include "itkLBFGSBOptimizer.h"
#include "itkLineSpatialObject.h"
#include "itkLineSpatialObjectPoint.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkResampleImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSpatialObjectReader.h"
#include "itkSpatialObjectWriter.h"
#include "itkSubtractImageFilter.h"
// General includes
#include<fstream>
using namespace std;
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::LBFGSBOptimizer 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( !(itk::IterationEvent().CheckEvent( &event )) )
{
return;
}
ofstream logfile("deform2.log", ios::app);
logfile << optimizer->GetCurrentIteration() << " ";
logfile << optimizer->GetValue() << " ";
logfile << optimizer->GetInfinityNormOfProjectedGradient() <<
endl;
logfile.close();
}
};
// nonlinear transformation in 2 dimensions
__declspec(dllexport) void deformtrans(float *fim, unsigned short *fdim,
double *fspac, double *forg, float *mim, unsigned short *mdim, double*
mspac, double *morg, double cfcf, double pgt, unsigned int mni, unsigned
int mne, unsigned int mnc)
{
ofstream logfile("deform2.log", ios::app);
const unsigned int ImageDimension = 2;
typedef float PixelType;
typedef itk::Image< PixelType, ImageDimension > FixedImageType;
typedef itk::Image< PixelType, ImageDimension > MovingImageType;
typedef itk::LineSpatialObject<ImageDimension> LineType;
typedef LineType::Pointer LinePointer;
typedef itk::LineSpatialObjectPoint<ImageDimension> LinePointType;
LinePointer Line = LineType::New();
LineType::PointListType list;
for(unsigned int i=0; i<3; i++)
{
LinePointType p;
p.SetPosition(100+10*i,100+10*i);
p.SetColor(1,0,0,1);
// VectorType normal1;
// VectorType normal2;
// for(unsigned int j=0;j<3;j++)
// {
// normal1[j]=j;
// normal2[j]=j*2;
// }
//
// p.SetNormal(normal1,0);
// p.SetNormal(normal2,1);
list.push_back(p);
}
Line->GetProperty()->SetName("Line1");
Line->SetId(1);
Line->SetPoints(list);
LineType::PointListType pointList = Line->GetPoints();
logfile << "Number of points representing the line: ";
logfile << pointList.size() << endl;
LineType::PointListType::const_iterator it =
Line->GetPoints().begin();
while(it != Line->GetPoints().end())
{
logfile << "Position = " << (*it).GetPosition() << endl;
logfile << "Color = " << (*it).GetColor() << endl;
// logfile << "First normal = " << (*it).GetNormal(0) << endl;
// logfile << "Second normal = " << (*it).GetNormal(1) << endl;
logfile << endl;
it++;
}
typedef itk::SpatialObjectWriter<ImageDimension> LineWriterType;
LineWriterType::Pointer linewriter = LineWriterType::New();
linewriter->SetInput(Line);
linewriter->SetFileName("lijn.meta");
linewriter->Update();
const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineDeformableTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > TransformType;
typedef itk::LBFGSBOptimizer OptimizerType;
typedef itk::MeanSquaresImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double >
InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType >
RegistrationType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
// Obtain moving image from function arguments (Matlab input)
typedef itk::ImportImageFilter<PixelType,ImageDimension>
ImportFilterType;
ImportFilterType::Pointer importFilterM =
ImportFilterType::New();
ImportFilterType::Pointer importFilterF =
ImportFilterType::New();
ImportFilterType::SizeType sizeM;
ImportFilterType::SizeType sizeF;
double originM[ImageDimension], originF[ImageDimension],
spacingM[ImageDimension], spacingF[ImageDimension];
for (int i=0;i<ImageDimension;i++)
{
sizeM[i]=mdim[i];
sizeF[i]=fdim[i];
originM[i]=morg[i];
originF[i]=forg[i];
spacingM[i]=mspac[i];
spacingF[i]=fspac[i];
}
const unsigned int numberOfPixelsM=mdim[0]*mdim[1]*mdim[2];
const unsigned int numberOfPixelsF=fdim[0]*fdim[1]*fdim[2];
ImportFilterType::IndexType start;
start.Fill(0);
ImportFilterType::RegionType regionM;
ImportFilterType::RegionType regionF;
regionM.SetIndex(start);
regionM.SetSize(sizeM);
regionF.SetIndex(start);
regionF.SetSize(sizeF);
importFilterM->SetRegion(regionM);
importFilterF->SetRegion(regionF);
importFilterM->SetOrigin(originM);
importFilterF->SetOrigin(originF);
importFilterM->SetSpacing(spacingM);
importFilterF->SetSpacing(spacingF);
importFilterM->SetImportPointer(mim,numberOfPixelsM,false);
importFilterF->SetImportPointer(fim,numberOfPixelsF,false);
FixedImageType::ConstPointer fixedImage = importFilterF->GetOutput();
registration->SetFixedImage( fixedImage );
registration->SetMovingImage( importFilterM->GetOutput() );
importFilterF->Update();
FixedImageType::RegionType fixedRegion =
fixedImage->GetBufferedRegion();
registration->SetFixedImageRegion( fixedRegion );
typedef TransformType::RegionType RegionType;
RegionType bsplineRegion;
RegionType::SizeType gridSizeOnImage;
RegionType::SizeType gridBorderSize;
RegionType::SizeType totalGridSize;
gridSizeOnImage.Fill( 5 );
gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower,
2 upper )
totalGridSize = gridSizeOnImage + gridBorderSize;
bsplineRegion.SetSize( totalGridSize );
typedef TransformType::SpacingType SpacingType;
SpacingType spacing = fixedImage->GetSpacing();
typedef TransformType::OriginType OriginType;
OriginType origin = fixedImage->GetOrigin();;
FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
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];
}
transform->SetGridSpacing( spacing );
transform->SetGridOrigin( origin );
transform->SetGridRegion( bsplineRegion );
typedef TransformType::ParametersType ParametersType;
const unsigned int numberOfParameters =
transform->GetNumberOfParameters();
ParametersType parameters( numberOfParameters );
parameters.Fill( 0.0 );
transform->SetParameters( parameters );
registration->SetInitialTransformParameters(
transform->GetParameters() );
OptimizerType::BoundSelectionType boundSelect(
transform->GetNumberOfParameters() );
OptimizerType::BoundValueType upperBound(
transform->GetNumberOfParameters() );
OptimizerType::BoundValueType lowerBound(
transform->GetNumberOfParameters() );
boundSelect.Fill( 0 );
upperBound.Fill( 0.0 );
lowerBound.Fill( 0.0 );
optimizer->SetBoundSelection( boundSelect );
optimizer->SetUpperBound( upperBound );
optimizer->SetLowerBound( lowerBound );
optimizer->SetCostFunctionConvergenceFactor( cfcf );
optimizer->SetProjectedGradientTolerance( pgt );
optimizer->SetMaximumNumberOfIterations( mni );
optimizer->SetMaximumNumberOfEvaluations( mne );
optimizer->SetMaximumNumberOfCorrections( mnc );
CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
logfile << endl << "start registratie" << endl;
registration->StartRegistration();
logfile << endl << "einde registratie" << endl;
}
catch( itk::ExceptionObject & err )
{
logfile << "ExceptionObject caught 1" << endl;
logfile << err << endl;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
logfile << "Last Transform Parameters" << endl;
logfile << finalParameters << endl;
transform->SetParameters( finalParameters );
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( transform );
resample->SetInput( importFilterM->GetOutput() );
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize()
);
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, ImageDimension >
OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
char *fnaam = "c.dcm";
writer->SetFileName( fnaam );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
try
{
writer->Update();
}
catch( itk::ExceptionObject & err )
{
logfile << "ExceptionObject caught 2" << endl;
logfile << err << endl;
}
linewriter->SetTransform( transform );
linewriter->SetInput(Line);
linewriter->SetFileName("lijn2.meta");
try
{
linewriter->Update();
}
catch( itk::ExceptionObject & err )
{
logfile << "ExceptionObject caught 2b" << endl;
logfile << err << endl;
}
typedef itk::SubtractImageFilter<
FixedImageType,
FixedImageType,
FixedImageType > DifferenceFilterType;
DifferenceFilterType::Pointer difference =
DifferenceFilterType::New();
difference->SetInput1( importFilterF->GetOutput() );
difference->SetInput2( resample->GetOutput() );
typedef itk::RescaleIntensityImageFilter<
FixedImageType,
OutputImageType > RescalerType;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput( difference->GetOutput() );
intensityRescaler->SetOutputMinimum( 0 );
intensityRescaler->SetOutputMaximum( 255 );
resample->SetDefaultPixelValue( 1 );
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput( intensityRescaler->GetOutput() );
fnaam = "d.dcm";
writer2->SetFileName( fnaam );
try
{
writer2->Update();
}
catch( itk::ExceptionObject & err )
{
logfile << "ExceptionObject caught 3" << endl;
logfile << err << endl;
}
logfile.close();
}
########################################################################
############################################
error message form Visual C++ 2005 Express Edition
########################################################################
############################################
------ Build started: Project: deform2, Configuration: Release Win32
------
Compiling...
deform2.cpp
.\deform2.cpp(351) : error C2039: 'SetTransform' : is not a member of
'itk::SpatialObjectWriter<NDimensions>'
with
[
NDimensions=2
]
Build log was saved at
"file://d:\Igor\ITKapps\deform2\deform2.dir\Release\BuildLog.htm"
deform2 - 1 error(s), 0 warning(s)
========== Build: 0 succeeded, 1 failed, 1 up-to-date, 0 skipped
==========
De inhoud van dit bericht is vertrouwelijk en alleen bestemd voor de geadresseerde(n). Anderen dan de geadresseerde mogen geen gebruik maken van dit bericht, het openbaar maken of op enige wijze verspreiden of vermenigvuldigen. Het UMCG kan niet aansprakelijk gesteld worden voor een incomplete aankomst of vertraging van dit verzonden bericht.
The contents of this message are confidential and only intended for the eyes of the addressee(s). Others than the addressee(s) are not allowed to use this message, to make it public or to distribute or multiply this message in any way. The UMCG cannot be held responsible for incomplete reception or delay of this transferred message.
More information about the Insight-users
mailing list