[Insight-users] Re:My Registration problem
Jessica de Ryk
jessica-deryk at uiowa . edu
Tue, 9 Dec 2003 14:45:44 -0600
This is a multi-part message in MIME format.
------=_NextPart_000_00AB_01C3BE63.20256FF0
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
Attached is the code I am having trouble with and the error file. I wish to
take two grayscale images which share a 20% overlap region, segment out the
common overlaping region, register the overlaped sections and apply the
transform to the original color images such that they are stitched together
to form one image. As mentioned previously I have tested the sections
seperately and the ImageCrop filter produced a blurr issue do to the spacing
(The images are originally captured from the microscope+CCD as Bitmaps and I
convert them to PNG format, I am currently altering the capture to PNG
format so hopefully that will help with the spacing issue)
For the final joining of the two color images, I am currently reverting to
MATLAB as I get errors about converting from two to three dimensions (Can
resample filter only be applied to grayscale images or am I not initializing
something correctly?)
Thanks
Jessica
----- Original Message -----
From: "Luis Ibanez" <luis . ibanez at kitware . com>
To: "Jessica de Ryk" <jessica-deryk at uiowa . edu>
Cc: <Insight-users at public . kitware . com>
Sent: Monday, December 08, 2003 9:47 PM
Subject: Re: [Insight-users] An overlap invariant mutual information metric
>
> Hi Jessica,
>
> You declarations for the RegionOfInterestImageFilter
> look fine.
>
> You don't have any variable named sizeF, or startF though...
>
> Are you still having trouble with this code ?
>
> If so, could you please post the current file along
> with the error messages ?
>
>
>
>
------=_NextPart_000_00AB_01C3BE63.20256FF0
Content-Type: application/octet-stream;
name="Reg2ImageConstruct.cxx"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="Reg2ImageConstruct.cxx"
//MyImageRegistration.cxx
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkGradientDescentOptimizer.h"
#include "itkImage.h"
#include "itkNormalizeImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRegionOfInterestImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include <sstream>
//
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
//
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command=20
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
//typedef itk::GradientDescentOptimizer 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)
{
#if 0
OptimizerPointer optimizer =3D=20
dynamic_cast< OptimizerPointer >( object );
if( typeid( event ) !=3D typeid( itk::IterationEvent ) )
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition()[0] << " ";
std::cout << optimizer->GetCurrentPosition()[1] << std::endl;
#endif
}
};
int main( int argc, char **argv )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile MovingImageLocation";
std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] =
";
std::cerr << " [TotalfixedImage] [TotalMovingImage] ";
std::cerr << " [TotalOutputImage] "<< std::endl;
return 1;
}
=20
// For MovingImageLocation (assume 20% overlap)
// 1 for moving image right of fixed image
// 2 for moving image above fixed image
=20
// The types of images should be declared first.=20
=20
const unsigned int Dimension =3D 2;
const unsigned int TotalDimension =3D 3;
typedef unsigned short PixelType;
typedef float InternalPixelType;
typedef float ROIPixelType;
typedef float TotalPixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
typedef itk::Image< ROIPixelType, Dimension > ROIFixedImageType;
typedef itk::Image< ROIPixelType, Dimension > ROIMovingImageType;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
typedef itk::Image< TotalPixelType, TotalDimension > =
TotalFixedImageType;
typedef itk::Image< TotalPixelType, TotalDimension > =
TotalMovingImageType;
typedef itk::TranslationTransform< double, Dimension > TransformType;
typedef itk::GradientDescentOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction<=20
InternalImageType,
double > =
InterpolatorType;
typedef itk::MutualInformationImageToImageMetric<=20
InternalImageType,=20
InternalImageType > =
MetricType;
TransformType::Pointer transform =3D TransformType::New();
InterpolatorType::Pointer interpolator =3D InterpolatorType::New();
MetricType::Pointer metric =3D MetricType::New();
=20
metric->SetFixedImageStandardDeviation( 0.4 );
metric->SetMovingImageStandardDeviation( 0.4 );
metric->SetNumberOfSpatialSamples( 100 );
metric->SetInterpolator( interpolator );
metric->SetTransform( transform );
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
=20
// Settings for the region of interest filter for the moving and fixed =
images
typedef itk::RegionOfInterestImageFilter< FixedImageType, =
ROIFixedImageType > FixedROIFilterType;
typedef itk::RegionOfInterestImageFilter< MovingImageType, =
ROIMovingImageType > MovingROIFilterType;
FixedROIFilterType::Pointer fixedROI =3D FixedROIFilterType::New();
MovingROIFilterType::Pointer movingROI =3D =
MovingROIFilterType::New();
// If moving image is to right of fixed image
if(atoi(argv[3]) =3D=3D 1){
ROIFixedImageType::IndexType startF;
startF[0] =3D 1040; =20
startF[1] =3D 0;
=20
ROIFixedImageType::SizeType sizeF;
sizeF[0] =3D 260; =20
sizeF[1] =3D 1030;
=20
ROIMovingImageType::IndexType startM;
startM[0] =3D 0;=20
startM[1] =3D 0;
ROIMovingImageType::SizeType sizeM;
sizeM[0] =3D 260; =20
sizeM[1] =3D 1030;
}
/*
// If moving image is above fixed image
if(atoi(argv[3]) =3D=3D 2){
ROIFixedImageType::IndexType startF;
startF[0] =3D 0; =20
startF[1] =3D 824;
ROIFixedImageType::SizeType sizeF;
sizeF[0] =3D 1300; =20
sizeF[1] =3D 206; =20
=20
ROIMovingImageType::IndexType startM;
startM[0] =3D 0;=20
startM[1] =3D 0;
ROIMovingImageType::SizeType sizeM;
sizeM[0] =3D 1300; =20
sizeM[1] =3D 206;=20
}
*/
ROIFixedImageType::RegionType desiredFixedRegion;
desiredFixedRegion.SetIndex( startF );
desiredFixedRegion.SetSize( sizeF );
=20
ROIMovingImageType::RegionType desiredMovingRegion;
desiredMovingRegion.SetIndex( startM );
desiredMovingRegion.SetSize( sizeM );
=20
FixedImageReaderType::Pointer fixedImageReader =3D =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =3D =
MovingImageReaderType::New();
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
fixedROI->SetInput( fixedImageReader->GetOutput() );
movingROI->SetInput( movingImageReader->GetOutput() );
// The normalization filters are declared using the fixed and moving =
image
// types as input and the internal image type as output.
typedef itk::NormalizeImageFilter<=20
ROIFixedImageType, InternalImageType > =
FixedNormalizeFilterType;
=20
typedef itk::NormalizeImageFilter<=20
ROIMovingImageType, InternalImageType > =
MovingNormalizeFilterType;
FixedNormalizeFilterType::Pointer fixedNormalizer =3D =
FixedNormalizeFilterType::New();
MovingNormalizeFilterType::Pointer movingNormalizer =3D =
MovingNormalizeFilterType::New();
// The output of the region of interest filters is connected as input =
to the=20
// normalization filters. The inputs to the registration method are =
taken
// from the normalization filters.
fixedNormalizer->SetInput( fixedROI->GetOutput() );
movingNormalizer->SetInput( movingROI->GetOutput() );
fixedNormalizer->Update();
metric->SetFixedImage( fixedNormalizer->GetOutput() );
metric->SetMovingImage( movingNormalizer->GetOutput() );
metric->SetFixedImageRegion(=20
fixedNormalizer->GetOutput()->GetBufferedRegion() );
=20
metric->Initialize();
FixedImageType::Pointer fixedImage =3D fixedImageReader->GetOutput();
const double* fixedSpacing =3D fixedImage->GetSpacing();
const double* fixedOrigin =3D fixedImage->GetOrigin();
FixedImageType::SizeType fixedSize =3D =
fixedImage->GetLargestPossibleRegion().GetSize();
=20
MetricType::TransformParametersType transformParameters( 2 );
transformParameters[0] =3D 0; //translation in x
transformParameters[1] =3D 0; //translation in y
double best_translation_x =3D transformParameters[0];
double best_translation_y =3D transformParameters[1];
double best_mi_value =3D metric->GetValue( transformParameters );
std::ofstream file_stream("Iterations.txt");
for ( double x =3D 41.0; x <=3D 42.0; x +=3D 0.125 ) {=20
for ( double y =3D -2; y <=3D0.0; y +=3D 0.125 ) {
transformParameters[0] =3D x;
transformParameters[1] =3D y;
double sum =3D 0.0;
for ( int i =3D 0; i < 100; i++ ) {
double mi_value =3D metric->GetValue( transformParameters );
sum +=3D mi_value;
#if 0
if ( mi_value > best_mi_value ) {
best_translation_x =3D x;
best_translation_y =3D y;
best_mi_value =3D mi_value;
}
#endif
}
double avg =3D sum / 100.0;
std::cout << x << " " << y << " " << avg << '\n';
file_stream << x << " " << y << " " << avg << '\n';
=09
if ( avg > best_mi_value ) {
best_mi_value =3D avg;
best_translation_x =3D x;
best_translation_y =3D y;
}
}
}
=09
=20
double TranslationAlongX =3D best_translation_x;
double TranslationAlongY =3D best_translation_y;
double bestValue =3D best_mi_value;
// Print out results
std::cout << "# Result =3D " << std::endl;
std::cout << "# Translation X =3D " << TranslationAlongX << =
std::endl;
std::cout << "# Translation Y =3D " << TranslationAlongY << =
std::endl;
std::cout << "# Metric value =3D " << bestValue << =
std::endl;
file_stream << "# Translation X =3D " << TranslationAlongX <<'\n';
file_stream << "# Translation Y =3D " << TranslationAlongY <<'\n';
file_stream << "# Metric value =3D " << bestValue;
file_stream.close();
typedef itk::ResampleImageFilter<=20
MovingImageType,=20
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform =3D TransformType::New();
transformParameters[0] =3D best_translation_x;
transformParameters[1] =3D best_translation_y;
finalTransform->SetParameters( transformParameters );
ResampleFilterType::Pointer resample =3D ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( movingImageReader->GetOutput() );
=20
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() =
);
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
=20
typedef itk::CastImageFilter<=20
FixedImageType,
OutputImageType > CastFilterType;
=20
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer =3D WriterType::New();
CastFilterType::Pointer caster =3D CastFilterType::New();
writer->SetFileName( argv[3] );
=20
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
typedef itk::SquaredDifferenceImageFilter<=20
FixedImageType,=20
FixedImageType,=20
OutputImageType > =
DifferenceFilterType;
DifferenceFilterType::Pointer difference =3D =
DifferenceFilterType::New();
WriterType::Pointer writer2 =3D WriterType::New();
writer2->SetInput( difference->GetOutput() ); =20
=20
if( argc >=3D 5 )
{
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resample->GetOutput() );
writer2->SetFileName( argv[4] );
writer2->Update();
}
// Compute the difference image between the=20
// fixed and moving image before registration.
if( argc >=3D 6 )
{
writer2->SetFileName( argv[5] );
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( movingImageReader->GetOutput() );
writer2->Update();
}
=20
//to generate the stitched color images
if( argc>=3D 7)
{
typedef itk::ImageFileReader< TotalFixedImageType > =
TotalFixedImageReaderType;
typedef itk::ImageFileReader< TotalMovingImageType > =
TotalMovingImageReaderType;
=20
TotalFixedImageReaderType::Pointer TotalfixedImageReader =3D =
TotalFixedImageReaderType::New();
TotalMovingImageReaderType::Pointer TotalmovingImageReader =3D =
TotalMovingImageReaderType::New();
=20
TotalfixedImageReader->SetFileName( argv[6] );
TotalmovingImageReader->SetFileName( argv[7] );
typedef itk::ResampleImageFilter<=20
TotalMovingImageType,=20
TotalFixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform =3D TransformType::New();
transformParameters[0] =3D best_translation_x + 260; //translate =
transform plus amount of overlap
transformParameters[1] =3D best_translation_y;
finalTransform->SetParameters( transformParameters );
ResampleFilterType::Pointer resample =3D ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( TotalmovingImageReader->GetOutput() );
=20
FixedImageType::Pointer TotalfixedImage =3D =
TotalfixedImageReader->GetOutput();
=20
=20
resample->SetSize( =
TotalfixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( TotalfixedImage->GetOrigin() );
resample->SetOutputSpacing( TotalfixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char TotalOutputPixelType;
typedef itk::Image< TotalOutputPixelType, TotalDimension > =
TotalOutputImageType;
=20
typedef itk::CastImageFilter<=20
TotalFixedImageType,
TotalOutputImageType > CastFilterType;
=20
typedef itk::ImageFileWriter< TotalOutputImageType > WriterType;
WriterType::Pointer writer =3D WriterType::New();
CastFilterType::Pointer caster =3D CastFilterType::New();
writer->SetFileName( argv[8] );
=20
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
}
return 0;
}
------=_NextPart_000_00AB_01C3BE63.20256FF0
Content-Type: text/plain;
name="Output - Build.txt"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="Output - Build.txt"
--------------------Configuration: Reg2ImageConstruct - Win32 =
Debug--------------------
Compiling...
Reg2ImageConstruct.cxx
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(179) =
: error C2065: 'startF' : undeclared identifier
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(180) =
: error C2065: 'sizeF' : undeclared identifier
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(184) =
: error C2065: 'startM' : undeclared identifier
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(185) =
: error C2065: 'sizeM' : undeclared identifier
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(400) =
: error C2664: 'SetTransform' : cannot convert parameter 1 from 'class =
itk::SmartPointer<class itk::TranslationTransform<double,2> >' to 'class =
itk::Transform<double,3,3> *'
No user-defined-conversion operator available that can perform =
this conversion, or the operator cannot be called
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(403) =
: error C2440: 'initializing' : cannot convert from 'class =
itk::Image<float,3> *' to 'class itk::SmartPointer<class =
itk::Image<unsigned short,2> >'
No constructor could take the source type, or constructor =
overload resolution was ambiguous
D:\Projects\Registration\TotalImageConstruct\Reg2ImageConstruct.cxx(406) =
: error C2664: 'SetSize' : cannot convert parameter 1 from 'const class =
itk::Size<2>' to 'const class itk::Size<3>'
No constructor could take the source type, or constructor =
overload resolution was ambiguous
Error executing cl.exe.
ALL_BUILD - 7 error(s), 0 warning(s)
------=_NextPart_000_00AB_01C3BE63.20256FF0--