### Author:  [@rolof](https://discourse.itk.org/u/rolof)
#### Post date: [June 25, 2018, 5:28pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/1 "2018-06-25T17:28:38Z")
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](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.
---
### Author:  [@phcerdan](https://discourse.itk.org/u/phcerdan)
#### Post date: [June 26, 2018, 1:58pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/2 "2018-06-26T13:58:22Z")
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:
```auto
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.
```auto
/** 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](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/include/itkImageBase.hxx#L189-L214)
```auto
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! 🙂
---
### Author:  [@blowekamp](https://discourse.itk.org/u/blowekamp)
#### Post date: [June 27, 2018, 2:35pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/3 "2018-06-27T14:35:40Z")
Good references to the code, but it looks like a goof in this formula to me:
> [@phcerdan](#):
>
> PhysicalPoint = Origin + (Spacing \* IdentityMatrix) \* (RotationMatrix) \* Index
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}\")
---
### Author:  [@phcerdan](https://discourse.itk.org/u/phcerdan)
#### Post date: [June 27, 2018, 3:14pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/4 "2018-06-27T15:14:22Z")
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!
---
### Author:  [@blowekamp](https://discourse.itk.org/u/blowekamp)
#### Post date: [June 27, 2018, 3:24pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/5 "2018-06-27T15:24:00Z")
Yup, I made a mistake there too. I have corrected the post to hopefully avoid future confusion
---
### Author:  [@phcerdan](https://discourse.itk.org/u/phcerdan)
#### Post date: [June 27, 2018, 3:33pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/6 "2018-06-27T15:33:13Z")
> 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](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.
---
### Author:  [@blowekamp](https://discourse.itk.org/u/blowekamp)
#### Post date: [June 27, 2018, 3:47pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/7 "2018-06-27T15:47:27Z")
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.
---
### Author:  [@phcerdan](https://discourse.itk.org/u/phcerdan)
#### Post date: [June 27, 2018, 3:53pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/8 "2018-06-27T15:53:29Z")
Nice, crystal clear now!. The `diag` is rigorous, but it doesn’t transform a vector to a diagonal matrix.
---
### Author:  [@rolof](https://discourse.itk.org/u/rolof)
#### Post date: [June 29, 2018, 10:10am UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/9 "2018-06-29T10:10:26Z")
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.
---
### Author:  [@phcerdan](https://discourse.itk.org/u/phcerdan)
#### Post date: [January 9, 2019, 5:59pm UTC](https://discourse.itk.org/t/solved-transformindextophysicalpoint-manually/1031/10 "2019-01-09T17:59:03Z")
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!
