polygon.cpp

Go to the documentation of this file.

/* Copyright, 2000 Nevrax Ltd.
 *
 * This file is part of NEVRAX NEL.
 * NEVRAX NEL 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, or (at your option)
 * any later version.

 * NEVRAX NEL 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 NEVRAX NEL; see the file COPYING. If not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

#include "stdmisc.h"

#include "nel/misc/polygon.h"
#include "nel/misc/plane.h"
#include "nel/misc/triangle.h"


using namespace std;
using namespace NLMISC;


namespace NLMISC
{


//==================================//
//      CPolygon implementation     //
//==================================//

// ***************************************************************************
CPolygon::CPolygon(const CVector &a, const CVector &b, const CVector &c)
{
    Vertices.reserve(3);
    Vertices.push_back(a);
    Vertices.push_back(b);
    Vertices.push_back(c);
}

// ***************************************************************************
void CPolygon::toTriFan(std::vector<NLMISC::CTriangle> &dest) const
{
    sint count = (sint) Vertices.size() - 2;
    for(sint k = 0; k < count; ++k)
    {
        dest.push_back(CTriangle(Vertices[0], Vertices[k + 1], Vertices[k + 2]));
    }
}

// ***************************************************************************
float CPolygon::computeArea() const
{
    float area = 0.f;
    sint numVerts = (sint) Vertices.size();
    for(sint k = 0; k < numVerts - 2; ++k)
    {
        CVector v0 = Vertices[k + 1] - Vertices[0];
        CVector v1 = Vertices[k + 2] - Vertices[0];
        area += (v0 ^ v1).norm();
    }
    return 0.5f * fabsf(area);
}

// ***************************************************************************
void            CPolygon::clip(const CPlane  *planes, uint nPlanes)
{
    if(nPlanes==0 || getNumVertices()==0)
        return;

    // The final polygon has at maximum currentVertices+number of clipping planes.
    // For performance, the vectors are static, so reallocation rarely occurs.
    static  vector<CVector>     tab0, tab1;
    tab0.resize(getNumVertices()+nPlanes);
    tab1.resize(getNumVertices()+nPlanes);
    // Init tab0 with Vertices.
    copy(Vertices.begin(), Vertices.end(), tab0.begin());
    CVector             *in=&(*tab0.begin()), *out= &(*tab1.begin());
    sint                nin= getNumVertices(), nout;
    for(sint i=0;i<(sint)nPlanes;i++)
    {
        nout= planes[i].clipPolygonBack(in, out, nin);
        swap(in, out);
        nin= nout;
        if(nin==0)
            break;
    }

    // Final result in "in".
    Vertices.resize(nin);
    if(nin>0)
    {
        memcpy(&(*Vertices.begin()), in, nin*sizeof(CVector));
    }
}


// ***************************************************************************
void            CPolygon::clip(const std::vector<CPlane> &planes)
{
    if(planes.size()==0)
        return;
    clip(&(*planes.begin()), planes.size());
}



// ***************************************************************************
void CPolygon::serial(NLMISC::IStream &f) throw(NLMISC::EStream)
{
    f.serialVersion(0);
    f.serialCont(Vertices);
}

// ***************************************************************************
void CPolygon::getBestTriplet(uint &index0,uint &index1,uint &index2)
{
    nlassert(Vertices.size() >= 3);
    uint i, j, k;
    float bestArea = 0.f;
    const uint numVerts = Vertices.size();
    for (i = 0; i < numVerts; ++i)
    {
        for (j = 0; j < numVerts; ++j)
        {
            if (i != j)
            {
                for (k = 0; k < numVerts; ++k)
                {
                    if (k != i && k != j)
                    {
                        CVector v0 = Vertices[j] - Vertices[i];
                        CVector v1 = Vertices[k] - Vertices[i];
                        float area = (v0 ^ v1).norm();
                        if (area > bestArea)
                        {
                            bestArea = area;
                            index0 = i;
                            index1 = j;
                            index2 = k;
                        }
                    }
                }
            }
        }
    }

}

// ***************************************************************************
void CPolygon::buildBasis(CMatrix &dest)
{
    nlassert(Vertices.size() > 3);
    uint i1, i2, i3;
    getBestTriplet(i1, i2, i3);
    CVector v1 = (Vertices[i2] - Vertices[i1]).normed();
    CVector v2 = (Vertices[i3] - Vertices[i1]).normed();
    CVector K = v2 ^ v1;
    CVector I = v1 - (v1 * K) * v1;
    CVector J = K ^ I;
    dest.setRot(I, J, K);
    dest.setPos(Vertices[i1]);
}

// ***************************************************************************

class CConcavePolygonsVertexDesc
{
public:

    CConcavePolygonsVertexDesc (float length, uint index)
    {
        Length = length;
        Index = index;
    }

    // Length > 0
    float   Length;

    // First vertex index
    uint    Index;
};
typedef std::map<float, CConcavePolygonsVertexDesc> TCConcavePolygonsVertexMap;


// ***************************************************************************
bool CPolygon::toConvexPolygonsEdgeIntersect (const CVector2f& a0, const CVector2f& a1, const CVector2f& b0, const CVector2f& b1)
{
    // both vertical?
    if( a0.x-a1.x==0 && b0.x-b1.x==0 )
        return false;

    // compute intersection of both lines
    CVector2f intersection;

    // first edge vertical?
    if(a0.x - a1.x==0)
    {
        float Ab = (b0.y - b1.y) / (b0.x - b1.x);

        // Intersection
        intersection.x = a0.x;
        intersection.y = b0.y + (a0.x-b0.x) * Ab;
    }
    // second edge vertical?
    else if(b0.x - b1.x==0)
    {
        float Aa = (a0.y - a1.y) / (a0.x - a1.x);

        // Intersection
        intersection.x = b0.x;
        intersection.y = a0.y + (b0.x-a0.x) * Aa;
    }
    // standard case
    else
    {
        float Aa = (a0.y - a1.y) / (a0.x - a1.x);
        float Ba = a0.y - a0.x * Aa;
        float Ab = (b0.y - b1.y) / (b0.x - b1.x);
        float Bb = b0.y - b0.x * Ab;

        // colinear?
        if(Aa==Ab)
            return false;

        // Intersection
        intersection.x = (Bb - Ba) / (Aa - Ab);
        intersection.y = Aa * intersection.x + Ba;
    }

    // In it ?
    return ( ( (a0-intersection)*(a1-intersection) < 0 ) && ( (b0-intersection)*(b1-intersection) < 0 ) );
}

// ***************************************************************************

class CBSPNode2v
{
public:
    CBSPNode2v ()
    {
        Back = NULL;
        Front = NULL;
    }
    CBSPNode2v ( const CPlane &plane, CVector p0, CVector p1, uint v0, uint v1 ) : Plane (plane), P0 (p0), P1 (p1)
    {
        Back = NULL;
        Front = NULL;
        Parent = NULL;
        V0 = v0;
        V1 = v1;
    }
    ~CBSPNode2v ()
    {
        if (Front)
            delete Front;
        if (Back)
            delete Back;
    }

    void insert (CBSPNode2v *node)
    {
        // Front ?
        bool p0Front = (Plane * node->P0) > 0;
        bool p1Front = (Plane * node->P1) > 0;
        if (p0Front && p1Front)
        {
            // Front child ?
            if (Front)
                Front->insert (node);
            else
            {
                // Link left
                Front = node;
                node->Parent = this;
            }
        }
        else if ((!p0Front) && (!p1Front))
        {
            // Back child ?
            if (Back)
                Back->insert (node);
            else
            {
                // Link left
                Back = node;
                node->Parent = this;
            }
        }
        else
        {
            // Split vertex
            CVector newVertex = Plane.intersect (node->P0, node->P1);

            // New node
            CBSPNode2v *newNode = new CBSPNode2v (node->Plane, node->P0, newVertex, node->V0, node->V1);

            // Old node
            node->P0 = newVertex;

            // Insert child
            CBSPNode2v **p0Parent;
            CBSPNode2v **p1Parent;

            // Get insertion pointer
            if (p0Front)
            {
                p0Parent = &Front;
                p1Parent = &Back;
            }
            else
            {
                p0Parent = &Back;
                p1Parent = &Front;
            }

            // Insert children
            if (*p0Parent)
            {
                (*p0Parent)->insert (newNode);
            }
            else
            {
                *p0Parent = newNode;
                newNode->Parent = this;
            }

            // Insert children
            if (*p1Parent)
            {
                (*p1Parent)->insert (node);
            }
            else
            {
                *p1Parent = node;
                node->Parent = this;
            }
        }
    }

    bool intersect (const CVector &p0, const CVector &p1, uint v0, uint v1) const
    {
        // Front ?
        bool p0Front = (Plane * p0) > 0;
        bool p1Front = (Plane * p1) > 0;

        if (p0Front != p1Front)
            if ( (v0 != V0) && (v0 != V1) && (v1 != V0) && (v1 != V1) )
                if (CPolygon::toConvexPolygonsEdgeIntersect ((CVector2f) P0, (CVector2f) P1, (CVector2f) p0, (CVector2f) p1))
                    return true;

        if (p0Front || p1Front)
        {
            if (Front)
                if (Front->intersect (p0, p1, v0, v1))
                    return true;
        }

        if ((!p0Front) || (!p1Front))
        {
            if (Back)
                if (Back->intersect (p0, p1, v0, v1))
                    return true;
        }

        return false;
    }

    CBSPNode2v  *Back, *Front, *Parent;
    CPlane      Plane;
    CVector     P0;
    CVector     P1;
    uint        V0;
    uint        V1;
};

// ***************************************************************************

bool CPolygon::toConvexPolygonsLeft (const std::vector<CVector> &vertex, uint a, uint b, uint c)
{
    return ( (vertex[b].x - vertex[a].x) * (vertex[c].y - vertex[a].y) - (vertex[c].x - vertex[a].x) * (vertex[b].y - vertex[a].y) ) < 0;
}

// ***************************************************************************

bool CPolygon::toConvexPolygonsLeftOn (const std::vector<CVector> &vertex, uint a, uint b, uint c)
{
    return ( (vertex[b].x - vertex[a].x) * (vertex[c].y - vertex[a].y) - (vertex[c].x - vertex[a].x) * (vertex[b].y - vertex[a].y) ) <= 0;
}

// ***************************************************************************

bool CPolygon::toConvexPolygonsInCone (const std::vector<CVector> &vertex, uint a, uint b)
{
    // Prev and next
    uint a0 = a+1;
    if (a0==vertex.size())
        a0=0;
    uint a1;
    if (a==0)
        a1=vertex.size()-1;
    else
        a1= a-1;

    if (toConvexPolygonsLeftOn (vertex, a, a1, a0) )
    {
        return toConvexPolygonsLeft ( vertex, a, b, a0) && toConvexPolygonsLeft ( vertex, b, a, a1);
    }
    else
    {
        return !( toConvexPolygonsLeft ( vertex, a, b, a1) && toConvexPolygonsLeft ( vertex, b, a, a0) );
    }
}

// ***************************************************************************

bool CPolygon::toConvexPolygonsDiagonal (const std::vector<CVector> &vertex, const CBSPNode2v &bsp, uint a, uint b)
{
    // Check it is a border
    if ( ( (b - a) == 1) || ( (a - b) == 1) || ( (a==0) && (b ==(vertex.size()-1))) || ( (b==0) && (a ==(vertex.size()-1))) )
        return true;

    // Check visibility
    if (toConvexPolygonsInCone (vertex, a, b) && toConvexPolygonsInCone (vertex, b, a))
    {
        // Intersection ?
        return !bsp.intersect (vertex[a], vertex[b], a, b);
    }
    return false;
}

// ***************************************************************************

void CPolygon::toConvexPolygonsLocalAndBSP (std::vector<CVector> &localVertices, CBSPNode2v &root, const CMatrix &basis) const
{
    // Invert matrix
    CMatrix invert = basis;
    invert.invert ();

    // Insert vertices in an ordered table
    uint vertexCount = Vertices.size();
    TCConcavePolygonsVertexMap vertexMap;
    localVertices.resize (vertexCount);
    uint i, j;

    // Transform the vertex
    for (i=0; i<vertexCount; i++)
    {
        CVector local = invert*Vertices[i];
        localVertices[i] = CVector (local.x, local.y, 0);
    }

    // Plane direction
    i=0;
    j=Vertices.size()-1;
    CVector normal = localVertices[i] - localVertices[j];
    normal = normal ^ CVector::K;
    CPlane clipPlane;
    clipPlane.make(normal, localVertices[i]);

    // Build the BSP root
    root = CBSPNode2v (clipPlane, localVertices[i], localVertices[j], i, j);

    // Insert all others edges
    j=i++;
    for (; i<Vertices.size(); i++)
    {
        // Plane direction
        normal = localVertices[i] - localVertices[j];
        normal = normal ^ CVector::K;
        clipPlane.make(normal, localVertices[i]);

        // Build the BSP root
        root.insert ( new CBSPNode2v (clipPlane, localVertices[i], localVertices[j], i, j) );

        j=i;
    }
}


// ***************************************************************************
bool CPolygon::toConvexPolygons (std::list<CPolygon>& outputPolygons, const CMatrix& basis) const
{
    // Some vertices ?
    if (Vertices.size()>2)
    {
        // Local vertices
        std::vector<CVector>    localVertices;

        // Build the BSP root
        CBSPNode2v root;

        // Build the local array and the BSP
        toConvexPolygonsLocalAndBSP (localVertices, root, basis);

        // Build a vertex list
        std::list<uint> vertexList;
        uint i;
        for (i=0; i<Vertices.size(); i++)
            vertexList.push_back (i);

        // Clip ears while there is some polygons
        std::list<uint>::iterator current=vertexList.begin();
        std::list<uint>::iterator begin=vertexList.begin();
        do
        {
again:;
            // Search for a diagonal
            bool found = false;

            // Get next vertex
            std::list<uint>::iterator first = current;
            std::list<uint>::iterator lastPreviousPrevious=current;
            std::list<uint>::iterator lastPrevious=current;
            lastPrevious++;
            if (lastPrevious==vertexList.end())
                lastPrevious = vertexList.begin();
            std::list<uint>::iterator currentNext = lastPrevious;
            std::list<uint>::iterator last = lastPrevious;
            last++;
            if (last==vertexList.end())
                last = vertexList.begin();
            while (last != current)
            {
                // Is a diagonal ?
                if (
                    (toConvexPolygonsDiagonal (localVertices, root, *lastPreviousPrevious, *last)) &&
                    (toConvexPolygonsDiagonal (localVertices, root, *currentNext, *last)) &&
                    (toConvexPolygonsDiagonal (localVertices, root, *last, *current))
                    )
                {
                    // Find one
                    found = true;
                }
                else
                {
                    // Come back
                    last = lastPrevious;
                    lastPrevious = lastPreviousPrevious;
                    break;
                }

                // Next vertex
                lastPreviousPrevious = lastPrevious;
                lastPrevious = last++;
                if (last==vertexList.end())
                    last = vertexList.begin();
            }

            // Last polygon ?
            if (last==current)
            {
                // Add a polygon
                outputPolygons.push_back (CPolygon());
                CPolygon &back = outputPolygons.back ();
                back.Vertices.reserve (vertexList.size());

                // Add each vertex in the new polygon
                current=vertexList.begin();
                while (current!=vertexList.end())
                {
                    back.Vertices.push_back (Vertices[*current]);
                    current++;
                }

                // Exit
                return true;
            }
            else
            {
                std::list<uint>::iterator firstNext = current;
                std::list<uint>::iterator firstNextNext = currentNext;
                if (first != vertexList.begin())
                    first--;
                else
                {
                    first = vertexList.end();
                    first--;
                }

                while (current != first)
                {
                    // Is a diagonal ?
                    if (
                        (toConvexPolygonsDiagonal (localVertices, root, *firstNextNext, *first)) &&
                        (toConvexPolygonsDiagonal (localVertices, root, *lastPrevious, *first)) &&
                        (toConvexPolygonsDiagonal (localVertices, root, *last, *first))
                        )
                    {
                        // Find one
                        found = true;
                    }
                    else
                    {
                        // Come back
                        first = firstNext;
                        break;
                    }

                    // Next vertex
                    firstNextNext = firstNext;
                    firstNext = first;
                    if (first==vertexList.begin())
                    {
                        first = vertexList.end();
                        first--;
                    }
                    else
                        first--;
                }
            }

            // Found ?
            if (found)
            {
                // Count vertex
                outputPolygons.push_back (CPolygon());
                CPolygon &back = outputPolygons.back ();

                // Vertex count
                uint vertexCount = 1;
                current = first;
                while (current != last)
                {
                    vertexCount++;
                    current++;
                    if (current == vertexList.end())
                        current = vertexList.begin();
                }

                // Alloc vertices
                back.Vertices.reserve (vertexCount);

                // Copy and remove vertices
                back.Vertices.push_back (Vertices[*first]);
                first++;
                if (first == vertexList.end())
                    first = vertexList.begin();
                while (first != last)
                {
                    back.Vertices.push_back (Vertices[*first]);

                    // Remove from list
                    first = vertexList.erase (first);
                    if (first == vertexList.end())
                        first = vertexList.begin();
                    nlassert (first != vertexList.end());
                }
                back.Vertices.push_back (Vertices[*first]);
                current = begin = last;
                goto again;
            }

            // Next current
            current++;
            if (current == vertexList.end())
                current = vertexList.begin ();
        }
        while (current != begin);
    }
    return false;
}

// ***************************************************************************

bool CPolygon::chain (const std::vector<CPolygon> &other, const CMatrix& basis)
{
    // Local vertices
    std::vector<CVector>    localVertices;

    // Build the BSP root
    CBSPNode2v root;

    // Build the local array and the BSP
    toConvexPolygonsLocalAndBSP (localVertices, root, basis);

    // Local vertices
    std::vector<std::vector<CVector> >  localVerticesOther (other.size());

    // Build the BSP root
    std::vector<CBSPNode2v> rootOther (other.size());

    // Build a copy of the polygons
    std::vector<CPolygon> copy = other;

    // Main copy
    CPolygon mainCopy = *this;

    // For each other polygons
    uint o;
    for (o=0; o<other.size(); o++)
    {
        // Build the local array and the BSP
        other[o].toConvexPolygonsLocalAndBSP (localVerticesOther[o], rootOther[o], basis);
    }

    // Look for a couple..
    uint thisCount = Vertices.size();
    uint i, j;
    for (o=0; o<other.size(); o++)
    {
        uint otherCount = other[o].Vertices.size();

        // Try to link in the main polygon
        for (i=0; i<thisCount; i++)
        {
            for (j=0; j<otherCount; j++)
            {
                // Test this segement
                if (!root.intersect (localVertices[i], localVerticesOther[o][j], i, 0xffffffff))
                {
                    // Test each other polygons
                    uint otherO;
                    for (otherO=0; otherO<other.size(); otherO++)
                    {
                        // Intersect ?
                        if (rootOther[otherO].intersect (localVertices[i], localVerticesOther[o][j], 0xffffffff, (otherO == o)?j:0xffffffff))
                            break;
                    }

                    // Continue ?
                    if (otherO==other.size())
                    {
                        // Insert new vertices
                        mainCopy.Vertices.insert (mainCopy.Vertices.begin()+i, 2+otherCount, CVector());

                        // Copy the first vertex
                        mainCopy.Vertices[i] = mainCopy.Vertices[i+otherCount+2];

                        // Copy the new vertices
                        uint k;
                        for (k=0; k<otherCount; k++)
                        {
                            uint index = j+k;
                            if (index>=otherCount)
                                index -= otherCount;
                            mainCopy.Vertices[i+k+1] = copy[o].Vertices[index];
                        }

                        // Copy the last one
                        mainCopy.Vertices[i+otherCount+1] = copy[o].Vertices[j];
                        break;
                    }
                }
            }
            if (j!=otherCount)
                break;
        }

        // Not found ?
        if (i==thisCount)
        {
            // Try to link in the sub polygons
            uint otherToCheck;
            for (otherToCheck=o+1; otherToCheck<other.size(); otherToCheck++)
            {
                uint otherToCheckCount = other[otherToCheck].Vertices.size();
                for (i=0; i<otherToCheckCount; i++)
                {
                    for (j=0; j<otherCount; j++)
                    {
                        // Test this segement
                        if (!rootOther[otherToCheck].intersect (localVerticesOther[otherToCheck][i], localVerticesOther[o][j], i, 0xffffffff))
                        {
                            // Test each other polygons
                            uint otherO;
                            for (otherO=0; otherO<other.size(); otherO++)
                            {
                                // Intersect ?
                                if (rootOther[otherO].intersect (localVerticesOther[otherToCheck][i], localVerticesOther[o][j],  (otherToCheck == otherO)?i:0xffffffff,  (otherO == o)?j:0xffffffff))
                                    break;
                            }

                            // Continue ?
                            if (otherO==other.size())
                            {
                                // Insert new vertices
                                copy[otherToCheck].Vertices.insert (copy[otherToCheck].Vertices.begin()+i, 2+otherCount, CVector());

                                // Copy the first vertex
                                copy[otherToCheck].Vertices[i] = copy[otherToCheck].Vertices[i+otherCount+2];

                                // Copy the new vertices
                                uint k;
                                for (k=0; k<otherCount; k++)
                                {
                                    uint index = j+k;
                                    if (index>=otherCount)
                                        index -= otherCount;
                                    copy[otherToCheck].Vertices[i+k+1] = copy[otherO].Vertices[index];
                                }

                                // Copy the last one
                                copy[otherToCheck].Vertices[i+otherCount+1] = copy[otherO].Vertices[j];
                                break;
                            }
                        }
                    }
                    if (j!=otherCount)
                        break;
                }
                if (i!=otherToCheckCount)
                    break;
            }
            if (otherToCheck==other.size())
            {
                // Not ok
                return false;
            }
        }
    }

    // Ok
    *this = mainCopy;
    return true;
}

// ***************************************************************************


//====================================//
//      CPolygon2d implementation     //
//====================================//



// ***************************************************************************
CPolygon2D::CPolygon2D(const CPolygon &src, const CMatrix &projMat)
{
    fromPolygon(src, projMat);
}

// ***************************************************************************
void CPolygon2D::fromPolygon(const CPolygon &src, const CMatrix &projMat /*=CMatrix::Identity*/)
{
    uint size = src.Vertices.size();
    Vertices.resize(size);
    for (uint k = 0; k < size; ++k)
    {
        CVector proj = projMat * src.Vertices[k];
        Vertices[k].set(proj.x, proj.y);
    }
}

// ***************************************************************************
bool        CPolygon2D::isConvex()
{
    bool Front  = true, Back = false;
    // we apply a dummy algo for now : check whether every vertex is in the same side
    // of every plane defined by a segment of this poly
    uint numVerts = Vertices.size();
    if (numVerts < 3) return true;
    CVector     segStart, segEnd;
    CPlane      clipPlane;
    for (TVec2fVect::const_iterator it = Vertices.begin(); it != Vertices.end(); ++it)
    {
        segStart.set(it->x, it->y, 0);            // segment start
        segEnd.set((it + 1)->x, (it + 1)->y, 0);  // segment end
        float n = (segStart - segEnd).norm();     // segment norm
        if (n != 0)
        {
            clipPlane.make(segStart, segEnd, (n > 10 ? n : 10) * CVector::K + segStart); // make a plane, with this segment and the poly normal
            // check each other vertices against this plane
            for (TVec2fVect::const_iterator it2 = Vertices.begin(); it2 != Vertices.end(); ++it2)
            {
                if (it2 != it && it2 != (it + 1)) // the vertices must not be part of the test plane (because of imprecision)
                {

                    float dist  = clipPlane * CVector(it2->x, it2-> y, 0);
                    if (dist != 0) // midlle pos
                    {
                        if (dist > 0) Front = true; else Back = true;
                        if (Front && Back) return false; // there are both front end back vertices -> failure
                    }
                }
            }
        }
    }
    return true;
}

// ***************************************************************************

void        CPolygon2D::buildConvexHull(CPolygon2D &dest) const
{
    nlassert(&dest != this);

    if (this->Vertices.size() == 3) // with 3 points it is always convex
    {
        dest = *this;
        return;
    }
    uint k, l;
    uint numVerts = Vertices.size();
    CVector2f p, curr, prev;
    uint      pIndex, p1Index, p2Index, pCurr, pPrev;
    // this is not optimized, but not used in realtime.. =)
    nlassert(numVerts >= 3);
    dest.Vertices.clear();

    typedef std::set<uint> TIndexSet;
    TIndexSet leftIndex;
    for (k = 0; k < Vertices.size(); ++k)
    {
        leftIndex.insert(k);
    }


    // 1) find the highest point p of the set. We are sure it belongs to the hull
    pIndex = 0;
    p = Vertices[0];
    for (k = 1; k < numVerts; ++k)
    {
        if (Vertices[k].y < p.y)
        {
            pIndex = k;
            p = Vertices[k];
        }
    }

    leftIndex.erase(pIndex);


    float bestCP = 1.1f;
    p1Index = p2Index = pIndex;

    for (k = 0; k < numVerts; ++k)
    {
        if (k != pIndex)
        {
            for (l = 0; l < numVerts; ++l)
            {
                if (l != pIndex && l != k)
                {
                    CVector2f seg1 = (Vertices[l] - p).normed();
                    CVector2f seg2 = (Vertices[k] - p).normed();

                    //CVector cp = CVector(seg1.x, seg1.y, 0) ^ CVector(seg2.x, seg2.y, 0);
                    //float n = fabsf(cp.z);
                    float n = seg1 * seg2;
                    if (n < bestCP)
                    {
                        p1Index = l;
                        p2Index = k;
                        bestCP  = n;
                    }
                }
            }
        }
    }


    leftIndex.erase(p2Index);



    // start from the given triplet, and complete the poly until we reach the first point
    pCurr = p2Index;
    pPrev = pIndex;

    curr = Vertices[pCurr];
    prev = Vertices[pPrev];

    // create the first triplet vertices
    dest.Vertices.push_back(Vertices[p1Index]);
    dest.Vertices.push_back(prev);
    dest.Vertices.push_back(curr);

    uint step = 0;

    for(;;)
    {
        bestCP = 1.1f;
        CVector2f seg2 = (prev - curr).normed();
        TIndexSet::iterator bestIt = leftIndex.end();
        for (TIndexSet::iterator it =  leftIndex.begin(); it != leftIndex.end(); ++it)
        {
            if (step == 0 && *it == p1Index) continue;
            CVector2f seg1 = (Vertices[*it] - curr).normed();
            float n = seg1 * seg2;
            if (n < bestCP)
            {
                bestCP = n;
                bestIt = it;
            }
        }

        nlassert(bestIt != leftIndex.end());
        if (*bestIt == p1Index)
        {
            return; // if we reach the start point we have finished
        }
        prev = curr;
        curr = Vertices[*bestIt];
        pPrev = pCurr;
        pCurr = *bestIt;
        // add new point to the destination
        dest.Vertices.push_back(curr);
        ++step;
        leftIndex.erase(bestIt);
    }
}

// ***************************************************************************


void CPolygon2D::serial(NLMISC::IStream &f) throw(NLMISC::EStream)
{
    (void)f.serialVersion(0);
    f.serialCont(Vertices);
}

// ***************************************************************************
void        CPolygon2D::getBestTriplet(uint &index0, uint &index1, uint &index2)
{
    nlassert(Vertices.size() >= 3);
    uint i, j, k;
    float bestArea = 0.f;
    const uint numVerts = Vertices.size();
    for (i = 0; i < numVerts; ++i)
    {
        for (j = 0; j < numVerts; ++j)
        {
            if (i != j)
            {
                for (k = 0; k < numVerts; ++k)
                {
                    if (k != i && k != j)
                    {
                        CVector2f v0 = Vertices[j] - Vertices[i];
                        CVector2f v1 = Vertices[k] - Vertices[i];
                        float area = fabsf((CVector(v0.x, v0.y, 0) ^ CVector(v1.x, v1.y, 0)).norm());
                        if (area > bestArea)
                        {
                            bestArea = area;
                            index0 = i;
                            index1 = j;
                            index2 = k;
                        }
                    }
                }
            }
        }
    }
}


// scan a an edge of a poly and write it into a table
static void ScanEdge(CPolygon2D::TRasterVect &outputVect, sint topY, const CVector2f &v1, const CVector2f &v2, bool rightEdge = true)
{
     const uint rol16 = 65536;
     sint ceilY1 = (sint) ceilf(v1.y);
     sint height;
     float deltaX, deltaY;
     float fInverseSlope;
     sint  iInverseSlope, iPosX;

     // check whether this segment gives a contribution to the final poly
     height = (sint) (ceilf(v2.y) - ceilY1);
     if (height <= 0) return;

     // compute slope
     deltaY = v2.y - v1.y;
     deltaX = v2.x - v1.x;
     fInverseSlope = deltaX / deltaY;


     CPolygon2D::TRasterVect::iterator  outputIt = outputVect.begin() + (ceilY1 - topY);

     // slope with ints
     iInverseSlope = (sint) (rol16 * fInverseSlope);

     // sub-pixel accuracy
     iPosX = (int) (rol16 * (v1.x + fInverseSlope * (ceilY1 - v1.y)));

     const CPolygon2D::TRasterVect::iterator endIt = outputIt + height;
     if (rightEdge)
     {
         do
         {
           outputIt->second =  iPosX >> 16;
           iPosX += iInverseSlope;
           ++outputIt;
         }
         while (outputIt != endIt);
     }
     else
     {
         iPosX += (rol16 - 1);
         do
         {
           outputIt->first =  iPosX >> 16;
           iPosX += iInverseSlope;
           ++outputIt;
         }
         while (outputIt != endIt);
     }
}


// *******************************************************************************
// This function alow to cycle forward through a vertex vector like if it was a circular list
static inline CPolygon2D::TVec2fVect::const_iterator Next(const CPolygon2D::TVec2fVect::const_iterator &it, const CPolygon2D::TVec2fVect &cont)
{
    nlassert(cont.size() != 0);
    if ((it + 1) == cont.end()) return cont.begin();
    return (it + 1);
}


// *******************************************************************************
// This function alow to cycle backward through a (non null) vertex vector like if it was a circular list
static inline CPolygon2D::TVec2fVect::const_iterator Prev(const CPolygon2D::TVec2fVect::const_iterator &it, const CPolygon2D::TVec2fVect &cont)
{
    nlassert(cont.size() != 0);
    if (it == cont.begin()) return cont.end() - 1;
    return (it - 1);
}


// *******************************************************************************
bool CPolygon2D::isCCWOriented() const
{
    const TVec2fVect &V = Vertices;
    nlassert(Vertices.size() >= 3);
    // compute highest and lowest pos of the poly
    float fHighest = V[0].y;
    float fLowest = fHighest;
    // iterators to the highest and lowest vertex
    TVec2fVect::const_iterator it = V.begin() ;
    const TVec2fVect::const_iterator endIt = V.end();
    TVec2fVect::const_iterator pHighest = V.begin();
    do
    {
        if (it->y < fHighest)
        {
            fHighest = it->y;
            pHighest = it;
        }
        fLowest = std::max(fLowest, it->y);
        ++it;
    }
    while (it != endIt);
    // we seek this vertex
    TVec2fVect::const_iterator pHighestRight = pHighest;
    if (fLowest == fHighest)
    {
        // special case : degenerate poly
        while (pHighestRight->x == pHighest->x)
        {
            pHighestRight = Next(pHighestRight, V);
            if (pHighestRight == pHighest) return false; // the poly is reduced to a point, returns an abritrary value
        }
        return pHighest->x <= pHighestRight->x;
    }
    // iterator to the first vertex that has an y different from the top vertex
    while (Next(pHighestRight, V)->y == fHighest)
    {
        pHighestRight = Next(pHighestRight, V);
    }

    // iterator to the first vertex after pHighestRight, that has the same y than the highest vertex
    TVec2fVect::const_iterator pHighestLeft = Next(pHighestRight, V);
    // seek the vertex
    while (pHighestLeft->y != fHighest)
    {
        pHighestLeft = Next(pHighestLeft, V);
    }
    TVec2fVect::const_iterator pPrevHighestLeft = Prev(pHighestLeft, V);
    // we need to get the orientation of the polygon
    // There are 2 case : flat, and non-flat top
    // check for flat top
    if (pHighestLeft->x != pHighestRight->x)
    {
        // compare right and left side
        return pHighestLeft->x <= pHighestRight->x;
    }
    // The top of the poly is sharp
    // We perform a cross product of the 2 highest vect to get its orientation
     float deltaXN = Next(pHighestRight, V)->x - pHighestRight->x;
     float deltaYN = Next(pHighestRight, V)->y - pHighestRight->y;
     float deltaXP = pPrevHighestLeft->x - pHighestLeft->x;
     float deltaYP = pPrevHighestLeft->y - pHighestLeft->y;
     return (deltaXN * deltaYP - deltaYN * deltaXP) >= 0;
}

// *******************************************************************************
void    CPolygon2D::computeBorders(TRasterVect &borders, sint &highestY) const
{
    #ifdef NL_DEBUG
        checkValidBorders();
    #endif
    // an 'alias' to the vertices
    const TVec2fVect &V = Vertices;
    if (Vertices.size() < 3)
    {
        borders.clear();
        return;
    }
    bool    ccw;  // set to true when it has a counter clock wise orientation

    // compute highest and lowest pos of the poly
    float fHighest = V[0].y;
    float fLowest  = fHighest;

    // iterators to the thighest and lowest vertex
    TVec2fVect::const_iterator pLowest = V.begin(), pHighest = V.begin();
    TVec2fVect::const_iterator it = V.begin() ;
    const TVec2fVect::const_iterator endIt = V.end();
    do
    {
        if (it->y > fLowest)
        {
            fLowest = it->y;
            pLowest = it;
        }
        else
        if (it->y < fHighest)
        {
            fHighest = it->y;
            pHighest = it;
        }
        ++it;
    }
    while (it != endIt);


    sint iHighest = (sint) ceilf(fHighest) ;
    sint iLowest  = (sint) ceilf(fLowest) ;

    highestY = iHighest;


    uint polyHeight = iLowest - iHighest;
    if (polyHeight <= 0)
    {
        borders.clear();
        return;
    }

    borders.resize(polyHeight);

    // iterator to the first vertex that has an y different from the top vertex
    TVec2fVect::const_iterator pHighestRight = pHighest;
    // we seek this vertex
    while (Next(pHighestRight, V)->y == fHighest)
    {
        pHighestRight = Next(pHighestRight, V);
    }

    // iterator to the first vertex after pHighestRight, that has the same y than the highest vertex
    TVec2fVect::const_iterator pHighestLeft = Next(pHighestRight, V);
    // seek the vertex
    while (pHighestLeft->y != fHighest)
    {
        pHighestLeft = Next(pHighestLeft, V);
    }

    TVec2fVect::const_iterator pPrevHighestLeft = Prev(pHighestLeft, V);

    // we need to get the orientation of the polygon
    // There are 2 case : flat, and non-flat top

    // check for flat top
    if (pHighestLeft->x != pHighestRight->x)
    {
        // compare right and left side
        if (pHighestLeft->x > pHighestRight->x)
        {
            ccw = true;  // the list is CCW oriented
            std::swap(pHighestLeft, pHighestRight);
        }
        else
        {
            ccw = false; // the list is CW oriented
        }
    }
    else
    {
        // The top of the poly is sharp
        // We perform a cross product of the 2 highest vect to get its orientation

         const float deltaXN = Next(pHighestRight, V)->x - pHighestRight->x;
         const float deltaYN = Next(pHighestRight, V)->y - pHighestRight->y;
         const float deltaXP = pPrevHighestLeft->x - pHighestLeft->x;
         const float deltaYP = pPrevHighestLeft->y - pHighestLeft->y;
         if ((deltaXN * deltaYP - deltaYN * deltaXP) < 0)
         {
            ccw = true;  // the list is CCW oriented
            std::swap(pHighestLeft, pHighestRight);
         }
         else
         {
            ccw = false; // the list is CW oriented
         }
    }


    // compute borders
    TVec2fVect::const_iterator currV, nextV; // current and next vertex
    if (!ccw) // clock wise order ?
    {
        currV = pHighestRight ;
        // compute right edge from top to bottom
        do
        {
            nextV = Next(currV, V);
            ScanEdge(borders, iHighest, *currV, *nextV, true);
            currV = nextV;
        }
        while (currV != pLowest); // repeat until we reach the bottom vertex

        // compute left edge from bottom to top
        do
        {
            nextV = Next(currV, V);
            ScanEdge(borders, iHighest, *nextV, *currV, false);
            currV = nextV;
        }
        while (currV != pHighestLeft);
    }
    else // ccw order
    {
        currV = pHighestLeft;
        // compute left edge from top to bottom
        do
        {
            nextV = Next(currV, V);
            ScanEdge(borders, iHighest, *currV, *nextV, false) ;
            currV = nextV;
        }
        while (currV != pLowest) ;

        // compute right edge from bottom to top
        do
        {
            nextV = Next(currV, V);
            ScanEdge(borders, iHighest, *nextV, *currV, true);
            currV = nextV;
        }
        while (currV != pHighestRight)  ;
    }
}

//=========================================================================
// scan outer right edge of a poly
static void ScanOuterEdgeRight(CPolygon2D::TRaster *r, float x1, float y1, float x2, float y2, sint minY)
{
    CPolygon2D::TRaster *currRaster;
    float deltaX, deltaY;
    float inverseSlope;
    sint32  iInverseSlope, iposx;
    sint  height;
    deltaX = x2 - x1;
    height = (sint) (ceilf(y2) - floorf(y1)) ;
    if (height <= 0) return;
    if (deltaX >= 0.f)
    {
        if (height == 1)
        {
            currRaster = r + ((sint) floorf(y1) - minY);
            currRaster->second = std::max((sint) floorf(x2), currRaster->second);
        }
        else
        {
            deltaY = y2 - y1;
            if(deltaY)
                inverseSlope = deltaX / deltaY;
            else
                inverseSlope = 0;
            iInverseSlope = (sint32) (65536.0 * inverseSlope);
            currRaster = r + ((sint) floorf(y1) - minY);
            iposx = (sint32) (65536.0 * (x1 + inverseSlope * (ceilf(y1) - y1))); // sub-pixel accuracy
            if (ceilf(y1) == y1)
            {
                iposx += iInverseSlope;
            }
            do
            {
                currRaster->second = std::max((sint) (iposx >> 16), currRaster->second);
                iposx += iInverseSlope;
                ++ currRaster;
                -- height;
            }
            while (height != 1);
            // correction for last line
            currRaster->second = std::max((sint) floorf(x2), currRaster->second);
        }
    }
    else
    {
        deltaY = y2 - y1;
        if(deltaY)
            inverseSlope = deltaX / deltaY;
        else
            inverseSlope = 0;
        iInverseSlope = (sint32) (65536.0 * inverseSlope);
        currRaster = r + ((sint) floorf(y1) - minY);
        currRaster->second = std::max((sint) floorf(x1), currRaster->second);
        ++ currRaster;
        iposx = (sint32) (65536.0 * (x1 + inverseSlope * (ceilf(y1) - y1))); // sub-pixel accuracy
        if (ceilf(y1) == y1)
        {
            iposx += iInverseSlope;
        }
        while (--height)
        {
            currRaster->second = std::max((sint) (iposx >> 16), currRaster->second);
            iposx += iInverseSlope;
            ++ currRaster;
        }
    }
}

//=========================================================================
// scan outer left edge of a poly
static void ScanOuterEdgeLeft(CPolygon2D::TRaster *r, float x1, float y1, float x2, float y2, sint minY)
{
    CPolygon2D::TRaster *currRaster;
    float deltaX, deltaY;
    float inverseSlope;
    sint32   iInverseSlope, iposx;
    sint  height;
    deltaX = x2 - x1;
    height = (sint) (ceilf(y2) - floorf(y1)) ;
    if (height <= 0) return;
    if (deltaX < 0.f)
    {
        if (height == 1)
        {
            currRaster = r + ((sint) floorf(y1) - minY);
            currRaster->first = std::min((sint) floorf(x2), currRaster->first);
        }
        else
        {
            deltaY = y2 - y1;
            if(deltaY)
                inverseSlope = deltaX / deltaY;
            else
                inverseSlope = 0;
            iInverseSlope = (sint32) (65536.0 * inverseSlope);
            currRaster = r + ((sint) floorf(y1) - minY);
            iposx = (sint32) (65536.0 * (x1 + inverseSlope * (ceilf(y1) - y1))); // sub-pixel accuracy
            if (ceilf(y1) == y1)
            {
                iposx += iInverseSlope;
            }
            do
            {
                currRaster->first = std::min((sint) (iposx >> 16), currRaster->first);
                iposx += iInverseSlope;
                ++ currRaster;
                -- height;
            }
            while (height != 1);
            // correction for last line
            currRaster->first = std::min((sint) floorf(x2), currRaster->first);
        }
    }
    else
    {
        deltaY = y2 - y1;
        if(deltaY)
            inverseSlope = deltaX / deltaY;
        else
            inverseSlope = 0;
        iInverseSlope = (sint32) (65536.0 * inverseSlope);
        currRaster = r + ((sint) floorf(y1) - minY);
        currRaster->first = std::min((sint) floorf(x1), currRaster->first);
        ++ currRaster;
        iposx = (sint32) (65536.0 * (x1 + inverseSlope * (ceilf(y1) - y1))); // sub-pixel accuracy
        if (ceilf(y1) == y1)
        {
            iposx += iInverseSlope;
        }
        while (--height)
        {
            currRaster->first = std::min((sint) (iposx >> 16), currRaster->first);
            iposx += iInverseSlope;
            ++ currRaster;
        }
    }
}


// *******************************************************************************
void CPolygon2D::computeOuterBorders(TRasterVect &borders, sint &minimumY) const
{
    #ifdef NL_DEBUG
        checkValidBorders();
    #endif
    borders.clear();
    // NB : this version is not much optimized, because of the min/max test
    // during rasterization.
    // TODO : optimize if needed ...

    if (Vertices.empty())
    {
        minimumY = -1;
        return;
    }
    const CVector2f *first = &Vertices[0];
    const CVector2f *last  = first + Vertices.size();

    const CVector2f *curr = first, *next, *plowest ,*phighest;
    const CVector2f *pHighestRight, *pHighestRightNext, *pHighestLeft;
    const CVector2f *pPrevHighestLeft;
    double          deltaXN, deltaYN, deltaXP, deltaYP;
    bool            ccw;  // true if CCW oriented
    sint            polyHeight;
    sint            highest, lowest;

    float fright   = curr->x;
    float fleft    = curr->x;
    float fhighest = curr->y;
    float flowest  = curr->y;
    plowest = phighest = curr;

    // compute highest and lowest pos of the poly
    do
    {
        fright = std::max(fright, curr->x);
        fleft  = std::min(fleft, curr->x);
        if (curr->y > flowest)
        {
            flowest = curr->y;
            plowest = curr;
        }
        if (curr->y < fhighest)
        {
            fhighest = curr->y;
            phighest = curr;
        }
        ++curr;
    }
    while (curr != last);


    highest = (sint) floorf(fhighest);
    lowest = (sint) floorf(flowest);

    polyHeight = lowest - highest + 1;
    nlassert(polyHeight > 0);

    // make room for rasters
    borders.resize(polyHeight);
    // fill with xmin / xman
    sint ileft = (sint) floorf(fleft);
    sint iright = (sint) ceilf(fright);
    minimumY = highest;
    if (flowest == fhighest) // special case : degenerate poly
    {

        borders.resize(1);
        borders.front().first =  ileft;
        borders.front().second =  ileft;
        return;
    }
    //
    for(TRasterVect::iterator it = borders.begin(); it != borders.end(); ++it)
    {
        it->second = ileft;
        it->first = iright;
    }



    pHighestRight = phighest;
    for (;;)
    {
        pHighestRightNext  = pHighestRight + 1;
        if (pHighestRightNext == last) pHighestRightNext = first;
        if (pHighestRightNext->y != pHighestRight->y) break;
        pHighestRight = pHighestRightNext;
    }

    pPrevHighestLeft = pHighestRight;
    pHighestLeft = pHighestRight;
    ++pHighestLeft;
    if (pHighestLeft == last) pHighestLeft = first;

    while (pHighestLeft->y != fhighest)
    {
        pPrevHighestLeft = pHighestLeft;
        ++pHighestLeft;
        if (pHighestLeft == last) pHighestLeft = first;
    }


    // we need to get the orientation of the polygon
    // There are 2 case : flat, and non-flat top

    // check for flat top
    if (pHighestLeft->x != pHighestRight->x)
    {
        // compare right and left side
        if (pHighestLeft->x > pHighestRight->x)
        {
            ccw = true;  // the list is CCW oriented
            std::swap(pHighestLeft, pHighestRight);
        }
        else
        {
            ccw = false; // the list is CW oriented
        }
    }
    else
    {
        pHighestRightNext = pHighestRight + 1;
        if (pHighestRightNext == last) pHighestRightNext = first;
         deltaXN = pHighestRightNext->x - pHighestRight->x;
         deltaYN = pHighestRightNext->y - pHighestRight->y;
         deltaXP = pPrevHighestLeft->x - pHighestLeft->x;
         deltaYP = pPrevHighestLeft->y - pHighestLeft->y;
         if ((deltaXN * deltaYP - deltaYN * deltaXP) < 0)
         {
            ccw = true;
            std::swap(pHighestLeft, pHighestRight);
         }
         else
         {
            ccw = false;
         }
    }

    if (!ccw)
    {
        // clock wise oriented list
        curr = pHighestRight;
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanOuterEdgeRight(&borders[0], curr->x, curr->y, next->x, next->y, minimumY);
            curr = next;
        }
        while (curr != plowest);
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanOuterEdgeLeft(&borders[0], next->x, next->y, curr->x, curr->y, minimumY);
            curr = next;
        }
        while (curr != pHighestLeft);
    }
    else
    {
        // ccw oriented
        curr = pHighestLeft;
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanOuterEdgeLeft(&borders[0], curr->x, curr->y, next->x, next->y, minimumY);
            curr = next;
        }
        while (curr != plowest);
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanOuterEdgeRight(&borders[0], next->x, next->y, curr->x, curr->y, minimumY);
            curr = next;
        }
        while (curr != pHighestRight);
    }
}


//=========================================================================
// scan inner right edge of a poly
static void ScanInnerEdge(CPolygon2D::TRaster *r, float x1, float y1, float x2, float y2, sint minY, bool rightEdge)
{
    const uint rol16 = 65536;
    CPolygon2D::TRaster *currRaster;
    float deltaX, deltaY;
    float inverseSlope;
    sint32  iInverseSlope, iposx;
    sint  height;
    deltaX = x2 - x1;
    height = (sint) (ceilf(y2) - floorf(y1));
    if (height <= 0) return;
    deltaY = y2 - y1;
    if(deltaY)
        inverseSlope = deltaX / deltaY;
    else
        inverseSlope = 0;
    iInverseSlope = (sint32) (rol16 * inverseSlope);
    currRaster = r + ((sint) floorf(y1) - minY);
    //
    iposx = (sint32) (rol16 * (x1 + inverseSlope * (ceilf(y1) - y1))); // sub-pixel accuracy
    if (rightEdge)
    {
        iposx -= rol16 - 1;
        if (deltaX >= 0.f)
        {
            // start of segment
            if (floorf(y1) != y1)
            {
                currRaster->second = std::min((sint) floorf(x1) - 1, currRaster->second);
                ++ currRaster;
                -- height;
                if (height == 0) return;
            }
            do
            {
                currRaster->second = std::min((sint) (iposx >> 16), currRaster->second);
                iposx += iInverseSlope;
                ++ currRaster;
            }
            while (--height);
        }
        else
        {
            // start of segment
            if (floorf(y1) != y1)
            {
                currRaster->second = std::min((sint) (iposx >> 16), currRaster->second);
                ++ currRaster;
                -- height;
                if (height == 0) return;
            }
            while (--height)
            {
                iposx += iInverseSlope;
                currRaster->second = std::min((sint) (iposx >> 16), currRaster->second);
                ++ currRaster;
            }
            // fill bottom of segment
            currRaster->second = std::min((sint) floorf(x2) - 1, currRaster->second);
        }
    }
    else
    {
        iposx += rol16 - 1;
        if (deltaX < 0.f)
        {
            // start of segment
            if (floorf(y1) != y1)
            {
                currRaster->first = std::max((sint) ceilf(x1), currRaster->first);
                ++ currRaster;
                -- height;
                if (height == 0) return;
            }
            do
            {
                currRaster->first = std::max((sint) (iposx >> 16), currRaster->first);
                iposx += iInverseSlope;
                ++ currRaster;
            }
            while (--height);
        }
        else
        {
            // start of segment
            if (floorf(y1) != y1)
            {
                currRaster->first = std::max((sint) (iposx >> 16), currRaster->first);
                ++ currRaster;
                -- height;
                if (height == 0) return;
            }
            while (--height)
            {
                iposx += iInverseSlope;
                currRaster->first = std::max((sint) (iposx >> 16), currRaster->first);
                ++ currRaster;
            }
            // fill bottom of segment
            currRaster->first = std::max((sint) ceilf(x2), currRaster->first);
        }
    }
}

// *******************************************************************************
void CPolygon2D::computeInnerBorders(TRasterVect &borders, sint &minimumY) const
{
    #ifdef NL_DEBUG
        checkValidBorders();
    #endif
    borders.clear();
    if (Vertices.empty())
    {
        minimumY = -1;
        return;
    }
    const CVector2f *first = &Vertices[0];
    const CVector2f *last  = first + Vertices.size();

    const CVector2f *curr = first, *next, *plowest ,*phighest;
    const CVector2f *pHighestRight, *pHighestRightNext, *pHighestLeft;
    const CVector2f *pPrevHighestLeft;
    double          deltaXN, deltaYN, deltaXP, deltaYP;
    bool            ccw;  // true if CCW oriented
    sint            polyHeight;
    sint            highest, lowest;

    float fright   = curr->x;
    float fleft    = curr->x;
    float fhighest = curr->y;
    float flowest  = curr->y;
    plowest = phighest = curr;

    // compute highest (with lowest y) and lowest (with highest y) points of the poly
    do
    {
        fright = std::max(fright, curr->x);
        fleft  = std::min(fleft, curr->x);
        if (curr->y > flowest)
        {
            flowest = curr->y;
            plowest = curr;
        }
        if (curr->y < fhighest)
        {
            fhighest = curr->y;
            phighest = curr;
        }
        ++curr;
    }
    while (curr != last);
    if (flowest == fhighest)
    {
        minimumY = -1;
        return;
    }
    highest = (sint) floorf(fhighest);
    lowest = (sint) ceilf(flowest);

    polyHeight = lowest - highest;
    minimumY = highest;
    if (polyHeight == 0)
    {
        minimumY = -1;
        return;
    }
    // make room for rasters
    borders.resize(polyHeight);
    // fill with xmin / xman
    sint ileft = (sint) floorf(fleft) - 1;
    sint iright = (sint) ceilf(fright);
    for(TRasterVect::iterator it = borders.begin(); it != borders.end(); ++it)
    {
        it->second = iright;
        it->first = ileft;
    }
    pHighestRight = phighest;
    for (;;)
    {
        pHighestRightNext  = pHighestRight + 1;
        if (pHighestRightNext == last) pHighestRightNext = first;
        if (pHighestRightNext->y != pHighestRight->y) break;
        pHighestRight = pHighestRightNext;
    }

    pPrevHighestLeft = pHighestRight;
    pHighestLeft = pHighestRight;
    ++pHighestLeft;
    if (pHighestLeft == last) pHighestLeft = first;

    while (pHighestLeft->y != fhighest)
    {
        pPrevHighestLeft = pHighestLeft;
        ++pHighestLeft;
        if (pHighestLeft == last) pHighestLeft = first;
    }


    // we need to get the orientation of the polygon
    // There are 2 case : flat, and non-flat top

    // check for flat top
    if (pHighestLeft->x != pHighestRight->x)
    {
        // compare right and left side
        if (pHighestLeft->x > pHighestRight->x)
        {
            ccw = true;  // the list is CCW oriented
            std::swap(pHighestLeft, pHighestRight);
        }
        else
        {
            ccw = false; // the list is CW oriented
        }
    }
    else
    {
        pHighestRightNext = pHighestRight + 1;
        if (pHighestRightNext == last) pHighestRightNext = first;
         deltaXN = pHighestRightNext->x - pHighestRight->x;
         deltaYN = pHighestRightNext->y - pHighestRight->y;
         deltaXP = pPrevHighestLeft->x - pHighestLeft->x;
         deltaYP = pPrevHighestLeft->y - pHighestLeft->y;
         if ((deltaXN * deltaYP - deltaYN * deltaXP) < 0)
         {
            ccw = true;
            std::swap(pHighestLeft, pHighestRight);
         }
         else
         {
            ccw = false;
         }
    }

    if (!ccw)
    {
        // cw oriented
        curr = pHighestRight;
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanInnerEdge(&borders[0], curr->x, curr->y, next->x, next->y, minimumY, true);
            curr = next;
        }
        while (curr != plowest);
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanInnerEdge(&borders[0], next->x, next->y, curr->x, curr->y, minimumY, false);
            curr = next;
        }
        while (curr != pHighestLeft);
    }
    else
    {
        // ccw oriented
        curr = pHighestLeft;
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanInnerEdge(&borders[0], curr->x, curr->y, next->x, next->y, minimumY, false);
            curr = next;
        }
        while (curr != plowest);
        do
        {
            next = curr + 1;
            if (next == last) next = first;
            ScanInnerEdge(&borders[0], next->x, next->y, curr->x, curr->y, minimumY, true);
            curr = next;
        }
        while (curr != pHighestRight);
    }
    // fix for top
    if (floorf(fhighest) != fhighest)
    {
        borders[0].first = 1;
        borders[0].second = 0;
    }
    // fix for bottom
    if (floorf(flowest) != flowest)
    {
        borders.back().first = 1;
        borders.back().second = 0;
    }
}

// *******************************************************************************
void CPolygon2D::checkValidBorders() const
{
    for (uint k = 0; k < Vertices.size(); ++k)
    {
        nlassert(Vertices[k].x >= -32000.f); // coordinate too big !
        nlassert(Vertices[k].x < 32000.f);   // coordinate too big !
        nlassert(Vertices[k].y >= -32000.f); // coordinate too big !
        nlassert(Vertices[k].y < 32000.f);   // coordinate too big !
    }
}

// *******************************************************************************
float       CPolygon2D::sumDPAgainstLine(float a, float b, float c) const
{
    float sum = 0.f;
    for (uint k = 0; k < Vertices.size(); ++k)
    {
        const CVector2f &p = Vertices[k];
        sum += a * p.x + b * p.y + c;
    }
    return sum;
}


// *******************************************************************************
bool  CPolygon2D::getNonNullSeg(uint &index) const
{
    nlassert(Vertices.size() > 0);
    float bestLength = 0.f;
    sint  bestIndex = -1;
    for (uint k = 0; k < Vertices.size() - 1; ++k)
    {
        float norm2 = (Vertices[k + 1] - Vertices[k]).sqrnorm();
        if ( norm2 > bestLength)
        {
            bestLength = norm2;
            bestIndex = (int) k;
        }
    }
    float norm2 = (Vertices[Vertices.size() - 1] - Vertices[0]).sqrnorm();
    if ( norm2 > bestLength)
    {
        index = Vertices.size() - 1;
        return true;
    }

    if (bestIndex != -1)
    {
        index = bestIndex;
        return true;
    }
    else
    {
        return false;
    }
}


// *******************************************************************************
void  CPolygon2D::getLineEquation(uint index, float &a, float &b, float &c) const
{
    nlassert(index < Vertices.size());
    const CVector2f &v0 = getSegRef0(index);
    const CVector2f &v1 = getSegRef1(index);

    NLMISC::CVector2f seg = v0 - v1;
    a = seg.y;
    b = - seg.x;
    c = - v0.x * a - v0.y * b;
}

// *******************************************************************************
bool        CPolygon2D::intersect(const CPolygon2D &other) const
{
    nlassert(other.Vertices.size() > 0);
    uint nonNullSegIndex;
    if (getNonNullSeg(nonNullSegIndex))
    {
        float a0, b0, c0; 
        getLineEquation(nonNullSegIndex, a0, b0, c0);
        float orient = sumDPAgainstLine(a0, b0, c0);

        for (uint k = 0; k < Vertices.size(); ++k)
        {
            if ( (getSegRef0(k) - getSegRef1(k)).sqrnorm() == 0.f) continue;

            float a, b, c; 
            getLineEquation(k, a, b, c);
            uint l;
            for (l = 0; l < other.Vertices.size(); ++l)
            {
                const CVector2f &ov = other.Vertices[l];
                if ( orient * (ov.x * a + ov.y * b +c) > 0.f) break;
            }
            if (l == other.Vertices.size()) // all point on the outside
            {
                return false; // outside
            }
        }
        return true;
    }
    else // this poly is just a point
    {
        return other.contains(Vertices[0]);
    }
}

// *******************************************************************************
bool        CPolygon2D::contains(const CVector2f &p, bool hintIsConvex /*= true*/) const
{
    if (hintIsConvex)
    {
        uint numVerts = Vertices.size();
        nlassert(numVerts >= 0.f);
        for (uint k = 0; k < numVerts; ++k)
        {
            if (getSegRef0(k) != getSegRef1(k))
            {
                float a, b, c; 
                getLineEquation(k, a, b, c);
                float orient = a * p.x + b * p.y + c;
                for(uint l = k + 1; l < numVerts; ++l)
                {
                    getLineEquation(l, a, b, c);
                    float newOrient = a * p.x + b * p.y + c;
                    if (newOrient * orient < 0.f) return false;
                }
                return true;
            }
        }
        // the poly reduces to a point
        return p == Vertices[0];
    }
    else
    {
        // concave case
        static std::vector<float> xInter;
        xInter.clear();
        for(uint k = 0; k < Vertices.size(); ++k)
        {
            const CVector2f &p0 = getSegRef0(k);
            const CVector2f &p1 = getSegRef1(k);
            if (p0.y == p1.y)
            {
                if (p.y == p0.y)
                {
                    if ((p.x >= p0.x && p.x <= p1.x)
                        || (p.x >= p1.x && p.x <= p0.x))
                    {
                        return true;
                    }
                }
            }
            if ((p.y >= p0.y && p.y < p1.y) ||
                (p.y >= p1.y && p.y < p0.y)
               )
            {

                float inter = p0.x + (p.y - p0.y) * (p1.x - p0.x) / (p1.y- p0.y);
                xInter.push_back(inter);
            }
        }
        if (xInter.size() < 2) return false;
        std::sort(xInter.begin(), xInter.end());
        for(uint k = 0; k < xInter.size() - 1; ++k)
        {
            if (p.x >= xInter[k] && p.x <= xInter[k + 1])
            {
                return (k & 1) == 0;
            }
        }
        return false;
    }
}


// *******************************************************************************
CPolygon2D::CPolygon2D(const CTriangle &tri, const CMatrix &projMat)
{
    Vertices.resize(3);
    NLMISC::CVector proj[3] = { projMat * tri.V0, projMat * tri.V1, projMat * tri.V2 };
    Vertices[0].set(proj[0].x, proj[0].y);
    Vertices[1].set(proj[1].x, proj[1].y);
    Vertices[2].set(proj[2].x, proj[2].y);
}

// *******************************************************************************
void CPolygon2D::getBoundingRect(CVector2f &minCorner, CVector2f &maxCorner) const
{
    nlassert(!Vertices.empty())
    minCorner = maxCorner = Vertices[0];
    uint numVertices = Vertices.size();
    for(uint k = 0; k < numVertices; ++k)
    {
        minCorner.minof(minCorner, Vertices[k]);
        maxCorner.maxof(minCorner, Vertices[k]);
    }
}

// *******************************************************************************
bool operator ==(const CPolygon2D &lhs,const CPolygon2D &rhs)
{
    if (lhs.Vertices.size() != rhs.Vertices.size()) return false;
    return std::equal(lhs.Vertices.begin(), lhs.Vertices.end(), rhs.Vertices.begin());
}

// *******************************************************************************
bool operator < (const CPolygon2D &lhs, const CPolygon2D &rhs)
{
    if (lhs.Vertices.size() != rhs.Vertices.size()) return lhs.Vertices.size() < rhs.Vertices.size();
    for(uint k = 0; k < lhs.Vertices.size(); ++k)
    {
        if (lhs.Vertices[k] != rhs.Vertices[k]) return lhs.Vertices[k] < rhs.Vertices[k];
    }
    return false;
}


// *******************************************************************************
static inline bool testSegmentIntersection(const CVector2f &a, const CVector2f &b,
                                           const CVector2f &c, const CVector2f &d)
{
    double denom = a.x * double(d.y - c.y) +
            b.x * double(c.y - d.y) +
            d.x * double(b.y - a.y) +
            c.x * double(a.y - b.y);

    if (denom == 0) return false;
    //
    double num = a.x * double(d.y - c.y) +
                 c.x * double(a.y - d.y) +
                 d.x * double(c.y - a.y);

    if (num == 0 || (num == denom)) return false;
    double lambda = num / denom;
    if (lambda <= 0 || lambda >= 1) return false;
    //
    num = - (a.x * double(c.y - b.y) +
          b.x * double(a.y - c.y) +
          c.x * double(b.y - a.y));

    if (num == 0 || (num == denom)) return false;
    lambda = num / denom;
    if (lambda <= 0 || lambda >= 1) return false;
    return true;
}


// *******************************************************************************
bool CPolygon2D::selfIntersect() const
{
    if (Vertices.size() < 3) return false;
    uint numEdges = Vertices.size();
    for(uint k = 0; k < numEdges; ++k)
    {
        // test intersection with all other edges that don't share a vertex with this one
        const CVector2f &p0 = getSegRef0(k);
        const CVector2f &p1 = getSegRef1(k);
        for(uint l = 0; l < k; ++l)
        {
            const CVector2f &v0 = getSegRef0(l);
            const CVector2f &v1 = getSegRef1(l);
            if (v0 == p0 || v0 == p1 || v1 == p0 || v1 == p1) continue;
            //
            if (testSegmentIntersection(p0, p1, v0, v1)) return true;
        }
    }
    return false;
}



} // NLMISC