Numpy vnl bridge unexpected behavior

(Simon Rit) #1

Hi,
I tend to write too many things on the same line in Python and this caused an unexpected behavior for the numpy vnl bridge. It occurred in a different context but I have the following minimal example

import itk
                                                                                                                                                                                                                                                                                                                                                                                                                                         
m = itk.Matrix[itk.D,4,4]()
m.SetIdentity()
m2 = m+m
print(itk.GetArrayFromVnlMatrix((m+m).GetVnlMatrix().as_matrix()))
print(itk.GetArrayFromVnlMatrix(m2.GetVnlMatrix().as_matrix()))

which should show the same two matrices but actually shows

[[4.68850067e-310 0.00000000e+000 0.00000000e+000 0.00000000e+000]
 [0.00000000e+000 2.00000000e+000 0.00000000e+000 0.00000000e+000]
 [0.00000000e+000 0.00000000e+000 2.00000000e+000 0.00000000e+000]
 [0.00000000e+000 0.00000000e+000 0.00000000e+000 2.00000000e+000]]
[[2. 0. 0. 0.]
 [0. 2. 0. 0.]
 [0. 0. 2. 0.]
 [0. 0. 0. 2.]]

Any clue what’s going on and if this can and should be corrected? It seems to be some memory release before the deep copy of the numpy array but I could not see what was wrong in the bridge code.
Thanks in advance!
Simon

2 Likes
(Matt McCormick) #2

@fbudin any ideas?

(Francois Budin) #3

No idea, but it is worth the look. Interestingly enough, since the implementation between images and VNL matrices are similar for the NumPy bridge, this problem could potentially arise with images…

(Francois Budin) #4

As a follow up, trying:

In [7]: a=(m+m).GetVnlMatrix()

In [8]: print(itk.GetArrayFromVnlMatrix(a.as_matrix()))
[[4.65803687e-310 6.90632463e-310 6.90632461e-310 6.90632465e-310]
 [6.90632463e-310 6.90632463e-310 6.90632464e-310 6.90632464e-310]
 [6.90632463e-310 6.90632464e-310 6.90632463e-310 6.90632463e-310]
 [6.90632462e-310 6.90632462e-310 6.90632461e-310 6.90632464e-310]]

also returns garbage, which makes me think that the problem is in GetVnlMatrix(), not in the bridge part.

This is most likely due to the return by reference here.

1 Like
(Francois Budin) #5

This PR solves the problem on my computer. @simon.rit: Could you test it on your machine? It is not the solution I wanted to implement (I wanted to use a SWIG typemap), but I couldn’t get the other solution to work. So I think that for now this could be a reasonable patch.

If this works for you, would you mind improving the PR with a test?

Thanks!

2 Likes
(Simon Rit) #6

Great, thanks a lot for your efforts. Yes, I’ll test and provide a test. I now understand the reference issue and can provide a simpler example for the problem, without some one line operation:

>>> m = itk.Matrix[itk.D,3,3]()
>>> m.SetIdentity()
>>> v = m.GetVnlMatrix()
>>> del m
>>> print(itk.GetArrayFromVnlMatrix(v.as_matrix()))
[[4.65884433e-310 0.00000000e+000 0.00000000e+000]
 [0.00000000e+000 1.00000000e+000 0.00000000e+000]
 [0.00000000e+000 0.00000000e+000 1.00000000e+000]]

I’ll base my test on this, just the time to recompile the wrapping with your master branch and your PR.

2 Likes
(Francois Budin) #7

I created another PR to add helper functions to convert NumPy array to ITK matrices (and back) directly, without having to worry about VNL matrices:

1 Like
(Simon Rit) #8

See my PR


I wasn’t sure if a new commit was in order. Let me know if this is not how it should be done, I’m happy to correct the PR if you tell me how.

1 Like