# Obtain new point after reslice a NIFTII image

I’m a new user of VTK. And excuse my english.
I have a niftii image volume with many landmarks. I want to reslice this image with an affine transformation matrix, at this point i use this program and that’s work well.

``````import vtk
import numpy as np
import random
import pandas as pd
import scipy.stats
import numpy.linalg as npl
from datetime import datetime
from functions import *

random.seed(datetime.now())

def samples_pdf(df):
samples = []
return samples

def generate_mat_9var(s):
getsamp = lambda samples, i: samples[i].resample(1).T[:,0][0]
MTrans = np.array([
[getsamp(s,0), getsamp(s,1), getsamp(s,2), 0.],
[getsamp(s,3), getsamp(s,4), getsamp(s,5), 0.],
[getsamp(s,6), getsamp(s,7), getsamp(s,8), 0.],
[0.,  0.,  0., 1.]
])
return MTrans

""" Return mm point (world) to voxel point """
def mmtovox(pointmm, target_affine):
return (nib.affines.apply_affine(npl.inv(target_affine), pointmm))

""" Return voxel point to mm point (world) """
def voxtomm(pointvox, target_affine):
return nib.affines.apply_affine(target_affine, pointvox)

samples = samples_pdf(df_9var)

print("Open image")
file = "../../DATA/100RAS.nii.gz"

print("Make transform")
bounds = np.array(Output.GetBounds())
transform = vtk.vtkTransform()
center = np.zeros(3)
transformCenter = np.zeros(3)
for i in range( len(center) ):
center[i] = 0.5 * (bounds[i*2]+bounds[i*2+1])

Mat = generate_mat_9var(samples)
transform.SetMatrix(Mat.flatten())

print("Reslice Image")
reslice = vtk.vtkImageReslice()
reslice.SetInputData( Output )
reslice.AutoCropOutputOn();
reslice.SetResliceTransform( transform )
reslice.SetBackgroundLevel( Output.GetPointData().GetScalars().GetRange()[0] )
reslice.Update()

print("Save Image")
savevol(reslice)
``````

But when I want to obtain the new point coordinates, that’s not work.
I’ve tried some calcul like
`MyNewPoint = np.Dot( MyActualPoint, Mat) + origin`

Anyone can help me pleaz ?

Could you please put your code between three bactick characters? It is now formatted as markdown and really hard to read. [quote=“LoiseauNicolas, post:1, topic:4054”]
when I want to obtain the new point coordinates, that’s not work.
I’ve tried some calcul like
`MyNewPoint = np.Dot( MyActualPoint, Mat) + origin`
[/quote]

Normally translation is included in the homogeneous transformation matrix and you multiply point vector from left:

``MyPoint_LPS = np.Dot(Mat_IJKToLPS, MyPoint_IJK)``
1 Like

I’ve tried this :

MyPoint_LPS = np.Dot(Mat_IJKToLPS, MyPoint_IJK)

As you can see in the first image, I have a point in [253, 390, 254] in voxel space and [-2.5535,-40.5299,-319.125] in world units.
I have applied this matrix to the volume

array([[ 1.19189022, -0.08946897, 0.10081192, 0. ],
[ 0.04594865, 0.8423345 , 0.02081823, 0. ],
[-0.18746262, -0.35562263, 1.1700583 , 0. ],
[ 0. , 0. , 0. , 1. ]])

and when i apply

MyPoint_LPS = np.Dot(Mat_IJKToLPS, MyPoint_IJK)
I obtain

array([ -31.58892633, -40.90068131, -358.50281916])

in world units and

array([292.26155525, 345.42329302, 111.07393835])

in voxel units.

But the correct point is in [136.5, 312.7, 444] for world units or [253, 475, 297] for voxel units, as we can see in the next image.

Probably due to this line ?

reslice.AutoCropOutputOn();

I see that you work in medical imaging domain - then I would recommend to simply use VTK in 3D Slicer’s Python environment. It simplifies everything - from image reslicing to coordinate system conversions.

For example, you can easily convert between anatomical<->voxel coordinates of even non-linearly transformed (warped) volumes as shown in these Python scripts: anatomical to voxel, voxel to anatomical.

Here is an example of how to reslice a volume using a plane defined by 3 points (that you can interactively place and move): https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Set_slice_position_and_orientation_from_3_markup_fiducials

You can give 3D Slicer a quick try online via Jupyter notebooks (see more details here).

1 Like

Thank you again, I will read those articles and then I come back !

Hello Lassoan,
I have take a look on slicer3D and that seems really interesting, I think i will use it insteed of itksnap !

But I have some problem with slicer 3D … I can’t just use with python without GUI/interface. I want to affine transform some images automatically and find new positions of landmarks in a program.

I’ve tried to transform a niftii volume with affine transform and i’ve used AutoCropOutputOn() to crop all new reslice image without loose anything. That work really good but i really have a problem to get new landmarks positions.

Anyone can help me ?

You can choose to use 3D Slicer with or without GUI. You can even use it as a Jupyter notebook kernel (you can try it in your web browser, without installing anything - see here: https://github.com/Slicer/SlicerJupyter#option-1-run-using-binder).

You should be able to apply the same transform to your landmark points as you applied to the image (simply multiply the point’s homogeneous coordinates by the transformation matrix).