<br>Hi Toni,<br><br>How big are these images ?<br><br>Do you have a computer with enough RAM to load these images twice ?<br><br>The process may be simpler if you don&#39;t have memory constrains.<br><br>Several comments.<br>
<br>1) It doesn&#39;t seem that the TileImageFilter is what you need here.<br>    The PastImageFilter would look like a better match.<br><br>2) Have you verified the registration between the first two slices ?<br>    It is a good idea to solve the registration process in small<br>
    chunks of work. You are certainly doing a lot in this code, and<br>    may have skipped to fine tune the intermediate components.<br><br><br><br>    Regards,<br><br><br>         Luis<br><br><br>---------------------------------------------------------------------------------------------------<br>
<div class="gmail_quote">On Tue, Jul 14, 2009 at 11:36 AM, Toni Uusimäki <span dir="ltr">&lt;<a href="mailto:tonimuusimaki@gmail.com">tonimuusimaki@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<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 &quot;something&quot;. 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 &quot;itkImageFileReader.h&quot;<br>
#include &quot;itkImageFileWriter.h&quot;<br>
#include &quot;itkExtractImageFilter.h&quot;<br>
#include &quot;itkImage.h&quot;<br>
#include &quot;itkImageRegistrationMethod.h&quot;<br>
#include &quot;itkMeanSquaresImageToImageMetric.h&quot;<br>
#include &quot;itkLinearInterpolateImageFunction.h&quot;<br>
#include &quot;itkRegularStepGradientDescentOptimizer.h&quot;<br>
#include &quot;itkCenteredRigid2DTransform.h&quot;<br>
#include &quot;itkCenteredTransformInitializer.h&quot;<br>
#include &quot;itkResampleImageFilter.h&quot;<br>
#include &quot;itkCastImageFilter.h&quot;<br>
#include &quot;itkRescaleIntensityImageFilter.h&quot;<br>
#include &quot;itkCommand.h&quot;<br>
#include &lt;itkTileImageFilter.h&gt;<br>
class CommandIterationUpdate : public itk::Command <br>
{<br>
public:<br>
  typedef  CommandIterationUpdate   Self;<br>
  typedef  itk::Command             Superclass;<br>
  typedef itk::SmartPointer&lt;Self&gt;  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 &amp; event)<br>
    {<br>
      Execute( (const itk::Object *)caller, event);<br>
    }<br>
<br>
  void Execute(const itk::Object * object, const itk::EventObject &amp; event)<br>
    {<br>
      OptimizerPointer optimizer = <br>
        dynamic_cast&lt; OptimizerPointer &gt;( object );<br>
      if( ! itk::IterationEvent().CheckEvent( &amp;event ) )<br>
        {<br>
        return;<br>
        }<br>
      std::cout &lt;&lt; optimizer-&gt;GetCurrentIteration() &lt;&lt; &quot;   &quot;;<br>
      std::cout &lt;&lt; optimizer-&gt;GetValue() &lt;&lt; &quot;   &quot;;<br>
      std::cout &lt;&lt; optimizer-&gt;GetCurrentPosition() &lt;&lt; std::endl;<br>
    }<br>
};<br>
<br>
<br>
int main( int argc, char ** argv )<br>
{<br>
<br>
  if( argc &lt; 3 )<br>
    {<br>
    std::cerr &lt;&lt; &quot;Usage: &quot; &lt;&lt; std::endl;<br>
    std::cerr &lt;&lt; argv[0] &lt;&lt; &quot; input3DImageFile  &quot; &lt;&lt; std::endl;<br>
    std::cerr &lt;&lt; &quot; sliceNumber &quot; &lt;&lt; std::endl;<br>
    return EXIT_FAILURE;<br>
    }<br>
<br>
<br>
<br>
  typedef signed short        InputPixelType;<br>
  <br>
<br>
  typedef itk::Image&lt; InputPixelType,  3 &gt;    InputImageType;<br>
  typedef itk::Image&lt; InputPixelType, 2 &gt;  FixedImageType;<br>
  typedef itk::Image&lt; InputPixelType, 2 &gt;  MovingImageType;<br>
<br>
  typedef itk::ImageFileReader&lt; InputImageType  &gt;  ReaderType;<br>
  typedef itk::ImageFileWriter&lt; FixedImageType &gt;  WriterType;<br>
  typedef itk::ImageFileWriter&lt; MovingImageType &gt;  WriterType1;<br>
<br>
  const char * inputFilename  = argv[1];<br>
  <br>
  ReaderType::Pointer reader = ReaderType::New();<br>
  reader-&gt;SetFileName( inputFilename  );<br>
<br>
  typedef itk::CenteredRigid2DTransform&lt; double &gt; TransformType;<br>
<br>
<br>
<br>
  typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;<br>
  typedef itk::MeanSquaresImageToImageMetric&lt; <br>
                                    FixedImageType, <br>
                                    MovingImageType &gt;    MetricType;<br>
  typedef itk:: LinearInterpolateImageFunction&lt; <br>
                                    MovingImageType,<br>
                                    double          &gt;    InterpolatorType;<br>
  typedef itk::ImageRegistrationMethod&lt; <br>
                                    FixedImageType, <br>
                                    MovingImageType &gt;    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-&gt;SetMetric(        metric        );<br>
  registration-&gt;SetOptimizer(     optimizer     );<br>
  registration-&gt;SetInterpolator(  interpolator  );<br>
<br>
  TransformType::Pointer  transform = TransformType::New();<br>
  registration-&gt;SetTransform( transform );<br>
 <br>
  typedef itk::ImageFileReader&lt; FixedImageType  &gt; FixedImageReaderType;<br>
  typedef itk::ImageFileReader&lt; MovingImageType &gt; MovingImageReaderType;<br>
  FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();<br>
  MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();<br>
  <br>
<br>
  typedef itk::ExtractImageFilter&lt; InputImageType, FixedImageType &gt; FilterType;<br>
  FilterType::Pointer filter = FilterType::New();<br>
  typedef itk::ExtractImageFilter&lt; InputImageType, MovingImageType &gt; FilterType1;<br>
  FilterType1::Pointer filter1 = FilterType1::New();<br>
<br>
  reader-&gt;Update();<br>
  InputImageType::RegionType inputRegion =<br>
           reader-&gt;GetOutput()-&gt;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&lt;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-&gt;SetExtractionRegion( desiredRegion );<br>
  filter1-&gt;SetExtractionRegion( desiredRegion1 );<br>
<br>
  filter-&gt;SetInput( reader-&gt;GetOutput() );<br>
  filter1-&gt;SetInput( reader-&gt;GetOutput() );<br>
  <br>
  registration-&gt;SetFixedImage(    filter-&gt;GetOutput()    );<br>
  registration-&gt;SetMovingImage(   filter1-&gt;GetOutput()   );<br>
  filter-&gt;Update();<br>
  filter1-&gt;Update();<br>
<br>
  registration-&gt;SetFixedImageRegion( <br>
  filter-&gt;GetOutput()-&gt;GetBufferedRegion() );<br>
<br>
 typedef itk::CenteredTransformInitializer&lt; <br>
                                    TransformType, <br>
                                    FixedImageType, <br>
                                    MovingImageType &gt;  TransformInitializerType;<br>
<br>
  TransformInitializerType::Pointer initializer = TransformInitializerType::New();<br>
  <br>
  initializer-&gt;SetTransform(   transform );<br>
  initializer-&gt;SetFixedImage(  filter-&gt;GetOutput() );<br>
  initializer-&gt;SetMovingImage( filter1-&gt;GetOutput() );<br>
<br>
  initializer-&gt;MomentsOn();<br>
 <br>
  initializer-&gt;InitializeTransform();<br>
<br>
  transform-&gt;SetAngle( 0.0 );<br>
<br>
  registration-&gt;SetInitialTransformParameters( transform-&gt;GetParameters() );<br>
 <br>
  typedef OptimizerType::ScalesType       OptimizerScalesType;<br>
  OptimizerScalesType optimizerScales( transform-&gt;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-&gt;SetScales( optimizerScales );<br>
<br>
  optimizer-&gt;SetMaximumStepLength( 0.1    ); <br>
  optimizer-&gt;SetMinimumStepLength( 0.001 );<br>
  optimizer-&gt;SetNumberOfIterations( 200 );<br>
<br>
 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br>
  optimizer-&gt;AddObserver( itk::IterationEvent(), observer );<br>
<br>
  try <br>
    { <br>
    registration-&gt;StartRegistration(); <br>
    } <br>
  catch( itk::ExceptionObject &amp; err ) <br>
    { <br>
    std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl; <br>
    std::cerr &lt;&lt; err &lt;&lt; std::endl; <br>
    return EXIT_FAILURE;<br>
    } <br>
  <br>
  OptimizerType::ParametersType finalParameters = <br>
                    registration-&gt;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-&gt;GetCurrentIteration();<br>
  const double bestValue = optimizer-&gt;GetValue();<br>
<br>
  const double finalAngleInDegrees = finalAngle * 180.0 / vnl_math::pi;<br>
<br>
  std::cout &lt;&lt; &quot;Result = &quot; &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Angle (radians) &quot; &lt;&lt; finalAngle  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Angle (degrees) &quot; &lt;&lt; finalAngleInDegrees  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Center X      = &quot; &lt;&lt; finalRotationCenterX  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Center Y      = &quot; &lt;&lt; finalRotationCenterY  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Translation X = &quot; &lt;&lt; finalTranslationX  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Translation Y = &quot; &lt;&lt; finalTranslationY  &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Iterations    = &quot; &lt;&lt; numberOfIterations &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot; Metric value  = &quot; &lt;&lt; bestValue          &lt;&lt; std::endl;<br>
<br>
  transform-&gt;SetParameters( finalParameters );<br>
<br>
  TransformType::MatrixType matrix = transform-&gt;GetRotationMatrix();<br>
  TransformType::OffsetType offset = transform-&gt;GetOffset();<br>
<br>
  std::cout &lt;&lt; &quot;Matrix = &quot; &lt;&lt; std::endl &lt;&lt; matrix &lt;&lt; std::endl;<br>
  std::cout &lt;&lt; &quot;Offset = &quot; &lt;&lt; std::endl &lt;&lt; offset &lt;&lt; std::endl;<br>
 <br>
  typedef itk::ResampleImageFilter&lt; <br>
                            MovingImageType, <br>
                            FixedImageType &gt;    ResampleFilterType;<br>
  TransformType::Pointer finalTransform = TransformType::New();<br>
  finalTransform-&gt;SetParameters( finalParameters );<br>
  finalTransform-&gt;SetFixedParameters( transform-&gt;GetFixedParameters() );<br>
<br>
  ResampleFilterType::Pointer resample = ResampleFilterType::New();<br>
<br>
  resample-&gt;SetTransform( finalTransform );<br>
  resample-&gt;SetInput( filter1-&gt;GetOutput() );<br>
<br>
  FixedImageType::Pointer fixedImage = filter-&gt;GetOutput();<br>
<br>
  resample-&gt;SetSize(    fixedImage-&gt;GetLargestPossibleRegion().GetSize() );<br>
  resample-&gt;SetOutputOrigin(  fixedImage-&gt;GetOrigin() );<br>
  resample-&gt;SetOutputSpacing( fixedImage-&gt;GetSpacing() );<br>
  resample-&gt;SetOutputDirection( fixedImage-&gt;GetDirection() );<br>
  resample-&gt;SetDefaultPixelValue( 100 );<br>
  <br>
  <br>
<br>
<br>
  typedef itk::TileImageFilter&lt; MovingImageType, InputImageType &gt;     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-&gt;SetLayout( layout );<br>
<br>
<br>
  resample-&gt;Update();<br>
  tileFilter-&gt;SetInput(resample-&gt;GetOutput());<br>
<br>
  //tileFilter-&gt;PushBackInput(someFilter-&gt;GetOutput()));<br>
  tileFilter-&gt;Update();<br>
<br>
<br>
  }<br>
  return EXIT_SUCCESS;<br>
}<br>
<br>_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at: <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
<br></blockquote></div><br>