Advanced Topics

Who This Vignette Is For

Audience: Advanced users, researchers, algorithm developers, curious minds

Prerequisites:

• Familiarity with lap_solve() (vignette("getting-started"))
• Basic complexity analysis (Big-O notation)
• Interest in algorithm design or scientific computing

What You’ll Learn:

• Why exact LAP becomes infeasible for large n
• Three approximation strategies and their trade-offs
• How matching problems appear in ecology, physics, and chemistry
• Mathematical connections to optimal transport theory

Time to complete: 45-60 minutes (conceptual reading)

Documentation Roadmap

You are here: Pixel Morphing (Advanced)

Why Pixels?

Pixels provide an ideal testbed for understanding assignment problems:

• Each pixel is an entity with measurable properties
• Color = feature (what it looks like)
• Position = spatial location (where it is)
• The matching is visually verifiable: you can see if it worked

The same algorithms that morph images smoothly also track particles in physics, align molecules in chemistry, and match vegetation plots in ecology.

Scientific Applications at a Glance

Problem Formulation

Given two sets of entities \(A = \{a_1, \ldots, a_n\}\) and \(B = \{b_1, \ldots, b_n\}\), find the optimal one-to-one correspondence by minimizing

\[ \min_{\pi} \sum_{i=1}^{n} c_{i,\pi(i)} \]

where the cost combines feature similarity and spatial proximity:

\[ c_{ij} = \alpha \, d_{\text{feature}}(a_i, b_j) + \beta \, d_{\text{spatial}}(\mathbf{x}_i, \mathbf{x}_j). \]

Feature distance \(d_{\text{feature}}\): domain-specific similarity

• Ecology: Bray-Curtis dissimilarity between species vectors
• Physics: difference in particle intensity or size
  
• Chemistry: penalty for mismatched atom types
  
• Images: Euclidean distance in RGB color space

Spatial distance \(d_{\text{spatial}}\): physical proximity

• Ecology: geographic distance between plot centers
  
• Physics: Euclidean distance accounting for predicted motion
  
• Chemistry: 3D distance between atomic coordinates
  
• Images: 2D pixel position distance

Weights \(\alpha, \beta \ge 0\) balance feature matching vs. spatial coherence.

Computational Challenge

Exact solution: solve the full \(n \times n\) LAP.

• Complexity: \(O(n^3)\) using Jonker-Volgenant
• Feasible: up to \(n \approx 1000\) (about \(30 \times 30\) images, or \(1000\) plots/particles/atoms)
  
• Prohibitive: for \(n = 10\,000\) (\(100 \times 100\) images), runtime and memory become expensive

Real applications often involve

• High-resolution images: \(200 \times 200 = 40\,000\) pixels
  
• Large ecological surveys: \(5000+\) plots
  
• Particle tracking: \(10\,000+\) particles per frame
  
• Molecular dynamics: \(100\,000+\) atoms

We therefore need approximations that are much faster but still produce high-quality matchings.

Exact Pixel Matching

The exact pixel morph uses a full LAP solution on a \(1600 \times 1600\) cost matrix. For each pair of pixels \((i, j)\) we compute

\[ c_{ij} = \alpha \,\lVert \text{RGB}_i^A - \text{RGB}_j^B \rVert_2 + \beta \,\lVert (x_i, y_i) - (x_j, y_j) \rVert_2, \]

where color distances are normalized to \([0, \sqrt{3}]\) (RGB in \([0,1]\)) and spatial distances to \([0,1]\) using the image diagonal.

This yields an optimal one-to-one assignment of pixels. The resulting animations are smooth and artifact-free but require solving the full LAP.

Feature Quantization Morph (Strategy 1)

In the feature quantization morph, similar colors are grouped, and groups are matched rather than individual pixels. Colors move as coherent “bands,” preserving global color structure but losing fine-grained per-pixel detail.

Hierarchical Morph (Strategy 2)

The hierarchical morph first matches large patches, then refines within patches. The motion is locally coherent and scales well to large problems, at the price of potentially missing some globally optimal cross-patch matches.

Strategy 1: Feature Quantization

Core idea: reduce problem size by grouping entities with similar features, then match groups.

Strategy 2: Hierarchical Decomposition

Core idea: split the domain into smaller subproblems by spatial partitioning, solve subproblems, and combine.

// Pseudocode for hierarchical LAP matching

FUNCTION match_hierarchical(region_A, region_B, threshold, level):

  // Base case: region small enough for exact LAP
  IF size(region_A) <= threshold THEN
    cost ← compute_cost_matrix(region_A, region_B, α, β)
    RETURN lap_solve(cost)
  END IF

  // Divide into 2×2 spatial grid (4 patches)
  patches_A ← spatial_partition(region_A, grid = 2×2)
  patches_B ← spatial_partition(region_B, grid = 2×2)

  // Compute patch representatives
  FOR each patch p DO
    centroid[p]      ← mean(positions in p)
    mean_feature[p]  ← mean(features in p)
  END FOR

  // Match patches using 4×4 LAP
  patch_cost ← matrix(4, 4)
  FOR i = 1 TO 4 DO
    FOR j = 1 TO 4 DO
      patch_cost[i, j] ← α·distance(mean_feature_A[i], mean_feature_B[j]) +
                         β·distance(centroid_A[i], centroid_B[j])
    END FOR
  END FOR

  patch_assignment ← lap_solve(patch_cost)

  // Recursively solve within matched patches
  assignments ← []
  FOR i = 1 TO 4 DO
    j ← patch_assignment[i]
    sub_assignment ← match_hierarchical(
      patches_A[i],
      patches_B[j],
      threshold,
      level + 1
    )
    assignments ← append(assignments, sub_assignment)
  END FOR

  RETURN concatenate(assignments)
END FUNCTION

Strategy 3: Resolution Reduction

Core idea: solve the LAP on a coarse grid, then lift/upscale the assignment to the full-resolution grid.

Strategy Comparison

Practical rules of thumb

• \(n < 1000\): use the exact LAP.
  
• \(1000 < n < 5000\): feature quantization or a shallow hierarchy.
  
• \(n > 5000\): hierarchical decomposition with 2-3 levels.
• \(n > 50\,000\): combine \(s = 2\) resolution reduction with a hierarchical method.

Ecology: Vegetation Plot Matching

Problem: match \(n\) vegetation plots surveyed at time \(t\) to \(n\) plots at time \(t + \Delta t\) to track community dynamics.

Feature distance: Bray-Curtis dissimilarity between species abundance vectors

\[ d_{\text{BC}}(a, b) = \frac{\sum_s \lvert a_s - b_s \rvert} {\sum_s (a_s + b_s)}, \]

where \(a_s, b_s\) are abundances of species \(s\) in plots \(a\) and \(b\).

Spatial distance: geographic distance (e.g. in kilometers) between plot centers.

Exact solution for small studies (\(n < 100\)):

// Pseudocode for ecological plot matching
FOR i = 1 TO n_plots_t DO
  FOR j = 1 TO n_plots_tplus DO

    // Bray-Curtis dissimilarity for species composition
    numerator   ← sum over species s of |abundance_t[i, s] - abundance_tplus[j, s]|
    denominator ← sum over species s of (abundance_t[i, s] + abundance_tplus[j, s])
    bc_distance ← numerator / denominator

    // Geographic distance (kilometers)
    geo_distance ← sqrt((x_t[i] - x_tplus[j])^2 +
                        (y_t[i] - y_tplus[j])^2)

    // Combined cost (α = 0.7 emphasizes species composition)
    cost[i, j] ← 0.7 · bc_distance + 0.3 · (geo_distance / max_distance)

  END FOR
END FOR

plot_correspondence ← lap_solve(cost)

For large studies (\(n > 1000\)) a hierarchical approach by region is more practical:

// Hierarchical decomposition by geographic region

// 1. Divide landscape into spatial grid (e.g. 10 km × 10 km cells)
regions_t     ← spatial_partition(plots_t,     grid_size = 10 km)
regions_tplus ← spatial_partition(plots_tplus, grid_size = 10 km)

// 2. Compute region representatives
FOR each region r DO
  mean_composition[r] ← average species vector across plots in r
  centroid[r]         ← geographic center of r
END FOR

// 3. Match regions (small LAP: ~100 regions)
region_cost       ← compute_cost(mean_composition, centroids, α = 0.7, β = 0.3)
region_assignment ← lap_solve(region_cost)

// 4. Within matched regions, solve plot-level LAP
full_assignment ← []
FOR r = 1 TO n_regions DO
  r_matched ← region_assignment[r]
  plots_A   ← plots in regions_t[r]
  plots_B   ← plots in regions_tplus[r_matched]

  // Local LAP (smaller problem, e.g. 50 × 50)
  cost_local       ← compute_plot_cost(plots_A, plots_B, α = 0.7, β = 0.3)
  local_assignment ← lap_solve(cost_local)

  full_assignment ← append(full_assignment, local_assignment)
END FOR

RETURN full_assignment

This allows tracking individual plot trajectories across time, distinguishing stable communities, successional trends, and invasion fronts.

Physics: Particle Tracking

Problem: track \(n\) particles between frame \(t\) and \(t + \Delta t\) in experimental video.

Feature distance: differences in intensity, size, or shape.

Spatial distance: displacement relative to predicted motion:

\[ d_{\text{spatial}}(i, j) = \bigl\| \mathbf{x}_i + \mathbf{v}_i \Delta t - \mathbf{x}_j \bigr\|_2, \]

where \(\mathbf{v}_i\) is the estimated velocity from previous frames.

We also impose a maximum displacement \(d_{\max}\) beyond which matches are physically implausible.

Exact solution (moderate \(n\)):

// Pseudocode for particle tracking with velocity prediction

// Initialize cost matrix as forbidden everywhere
cost ← matrix(Inf, n_particles_t, n_particles_tplus)

FOR i = 1 TO n_particles_t DO
  // Predict position using previous velocity
  x_predicted ← x_t[i] + v_x_t[i] · Δt
  y_predicted ← y_t[i] + v_y_t[i] · Δt

  FOR j = 1 TO n_particles_tplus DO
    // Distance from predicted position
    dx ← x_predicted - x_tplus[j]
    dy ← y_predicted - y_tplus[j]
    spatial_distance ← sqrt(dx^2 + dy^2)

    // Only consider physically plausible matches
    IF spatial_distance <= max_displacement THEN
      // Feature similarity (intensity, size, etc.)
      feature_distance ← |intensity_t[i] - intensity_tplus[j]|

      // Combined cost
      cost[i, j] ← α · feature_distance + β · spatial_distance
    END IF
  END FOR
END FOR

// Solve assignment (Inf entries are forbidden)
particle_tracks ← lap_solve(cost)

// Update velocities from assignments
FOR i = 1 TO n_particles_t DO
  j ← particle_tracks[i]
  velocity_new[i] ← (position_tplus[j] - position_t[i]) / Δt
END FOR

For dense tracking (\(n > 5000\)), we can first cluster particles:

// Two-stage: clustering then local matching

// Stage 1: spatial clustering
clusters_t     ← spatial_cluster(particles_t,     radius = 2 · pixel_size)
clusters_tplus ← spatial_cluster(particles_tplus, radius = 2 · pixel_size)

// Compute cluster representatives
FOR each cluster c DO
  centroid[c]       ← mean position of particles in c
  mean_intensity[c] ← mean intensity
  mean_velocity[c]  ← mean velocity (if available)
END FOR

// Match clusters
cluster_cost   ← compute_cluster_similarity(clusters_t, clusters_tplus)
cluster_tracks ← lap_solve(cluster_cost)

// Stage 2: within matched clusters, track individual particles
full_tracks ← []
FOR c = 1 TO n_clusters DO
  c_matched   ← cluster_tracks[c]
  particles_A ← particles in clusters_t[c]
  particles_B ← particles in clusters_tplus[c_matched]

  cost_local ← compute_particle_distance(
    particles_A, particles_B,
    max_displacement = 5,
    α = 0.3,
    β = 0.7
  )

  local_tracks ← lap_solve(cost_local)
  full_tracks  ← append(full_tracks, local_tracks)
END FOR

RETURN full_tracks

This yields efficient and robust trajectories even for very dense particle fields.

Chemistry: Molecular Conformation Alignment

Problem: align two conformations of the same molecule (e.g. a protein) with \(n\) atoms to compute RMSD and analyze structural change.

Feature distance: strict element matching

\[ d_{\text{element}}(i, j) = \begin{cases} 0, & \text{if } \text{element}_i = \text{element}_j, \\ \infty, & \text{otherwise.} \end{cases} \]

Spatial distance: 3D Euclidean distance between atomic coordinates.

Exact LAP for small molecules:

// Pseudocode for molecular conformation alignment

n_atoms ← number of atoms in molecule
cost    ← matrix(0, n_atoms, n_atoms)

FOR i = 1 TO n_atoms DO
  FOR j = 1 TO n_atoms DO

    // Enforce strict element type matching
    IF element_type_A[i] ≠ element_type_B[j] THEN
      cost[i, j] ← Inf
    ELSE
      dx ← x_A[i] - x_B[j]
      dy ← y_A[i] - y_B[j]
      dz ← z_A[i] - z_B[j]
      cost[i, j] ← sqrt(dx^2 + dy^2 + dz^2)
    END IF

  END FOR
END FOR

// Solve alignment
alignment ← lap_solve(cost)

// Compute RMSD
sum_sq_dist ← 0
FOR i = 1 TO n_atoms DO
  j            ← alignment[i]
  sum_sq_dist ← sum_sq_dist + cost[i, j]^2
END FOR

rmsd ← sqrt(sum_sq_dist / n_atoms)

For large biomolecules, we again use a hierarchical strategy, this time by secondary structure elements (helices, sheets, loops, etc.), aligning segments first and then atoms within matched segments.

Customizing Morph Duration

The morphing examples use default settings, but you can customize the number of frames and speed:

# From inst/scripts/generate_examples.R
generate_morph <- function(assignment, pixels_A, pixels_B,
                           n_frames    = 30,   # number of frames
                           frame_delay = 0.1)  # delay between frames (seconds)
{
  frames <- lapply(seq(0, 1, length.out = n_frames), function(t) {
    interpolate_frame(t, assignment, pixels_A, pixels_B)
  })

  save_gif(frames, delay = frame_delay)
}

Total animation duration is n_frames * frame_delay seconds.

Using the Example Code

The morphing implementation is provided in inst/scripts/generate_examples.R:

# View the source
example_script <- system.file("scripts", "generate_examples.R", package = "couplr")
file.show(example_script)

# Or source to use its helpers
source(example_script)

# Apply to your own data
my_cost       <- build_cost_matrix(my_data_A, my_data_B)
my_assignment <- lap_solve(my_cost)

Regenerating Examples

To regenerate all demo GIFs and PNGs:

source("inst/scripts/generate_examples.R")

This will write the assets under inst/extdata.

Monge Problem

The original Monge formulation (1781) seeks a transport map \(T: A \to B\) minimizing

\[ \int_A c(\mathbf{x}, T(\mathbf{x})) \,\mathrm{d}\mu(\mathbf{x}). \]

Kantorovich Relaxation

Kantorovich (1942) relaxed this to a transport plan \(\gamma\) on \(A \times B\):

\[ \min_{\gamma} \int_{A \times B} c(\mathbf{x}, \mathbf{y}) \,\mathrm{d}\gamma(\mathbf{x}, \mathbf{y}) \]

subject to marginal constraints on \(\gamma\).

Discrete Linear Assignment

For discrete uniform distributions with \(n\) points in \(A\) and \(B\) we obtain exactly the linear assignment problem:

\[ \min_{\pi \in S_n} \sum_{i=1}^n c_{i,\pi(i)}, \]

which is what couplr solves efficiently.

Wasserstein Distance

With \(c_{ij} = d(\mathbf{x}_i, \mathbf{x}_j)\) (often Euclidean distance), the optimal cost defines the \(1\)‑Wasserstein distance:

\[ W_1(\mu, \nu) = \min_{\pi \in S_n} \sum_{i=1}^n c_{i,\pi(i)}. \]

This appears in

• Earth mover’s distance for image retrieval
  
• Distributional similarity in statistics
  
• Generative modeling (e.g. Wasserstein GANs)

Optimal Transport Theory

• Peyré, G., & Cuturi, M. (2019). Computational Optimal Transport. Foundations and Trends in Machine Learning.
  
• Villani, C. (2008). Optimal Transport: Old and New. Springer.

Scientific Applications

Ecology

• Anderson, M. J. et al. (2011). Navigating the multiple meanings of beta diversity. Ecology Letters.
  
• Legendre, P., & Legendre, L. (2012). Numerical Ecology. Elsevier.

Physics

• Adrian, R. J., & Westerweel, J. (2011). Particle Image Velocimetry. Cambridge University Press.
  
• Crocker, J. C., & Grier, D. G. (1996). Methods of digital video microscopy. Journal of Colloid and Interface Science.

Chemistry

• Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Crystallographica.
  
• Coutsias, E. A. et al. (2004). Using quaternions to calculate RMSD. Journal of Computational Chemistry.

Assignment Algorithms

• Burkard, R., Dell’Amico, M., & Martello, S. (2009). Assignment Problems. SIAM.
  
• For implementation details in this package see vignette("algorithms").

When Do These Strategies Apply?

Practical Recommendations

For n < 3,000: Use match_couples() with exact algorithms:

result <- match_couples(left, right, vars = c("x", "y", "z"), auto_scale = TRUE)

For 3,000 < n < 10,000: Use blocking to create smaller subproblems:

blocks <- matchmaker(left, right, block_type = "cluster", n_blocks = 10)
result <- match_couples(blocks$left, blocks$right, vars = vars, block_id = "block_id")

For n > 10,000: Use greedy matching:

result <- greedy_couples(left, right, vars = vars, strategy = "sorted")

For n > 50,000: Combine strategies (blocking + greedy within blocks), or implement custom approximations using the techniques in this vignette.

Limitations of Approximation Strategies

Each approximation trades accuracy for speed. Know the failure modes: