GEO_ImplicitSurface.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:        GEO_ImplicitSurface.h
 *
 * COMMENTS:    Compound implicit surface evaluation with recursive operations.
 */

#pragma once


#ifndef __GEO_ImplicitSurface_h__
#define __GEO_ImplicitSurface_h__

#include "GEO_ImplicitSurfaceOp.h"
#include "GEO_ImplicitBasicSurface.h"

#include <UT/UT_Vector.h>
#include <UT/UT_Tracing.h>
#include <UT/UT_BoundingBox.h>
#include <GA/GA_Handle.h>
#include <GA/GA_PrimVolumeXform.h>
#include <GEO/GEO_Detail.h>


//////////////////////////////////////////////////////////////////////////////
// GEO_ImplicitSurface
//////////////////////////////////////////////////////////////////////////////

/// Compound Implicit Surface - evaluates compound SDFs built from basic shapes
/// combined with boolean/blending operations.
///
/// Constructed from a GEO_Detail whose points carry implicit surface attributes.
/// The evaluation uses a recursive algorithm to handle nested Begin/End groups.
///
/// Shape/operation sequence structure:
///   - Begin op: starts a new group (recursive call)
///   - End shape + Binary op: ends group, returns to caller who applies the op
///   - Unary op: modifies current SDF in place
///   - Binary op + regular shape: evaluates shape, combines with current SDF
template <VEX_Precision PREC>
class GEO_ImplicitSurface
{
public:
    using Float = VEXfloat<PREC>;
    using Vec3 = VEXvec3<PREC>;
    using Vec4 = UT_Vector4T<Float>;
    using Quat = UT_QuaternionT<Float>;
    using Shape = GEO_ImplicitBasicSurface<PREC>;
    using Op = GEO_ImplicitSurfaceOp<PREC>;

    /// Time derivatives of shape/operation parameters for computing dparm
    struct DParms
    {
        /// Linear velocity (time derivative of position)
        Vec3 v = {0, 0, 0};
        /// Angular velocity in world space: vel(x) = cross(w, x - pos) + v
        Vec3 w = {0, 0, 0};
        /// Time derivative of shape parameters (e.g., change in box size)
        Vec4 shape_dparms = {0, 0, 0, 0};
        /// Time derivative of operation parameters (e.g., change in smooth k)
        Vec4 op_dparms = {0, 0, 0, 0};
    };

    using IntervalF = UT_IntervalT<Float>;
    using BBox3 = UT_BoundingBoxT<Float>;

    // Data members
    UT_ValArray<Shape> myShapes;
    UT_ValArray<Op> myOps;
    UT_ValArray<DParms> myDParms;
    int myNPoints = 0;

    GEO_ImplicitSurface() = default;

    /// Construct from geometry, reading all implicit surface attributes
    /// from point data into the internal arrays.
    GEO_ImplicitSurface(const GEO_Detail &geo)
    {
        utZoneScopedN("implicit surface - init");
        
        myNPoints = (int)geo.getNumPoints();
        if (myNPoints == 0)
            return;

        const GA_IndexMap &ptmap = geo.getPointMap();

        // Bind attribute handles
        GA_ROHandleV3 h_P(&geo, GA_ATTRIB_POINT, "P");
        GA_ROHandleM3 h_transform(&geo, GA_ATTRIB_POINT, "transform");
        
        GA_ROHandleS h_shape(&geo, GA_ATTRIB_POINT, "shape");
        GA_ROHandleV4 h_shapevals(&geo, GA_ATTRIB_POINT, "shapevals");

        GA_ROHandleS h_op(&geo, GA_ATTRIB_POINT, "implicitsurface");
        GA_ROHandleV4 h_opvals(&geo, GA_ATTRIB_POINT, "implicitsurfacevals");

        // Optional velocity attributes
        GA_ROHandleV3 h_v(&geo, GA_ATTRIB_POINT, "v");
        GA_ROHandleV3 h_w(&geo, GA_ATTRIB_POINT, "w");
        GA_ROHandleV4 h_shapedvals(&geo, GA_ATTRIB_POINT, "shapedvals");
        GA_ROHandleV4 h_opdvals(&geo, GA_ATTRIB_POINT, "implicitsurfacedvals");

        if (!h_shape.isValid() || !h_op.isValid())
        {
            myNPoints = 0;
            return;
        }

        myShapes.setSizeNoInit(myNPoints);
        myOps.setSizeNoInit(myNPoints);
        myDParms.setSizeNoInit(myNPoints);

        for (int i = 0; i < myNPoints; i++)
        {
            GA_Offset off = ptmap.offsetFromIndex(i);

            // Shape
            VEXmat3<PREC> xform;
            if (h_transform.isValid())
                xform = h_transform.get(off);
            else
                xform.identity();
            myShapes[i] = Shape(h_shape.get(off),
                                h_P.isValid() ? Vec3(h_P.get(off)) : Vec3(0, 0, 0),
                                xform,
                                h_shapevals.isValid() ? Vec4(h_shapevals.get(off)) : Vec4(1, 1, 1, 0));

            // Op
            myOps[i] = Op(h_op.get(off),
                          xform,
                          h_opvals.isValid() ? Vec4(h_opvals.get(off)) : Vec4(0, 0, 0, 0));

            // DParms
            DParms &dp = myDParms[i];
            dp.v = h_v.isValid() ? Vec3(h_v.get(off)) : Vec3(0, 0, 0);
            dp.w = h_w.isValid() ? Vec3(h_w.get(off)) : Vec3(0, 0, 0);
            dp.shape_dparms = h_shapedvals.isValid() ? Vec4(h_shapedvals.get(off)) : Vec4(0, 0, 0, 0);
            dp.op_dparms = h_opdvals.isValid() ? Vec4(h_opdvals.get(off)) : Vec4(0, 0, 0, 0);
        }
    }

    // ========================================
    // Transform
    // ========================================

    /// Apply translation, rotation, and uniform scale to all shapes.
    /// Positions: new_pos = rotate(old_pos * uniform_scale) + translate
    /// Orientations: composed with the rotation quaternion
    /// Distance parameters: scaled by uniform_scale on shapes and ops
    /// NOTE: DParms (v, w) are not transformed by this method.
    void applyTransform(const Vec3 &translate, const Quat &rotate,
                        Float uniform_scale)
    {
        VEXmat3<PREC> rot_mat;
        rotate.getRotationMatrix(rot_mat);
        for (int i = 0; i < myNPoints; i++)
        {
            myShapes[i].myPos = rotate.rotate(myShapes[i].myPos * uniform_scale)
                              + translate;
            myShapes[i].myTransform = rot_mat * myShapes[i].myTransform;

            myShapes[i].scale(uniform_scale);
            myOps[i].scale(uniform_scale);
        }
    }

    /// Transform the implicit surface from world space into the index space
    /// of a VDB defined by the given GEO_PrimVolumeXform.
    /// After this call, evaluating positions in index space will produce
    /// SDF values in index-space units (multiply by voxelsize for world).
    void applyIndexSpaceTransform(const GA_PrimVolumeXform &index_xform)
    {
        // Decompose world->index into translate, rotate, uniform_scale.
        // toVoxelSpace(pos) - 0.5 = (pos - center) * myInverseXform * 0.5
        UT_Matrix3D rot_mat(index_xform.myInverseXform);
        UT_Vector3D scales;
        rot_mat.extractScalesT(scales, (UT_Vector3D *)nullptr);
        Float uniform_scale = Float(scales.x()) * Float(0.5);
        Quat rotate;
        rotate.updateFromRotationMatrix(rot_mat);
        Vec3 center(index_xform.myCenter);
        Vec3 translate = -rotate.rotate(center * uniform_scale);

        applyTransform(translate, rotate, uniform_scale);
    }

    /// Apply offset to the surface
    void applyOffset(const Float offset);

    // ========================================
    // Evaluation methods
    // ========================================

    Float operator()(const Vec3 &P) const;
    Float operator()(const Vec3 &P, Vec3 &grad) const;
    Float operator()(const Vec3 &P, Vec3 &grad, Float &dparm) const;

    // ========================================
    // Batched evaluation (surfaces outer, points inner)
    // ========================================

    void evalBatched(int npts, const Vec3 *P, Float *sdf) const;
    void evalBatched(int npts, const Vec3 *P, Float *sdf, Vec3 *grad) const;
    void evalBatched(int npts, const Vec3 *P, Float *sdf, Vec3 *grad,
                     Float *dparm) const;

    /// Same as above, but also outputs the winning shape id per point.
    /// `grad`/`dparm` may be nullptr if not needed.
    void evalBatched(int npts, const Vec3 *P, Float *sdf, int *shape_id,
                     Vec3 *grad, Float *dparm) const;

    /// Batched evaluation that also interpolates an arbitrary point
    /// attribute (float or vector) across the blended shapes.
    void evalBatchedAttr(int npts, const Vec3 *P, Float *sdf, int *shape_id,
                         Vec3 *grad, Float *dparm,
                         const GA_ROHandleF &h_attr, const GA_IndexMap &ptmap,
                         Float *attr_out) const;
    void evalBatchedAttr(int npts, const Vec3 *P, Float *sdf, int *shape_id,
                         Vec3 *grad, Float *dparm,
                         const GA_ROHandleV3 &h_attr, const GA_IndexMap &ptmap,
                         Vec3 *attr_out) const;

    /// Batched evaluation that also blends the per-shape rigid velocity
    /// (v + w x (P - pos)) across the blended shapes.
    void evalBatchedVel(int npts, const Vec3 *P, Float *sdf, int *shape_id,
                        Vec3 *grad, Float *dparm, Vec3 *rigid_vel) const;

    // ========================================
    // Interval evaluation (for spatial pruning)
    // ========================================

    /// Evaluate SDF over an interval region, returning relevant shape IDs.
    IntervalF evalInterval(const BBox3 &P, UT_Array<int> &relevant_ids) const;

    /// Evaluate SDF over an interval region
    IntervalF evalInterval(const BBox3 &P) const;

    /// Evaluate SDF over an interval region
    IntervalF evalInterval(const Vec3 &center, const Float radius) const;

    /// Evaluate interval over a sphere and build a pruned surface expression.
    /// Pruning follows operation interval optimization flags and emits a new
    /// shape/op sequence equivalent on the queried sphere.
    IntervalF outputBoundAndPrune(
            const Vec3 &center,
            const Float radius,
            GEO_ImplicitSurface<PREC> &pruned) const;

    /// Return a conservative estimate on the isoline offset.
    /// This is used whe we get bounding box based on the visualization geometry,
    /// we expand this bounding box by this offset to get a better estimate
    /// on a region where implicut surface is contained.
    Float offsetEstimate(){
        Float offset = 0.0;
        for (auto& op : myOps)
        {
            if (op.myId == GEO_ImplicitSurfaceOpID::Offset)
            {
                offset += SYSabs(op.myParms[0]);
            }
            if (op.myId == GEO_ImplicitSurfaceOpID::Shell)
            {
                offset += SYSabs(op.myParms[0]) + SYSabs(op.myParms[1])/2;
            }
        }
        return offset;
    }

private:
    // ========================================
    // Internal recursive evaluation helpers
    // ========================================
    // Not part of the public interface: these are implemented and
    // implicitly instantiated only inside GEO_ImplicitSurface.C, so their
    // definitions (and the lambda closure types some of them are called
    // with) never need to be visible outside that translation unit.

    /// Recursive single-point evaluation of a shape/operation group.
    /// Returns the index of the End that terminated the group, or
    /// myNPoints if no End was found.
    template <bool WantGrad, bool WantDParm>
    int evalGroup(int start, const Vec3 &P, Float &sdf, int &shape_id,
                 Vec3 *grad, Float *dparm) const;

    /// Recursive batched evaluation of a shape/operation group. `valFn`
    /// computes the interpolated `Val` contributed by shape `i` at point
    /// P when WantInterp is true; ignored otherwise.
    template <bool WantGrad, bool WantDParm, bool WantInterp,
             typename Val, typename ValFn>
    int evalGroupBatched(int start, int npts, const Vec3 *P, Float *sdf,
                         int *shape_id, Vec3 *grad, Float *dparm,
                         const ValFn &valFn, Val *interp_val) const;

    /// Recursive interval evaluation with relevant_ids pruning.
    int evalGroupInterval(int start, const BBox3 &P, IntervalF &sdf,
                          UT_Array<int> &relevant_ids) const;

    /// Recursive interval evaluation over a bounding sphere.
    int evalGroupIntervalSphere(int start, const Vec3 &center, Float radius,
                                IntervalF &sdf) const;
};

#endif // __GEO_ImplicitSurface_h__