# Matplotlib with invalid triangulations

 Classic List Threaded
4 messages
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

## 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) # +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- _______________________________________________ Matplotlib-users mailing list [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Matplotlib with invalid triangulations

 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 mtriimport numpy as np# manually construct an invalid triangulation having a zero area elementx = 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 Triangulationtriang = 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 interpolationinterpolator = 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.IanOn 19 November 2016 at 12:24, Pat Prodanovic 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 [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users _______________________________________________ Matplotlib-users mailing list [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Matplotlib with invalid triangulations

 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 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 [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users _______________________________________________ Matplotlib-users mailing list [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Matplotlib with invalid triangulations

 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.IanOn 21 November 2016 at 11:45, Pat Prodanovic 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 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 [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users _______________________________________________ Matplotlib-users mailing list [hidden email] https://mail.python.org/mailman/listinfo/matplotlib-users
Loading...