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.
@rolof, take advantage of the open source nature of ITK and study the source code of ITK to find out the details:
/** Get a physical point (in the space which
* the origin and spacing information comes from)
* from a discrete index (in the index space)
* \sa Transform */
template< typename TCoordRep >
const IndexType & index,
Point< TCoordRep, VImageDimension > & point) const
for ( unsigned int i = 0; i < VImageDimension; ++i )
point[i] = this->m_Origin[i];
for ( unsigned int j = 0; j < VImageDimension; ++j )
point[i] += m_IndexToPhysicalPoint[i][j] * index[j];
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:
template< unsigned int VImageDimension >
ImageBase< VImageDimension >
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!
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:
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
I don’t think the software guide notation is conventional either, but at least specify the 3x3 shape of the operation.
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();
const ImageType::SpacingType & ImageSpacing = image->GetSpacing();
S11 = ImageSpacing;
S22 = ImageSpacing;
S33 = ImageSpacing;
const ImageType::DirectionType& direct = image->GetDirection();
D11 = direct;
D12 = direct;
D13 = direct;
D21 = direct;
D22 = direct;
D23 = direct;
D31 = direct;
D32 = direct;
D33 = direct;
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.
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!