[Insight-users] Re: My Registration problem : 3D / 2D mix of types
Luis Ibanez
luis . ibanez at kitware . com
Wed, 10 Dec 2003 11:04:35 -0500
This is a multi-part message in MIME format.
--------------090400020301000000040606
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
Hi Jessica,
Thanks for posting your code,
it made things much easier.
There were two major problems in your code:
A) This first one is just C++:
You were declaring the variables
sizeF,sizeM,startF,startM inside
the block of an "if()" statement.
Variables declared inside any "{}"
block are only valid in the scope
of the block. Therefore the variable
do not exist outside the brackets.
The solution is to simply declare
the variables outside the brackets
before the "if" statement, and then
inside the "{}" just do the assignment.
B) This is a conceptual problem:
You are dealing with 3D and 2D images
in tha same process.
You use
Dimension = 2
TotalDimension = 3
what I understand from your code is that
your real problem is 3D and you are using
a 2D shortcut to solve it.
You seem to have two 3D datasets to stitch
together, and in order to find the translation
transform that will align them with the proper
overlap you are first registering only two 2D
slices taken from the 3D datasets.
That looks fine. The only difficulty is that the
2D registration process that you are performing
produces a 2D translation transform, while you
need a 3D translation transform for stitching the
3D datasets. No big deal, you can simply create
a 3D translation transform, copy the first two
components from the 2D transform resulting from
the registration and set the third component to
zero.
Please find attached the modified version of your file.
Both changes (A) and (B) were made in the attached file.
It is now compiling fine.
NOTE that PNG and BMP are very poor formats for images.
spacing and origin are note represented correctly in
this two file formats. You may want to convert your images
to a format like MetaImage or Analyze where there is better
support for the technical information about the image.
I doubt very much that you microscope will put the right
pixel spacing in the PNG file. You may have to check the
correct values in millimeters or microns and introduce this
values in the file format of your choice. This is *fundamental*
for the registration process to work correctly.
Let us know if you have further questions,
Thanks
Luis
--------------------------
Jessica de Ryk wrote:
>
> 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 ?
>>
>>
>
>
>>
>
>
> ------------------------------------------------------------------------
>
> --------------------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)
--------------090400020301000000040606
Content-Type: text/plain;
name="Reg2ImageConstruct.cxx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
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
{
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 =
dynamic_cast< OptimizerPointer >( object );
if( typeid( event ) != 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;
}
// For MovingImageLocation (assume 20% overlap)
// 1 for moving image right of fixed image
// 2 for moving image above fixed image
// The types of images should be declared first.
const unsigned int Dimension = 2;
const unsigned int TotalDimension = 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::TranslationTransform< double, TotalDimension > TotalTransformType;
typedef itk::GradientDescentOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction<
InternalImageType,
double > InterpolatorType;
typedef itk::MutualInformationImageToImageMetric<
InternalImageType,
InternalImageType > MetricType;
TransformType::Pointer transform = TransformType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
MetricType::Pointer metric = MetricType::New();
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;
// 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 = FixedROIFilterType::New();
MovingROIFilterType::Pointer movingROI = MovingROIFilterType::New();
ROIFixedImageType::IndexType startF;
ROIFixedImageType::SizeType sizeF;
ROIMovingImageType::IndexType startM;
ROIMovingImageType::SizeType sizeM;
// If moving image is to right of fixed image
if(atoi(argv[3]) == 1)
{
startF[0] = 1040;
startF[1] = 0;
sizeF[0] = 260;
sizeF[1] = 1030;
startM[0] = 0;
startM[1] = 0;
sizeM[0] = 260;
sizeM[1] = 1030;
}
// If moving image is above fixed image
if(atoi(argv[3]) == 2)
{
startF[0] = 0;
startF[1] = 824;
sizeF[0] = 1300;
sizeF[1] = 206;
startM[0] = 0;
startM[1] = 0;
sizeM[0] = 1300;
sizeM[1] = 206;
}
ROIFixedImageType::RegionType desiredFixedRegion;
desiredFixedRegion.SetIndex( startF );
desiredFixedRegion.SetSize( sizeF );
ROIMovingImageType::RegionType desiredMovingRegion;
desiredMovingRegion.SetIndex( startM );
desiredMovingRegion.SetSize( sizeM );
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = 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<
ROIFixedImageType, InternalImageType > FixedNormalizeFilterType;
typedef itk::NormalizeImageFilter<
ROIMovingImageType, InternalImageType > MovingNormalizeFilterType;
FixedNormalizeFilterType::Pointer fixedNormalizer = FixedNormalizeFilterType::New();
MovingNormalizeFilterType::Pointer movingNormalizer = MovingNormalizeFilterType::New();
// The output of the region of interest filters is connected as input to the
// 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(
fixedNormalizer->GetOutput()->GetBufferedRegion() );
metric->Initialize();
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
const double* fixedSpacing = fixedImage->GetSpacing();
const double* fixedOrigin = fixedImage->GetOrigin();
FixedImageType::SizeType fixedSize = fixedImage->GetLargestPossibleRegion().GetSize();
MetricType::TransformParametersType transformParameters( 2 );
transformParameters[0] = 0; //translation in x
transformParameters[1] = 0; //translation in y
double best_translation_x = transformParameters[0];
double best_translation_y = transformParameters[1];
double best_mi_value = metric->GetValue( transformParameters );
std::ofstream file_stream("Iterations.txt");
for ( double x = 41.0; x <= 42.0; x += 0.125 ) {
for ( double y = -2; y <=0.0; y += 0.125 ) {
transformParameters[0] = x;
transformParameters[1] = y;
double sum = 0.0;
for ( int i = 0; i < 100; i++ ) {
double mi_value = metric->GetValue( transformParameters );
sum += mi_value;
#if 0
if ( mi_value > best_mi_value ) {
best_translation_x = x;
best_translation_y = y;
best_mi_value = mi_value;
}
#endif
}
double avg = sum / 100.0;
std::cout << x << " " << y << " " << avg << '\n';
file_stream << x << " " << y << " " << avg << '\n';
if ( avg > best_mi_value ) {
best_mi_value = avg;
best_translation_x = x;
best_translation_y = y;
}
}
}
double TranslationAlongX = best_translation_x;
double TranslationAlongY = best_translation_y;
double bestValue = best_mi_value;
// Print out results
std::cout << "# Result = " << std::endl;
std::cout << "# Translation X = " << TranslationAlongX << std::endl;
std::cout << "# Translation Y = " << TranslationAlongY << std::endl;
std::cout << "# Metric value = " << bestValue << std::endl;
file_stream << "# Translation X = " << TranslationAlongX <<'\n';
file_stream << "# Translation Y = " << TranslationAlongY <<'\n';
file_stream << "# Metric value = " << bestValue;
file_stream.close();
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
transformParameters[0] = best_translation_x;
transformParameters[1] = best_translation_y;
finalTransform->SetParameters( transformParameters );
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( movingImageReader->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, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
typedef itk::SquaredDifferenceImageFilter<
FixedImageType,
FixedImageType,
OutputImageType > DifferenceFilterType;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput( difference->GetOutput() );
if( argc >= 5 )
{
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resample->GetOutput() );
writer2->SetFileName( argv[4] );
writer2->Update();
}
// Compute the difference image between the
// fixed and moving image before registration.
if( argc >= 6 )
{
writer2->SetFileName( argv[5] );
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( movingImageReader->GetOutput() );
writer2->Update();
}
//to generate the stitched color images
if( argc>= 7)
{
typedef itk::ImageFileReader< TotalFixedImageType > TotalFixedImageReaderType;
typedef itk::ImageFileReader< TotalMovingImageType > TotalMovingImageReaderType;
TotalFixedImageReaderType::Pointer TotalfixedImageReader = TotalFixedImageReaderType::New();
TotalMovingImageReaderType::Pointer TotalmovingImageReader = TotalMovingImageReaderType::New();
TotalfixedImageReader->SetFileName( argv[6] );
TotalmovingImageReader->SetFileName( argv[7] );
typedef itk::ResampleImageFilter<
TotalMovingImageType,
TotalFixedImageType > ResampleFilterType;
TotalTransformType::Pointer totalFinalTransform = TotalTransformType::New();
TotalTransformType::ParametersType totalTransformParameters(3);
totalTransformParameters[0] = best_translation_x + 260; //translate transform plus amount of overlap
totalTransformParameters[1] = best_translation_y;
totalTransformParameters[2] = 0;
totalFinalTransform->SetParameters( totalTransformParameters );
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( totalFinalTransform );
resample->SetInput( TotalmovingImageReader->GetOutput() );
TotalFixedImageType::ConstPointer TotalfixedImage = TotalfixedImageReader->GetOutput();
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;
typedef itk::CastImageFilter<
TotalFixedImageType,
TotalOutputImageType > CastFilterType;
typedef itk::ImageFileWriter< TotalOutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[8] );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
}
return 0;
}
--------------090400020301000000040606--