Matplotlib with invalid triangulations

Hello,

I am using JR Shewchuk's Triangle to create triangulations for use in
floodplain modeling. I am using a city's digital terrain model with
hundreds of thousands of breaklines that constrain where triangles can
form in the triangulations (streets, rivers, etc). Triangle does this
very efficiently.

Sometimes the input topology I am using has bad inputs which makes
Triangle create zero area elements. When I import these triangulations
to Matplotlib I get the error that such triangulations are invalid (when
using the LinearTriInterpolator() method). I understand the trapezoid
map algorithm implemented requires only valid triangulations. So far, so
good.

The option of cleaning the input topology before using Matplotlib
exists, but it is really cumbersome. Rather than topology cleaning, am I
am able to somehow throw out the zero area elements from the call to
LinearTriInterpolator() method, and still use the triangulation in
Matplotlib? My other alternative is to use something other than
trapezoidal map algorithm, but this is just not computationally efficient.

I've reproduced the following example that illustrates the problem in a
small code snippet. Any suggestions?

Thanks,

Pat Prodanovic

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

Hello Pat,

The solution is to use the function Triangulation.set_mask() to mask out
the zero-area triangles. The masked-out triangles will be ignored in
subsequent calls to LinearTriInterpolator, tricontourf, etc. For example:

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# ---------- start of new code ----------
xy = np.dstack((triang.x[triang.triangles], triang.y[triang.triangles])) #
shape (ntri,3,2)
twice_area = np.cross(xy[:,1,:] - xy[:,0,:], xy[:,2,:] - xy[:,0,:]) #
shape (ntri)
mask = twice_area < 1e-10 # shape (ntri)

if np.any(mask):
    triang.set_mask(mask)
# ---------- end of new code ----------

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

Note that I have used a small positive number to test the triangle areas
against rather than zero. This is to avoid problems with rounding errors.
You may need to alter this threshold.

Ian

···

On 19 November 2016 at 12:24, Pat Prodanovic <pprodano at gmail.com> wrote:

Hello,

I am using JR Shewchuk's Triangle to create triangulations for use in
floodplain modeling. I am using a city's digital terrain model with
hundreds of thousands of breaklines that constrain where triangles can form
in the triangulations (streets, rivers, etc). Triangle does this very
efficiently.

Sometimes the input topology I am using has bad inputs which makes
Triangle create zero area elements. When I import these triangulations to
Matplotlib I get the error that such triangulations are invalid (when using
the LinearTriInterpolator() method). I understand the trapezoid map
algorithm implemented requires only valid triangulations. So far, so good.

The option of cleaning the input topology before using Matplotlib exists,
but it is really cumbersome. Rather than topology cleaning, am I am able to
somehow throw out the zero area elements from the call to
LinearTriInterpolator() method, and still use the triangulation in
Matplotlib? My other alternative is to use something other than trapezoidal
map algorithm, but this is just not computationally efficient.

I've reproduced the following example that illustrates the problem in a
small code snippet. Any suggestions?

Thanks,

Pat Prodanovic

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users at python.org
https://mail.python.org/mailman/listinfo/matplotlib-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20161121/e966b5b7/attachment.html>

Hi Ian,

Thank you for your reply.

The modification you provided correctly finds the zero area element, and
masks it from the triangulation. In the example from the previous post,
masking the zero area element works.

When I try and make a slightly different triangulation (see below), and
try to mask the zero area elements, I still get an invalid
triangulation. I am using v1.4.2. Do you have a sense as to what could
be going on here?

Thanks,

Pat

import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

# slightly modified from what I originally posted
triangles = np.array( [[0, 1, 4], [2, 3, 4], [0, 3, 2], [0, 4, 3]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# ---------- start of new code ----------
xy = np.dstack((triang.x[triang.triangles], triang.y[triang.triangles]))
#shape (ntri,3,2)
twice_area = np.cross(xy[:,1,:] - xy[:,0,:], xy[:,2,:] - xy[:,0,:]) #
shape (ntri)
mask = twice_area < 1e-10 # shape (ntri)

if np.any(mask):
     triang.set_mask(mask)
# ---------- end of new code ----------

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)

···

On 11/21/2016 03:47 AM, Ian Thomas wrote:

Hello Pat,

The solution is to use the function Triangulation.set_mask() to mask
out the zero-area triangles. The masked-out triangles will be ignored
in subsequent calls to LinearTriInterpolator, tricontourf, etc. For
example:

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# ---------- start of new code ----------
xy = np.dstack((triang.x[triang.triangles],
triang.y[triang.triangles])) # shape (ntri,3,2)
twice_area = np.cross(xy[:,1,:] - xy[:,0,:], xy[:,2,:] - xy[:,0,:]) #
shape (ntri)
mask = twice_area < 1e-10 # shape (ntri)

if np.any(mask):
    triang.set_mask(mask)
# ---------- end of new code ----------

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

Note that I have used a small positive number to test the triangle
areas against rather than zero. This is to avoid problems with
rounding errors. You may need to alter this threshold.

Ian

On 19 November 2016 at 12:24, Pat Prodanovic <pprodano at gmail.com > <mailto:pprodano at gmail.com>> wrote:

    Hello,

    I am using JR Shewchuk's Triangle to create triangulations for use
    in floodplain modeling. I am using a city's digital terrain model
    with hundreds of thousands of breaklines that constrain where
    triangles can form in the triangulations (streets, rivers, etc).
    Triangle does this very efficiently.

    Sometimes the input topology I am using has bad inputs which makes
    Triangle create zero area elements. When I import these
    triangulations to Matplotlib I get the error that such
    triangulations are invalid (when using the LinearTriInterpolator()
    method). I understand the trapezoid map algorithm implemented
    requires only valid triangulations. So far, so good.

    The option of cleaning the input topology before using Matplotlib
    exists, but it is really cumbersome. Rather than topology
    cleaning, am I am able to somehow throw out the zero area elements
    from the call to LinearTriInterpolator() method, and still use the
    triangulation in Matplotlib? My other alternative is to use
    something other than trapezoidal map algorithm, but this is just
    not computationally efficient.

    I've reproduced the following example that illustrates the problem
    in a small code snippet. Any suggestions?

    Thanks,

    Pat Prodanovic

    #
    +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
    import matplotlib.tri as mtri
    import numpy as np

    # manually construct an invalid triangulation having a zero area
    element
    x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
    y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
    z = np.zeros(5)

    triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

    # create a Matplotlib Triangulation
    triang = mtri.Triangulation(x,y,triangles)

    # to perform the linear interpolation
    interpolator = mtri.LinearTriInterpolator(triang, z)
    m_z = interpolator(1.0, 1.0)
    #
    +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

    _______________________________________________
    Matplotlib-users mailing list
    Matplotlib-users at python.org <mailto:Matplotlib-users at python.org>
    https://mail.python.org/mailman/listinfo/matplotlib-users
    <https://mail.python.org/mailman/listinfo/matplotlib-users>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20161121/f736ab39/attachment-0001.html>

Hi Pat,

That's really unlucky!

Point 3 at (1,1) lies along one edge of triangle (0, 1, 4). However, due
to rounding error caused by finite-precision maths, there is some test when
creating LinearTriInterpolator for which point 3 lies just inside triangle
(0, 1, 4). Hence the triangulation is invalid. You can prove this by
changing the line
    y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
to
    y = np.array([1.0, 0.0, 2.0, 1.0+1e-10, 1.0])
instead.

Really the only solution if you want to use LinearTriInterpolator is to
ensure that you don't have any very thin triangles before going anywhere
near matplotlib. This is down to whatever you use to generate your
triangulation and/or whatever preprocessing you perform on it.

Ian

···

On 21 November 2016 at 11:45, Pat Prodanovic <pprodano at gmail.com> wrote:

Hi Ian,

Thank you for your reply.

The modification you provided correctly finds the zero area element, and
masks it from the triangulation. In the example from the previous post,
masking the zero area element works.

When I try and make a slightly different triangulation (see below), and
try to mask the zero area elements, I still get an invalid triangulation. I
am using v1.4.2. Do you have a sense as to what could be going on here?

Thanks,

Pat

import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

# slightly modified from what I originally posted
triangles = np.array( [[0, 1, 4], [2, 3, 4], [0, 3, 2], [0, 4, 3]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# ---------- start of new code ----------
xy = np.dstack((triang.x[triang.triangles], triang.y[triang.triangles]))
#shape (ntri,3,2)
twice_area = np.cross(xy[:,1,:] - xy[:,0,:], xy[:,2,:] - xy[:,0,:]) #
shape (ntri)
mask = twice_area < 1e-10 # shape (ntri)

if np.any(mask):
    triang.set_mask(mask)
# ---------- end of new code ----------

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
On 11/21/2016 03:47 AM, Ian Thomas wrote:

Hello Pat,

The solution is to use the function Triangulation.set_mask() to mask out
the zero-area triangles. The masked-out triangles will be ignored in
subsequent calls to LinearTriInterpolator, tricontourf, etc. For example:

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# ---------- start of new code ----------
xy = np.dstack((triang.x[triang.triangles], triang.y[triang.triangles]))
# shape (ntri,3,2)
twice_area = np.cross(xy[:,1,:] - xy[:,0,:], xy[:,2,:] - xy[:,0,:]) #
shape (ntri)
mask = twice_area < 1e-10 # shape (ntri)

if np.any(mask):
    triang.set_mask(mask)
# ---------- end of new code ----------

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

Note that I have used a small positive number to test the triangle areas
against rather than zero. This is to avoid problems with rounding errors.
You may need to alter this threshold.

Ian

On 19 November 2016 at 12:24, Pat Prodanovic <pprodano at gmail.com> wrote:

Hello,

I am using JR Shewchuk's Triangle to create triangulations for use in
floodplain modeling. I am using a city's digital terrain model with
hundreds of thousands of breaklines that constrain where triangles can form
in the triangulations (streets, rivers, etc). Triangle does this very
efficiently.

Sometimes the input topology I am using has bad inputs which makes
Triangle create zero area elements. When I import these triangulations to
Matplotlib I get the error that such triangulations are invalid (when using
the LinearTriInterpolator() method). I understand the trapezoid map
algorithm implemented requires only valid triangulations. So far, so good.

The option of cleaning the input topology before using Matplotlib exists,
but it is really cumbersome. Rather than topology cleaning, am I am able to
somehow throw out the zero area elements from the call to
LinearTriInterpolator() method, and still use the triangulation in
Matplotlib? My other alternative is to use something other than trapezoidal
map algorithm, but this is just not computationally efficient.

I've reproduced the following example that illustrates the problem in a
small code snippet. Any suggestions?

Thanks,

Pat Prodanovic

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
import matplotlib.tri as mtri
import numpy as np

# manually construct an invalid triangulation having a zero area element
x = np.array([0.0, 1.0, 1.0, 1.0, 2.0])
y = np.array([1.0, 0.0, 2.0, 1.0, 1.0])
z = np.zeros(5)

triangles = np.array( [[0, 1, 2], [1, 3, 2], [1, 4, 2], [0, 4, 1]])

# create a Matplotlib Triangulation
triang = mtri.Triangulation(x,y,triangles)

# to perform the linear interpolation
interpolator = mtri.LinearTriInterpolator(triang, z)
m_z = interpolator(1.0, 1.0)
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users at python.org
https://mail.python.org/mailman/listinfo/matplotlib-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20161122/2afcf7aa/attachment.html>