[Insight-users] 3D->2D->3D problematics
Luis Ibanez
luis.ibanez at kitware.com
Tue Jul 14 12:19:51 EDT 2009
Hi Toni,
How big are these images ?
Do you have a computer with enough RAM to load these images twice ?
The process may be simpler if you don't have memory constrains.
Several comments.
1) It doesn't seem that the TileImageFilter is what you need here.
The PastImageFilter would look like a better match.
2) Have you verified the registration between the first two slices ?
It is a good idea to solve the registration process in small
chunks of work. You are certainly doing a lot in this code, and
may have skipped to fine tune the intermediate components.
Regards,
Luis
---------------------------------------------------------------------------------------------------
On Tue, Jul 14, 2009 at 11:36 AM, Toni Uusimäki <tonimuusimaki at gmail.com>wrote:
>
> 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;
> }
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090714/4053f8f8/attachment-0001.htm>
More information about the Insight-users
mailing list