polytri.cpp

来自「算断裂的」· C++ 代码 · 共 865 行 · 第 1/2 页

// ------------------------------------------------------------------
// polytri.cpp
//
// This file contains a function for triangulating a planar polygon
// or PSLG using CDT.
// ------------------------------------------------------------------
// Author: Stephen A. Vavasis
// Copyright (c) 1999 by Cornell University.  All rights reserved.
// 
// See the accompanying file 'Copyright' for authorship information,
// the terms of the license governing this software, and disclaimers
// concerning this software.
// ------------------------------------------------------------------
// This file is part of the QMG software.  
// Version 2.0 of QMG, release date September 3, 1999.
// ------------------------------------------------------------------

#ifdef _MSC_VER
#if _MSC_VER == 1200
#include "qmatvec.h"
#endif
#endif

#include "qpolytri.h"
#include "qmatvec.h"

#include <fstream>


namespace QMG {
  namespace PolyTri {
    Triangulation
      pslg_constrained_delaunay_triangulate(VertexList& vtxlist,
      const vector<Edge>& elist);
    Triangulation poly_cdt(const Matrix& vtxlist, 
      const vector<Edge>& elist,
      double tol);
  }
}

namespace {
  using namespace QMG;
  using namespace QMG::PolyTri;


  // The triangulation algorithm maintains a front, which is a
  // heap of labeled distances.  The label is an int.

  struct Labeled_dist {
    Real dist;
    int label;
    Labeled_dist(Real d, int l) : dist(d), label(l) { }
    Labeled_dist(const Labeled_dist& o) : dist(o.dist), label(o.label) { }
    Labeled_dist() : dist(0), label(-1) { }
    bool operator==(const Labeled_dist& other) const;
    bool operator<(const Labeled_dist& other) const;
  };



  bool Labeled_dist::operator==(const Labeled_dist& other) const {
    if (dist == other.dist && label == other.label) return true;
    return false;
  }

  bool Labeled_dist::operator<(const Labeled_dist& other) const {
    if (dist < other.dist || (dist == other.dist && label < other.label))
      return true;
    return false;
  }
  Labeled_dist null_labeled_dist;

  // This routine projects three-dimensional vertices into 2D.
  // orthogonally to the plane containing the vertices.
  // The plane is computed via QR factorization with pivoting.


  void make_2d_vertices(const Matrix& vtxlist, VertexList& vtxlist2,
    double tol) {
    int m = vtxlist.numrows();
    int n = vtxlist.numcols();
    if (n < 2 || m < 3) {
      throw_error("Range error in make_2d_vertices");
    }
    if (n == 2) {
      for (int k = 0; k < m; ++k) {
        vtxlist2.add_vertex(vtxlist(k,0), vtxlist(k,1));
      }
      return;
    }
    Matrix A(n,m - 1);
    vector<Real> hhtrans(n);
    Real Afrob = 0.0;
    {
      for (int i = 1; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
          A(j, i - 1) = vtxlist(i,j) - vtxlist(0,j);
          Afrob += A(j,i-1) * A(j,i-1);
        }
      }

    }
    Afrob = sqrt(Afrob);

    {
      for (int l = 0; l < 2; ++l) {
        Real maxnorm = 0.0;
        int pivcol = 0;
        for (int p = 0; p < m - 1; ++p) {
          Real thiscolnorm = 0.0;
          for (int i = l; i < n; ++i) {
            thiscolnorm += A(i,p) * A(i,p);
          }
          if (thiscolnorm > maxnorm) {
            maxnorm = thiscolnorm;
            pivcol = p;
          }
        }
        Real beta = Matrix::compute_hh_transform(&hhtrans[l], 
          &A(l,pivcol), m - 1, n - l);


        {
          for (int j = 0; j < m - 1; ++j) {
            Real ip = 0.0;
            {
              for (int i = l; i < n; ++i) {
                ip += hhtrans[i] * A(i,j);
              }
            }
            Real beta1 = beta * ip;
            {
              for (int i = l; i < n; ++i) {
                A(i,j) -= beta1 * hhtrans[i];
              }
            }
          }
        }
      }
    }
    {
      Real resdfro = 0.0;
      for (int j = 0; j < m - 1; ++j) {
        for (int i = 2; i < n; ++i) {
          resdfro += A(i,j) * A(i,j);
        }
      }
      resdfro = sqrt(resdfro);
      if (resdfro >= tol * Afrob)
        throw_error("Polygon vertices exceed non-coplanarity tolerance");
    }
    vtxlist2.add_vertex(0.0, 0.0);
    {
      for (int j = 0; j < m - 1; ++j) {
        vtxlist2.add_vertex(A(0,j), A(1,j));
      }
    }
  }

  // These two structs are used by search tri for walking from
  // triangle to triangle.

  enum Triangle_Label {NOT_YET_SEARCHED, KEEP, DELETE};
  
  struct Two_Int {
    int first;
    int second;
    Two_Int() : first(-1), second(-1) { }
    Two_Int(const Two_Int& other) : first(other.first), second(other.second) { }
    Two_Int& operator=(const Two_Int& o) {
      first = o.first;
      second = o.second;
      return *this;
    }
    ~Two_Int() { }
  };

  
  // This routine takes as input a triangulation and labels
  // the triangles according to whether they are inside
  // or outside the polygon.  This is determined by
  // starting from the outer boundary (whose triangles
  // are not in the polygon) and then walking from
  // triangle to neighboring triangle.  Each time
  // a boundary is crossed, the disposition changes.


  void search_tri(const vector<Edge>& elist, 
    const Triangulation& init_tri, 
    vector<Triangle_Label>& disposition) {

    map<Edge,int> countmap;
    int nume = elist.size();

    {
      for (int k = 0; k < nume; ++k) {
        ++countmap[elist[k].sort()];
      }
    }


    int nt = init_tri.numtri();

    // Map edges to the two triangles that contain them.
    map<Edge,Two_Int> tri_occur;
    {
      for (int l = 0; l < nt; ++l) {
        for (int j = 0; j < 3; ++j) {
          Edge e(init_tri.tri_vertex(l,j), init_tri.tri_vertex(l, (j+1) % 3));
          e = e.sort();
          Two_Int twoint = tri_occur[e];
          if (twoint.first == -1)
            twoint.first = l;
          else if (twoint.second == -1)
            twoint.second = l;
          else
            throw_error("edge adjacent to >2 triangles");
          tri_occur[e] = twoint;
        }
      }
    }

    // Pick a starting triangle that was labeled by the caller.

    int starttri;
    for (starttri = 0; starttri < nt; ++starttri) {
      if (disposition[starttri] != NOT_YET_SEARCHED)
        break;
    }
    if (starttri == nt)
      throw_error("No starting triangle found");

    vector<int> search_stack;
    search_stack.push_back(starttri);

    // Now search all triangles with a stack.

    while (search_stack.size() > 0) {
      int curtri = search_stack.back();
      search_stack.pop_back();
      Triangle_Label thislab = disposition[curtri];
      if (thislab == NOT_YET_SEARCHED)
        throw_error("Stack inconsistency");
      for (int j = 0; j < 3; ++j) {
        Edge e(init_tri.tri_vertex(curtri, j), 
          init_tri.tri_vertex(curtri, (j+1) % 3));
        e = e.sort();
        Two_Int twoint = tri_occur[e];
        int nbrtri = -1;
        if (twoint.first == curtri)
          nbrtri = twoint.second;
        else if (twoint.second == curtri)
          nbrtri = twoint.first;
        if (nbrtri < 0)
          throw_error("Inconsistent map state");
        int occ_cnt = countmap[e];
        Triangle_Label newlab;
        if ((occ_cnt % 2) ^ (thislab == KEEP)) {
          newlab = KEEP;
        }
        else {
          newlab = DELETE;
        }
        if (disposition[nbrtri] == NOT_YET_SEARCHED) {
          disposition[nbrtri] = newlab;
          search_stack.push_back(nbrtri);
        }
        else {
          if (disposition[nbrtri] != newlab)
            throw_error("Illegal polygon -- boundary does not close");
        }
      }
    }

    {
      for (int l = 0; l < nt; ++l) {
        if (disposition[l] == NOT_YET_SEARCHED)
          throw_error("Disconnected triangulation?");
      }
    }
  }
}


namespace QMG {
  namespace PolyTri {
    ostream& operator<<(ostream& os, const Edge& e) {
      os << "Ed:v" << e.endpoint(0) << "/v" << e.endpoint(1) << '/';
      return os;
    }
  }
}


// #define POLYTRILOG

// ------------------------------------------------------------------
// This routine computes the constrained Delaunay triangulation (CDT)
// of a planar straight line graph (PSLG).

QMG::PolyTri::Triangulation 
QMG::PolyTri::pslg_constrained_delaunay_triangulate(VertexList& vtxlist,
                                                    const vector<Edge>& elist) {


  // Make a list of all constraint edges.
  // Use a map to prevent duplicates.

  vector<Edge> constraint;
  vector<bool> top_done;
  vector<bool> bottom_done;
  map<Edge,int> checkedge;

  using std::advance;



#ifdef POLYTRILOG

  using std::endl;
  std::ofstream dbstream("pslglog");
  ostream* debug_str = &dbstream;
  for (int nv1 = 0; nv1 < vtxlist.size(); ++nv1) {
    *debug_str << "vertex v" << nv1 << "/ at " << vtxlist.coord(nv1,0) 
      << ',' << vtxlist.coord(nv1,1) << endl;
  }
#endif


  {
    int numv = vtxlist.size();
    for (int j = 0; j < elist.size(); ++j) {
      Edge e = elist[j];
      for (int k = 0; k < 2; ++k) {
      if (e.endpoint(k) < 0 || e.endpoint(k) >= numv)
        throw_error("Entry in edge list out of range");
      }
      if (e.endpoint(0) == e.endpoint(1)) continue;
      if (checkedge.find(e.sort()) == checkedge.end()) {

        int edgeind = constraint.size();
        constraint.push_back(e);
        top_done.push_back(false);
        bottom_done.push_back(false);
        checkedge[e.sort()] = edgeind;
#ifdef POLYTRILOG 
        *debug_str << " adding edge " << e << " at pos " << edgeind << endl;
#endif
      }
      else {
#ifdef POLYTRILOG
        *debug_str << " skipping edge " << e << endl;
#endif
      }
    }
  }

  // Add the bounding box to the stack, and four new vertices.

  {
    Real minx = BIG_REAL;
    Real maxx = -BIG_REAL;
    Real miny = BIG_REAL;
    Real maxy = -BIG_REAL;
    for (int j = 0; j < vtxlist.size(); ++j) {
      Real x = vtxlist.coord(j,0);
      if (x < minx)
        minx = x;
      if (x > maxx)
        maxx = x;
      Real y = vtxlist.coord(j,1);
      if (y < miny)
        miny = y;
      if (y > maxy)
        maxy = y;
    }

    if (vtxlist.size() == 0) {
      minx = 0;
      maxx = 1;
      miny = 0;
      maxy = 1;
    }

    Real deltax = maxx - minx;
    if (deltax == 0.0)
      deltax = fabs(minx) + 1;
    Real deltay = maxy - miny;
    if (deltay == 0.0)
      deltay = fabs(miny) + 1;
    Real lbx = minx - deltax;
    Real ubx = maxx + deltax;
    Real lby = miny - deltay;
    Real uby = maxy + deltay;
    VertexIndex newv1 = vtxlist.add_vertex(lbx,lby);
    VertexIndex newv2 = vtxlist.add_vertex(ubx,lby);
    VertexIndex newv3 = vtxlist.add_vertex(ubx,uby);
    VertexIndex newv4 = vtxlist.add_vertex(lbx,uby);
#ifdef POLYTRILOG
    *debug_str << " lbx, ubx, lby, uby = " << lbx << ',' << ubx << ',' 
      << lby << ',' <<  uby << endl;
    *debug_str << " new verts = v" << newv1 << "/ v" << newv2 << "/ v"
      << newv3 << "/ v" << newv4 << "/" << endl;
#endif
    int sz = constraint.size();
    constraint.push_back(Edge(newv1,newv2));
    top_done.push_back(false);
    bottom_done.push_back(true);
    checkedge[Edge(newv1,newv2).sort()] = sz++;
    constraint.push_back(Edge(newv2,newv3));
    top_done.push_back(false);
    bottom_done.push_back(true);
    checkedge[Edge(newv2,newv3).sort()] = sz++;
    constraint.push_back(Edge(newv3,newv4));
    top_done.push_back(false);
    bottom_done.push_back(true);
    checkedge[Edge(newv3,newv4).sort()] = sz++;
    constraint.push_back(Edge(newv4,newv1));
    top_done.push_back(false);
    bottom_done.push_back(true);
    checkedge[Edge(newv4,newv1).sort()] = sz++;
  }

  Matrix transf(2,2);
  Matrix smalleq(2,2);
  Real smallrhs[2];
  Real smallsoln[2];
  vector<Labeled_dist> vtxqueue;
  vector<Labeled_dist> edgequeue;
  VertexList tvtxlist;
  vector<Real> tan_point_x;
  vector<Real> tan_point_y;
  vector<Real> slope_c;