polyhedron.cpp

/***************************************************************************
 *   Copyright (C) 2007 by Pablo Diaz-Gutierrez   *
 *   pablo@ics.uci.edu   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 2 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <algorithm>
#include <float.h>
#include "polyhedron.h"
#include "polyhedron_gl.h"
#include "util.h"

using namespace std;
using namespace HE;
using namespace Geometry;

Polyhedron::Polyhedron(const char* filename)
  : _gl(0), _type(FT_OFF)
{
  if (!load(filename))
    throw HE_exception("Cannot load input file", filename);
}

Polyhedron::Polyhedron()
   : _gl(0), _type(FT_OFF)
{
}

Polyhedron::Polyhedron(const std::vector<Vertex*>& verts, const std::vector<std::vector<int> >& faces)
   : _gl(0), _type(FT_OFF)
{
  _vertices.resize(verts.size());
  for (unsigned v=0 ; v<verts.size() ; v++)
    _vertices[v] = new Vertex(verts[v]->position(), v);
  loadFaces( faces );
}

Polyhedron::~Polyhedron()
{
        for (face_iterator fit=fBegin() ; fit!=fEnd() ; ++fit)
                delete *fit;
        _faces.clear();

        for (edge_iterator eit=eBegin() ; eit!=eEnd() ; ++eit)
                delete *eit;
        _halfEdges.clear();

        for (vertex_iterator vit=vBegin() ; vit!=vEnd() ; ++vit)
                delete *vit;
        _vertices.clear();
}


Polyhedron* Polyhedron::refine() const
{
  Polyhedron* P = new Polyhedron;
  cerr << P->numFaces() << " faces in initial mesh\n";
  P->_vertices.reserve(_vertices.size() + _halfEdges.size()/2 + 1);
  for (unsigned i=0 ; i<_vertices.size() ; i++)
     P->_vertices.push_back( new Vertex(_vertices[i]->position(), _vertices[i]->index()) );

  map<const HalfEdge*, const Vertex*> edge2midpoints;
  for (const_face_iterator it=fBegin() ; it!=fEnd() ; ++it)
  {
    const Face* f = *it;
    //cerr << "***\nFace " << *f << "\n***\n";
    if (f->hole())
      continue;
    const HalfEdge* sentinel = f->edge();
    const HalfEdge* he = sentinel;
    const Vertex* V(0), *V2(0);
    vector<int> verts;
    //cerr << "[ Entering the loop...\n";
    do
    {
      if (edge2midpoints.end() == edge2midpoints.find(he))
      {
        V2 = P->addVertex(he->midpoint());
        edge2midpoints[he] = edge2midpoints[he->twin()] = V2;
      }
      else
      {
        V2 = edge2midpoints[he];
      }
      verts.push_back(V2->index());

      // Is this NOT the first halfedge we traverse in the face?
      if (V)
      {
        Face* f2 = P->addFace(V->index(), he->src()->index(), V2->index());
        f2->hole(f->hole());
      }

      V = V2;
      he = he->next();
    } while (he != sentinel) ;
    //cerr << "] Leaving the loop...\n";
    assert(verts.size() >= 3);
    Face* f2 = P->addFace(verts.back(), sentinel->src()->index(), verts.front());
    f2->hole(f->hole());
    f2 = P->addFace(verts);
    f2->hole(f->hole());
  }

  cerr << P->numFaces() << " in new Polyhedron\n";
  P->finalize();
  cerr << P->numFaces() << " in refined Polyhedron\n";
  cerr << numVertices() << " vertices in original Polyhedron\n";
  cerr << P->numVertices() << " vertices in refined Polyhedron\n";
  return P;
}


Vertex* Polyhedron::addVertex(const Vector3Df& v)
{
        Vertex* V = new Vertex(v, _vertices.size());
        _vertices.push_back(V);
        return V;
}

void Polyhedron::clearData()
{
   _faces.clear();
   _holes.clear();
   _halfEdges.clear();
   _vertices.clear();
}

bool Polyhedron::load(const char* filename)
{
        filetype_e format = fileFormat(filename);

   clearData();
        bool success = false;
        switch( format )
        {
                case FT_PLY_ASCII:
                        success = readPly(filename);
                        type(FT_PLY_ASCII);
                        break;

    case FT_OFF:
      success = readOff(filename);
      type(FT_OFF);
      break;

    case FT_TRI:
      success = readTri(filename);
      type(FT_TRI);
      break;

    case FT_PLY_BIN:
                        throw HE_exception("Cannot read binary PLY format yet", filename);
                        type(FT_PLY_BIN);
                        break;

                case FT_OBJ:
                        throw HE_exception("Cannot read OBJ format yet", filename);
                        type(FT_OBJ);
                        break;

                case FT_NONE:
                        throw HE_exception("Unknown input file format", filename);
                        type(FT_NONE);
                        break;
        }

        packVertices();
        return success;
}


void Polyhedron::finalize()
{
  assert(check(false));

  // First, find the connections that need finalization. During construction, this is all of them.
  map<pair<int,int>, HalfEdge*> connections;
  for (edge_iterator it=eBegin() ; it!=eEnd() ; ++it)
  {
    HalfEdge* he = *it;
    if (!he->twin())
    {
      int a = he->prev()->dst()->index();
      int b = he->dst()->index();
      pair<int,int> idx(a, b);
      if (connections.end() != connections.find(idx))
         throw HE_exception("Repeated half-edge: This is a non-manifold mesh, not supported by half-edge data structures.");

      pair<int,int> idx2(b, a);
      map<pair<int,int>, HalfEdge*>::iterator cit = connections.find(idx2);
      if (connections.end() == cit)
      {
         connections[idx] = he;
      }
      else
      {
         HalfEdge*& he2 = cit->second;
         he->twin(he2);
         he2->twin(he);
         connections.erase(cit);
      }
      //connections[make_pair(b, a)] = he;
    }
    else
    {
       assert(he->twin()->twin() == he);
    }
  }

  assert(check(false));
  closeSurface(connections);
  assert(check());
}


void Polyhedron::closeSurface(map<pair<int,int>, HalfEdge*>& connections)
{
  //cerr << "closeSurface\n";
   map<pair<int,int>, HalfEdge*>& endPoints = connections; // Which vertex points are connected by which halfedge?
   map<int,vector<pair<int,int> > >v2e;    // Which edges point to each vertex?
   for (map<pair<int,int>, HalfEdge*>::const_iterator it=endPoints.begin() ; it!=endPoints.end() ; ++it)
   {
      HalfEdge* e = it->second;
      assert(!e->twin());

      int src = e->prev()->dst()->index();
      int dst = e->dst()->index();
      pair<int,int> pts(src, dst);
      v2e[dst].push_back(pts);
      assert(find(v2e[pts.second].begin(), v2e[pts.second].end(), pts) != v2e[pts.second].end());
      //if (v2e[dst].size() > 1)
         //throw HE_exception("This mesh has a pinched vertex, which we do not support yet.");
   }

        // Add half-edges around the detected boundaries, each connected set being assigned to a face/hole
   while (!endPoints.empty())
   {
      map<pair<int,int>, HalfEdge*>::iterator pit = endPoints.begin();
      pair<int,int> pts = pit->first;
      HalfEdge* e = pit->second;
      HalfEdge* e2;
      Face* f = new Face(0, true);
      _faces.push_back(f);
      _holes.insert(f);

      //cerr << "processing edge <" << pts.first << ',' << pts.second << ">\n";
      int start = pts.second; // pts.second is the dst() of a surface half-edge
      int v = start;
      vector<HalfEdge*> bdry;
      do
      {
         e2 = new HalfEdge(_vertices[pts.first], f, e, 0);
         _halfEdges.push_back(e2);
         e->twin(e2);
         assert(e2->twin() == e);
         assert(e->dst() == e2->src());
         assert(e->src() == e2->dst());
         bdry.push_back(e2);
         assert(endPoints.end() != endPoints.find(pts));
         endPoints.erase(pit);
         assert(endPoints.end() == endPoints.find(pts));
         //assert(! v2e[pts.first].empty() );

         // Move on to the next halfedge in the boundary sequence.
         // Notice here that due to the edge shift, pts.first becomes pts.second... This caused a bug!
         int idx = pts.first;
         pts = v2e[idx].back();
         v2e[idx].pop_back();
         assert(pts.second == idx);
         //cerr << "trying to move to edge <" << pts.first << ',' << pts.second << ">\n";
         pit = endPoints.find(pts);
         if(endPoints.end() == pit)
            break;
         assert(pit->first == pts);
         //cerr << "moved to edge <" << pts.first << ',' << pts.second << ">\n";
         e = pit->second;
         v = pts.second;
         assert(start != v);
      } while(start != v);
      f->edge(e2);

                // Connect each new half-edge with its immediately next around the hole
      //cerr << "closing a loop w/ " << bdry.size() << " elements:\n";
      for (unsigned i=0 ; i<bdry.size() ; i++)
      {
        //cerr << '<' << bdry[i]->src()->index() << ',' << bdry[i]->dst()->index() << "> is followed by <"
          //  << bdry[(i+1)%bdry.size()]->src()->index() << ',' << bdry[(i+1)%bdry.size()]->dst()->index() << ">\n";
         bdry[i]->next( bdry[(i+1)%bdry.size()] );
      }

      // Check the connections
#ifndef NDEBUG
      for (unsigned i=0 ; i<bdry.size() ; i++)
      {
         const HalfEdge* e = bdry[i];
         const HalfEdge* e2 = e->twin();
         assert(e->prev()->dst() == e->src());
         assert(e2->prev()->dst() == e2->src());
      }
#endif
   }
}



HalfEdge* Polyhedron::edge(int from, int to) const
{
        return _vertices[from]->edgeTo(to);
}


void Polyhedron::boundingBox(Geometry::Vector3Df& m, Geometry::Vector3Df& M) const
{
        m = Vector3Df(FLT_MAX, FLT_MAX, FLT_MAX);
        M = -m;

        // Find the bounding min and max values for each coordinate, each point
        for (const_vertex_iterator vit=vBegin() ; vit!=vEnd() ; ++vit)
        {
                const Vector3Df& v = (*vit)->position();
                for (int i=0 ; i<3 ; i++)
                {
                        if (v[i] > M[i])
                                M[i] = v[i];
                        if (v[i] < m[i])
                                m[i] = v[i];
                }
        }
}


void Polyhedron::scaleAndCenter(float radius)
{
        Geometry::Vector3Df m, M;
        boundingBox(m, M);
        Geometry::Vector3Df ctr = (m+M) / 2;
        float diff = 0;
        for (int i=0 ; i<3 ; i++)
        {
                assert(M[i] >= m[i]);
                if (diff < M[i]-m[i])
                        diff = M[i]-m[i];
        }

        for (vertex_iterator vit=vBegin() ; vit!=vEnd() ; ++vit)
        {
                Vector3Df v = (*vit)->position() - ctr;
                v *= radius / diff;
                (*vit)->position(v);
        }

        if (_gl)
                _gl->recomputeNormals();
}


void Polyhedron::invertNormals()
{
        // Invert the order of the edges around the faces
        map<const HalfEdge*, HalfEdge*> next;
        map<const HalfEdge*, Vertex*> dst;
        for (edge_iterator eit=eBegin() ; eit!=eEnd() ; ++eit)
        {
                HalfEdge* e = *eit;

                assert(next.end() == next.find(e->next()));
                next[e->next()] = e;

                assert(dst.end() == dst.find(e));
                dst[e] = e->src();
        }

        // Apply such changes
        for (edge_iterator eit=eBegin() ; eit!=eEnd() ; ++eit)
        {
                HalfEdge* e = *eit;
                e->next( next[e] );

                Vertex* v = e->dst();
                e->dst( dst[e] );
                v->edge(e);
        }

        if (_gl)
                _gl->recomputeNormals();

        assert(check());
}


bool faceIsHole(const Face* f) { return f->hole(); }

int Polyhedron::numHoles() const
{
        return count_if(fBegin(), fEnd(), faceIsHole);
}


bool Polyhedron::check(bool holesFilled) const
{
#if 0
        for (Polyhedron::const_vertex_iterator it=vBegin() ; it!=vEnd() ; ++it)
        {
                assert((*it)->edge() != 0);
        }
#endif

        for (Polyhedron::const_edge_iterator it=eBegin() ; it!=eEnd() ; ++it)
        {
                if (holesFilled)
                {
                        assert((*it)->twin() != 0);
                        assert((*it)->twin()->twin() == (*it));
                }
                assert((*it)->next() != 0);
                assert((*it)->face() != 0);
      assert((*it)->dst() != 0);
      assert(!(*(*it)->dst() == *(*it)->prev()->dst()));
                
                // Make sure the half edge is correctly linked in a face
                {
                        set<const HalfEdge*> visited;
                        const HalfEdge* e = *it;
                        const HalfEdge* sentinel = e;
                        do
                        {
                                assert(visited.end() == visited.find(e));
                                visited.insert(e);
                                e = e->next();
                        } while(e != sentinel);
                }

                assert((*it)->prev()->next() == *it);
                if (holesFilled)
      {
        assert((*it)->prev()->dst() == (*it)->src());
      }
        }

        for (Polyhedron::const_face_iterator it=fBegin() ; it!=fEnd() ; ++it)
        {
                set<const HalfEdge*> visited;
                const Face* f = *it;
                Face::const_edge_circulator eit = f->begin();
                Face::const_edge_circulator sentinel = eit;
                do
                {
                        assert(visited.end() == visited.find(*eit));
                        visited.insert(*eit);
                        ++eit;
                } while(eit != sentinel);
        }

        return true;
}


void Polyhedron::packVertices()
{
        int n = 0;
        for (int v=0 ; v<numVertices() ; v++)
        {
                if (_vertices[v]->edge())
                {
                        _vertices[v]->index(n);
                        _vertices[n++] = _vertices[v];
                }
        }

        //cerr << (numVertices()-n) << " vertices have been eliminated\n";
        _vertices.resize(n);
}


bool Polyhedron::adjacent(const Vertex* v1, const Vertex* v2) const
{
   Vertex::const_edge_circulator eit = v1->begin();
   Vertex::const_edge_circulator sentinel = eit;
   do
   {
      if ((*eit)->dst() == v2)
         return true;
      ++eit;
   } while(eit != sentinel);
   return false;
}

bool Polyhedron::adjacent(const Face* f1, const Face* f2) const
{
   Face::const_edge_circulator eit = f1->begin();
   Face::const_edge_circulator sentinel = eit;
   do
   {
      if ((*eit)->twin()->face() == f2)
         return true;
      ++eit;
   } while(eit != sentinel);
   return false;
}

---