DelaunayMeshGenerator.cpp

//
//  file    : DelaunayMeshGenerator.cpp
//  author  : jh
//  date    : 15 jul 2002
//  changed : 29 jul 2002
//

#include <stdlib.h>

#include "DelaunayMeshGenerator.h"
#include "geometry.h"

#include <qdatetime.h>


extern char* __delaunay_lib_version;



DelaunayMeshGenerator::DelaunayMeshGenerator(const QString& filename, bool verbose)
  : QXmlDefaultHandler(), _verbose(verbose)
{
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::DelaunayMeshGenerator()>\n");
#endif

  reset();
  QFile xmlfile(filename);
  open(xmlfile);


#ifdef _DEBUG_ 
  fprintf(stderr, "numer of vertices : %d \n", _vertices.count() );
  fprintf(stderr, "</DelaunayMeshGenerator::DelaunayMeshGenerator()>\n");
#endif

}



DelaunayMeshGenerator::DelaunayMeshGenerator(QIODevice& iodevice, bool verbose)
  : QXmlDefaultHandler(), _verbose(verbose)
{
  reset();
  open(iodevice);
}




void
DelaunayMeshGenerator::reset()
{
  _mesh     = 0;
  _count_ID = 1;
}



DelaunayMeshGenerator::~DelaunayMeshGenerator()
{
  // _vertices.setAutoDelete(true);
  // _active_vertices.setAutoDelete(true);
  
  _edges.setAutoDelete(true);
  _active_edges.setAutoDelete(true);
  
  // _triangles.setAutoDelete(true);
  
  if(_mesh != 0)
    delete _mesh;
}



void
DelaunayMeshGenerator::add_vertex(Vertex* v) 
{
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::addVertex()>\n");
#endif

  _vertices.append(v); 
  
#ifdef _DEBUG_
  fprintf(stderr, "</DelaunayMeshGenerator::addVertex()>\n");
#endif

}



void
DelaunayMeshGenerator::generate()
{
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::generate()>\n");
#endif

  init();

#ifdef _DEBUG_
  fprintf(stderr, "[DelaunayMeshGenerator::generate()] bbox (%f/%f)-(%f/%f) \n", 
          minX, minY, maxX, maxY );
#endif

  _active_edges    = _edges;
  _active_vertices = _vertices;

  Vertex* new_vertex = 0;


  int count = 0;


  while(!_active_edges.isEmpty()) {
    
    Edge* edge = _active_edges.take(0);

    /*
      fprintf(stderr, "<DelaunayMeshGenerator::generate()> [%d] next active edge [%s]-[%s] ... \n", 
      count++,
      edge->get_vertex(0)->get_ID().latin1(),
      edge->get_vertex(1)->get_ID().latin1() );
    */

    if(_verbose) {
      count++;
      if(count % 100 == 0) {
        fprintf(stderr, "number of active vertices = %d \n", _active_vertices.count() );
      }
    }

    ((Vertex*) edge->vertex(0))->grad -= 1;
    ((Vertex*) edge->vertex(1))->grad -= 1;

    box_minX = minX;
    box_maxX = maxX;
    box_minY = minY;
    box_maxY = maxY;

    new_vertex = 0;
    u_alt      = 1e15;


    QPtrListIterator<Vertex> it(_active_vertices);

    Vertex* vertex;

    while( (vertex = it.current()) != 0 ) {
      
      ++it;

#ifdef _DEBUG_
      fprintf(stderr, "check vertex [%s] (%d active vertices) \n", 
              vertex->get_ID().latin1(), _active_vertices.count() );
#endif
    

      if(vertex == edge->vertex(0) || vertex == edge->vertex(1) ) {
#ifdef _DEBUG_
        fprintf(stderr, "continue [%s] (1) \n", vertex->get_ID().latin1() );
#endif
        continue;     
      }


      if(vertex->x() < box_minX || vertex->x() > box_maxX ||
         vertex->y() < box_minY || vertex->y() > box_maxY  ) {

#ifdef _DEBUG_
        fprintf(stderr, "continue [%s] (2) \n", vertex->ID().latin1() );
#endif
        
        continue;   

      }

      
      if( geometry::area(edge->vertex(0),
                         edge->vertex(1),
                         vertex               ) <= CEps ) {
#ifdef _DEBUG_
        fprintf(stderr, "continue [%s] (3) \n", vertex->ID().latin1() );
#endif
        continue;
      }



      if( test_intersection(vertex, edge->vertex(0)) ) {
#ifdef _DEBUG_
        fprintf(stderr, "continue [%s] (4) \n", vertex->ID().latin1() );
#endif
        continue;
      }

        

      if( test_intersection(edge->vertex(1), vertex) ) {

#ifdef _DEBUG_
        fprintf(stderr, "continue [%s] (5) \n", vertex->ID().latin1() );
#endif
        continue;            

      }
    

      if( is_voronoi_triangle(edge, vertex) )
        new_vertex = vertex;


    }
    
    if(new_vertex != 0) {
      create_triangle(edge, new_vertex);
    } else {
      fprintf(stderr, "ERROR : newpoint = null! \n");
      exit(1);
    }

  }

  /*
    const QDict<Vertex>& v_list = _mesh->get_vertices();    
    fprintf(stderr, "</generate()> numer of vertices (mesh) : %d \n", v_list.count() );
  */

#ifdef _DEBUG_
  fprintf(stderr, "</DelaunayMeshGenerator::generate()>\n");
#endif

}



void
DelaunayMeshGenerator::create_triangle(Edge* basisedge, Vertex* vertex)
{
#ifdef _DEBUG_
  fprintf(stderr, "</DelaunayMeshGenerator::createTriangle()>\n");
  fprintf(stderr, "[/DelaunayMeshGenerator::createTriangle()] basisedge [%s]-[%s], vertex [%s]\n",
          basisedge->vertex(0)->ID().latin1(), 
          basisedge->vertex(1)->ID().latin1(),
          vertex->ID().latin1() );
#endif


  Triangle* tri;
  Edge* edge;
  Edge* e[2];
   
  tri = new Triangle(createID(),
                     basisedge->vertex(0),
                     basisedge->vertex(1),
                     vertex);

  vertex->used = true;


  //  check if new edge already exists

  e[0] = e[1] = 0;


  QPtrListIterator<Edge> it(_active_edges);

  while( (edge = it.current()) != 0 ) {

    ++it;

    for(int i=0; i<2; i++) {

      if(vertex == edge->vertex(i) &&
         basisedge->vertex(i) == edge->vertex((i+1)%2)) {
          
        e[i] = edge;
          
        ((Vertex*) edge->vertex(0))->grad -= 1;
        ((Vertex*) edge->vertex(1))->grad -= 1;   

        if(((Vertex*) edge->vertex(0))->grad == 0 && 
           ((Vertex*) edge->vertex(0))->used == true )
          _active_vertices.removeRef( ((Vertex*) edge->vertex(0)) );
               
        if(((Vertex*) edge->vertex(1))->grad == 0 && 
           ((Vertex*) edge->vertex(1))->used == true )
          _active_vertices.removeRef( ((Vertex*) edge->vertex(1)) ); 
        
      }

      if(e[0] != 0 && e[1] != 0)
        break;

    }

  }  // while(iter.hasNext()) 

          
  if(e[0] == 0) {
    
#ifdef _DEBUG_
    fprintf(stderr, "[DelaunayMeshGenerator::createTriangle()] create new edge \n");
    fprintf(stderr, " new edge [%s][%s] \n", 
            basisedge->vertex(0)->ID().latin1(), 
            vertex->ID().latin1() );
#endif

    e[0] = new Edge(basisedge->vertex(0), vertex);

    vertex->grad += 1;
    ((Vertex*) basisedge->vertex(0))->grad += 1;

    _active_edges.append( e[0] );
    
  } else {
    _active_edges.removeRef( e[0] );
  }
  
  
  if(e[1] == 0) {    

#ifdef _DEBUG_
    fprintf(stderr, "[DelaunayMeshGenerator::createTriangle()] create new edge \n");
    fprintf(stderr, " new edge [%s][%s] \n", 
            vertex->ID().latin1(),
            basisedge->vertex(1)->ID().latin1() );
#endif
    
    e[1] = new Edge(vertex, basisedge->vertex(1));

    vertex->grad += 1;
    ((Vertex*) basisedge->vertex(1))->grad += 1; 

    _active_edges.append( e[1] );
    
  } else {
    _active_edges.removeRef( e[1] );
  }
  
    
  /*
   *  finally add new triangle to array
   */
  _mesh->add_triangle(tri);

#ifdef _DEBUG_
  fprintf(stderr, "</DelaunayMeshGenerator::createTriangle(): triangle added.\n");
#endif


}



bool
DelaunayMeshGenerator::test_intersection(Vertex* v1, Vertex* v2)
{
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::testIntersection(%s - %s)\n",
          v1->ID().latin1(), v2->ID().latin1() );
  fprintf(stderr, "          (%f/%f)-(%f/%f))>\n",
          v1->x(), v1->y(), v2->x(), v2->y() );
#endif

  /*
    QPtrListIterator<Edge> it(_active_edges );
  */

  Edge* edge;

  geometry::line l1, l2;
  
  l1.p1.x = v1->x();
  l1.p1.y = v1->y();
    
  l1.p2.x = v2->x();
  l1.p2.y = v2->y();


  // while( (edge = it.current()) != 0 ) {

  int size = _active_edges.count();

#ifdef _DEBUG_
  fprintf(stderr, "[DelaunayMeshGenerator::testIntersection()] %d active edges\n",
          size );
#endif

  for(int i=0; i<size; i++) {

    edge = _active_edges.at(i);

#ifdef _DEBUG_
    fprintf(stderr, "[DelaunayMeshGenerator::testIntersection()] test [%d/%d] active edge [%s]-[%s]\n",
            i, size, 
            edge->vertex(0)->ID().latin1(),
            edge->vertex(1)->ID().latin1() );
    fprintf(stderr, "          edge [%s - %s] \n",
            edge->vertex(0)->ID().latin1(), 
            edge->vertex(1)->ID().latin1() );
#endif
    
    if( edge->vertex(0) == v1 || edge->vertex(1) == v1 ) {
      continue;
    }

    if( edge->vertex(0) == v2 || edge->vertex(1) == v2 ) {
      continue;
    }

    l2.p1.x = edge->vertex(0)->x();
    l2.p1.y = edge->vertex(0)->y();
    l2.p2.x = edge->vertex(1)->x();
    l2.p2.y = edge->vertex(1)->y();

    if(geometry::intersect(l1,l2)) {

#ifdef _DEBUG_
      fprintf(stderr, "<DelaunayMeshGenerator::testIntersection()> return <true>\n");
#endif

      return true;

    }

  }
  
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::testIntersection()> return <false>\n");
#endif

  return false;

}



bool
DelaunayMeshGenerator::is_voronoi_triangle(Edge* edge, Vertex* v)
{
  double u;
      
  double e2[2];       
    
  double xm,ym;       
  double dx,dy;       
  double radius;      


  double x1 = edge->vertex(0)->x();
  double y1 = edge->vertex(0)->y();
    
  double x2 = edge->vertex(1)->x();
  double y2 = edge->vertex(1)->y();
    
  double x3 = v->x();
  double y3 = v->y();
    

  e2[0] = -(y2 - y1) / sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1));
  e2[1] =  (x2 - x1) / sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1));
   
   
  u = ((x3 - x2)*(x3 - x1) + (y3 - y2)*(y3 - y1)) /
    (2.0 * (e2[0] * (x3 - x1) +   e2[1] * (y3 - y1)));
  
      
  if(u_alt > u) {
    
    u_alt = u;
    
    xm        = 0.5 * (x1 + x2) + u_alt * e2[0];
    ym        = 0.5 * (y1 + y2) + u_alt * e2[1];
    dx        = xm - x1;
    dy        = ym - y1;
    radius    = sqrt(dx*dx+dy*dy);
    
    box_minX = xm - radius;   
    box_minY = ym - radius;
    box_maxX = xm + radius;
    box_maxY = ym + radius;
    
    return true;
    
  }
  
  return false;

}



void
DelaunayMeshGenerator::init() 
{
  _mesh = new Mesh();

  minX = minY =  1e15;
  maxX = maxY = -1e15;

  for(Vertex* v = _vertices.first(); v; v = _vertices.next()) {

    _mesh->add_vertex(v);

    v->grad = 0;
    v->used = false;
    
    if(minX > v->x()) minX = v->x();
    if(maxX < v->x()) maxX = v->x();
    if(minY > v->y()) minY = v->y();
    if(maxY < v->y()) maxY = v->y();

  }


  /*
    const QDict<Vertex>& v_list = _mesh->get_vertices();  
    fprintf(stderr, "numer of vertices (mesh) : %d \n", v_list.count() );
  */


  //  set state of all edges to contour

  for(Edge* e = _edges.first(); e; e = _edges.next()) {

    for(int i=0; i<2; i++) {
      Vertex* v = (Vertex*) e->vertex(i);
      v->grad += 1;
      v->used = true;
    }

  }

}



void
DelaunayMeshGenerator::write_XML(QTextStream& cout)
{
  cout << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>\n";
  cout << "<!-- created by delaunay mesh generator \n";
  cout << "     " << __delaunay_lib_version << endl;
  cout << "     " << QTime::currentTime().toString() << " -->\n";

  cout << "<objects>\n";

  for(Vertex* v = _vertices.first(); v; v = _vertices.next()) {
    v->write_XML(cout);
  }

  for(Edge* e = _edges.first(); e; e = _edges.next()) {
    e->write_XML(cout);
  }

  cout << "</objects>\n";
}



bool
DelaunayMeshGenerator::startElement(const QString&, const QString&, 
                                    const QString& qName, 
                                    const QXmlAttributes& xmlattr)
{
  
  QString qNameL = qName.stripWhiteSpace().lower();  
  
#ifdef _DEBUG_
  fprintf(stderr, "<DelaunayMeshGenerator::startElement()> elem = %s \n", 
          qNameL.latin1() );
#endif
  
  if(qNameL == "point") {
    
    QString id = xmlattr.value("id");
    QString x  = xmlattr.value("x");
    QString y  = xmlattr.value("y");

    Vertex* v = new Vertex(id, x.toDouble(), y.toDouble());

    // fprintf(stderr, "created vertex : [%s] \n", v->ID().latin1() );

    add_vertex(v);

    return true;
    
  }   

  if(qNameL == "edge") {
    
    QString from  = xmlattr.value("from");
    QString to    = xmlattr.value("to");
    
    Vertex* v1 = vertex(from);
    Vertex* v2 = vertex(to);
    
    if(v1 == 0) {
      fprintf(stderr, "<DelaunayMeshGenerator::startElement()> : ERROR = \n");
      fprintf(stderr, "  could not find point [%s] !!\n", from.latin1());
      exit(1);
    }

    if(v2 == 0) {
      fprintf(stderr, "<DelaunayMeshGenerator::startElement()> : ERROR = \n");
      fprintf(stderr, "  could not find point [%s] !!\n", to.latin1());
      exit(1);
    }
    
    Edge* e = new Edge(v1,v2);

    add_edge(e);
    
    return true;

  }
  
  return true;

}

---