[Insight-users] 3D->2D->3D problematics
Toni Uusimäki
tonimuusimaki at gmail.com
Tue Jul 14 11:36:11 EDT 2009
Dear ITK Users,
i have tried for a long time now to connect certain example scripts
(ImageReadExtractWrite.cxx, ImageRegistration5.cxx) without any success.
I have an imagestack of electron tomography projections which i have to
align. So i set up a for loop to extract two sequential images and
pass them to the registration as fixed and moving image from the
extractfilter. Then after the registration i want to rewrite the moving
image from the stack
with the transformed one and use this as a fixed image to the second
registration and so on... Below is the tryout to combine the above examples
with tileimagefilter. I am
pretty new to ITK so it has propably lots of errors on it but the building
worked with visual c++ 2008 express and the program is doing "something".
Please give
me some feedback if this is a plausible way of doing this and if so, where i
am going wrong!!!
Cheers and many thanks for any reply!!!
p.s. the imagestack is a signed 16bit .tif file saved with Digital
Micrograph and it gives me some errors about the header. Is there some way
to read .mrc formats?
I also wanted to thank Luis and Sharath for the help on the last
issue i had but found it difficult to send reply messages to the same
topic...
T: Don I
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkExtractImageFilter.h"
#include "itkImage.h"
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkCenteredRigid2DTransform.h"
#include "itkCenteredTransformInitializer.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkCommand.h"
#include <itkTileImageFilter.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::RegularStepGradientDescentOptimizer 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;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
};
int main( int argc, char ** argv )
{
if( argc < 3 )
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " input3DImageFile " << std::endl;
std::cerr << " sliceNumber " << std::endl;
return EXIT_FAILURE;
}
typedef signed short InputPixelType;
typedef itk::Image< InputPixelType, 3 > InputImageType;
typedef itk::Image< InputPixelType, 2 > FixedImageType;
typedef itk::Image< InputPixelType, 2 > MovingImageType;
typedef itk::ImageFileReader< InputImageType > ReaderType;
typedef itk::ImageFileWriter< FixedImageType > WriterType;
typedef itk::ImageFileWriter< MovingImageType > WriterType1;
const char * inputFilename = argv[1];
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( inputFilename );
typedef itk::CenteredRigid2DTransform< double > TransformType;
typedef itk::RegularStepGradientDescentOptimizer 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 );
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
typedef itk::ExtractImageFilter< InputImageType, FixedImageType >
FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::ExtractImageFilter< InputImageType, MovingImageType >
FilterType1;
FilterType1::Pointer filter1 = FilterType1::New();
reader->Update();
InputImageType::RegionType inputRegion =
reader->GetOutput()->GetLargestPossibleRegion();
InputImageType::SizeType size = inputRegion.GetSize();
size[2] = 0;
const unsigned int sliceNumber = atoi( argv[2] );
InputImageType::IndexType start = inputRegion.GetIndex();
InputImageType::IndexType start1 = inputRegion.GetIndex();
for(int i=0;i<sliceNumber-1;i++)
{
start[2] = i;
start1[2] = i+1;
InputImageType::RegionType desiredRegion;
desiredRegion.SetSize( size );
desiredRegion.SetIndex( start );
InputImageType::RegionType desiredRegion1;
desiredRegion1.SetSize( size );
desiredRegion1.SetIndex( start1 );
filter->SetExtractionRegion( desiredRegion );
filter1->SetExtractionRegion( desiredRegion1 );
filter->SetInput( reader->GetOutput() );
filter1->SetInput( reader->GetOutput() );
registration->SetFixedImage( filter->GetOutput() );
registration->SetMovingImage( filter1->GetOutput() );
filter->Update();
filter1->Update();
registration->SetFixedImageRegion(
filter->GetOutput()->GetBufferedRegion() );
typedef itk::CenteredTransformInitializer<
TransformType,
FixedImageType,
MovingImageType >
TransformInitializerType;
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform( transform );
initializer->SetFixedImage( filter->GetOutput() );
initializer->SetMovingImage( filter1->GetOutput() );
initializer->MomentsOn();
initializer->InitializeTransform();
transform->SetAngle( 0.0 );
registration->SetInitialTransformParameters( transform->GetParameters() );
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = translationScale;
optimizerScales[2] = translationScale;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizer->SetScales( optimizerScales );
optimizer->SetMaximumStepLength( 0.1 );
optimizer->SetMinimumStepLength( 0.001 );
optimizer->SetNumberOfIterations( 200 );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
const double finalAngle = finalParameters[0];
const double finalRotationCenterX = finalParameters[1];
const double finalRotationCenterY = finalParameters[2];
const double finalTranslationX = finalParameters[3];
const double finalTranslationY = finalParameters[4];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
const double finalAngleInDegrees = finalAngle * 180.0 / vnl_math::pi;
std::cout << "Result = " << std::endl;
std::cout << " Angle (radians) " << finalAngle << std::endl;
std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;
std::cout << " Center X = " << finalRotationCenterX << std::endl;
std::cout << " Center Y = " << finalRotationCenterY << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
transform->SetParameters( finalParameters );
TransformType::MatrixType matrix = transform->GetRotationMatrix();
TransformType::OffsetType offset = transform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters( finalParameters );
finalTransform->SetFixedParameters( transform->GetFixedParameters() );
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( filter1->GetOutput() );
FixedImageType::Pointer fixedImage = filter->GetOutput();
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetOutputDirection( fixedImage->GetDirection() );
resample->SetDefaultPixelValue( 100 );
typedef itk::TileImageFilter< MovingImageType, InputImageType >
TilerType;
TilerType::Pointer tileFilter = TilerType::New();
TilerType::LayoutArrayType layout;
layout[0] = 1;
layout[1] = 1;
layout[2] = 0;
tileFilter->SetLayout( layout );
resample->Update();
tileFilter->SetInput(resample->GetOutput());
//tileFilter->PushBackInput(someFilter->GetOutput()));
tileFilter->Update();
}
return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090714/211930f6/attachment-0001.htm>
More information about the Insight-users
mailing list