[Insight-users] Resampler filter giving blank output

[Luis Ibanez]

On Sat, Oct 10, 2009 at 4:09 PM, darkflame2000 at libero.it
<darkflame2000 at libero.it> wrote:
>> Hi Andrea,
>>
>> Well, ...
>> a multiple of 360 degrees is equivalent to zero degrees,
>> and therefore, to not rotating at all.
>>
>> Something fundamental is wrong here...
>>
>> Could you please do the following:
>>
>> a) Upload your input image to:
>> http://www.kitware.com/cgi-bin/uploadfile.cgi
>> and indicate that it is intended for me.
>>
>> b) Send me you full program
>>     (assuming that it is a single .cxx file)
>>
>> In that way I will be able to run in it here
>> an attempt to replicate what you are seeing.
>>
>>
>>      Thanks
>>
>>
>>            Luis
>>
>>
> Hi Luis,
>
> thank you very much for the help! I've sent you a .ZIP file where you can find the entire project that includes:
>
> - the atlas 3d image (aa.hdr and aal.img)
>
> - a .TXT file that contains all the anatomical part extracted from Tzourio-Mazoyer's Brain Atlas
>
> - a CMake file
>
> The goal of my project is to extract one of the 116 parts (or more than one part) of the Tzourio-Mazoyer's AAL atlas from it, and then save such a part on an arbitrary sized 3D image. This 3d image will be visualized with a program ITKSNAP...
> The part of the project that extracts from the .TXT file a number (that concerns a determinate color) works correctly. I've sent a version of project that extract all the 116 parts ( with the instruction inside aalextractor.cpp: ThresholdFilter->SetLowerThreshold(1); ThresholdFilter->SetUpperThreshold(116);  )
>
> The class aalImageProcessor is the problem!
> It implements the rotation, the resize, and the conversion from itk to vtk image.
> I've tried to resize the atlas  3d image in two differents formats:
>
>        output_dim[0] = 181;//78;
>        output_dim[1] = 217;//110;
>        output_dim[2] = 181;//64;
>
> 181 x 217 x 181 are the same dimensions of the atlas  3d image (scaling=1 for every directions)
> 78 x 110 x 64 are the dimensions that I want for the output 3d image.
>
> Inside the main you can see these instructions:
>
>    imageprocessor -> Rotation (itkPart, itkRotate);
>    imageprocessor -> Process (itkRotate, itkResize);
>    imageprocessor -> Paste (itkResize, itkFinal);
>    imageprocessor -> Convert (itkFinal, vtkPart);
>
> that represent the image processing that I want.
>
> If I try to rotate 360 degrees without resizing the 3d image (using the dimension 181 x 217 x 181) the program produces a correct output, that I visualize with ITKSNAP...
>
> If I try to rotate 180 degrees without resizing the 3d image (using the dimension 181 x 217 x 181) the program produces an incorrect output: the output file is blank.
>
> If I try to rotate 360 degrees resizing the 3d image (using the dimension 78 x 110 x 64) the program produces an incorrect output, the output file is blank.
>
> If I try to rotate 180 degrees resizing the 3d image (using the dimension 78 x 110 x 64) the program produces an incorrect output: the output file is blank.
>
> The last case is the most important for me... It's the resolution of the problem.
>
> At this project partecipates also Arturo Caissut. We are two students that work together at this project.
>
> I hope that your experience can help us...
>
> Thank you very much,
>
> Andrea
>
>



-----


Hi Andrea

Thanks for sending your code.
It was very helpful.

Here is what happens:

You input image has a Direction Cosines matrix
that is not an Identity.

Instead it is:

  Direction:
              1 0 0
              0 -1 0
              0 0 1

However, your computation of the image
center in the function that performs the
rotation was done as:

	double imageCenter[3];
	imageCenter[0] = origin[0] + spacing[0] * size[0] / 2.0;
	imageCenter[1] = origin[1] + spacing[1] * size[1] / 2.0;
	imageCenter[2] = origin[2] + spacing[2] * size[2] / 2.0;
	
This way of computing the center physical coordinates
is only valid with the Direction Cosines matrix is an identity.

In general the computation should be (in vectors):

          center = origin + DM * SP * Index

where

DM = Direction cosines
SP = Spacing (diagonal) matrix

The simple fix, is to replace the computation of the
center above with:

 AtlasImageType::IndexType centralIndex;

  centralIndex[0] = size[0] / 2.0;
  centralIndex[1] = size[1] / 2.0;
  centralIndex[2] = size[2] / 2.0;

  AtlasImageType::PointType centralPoint;

  InputImage->TransformIndexToPhysicalPoint( centralIndex, centralPoint );


The image method TransformIndexToPhysicalPoint() will take
into account the image: Origin, Spacing, and Direction Cosines.

The problem that you were observing was due to the fact
that the previous computation of the image center was actually
setting the center physical coordinates outside of the physical
extent of the image.

Therefore, as you rotated the image, it was moved away from
the extent of the grid where you could see any of its pixels.

---

Please note also that you can simplify / improve the API
of your functions by passing the second smart pointer by
reference. In this way you don't have to Graft() the output
image in it.

Please find attached the modified files

   * aalImageProcessor.h
   * aalImageProcessor.cpp


Please give them a try and let us know if you find
any residual problem.


     Thanks


            Luis
-------------- next part --------------
A non-text attachment was scrubbed...
Name: aalImageProcessor.cpp
Type: text/x-c++src
Size: 6260 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091012/e69dde57/attachment.cpp>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: aalImageProcessor.h
Type: text/x-chdr
Size: 1847 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091012/e69dde57/attachment.h>