Here we will implement an algorithm for breaking a polygon into triangles which can be used for any planar polygon (not necessarily convex).

First we define several auxilary functions. Two functions below are for comparing machine numbers. feq is defined in such a way so as to work for real numbers, complex numbers and vectors (i.e. lists).

\(feq[x_, y_, eps_:  10^\(-10\)] := 
	Plus\ @@\ Abs[x - y] < eps\)

fsign[x_] := If[feq[x, 0], 0, Sign[x]]

interq checks if two given segments (not lines!) intersect; onsegq determines if a point lies on a segment. The general idea is to construct the equation of a line and then see on which side of this line certain points are. interq can give incorrect results if an end of one segment lies on the line containing the other segment but we'll get around it.

\(interq[{{x11_, y11_}, {x12_, y12_}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{x21_, y21_}, {x22_, y22_}}] := 
      Module[
		{line1, line2}, 
		
        line1[x_, y_] = 
          x \((y12 - y11)\) + y \((x11 - x12)\) - \((x11\ y12 - x12\ y11)\); 
        
		line2[x_, y_] = 
          x \((y22 - y21)\) + y \((x21 - x22)\) - \((x21\ y22 - x22\ y21)\); 
        
		Return[
			
          fsign[line1[x21, y21] line1[x22, y22]] < 0 && 
				
            fsign[line2[x11, y11] line2[x12, y12]] < 0]]\)

\(onsegq[p : {x_, y_}, {p1 : {x1_, y1_}, p2 : {x2_, y2_}}] := 
	
      fsign[x \((y2 - y1)\) + y \((x1 - x2)\) - \((x1\ y2 - x2\ y1)\)] == 0
           && 
		feq[\(Sort[
				{p, p1, p2}]\)\ [\([2]\)], p]\)

Now, the important fact is that we can compute the area of a polygon by coordinates of its vertices. It is a particular case of the general formula for the oriented area surrounded by a closed contour C:

S = (1/2) (∮_C) x dy - 
              y dx

If we break this integral into sum of integrals along each side of the polygon and write parametric equations for each side (segment), we'll get the following formula:

\(TraditionalForm
      \`S = \(1\/2\) 
          \(∑\+\(k = 1\)\%n\((
              \(x\_k\) y\_\(k + 1\) - \(x\_\(K + 1\)\) y\_k)\)\)\),

where ((Subscript[x, i, ]) (Subscript[y, i])) are apexes of the n-gon and we define \(TraditionalForm
      \`\((x\_\(n + 1\), y\_\(n + 1\))\) = \((x\_1, y\_1)\)\).

\(area[poly_] := 
	
      Abs[Plus\ @@
			
            MapThread[Dot, 

              				{poly, 
					
                \({#[\([2]\)], \(-#[\([1]\)]\)}&\)\ /@
						
                  RotateLeft[poly]}\ ]\ ]/2\)

Now for the algorithm itself. On each step we choose a diagonal of our polygon P to split it into two polygons (with smaller number of vertices) We need to check for two things to determine the validity of the diagonal. First, we make sure that the diagonal doesn't intersect with the sides of the polygon  If it doesn't, then it breaks P into two correct polygons (Subscript[P, 1]) and (Subscript[P, 2]) (i.e. closed polylines without self-intersections). Also, such diagonal lies either entirely inside of P or entirely outside. And then we only need to check if the area of P is equal to the sum of areas of (Subscript[P, 1]) and (Subscript[P, 2]) (if it's not, then the diagonal lies outside of the polygon P and one of the polygons (Subscript[P, 1]), (Subscript[P, 2]) lies inside the other.

\(split[poly_] := 
      Module[
		{n, s, sides, 
			diag, inters, apxs, 
			poly1, 
          poly2}, 
		
        If[\((n = Length[poly])\) == 3, 
			Return[{poly}]\ ]; 
		
        s = area[poly]; 
		sides = Transpose[{poly, RotateLeft[poly]}]; 

        		Do[
			If[{ap1, ap2} == {1, n}, 
				Continue[]\ ]; 

          			diag = {poly[\([ap1]\)], poly[\([ap2]\)]}; 
			
          inters = Select[sides, 
					interq[#, diag]&, 1]; 
			
          If[inters != {}, 
				Continue[]\ ]; 
			
          apxs = Delete[poly, 
					{{ap1}, {ap2}}\ ]; 
			
          apxs = Select[apxs, 
					onsegq[#, diag]&, 1]; 
			
          If[apxs != {}, 
				Continue[]\ ]; 
			
          poly1 = Take[
					RotateLeft[poly, ap1 - 1], 
					
              ap2 - ap1 + 1]; 
			
          poly2 = Take[
					RotateLeft[poly, ap2 - 1], 
					
              n + 2 - \((ap2 - ap1 + 1)\)]; 
			
          If[\(! feq[area[poly1] + area[poly2], s]\), 
				Continue[]]; 
          
			Return[Join[
					split[poly1], split[poly2]]\ ], 

          			{ap1, n - 2}, {ap2, ap1 + 2, n}]\ ]\)

In the same way we can solve the problem of finding the convex hull of a polygon. In this case we have to look for a diagonal that lies outside of the polygon. If there isn't such a diagonal, P is convex. If we find one, we replace P with the largest one of (Subscript[P, 1]) and (Subscript[P, 2]), and so our polygon will 'grow' until we construct a convex polygon.

\(hull[poly_] := 
      Module[
		{n, s, s1, s2, sides, 
			diag, inters, apxs, 
			
          poly1, poly2, res = {}}, 
		
        If[\((n = Length[poly])\) == 3, 
			Return[poly]]; 
		
        s = area[poly]; 
		sides = Transpose[{poly, RotateLeft[poly]}]; 

        		Do[
			If[{ap1, ap2} == {1, n}, 
				Continue[]\ ]; 

          			diag = {poly[\([ap1]\)], poly[\([ap2]\)]}; 
			
          inters = Select[sides, 
					interq[#, diag]&, 1]; 
			
          If[inters != {}, 
				Continue[]\ ]; 
			
          apxs = Delete[poly, 
					{{ap1}, {ap2}}\ ]; 
			
          apxs = Select[apxs, 
					onsegq[#, diag]&, 1]; 
			
          If[apxs != {}, 
				Continue[]\ ]; 
			
          poly1 = Take[
					RotateLeft[poly, ap1 - 1], 
					
              ap2 - ap1 + 1]; 
			
          poly2 = Take[
					RotateLeft[poly, ap2 - 1], 
					
              n + 2 - \((ap2 - ap1 + 1)\)]; 

          			{s1, s2} = area\ /@\ {poly1, poly2}; 
			
          If[feq[s1 + s2, s], 
				Continue[]]; 
			
          Return[res = hull[
					If[s1 > s2, poly1, poly2]\ ]\ ], 

          			{ap1, n - 2}, {ap2, ap1 + 2, n}]; 
		
        If[res == {}, poly, res]\ ]\)

Finally, this algorithm can be applied to a planar polygon in space. In order to do this, we may find a vector normal to the plane containing the polygon P, rotate P so as to make it parallel to horizontal plane and then discard the third coordinate of P's apexes to reduce the problem to the previous one.

rotm[{x, y, z}] returns a rotation matrix A such that A . {x, y, z} == {0, 0, z'} (we need the normal vector to be parallel to z axis).

rotm[{x_, y_, z_}]:=Module[
		{α, β},
		{α, β} = {ArcTan[z, y], ArcTan[(Sqrt[(y^2) + (z^2)]), x]}
                ;
		(({{Cos[β], 0, -Sin[β]}, {0, 1, 0}, {Sin[β], 0, Cos[β]}})).(({{1, 0, 0}, {0, Cos[α], -Sin[α]}, {0, Sin[α], Cos[α]}}))]

\(split3d[poly_] := 
      Module[
		{vec, zs, 
			rot, flat, triangs}, 
		
        vec = \((poly[\([3]\)] - poly[\([1]\)])\)\[Cross]\((
              poly[\([2]\)] - poly[\([1]\)])\); 
		
        If[feq[vec, 0], 
			Print["Unable to find normal vector"]; 

          			Return[]\ ]; 
		rot = \(rotm[vec].#&\)\ /@\ poly; 
		
        zs = \(#[\([3]\)]&\)\ /@\ rot; 
		
        If[Length[
				Split[zs, feq]\ ] != 1, 
			
          Print["Polygon is not planar"]; 
			Return[]\ ]; 
		
        flat = \(Take[#, 2]&\)\ /@\ rot; 
		triangs = split[flat]; 
		
        triangs = 
          Map[
				Append[#, zs[\([1]\)]\ ]&, 
				triangs, {2}]; 

        		Map[Transpose[rotm[vec]].#&, 
			triangs, {2}]\ ]\)

Here is a test. Polygon poly is a star rotated by a random angle in space. It can be visualized incorrectly but triangles are guaranteed to be displayed right.

\(poly = 
      Table[{Sin[\(π\ k\)\/7], Cos[\(π\ k\)\/7], 0} 
          \((2 + Mod[k, 2])\), {k, 14}]; 

    poly = \(rotm[{1, 2, 3}].#&\)\ /@\ N[poly]\)

Evaluate the cells below to see the results.

\(Show[Graphics3D[
		Line[Append[poly, poly[\([1]\)]\ ]]\ ]]\)

\(Show[Graphics3D[
		Polygon[poly]\ ]]\)

poly2 = split3d[poly];

\(Show[Graphics3D[
		Polygon\ /@\ poly2\ ]]\)

\(Show[Graphics3D[
		{EdgeForm[], Polygon\ /@\ poly2}\ ]]\)

Copyright © 1998, Maxim Rytin