Use Delauney triangulation to minimize minimum angle

src/dtm-surface.adb

@@ -1,5 +1,6 @@
 with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
 with Ada.Text_IO; use Ada.Text_IO;
+with Ada.Numerics.Generic_Elementary_Functions; use Ada.Numerics;
 
 package body Dtm.Surface is
    Zero : constant Dtm_Value := Dtm_Value(0);
@@ -123,7 +124,6 @@ package body Dtm.Surface is
       Add_Left(SV, G);
       Add_Right(SV, G);
    end Add_Sides;
-
    
    procedure Constrain(SV : in out Surface_Vector; Scale_Factor : Dtm_Value) is
    begin
@@ -134,13 +134,20 @@ package body Dtm.Surface is
          end loop;
       end loop;
    end Constrain;
-   
+
+   function Vector_Length(V : Point) return Dtm_Value is
+      package Dtm_Math is new Generic_Elementary_Functions(Dtm.Dtm_Value); use Dtm_Math;
+   begin
+      return Sqrt(V.X**2 + V.Y**2 + V.Z**2);
+   end Vector_Length;
+        
    procedure Surface_From_Grid(Triangles : in out Surface_Vector; Data : Dtm_Data) is
       Grid : Dtm_Grid := Data.Grid;
    begin
       for X in Grid'First(1) .. Grid'Last(1) - 1 loop
          for Y in Grid'First(2) .. Grid'Last(2) - 1 loop
             declare
+               -- Option 1
                T1_P1 : Point := (Dtm_Value(X), Dtm_Value(Y), Grid(X, Y));
                T1_P2 : Point := (Dtm_Value(X+1), Dtm_Value(Y), Grid(X+1, Y));
                T1_P3 : Point := (Dtm_Value(X+1), Dtm_Value(Y+1), Grid(X+1, Y+1));
@@ -150,10 +157,29 @@ package body Dtm.Surface is
                T2_P2 : Point := (Dtm_Value(X+1), Dtm_Value(Y+1), Grid(X+1, Y+1));
                T2_P3 : Point := (Dtm_Value(X), Dtm_Value(Y+1), Grid(X, Y+1));
                Tri2 : Triangle := (T2_P1, T2_P2, T2_P3);
+
+               O1_Length : Dtm_Value := Vector_Length(T1_P1 - T1_P3);
+               
+               -- Option 2
+               T3_P1 : Point := (Dtm_Value(X), Dtm_Value(Y), Grid(X, Y));
+               T3_P2 : Point := (Dtm_Value(X+1), Dtm_Value(Y), Grid(X+1, Y));
+               T3_P3 : Point := (Dtm_Value(X), Dtm_Value(Y+1), Grid(X, Y+1));
+               Tri3 : Triangle := (T3_P1, T3_P2, T3_P3);
+               
+               T4_P1 : Point := (Dtm_Value(X+1), Dtm_Value(Y), Grid(X+1, Y));
+               T4_P2 : Point := (Dtm_Value(X+1), Dtm_Value(Y+1), Grid(X+1, Y+1));
+               T4_P3 : Point := (Dtm_Value(X), Dtm_Value(Y+1), Grid(X, Y+1));
+               Tri4 : Triangle := (T4_P1, T4_P2, T4_P3);
+               
+               O2_Length : Dtm_Value := Vector_Length(T3_P2 - T3_P3);
             begin
-               null;
-               Triangles.Append(Tri1);
-               Triangles.Append(Tri2);
+               if O1_Length < O2_Length then
+                  Triangles.Append(Tri1);
+                  Triangles.Append(Tri2);
+               else
+                  Triangles.Append(Tri3);
+                  Triangles.Append(Tri4);
+               end if;
             end;
          end loop;
       end loop;

src/stl.adb

@@ -14,9 +14,6 @@ package body Stl is
       R.X := A.Y * B.Z - A.Z * B.Y;
       R.Y := A.Z * B.X - A.X * B.Z;
       R.Z := A.X * B.Y - A.Y * B.X;
-      if R.Z < 0.0 then
-         R.Z := -R.Z;
-      end if;
       return R;
    end Cross;
    

src/tesselada.adb

@@ -12,7 +12,7 @@ procedure Tesselada is
    Input : File_Type;
    Minimum_Thickness : constant := 1.0;
    Z_Scaling_Value : constant := 10.0;
-   Max_X_Y_Dimension : constant := 100.0;
+   Max_X_Y_Dimension : constant := 244.0;
 begin
    Ada.Text_IO.Put_Line("Tess started.");
    Open(Input, In_File, "sample-data/katie.dat");

tesselada.gpr

@@ -0,0 +1,5 @@
+project Tesselada is
+    for Source_Dirs use ("src");
+    for Object_Dir use "obj";
+    for Main use ("tesselada.adb");
+end Tesselada;