[Insight-users] Attempt to rotate image resulting in black image.

[Mike Jackson]

Hi Luis,
  Thanks for the detailed explanation. In our case the origin does
correspond to a true physical space. The data comes from a Zeiss
Microscope which is collecting a mosaic of images. Yes, I can "stitch"
the images together just by simply using the coordinates stored in the
image data and turn that into a rather large mosaic (about 10,000
pixels on a side.).

  As far as your example code below I will try it out in the morning
(at home now) but that code does not really correspond to what is in
the itkManual just FYI. Seeming as that code does not seem to work
I'll definitely give yours a try (it also seems more intuitive) and
report back.

 Again, My deepest thanks for the detailed explanation.

Mike Jackson

On Wed, Jun 3, 2009 at 8:42 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
>
> Hi Michael,
>
>
> We define as "Origin" of an image the physical coordinates of its
> pixel with index [0,0,0].
>
> Rotations, are performed by default, with respect to the origin
> of the coordinate system (the point with physical coordinates
> (0.0, 0.0, 0.0), which doesn't have to coincide with the first
> pixel of the image, nor with the central pixel of the image.
>
> In ITK we assume that your image is a representation of some
> physical reality (for example a patient's body), and therefore
> all transformations should be performed in the context of the
> physical coordinates, not the grid of pixels.
>
> The purpose of the "Origin" is that if you use your physical
> image acquisition device (microscope, telescope, CT scanner,
> MRI scanner, ultrasound device, lidar...) and you acquire
> images of different sections of an object that has a physical
> manifestation in the real world, then you should be able to
> recreate a mosaic of such real object by simply placing the
> images in a common coordinate system according to the coordinates
> of their origin.
>
>
> ----
>
>
> In order to understand the behavior of the transformations,
> you should think of your image in the context of the physical
> coordinate system.
>
> If you want your image to rotate around the one of the pixel
> in one of the image corners the you should do the following:
>
>
>  A) Compute the Index of that pixel
>
>  B) Call image->TransformIndexToPhysicalPoint( index, point )
>     to compute the physical coordinates corresponding to this
>     pixel. This transformation will take into account the
>     image origin, pixel spacing and image orientation.
>
>  C) Call transform->SetCenter( point )
>     to tell the transform that you want to use this point
>     as the center of rotation.
>
>  D) Connect the transform to the resample filter.
>
>  E) call Update() in the resample filter.
>
>
>
> Please let us know if you still find any problem,
>
>
>    Thanks,
>
>
>       Luis
>
>
>
> ---------------------------
> Michael Jackson wrote:
>>
>> Just to follow up with this, If I force my origins to be (0,0) then  the
>> rotation seems to work. So I guess I am not understanding what  exactly the
>> origins are? I have reread through several sections of the  ITK guide and I
>> thought I had it figured out. The origins of where the  image is taken is
>> actually stored in our data files so I just thought  I would use that. The
>> values are in Microns and are typically around  the 48,000 to 50,000 range.
>>
>> Is there an implicit units (like mm) in vtkImage that I don't know  about?
>> I just can not get this figured out.
>>
>>  Explanations are truly appreciated at this point.
>>
>> _________________________________________________________
>> Mike Jackson                  mike.jackson at bluequartz.net
>> BlueQuartz Software                    www.bluequartz.net
>> Principal Software Engineer                  Dayton, Ohio
>>
>>
>>
>> On Jun 3, 2009, at 5:50 PM, Michael Jackson wrote:
>>
>>> I am trying to rotate an image about the "upper left" corner. I had  to
>>> write my own ImageSource derived class to read the data from and  HDF5 file.
>>> I am pretty sure the "Hdf5 reader" works as I can create  the ImageReader
>>> instance from it, read an image and write the image  to disk as a tiff file.
>>>
>>>  Where my code seems to be failing is in the application of an
>>>  AffineTransformFilter. I am basically taking the code straight from  one of
>>> the examples adapting it to use my ImageSource derived reader  and running
>>> it. And I get the dreaded "black image". Here is the  code and the output
>>> from the program. Oddly, if I switch the reader  back to the usual
>>> ImageFileReader and read up a tiff image I get the  correct result, so maybe
>>> my custom reader class is not quite correct?
>>>
>>>
>>> #define USE_MXA_FILE 1
>>> int main(int argc, char **argv) {
>>>  int exampleAction =0;
>>>  const     unsigned int   Dimension = 2;
>>>  typedef   unsigned char  InputPixelType;
>>>  typedef   unsigned char  OutputPixelType;
>>>
>>>  typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
>>>  typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
>>>
>>> #if USE_MXA_FILE
>>>  typedef itkR3DImageIO  ReaderType;
>>> #else
>>>  typedef itk::ImageFileReader< InputImageType >  ReaderType;
>>> #endif
>>>  typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>>>
>>>  ReaderType::Pointer reader = ReaderType::New();
>>>  WriterType::Pointer writer = WriterType::New();
>>>
>>>  reader->SetFileName( argv[1] );
>>> #if USE_MXA_FILE
>>>  reader->SetDatasetPath(argv[2]);
>>> #endif
>>>  reader->Update();
>>>  writer->SetFileName( argv[3]);
>>>  typedef itk::ResampleImageFilter<InputImageType, OutputImageType >
>>> FilterType;
>>>  FilterType::Pointer filter = FilterType::New();
>>>  typedef itk::AffineTransform< double, Dimension >  TransformType;
>>>  TransformType::Pointer transform = TransformType::New();
>>>
>>>  typedef  itk::NearestNeighborInterpolateImageFunction<InputImageType,
>>> double  >  InterpolatorType;
>>>  InterpolatorType::Pointer interpolator = InterpolatorType::New();
>>>  filter->SetInterpolator( interpolator );
>>>  filter->SetDefaultPixelValue( 0 );
>>>  const InputImageType::SpacingType& spacing = reader->GetOutput()-
>>> >GetSpacing();
>>>  const InputImageType::PointType& origin = reader->GetOutput()-
>>> >GetOrigin();
>>>  InputImageType::SizeType size = reader->GetOutput()-
>>> >GetLargestPossibleRegion().GetSize();
>>>  filter->SetOutputOrigin( origin );
>>>  filter->SetOutputSpacing( spacing );
>>>  filter->SetSize( size );
>>>
>>>  filter->SetInput( reader->GetOutput() );
>>>  writer->SetInput( filter->GetOutput() );
>>>
>>>  std::cout << "Input Size: " << size[0] << "," << size[1] <<  std::endl;
>>>  std::cout << "Input Origin: " << origin[0] << "," << origin[1] <<
>>>  std::endl;
>>>  std::cout << "Input Spacing: " << spacing[0] << "," << spacing[1]  <<
>>> std::endl;
>>>
>>>  TransformType::OutputVectorType translation1;
>>>  translation1[0] =   -origin[0];
>>>  translation1[1] =   -origin[1];
>>>  transform->Translate( translation1 );
>>>
>>>  const double degreesToRadians = atan(1.0) / 45.0;
>>>  transform->Rotate2D( -10.0 * degreesToRadians, false );
>>>
>>>  TransformType::OutputVectorType translation2;
>>>  translation2[0] =   origin[0];
>>>  translation2[1] =   origin[1];
>>>  transform->Translate( translation2, false );
>>>  filter->SetTransform( transform );
>>>  try
>>>   {
>>>   writer->Update();
>>>   }
>>>  catch( itk::ExceptionObject & excep )
>>>   {
>>>   std::cerr << "Exception caught !" << std::endl;
>>>   std::cerr << excep << std::endl;
>>>   }
>>>
>>>  return EXIT_SUCCESS;
>>> }
>>>
>>>
>>> output:
>>> GenerateData()
>>> GenerateData()
>>> Input Size: 1292,968
>>> Input Origin: 46978.9,48347.8
>>> Input Spacing: 0.207987,0.207987
>>> GenerateData()
>>>
>>> Any help would be great.
>>> _________________________________________________________
>>> Mike Jackson                  mike.jackson at bluequartz.net
>>> BlueQuartz Software                    www.bluequartz.net
>>> Principal Software Engineer                  Dayton, Ohio
>>>
>>>
>>>
>>
>> _____________________________________
>> 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
>>
>