Virtual Support Function Dispatch for Coal

[rjoomen]

This is a design brief, created for discussion as an alternative option to PR #822.

The current dispatch is manual virtual dispatch

Coal's GJK support function dispatch uses NODE_TYPE switches, template specializations, and function pointers to route each shape to its support function body. This is virtual dispatch — implemented with templates and switches instead of a vtable.

The call chain from GJK to the actual support body is 4 levels deep:

MinkowskiDiff::support()          → getSupportFunc(...)           [function pointer]
  → getSupportFuncTpl<S0,S1,I,O>  → getSupportTpl<S0,S1,I,O>    [template trampoline]
    → getShapeSupport<O>(S0*, ...) + getShapeSupport<O>(S1*, ...) [actual body, cross-TU]

The function pointer is resolved once per collision pair in MinkowskiDiff::set() via two nested switches (makeGetSupportFunction0 → makeGetSupportFunction1), producing ~400 template instantiations (10 shapes x 10 x 2 identity x 2 SupportOptions).

This manual dispatch has three structural problems:

1. Per-shape logic is separated from the shapes. Each shape's support function lives in support_functions.cpp as a template specialization of getShapeSupport<T>, not on the shape class itself. To understand how Box computes its support point, you read support_functions.cpp, not Box.

2. Per-shape configuration is embedded in dispatch machinery. makeGetSupportFunction0/1 doesn't just dispatch — it also handles per-shape setup that only those shapes know about:

• Sphere/Capsule radius promotion (swept_sphere_radius += radius) — duplicated in both functions
• ConvexBase visited-array allocation and Small/Large algorithm selection — duplicated in both functions
• Nesterov heuristic — a 35-line switch that reads a compile-time trait to produce a runtime bool, for a flag that is true for exactly one shape type

3. The shape list is duplicated in 6 switch sites across 4 files. makeGetSupportFunction0, makeGetSupportFunction1, getSupport, getSupportSet, getNormalizeSupportDirection, makeSupportSetFunction — totaling ~350 lines of boilerplate. Adding a shape means editing all six and matching the existing patterns for radius promotion, data initialization, and traits.

Proposed design: use actual virtual dispatch

Replace the manual dispatch with virtual methods on ShapeBase. Each shape owns its support logic, its inflation semantics, and its setup requirements. The dispatch machinery disappears.

Virtual methods on ShapeBase

/// @brief Compute the support point of this shape in the given direction.
///
/// Returns the core support point in the shape's local frame, WITHOUT
/// swept-sphere inflation.
///
/// @param[in] dir Support direction.  May be non-unit-length — the existing
///   template specializations already receive raw directions from GJK.
///   Shapes whose formula requires |dir| = 1 (Ellipsoid, Cone) normalize
///   internally, same as they do today.
/// @param[out] support The computed support point.
/// @param[in,out] hint Warm-start index (ConvexBase).
/// @param[in,out] data Scratch space (ConvexBase uses data.visited,
///   data.last_dir; other shapes ignore it).
virtual void getShapeSupport(const Vec3s& dir, Vec3s& support,
                             int& hint, ShapeSupportData& data) const = 0;

/// @brief Total inflation radius for GJK swept-sphere correction.
///
/// Under NoSweptSphere mode, GJK operates on the bare core shape and
/// stores this radius for post-convergence distance/witness correction.
/// Default: getSweptSphereRadius().
/// Sphere overrides: radius + getSweptSphereRadius() (point core).
/// Capsule overrides: radius + getSweptSphereRadius() (segment core).
virtual Scalar getInflationRadius() const { return getSweptSphereRadius(); }

/// @brief Initialize per-shape scratch data before a GJK/EPA solve.
///
/// Called once by MinkowskiDiff::set(), not per iteration.
/// Default: no-op.  ConvexBaseTpl overrides to pre-size data.visited
/// for large convex shapes.
virtual void initSupportData(ShapeSupportData& data) const {
  COAL_UNUSED_VARIABLE(data);
}

/// @brief Whether this shape needs the Nesterov normalize heuristic.
///
/// True for ConvexBase (polytopes where support is non-smooth).
/// MinkowskiDiff::set() ANDs both shapes' values.
virtual bool needNesterovNormalizeHeuristic() const { return false; }

/// @brief Whether this shape is strictly convex (single support point
/// in every direction).
///
/// True for Sphere, Ellipsoid.  Used by contact patch computation to
/// skip support set queries.  Replaces shape_traits<T>::IsStrictlyConvex.
virtual bool isStrictlyConvex() const { return false; }

/// @brief Compute the support set (contact patch vertices).
///
/// Returns all points within @p tol of the extremal support value in
/// direction support_set.getNormal().  Default: single-point fallback
/// via getShapeSupport().  Shapes with flat faces override for richer
/// contact patches.
virtual void getShapeSupportSet(SupportSet& support_set, int& hint,
                                ShapeSupportData& data,
                                size_t num_sampled_supports,
                                Scalar tol) const;

getShapeSupport is pure virtual — every concrete ShapeBase must provide a support function. The compiler catches missing overrides. Shapes that should never reach GJK (Plane, Halfspace) implement an override that throws — these shapes have analytical specializations for every pair, so a GJK call indicates a bug.

What moves where

Each built-in shape moves its getShapeSupport<NoSweptSphere> template specialization body into a getShapeSupport override marked final. Per-shape configuration moves from the dispatch machinery onto the shapes:

Shape	Support body (from support_functions.cpp)	Configuration (from makeGetSupportFunction0/1)
Box	Per-axis sign selection with halfSide	—
Sphere	support.setZero() (point core)	getInflationRadius → radius + getSweptSphereRadius()
Ellipsoid	v / sqrt(v . dir) (normalizes internally)	—
Capsule	±halfLength on Z (segment core)	getInflationRadius → radius + getSweptSphereRadius()
Cone	Tip-or-rim conditional (normalizes internally)	—
Cylinder	±halfLength on Z + radial	—
TriangleP	Max dot of 3 vertices	—
ConvexBaseTpl	Linear scan or hill-climbing (decides internally)	initSupportData (visited array), needNesterovNormalizeHeuristic
Plane	Throws (should never reach GJK)	—
Halfspace	Throws (should never reach GJK)	—

ConvexBaseTpl deserves a note: currently the dispatch machinery selects between SmallConvex and LargeConvex tag types based on vertex count. With virtual dispatch, ConvexBaseTpl::getShapeSupport makes this decision internally. The tag types become a private implementation detail.

Marking overrides final allows the compiler to devirtualize when the concrete type is statically known.

Simplified MinkowskiDiff

The ~250 lines of dispatch machinery in minkowski_difference.cpp reduce to:

struct MinkowskiDiff {
  const ShapeBase* shapes[2];
  Matrix3s oR1;           // relative rotation
  Vec3s ot1;              // relative translation
  bool identity;          // true when oR1 ≈ I, ot1 ≈ 0
  bool normalize_support_direction;  // Nesterov heuristic
  Eigen::Array<Scalar, 1, 2> swept_sphere_radius;

  // Mutable: support data is internal cache state (visited-vertex sets,
  // warm-start hints) mutated during const support() queries.
  // Replaces the current const_cast<ShapeSupportData*>(data).
  mutable ShapeSupportData data[2];

  template <int _SupportOptions>
  void set(const ShapeBase* s0, const ShapeBase* s1,
           const Transform3s& tf0, const Transform3s& tf1) {
    shapes[0] = s0; shapes[1] = s1;

    // Relative transform
    oR1.noalias() = tf0.getRotation().transpose() * tf1.getRotation();
    ot1.noalias() = tf0.getRotation().transpose() *
                    (tf1.getTranslation() - tf0.getTranslation());
    identity = (oR1.isIdentity() && ot1.isZero());

    // Swept sphere radius.
    // WithSweptSphere: inflation is baked into the support return value,
    // so the external field is zeroed.
    // NoSweptSphere: GJK operates on the bare core shape; the radius is
    // stored here for post-convergence correction.
    if (_SupportOptions == SupportOptions::WithSweptSphere) {
      swept_sphere_radius.setZero();
    } else {
      swept_sphere_radius[0] = s0->getInflationRadius();
      swept_sphere_radius[1] = s1->getInflationRadius();
    }

    // Per-shape scratch data
    s0->initSupportData(data[0]);
    s1->initSupportData(data[1]);

    // Nesterov heuristic
    normalize_support_direction =
        s0->needNesterovNormalizeHeuristic() &&
        s1->needNesterovNormalizeHeuristic();
  }

  inline void support(const Vec3s& dir, Vec3s& supp0, Vec3s& supp1,
                      support_func_guess_t& hint) const {
    shapes[0]->getShapeSupport(dir, supp0, hint[0], data[0]);

    if (identity) {
      shapes[1]->getShapeSupport(-dir, supp1, hint[1], data[1]);
    } else {
      const Vec3s dir1 = -(oR1.transpose() * dir);
      shapes[1]->getShapeSupport(dir1, supp1, hint[1], data[1]);
      supp1.noalias() = oR1 * supp1 + ot1;
    }
  }
};

Similarly, contact patch dispatch simplifies from a switch + function pointer (makeSupportSetFunction) to a direct virtual call (shape->getShapeSupportSet(...)).

Removed from minkowski_difference.cpp: GetSupportFunction typedef, getSupportFunc pointer, getSupportTpl, getSupportFuncTpl, makeGetSupportFunction0, makeGetSupportFunction1, getNormalizeSupportDirection, getNormalizeSupportDirectionFromShapes.

Consequence for function matrices

MinkowskiDiff::set() no longer needs concrete types — it takes ShapeBase* and calls virtuals. This affects all three function matrices:

Collision and distance matrices. (CollisionFunctionMatrix/DistanceFunctionMatrix) Each has 121 explicit shape-shape entries (11 x 11) using ShapeShapeCollide<S0, S1> / ShapeShapeDistance<S0, S1> — the template types exist solely to thread concrete types through to MinkowskiDiff::set(). Of the 55 unique pairs, 26 have analytical specializations (Sphere-Sphere, Sphere-Capsule, *-Plane, *-Halfspace, etc.) that bypass GJK and need concrete types. The remaining 29 go through GJK and can use a generic <ShapeBase, ShapeBase> fallback. Each matrix shrinks from 121 explicit entries to ~26.

Contact patch matrix. (ContactPatchFunctionMatrix) The contact patch solver calls makeSupportSetFunction, which is one of the 6 switch sites eliminated by this refactoring. With virtual getShapeSupportSet, the solver doesn't need concrete types at all — a single generic entry covers all shape-shape pairs.

In all three matrices, adding a shape no longer requires adding matrix entries — GJK collision, distance, and contact patches work automatically via the virtual interface. This is a follow-up change enabled by the refactoring.

What the current design does well

The current approach has genuine performance advantages that the refactoring trades away:

• One indirect call per iteration, not two. The pair function pointer (getSupportFuncTpl<Box, Sphere, false, NoSS>) dispatches both shapes in a single indirect call. Inside, the two getShapeSupport calls are direct (concrete type baked in via template parameter). Virtual dispatch replaces this with two indirect calls. The total call count drops (2 vs 3), but the indirect count rises (2 vs 1) — indirect calls have slightly higher latency due to branch prediction pipeline effects.

• Identity transform compiled out. TransformIsIdentity is a template parameter that selects between two function pointer instantiations at set() time. The identity instantiation contains no rotation/translation instructions at all — no oR1.transpose() * dir, no oR1 * supp1 + ot1. Virtual dispatch replaces this with a runtime if (identity) branch inside support(). The branch is well-predicted, but the non-identity code (matrix multiply + transform) is still compiled into the function body and occupies instruction cache.

• Inlinable under LTO. The concrete type is statically known throughout the call chain (getSupportFuncTpl<Box, ...> → getShapeSupport(const Box*, ...)). Under LTO, the linker can inline across TUs because the types are explicit. With virtual dispatch through ShapeBase*, LTO can only inline if it proves the dynamic type — harder even with final, since MinkowskiDiff stores ShapeBase* that could point to any concrete shape.

These are real tradeoffs. They're small — GJK iterations are dominated by the support body math (dot products, comparisons, hill-climbing), not call overhead — but should be benchmarked rather than assumed negligible.

Performance summary

Metric	Current	Proposed
Indirect calls per iteration	1	2
Direct non-inlinable calls per iteration	2	0
Identity transform	Compiled out (no rotation/translation instructions)	Runtime branch (predicted, but non-identity code in icache)
Inlinable with LTO	Yes (concrete types known)	Only with final + type proof
Template instantiations	~400	0
Binary size	~400 monomorphized bodies	~11 virtual method bodies

Risk areas

• Sphere/Capsule radius promotion. getInflationRadius() must reproduce the current accounting exactly. Wrong values silently produce incorrect collision results.
• ConvexBaseTpl data initialization. A missing initSupportData override causes silent hill-climbing failure. Add a debug assert in ConvexBaseTpl::getShapeSupport that data.visited is correctly sized for large meshes.
• Ellipsoid/Cone normalization. These shapes must normalize dir internally. The formulas are moved unchanged from the existing specializations, but should be verified.
• Identity transform tolerances. Must match the current isIdentity() / isZero() behavior.

Files to modify

File	Changes
include/coal/shape/geometric_shapes.h	Add virtuals to ShapeBase. Add final overrides to each built-in shape. Add #include "coal/narrowphase/support_data.h".
src/narrowphase/support_functions.cpp	Move each getShapeSupport<T> body into T::getShapeSupport. Simplify getSupport free function from 40-line switch to convenience wrapper (creates local ShapeSupportData, calls virtual). Remove per-shape template specializations.
include/coal/narrowphase/support_functions.h	Remove per-shape getShapeSupport declarations. Keep getSupport free function.
include/coal/narrowphase/minkowski_difference.h	Remove GetSupportFunction typedef and getSupportFunc pointer. Make data[2] mutable. Add identity member. Replace support() with virtual calls.
src/narrowphase/minkowski_difference.cpp	Remove all dispatch machinery. Simplify set() to ~20 lines.
src/contact_patch/contact_patch_solver.cpp	Remove makeSupportSetFunction and related templates. Call getShapeSupportSet directly.
include/coal/shape/geometric_shapes_traits.h	Removed entirely. NeedNesterovNormalizeHeuristic and IsStrictlyConvex replaced by virtuals. Three other traits (NeedNormalizedDir, IsInflatable, HasInflatedSupportFunction) are dead code.

All existing GJK/EPA/contact-patch tests must pass unchanged.

Follow-up: virtual computeAABB

computeBV<AABB, ShapeBase> has the same structural problem as the support dispatch: it's a template specialization that custom shapes can't override. After the fork's change to use aabb_local instead of support queries, it also has a precondition (computeLocalAABB() must have been called first) that makes it unsuitable for bootstrapping a shape's own AABB.

With virtual dispatch on ShapeBase, this becomes a virtual method with a default implementation:

/// @brief Compute the AABB of this shape at the given transform.
///
/// Default: rotated-AABB formula on aabb_local (|R|*half-extents).
/// Requires computeLocalAABB() to have been called.
/// Shapes that need exact AABBs or want to bootstrap from geometry
/// can override.
virtual void computeAABB(const Transform3s& tf, AABB& bv) const {
  const Matrix3s& R = tf.getRotation();
  const Vec3s& T = tf.getTranslation();
  const Vec3s half = (aabb_local.max_ - aabb_local.min_) * 0.5;
  const Vec3s center = R * aabb_local.center() + T;
  const Vec3s delta(R.cwiseAbs() * half);
  bv.min_ = center - delta;
  bv.max_ = center + delta;
}

Built-in shapes like Cylinder can override with their closed-form (R.cwiseAbs() * (radius, radius, halfLength)). Custom shapes override to compute from their own geometry. The precondition on the default is still there, but shapes that need to bootstrap (like DeformedCylinder) simply override instead of being forced through computeShapeSupport.

This also unifies computeBV<AABB, T> (currently one specialization per shape) with CollisionObject::computeAABB() (which already uses the aabb_local formula). Both would delegate to the same virtual.

Not in scope for the initial refactoring — the support dispatch is the priority. But it's a natural follow-up that completes the move from template specializations to virtual methods.

All existing GJK/EPA/contact-patch tests must pass unchanged.

API impact

Top-level API unchanged. coal::collide(), coal::distance(), coal::computeContactPatch(), and all request/result types are unaffected.

ShapeBase gains new virtuals. Most notably a pure virtual getShapeSupport. This does not change ShapeBase's instantiability — it is already abstract (inherits 4 pure virtuals from CollisionGeometry: clone, computeLocalAABB, getNodeType, isEqual). External ShapeBase subclasses must add a getShapeSupport override; the remaining new virtuals (getInflationRadius, initSupportData, needNesterovNormalizeHeuristic, isStrictlyConvex, getShapeSupportSet) have defaults.

MinkowskiDiff internals change, public methods don't. GetSupportFunction typedef and getSupportFunc pointer are removed. support(), support0(), support1() keep their signatures — Python bindings are preserved. MinkowskiDiff lives in coal::details, signaling internal use.

getSupport free function unchanged. Signature preserved, implementation simplified. Also in coal::details.

Removed from public headers:

• Per-shape getShapeSupport<T> declarations from support_functions.h
• ContactPatchSolver::SupportSetFunction typedef and makeSupportSetFunction
• geometric_shapes_traits.h (entire file)

All removed items are in coal::details or internal headers. Code that depends on them is reaching into implementation details.

Implementing a custom shape outside of Coal

A custom shape subclasses ShapeBase, overrides virtuals, and works with Coal's collision/distance pipeline without modifying any Coal source files.

Minimal example

#include <coal/shape/geometric_shapes.h>

/// A rounded box: a box with filleted edges (Minkowski sum of Box and Sphere).
class RoundedBox : public coal::ShapeBase {
 public:
  coal::Vec3s halfSide;
  coal::Scalar fillet_radius;

  RoundedBox(const coal::Vec3s& halfSide, coal::Scalar fillet_radius)
      : halfSide(halfSide), fillet_radius(fillet_radius) {}

  // --- Required overrides (inherited pure virtuals) ---

  RoundedBox* clone() const override { return new RoundedBox(*this); }
  void computeLocalAABB() override {
    coal::Vec3s inflated = halfSide + coal::Vec3s::Constant(fillet_radius);
    aabb_local = coal::AABB(-inflated, inflated);
    aabb_center = aabb_local.center();
    aabb_radius = (inflated).norm();
  }
  coal::NODE_TYPE getNodeType() const override { return coal::GEOM_CUSTOM; }
  bool isEqual(const coal::CollisionGeometry& other) const override {
    const auto* o = dynamic_cast<const RoundedBox*>(&other);
    return o && o->halfSide == halfSide && o->fillet_radius == fillet_radius;
  }

  // --- Support function ---

  void getShapeSupport(const coal::Vec3s& dir, coal::Vec3s& support,
                       int& /*hint*/,
                       coal::details::ShapeSupportData& /*data*/) const final {
    // Box support: per-axis sign selection
    support = coal::Vec3s(
        dir[0] > 0 ? halfSide[0] : -halfSide[0],
        dir[1] > 0 ? halfSide[1] : -halfSide[1],
        dir[2] > 0 ? halfSide[2] : -halfSide[2]);
    // No fillet inflation here — handled by getInflationRadius()
  }

  coal::Scalar getInflationRadius() const final {
    return fillet_radius + getSweptSphereRadius();
  }
};

What to override

Virtual	Required?	Purpose
getShapeSupport	Yes (pure virtual)	Core support point without inflation.
getNodeType	Yes (pure virtual)	Return GEOM_CUSTOM.
clone	Yes (pure virtual)	Deep copy.
computeLocalAABB	Yes (pure virtual)	Bounding box for broadphase.
isEqual	Yes (pure virtual)	Equality comparison.
getInflationRadius	If core shape is smaller than actual geometry	Total radius for GJK post-correction. Default: getSweptSphereRadius().
initSupportData	If shape needs warm-start data	Pre-allocate scratch space. Default: no-op.
needNesterovNormalizeHeuristic	If shape is a polytope	Nesterov acceleration flag. Default: false.
isStrictlyConvex	If shape has unique support in every direction	Skips contact patch support set. Default: false.
getShapeSupportSet	If shape has flat faces	Multi-point contact patches. Default: single-point fallback.

Contrast with the current GEOM_CUSTOM fork

The current fork (13 commits on top of 14b40842f) adds custom shape support by bolting GEOM_CUSTOM onto the existing dispatch machinery. The virtual dispatch refactoring makes this machinery unnecessary.

Aspect	Current GEOM_CUSTOM fork	Virtual dispatch
Approach	Adds extension point to existing machinery	Replaces the machinery
Coal modifications	~15 files: adds GEOM_CUSTOM cases to 6 switches, 3 matrices, NODE_TYPE, ShapeBase	~10 files: moves support bodies to shapes, removes switches and dispatch machinery
Dispatch path for custom shapes	Function pointer → normalization wrapper (sqrt + 3 div) → virtual call	Direct virtual call (same path as built-in shapes)
dir contract for custom shapes	Guaranteed unit-length (wrapper normalizes)	Raw direction (same as built-in shapes receive)
Indirections per custom shape per iteration	2 (function pointer + virtual)	1 (virtual)
Contact patches for custom shapes	Single-point only (no override point for support sets)	Override getShapeSupportSet for multi-point patches
Composition (e.g., swept volume wrapping Box)	Re-enters NODE_TYPE switch for inner shape	Calls inner->getShapeSupport() directly
Built-in shape dispatch	Unchanged (template specializations)	Changed (virtual methods)

Both approaches require significant modifications to Coal. The difference is structural: the current fork adds a second-class extension point alongside the existing dispatch machinery, leaving the machinery itself intact. The virtual dispatch refactoring replaces the dispatch machinery, making extensibility a consequence of the cleaner internal structure.

What works automatically

Once the shape class exists, it participates in Coal's full pipeline without registration:

• GJK collision and distance — MinkowskiDiff::set() calls the virtuals; no matrix entry needed for GJK pairs.
• Contact patches — getShapeSupportSet is called directly; default produces a single contact point.
• Broadphase — works via computeLocalAABB().
• Composability — a custom shape can delegate to another shape's getShapeSupport (e.g., a swept volume wrapping a Box).

[lmontaut]

Hi @rjoomen, I think it's a bad idea to use virtualization for the support functions. The reason why coal is fast is that it efficiently dispatches the support function. The tag system in coal is ugly, I agree. But the support function calls need to happen as fast as possible, as they are called many times during collision detection.

Currently, when the GJK solver is created, it automatically resolves the exact support function of the minkowski difference it needs to run. When the solver runs, it never goes through the virtual table of the shapes.

[rjoomen]

That makes sense, although, as you can see above, question is how much this actually will hurt performance. I started this discussion just to pose an alternative solution to the GEOM_CUSTOM PR. I'd be happy with either solution, as long as Coal is extensible.