GU_Watershed.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:        SOP library (C++)
 *
 * COMMENTS:
 *
 * The GU_Watershed class implements a Morse-Smale complex on a polygonal
 * mesh from a geometry detail, by constructing a discrete gradient field
 * using a scalar function defined on points, primitives, or both. The
 * definition of the complex and discrete gradient field follow the work
 * of Robin Forman in discrete Morse theory.
 *
 * In addition, we implement a persistence-based simplification of the complex
 * (Cancellation of critical pairs based on the persistent homology of sub- or
 * super-level sets). This is done using the algorithm of:
 *
 *   Optimal topological simplification of discrete functions on surfaces
 *   Ulrich Bauer, Carsten Lange, and Max Wardetzky. DCG 47:2 (2012), 347–377
 *
 * Note that this simplification does not alter function values
 * and merely merges complex cells.
 */

#ifndef __GU_Watershed_h__
#define __GU_Watershed_h__

#include "GU_API.h"

#include "GU_Detail.h"
#include <GEO/GEO_HedgeInterface.h>
#include <UT/UT_Classifier.h>
#include <UT/UT_Function.h>
#include <UT/UT_Map.h>
#include <UT/UT_Set.h>
#include <UT/UT_UniquePtr.h>


class GU_API GU_Watershed
{
public:
    using HedgeInterface = GEO_DetachedHedgeInterface;
    using HedgeArray = UT_Array<GEO_Hedge>;
    using GetFunctor = UT_Function<fpreal(GA_Offset)>;

                         GU_Watershed(const GU_Detail *gdp,
                                const GA_PrimitiveGroup *prims,
                                const GetFunctor &point_fn,
                                GetFunctor face_fn = nullptr,
                                const HedgeInterface *hip = nullptr);

    void                 simplify(fpreal tol);

    const
    GA_OffsetArray      &minima() const { return myMins; }
    int                  numMinima() const { return int(myMins.size()); }

    const
    GA_OffsetArray      &maxima() const { return myMaxs; }
    int                  numMaxima() const { return int(myMaxs.size()); }

    const HedgeArray    &saddles() const { return mySaddles; }

    int                  numSaddles() const { return int(mySaddles.size()); }
    int                  saddleSign(int i) const { return mySaddleSign(i); }
    GEO_Hedge            saddleHedge(int i) const { return mySaddles(i); }

    bool                 isSaddleCanceled(int i) const
                            { return mySaddleSign(i) == 0; }

    void                 getPosSaddles(HedgeArray &saddles) const;
    void                 getNegSaddles(HedgeArray &saddles) const;

    void                 orientMinTree();
    void                 orientMaxTree();

    int                  faceMax(GA_Offset poly) const
                            { return maxTreeRoot(poly); }

    int                  pointMin(GA_Offset pt) const
                            { return minTreeRoot(pt); }

    fpreal               saddleFn(GEO_Hedge h) const
                            { return SYSmax(myPointFn(srcPoint(h)),
                                            myPointFn(dstPoint(h))); }

    void                 dump();

    int                  faceMin(GA_Offset poly) const;
    int                  pointMax(GA_Offset pt) const;

    SYS_FORCE_INLINE
    GEO_Hedge            pointMinHedge(GA_Offset pt) const
    {
        return GEO_Hedge(GA_Offset(myParent.get(pt)));
    }

    SYS_FORCE_INLINE
    GEO_Hedge            faceMaxHedge(GA_Offset poly) const
    {
        return GEO_Hedge(GA_Offset(myCoparent.get(poly)));
    }


    GA_Offset            pointDescentSucc(GA_Offset pt) const
    {
        auto h = pointMinHedge(pt);
        if (h == GEO_INVALID_HEDGE)
            return GA_INVALID_OFFSET;
        return srcPoint(h) == pt ? dstPoint(h) : srcPoint(h);
    }

    GA_Offset            faceAscentSucc(GA_Offset poly) const
    {
        auto h = faceMaxHedge(poly);
        if (h == GEO_INVALID_HEDGE)
            return GA_INVALID_OFFSET;
        return hedgePoly(sym(h));
    }

    bool                 isInMinTree(GEO_Hedge h) const;
    bool                 isInMaxTree(GEO_Hedge h) const;

    SYS_FORCE_INLINE
    bool                 isSaddle(GEO_Hedge h) const
                            { return !isInMinTree(h) && !isInMaxTree(h); }

    int                  negSaddlePairedMin(GEO_Hedge h) const;
    int                  posSaddlePairedMax(GEO_Hedge h) const;


    fpreal               posSaddlePersistence(GEO_Hedge h) const;
    fpreal               negSaddlePersistence(GEO_Hedge h) const;


    bool                 isRidgeEdge(GEO_Hedge h) const;
    bool                 isValleyEdge(GEO_Hedge h) const;

    void                 pointDescentPath(GA_Offset pt,
                                HedgeArray &path) const;

    void                 faceAscentDualPath(GA_Offset prim,
                                HedgeArray &path) const;

    bool                 isMinCanceled(int i) const
                            { return pointMin(myMins(i)) != i; }

    bool                 isMaxCanceled(int i) const
                            { return faceMax(myMaxs(i)) != i; }

private:
    SYS_FORCE_INLINE
    GA_OffsetListRef     vertexList(GA_Offset poly) const
                            { return myGdp->getPrimitiveVertexList(poly); }

    SYS_FORCE_INLINE
    GA_Offset            vertexPoint(GA_Offset vtx) const
                            { return myGdp->vertexPoint(vtx); }

    SYS_FORCE_INLINE
    GA_Offset            vertexPoly(GA_Offset vtx) const
                            { return myGdp->vertexPrimitive(vtx); }

    SYS_FORCE_INLINE
    GA_Offset            hedgePoly(GEO_Hedge h) const
                            { return vertexPoly(myHip->srcVertex(h)); }

    SYS_FORCE_INLINE
    GEO_Hedge            polyHedge(GA_Offset poly) const
    {
        return GEO_Hedge(myGdp->getPrimitiveVertexList(poly)(0));
    }

    SYS_FORCE_INLINE
    GA_Offset            srcVertex(GEO_Hedge h) const
                            { return myHip->srcVertex(h); }

    SYS_FORCE_INLINE
    GA_Offset            dstVertex(GEO_Hedge h) const
                            { return myHip->dstVertex(h); }

    SYS_FORCE_INLINE
    GA_Offset            srcPoint(GEO_Hedge h) const
                            { return myHip->srcPoint(h); }

    SYS_FORCE_INLINE
    GA_Offset            dstPoint(GEO_Hedge h) const
                            { return myHip->dstPoint(h); }

    SYS_FORCE_INLINE
    GEO_Hedge            sym(GEO_Hedge h) const
                            { return myHip->sym(h); }


    SYS_FORCE_INLINE
    GEO_Hedge            lnext(GEO_Hedge h) const
                            { return myHip->lnext(h); }

    SYS_FORCE_INLINE
    GEO_Hedge            lprev(GEO_Hedge h) const
                            { return myHip->lprev(h); }

    using ClassRootMap = UT_Map<GA_Offset, int>;

    SYS_FORCE_INLINE
    void                 setClassRoot(ClassRootMap &root_map,
                                const UT_Classifier &classes,
                                GA_Offset elem, int root_id)
    {
        root_map[GA_Offset(classes.classRoot(elem))] = root_id;
    };

    int                  classRoot(const ClassRootMap &root_map,
                                const UT_Classifier &classes,
                                GA_Offset elem) const
    {
        auto it = root_map.find(GA_Offset(classes.classRoot(elem)));
        if (it == root_map.end())
            return -1;
        return it->second;
    }

    SYS_FORCE_INLINE
    void                 setMaxTreeRoot(GA_Offset poly, int root_id)
    {
        return setClassRoot(myMaxTreeRoot, myMaxTrees, poly, root_id);
    };

    SYS_FORCE_INLINE
    int                  maxTreeRoot(GA_Offset poly) const
    {
        return classRoot(myMaxTreeRoot, myMaxTrees, poly);
    };

    SYS_FORCE_INLINE
    void                 setMinTreeRoot(GA_Offset pt, int root_id)
    {
        return setClassRoot(myMinTreeRoot, myMinTrees, pt, root_id);
    };

    SYS_FORCE_INLINE
    int                  minTreeRoot(GA_Offset rep_pt) const
    {
        return classRoot(myMinTreeRoot, myMinTrees, rep_pt);
    };

    bool                 mergeMinTrees(GEO_Hedge h);
    bool                 mergeMaxTrees(GEO_Hedge h);

    using HedgeInterfaceUptr = UT_UniquePtr<HedgeInterface>;


    ClassRootMap         myMaxTreeRoot, myMinTreeRoot;
    UT_Classifier        myMaxTrees, myMinTrees;

    GA_VertexGroupUPtr   myMaxForest;
    GA_VertexGroupUPtr   myMinForest;

    HedgeInterfaceUptr   myOwnHip;
    const GU_Detail     *myGdp;

    const
    HedgeInterface      *myHip;

    const GetFunctor     myPointFn;
    GetFunctor           myFaceFn;

    GA_AttributeUPtr     myParentOwner;
    GA_AttributeUPtr     myCoparentOwner;

    GA_RWHandleI         myParent;
    GA_RWHandleI         myCoparent;

    GA_OffsetArray       myMins;
    GA_OffsetArray       myMaxs;

    HedgeArray           mySaddles;
    UT_IntArray          mySaddleSign;
};

#endif