Numpy Area of Triangle and equation of a plane on which triangle lies on-CodePudding

I have written code to find area of triangle, and the plane of which it lies on. But it feels pretty slow. Is there a way to use numpy array to find area of the triangle and plane equation with out using For loops and store it into a new array.

i have stored my coordinates in 3D array.

array([[[405.05, 805.52, 345.1 ],
        [406.78, 805.73, 314.15],
        [405.84, 805.62, 345.14]],

       [[405.84, 805.62, 345.14],
        [406.78, 805.73, 314.15],
        [407.57, 805.83, 314.19]],

       [[407.57, 805.83, 314.19],
        [406.78, 805.73, 314.15],
        [407.9 , 805.19, 309.75]]])

def area(self):
        side1=length(self.a,self.b)
        side2=length(self.b,self.c)
        side3=length(self.c,self.a)
        try:
                 s=(side1 side2 side3)/2
                 self.area1=sqrt(s*(s-side1)*(s-side2)*(s-side3))
        except ValueError:
                self.area1=0
        return self.area1

def length(a,b):
         l1=sqrt(((a[0]-b[0])**2) ((a[1]-b[1])**2) ((a[2]-b[2])**2))
         return l1

def Vector_values(Coordinates):
         """For calculation of Vector PQ and PR for the plane equation """
         arr1=np.asarray(Coordinates[0])
         arr2=np.asarray(Coordinates[1])
         arr3=arr1-arr2
         return arr3

def plane_equation(Vertex_triangle):
         """For calculation the plane equation """
         PQ=Vector_values([Vertex_triangle[0],Vertex_triangle[1]])
         PR=Vector_values([Vertex_triangle[0],Vertex_triangle[2]])
         ijk=np.vstack((PQ,PR))
         i=round(np.linalg.det(ijk[0:2,1:3]),2)
         j=-round(np.linalg.det(ijk[0:2,::2]),2)
         k=round(np.linalg.det(ijk[0:2,0:2]),2)
         d=-round((i*Vertex_triangle[0][0] j*Vertex_triangle[0][1] k*Vertex_triangle[0][2]),2)
         return [i,j,k,d]

CodePudding user response：

You can calculate areas for multiple triangles this way:

def normal(triangles):
    # The cross product of two sides is a normal vector
    return np.cross(triangles[:,1] - triangles[:,0], 
                    triangles[:,2] - triangles[:,0], axis=1)

def area(triangles):
    # The norm of the cross product of two sides is twice the area
    return np.linalg.norm(normal(triangles), axis=1) / 2

In your example:

>>> area(triangles)
array([12.35765785, 12.35765785,  1.79312445])

And for the planes:

def plane(triangles):
    # The unit normal u can be used as the coefficients a,b,c
    n = normal(triangles)
    u = n / np.linalg.norm(n, axis=1, keepdims=True)
    # Coefficient d is the dot product between a point on each plane
    # and each unit normal (with opposite sign)
    d = -np.einsum('ij,ij->i', triangles[:, 0], u)
    return np.hstack((u, d[:, None]))

In your example:

>>> plane(triangles)
array([[ 1.25565865e-01, -9.92085244e-01,  2.87271265e-04, 7.48184915e 02],
       [ 1.25565865e-01, -9.92085244e-01,  2.87271265e-04, 7.48184915e 02],
       [ 1.16667864e-01, -9.81750040e-01,  1.50184779e-01, 6.96386758e 02]])