Delaunay triangulation

Finally I have a little time to implement my project. I have implemented Graphics Device Abstraction Layer, XNA Implementations for WinForms and WPF. Than I have implementing Mathematics Subsystem and now I have reached Topology Subsystem.
There are two kind of Delaunay triangulation. This is 3D and 2D triangulation. In my case I need to tessellate a 3D polygon which locates on the plane. Well I decided implement 2d case for a start.
Let’s see what has been done. The essence of the algorithm is adding step by step points and edges. An important thing, adding new point we have to broke some edges to push new point in optimal way. Let’s see the code.

float xmin = points[0].X;
float ymin = points[0].Y;
float xmax = xmin;
float ymax = ymin;
for (int i = 1; i < points.Length; i++)
{
     if (points[i].X < xmin) xmin = points[i].X;
     if (points[i].X > xmax) xmax = points[i].X;
     if (points[i].Y < ymin) ymin = points[i].Y;
     if (points[i].Y > ymax) ymax = points[i].Y;
}

float dx = xmax - xmin;
float dy = ymax - ymin;
float dmax = (dx > dy) ? dx : dy;

float xmid = (xmax + xmin) * 0.5f;
float ymid = (ymax + ymin) * 0.5f;

In this fragment we calculate the maximum and minimum coordinate. This will be useful for calculating the bounding triangle. Bounding triangle (or supertriangle) includes all input points. Now we can construct the supertriangle and add to the list.

Vector2[] vertices = new Vector2[points.Length + 3];
points.CopyTo(vertices, 0);

vertices[points.Length + 0] = new Vector2(
     xmid - 2.0f * dmax, ymid - dmax);
vertices[points.Length + 1] = new Vector2(
     xmid, ymid + 2.0f * dmax);
vertices[points.Length + 2] = new Vector2(
     xmid + 2 * dmax, ymid - dmax);

List<IndexedTriangle> triangles =
     new List<IndexedTriangle>();
triangles.Add(new IndexedTriangle(
     points.Length,
     points.Length + 1,
     points.Length + 2));

Now we have all points in the "vertices" and the supertriangle in the "triangles". The algorithm begins here...

for (int i = 0; i < points.Length; i++)
{
   List<IndexedEdge> edges = new List<IndexedEdge>();

   for (int j = 0; j < triangles.Count; j++)
   {
      Circle2 circumcircle = new Circle2();
      try
      {
            circumcircle = new Circle2(
            vertices[triangles[j].First],
            vertices[triangles[j].Second],
            vertices[triangles[j].Third]);
      }
      catch
      {
         // Points are collinear
         continue;
      }

      if (circumcircle.IsInside(vertices[i]))
      {
         edges.Add(new IndexedEdge(
            triangles[j].First, triangles   [j].Second));
         edges.Add(new IndexedEdge(
            triangles[j].Second, triangles[j].Third));
         edges.Add(new IndexedEdge(
            triangles[j].Third, triangles[j].First));
         triangles.RemoveAt(j);
         j--;
      }
   }
   // Clear dublicates
   for (int j = edges.Count - 2; j >= 0; j--)
   {
      for (int k = edges.Count - 1; k >= j + 1; k--)
      {
         if (edges[j] == edges[k])
         {
            edges.RemoveAt(k);
            edges.RemoveAt(j);
            k--;
            continue;
         }
      }
   }

   // Form triangles
   for (int j = 0; j < edges.Count; j++)
   {
      triangles.Add(MakeClockwise(
      new IndexedTriangle(
         edges[j].First,
         edges[j].Second, i), vertices));
   }
}

I have written a lot of modul tests and got "green bar". If you need the sources - ask me daVinci@mail.ru