Obtain new point after reslice a NIFTII image

Support
1

Hello everyone and thank you for your answers,

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 = []
	for dfhead in df.columns:
		samples.append(scipy.stats.gaussian_kde(df_9var[str(dfhead)].to_numpy()))
	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)


df_9var = pd.read_csv('9var.csv',   sep=',', error_bad_lines=False)
samples = samples_pdf(df_9var)


print("Open image")
file = "../../DATA/100RAS.nii.gz"
reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(file)
reader.Update()
Output = reader.GetOutput()

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 ?

2

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
3

Hello @lassoan and thank you for your answer !

I’ve tried this :

MyPoint_LPS = np.Dot(Mat_IJKToLPS, MyPoint_IJK)

original
original899×623 148 KB

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.

newApproximatly
newApproximatly898×661 179 KB

4

Probably due to this line ?

reslice.AutoCropOutputOn();

5

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
6

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

7

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 ?

8

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: GitHub - Slicer/SlicerJupyter: Extension for 3D Slicer that allows the application to be used from Jupyter notebook).

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).

1 Like