[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