Thanks @T4mmi.
It looks like, at least for now, itk.GetVnlMatrixFromArray requires a C-order array and does not handle a Fortran-order array. We will work on fixing that, or at least improving the error message.
To convert it, try:
rotationMatrix = itk.GetVnlMatrixFromArray(np.ascontiguousarray(eigenVectors))
Well, again thanks a lot …
Hi,
I tried to reproduce this code but I get
AttributeError: ‘SwigPyObject’ object has no attribute ‘copy_in’
with the master HEAD. Has something changed? Is there any other way to modify the values of an itk.Matrix than to use a numpy array?
Thanks,
Simon
There is currently no easy way to change the values in an ITK matrix. You can query the value with [ ], but to change the value you have to use a VNL matrix. That could be improved…
I haven’t tried yet to reproduce your error with master, but in theory nothing should have changed.
Thanks. Actually, I had wrapped new itk::Matrix templates but I had not wrapped the corresponding vnl_matrix_fixed. Now that this is fixed, I can easily change the values using (here with a matrix wrapped by default),
>>> import itk
>>> m = itk.Matrix[itk.D,3,3]()
>>> print(m(2,1))
0.0
>>> m.GetVnlMatrix().put(2,1,3.14)
>>> print(m(2,1))
3.14
Be aware that my previous solution is no longer valid after this commit (after version 5rc01):
If you want to modify an ITK matrix from python, the simplest in my opinion is to go to NumPy with itk.GetArrayFromVnlMatrix and itk.GetVnlMatrixFromArray (or soon with GetArrayFromMatrix and GetMatrixFromArray see
).
Hi simon.rit,
I found m.GetVnlMatrix().set(2, 1, 3.14) can assign value to the matrix as well. What is difference between “set” and “put”?
Regards,
Zhuangming Shen
They both do the same job but set returns a reference to the matrix. See the code here. I’m not sure when set would be useful to be honest. Both put and set won’t change the elements of the itk.Matrix after this commit.
Thanks for your explanation.
Hi,
I’m having trouble updating the direction on an existing 3d itkImage.
Basically I want to reset the origin to [0,0,0] and direction to eye(3).
This first step seems to work but then when I use TransformPhysicalPointToContinuousIndex it returns the wrong index position. Here is a minimal example reproducing the error.
from __future__ import print_function
import itk
import sys
import numpy as np
def saveDirection(image):
''' copy direction to numpy array '''
direction = np.zeros((3,3))
for i in range(9):
direction[i//3, i%3] = image.GetDirection().GetVnlMatrix().get(i//3, i%3)
return direction
def readImageStack(folder):
''' load a single dicom series from a folder and reset origin and direction '''
ImageType = itk.Image[itk.SS, 3]
namesGenerator = itk.GDCMSeriesFileNames.New()
namesGenerator.SetUseSeriesDetails(True)
namesGenerator.SetGlobalWarningDisplay(False)
namesGenerator.SetDirectory(folder)
seriesUID = namesGenerator.GetSeriesUIDs()
if len(seriesUID) != 1:
print("Could not load images in folder %s. No image or multiple series."%folder)
sys.exit(1)
uid = seriesUID[0]
fileNames = namesGenerator.GetFileNames(uid)
reader = itk.ImageSeriesReader[ImageType].New()
dicomIO = itk.GDCMImageIO.New()
reader.SetImageIO(dicomIO)
reader.SetFileNames(fileNames)
reader.ForceOrthogonalDirectionOff()
reader.Update()
image = reader.GetOutput()
oldDirection = saveDirection(image)
print("old direction :\n", oldDirection)
image.SetOrigin((0,0,0)) # reset origin : works
# reset direction
np_dir_vnl = itk.GetVnlMatrixFromArray(np.eye(3))
direction = image.GetDirection()
direction.GetVnlMatrix().copy_in(np_dir_vnl.data_block())
image.SetDirection(direction)
# old code, gives the same bad result
#tmp = image.GetDirection().GetVnlMatrix()
#for i in range(9):
# tmp.set(i//3, i%3, 0)
#for i in range(3):
# tmp.set(i, i, 1)
#direction = image.GetDirection()
#direction.SetVnlMatrix(tmp)
#image.SetDirection(direction)
return image, oldDirection
if __name__ == "__main__":
DATA = "F:/Data/Images/Pat_John Doe/Study_1.3.12.2.1107.5.4.5.164089.30000018061106463446200000037/Series_1.3.12.2.1107.5.4.5.164089.30000018061106500168800005613"
image, oldDirection = readImageStack(DATA)
print(image)
# prob intensity at a given 3d physical point
p = np.array([11.44909474, 23.79051299, 40.75303588])
# since Direction is now identity and origin is [0,0,0], continuous index should be p/spacing
spacing = image.GetSpacing()
print("expected continuous index", [c/s for c,s in zip(p, spacing)])
# using TransformPhysicalPointToContinuousIndex -> off by a few voxels
pixelIndex = image.TransformPhysicalPointToContinuousIndex(p)
print(pixelIndex)
#using saved old direction :
#S-1*D-1*p
print("using old direction", np.linalg.inv(np.diag(spacing)).dot(oldDirection.T).dot(p))
The code above gives the following output :
old direction :
[[ 0.99998595 0.00174611 -0.00500522]
[-0.00175405 0.99999721 -0.00158081]
[ 0.00500245 0.00158957 0.99998622]]
Image (0000020CF83CC080)
RTTI typeinfo: class itk::Image<short,3>
Reference Count: 1
Modified Time: 326554
Debug: Off
Object Name:
Observers:
none
Source: (none)
Source output name: (none)
Release Data: Off
Data Released: False
Global Release Data: Off
PipelineMTime: 310064
UpdateMTime: 326552
RealTimeStamp: 0 seconds
LargestPossibleRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
BufferedRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
RequestedRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
Spacing: [0.108077, 0.108077, 0.108077]
Origin: [0, 0, 0]
Direction:
1 0 0
0 1 0
0 0 1
IndexToPointMatrix:
0.108076 0.000188715 -0.000540951
-0.000189573 0.108077 -0.00017085
0.000540651 0.000171796 0.108076
PointToIndexMatrix:
9.25251 -0.0162296 0.0462858
0.0161562 9.25261 0.0147077
-0.0463115 -0.0146267 9.25251
Inverse Direction:
0.999986 -0.00175405 0.00500245
0.00174611 0.999997 0.00158957
-0.00500522 -0.00158081 0.999986
PixelContainer:
ImportImageContainer (0000020D24629F90)
RTTI typeinfo: class itk::ImportImageContainer<unsigned __int64,short>
Reference Count: 1
Modified Time: 310132
Debug: Off
Object Name:
Observers:
none
Pointer: 0000020D09DF4040
Container manages memory: true
Size: 122945536
Capacity: 122945536
expected continuous index [105.93429698107786, 220.1249378791367, 377.0729742247302]
itkContinuousIndexD3 ([107.433, 220.909, 376.19])
using old direction [107.43298801 220.9086811 376.18957873]
We can see that origin and direction of the image are properly set to the expected values. But somehow IndexToPointMatrix, PointToIndexMatrix and Inverse Direction are not.
This then leads to incorrect behavior of image.TransformPhysicalPointToContinuousIndex() as it uses this out-of-date matrices.
Is this a bug or am I doing something wrong when updating the direction ? I did not try this in C++ but would the result be the same ?
I’m using itk for python version 5.0b3.
I remember there were some related fixed to Python wrapping recently. Could you try ITK for Python version 5.0.0.post1, or soon to be released 5.0.1?
Just upgraded to 5.0.0post1 and I get the same exact output.
While this patch, included in 5.0.0,:
fixed one issue, it prevented the above approach from working.
A different approach, which creates a new Matrix, instead of replacing the buffer contents of an existing on, is to use
np_dir_vnl = itk.vnl_matrix_from_array(np.eye(3))
DirectionType = type(image.GetDirection())
direction = DirectionType(np_dir_vnl)
image.SetDirection(direction)
That said, this is all not very Pythonic, and we can improve the API by supporting:
image.SetDirection(np.eye(3))
i.e., passing a NumPy array directly.
And support the more Pythonic:
image.direction = np.eye(3)
Tracking issues:
Thanks Matt, it’s working now.
It is worth noting that ITK 5.0 already includes a much cleaner approach here that avoids interacting with vnl_matrix:
new_direction = np.eye(2)
image.SetDirection(itk.matrix_from_array(new_direction))
@matt.mccormick — this is very nice indeed. I’m wondering, where can I have a look at the rest of Python API (list of what’s what) for my general education on the matter of Python ITK for functionality discovery?
I’m aware there’s always python online help such as dir(itk) and itk.*matrix*? and friends, along with https://docs.itk.org/en/latest/learn/python_quick_start.html, but I can’t seem to be able an equivalent to https://itk.org/Doxygen/html/index.html, which facilitates functionality discovery very much.
@constantine yes, the quick start is helpful, and the embedded docs, accessed with the Python built-ins, dir() or help() are informative, albeit verbose. One motivation for ITK-Wasm is that it enables Pythonic interfaces along with Markdown/Sphinx documentation, linkable via InterSphinx – example.
Yes! I see ITK-Wasm (Python) API reference — it’s great.
Is there ITK Python API reference like that?
Not yet – it would come with ITK 6.
ITK 6 — this is exciting!