<br>Dear ITK Users,<br><br>i have tried for a long time now to connect certain example scripts (ImageReadExtractWrite.cxx, ImageRegistration5.cxx) without any success.<br>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<br>
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<br>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<br>
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<br>me some feedback if this is a plausible way of doing this and if so, where i am going wrong!!!<br>
<br>Cheers and many thanks for any reply!!!<br><br>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?<br> 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...<br>
<br>T: Don I<br><br><br>
#if defined(_MSC_VER)<br>
#pragma warning ( disable : 4786 )<br>
#endif<br>
<br>
#ifdef __BORLANDC__<br>
#define ITK_LEAN_AND_MEAN<br>
#endif<br>
<br>
#include "itkImageFileReader.h"<br>
#include "itkImageFileWriter.h"<br>
#include "itkExtractImageFilter.h"<br>
#include "itkImage.h"<br>
#include "itkImageRegistrationMethod.h"<br>
#include "itkMeanSquaresImageToImageMetric.h"<br>
#include "itkLinearInterpolateImageFunction.h"<br>
#include "itkRegularStepGradientDescentOptimizer.h"<br>
#include "itkCenteredRigid2DTransform.h"<br>
#include "itkCenteredTransformInitializer.h"<br>
#include "itkResampleImageFilter.h"<br>
#include "itkCastImageFilter.h"<br>
#include "itkRescaleIntensityImageFilter.h"<br>
#include "itkCommand.h"<br>
#include <itkTileImageFilter.h><br>
class CommandIterationUpdate : public itk::Command <br>
{<br>
public:<br>
typedef CommandIterationUpdate Self;<br>
typedef itk::Command Superclass;<br>
typedef itk::SmartPointer<Self> Pointer;<br>
itkNewMacro( Self );<br>
protected:<br>
CommandIterationUpdate() {};<br>
public:<br>
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br>
typedef const OptimizerType * OptimizerPointer;<br>
<br>
void Execute(itk::Object *caller, const itk::EventObject & event)<br>
{<br>
Execute( (const itk::Object *)caller, event);<br>
}<br>
<br>
void Execute(const itk::Object * object, const itk::EventObject & event)<br>
{<br>
OptimizerPointer optimizer = <br>
dynamic_cast< OptimizerPointer >( object );<br>
if( ! itk::IterationEvent().CheckEvent( &event ) )<br>
{<br>
return;<br>
}<br>
std::cout << optimizer->GetCurrentIteration() << " ";<br>
std::cout << optimizer->GetValue() << " ";<br>
std::cout << optimizer->GetCurrentPosition() << std::endl;<br>
}<br>
};<br>
<br>
<br>
int main( int argc, char ** argv )<br>
{<br>
<br>
if( argc < 3 )<br>
{<br>
std::cerr << "Usage: " << std::endl;<br>
std::cerr << argv[0] << " input3DImageFile " << std::endl;<br>
std::cerr << " sliceNumber " << std::endl;<br>
return EXIT_FAILURE;<br>
}<br>
<br>
<br>
<br>
typedef signed short InputPixelType;<br>
<br>
<br>
typedef itk::Image< InputPixelType, 3 > InputImageType;<br>
typedef itk::Image< InputPixelType, 2 > FixedImageType;<br>
typedef itk::Image< InputPixelType, 2 > MovingImageType;<br>
<br>
typedef itk::ImageFileReader< InputImageType > ReaderType;<br>
typedef itk::ImageFileWriter< FixedImageType > WriterType;<br>
typedef itk::ImageFileWriter< MovingImageType > WriterType1;<br>
<br>
const char * inputFilename = argv[1];<br>
<br>
ReaderType::Pointer reader = ReaderType::New();<br>
reader->SetFileName( inputFilename );<br>
<br>
typedef itk::CenteredRigid2DTransform< double > TransformType;<br>
<br>
<br>
<br>
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br>
typedef itk::MeanSquaresImageToImageMetric< <br>
FixedImageType, <br>
MovingImageType > MetricType;<br>
typedef itk:: LinearInterpolateImageFunction< <br>
MovingImageType,<br>
double > InterpolatorType;<br>
typedef itk::ImageRegistrationMethod< <br>
FixedImageType, <br>
MovingImageType > RegistrationType;<br>
<br>
MetricType::Pointer metric = MetricType::New();<br>
OptimizerType::Pointer optimizer = OptimizerType::New();<br>
InterpolatorType::Pointer interpolator = InterpolatorType::New();<br>
RegistrationType::Pointer registration = RegistrationType::New();<br>
<br>
<br>
registration->SetMetric( metric );<br>
registration->SetOptimizer( optimizer );<br>
registration->SetInterpolator( interpolator );<br>
<br>
TransformType::Pointer transform = TransformType::New();<br>
registration->SetTransform( transform );<br>
<br>
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<br>
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;<br>
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();<br>
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();<br>
<br>
<br>
typedef itk::ExtractImageFilter< InputImageType, FixedImageType > FilterType;<br>
FilterType::Pointer filter = FilterType::New();<br>
typedef itk::ExtractImageFilter< InputImageType, MovingImageType > FilterType1;<br>
FilterType1::Pointer filter1 = FilterType1::New();<br>
<br>
reader->Update();<br>
InputImageType::RegionType inputRegion =<br>
reader->GetOutput()->GetLargestPossibleRegion();<br>
<br>
InputImageType::SizeType size = inputRegion.GetSize();<br>
size[2] = 0;<br>
<br>
const unsigned int sliceNumber = atoi( argv[2] );<br>
<br>
InputImageType::IndexType start = inputRegion.GetIndex();<br>
InputImageType::IndexType start1 = inputRegion.GetIndex();<br>
<br>
for(int i=0;i<sliceNumber-1;i++)<br>
{<br>
<br>
<br>
start[2] = i;<br>
start1[2] = i+1;<br>
<br>
InputImageType::RegionType desiredRegion;<br>
desiredRegion.SetSize( size );<br>
desiredRegion.SetIndex( start );<br>
InputImageType::RegionType desiredRegion1;<br>
desiredRegion1.SetSize( size );<br>
desiredRegion1.SetIndex( start1 );<br>
<br>
filter->SetExtractionRegion( desiredRegion );<br>
filter1->SetExtractionRegion( desiredRegion1 );<br>
<br>
filter->SetInput( reader->GetOutput() );<br>
filter1->SetInput( reader->GetOutput() );<br>
<br>
registration->SetFixedImage( filter->GetOutput() );<br>
registration->SetMovingImage( filter1->GetOutput() );<br>
filter->Update();<br>
filter1->Update();<br>
<br>
registration->SetFixedImageRegion( <br>
filter->GetOutput()->GetBufferedRegion() );<br>
<br>
typedef itk::CenteredTransformInitializer< <br>
TransformType, <br>
FixedImageType, <br>
MovingImageType > TransformInitializerType;<br>
<br>
TransformInitializerType::Pointer initializer = TransformInitializerType::New();<br>
<br>
initializer->SetTransform( transform );<br>
initializer->SetFixedImage( filter->GetOutput() );<br>
initializer->SetMovingImage( filter1->GetOutput() );<br>
<br>
initializer->MomentsOn();<br>
<br>
initializer->InitializeTransform();<br>
<br>
transform->SetAngle( 0.0 );<br>
<br>
registration->SetInitialTransformParameters( transform->GetParameters() );<br>
<br>
typedef OptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<br>
const double translationScale = 1.0 / 1000.0;<br>
<br>
optimizerScales[0] = 1.0;<br>
optimizerScales[1] = translationScale;<br>
optimizerScales[2] = translationScale;<br>
optimizerScales[3] = translationScale;<br>
optimizerScales[4] = translationScale;<br>
<br>
optimizer->SetScales( optimizerScales );<br>
<br>
optimizer->SetMaximumStepLength( 0.1 ); <br>
optimizer->SetMinimumStepLength( 0.001 );<br>
optimizer->SetNumberOfIterations( 200 );<br>
<br>
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br>
optimizer->AddObserver( itk::IterationEvent(), observer );<br>
<br>
try <br>
{ <br>
registration->StartRegistration(); <br>
} <br>
catch( itk::ExceptionObject & err ) <br>
{ <br>
std::cerr << "ExceptionObject caught !" << std::endl; <br>
std::cerr << err << std::endl; <br>
return EXIT_FAILURE;<br>
} <br>
<br>
OptimizerType::ParametersType finalParameters = <br>
registration->GetLastTransformParameters();<br>
<br>
<br>
const double finalAngle = finalParameters[0];<br>
const double finalRotationCenterX = finalParameters[1];<br>
const double finalRotationCenterY = finalParameters[2];<br>
const double finalTranslationX = finalParameters[3];<br>
const double finalTranslationY = finalParameters[4];<br>
<br>
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<br>
const double bestValue = optimizer->GetValue();<br>
<br>
const double finalAngleInDegrees = finalAngle * 180.0 / vnl_math::pi;<br>
<br>
std::cout << "Result = " << std::endl;<br>
std::cout << " Angle (radians) " << finalAngle << std::endl;<br>
std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;<br>
std::cout << " Center X = " << finalRotationCenterX << std::endl;<br>
std::cout << " Center Y = " << finalRotationCenterY << std::endl;<br>
std::cout << " Translation X = " << finalTranslationX << std::endl;<br>
std::cout << " Translation Y = " << finalTranslationY << std::endl;<br>
std::cout << " Iterations = " << numberOfIterations << std::endl;<br>
std::cout << " Metric value = " << bestValue << std::endl;<br>
<br>
transform->SetParameters( finalParameters );<br>
<br>
TransformType::MatrixType matrix = transform->GetRotationMatrix();<br>
TransformType::OffsetType offset = transform->GetOffset();<br>
<br>
std::cout << "Matrix = " << std::endl << matrix << std::endl;<br>
std::cout << "Offset = " << std::endl << offset << std::endl;<br>
<br>
typedef itk::ResampleImageFilter< <br>
MovingImageType, <br>
FixedImageType > ResampleFilterType;<br>
TransformType::Pointer finalTransform = TransformType::New();<br>
finalTransform->SetParameters( finalParameters );<br>
finalTransform->SetFixedParameters( transform->GetFixedParameters() );<br>
<br>
ResampleFilterType::Pointer resample = ResampleFilterType::New();<br>
<br>
resample->SetTransform( finalTransform );<br>
resample->SetInput( filter1->GetOutput() );<br>
<br>
FixedImageType::Pointer fixedImage = filter->GetOutput();<br>
<br>
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );<br>
resample->SetOutputOrigin( fixedImage->GetOrigin() );<br>
resample->SetOutputSpacing( fixedImage->GetSpacing() );<br>
resample->SetOutputDirection( fixedImage->GetDirection() );<br>
resample->SetDefaultPixelValue( 100 );<br>
<br>
<br>
<br>
<br>
typedef itk::TileImageFilter< MovingImageType, InputImageType > TilerType;<br>
<br>
TilerType::Pointer tileFilter = TilerType::New();<br>
<br>
TilerType::LayoutArrayType layout;<br>
<br>
layout[0] = 1;<br>
layout[1] = 1;<br>
layout[2] = 0;<br>
<br>
<br>
tileFilter->SetLayout( layout );<br>
<br>
<br>
resample->Update();<br>
tileFilter->SetInput(resample->GetOutput());<br>
<br>
//tileFilter->PushBackInput(someFilter->GetOutput()));<br>
tileFilter->Update();<br>
<br>
<br>
}<br>
return EXIT_SUCCESS;<br>
}<br>