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

Luis Ibanez luis.ibanez at kitware.com
Sun Nov 8 15:03:26 EST 2009


Hi Michael,

Please do the following:


Just after you call   writer->Update() do:


    reader->GetOutput()->Print( std::cout );
    std::cout << "Transform " << transform << std::endl;
    filter->GetOutput()->Print( std::cout );


and post the output to the mailing list.

This will help us track down any potential
issue with the coordinate system of your
images.


     Thanks


          Luis


---------------------------------------------------------
On Thu, Jun 4, 2009 at 9:22 AM, Michael Jackson
<mike.jackson at bluequartz.net> wrote:
> I took the code directly from the RescaleExample5.cxx and just substituted
> my reader class in for the "itkImageFileReader" class and I still get a
> black image. The code is pasted below. The general characteristics of my
> image are:
>
> Input Size: 1292,968 (pixels)
> Input Origin: 46978.9,48347.8  (Microns)
> Input Spacing: 0.207987,0.207987 (Microns/pixel)
> Image Center:  47113.3,48448.5 (Microns)
>
> Is it possible my itkR3DImageIO class is not setting something correctly?
>
>  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;
>
>  //typedef itk::ImageFileReader< InputImageType  >  ReaderType;
>  typedef itkR3DImageIO ReaderType;
>  typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>
>  ReaderType::Pointer reader = ReaderType::New();
>  WriterType::Pointer writer = WriterType::New();
>
>  reader->SetFileName( argv[1] );
>  reader->SetDatasetPath( argv[2]);
>  writer->SetFileName( argv[3] );
>
>  const double angleInDegrees = 45.0;// atof( argv[3] );
>  const double scale          = 0.5; //atof( argv[4] );
>
>  typedef itk::ResampleImageFilter<
>                  InputImageType, OutputImageType >  FilterType;
>
>  FilterType::Pointer filter = FilterType::New();
>
>  typedef itk::Similarity2DTransform< double >  TransformType;
>  TransformType::Pointer transform = TransformType::New();
>  typedef itk::LinearInterpolateImageFunction<
>                       InputImageType, double >  InterpolatorType;
>  InterpolatorType::Pointer interpolator = InterpolatorType::New();
>
>  filter->SetInterpolator( interpolator );
>
>  filter->SetDefaultPixelValue( 0 );
>
>  reader->Update();
>  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() );
>
>  TransformType::InputPointType rotationCenter;
>  rotationCenter[0] = origin[0] + spacing[0] * size[0] / 2.0;
>  rotationCenter[1] = origin[1] + spacing[1] * size[1] / 2.0;
>  transform->SetCenter( rotationCenter );
>  const double degreesToRadians = atan(1.0) / 45.0;
>  const double angle = angleInDegrees * degreesToRadians;
>  transform->SetAngle( angle );
>
>  transform->SetScale( scale );
>  TransformType::OutputVectorType translation;
>
>  translation[0] =   13.0;
>  translation[1] =   17.0;
>
>  transform->SetTranslation( translation );
>
>  filter->SetTransform( transform );
>
>  try
>    {
>    writer->Update();
>    }
>  catch( itk::ExceptionObject & excep )
>    {
>    std::cerr << "Exception catched !" << std::endl;
>    std::cerr << excep << std::endl;
>    }
>  return EXIT_SUCCESS;
> _________________________________________________________
> Mike Jackson                  mike.jackson at bluequartz.net
> BlueQuartz Software                    www.bluequartz.net
> Principal Software Engineer                  Dayton, Ohio
>
>
>
> On Jun 3, 2009, at 8:42 PM, Luis Ibanez 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
>
>


More information about the Insight-users mailing list