title: Delaunay Triangulation in C++
body: <p>The following steps assume you start with two points – A &amp; B – and are trying to determine the point C to be used to form a
  triangle:</p>
  
  <p>a. Create a member function that will determine if a given point C is left or right of the line formed by two points A and B. Hint: to do this, take the cross product of the vector
  between the two points A and B and the vector between A and C. Since the cross product
  is proportional to the sine of the angle between two vectors, it will be a positive value for an angle between 0 and 180 degrees (i.e., if point C is to the left of the line from A to B).</p>
  
  <p>b. Create a member function that will determine if a given point is inside the circle formed by three other points (in this case, the three points will be the points of a triangle and the resulting circle will be the circumcircle). Hint: a very elegant implementation of this
  function can be found on the Wikipedia entry for Delaunay Triangulation. If you borrow
  this implementation, make sure to test it thoroughly, reference it, and explain how it
  works in your lab report.</p>
  
  <p>c. Create a member function that, given two points, finds the point for the next Delaunay
  triangle. This function can search the whole list (not the most efficient implementation,
  but sufficient for this lab) and most likely will call the functions defined in 4a and 4b
  above.</p>
  
  <p>d. Create a recursive member function, Delaunay(Point A, Point B), that calls the function
  from 4c above the find the next point C. When point C is found, this function should
  insert a new triangle into the triangle vector and also update the Boolean variable in the
  point list. This function should then recursively call itself using points A and C, and again using points C and B. Make sure the function also has two base cases; one for when the point C found has already been used, in which case the triangle should be added, but the recursive call should not occur, and; another case where no point C is found (because A
  and B form an outside edge of the dataset)</p>
  
  <p>I have completed the code but seem to keep coming up with error in the debugging process and I believe I've narrowed it down to the read function so here is my code:</p>
  
  <pre><code>PointList::PointList()
  {
  totPt = 0;
  totTri = 0;
  }
  
  
  void PointList::read(const char* filename)
  {
  FILE* file;
  file = fopen ( filename, "r" );
  
  fscanf ( file, "%d  \n", &amp;totPt);
  datapoint.resize( totPt );
  
  for ( unsigned int i=0; i &lt; datapoint.size(); i++ )
  {
      fscanf ( file, " %d %lf %lf \n", &amp;datapoint.at(i).ptnum, &amp;datapoint.at(i).xCoord, &amp;datapoint.at(i).yCoord  );
  }
  
  }
  
  void PointList::TriPrint()
  {
  cout &lt;&lt; "# of triangles created: " &lt;&lt; triPoint.size() &lt;&lt; endl;
  cout &lt;&lt; "Triangle # : Vertice1   Vertice2  Vertice3" &lt;&lt; endl;
  
  for ( unsigned int i = 0; i &lt; triPoint.size() ; i++)
  {
      cout &lt;&lt; triPoint.at(i).triNum &lt;&lt; " : " &lt;&lt; triPoint.at(i).vert1.ptnum &lt;&lt; "   " &lt;&lt; triPoint.at(i).vert2.ptnum &lt;&lt; "   " &lt;&lt; triPoint.at(i).vert3.ptnum &lt;&lt; endl;
  }
  return;
   }
  </code></pre>
  
  <p>//TASK 4a: will determine if a given point C is left or right of the line forme by the two points a and b</p>
  
  <pre><code>bool PointList::isRight( double a, double b, double c )
  {
   Point A = datapoint.at(a-1);
  Point B = datapoint.at(b-1);
  Point C = datapoint.at(c-1);
  
  //Cross product of AB and AC
  Point AB;
  Point AC;
  
   AB.xCoord = B.xCoord - A.xCoord;
  AB.yCoord = B.yCoord - A.yCoord;
  
   AC.xCoord = C.xCoord - A.xCoord;
  AC.yCoord = C.yCoord - A.yCoord;
  
   double det = (AB.xCoord*AC.yCoord) - (AC.xCoord*AB.yCoord);
  
  if ( det &gt;= 0 )
  {  
      return true;
  }
  else
  {
      return false;
  }
  }
  </code></pre>
  
  <p>//Task 4b: determine if a given point is inside circle formed by three other points</p>
  
  <pre><code>bool PointList::inCircle ( double a, double b, double c, double d )
  {
  bool inCirc = false;
  Point A = datapoint.at(a-1);
  Point B = datapoint.at(b-1);
  Point C = datapoint.at(c-1);
  Point D = datapoint.at(d-1);
  
  vector &lt;vector&lt;double&gt;&gt; Matrix;
  Matrix.resize(3);
  for ( int i = 0; i &lt; 3 ; i++)
  {
      Matrix.at(i).resize(3);
  }
  
  Matrix.at(0).at(0) = A.xCoord - D.xCoord;
  Matrix.at(1).at(0) = B.xCoord - D.xCoord;
  Matrix.at(2).at(0) = C.xCoord - D.xCoord;
  
  Matrix.at(0).at(1) = A.yCoord - D.yCoord;
  Matrix.at(1).at(1) = B.yCoord - D.yCoord;
  Matrix.at(2).at(1) = C.yCoord - D.yCoord;
  
  Matrix.at(0).at(2) = ((A.xCoord*A.xCoord) - (D.xCoord*D.xCoord)) + ((A.yCoord*A.yCoord) - (D.yCoord*D.yCoord));
  Matrix.at(1).at(2) = ((B.xCoord*B.xCoord) - (D.xCoord*D.xCoord)) + ((B.yCoord*B.yCoord) - (D.yCoord*D.yCoord));
  Matrix.at(2).at(2) = ((C.xCoord*C.xCoord) - (D.xCoord*D.xCoord)) + ((C.yCoord*C.yCoord) - (D.yCoord*D.yCoord));
  
  
  double det = 0;
  
  det = (Matrix.at(0).at(0) * Matrix.at(1).at(1) * Matrix.at(2).at(2)) + (Matrix.at(0).at(1) * Matrix.at(1).at(2) * Matrix.at(2).at(0)) + (Matrix.at(0).at(2) * Matrix.at(1).at(0) * Matrix.at(2).at(1));
  det = det - (Matrix.at(2).at(0) * Matrix.at(1).at(1) * Matrix.at(0).at(2)) - (Matrix.at(2).at(1) * Matrix.at(1).at(2) * Matrix.at(0).at(0)) - (Matrix.at(2).at(2) * Matrix.at(1).at(0) * Matrix.at(0).at(1));
  
  if ( det &gt;= 0 )     // determinant is positive if and only if D lies inside the circumcircle
  {
      inCirc = false;
  }
  else
  {
      inCirc = true;
  }
  
  return inCirc;
  }
  </code></pre>
  
  <p>//Task 4c: Given two points, finds the point for the next Delaunay triangle</p>
  
  <pre><code>double PointList::nextDelaunay ( double a, double b )
  {
  bool incircle = false;
  bool isright = false;
   bool noRight = false;
   int f = -1;
  
  
  //check all points in file
  for ( int i = 0; i &lt; datapoint.size() ; i++)
  {
      if ( i == a || i == b)
      {
          continue;
      }
       else
      {
           isright = isRight( a, b, i);   //checks to see if its to the right
          if(isright)
          {
              noRight = false;
              for ( int j = 1; j &lt; datapoint.size(); j++ )
              {
                  if ( j == a || j == b || j == i )
                                     {
                      continue;
                  }
                  else
                  {
                      incircle = inCircle ( a, b, i, j );    
                                             //checks to see if point in circle
                       if(incircle)
                      {
                          break;
                      }
                  }
              }
              if ( !incircle)
              {
                  return i;
              }
  
          }
          else
          {
              continue;
          }
      }
  }
  if (noRight)
  {
      return f;
   }
  
  }
  </code></pre>
  
  <p>//Task 4d: a recursive member function, Delaunay, the calls from 4c above to find the next point C</p>
  
  <pre><code>void PointList::Delaunay( int a, int b )
  {
  Triangle x;
  Point A = datapoint.at(a - 1);
  Point B = datapoint.at(b - 1);
  
  int c = nextDelaunay(a,b);
  
  if ( c == -1)
  {
      return;
  }
  else
  {
      Point C = datapoint.at(c-1);
  
  
      if ( C.usedPt == false)
      {
          x.vert1.ptnum = a;
          x.vert1.xCoord = A.xCoord;
          x.vert1.yCoord = A.yCoord;
  
          x.vert2.ptnum = b;
          x.vert2.xCoord = B.xCoord;
          x.vert2.yCoord = B.yCoord;
  
          x.vert3.ptnum = c;
          x.vert3.xCoord = C.xCoord;
          x.vert3.yCoord = C.yCoord;
  
          x.triNum = triPoint.size()+1;
          triPoint.push_back(x);
          datapoint.at(c-1).usedPt = true;
  
          Delaunay( a, c );
          Delaunay( c, b );
  
      return;
      }
  
      if ( C.usedPt == true )
      {
          x.vert1.ptnum = a;
          x.vert1.xCoord = A.xCoord;
          x.vert1.yCoord = A.yCoord;
  
          x.vert2.ptnum = b;
          x.vert2.xCoord = B.xCoord;
          x.vert2.yCoord = B.yCoord;
  
          x.vert3.ptnum = c;
          x.vert3.xCoord = C.xCoord;
          x.vert3.yCoord = C.yCoord;
  
          x.triNum = triPoint.size()+1;
          triPoint.push_back(x);
  
      return;
      }
  }
  }
  </code></pre>