[Insight-users] Attempt to rotate image resulting in black image.
Michael Jackson
mike.jackson at bluequartz.net
Thu Jun 4 10:22:57 EDT 2009
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