UniformTriRefiner gives uncorrect results for transposed arrays

[miaocb]

The following code repeats the problem. Array ele1 is same as ele2 except that it is transpose of ele0, but the refined triangles using ele1 is incorrect.

!/bin/env python

import numpy as np
import matplotlib.tri as tri
import matplotlib.pyplot as plt

x = np.array([ 120.39299774, 120.59100342, 120.42900085, 120.31700134])
y = np.array([ 33.99900055, 34.00899887, 34.18799973, 34.18399811])

ele0 = np.array([[2, 2], [0, 1], [3, 0]])
ele1 = ele0.transpose() # ele1 is same as ele2 except that it is transpose of ele0
ele2 = np.array([[2, 0, 3], [2, 1, 0]])

triang1 = tri.Triangulation(x, y, ele1)
triang2 = tri.Triangulation(x, y, ele2)

refiner1 = tri.UniformTriRefiner(triang1)
refiner2 = tri.UniformTriRefiner(triang2)

fine_triang1 = refiner1.refine_triangulation(subdiv=1)
fine_triang2 = refiner2.refine_triangulation(subdiv=1)

fig = plt.figure()
ha1 = fig.add_subplot(121)
ha1.set_aspect('equal')
plt.triplot(fine_triang1, color='b', linewidth=0.5)
plt.triplot(triang1, color='k', linewidth=1)
plt.title('refine_triang1 is incorrect')

ha2 = fig.add_subplot(122)
ha2.set_aspect('equal')
plt.triplot(fine_triang2, color='b', linewidth=0.5)
plt.triplot(triang2, color='k', linewidth=1)
plt.title('refine_triang2 is correct')

plt.show()

[ianthomas23]

It looks fine for me using matplotlib master, python 2.7.6 and numpy 1.8.2:

What does your output image look like, and what versions of matplotlib, python and numpy are you using? It would also help if you could let us know what the output is if you add the following to the end of your script:

print('ele1', ele1)
print('ele2', ele2)
print('fine1', fine_triang1.triangles)
print('fine2', fine_triang2.triangles)

[miaocb]

Hi, Ian,
Thanks for your reply. The software version I used is as follows:
Python 2.7.8 :: Anaconda 2.1.0 (64-bit)
numpy-1.9.0-py27_0
matplotlib-1.4.0-np19py27_0

The output of print:
('ele1', array([[2, 0, 3], [2, 1, 0]]))
('ele2', array([[2, 0, 3], [2, 1, 0]]))
('ele1.shape', (2, 3))
('ele2.shape', (2, 3))
('ele1.strides', (8, 16))
('ele2.strides', (24, 8))
('fine1', array([[2, 4, 5],
[0, 4, 6],
[3, 6, 5],
[4, 5, 6],
[2, 8, 5],
[1, 5, 7],
[0, 7, 8],
[5, 8, 7]], dtype=int32))
('fine2', array([[2, 4, 5],
[3, 6, 4],
[0, 5, 6],
[4, 6, 5],
[2, 5, 8],
[0, 7, 5],
[1, 8, 7],
[5, 7, 8]], dtype=int32))

Obviously, fine1 is different with fine2, I use numpy 1.9.0 and python 2.7.8, while you used numpy 1.8.2 and python 2.7.6.

[tacaswell]

@miaocb The images did not make it through.

[miaocb]

[ianthomas23]

I have confirmed that the bug occurs with mpl 1.4.x and numpy 1.9.x. There is no problem with mpl 1.4.x and numpy 1.8.x, and mpl master is fine regardless of numpy version.

I recall that numpy 1.9 required quite a few changes to deal with slightly different python and C interfaces. This is something we must have missed at the time. The C++ triangulation code deals with contiguous numpy arrays and it assumes they are ordered as conventional C arrays, i.e. C-contiguous. A fortran-contiguous array (which is what this issue uses) is correctly converted to C-contiguous by our code for numpy < 1.9, but not for numpy >= 1.9. In the master branch, the rewrite to remove PyCXX has fortuitously removed the bug for us; this is good, but is probably why the bug has not been spotted before.

If I remember correctly, we won't be making any more releases from the 1.4.x branch as the next release will be 1.5 or 2.0 from the master branch. If so, there is no point in me fixing the bug in the 1.4.x branch as it will never be released. Can you confirm @tacaswell?

@miaocb: Thanks for reporting this bug. Until we make a new release you have a number of options to deal with this bug:

1. Use the mpl master branch, but this will involve compiling mpl from scratch.
2. Use mpl with an older version of numpy (1.8.x).
3. Whenever you use a transposed triangles array to create a tri.Triangulation, make sure it has been copied first as this will ensure the array is C-contiguous. For example, in your test script change line 10 to ele1 = ele0.transpose().copy()
4. The best solution is to alter the source code to eliminate the bug. Find the file tri/triangulation.py in your mpl installation directory. Around line 60 there is the following line:
   self.triangles = np.array(triangles, dtype=np.int32)
   Change this to
   self.triangles = np.array(triangles, dtype=np.int32, order='C')

[tacaswell]

2.0 (the color over haul) will be based off of the 1.4 branch rather than
master so it is worth fixing on the color_overhaul branch.

Once 2.0 is released the master branch will become 2.1

On Tue, Mar 3, 2015, 04:51 Ian Thomas notifications@github.com wrote:

I have confirmed that the bug occurs with mpl 1.4.x and numpy 1.9.x. There
is no problem with mpl 1.4.x and numpy 1.8.x, and mpl master is fine
regardless of numpy version.

I recall that numpy 1.9 required quite a few changes to deal with slightly
different python and C interfaces. This is something we must have missed at
the time. The C++ triangulation code deals with contiguous numpy arrays and
it assumes they are ordered as conventional C arrays, i.e. C-contiguous. A
fortran-contiguous array (which is what this issue uses) is correctly
converted to C-contiguous by our code for numpy < 1.9, but not for numpy >=
1.9. In the master branch, the rewrite to remove PyCXX has fortuitously
removed the bug for us; this is good, but is probably why the bug has not
been spotted before.

If I remember correctly, we won't be making any more releases from the
1.4.x branch as the next release will be 1.5 or 2.0 from the master branch.
If so, there is no point in me fixing the bug in the 1.4.x branch as it
will never be released. Can you confirm @tacaswell
https://github.com/tacaswell?

@miaocb https://github.com/miaocb: Thanks for reporting this bug. Until
we make a new release you have a number of options to deal with this bug:

1. Use the mpl master branch, but this will involve compiling mpl from
   scratch.
2. Use mpl with an older version of numpy (1.8.x).
3. Whenever you use a transposed triangles array to create a
   tri.Triangulation, make sure it has been copied first as this will
   ensure the array is C-contiguous. For example, in your test script change
   line 10 to ele1 = ele0.transpose().copy()
4. The best solution is to alter the source code to eliminate the bug.
   Find the file tri/triangulation.py in your mpl installation directory.
   Around line 60 there is the following line: self.triangles =
   np.array(triangles, dtype=np.int32) Change this to self.triangles =
   np.array(triangles, dtype=np.int32, order='C')

—
Reply to this email directly or view it on GitHub
#4180 (comment)
.

[jenshnielsen]

I guess it would make sense to turn this into a test too?

[miaocb]

Thanks. I changed the code as suggested.
"Around line 60 there is the following line: self.triangles = np.array(triangles, dtype=np.int32) Change this to self.triangles = np.array(triangles, dtype=np.int32, order='C')"