Set Image Direction from numpy array

Hello,

I’m struggling to use a numpy array in order to set the direction of an ITK Image using python bindings.

My attempts:

In [8]: image = itk.Image.New()

In [9]: np_dir = np.eye(2)

In [10]: image.SetDirection(itk.GetVnlMatrixFromArray(np_dir))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-7b1c85a3c519> in <module>()
----> 1 image.SetDirection(itk.GetVnlMatrixFromArray(np_dir))

TypeError: in method 'itkImageBase2_SetDirection', argument 2 of type 'itkMatrixD22 const &'

okay, so you want an itkMatrixD22, not a VnlMatrix.

Let’s try that:

In [11]: itk_dir = itk.Matrix[itk.D, 2, 2]

In [12]: dir(itk_dir)
Out[12]: 
[[...]
 'GetVnlMatrix',
 [...]]

No SetVnlMatrix (I would have liked that) but maybe I can GetVnlMatrix and then use this object to set the value…

In [13]: itk_dir.GetVnlMatrix()
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-13-83a493e54f83> in <module>()
----> 1 itk_dir.GetVnlMatrix()

NotImplementedError: Wrong number or type of arguments for overloaded function 'itkMatrixD22_GetVnlMatrix'.
  Possible C/C++ prototypes are:
    itkMatrixD22::GetVnlMatrix()
    itkMatrixD22::GetVnlMatrix() const

How do I set the values of my itkMatrixD22? I couldn’t find. :frowning:

I may be in a wrong direction, but as my post’s title states, I want to use a python object (a numpy array, but a tuple or a list are fine too) to set the direction matrix of my image. How do I do that?

Thanks

1 Like

Hello @nicoco,

It is a little tricky. Here are the calls:

np_dir_vnl = itk.GetVnlMatrixFromArray(np_dir)
direction = image.GetDirection()
direction.GetVnlMatrix().copy_in(np_dir_vnl.data_block())

HTH,
Matt

3 Likes

Thank you!

This is exactly what I needed and it is indeed a little tricky.

1 Like

Hi,

Having some trouble getting the VnlMatrix from a numpy array

I perform a PCA on the ‘positive’ pixels coordinates to approximate the ‘main axes’ ahead of a rotation. However, I got errors casting the Numpy.array to VnlMatrix:

In [56]: eigenVectors
Out[56]:
array([[ 0.        ,  1.        , -0.        ],
       [-0.70710678,  0.        , -0.70710678],
       [ 0.70710678,  0.        , -0.70710678]])

In [57]: rotationMatrix = itk.GetVnlMatrixFromArray(eigenVectors)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
RuntimeError: Cannot get an instance of NumPy array.

During handling of the above exception, another exception occurred:

SystemError                               Traceback (most recent call last)
<ipython-input-57-40745ee4ad3d> in <module>()
----> 1 rotationMatrix = itk.GetVnlMatrixFromArray(eigenVectors)

~/anaconda3/lib/python3.5/site-packages/itkExtras.py in GetVnlMatrixFromArray(arr)
    336     """Get a vnl matrix from a Python array.
    337     """
--> 338     return _GetVnlObjectFromArray(arr, "GetVnlMatrixFromArray")
    339
    340 # return an image

~/anaconda3/lib/python3.5/site-packages/itkExtras.py in _GetVnlObjectFromArray(arr, function)
    326     PixelType = _get_itk_pixelid(arr)
    327     templatedFunction = getattr(itk.PyVnl[PixelType], function)
--> 328     return templatedFunction(arr)
    329
    330 def GetVnlVectorFromArray(arr):

~/anaconda3/lib/python3.5/site-packages/itk/Configuration/../itkPyBufferPython.py in GetVnlMatrixFromArray(ndarr)
   6506             "Only arrays of 2 dimensions are supported."
   6507
-> 6508         mat = itkPyVnlD._GetVnlMatrixFromArray( ndarr, ndarr.shape)
   6509
   6510         return mat

~/anaconda3/lib/python3.5/site-packages/itk/Configuration/../itkPyBufferPython.py in _GetVnlMatrixFromArray(arr, shape)
   6391     def _GetVnlMatrixFromArray(arr: 'PyObject *', shape: 'PyObject *') -> "vnl_matrixD const":
   6392         """_GetVnlMatrixFromArray(PyObject * arr, PyObject * shape) -> vnl_matrixD"""
-> 6393         return _itkPyBufferPython.itkPyVnlD__GetVnlMatrixFromArray(arr, shape)
   6394
   6395     _GetVnlMatrixFromArray = staticmethod(_GetVnlMatrixFromArray)

SystemError: <built-in function itkPyVnlD__GetVnlMatrixFromArray> returned a result with an error set

@T4mmi could you please save eigenVectors and share them?

Here is the matrix

eigenVectors.npy (200 Bytes)

I obtain it from the following :

# Main axes computation
objectImageFilter.Update()
mainConnComp = itk.GetArrayViewFromImage(objectImageFilter.GetOutput())
pointsCoords = numpy.asarray(numpy.where(mainConnComp)).T
# rotationCenter = [round(i/2) for i in mainConnComp.shape]  # rotate around the image center
rotationCenter = [int(i) for i in numpy.mean(
    pointsCoords, axis=0)]  # rotate around the object center
covarianceMatrix = numpy.cov(pointsCoords, rowvar=False)
pointsCoords = None
eigenValues, eigenVectors = numpy.linalg.eigh(covarianceMatrix)
covarianceMatrix = None
descendantIndices = numpy.argsort(eigenValues)[::-1]
eigenVectors = eigenVectors[:, descendantIndices] 

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))
1 Like

Well, again thanks a lot

1 Like

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

3 Likes

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

).

2 Likes

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:

1 Like

Thanks Matt, it’s working now.
:+1:

1 Like