triangulator.cpp

Go to the documentation of this file.

/*
    Sheep - A Rigid Body Dynamics Engine
    Copyright (C) 2001-2004 Francois Beaune
    Contact: http://toxicengine.sourceforge.net/

    This file is part of Sheep.

    Sheep is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    Sheep 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 General Public License
    along with Sheep; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

// Code for the ComputePolygonOrientation2() function comes from:
// http://geometryalgorithms.com/Archive/algorithm_0101/#orientation2D_triangle()
// Copyright 2000, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.

// Code for the ComputePolygonArea2() and ComputePolygonNormal3() functions
// comes from a paper entitled "Fast Polygon Area and Newell Normal Computation",
// Daniel Sunday, journal of graphics tools, 7(2):9-13, 2002.
// http://www.acm.org/jgt/papers/Sunday02/FastArea.html

#include "triangulator.h"   // include first
#include "common/math/build_basis.h"
#include "common/math/point2.h"
#include "common/math/point3.h"
#include "common/math/real.h"
#include "common/math/vector2.h"
#include "common/math/vector3.h"

#include <cassert>
#include <list>

using namespace sheep;
using namespace std;

namespace {
    //! Computes the normal to the plane containing the specified 3D triangle.
    //! Computed vector *is not* unit length.
    Vector3 ComputeTrianglePlane3(const Point3 &v0,
                                  const Point3 &v1,
                                  const Point3 &v2)
    {
        const Vector3 e0 = v1 - v0;
        const Vector3 e1 = v2 - v0;

        assert(e0.SquareNorm() > 0.0);
        assert(e1.SquareNorm() > 0.0);

        return e0 ^ e1;
    }

    //! Returns the signed area of a 2D polygon.
    Real ComputePolygonArea2(vector<Point2> input) {
        const int n = input.size();

        // Duplicate the first two vertices at the end of the input vector.
        input.push_back(input[0]);
        input.push_back(input[1]);

        Real area = 0.0;

        for(int i = 1, j = 2, k = 0; i <= n; ++i, ++j, ++k) {
            area += input[i].m_x * (input[j].m_y - input[k].m_y);
        }

        return 0.5 * area;
    }

    // Outputs the approximate unit normal of a 3D nearly planar polygon.
    // Returns the area of the polygon.
    Real ComputePolygonNormal3(vector<Point3> input, Vector3 *normal) {
        assert(normal);

        const int n = input.size();

        // Duplicate the first two vertices at the end of the input vector.
        input.push_back(input[0]);
        input.push_back(input[1]);

        normal->m_x = normal->m_y = normal->m_z = 0.0;

        // Get the Newell normal.
        for(int i = 1, j = 2, k = 0; i <= n; ++i, ++j, ++k) {
            normal->m_x += input[i].m_y * (input[j].m_z - input[k].m_z);
            normal->m_y += input[i].m_z * (input[j].m_x - input[k].m_x);
            normal->m_z += input[i].m_x * (input[j].m_y - input[k].m_y);
        }

        // Area of polygon = length of Newell normal.
        const Real area = normal->Norm();

        // Normalize the polygon normal.
        assert(area > 0.0);
        *normal /= area;

        return area;
    }

    //! Projects an approximatly planar 3D polygon to a plane.
    void ProjectPlanarPolygon3(const vector<Point3> &input,
                               vector<Point2> *output)
    {
        assert(input.size() >= 3);
        assert(output);

        // Compute normal to the plane containing the polygon.
        Vector3 n;
        ComputePolygonNormal3(input, &n);

        // Build an orthonormal basis around the normal.
        Vector3 u, v;
        build_basis(n, &u, &v);

        // Project the 3D polygon to the plane.
        for(vector<Point3>::const_iterator i = input.begin(), e = input.end(); i != e; ++i) {
            const Point3 &p = *i;
            output->push_back(Point2(p * u, p * v));
        }

/*      
        assert(input.size() >= 3);
        assert(output);

        // Compute normal to the plane containing the polygon.
        Vector3 n;
        ComputePolygonNormal3(input, &n);

        const Real nx = fabs(n.m_x);
        const Real ny = fabs(n.m_y);
        const Real nz = fabs(n.m_z);

        enum { X, Y, Z } ignore;

        // Search the largest coordinate (in absolute value) of the normal vector.
        if(nx > ny) {
            if(nx > nz)
                ignore = X;
            else ignore = Z;
        } else {
            if(ny > nz)
                ignore = Y;
            else ignore = Z;
        }

        output->reserve(input.size());

        switch(ignore) {
        case X:
            assert(nx > 0.0);
            for(vector<Point3>::const_iterator i = input.begin(), e = input.end(); i != e; ++i)
                output->push_back(Point2(i->m_y, i->m_z));
            break;
        case Y:
            assert(ny > 0.0);
            for(vector<Point3>::const_iterator i = input.begin(), e = input.end(); i != e; ++i)
                output->push_back(Point2(i->m_x, i->m_z));
            break;
        case Z:
            assert(nz > 0.0);
            for(vector<Point3>::const_iterator i = input.begin(), e = input.end(); i != e; ++i)
                output->push_back(Point2(i->m_x, i->m_y));
            break;
        }*/
    }

    enum Orientation { CCW, CW };

    //! Computes the orientation of a 2D triangle.
    Orientation ComputeTriangleOrientation2(const Point2 &v0,
                                            const Point2 &v1,
                                            const Point2 &v2)
    {
        const Vector2 v0v1 = v1 - v0;
        const Vector2 v0v2 = v2 - v0;
        const Real cp = v0v1.m_x * v0v2.m_y - v0v2.m_x * v0v1.m_y;

        // Make sure the triangle is not degenerate.
//      assert(cp != 0.0);

        return cp > 0.0 ? CCW : CW;
    }

    //! Computes the orientation of a 2D polygon.
    Orientation ComputePolygonOrientation2(const vector<Point2> &vertices) {
        assert(vertices.size() >= 3);

        const int n = vertices.size();

        // First find rightmost lowest vertex of the polygon.
        int rmin = 0;
        Real xmin = vertices[0].m_x;
        Real ymin = vertices[0].m_y;

        for(int i = 1; i < n; ++i) {
            const Point2 &v = vertices[i];

            if(v.m_y > ymin)
                continue;

            if(v.m_y == ymin) {     // just as low
                if(v.m_x < xmin)    // and to left
                    continue;
            }

            rmin = i;   // a new rightmost lowest vertex
            xmin = v.m_x;
            ymin = v.m_y;
        }

        if(rmin == 0)
            return ComputeTriangleOrientation2(
                vertices[n - 1],
                vertices[0],
                vertices[1]);
        else if(rmin == n - 1)
            return ComputeTriangleOrientation2(
                vertices[n - 2],
                vertices[n - 1],
                vertices[0]);
        else
            return ComputeTriangleOrientation2(
                vertices[rmin - 1],
                vertices[rmin],
                vertices[rmin + 1]);
    }

    //! Returns true if point 'p' is inside the triangle (v0, v1, v2).
    //! Triangle vertices must be given in counter-clockwise order.
    bool IsPointInsideTriangle2(const Point2 &v0,
                                const Point2 &v1,
                                const Point2 &v2,
                                const Point2 &p)
    {
        assert(ComputeTriangleOrientation2(v0, v1, v2) == CCW);

        const Vector2 v0p = p - v0;
        const Vector2 v0v1 = v1 - v0;

        if(v0p.m_x * v0v1.m_y - v0v1.m_x * v0p.m_y > 0.0)
            return false;

        const Vector2 v1p = p - v1;
        const Vector2 v1v2 = v2 - v1;

        if(v1p.m_x * v1v2.m_y - v1v2.m_x * v1p.m_y > 0.0)
            return false;

        const Vector2 v2p = p - v2;
        const Vector2 v2v0 = v0 - v2;

        if(v2p.m_x * v2v0.m_y - v2v0.m_x * v2p.m_y > 0.0)
            return false;

        return true;
    }

    bool IsAnEar(const vector<Point2> &vertices,
                 const list<int> &active,
                 list<int>::const_iterator prev,
                 list<int>::const_iterator cur,
                 list<int>::const_iterator next)
    {
        const Point2 &v0 = vertices[*prev];
        const Point2 &v1 = vertices[*cur];
        const Point2 &v2 = vertices[*next];

        if(ComputeTriangleOrientation2(v0, v1, v2) != CCW)
            return false;

        for(list<int>::const_iterator i = active.begin(), e = active.end(); i != e; ++i) {
            if(i == prev || i == cur || i == next)
                continue;

            if(IsPointInsideTriangle2(v0, v1, v2, vertices[*i]))
                return false;
        }

        return true;
    }

    template<typename Iterator> inline
    Iterator succ(Iterator i) {
        return ++i;
    }
}

bool sheep::TriangulatePolygon3(const vector<Point3> &vertices,
                                vector<int> *triangles)
{
    assert(vertices.size() >= 3);
    assert(triangles);

    vector<Point2> vertices2;
    ProjectPlanarPolygon3(vertices, &vertices2);

    const int n = vertices2.size();

    list<int> active;

    switch(ComputePolygonOrientation2(vertices2)) {
    case CCW:
        // Keep vertices order.
        for(int i = 0; i < n; ++i)
            active.push_back(i);
        break;
    case CW:
        // Reverse vertices order.
        for(int i = n - 1; i >= 0; --i)
            active.push_back(i);
        break;
    }

    triangles->reserve(3 * (n - 2));

    list<int>::iterator prev = active.begin(), cur = succ(prev), next = succ(cur);
    int iter = 0;

    while(active.size() > 2) {
        if(IsAnEar(vertices2, active, prev, cur, next)) {
            // Create a new triangle.
            triangles->push_back(*prev);
            triangles->push_back(*cur);
            triangles->push_back(*next);

            // Remove the current vertex from the active list and
            // update 'cur' iterator.
            cur = active.erase(cur);
            if(cur == active.end())
                cur = active.begin();

            // Update 'next' iterator.
            next = succ(cur);
            if(next == active.end())
                next = active.begin();

            // Reset the iteration counter.
            iter = 0;
        } else {
            // Update all three iterators.
            prev = cur;
            cur = next;
            next = succ(next);
            if(next == active.end())
                next = active.begin();

            // Update the iteration counter and detect infinite loop.
            if(++iter > 2 * n)  // weak
                return false;
        }
    }

    return true;
}

---