polygon_mesher.C

/*
 * This file is part of the "Archon" framework.
 * (http://files3d.sourceforge.net)
 *
 * Copyright © 2002 by Kristian Spangsege and Brian Kristiansen.
 *
 * Permission to use, copy, modify, and distribute this software and
 * its documentation under the terms of the GNU General Public License is
 * hereby granted. No representations are made about the suitability of
 * this software for any purpose. It is provided "as is" without express
 * or implied warranty. See the GNU General Public License
 * (http://www.gnu.org/copyleft/gpl.html) for more details.
 *
 * The characters in this file are ISO8859-1 encoded.
 *
 * The documentation in this file is in "Doxygen" style
 * (http://www.doxygen.org).
 */

#include <cmath>
#include <limits>
#include <algorithm>

#include <archon/math/matrix.H>
#include <archon/math/functions.H>

#include <archon/util/polygon_mesher.H>

namespace Archon
{
  namespace Utilities
  {
    using namespace std;
    using namespace Math;

    namespace
    {
      struct PolygonException: MesherException
      {
        PolygonException(string l, string m): Utilities::Exception(l, m) {}
      };

      struct Vertex
      {
        Vector2 p;
        int i;
        Vertex() {}
        Vertex(const Vector2 &p, int i): p(p), i(i) {}
      };

      void splitX(const vector<Vertex> &polygon, double clip,
                  const Matrix3x3 &m, MesherContext *context,
                  vector<Vertex> &below, vector<Vertex> &above,
                  Vertex &reuseVertex, bool &reuse, double z)
      {
        Vertex rv = reuseVertex, v = polygon[polygon.size()-1];
        bool r = reuse;
        reuse = false;
        for(unsigned i=0; i<polygon.size(); ++i)
        {
          Vertex w = polygon[i];
          if(w.p[0] >= clip)
          {
            if(v.p[0] >= clip) above.push_back(w);
            else
            {
              if(!r)
              {
                Vector2 s(clip, linInterp(clip, v.p[0], w.p[0], v.p[1], w.p[1]));
                rv = Vertex(s, context->addVertex(m * Vector3(s[0], s[1], z)));
              }
              above.push_back(rv);
              above.push_back(w);
              below.push_back(rv);
            }
          }
          else
          {
            if(v.p[0] >= clip)
            {
              Vector2 s(clip, linInterp(clip, v.p[0], w.p[0], v.p[1], w.p[1]));
              Vertex t = Vertex(s, context->addVertex(m * Vector3(s[0], s[1], z)));
              below.push_back(t);
              below.push_back(w);
              above.push_back(t);
              reuseVertex = t;
              reuse = true;
            }
            else below.push_back(w);
          }
          v = w;
        }
      }

      void splitY(const vector<Vertex> &polygon, double clip,
                  const Matrix3x3 &m, MesherContext *context,
                  vector<Vertex> &below, vector<Vertex> &above,
                  double z)
      {
        Vertex v = polygon[polygon.size()-1];
        for(unsigned i=0; i<polygon.size(); ++i)
        {
          Vertex w = polygon[i];
          if(w.p[1] >= clip)
          {
            if(v.p[1] >= clip) above.push_back(w);
            else
            {
              Vector2 s(linInterp(clip, v.p[1], w.p[1], v.p[0], w.p[0]), clip);
              Vertex t = Vertex(s, context->addVertex(m * Vector3(s[0], s[1], z)));
              above.push_back(t);
              above.push_back(w);
              below.push_back(t);
            }
          }
          else
          {
            if(v.p[1] >= clip)
            {
              Vector2 s(linInterp(clip, v.p[1], w.p[1], v.p[0], w.p[0]), clip);
              Vertex t = Vertex(s, context->addVertex(m * Vector3(s[0], s[1], z)));
              below.push_back(t);
              below.push_back(w);
              above.push_back(t);
            }
            else below.push_back(w);
          }
          v = w;
        }
      }
    }

    void meshPolygon(const vector<Vector3> &polygon,
                     double maxPatchSize,
                     MesherContext *context)
    {
      // Determine a basis for the polygon
      Vector3 basisu = polygon[1] - polygon[0];
      Vector3 basisv = polygon[polygon.size()-1] - polygon[0];

      basisu.normalize();
      basisv.normalize();

      Vector3 basisw = basisu * basisv; // Normal
      basisv = basisw * basisu;

      basisv.normalize();
      basisw.normalize();

      Matrix3x3 mat(basisu, basisv, basisw);

      vector<Vertex> p1, p2, p3, p4;
      vector<Vertex> *poly = &p1, *stripe = &p2, *patch = &p3, *temp = &p4;

      // Transform to 2D and find AABB
      double
        minX = numeric_limits<double>::max(),
        minY = numeric_limits<double>::max(),
        maxX = numeric_limits<double>::min(),
        maxY = numeric_limits<double>::min();

      double z;    // Common z coordinate for all 2d-transformed vertices.
      Vector3 v;
      for(unsigned i=0; i<polygon.size(); ++i)
      {
        v = polygon[i] * mat;
        Vector2 w(v);
        if(i == 0) z = v[2];

        if(w[0] < minX) minX = w[0];
        if(w[0] > maxX) maxX = w[0];

        if(w[1] < minY) minY = w[1];
        if(w[1] > maxY) maxY = w[1];

        poly->push_back(Vertex(w, context->addVertex
                               (mat * Vector3(v[0], v[1], z))));
      }

      const unsigned numberOfPatchesX =
        static_cast<unsigned>(ceil((maxX-minX)/maxPatchSize));
      const unsigned numberOfPatchesY =
        static_cast<unsigned>(ceil((maxY-minY)/maxPatchSize));

      const double patchSizeX = (maxX-minX)/numberOfPatchesX;
      const double patchSizeY = (maxY-minY)/numberOfPatchesY;

      Vertex *reuseVertex = new Vertex[numberOfPatchesX-1];
      bool *reuse = new bool[numberOfPatchesX-1];

      for(unsigned i=0; i<numberOfPatchesX-1; ++i) reuse[i] = false;

      for(unsigned y=0; y<numberOfPatchesY; ++y)
      {
        if(y < numberOfPatchesY-1)
        {
          temp->clear();
          stripe->clear();
          splitY(*poly, minY + (y+1) * patchSizeY, mat,
                 context, *stripe, *temp, z);
          swap(poly, temp);
        }
        else swap(stripe, poly);

        for(unsigned x=0; x<numberOfPatchesX; ++x)
        {
          if(x < numberOfPatchesX-1)
          {
            temp->clear();
            patch->clear();
            splitX(*stripe, minX + (x+1) * patchSizeX, mat, context,
                   *patch, *temp, reuseVertex[x], reuse[x], z);
            swap(stripe, temp);
          }
          else swap(patch, stripe);

          if(patch->size() == 0) continue; // Patch and polygon are disjoint 

          if(patch->size() > 4)
          {
            int a = (*patch)[0].i;
            for(unsigned i=2; i<patch->size(); ++i)
            {
              vector<int> p;
              p.push_back(a);
              p.push_back((*patch)[i-1].i);
              p.push_back((*patch)[i].i);
              context->addPolygon(p);
            }
          }
          else
          {
            if(patch->size() < 3)
              ARCHON_THROW1(PolygonException, "Too few vertices");
            vector<int> p;
            p.push_back((*patch)[0].i);
            p.push_back((*patch)[1].i);
            p.push_back((*patch)[2].i);
            if(patch->size() == 4) p.push_back((*patch)[3].i);
            context->addPolygon(p);
          }
        }
      }

      delete[] reuseVertex;
      delete[] reuse;
    }

    namespace
    {
      struct Context: MesherContext
      {
        vector<Vector3> &vertices;
        vector<int> &polygons;
        Context(vector<Vector3> &vertices,
                vector<int> &polygons):
          vertices(vertices), polygons(polygons) {}

        int addVertex(const Vector3 &v)
        {
          int i = vertices.size();
          vertices.push_back(v);
          return i;
        }

        void addPolygon(const vector<int> &p)
        {
          polygons.insert(polygons.end(), p.begin(), p.end());
          polygons.push_back(-1);
        }
      };
    }

    void meshPolygon(const vector<Vector3> &polygon,
                     double maxPatchSize,
                     vector<Vector3> &resultingVertices,
                     vector<int> &resultingPolygons)
    {
      resultingVertices.clear();
      resultingPolygons.clear();
      Context context(resultingVertices, resultingPolygons);
      meshPolygon(polygon, maxPatchSize, &context);
    }
  }
}

---