UT_Convex.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME:        UT_Convex.h ( UT Library, C++)
 *
 * COMMENTS:    This code does general polygon convexing.  It's down here
 *              in UT since there are places other than the GEO lib which
 *              may need to use convexing code (i.e. the RE lib, or UI).
 *              This code convexes 2D polygons.
 */

#ifndef __UT_Convex__
#define __UT_Convex__

#include "UT_API.h"
#include "UT_NonCopyable.h"
#include "UT_Vector3.h"
#include "UT_Vector2.h"
#include "UT_Array.h"
#include <SYS/SYS_Types.h>
#include <stdlib.h>

class UT_ConvexPoint;
class UT_Interrupt;

class UT_API UT_Convex {
public:
             UT_Convex()
             {
                 myData = 0;
                 myConvexInfo = 0;
                 myConvexSize = 0;
                 myConvexEntries = 0;
                 myPointList = 0;
                 myPointSize = 0;
             }
    virtual ~UT_Convex()
             {
                 if (myConvexInfo) free(myConvexInfo);
                 if (myPointList) free(myPointList);
             }

             UT_NON_COPYABLE(UT_Convex);

    bool                isConcave();
    void                convex(exint maxsides, bool allowinterrupt = true,
                               bool avoiddegeneracy = false);

    virtual void        initialize();
    /// It's the sub-class responsibility to fill these out...
    virtual int         getPointCount() const = 0;
    virtual void        getPoint(int num, float &x, float &y) const = 0;

    /// beginPolygon will be told when the last polygon is being added
    virtual void        beginPolygon(int npoints, bool lastone) = 0;
    virtual void        polyVertex(int num) = 0;
    virtual void        endPolygon() = 0;

    /// If the last polygon to be added won't be added after all, e.g. if it is
    /// degenerate and avoiddegeneracy is true, the earlyTermination() method
    /// will be called.  This is so that the caller can clean up properly, since
    /// beginPolygon will *not* have been called with lastone==true.
    /// By default, this method does nothing.
    virtual void        earlyTermination();

    /// Returns true if splitting the specified quad along the line from
    /// point 0 to point 2 is possible
    static inline bool canSplitQuad02(const UT_Vector3 &p0,
                                      const UT_Vector3 &p1,
                                      const UT_Vector3 &p2,
                                      const UT_Vector3 &p3)
    {
        const UT_Vector3 v0 = p1 - p0;
        const UT_Vector3 v1 = p2 - p1;
        const UT_Vector3 v2 = p3 - p2;
        const UT_Vector3 v3 = p0 - p3;
        // Okay if normals at 1 and 3 point in the same direction
        return ((cross(v0, v1)).dot(cross(v2, v3)) > 0);
    }

    /// Returns true if the specified quad is convex.
    /// NOTE: This may return false negatives for highly non-planar quads,
    ///       since the definition of convex breaks down.
    static inline bool isQuadConvex(const UT_Vector3 &p0,
                                    const UT_Vector3 &p1,
                                    const UT_Vector3 &p2,
                                    const UT_Vector3 &p3)
    {
        const UT_Vector3 v0 = p1 - p0;
        const UT_Vector3 v1 = p2 - p1;
        const UT_Vector3 v2 = p3 - p2;
        const UT_Vector3 v3 = p0 - p3;
        // Convex if normals at 1 and 3 point in the same direction
        // and normals at 0 and 2 point in the same direction
        return (cross(v0, v1).dot(cross(v2, v3)) > 0)
            && (cross(v1, v2).dot(cross(v3, v0)) > 0);
    }

    /// Relative size below which a triangle or diagonal counts as degenerate.
    /// A template so that it is always inlined, never imported from the DLL.
    template <typename T = fpreal>
    static constexpr T smallAngleTolerance() { return T(1e-6); }

    /// Returns true if, for triangles (a,b,c) and (c,b,d) sharing edge b-c,
    /// the diagonal a-d gives the larger smallest angle. This is the "Avoid
    /// Small Angles" test of SOP_Divide and handles degenerate input.
    /*
     *       D                  D
     *      / \                /|\
     *     / 1 \              / | \
     *    /     \            /  |  \
     *   C-------B flips to C 3 | 2 B
     *    \     /            \  |  /
     *     \ 0 /              \ | /
     *      \ /                \|/
     *       A                  A
     */
    template <typename T>
    static bool shouldFlipEdge(const UT_Vector2T<T> &a,
                               const UT_Vector2T<T> &b,
                               const UT_Vector2T<T> &c,
                               const UT_Vector2T<T> &d)
    {
        const fpreal tol = smallAngleTolerance();

        // First check if a flip is even possible
        // If sign(ca X cd) != sign(bd x ba), can't flip
        UT_Vector2T<T> ca = a-c;
        UT_Vector2T<T> cd = d-c;
        UT_Vector2T<T> bd = d-b;
        UT_Vector2T<T> ba = a-b;
        fpreal ccross = cross(ca, cd);
        fpreal bcross = cross(bd, ba);
        fpreal lca = ca.length();
        fpreal lcd = cd.length();
        fpreal lbd = bd.length();
        fpreal lba = ba.length();
        fpreal cscale = lca * lcd;
        fpreal bscale = lbd * lba;
        // If c is definitely negative and b is definitely positive,
        // or c is definitely positive and b is definitely negative,
        // can't flip, since one angle is reflex.
        if ((ccross < -tol * cscale && bcross >  tol * bscale) ||
            (ccross >  tol * cscale && bcross < -tol * bscale))
            return false;

        fpreal perim = lca + lcd + lbd + lba;
        fpreal lsmall = tol * perim;

        // Compute diagonals
        UT_Vector2T<T> ad = d-a;
        UT_Vector2T<T> bc = c-b;
        fpreal lad = ad.length();
        fpreal lbc = bc.length();

        // Handle near-degenerate diagonals by
        // choosing the zero-length diagonal.
        // Having the = handles the case of exactly zero perimeter.
        if (lbc <= lsmall)
            return false;
        if (lad <= lsmall)
            return true;

        // In all other degenerate point cases, flipping will
        // result in just as bad a triangle, so don't bother.
        if (SYSmin(lba, lca, lbd, lcd) <= lsmall)
            return false;

        // Don't need to consider angles CAB or BDC,
        // because the condition is always determined by
        // ABC, CBD, DCB, BCA; CAD, DAB, BDA, ADC
        // Since these angles are always <= 90 degrees, (or close enough
        // in the case of roundoff with near-degeneracies), the (absolute)
        // sine of the angles suffices for this comparison.
        // Will use: |sin(XYZ)| = 2*area(XYZ)/(|XY|*|YZ|)
        fpreal area2x2overad = SYSabs(bcross) / lad;
        fpreal area3x2overad = SYSabs(ccross) / lad;
        // Compute both crosses for the other two triangles
        // One could be eliminated subtracting from the total area,
        // but need absolute areas even if negative area from roundoff.
        fpreal across = cross(ba, ca); // same as ab X ac
        fpreal dcross = cross(cd, bd); // same as dc X db
        fpreal area0x2overbc = SYSabs(across) / lbc;
        fpreal area1x2overbc = SYSabs(dcross) / lbc;

        fpreal min01 = SYSmin(area0x2overbc / lba,
                              area0x2overbc / lca,
                              area1x2overbc / lbd,
                              area1x2overbc / lcd);
        fpreal min23 = SYSmin(area2x2overad / lba,
                              area2x2overad / lbd,
                              area3x2overad / lca,
                              area3x2overad / lcd);
        return min23 > min01;
    }

    /// Flips edges of the triangulation tripts (three corner indices per
    /// triangle) until the smallest angles are maximized. pts(i) is the 2D
    /// position of corner i. on_flip(tri, tri2, i, j) is called before each
    /// flip of tri's edge i, shared with tri2's edge j.
    template <typename PTS, typename TRIS, typename ON_FLIP>
    static void flipEdgesToAvoidSmallAngles(const PTS &pts, TRIS &tripts,
                                            ON_FLIP &&on_flip)
    {
        const int ntripts = int(tripts.size());
        const int ntris = ntripts/3;

        // For a triangle ABC, the adjacent triangles are in the
        // order of the triangles on the opposite side of:
        // AB, then BC, then CA. -1 if no opposite triangle
        // The integer for a triangle's edge is the index
        // of the other triangle's edge, which points back.
        UT_Array<int> adjtrioff;
        adjtrioff.setSize(ntripts);
        for (int i = 0; i < ntripts; ++i)
            adjtrioff(i) = -1;
        for (int tri = 0; tri < ntris; ++tri)
        {
            int trioff = 3 * tri;
            for (int i = 0; i < 3; ++i)
            {
                // Don't revisit edges
                if (adjtrioff(trioff + i) >= 0)
                    continue;
                // Edge points are i and (i+1)%3
                int pti   = tripts(trioff + i);
                int ptip1 = tripts(trioff + (i + 1) % 3);
                bool found = false;
                // Search all following triangles' edges
                for (int tri2 = tri+1; tri2 < ntris && !found; ++tri2)
                {
                    int tri2off = 3 * tri2;
                    for (int j = 0; j < 3; ++j)
                    {
                        int ptj   = tripts(tri2off + j);
                        int ptjp1 = tripts(tri2off + (j + 1) % 3);
                        // Opposite edge is in opposite direction
                        if (ptjp1 == pti && ptj == ptip1)
                        {
                            adjtrioff(trioff  + i) = tri2off + j;
                            adjtrioff(tri2off + j) = trioff  + i;
                            found = true;
                            break;
                        }
                    }
                }
            }
        }

        // Loop through the edges until none need flipping
        bool edgesflipped;
        int64 count = 0;
        do
        {
            edgesflipped = false;
            for (int tri = 0; tri < ntris; ++tri)
            {
                int trioff = 3 * tri;
                for (int i = 0; i < 3; ++i)
                {
                    if (adjtrioff(trioff + i) < 0)
                        continue;

                    int j = adjtrioff(trioff + i) % 3;
                    int tri2off = adjtrioff(trioff + i) - j;
                    int tri2 = tri2off / 3;

                    int a = tripts(trioff + ((i + 2) % 3));
                    int b = tripts(trioff + i);
                    int c = tripts(trioff + ((i + 1) % 3));
                    int d = tripts(tri2off + ((j + 2) % 3));

                    if (!shouldFlipEdge(pts(a), pts(b), pts(c), pts(d)))
                        continue;

                    /* Flip the edge
                     *
                     *      j+2                j+2 i+1
                     *      / \                  /|\
                     *     /   \                / | \
                     *    / tri2\              / t|t \
                     *   j-------j+1 flips to j  r|r  i
                     * i+1\ tri /i             \ i|i /
                     *     \   /                \2| /
                     *      \ /                  \|/
                     *      i+2                j+1 i+2
                     */

                    // Update the caller's geometry
                    on_flip(tri, tri2, i, j);

                    // Vertex at i+1 must be changed to vertex at j+2
                    // Vertex at j+1 must be changed to vertex at i+2
                    tripts(trioff  + (i + 1) % 3) = tripts(tri2off + (j + 2) % 3);
                    tripts(tri2off + (j + 1) % 3) = tripts(trioff  + (i + 2) % 3);

                    // Edge i is where edge j+1 was
                    int nowacrossedgei = adjtrioff(tri2off + (j + 1) % 3);
                    adjtrioff(trioff  + i) = nowacrossedgei;
                    if (nowacrossedgei >= 0)
                        adjtrioff(nowacrossedgei) = trioff  + i;
                    // Edge j is where edge i+1 was
                    int nowacrossedgej = adjtrioff(trioff  + (i + 1) % 3);
                    adjtrioff(tri2off + j) = nowacrossedgej;
                    if (nowacrossedgej >= 0)
                        adjtrioff(nowacrossedgej) = tri2off + j;
                    // Edges i+1 and j+1 are between the two triangles
                    adjtrioff(trioff  + (i + 1) % 3) = tri2off + (j + 1) % 3;
                    adjtrioff(tri2off + (j + 1) % 3) = trioff  + (i + 1) % 3;

                    edgesflipped = true;

                    // Edge i+1 is now the one between these two triangles,
                    // so there's no need to check it next.
                    // All 3 edges will be checked again on the next pass anyway.
                    break;
                }
            }
            ++count;
        } while (edgesflipped && count < int64(ntris) * ntris);
    }

    /// As above, without a flip callback.
    template <typename PTS, typename TRIS>
    static void flipEdgesToAvoidSmallAngles(const PTS &pts, TRIS &tripts)
    {
        flipEdgesToAvoidSmallAngles(pts, tripts, [](int, int, int, int) {});
    }

    virtual int64 getMemoryUsage(bool inclusive) const;

private:
    UT_ConvexPoint      *polySetup(bool avoiddegeneracy);
    int                 leftOrRight(int v1, int v2, int v3, bool &double_back);
    int                 leftOrRight(int v1, int v2, int v3)
                        {
                            bool double_back;
                            return leftOrRight(v1, v2, v3, double_back);
                        }
    uint                leftMost(uint a, uint b, int ccwise);
    bool                testConcave(uint first, uint last);
    bool                arePointsSame(uint i, uint j);
    bool                isDegenerate(uint first, uint last);
    bool                addPolygon(
                                UT_Interrupt *boss, uint first,
                                uint last, bool lastone,
                                bool avoiddegeneracy);
    void                addTriangle(uint a, uint b, uint c,
                                    bool lastone);
    void                doSplit(uint first, uint newfirst,
                                uint newlast, uint last);
    bool                cornerTest(int side, uint la, uint lb,
                                   uint first, uint last);
    bool                insideOut(uint k, uint la,
                                  uint lb, uint first, uint last) const;
    uint                chooseCoincidentIntruder(
                                uint l, uint k1, uint k2,
                                uint first, uint last, uint ccwise);
    uint                findClosestInRange(uint lb, uint la,
                                 uint ls, uint first, uint last);
    uint                testIntersect(uint l, uint lb, uint la,
                                      uint first, uint last, int &ccwise);

protected:
    /// NOTE: You *must* set myData to NULL before
    ///       calling isConcave() or convex() !
    UT_ConvexPoint      *myData;

private:
    UT_ConvexPoint      *myConvexInfo;
    uint                 myConvexSize;
    uint                 myConvexEntries;
    int                 *myPointList;
    uint                 myPointSize;

    /// NOTE: The only purpose of this bool is because isConvex() must call
    ///       polySetup(true), so if there were degenerate points and
    ///       convex(__, __, false) is then called, it needs to know to call
    ///       polySetup(false) to get those omitted points.
    bool                 myAvoidDegeneracyPrev;
};

#endif