matryoshka

matryoshka(V, F, VB=None, FB=None, optimize='all', R=None, c=None, cut_point=None, cut_normal=None, a_plus=None, a_minus=None, n_samples=200, n_particles=30, max_iter=40, scale_tol=0.001, warm_start=False, verbose=False, seed=None)

Generalized Matryoshka: find a similarity transform of B that nests inside A such that A can be cut by a plane and pulled apart along a+/a- without colliding with the inner copy. Implements the algorithm of Jacobson (SGP 2017) on the CPU using particle swarm optimization.

Parameters:

Name	Type	Description	Default
V	(n,3) numpy double array	Vertex positions of the outer mesh A.	required
F	(m,3) numpy int array	Triangle indices of A.	required
VB	optional inner mesh B. If None, self-nesting is performed (B = A).		None
FB	optional inner mesh B. If None, self-nesting is performed (B = A).		None
optimize	str, optional (default 'all')	Which variables to optimize: * 'all' — scale, rotation, centroid, cut plane, removal directions * 'rigid' — scale, rotation, centroid (cut plane and removal directions fixed) * 'scale_only' — only the scale (everything else must be provided)	'all'
R	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
c	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
cut_point	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
cut_normal	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
a_plus	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
a_minus	optional fixed values used	either as defaults for non-optimized variables or, in 'scale_only' mode, as the entire configuration. If cut_point / cut_normal are not provided, a horizontal cut through the centroid of A is used.	None
n_samples	int, optional (default 200)	Number of random surface samples drawn from B (in addition to B's vertices) for the feasibility test.	200
n_particles	int	Particle swarm hyperparameters (kept small by default since each feasibility evaluation is expensive on the CPU).	30
max_iter	int	Particle swarm hyperparameters (kept small by default since each feasibility evaluation is expensive on the CPU).	30
scale_tol	float, optional (default 1e-3)	Binary-search tolerance for the inner scale search.	0.001
warm_start	bool or dict, optional (default False)	Only relevant for optimize='all'. If True, run a 'rigid' optimization first (with the same budget) and use its solution to seed the global-best tracker of the full optimization. Alternatively, pass a result dict from a previous matryoshka(...) call to use that as the seed (useful when you want to compare against a specific baseline). Either form guarantees that the returned scale is at least as large as the seed scale (modulo scale_tol).	False
verbose	bool, optional (default False)	Print particle-swarm progress.	False
seed	int or None	Seed for the surface-sampling RNG and the particle-swarm RNG.	None

Returns:

Name	Type	Description
result	dict with keys	s : float, the optimal scale. R : (3,3) rotation matrix applied to B. c : (3,) translation, the new centroid of T(B). B_center : (3,) original centroid of B used as the rotation pivot. cut_point : (3,) point on the cut plane. cut_normal : (3,) unit normal of the cut plane. a_plus : (3,) removal direction for A above the plane. a_minus : (3,) removal direction for A below the plane.

Notes

This is a CPU implementation of a method originally formulated with GPU depth peeling. It is therefore much slower than the original. Use modest n_samples, n_particles and max_iter for interactive experimentation.

Examples:

Inspect a result interactively with polyscope. The transformed inner copy is T(B) = c + s · R · (B − B_center); the cut plane is rendered as a square in the normal's tangent frame, and the two removal directions are drawn as a curve network:

import polyscope as ps
import numpy as np
import gpytoolbox as gpy

V, F = gpy.read_mesh("bunny.obj")
V = V - V.mean(0); V = V / np.max(np.abs(V))
res = gpy.matryoshka(V, F, optimize='rigid',
                     n_samples=80, n_particles=20, max_iter=20, seed=0)

# Inner copy.
T_B = (res['s'] * (V - res['B_center']) @ res['R'].T) + res['c']

# Cut plane as a square in the tangent frame of cut_normal.
diag = float(np.linalg.norm(V.max(0) - V.min(0)))
n = res['cut_normal']
e = np.array([1.,0.,0.]) if abs(n[0])<0.9 else np.array([0.,1.,0.])
t1 = e - np.dot(e, n) * n; t1 /= np.linalg.norm(t1)
t2 = np.cross(n, t1)
h, p0 = 0.75*diag, res['cut_point']
quad_V = np.stack([p0 + a*h*t1 + b*h*t2
                   for (a,b) in [(-1,-1),(1,-1),(1,1),(-1,1)]])
quad_F = np.array([[0,1,2],[0,2,3]], dtype=np.int32)

# Removal directions as a 2-edge curve network.
arrows = np.stack([p0, p0+0.5*diag*res['a_plus'],
                   p0, p0+0.5*diag*res['a_minus']])
edges = np.array([[0,1],[2,3]], dtype=np.int32)

ps.init()
ps.register_surface_mesh("A", V, F, transparency=0.35)
ps.register_surface_mesh("T(B)", T_B, F)
ps.register_surface_mesh("cut plane", quad_V, quad_F, transparency=0.4)
ps.register_curve_network("removal dirs", arrows, edges, radius=0.005)
ps.show()

Source code in src/gpytoolbox/matryoshka.py

def matryoshka(V, F,
               VB=None, FB=None,
               optimize='all',
               R=None, c=None,
               cut_point=None, cut_normal=None,
               a_plus=None, a_minus=None,
               n_samples=200,
               n_particles=30, max_iter=40,
               scale_tol=1e-3,
               warm_start=False,
               verbose=False,
               seed=None):
    """Generalized Matryoshka: find a similarity transform of `B` that nests
    inside `A` such that `A` can be cut by a plane and pulled apart along
    `a+`/`a-` without colliding with the inner copy. Implements the algorithm
    of Jacobson (SGP 2017) on the CPU using particle swarm optimization.

    Parameters
    ----------
    V : (n,3) numpy double array
        Vertex positions of the outer mesh A.
    F : (m,3) numpy int array
        Triangle indices of A.
    VB, FB : optional inner mesh B. If None, self-nesting is performed (B = A).
    optimize : str, optional (default 'all')
        Which variables to optimize:
        * 'all'        — scale, rotation, centroid, cut plane, removal directions
        * 'rigid'      — scale, rotation, centroid (cut plane and removal directions fixed)
        * 'scale_only' — only the scale (everything else must be provided)
    R, c, cut_point, cut_normal, a_plus, a_minus : optional fixed values used
        either as defaults for non-optimized variables or, in 'scale_only' mode,
        as the entire configuration. If `cut_point` / `cut_normal` are not
        provided, a horizontal cut through the centroid of A is used.
    n_samples : int, optional (default 200)
        Number of random surface samples drawn from B (in addition to B's
        vertices) for the feasibility test.
    n_particles, max_iter : int, optional
        Particle swarm hyperparameters (kept small by default since each
        feasibility evaluation is expensive on the CPU).
    scale_tol : float, optional (default 1e-3)
        Binary-search tolerance for the inner scale search.
    warm_start : bool or dict, optional (default False)
        Only relevant for `optimize='all'`. If True, run a `'rigid'`
        optimization first (with the same budget) and use its solution to
        seed the global-best tracker of the full optimization. Alternatively,
        pass a result dict from a previous `matryoshka(...)` call to use
        that as the seed (useful when you want to compare against a
        specific baseline). Either form guarantees that the returned scale
        is at least as large as the seed scale (modulo `scale_tol`).
    verbose : bool, optional (default False)
        Print particle-swarm progress.
    seed : int or None, optional
        Seed for the surface-sampling RNG and the particle-swarm RNG.

    Returns
    -------
    result : dict with keys
        s          : float, the optimal scale.
        R          : (3,3) rotation matrix applied to B.
        c          : (3,) translation, the new centroid of T(B).
        B_center   : (3,) original centroid of B used as the rotation pivot.
        cut_point  : (3,) point on the cut plane.
        cut_normal : (3,) unit normal of the cut plane.
        a_plus     : (3,) removal direction for A above the plane.
        a_minus    : (3,) removal direction for A below the plane.

    Notes
    -----
    This is a CPU implementation of a method originally formulated with GPU
    depth peeling. It is therefore much slower than the original. Use modest
    `n_samples`, `n_particles` and `max_iter` for interactive experimentation.

    Examples
    --------
    Inspect a result interactively with polyscope. The transformed inner
    copy is `T(B) = c + s · R · (B − B_center)`; the cut plane is rendered
    as a square in the normal's tangent frame, and the two removal
    directions are drawn as a curve network:
    ```python
    import polyscope as ps
    import numpy as np
    import gpytoolbox as gpy

    V, F = gpy.read_mesh("bunny.obj")
    V = V - V.mean(0); V = V / np.max(np.abs(V))
    res = gpy.matryoshka(V, F, optimize='rigid',
                         n_samples=80, n_particles=20, max_iter=20, seed=0)

    # Inner copy.
    T_B = (res['s'] * (V - res['B_center']) @ res['R'].T) + res['c']

    # Cut plane as a square in the tangent frame of cut_normal.
    diag = float(np.linalg.norm(V.max(0) - V.min(0)))
    n = res['cut_normal']
    e = np.array([1.,0.,0.]) if abs(n[0])<0.9 else np.array([0.,1.,0.])
    t1 = e - np.dot(e, n) * n; t1 /= np.linalg.norm(t1)
    t2 = np.cross(n, t1)
    h, p0 = 0.75*diag, res['cut_point']
    quad_V = np.stack([p0 + a*h*t1 + b*h*t2
                       for (a,b) in [(-1,-1),(1,-1),(1,1),(-1,1)]])
    quad_F = np.array([[0,1,2],[0,2,3]], dtype=np.int32)

    # Removal directions as a 2-edge curve network.
    arrows = np.stack([p0, p0+0.5*diag*res['a_plus'],
                       p0, p0+0.5*diag*res['a_minus']])
    edges = np.array([[0,1],[2,3]], dtype=np.int32)

    ps.init()
    ps.register_surface_mesh("A", V, F, transparency=0.35)
    ps.register_surface_mesh("T(B)", T_B, F)
    ps.register_surface_mesh("cut plane", quad_V, quad_F, transparency=0.4)
    ps.register_curve_network("removal dirs", arrows, edges, radius=0.005)
    ps.show()
    ```
    """
    V = np.asarray(V, dtype=np.float64)
    F = np.asarray(F, dtype=np.int32)
    if VB is None:
        VB = V
        FB = F
    VB = np.asarray(VB, dtype=np.float64)
    FB = np.asarray(FB, dtype=np.int32)

    rng = np.random.default_rng(seed)

    B_center = VB.mean(axis=0)
    A_center = V.mean(axis=0)
    A_min = V.min(axis=0)
    A_max = V.max(axis=0)
    A_extent = A_max - A_min

    # Defaults for fixed values (used when optimize != 'all' or values not
    # provided).
    if cut_normal is None:
        cut_normal = np.array([0.0, 0.0, 1.0])
    else:
        cut_normal = _normalize(np.asarray(cut_normal, dtype=np.float64))
    if cut_point is None:
        cut_point = A_center.copy()
    else:
        cut_point = np.asarray(cut_point, dtype=np.float64)
    if a_plus is None:
        a_plus = cut_normal.copy()
    else:
        a_plus = _normalize(np.asarray(a_plus, dtype=np.float64))
    if a_minus is None:
        a_minus = -cut_normal.copy()
    else:
        a_minus = _normalize(np.asarray(a_minus, dtype=np.float64))

    # Surface samples of B (canonical, before transformation).
    samples_B = _sample_B_surface(VB, FB, n_samples, rng)

    # Build the AABB tree + fast winding number BVH for A once and reuse them
    # across every feasibility evaluation. This is the single biggest CPU win
    # for the optimization — `signed_distance` is called hundreds of times
    # and would otherwise rebuild both trees on each call.
    A_aabb = squared_distance_precompute(V, F)
    A_fwn_bvh = fast_winding_number_precompute(V, F)
    A_rmi = ray_mesh_intersect_precompute(V, F)

    if optimize == 'scale_only':
        if R is None:
            R = np.eye(3)
        if c is None:
            c = A_center.copy()
        s = _largest_feasible_scale(samples_B, B_center, V, F, R, c,
                                    cut_point, cut_normal, a_plus, a_minus,
                                    cpp_aabb=A_aabb, fwn_bvh=A_fwn_bvh,
                                    intersector=A_rmi,
                                    tol=scale_tol)
        return dict(s=s, R=R, c=c, B_center=B_center,
                    cut_point=cut_point, cut_normal=cut_normal,
                    a_plus=a_plus, a_minus=a_minus)

    # Build (lb, ub) for the free parameters.
    # Always-free: 3 axis-angle + 3 centroid
    lb = [-np.pi, -np.pi, -np.pi,
          A_min[0], A_min[1], A_min[2]]
    ub = [np.pi, np.pi, np.pi,
          A_max[0], A_max[1], A_max[2]]
    fixed = dict(cut_point=cut_point, cut_normal=cut_normal,
                 a_plus=a_plus, a_minus=a_minus,
                 cut_anchor=A_center.copy())

    if optimize == 'all':
        # cut normal direction (3 free), offset along normal (1), a+ (3), a- (3)
        max_offset = 0.5 * np.linalg.norm(A_extent)
        lb += [-1.0, -1.0, -1.0, -max_offset, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0]
        ub += [1.0, 1.0, 1.0, max_offset, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

    lb = np.array(lb, dtype=np.float64)
    ub = np.array(ub, dtype=np.float64)

    # Use a deterministic rng feeding into particle_swarm via numpy global seed
    # (particle_swarm uses np.random internally for the Python path; the C++
    # path uses a non-seeded RNG, which is fine — the objective itself does
    # not depend on the swarm's RNG for reproducibility of the *result*).
    if seed is not None:
        np.random.seed(seed)

    # Track the best configuration encountered (the swarm only tracks the
    # best objective; we need the full configuration to return).
    best = {'s': 0.0, 'x': None}

    # Optional warm-start: seed the global-best tracker with a previously
    # computed configuration. The swarm itself starts from random particles,
    # but our external `best` tracker retains the seed whenever the swarm
    # doesn't beat it. The seed can be either a result dict from a prior
    # `matryoshka(...)` call, or True (run 'rigid' internally first).
    if warm_start and optimize == 'all':
        if isinstance(warm_start, dict):
            seed_res = warm_start
        else:
            seed_res = matryoshka(
                V, F, VB=VB, FB=FB, optimize='rigid',
                R=R, c=c,
                cut_point=cut_point, cut_normal=cut_normal,
                a_plus=a_plus, a_minus=a_minus,
                n_samples=n_samples,
                n_particles=n_particles, max_iter=max_iter,
                scale_tol=scale_tol, verbose=verbose, seed=seed)
        if seed_res['s'] > 0:
            seed_x = _encode_all(
                seed_res['R'], seed_res['c'],
                seed_res['cut_normal'], seed_res['cut_point'],
                seed_res['a_plus'], seed_res['a_minus'],
                fixed['cut_anchor'])
            # Re-evaluate the seed's scale with the outer call's samples_B,
            # which may differ from whatever sampling produced seed_res. This
            # keeps best['s'] self-consistent: it is always the maximum
            # feasible scale of best['x'] *as measured here*.
            R_s, c_s, cp_s, cn_s, ap_s, am_s = _decode(seed_x, 'all', fixed)
            best['s'] = _largest_feasible_scale(
                samples_B, B_center, V, F, R_s, c_s,
                cp_s, cn_s, ap_s, am_s,
                cpp_aabb=A_aabb, fwn_bvh=A_fwn_bvh, intersector=A_rmi,
                tol=scale_tol)
            best['x'] = seed_x

    def objective(x):
        R_x, c_x, cp, cn, ap, am = _decode(x, optimize, fixed)
        s = _largest_feasible_scale(samples_B, B_center, V, F, R_x, c_x,
                                    cp, cn, ap, am,
                                    cpp_aabb=A_aabb, fwn_bvh=A_fwn_bvh,
                                    intersector=A_rmi,
                                    tol=scale_tol)
        if s > best['s']:
            best['s'] = s
            best['x'] = x.copy()
        # We minimize: negative scale.
        return -s

    _x_best, _f_best = particle_swarm(
        objective, lb, ub,
        n_particles=n_particles, max_iter=max_iter,
        verbose=verbose, topology='full')

    if best['x'] is None:
        # Swarm never found a feasible scale; return zero-scale result.
        return dict(s=0.0,
                    R=np.eye(3), c=A_center.copy(), B_center=B_center,
                    cut_point=cut_point, cut_normal=cut_normal,
                    a_plus=a_plus, a_minus=a_minus)

    R_out, c_out, cp_out, cn_out, ap_out, am_out = _decode(
        best['x'], optimize, fixed)
    return dict(s=best['s'], R=R_out, c=c_out, B_center=B_center,
                cut_point=cp_out, cut_normal=cn_out,
                a_plus=ap_out, a_minus=am_out)