[Insight-users] Attempt to rotate image resulting in black image.
Darren Weber
darren.weber.lists at gmail.com
Wed Nov 4 16:54:41 EST 2009
Hi Mike,
I recently observed similar results using some microscope images (in tiff
format). The solution for me was to have a program input to specify the
"real" pixel spacing. Those input values were given to an
itkChangeInformationImageFilter and the output of that filter was fed into
the resampler (along with the transform and the spacing etc. from the
ChangeInformationFilter), so that everything was in "real world
coordinates". The output of the resampler was sent to the image writer.
http://www.itk.org/Doxygen316/html/classitk_1_1ChangeInformationImageFilter.html
Let me know if you want an example src file.
Take care,
Darren
On Thu, Jun 4, 2009 at 6: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
>>>
>>
> _____________________________________
> 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/20091104/c3965880/attachment-0001.htm>
More information about the Insight-users
mailing list