[SOLVED]TransformIndexToPhysicalPoint manually


Hi all,

I have a question How could I do the things that TransformIndexToPhysicalPoint does but manually?

I know that here (https://itk.org/Doxygen/html/Examples_2DataRepresentation_2Image_2Image4_8cxx-example.html) I have available an example.

But I need to transform form pixel index to physical position avoiding ITK templates and methods, I need to do that only with C basic data type. I saw the example and I am a bit confusing when I try to do the same in C.

Thanks in advance.

[SOLVED]Voxel to phisycal position in GPU
(Pablo Hernandez-Cerdan) #2

Hi @rolof, take advantage of the open source nature of ITK and study the source code of ITK to find out the details:

The transform is a simple equation, including a change of origin, scaling with the spacing, and a rotation matrix.

Without having direction into account this would be: PhysicalPoint = Origin + (Spacing * IdentityMatrix) * Index.

With Direction into account the final transform is:

PhysicalPoint = Origin + (Spacing * IdentityMatrix) * (RotationMatrix) * Index

Where everything is a vector/array with D components, except Identity and RotationMatrix that are a DxD Matrix. D being the image dimension.

In the code reference above m_IndexToPhysicalPoint is a matrix with value (Spacing * IdentityMatrix) * (RotationMatrix), the spacing (a D-vector) only influence the diagonal of the rotation matrix, think of it as a simple scaling transformation.

  /** Direction type alias support.  The Direction is a matrix of
   * direction cosines that specify the direction in physical space
   * between samples along each dimension. */
  using DirectionType = Matrix< SpacePrecisionType, VImageDimension, VImageDimension >;

Also here you have the code computing that matrix: https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/include/itkImageBase.hxx#L189-L214

template< unsigned int VImageDimension >
ImageBase< VImageDimension >
  DirectionType scale;

  for ( unsigned int i = 0; i < VImageDimension; i++ )
    if ( this->m_Spacing[i] == 0.0 )
      itkExceptionMacro("A spacing of 0 is not allowed: Spacing is " << this->m_Spacing);
    scale[i][i] = this->m_Spacing[i];

  if ( vnl_determinant( this->m_Direction.GetVnlMatrix() ) == 0.0 )
    itkExceptionMacro(<< "Bad direction, determinant is 0. Direction is " << this->m_Direction);

  this->m_IndexToPhysicalPoint = this->m_Direction * scale;
  this->m_PhysicalPointToIndex = m_IndexToPhysicalPoint.GetInverse();


Use the source! :slight_smile:

(Bradley Lowekamp) #3

Good references to the code, but it looks like a goof in this formula to me:

Here is a latex math version of the correct equation:

(Pablo Hernandez-Cerdan) #4

Thanks Brad, is s a vector there? wouldn’t the matrix shape of (D \dot s) be DIMx1?
I apply my own advice, and I should have better gone to the source, from the ITK Software guide:

I swapped the order of the RotationMatrix and the spacing, thanks for spoting the goof, and I learned a new word too!

(Bradley Lowekamp) #5

Yup, I made a mistake there too. I have corrected the post to hopefully avoid future confusion

(Pablo Hernandez-Cerdan) #6

D * diag(s) * i

The thing is that the diagonal of a Nx1 vector is not a NxN matrix but 1x1

We are missing a handy mathematical notation here to convert a vector into a diagonal matrix. See https://mathoverflow.net/questions/55820/vector-to-diagonal-matrix.

I don’t think the software guide notation is conventional either, but at least specify the 3x3 shape of the operation.

(Bradley Lowekamp) #7

Hmm… The diag function is common in Python numpy and Matlab, but I would a agree it is not rigorous. I have updated again with the values file in the matrix.

(Pablo Hernandez-Cerdan) #8

Nice, crystal clear now!. The diag is rigorous, but it doesn’t transform a vector to a diagonal matrix.


Thanks a lots @phcerdan and @blowekamp for the answers, so In C low level language programming it should be for 3D image the following Code If I didn’t make a mistake, With O(Origin), D(Direction), S(Spacing), I(index) and DxS(product between DxS):

float DxS11, DxS12, DxS13, DxS21, DxS22, DxS23, DxS31, DxS32, DxS33; 
float Ox, Oy, Oz;
float D11, D12, D13, D21, D22, D23, D31, D32, D33;
float S11, S22, S33;
float Ix, Iy, Iz;
float Px, Py, Pz;

const ImageType::PointType & origin = image->GetOrigin();

const ImageType::SpacingType & ImageSpacing = image->GetSpacing();
S11 = ImageSpacing[0];
S22 = ImageSpacing[1];
S33 = ImageSpacing[2];

const ImageType::DirectionType& direct = image->GetDirection();
D11 = direct[0][0];
D12 = direct[0][1];
D13 = direct[0][2];
D21 = direct[1][0];
D22 = direct[1][1];
D23 = direct[1][2];
D31 = direct[2][0];
D32 = direct[2][1];
D33 = direct[2][2];

DxS11= (D11*S11);
DxS12= (D12*S22);
DxS13= (D13*S33);
DxS21= (D21*S11);
DxS22= (D22*S22);
DxS23= (D23*S33);
DxS31= (D31*S11);
DxS32= (D32*S33);
DxS33= (D33*S33);

for (inputImageIterator.GoToBegin(); !inputImageIterator.IsAtEnd(); ++inputImageIterator)

  Px = Ox + ((DxS11*Ix) + (DxS12*Iy) + (DxS13*Iz));
  Py = Oy + ((DxS21*Ix) + (DxS22*Iy) + (DxS33*Iz));
  Pz = Oz + ((DxS31*Ix) + (DxS32*Iy) + (DxS33*Iz));

It could be improved with vectors instead of variables, but I do not have much time, and I wanted to test if it works and thanks to let me know how to get it working.

(Pablo Hernandez-Cerdan) #10

Pasting here a better formatted version of the transformation. I will contribute it to the ITK Software Guide to replace the non-standard math notation there. Hope is helpful!