@page edge_contraction Edge Contraction

# Edge Contraction Simplification for Dendritic Trees

## Overview

This document describes the implementation of a quadric error metric-based edge contraction algorithm for simplifying dendritic tree structures. The algorithm is inspired by Garland-Heckbert mesh simplification but adapted specifically for neuronal morphology trees.

## Algorithm Description

### Core Concept

The algorithm simplifies dendritic trees by iteratively contracting edges while minimizing geometric distortion, measured using quadric error metrics. Each vertex has an associated **quadric form** `f_v(x) = x^T · Q_v · x` that measures geometric error at position `x`.

### Vertex Classification

Vertices are classified by their **degree** (number of incident edges):

- **Leaf** (degree 1): Terminal vertices with no children
- **Regular** (degree 2): Linear continuation vertices with one parent and one child
- **Fork** (degree 3+): Branching points with multiple children

### Edge Contraction Operation

Contracting an edge `(a,b)` performs the following:

1. Removes the edge connecting vertices `a` and `b`
2. Replaces `a` and `b` with a new vertex `c`
3. Reconnects all edges incident to `a` or `b` to the new vertex `c`
4. Updates the quadric: `Q_c = Q_a + Q_b`
5. Updates the incident edge set: `pre(c) = pre(a) ∪ pre(b)`

**Topological Constraint**: An edge can only be contracted if **at least one endpoint is a regular vertex**. This preserves the topological type of the tree.

### Optimal Vertex Position

The position of the new vertex `c` is chosen to minimize geometric distortion:

- **Both endpoints regular**: `c` minimizes `f_a(c) + f_b(c)` (currently uses midpoint heuristic)
- **Only one endpoint regular**: `c` equals the position of the non-regular vertex

### Distortion Metric

The **distortion** caused by contracting edge `(a,b)` to position `c` is:

```
distortion(c) = f_a(c) + f_b(c)
```

where `f_v(x) = x^T · Q_v · x` evaluates the quadric error for vertex `v` at point `x`.

## Implementation Details

### Core Structures

#### `EdgeContractionSimplifier`

The main simplification engine that manages:
- Current set of vertices and edges
- Priority queue (heap) of edges sorted by distortion
- Incident edge sets for each vertex (`pre(v)`)
- Distortion threshold for halting

#### `EdgeWithDistortion`

Priority queue entry containing:
- Edge ID (hash-based, order-independent)
- Distortion value
- Vertex IDs of endpoints

Implements `Ord` with reverse ordering so minimum distortion edges have highest priority.

#### `Branch`

Represents a maximal sequence of edges connecting two non-regular vertices through only regular vertices.

### Key Algorithms

#### 1. Initialization (`initialize_heap`)

```rust
fn initialize_heap(&mut self)
```

- Builds incident edge sets for all vertices
- Computes distortion for each edge
- Populates priority queue with all edges

#### 2. Distortion Calculation (`compute_distortion`)

```rust
fn compute_distortion(&self, vertex_a: &Vertex, vertex_b: &Vertex) -> (f64, Vector3D)
```

Determines optimal contraction position based on vertex types and computes distortion value.

#### 3. Edge Contraction (`contract_edge`)

```rust
fn contract_edge(&mut self, edge: &Edge, new_position: Vector3D) -> u128
```

Performs the actual edge contraction:
1. Creates new vertex `c` at optimal position
2. Computes `Q_c = Q_a + Q_b`
3. Builds `pre(c) = pre(a) ∪ pre(b)`
4. Reconnects all incident edges to `c`
5. Removes vertices `a` and `b`
6. Returns ID of new vertex

#### 4. Main Loop (`simplify`)

```rust
pub fn simplify(&mut self) -> Dendrite
```

Executes the iterative simplification:

1. Initialize heap with all edges
2. While heap is not empty:
   - Pop minimum distortion edge
   - Check if distortion exceeds threshold --> halt
   - Verify at least one endpoint is regular --> skip if not
   - Compute optimal position
   - Contract edge
   - Update incident edge distortions
   - Check if regular vertices remain --> halt if none
3. Return simplified dendrite

#### 5. Branch Decomposition (`decompose_into_branches`)

```rust
pub fn decompose_into_branches(&self) -> Vec<Branch>
```

Identifies all branches in the simplified structure by tracing paths between non-regular vertices.

### Optimization Techniques

#### Lazy Heap Updates

The implementation uses a **lazy deletion** strategy for the priority queue:

- Old edge entries remain in the heap after modifications
- When popped, edges are checked for validity using `valid_edges` set
- Invalid entries are skipped
- New distortion values are added without removing old ones

This avoids expensive heap reorganization on every update.

#### Order-Independent Edge Hashing

Edges are identified by a hash computed from vertex IDs in canonical order:

```rust
fn compute_edge_id(v1: u128, v2: u128) -> u64
```

This ensures `Edge(a,b)` and `Edge(b,a)` have the same ID.

## Usage

### Basic Example

```rust
use quadrix::swc::load_swc_file;
use quadrix::simplification::EdgeContractionSimplifier;

// Load dendrite
let dendrite = load_swc_file("neuron.swc")?;

// Create simplifier with distortion threshold
let mut simplifier = EdgeContractionSimplifier::new(dendrite, 0.5);

// Simplify
let simplified = simplifier.simplify();

println!("Original: {} vertices, {} edges",
         dendrite.vertices.len(), dendrite.edges.len());
println!("Simplified: {} vertices, {} edges",
         simplified.vertices.len(), simplified.edges.len());
```

### Command-Line Example

```bash
# Run the example program
cargo run --example simplify_dendrite data/test_neuron.swc 0.5
```

### Branch Decomposition

```rust
// After simplification
let branches = simplifier.decompose_into_branches();

for (i, branch) in branches.iter().enumerate() {
    println!("Branch {}: {} edges, length = {:.2}",
             i, branch.num_edges(), branch.length());
}
```

## Algorithm Parameters

### Distortion Threshold

The `distortion_threshold` parameter controls when the algorithm halts:

- **Lower values** (e.g., 0.01): More conservative, preserves more detail
- **Higher values** (e.g., 10.0): More aggressive, creates simpler structures
- **Infinite**: Simplifies until no regular vertices remain

### Halting Conditions

The algorithm stops when either:

1. **Distortion threshold exceeded**: Next edge would cause too much distortion
2. **No regular vertices remain**: All remaining vertices are leaves or forks

## Mathematical Background

### Quadric Error Metrics

For an edge with direction `d`, the quadric matrix `Q` is constructed as:

```
Q = v_1·v_1ᵀ + v_2·v_2ᵀ
```

where `v_1` and `v_2` are orthonormal vectors perpendicular to `d`.

**Properties:**
- `Q` is symmetric (3×3)
- `Q` is positive semi-definite
- `Q · d = 0` (edge direction is in the null space)
- For point `p` on the edge: `pᵀ·Q·p` measures squared perpendicular distance

### Quadric Aggregation

When contracting edge `(a,b)` to vertex `c`:

```
Q_c = Q_a + Q_b
```

**Note**: This counts the contracted edge's contribution twice (once in `Q_a` and once in `Q_b`). The paper mentions this can be corrected by subtracting the edge's quadric once, but the current implementation keeps both contributions.

### Optimal Position (Theoretical)

The optimal position `c` should minimize:

```
f_c(x) = xᵀ·(Q_a + Q_b)·x
```

This requires solving:

```
(Q_a + Q_b)·c = 0
```

**Current Implementation**: Uses midpoint heuristic `c = (a + b)/2` when both vertices are regular. A full implementation would solve the linear system or compute the pseudoinverse.

## Implementation Status

### [DONE] Implemented Features

- [x] Vertex degree classification
- [x] Edge contraction operation
- [x] Quadric evaluation at arbitrary points
- [x] Distortion calculation
- [x] Priority queue management with lazy deletion
- [x] Main simplification loop
- [x] Branch decomposition
- [x] Comprehensive unit tests
- [x] Example program

### 🚧 Future Enhancements

- [ ] **True quadric minimization**: Solve `(Q_a + Q_b)·c = 0` using linear algebra library (e.g., nalgebra)
- [ ] **Pseudoinverse for degenerate cases**: Handle collinear or degenerate edge configurations
- [ ] **Edge quadric correction**: Subtract double-counted edge contribution from `Q_c`
- [ ] **Adaptive thresholding**: Automatically determine threshold based on geometry
- [ ] **Visualization**: Generate before/after plots of simplification
- [ ] **Progressive simplification**: Generate multiple levels of detail
- [ ] **Parallel processing**: Simplify independent subtrees in parallel

## Testing

The implementation includes comprehensive tests:

### `test_simple_linear_contraction`

Creates a linear chain `v1 -- v2 -- v3` where `v2` is regular (degree 2) and verifies it gets contracted.

### `test_fork_preservation`

Creates a Y-shaped structure with a fork vertex and verifies it is not contracted (forks are non-regular).

### `test_distortion_threshold`

Tests that the algorithm respects the distortion threshold and halts appropriately.

### `test_branch_decomposition`

Creates a structure with multiple branches and verifies branch decomposition works correctly.

Run tests with:

```bash
cargo test simplification
```

## Performance Considerations

### Time Complexity

- **Initialization**: O(E log E) for heap construction
- **Main loop**: O(E log E) worst case, typically better due to early termination
- **Per contraction**: O(log E) for heap operations, O(D) for incident edge updates where D is vertex degree

### Space Complexity

- O(V) for vertices
- O(E) for edges and heap
- O(V·D) for incident edge sets where D is average degree

### Optimization Opportunities

1. **Spatial indexing**: Use octree or k-d tree for nearest neighbor queries
2. **Incremental updates**: More efficient heap update strategies
3. **Parallelization**: Process independent subtrees concurrently
4. **Memory pooling**: Reuse vertex/edge allocations

## References

### Academic Papers

1. **Garland, M., & Heckbert, P. S. (1997)**. "Surface simplification using quadric error metrics." *SIGGRAPH '97*.
   - Original quadric error metric paper for mesh simplification

2. **Hoppe, H. (1996)**. "Progressive meshes." *SIGGRAPH '96*.
   - Edge collapse for progressive mesh representation

3. **Lindstrom, P., & Turk, G. (1998)**. "Fast and memory efficient polygonal simplification." *Visualization '98*.
   - Memory-efficient simplification techniques

### Related Work

- **Dendritic tree analysis**: Computational neuroscience methods for morphological characterization
- **Persistent homology**: Topological data analysis techniques adapted in `persistent_quadric.rs`
- **Tree simplification**: Algorithms for simplifying hierarchical structures

## License

MIT License - See LICENSE file for details

## Author

Yossi Bokor Bleile
December 2025

## Contributing

Contributions welcome! Areas for improvement:
- Implementing true quadric minimization
- Adding visualization capabilities
- Performance optimizations
- Additional test cases
- Documentation improvements
