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.
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 >
void
ImageBase< VImageDimension >
::ComputeIndexToPhysicalPointMatrices()
{
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();
this->Modified();
}
Use the source! 
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:![](https://discourse.itk.org/uploads/default/original/2X/f/f3fba852eae35532fecae6b68a4d7ab50ced7970.gif "\ D = \text{direction cosine matrix (NxN)} \ \mathbf{o} = \text{origin (Nx1)}\ \mathbf{s} = \text{spacing(Nx1)}\ \mathbf{i} = \text{index (Nx1)}\ \ IndexToPhysicalPoint(\mathbf{i})=\mathbf{o} +D \cdot \begin{pmatrix} s_{1} & \ & s_{2} & \ & & \ddots & \ & & & s_{n} \end{pmatrix} \cdot \mathbf{i}\")
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!
Yup, I made a mistake there too. I have corrected the post to hopefully avoid future confusion
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 linear algebra - vector to diagonal matrix - MathOverflow.
I don’t think the software guide notation is conventional either, but at least specify the 3x3 shape of the operation.
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.
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();
Ox=origin[0];
Oy=origin[1];
Oz=origin[2];
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)
{
Ix=(inputImageIterator.GetIndex())[0];
Iy=(inputImageIterator.GetIndex())[1];
Iz=(inputImageIterator.GetIndex())[2];
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.
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!
