Python Forum

Dear Forum memebers,

i want to calcualte the surface area of each triangle generated by the Delaunay Triangulation by the following python code.

could you please support with this issue?

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


axes=plt.axes()
x=np.random.random(5)
y=np.random.random(5)
axes.set_xlim([0,1])
axes.set_ylim([0,1])

# plt.scatter(x,y)
# plt.show()

plt.scatter(x,y)
plt.show()

# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

thank you in advance Big Grin

Vertices for each triangle (in a matrix) are provided after tesselation; just add the following code:

VerticesNodesNumbers = triang.triangles
print("Vertices of triangles;\n{}".format(VerticesNodesNumbers))

From vertices number, you've:

to get nodes coordinates to calculate the area of each triangle
the complete surface area is no more no less than the summ of all triangles areas

see here for examples of formulas

Thank you so much for your kind help.

I tried the following code. it seems to work for huge number of points.

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

np.random.seed(0)
axes=plt.axes()
x=np.random.random(10000)
y=np.random.random(10000)
axes.set_xlim([0,1])
axes.set_ylim([0,1])

plt.scatter(x,y)
plt.show()

# # Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
#plt.triplot(triang, 'ro:')

AZ=triang.neighbors
print(AZ)

# VerticesNodesNumbers = VNN
VNN = triang.triangles
print(VNN)
a,b=VNN.shape
print(a)
#print("x\n{}\n".format(x))
#print("y\n{}\n".format(y))

Arealiste=[]
# in order to calculate the area of each single triangle, we need to scann
# the VNN 2D Array by the number of raws only. 

for i in range(a):
    t1=(x[VNN[i][0]])*((y[VNN[i][1]])-(y[VNN[i][2]]))
    t2=(x[VNN[i][1]])*((y[VNN[i][2]])-(y[VNN[i][0]]))
    t3=(x[VNN[i][2]])*((y[VNN[i][0]])-(y[VNN[i][1]]))
    Arealiste.append(abs((t1+t2+t3)/2))

#print(Arealiste)

Areas=np.array(Arealiste)
#print(Areas)
print(np.sum(Areas))
k=(Areas/(np.sum(Areas)))*100
#print(k)
# we can now plot the histogram of the area distribution of the triangles:-)
plt.hist(k,1000)
plt.show()

Could you please check?
I appreciate your support for any further hints :-)

You can easily check using the example of the website I gave you.

My advices: use vectorization instead of a loop, and above all, do not use append (dynamic memory allocation is time consuming); have a look to the code here bellow as an example.

Of course by replacing and increasing the n value, the area tends to 1

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
 
axes=plt.axes()
n = 10_000
# x=np.random.random(n)
# y=np.random.random(n)
x= np.array( [ 15, 23, 50] ) # see https://www.mathopenref.com/coordtrianglearea.html 
y= np.array( [ 15, 30, 25] ) # Area = 222.5 for these 3 vertices

axes.set_xlim([0,1])
axes.set_ylim([0,1])
 
# plt.scatter(x,y)
# plt.show()
 
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: Area = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) | where a,b,c are 3 vertices
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Area = {}".format(Area))

Just for fun

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
R=2 
print("Circle of radius R=2")
x= np.array( [2, 1.8270909152852, 1.33826121271772, 0.618033988749899, \
              -0.209056926535301, -0.999999999999994, -1.61803398874989, \
              -1.95629520146761, -1.95629520146761, -1.6180339887499,  \
              -1.00000000000001, -0.209056926535314, 0.618033988749889, \
              1.33826121271771, 1.8270909152852, -0.856522911949342,   \
              0.635640538230626, -1.15448031396626, 1.00907475436824,  \
              0.11104325270382, -0.425564794983423, 1.41862257999745,  \
              0.150206110091391, 0.042749743825051, 0.717216890750844, \
              0.23402642946099, -1.47849160060644, -0.668068810539202, \
              0.862760612586134, -1.36476212412967, -0.996086145254554, \
              -0.457339674908644, 1.42321207043604, 1.22013478402029,  \
              -0.582216221137397, -0.285886702570627] )
    
y= np.array( [0, 0.813473286151599, 1.48628965095479, 1.90211303259031, \
              1.98904379073655, 1.73205080756888, 1.17557050458495,     \
              0.415823381635526, -0.415823381635509, -1.17557050458494, \
              -1.73205080756887, -1.98904379073655, -1.90211303259031,  \
              -1.48628965095479, -0.813473286151604, 0.938050975700476, \
              -1.05863282037826, -0.508934656032994, 0.717164036805174, \
              1.2944564466239, -1.22019130498736, -0.298702466121143,   \
              -1.39814565289391, -0.406736643075798, 1.23915884760709,  \
              0.405345666140318, 0.004667930500299, 0.149680059053867,  \
              -0.174425526495497, 0.594557386209038, -1.11495929089071, \
              1.42228369708499, 0.203021725197203, -0.75790529479237,   \
              -0.60844562355935, 0.800289654327748] )

axes=plt.axes()
axes.set_xlim([0,1])
axes.set_ylim([0,1])
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: A = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) |
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Circle area = {} vs theoritical value: {}".format(Area, np.pi*R**2))

(of course the accuracy depends on the mesh discretization Wink

)

(Jan-18-2021, 08:39 AM)paul18fr Wrote: [ -> ]You can easily check using the example of the website I gave you.

My advices: use vectorization instead of a loop, and above all, do not use append (dynamic memory allocation is time consuming); have a look to the code here bellow as an example.

Of course by replacing and increasing the n value, the area tends to 1

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
 
axes=plt.axes()
n = 10_000
# x=np.random.random(n)
# y=np.random.random(n)
x= np.array( [ 15, 23, 50] ) # see https://www.mathopenref.com/coordtrianglearea.html 
y= np.array( [ 15, 30, 25] ) # Area = 222.5 for these 3 vertices

axes.set_xlim([0,1])
axes.set_ylim([0,1])
 
# plt.scatter(x,y)
# plt.show()
 
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: Area = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) | where a,b,c are 3 vertices
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Area = {}".format(Area))

I do not understand the reason for saying "Of course by replacing and increasing the n value, the area tends to 1" and so for having n = 10_000. Could you elaborate more?

Thanks,
Valerio

Hi
Probably a typo. Anyway just focuss on the example and the code.

Paul

pythonformaterialsci

paul18fr

pythonformaterialsci

paul18fr

paul18fr

meValerio

paul18fr