WIP: feat(bilateral filter): implement bilateral filter for denoising before K-Means in RGB & CIELAB

docs/docs/reference/wasm/modules/image/bilateral_filter/api.md

@@ -2,20 +2,21 @@
 id: api
 title: Bilateral Filter — API & Reference
 sidebar_label: API / Usage
-sidebar_position: 5
+sidebar_position: 4
 ---
 
 # Bilateral Filter — API & Reference
 
 Quick reference for the function implemented in the header.
 
-```cpp title="Applies a bilateral filter to an RGBA image (modified in-place)."
-void bilateral_filter(uint8_t *image,
-                      size_t width, size_t height,
-                      double sigma_spatial,
-                      double sigma_range)
+```cpp title="Applies a bilateral filter to an RGBA uint8_t* image (modified in-place)."
+void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, double sigma_range, uint8_t color_space)
 ```
 
+:::important Alpha Channel Preservation
+The alpha channel, `image[i + 3]`, is left untouched - it is not part of the bilateral filter implementation.
+:::
+
 ## Parameters
 
 | Parameter | Type | Description |
@@ -25,6 +26,7 @@ void bilateral_filter(uint8_t *image,
 | `height` | `size_t` | Height of the image in pixels. |
 | `sigma_spatial` | `double` | Spatial standard deviation ($\sigma_s$). Controls how far pixels influence each other spatially. |
 | `sigma_range` | `double` | Range standard deviation ($\sigma_r$). Controls how much color definition is preserved (edge preservation). |
+| `color_space` | `uint8_t` | Toggle color space to use for range distance (0 - CIELAB, 1 - RGB). CIELAB produces perceptually better results but requires more computation. |
 
 :::info Implementation Details
 - **Namespace**: `bilateral` (C++)

docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md

@@ -0,0 +1,219 @@
+---
+id: color-spaces
+title: Color Space Selection — RGB vs CIELAB
+sidebar_label: Color Space Selection
+sidebar_position: 6
+---
+
+# Color Space Selection — RGB vs CIELAB
+
+The bilateral filter in Img2Num supports two color spaces for computing range (color) distances: **RGB** and **CIELAB**. This guide explains the differences, trade-offs, and when to use each.
+
+## Quick Comparison
+
+| Aspect | RGB | CIELAB |
+|:---|:---|:---|
+| **Perceptual accuracy** | Lower — equal Euclidean distances don't correspond to equal perceived color differences | Higher — designed to be perceptually uniform |
+| **Performance** | Faster — uses precomputed LUT | Slower — requires conversion and on-the-fly computation |
+| **Edge preservation** | Good for most images | Better for images with subtle color transitions |
+| **Best for** | General purpose, real-time applications | High-quality processing, perceptual accuracy |
+
+## When to Use Each Color Space
+
+### Use RGB when:
+- **Performance is critical** — RGB processing is significantly faster due to LUT optimization
+- **Working with high-contrast images** — where edge preservation is less sensitive to color space choice
+- **Real-time processing** — where milliseconds matter
+- **Sigma_range values are well-tuned** — and visual results are satisfactory
+
+### Use CIELAB when:
+- **Perceptual uniformity matters** — you want visually equal smoothing across different hues
+- **Working with skin tones or subtle gradients** — where human perception is sensitive
+- **Quality over speed** — when processing time is less critical than output quality
+- **Processing medical or scientific imagery** — where perceptual accuracy is important
+
+## Mathematical Differences
+
+### Distance Metrics
+
+Both color spaces compute the Euclidean distance between color vectors, but the ranges differ significantly.
+
+#### RGB Color Space
+RGB channels are bounded `[0, 255]` per channel:
+
+$$
+\text{distance}_{\text{RGB}} = \sqrt{\Delta R^2 + \Delta G^2 + \Delta B^2}
+$$
+
+Maximum possible distance:
+$$
+\text{max}_{\text{RGB}} = \sqrt{255^2 + 255^2 + 255^2} \approx 441.67
+$$
+
+#### CIELAB Color Space
+CIELAB channels have different ranges:
+- **L\***: `[0, 100]` (lightness)
+- **a\***: approximately `[-128, 127]` (green-red)
+- **b\***: approximately `[-128, 127]` (blue-yellow)
+
+$$
+\text{distance}_{\text{LAB}} = \sqrt{\Delta L^2 + \Delta a^2 + \Delta b^2}
+$$
+
+Maximum theoretical distance:
+$$
+\text{max}_{\text{LAB}} = \sqrt{100^2 + 255^2 + 255^2} \approx 373.56
+$$
+
+:::important Key Insight
+In practice, most real-world pixel differences are **much smaller** than the maximum possible distance. CIELAB distances for neighboring pixels are typically smaller than RGB distances due to:
+1. **Numerical compression** from the RGB→LAB conversion
+2. **Perceptual scaling** — LAB is designed to reflect human vision, which perceives smaller differences
+:::
+
+## Sigma_range Behavior Differences
+
+The `sigma_range` parameter controls edge preservation by weighting color similarity. However, the same `sigma_range` value produces **different visual results** in RGB vs CIELAB.
+
+### The Range Weight Formula
+
+The bilateral filter computes range weights using a Gaussian:
+
+$$
+w_{\text{range}} = \exp\left(-\frac{\text{distance}^2}{2\sigma_{\text{range}}^2}\right)
+$$
+
+- When distance is **small**, weight is **high** (≈1) → strong contribution
+- When distance is **large**, weight is **low** (≈0) → weak contribution
+
+### Why the Same Sigma Produces Different Results
+
+import RgbVsLabRangeKernel from '@site/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel';
+
+<RgbVsLabRangeKernel />
+
+**With `sigma_range = 50`**:
+- **RGB**: Typical neighboring pixel distances are small relative to 50, so many neighbors contribute significantly → **moderate blur**
+- **CIELAB**: Typical neighboring pixel distances are even smaller, so almost all neighbors contribute strongly → **stronger blur**
+
+### Sigma_range Scaling for Visual Consistency
+
+To achieve **visually similar** blur between RGB and CIELAB, you can scale `sigma_range`:
+
+```javascript
+// Example: Scaling RGB sigma_range to match CIELAB visual output
+const sigma_range_base = 50.0; // Target CIELAB sigma_range
+
+let sigma_range_actual;
+if (color_space === COLOR_SPACE_RGB) {
+    // Scale RGB sigma_range to match CIELAB perceptually
+    sigma_range_actual = sigma_range_base * 4.18;
+} else {
+    sigma_range_actual = sigma_range_base;
+}
+```
+
+:::important Scaling Factor
+The scaling factor of **~4.18** is empirically derived and works well for natural images. However:
+- It's **not universal** — depends on image statistics
+- It's **not mandatory** — the different behaviors are valid features of each color space
+- **Advanced users** may want different sigma_range values for each space
+:::
+
+### Visual Example
+
+Using the same `sigma_range = 50`:
+
+| Color Space | Visual Result |
+|:---|:---|
+| **CIELAB** | Stronger smoothing, better edge preservation in perceptually uniform manner |
+| **RGB** | Moderate smoothing, adequate edge preservation for most use cases |
+| **RGB (scaled)** | Similar smoothing to CIELAB when `sigma_range ≈ 209` |
+
+## Performance Considerations
+
+### RGB Performance
+- **Precomputed LUT**: All 195,075 possible squared distances are precomputed
+- **O(1) lookup**: Range weight retrieval is extremely fast
+- **Memory**: ~1.5 MB for LUT (acceptable for most applications)
+
+### CIELAB Performance
+- **Full image conversion**: RGB→LAB conversion for entire image upfront
+- **On-the-fly computation**: Range weights computed using `exp()` for each neighbor
+- **Slower but optimized**: Conversion is done once; only distance calculation repeated
+
+**Performance Impact**: CIELAB is typically **2-4× slower** than RGB, depending on image size and kernel radius.
+
+:::tip Optimization Note
+Future optimizations may include:
+- Taylor/Horner polynomial approximations for `exp(-x²)`
+- SIMD vectorization for distance calculations
+- Adaptive LUT for CIELAB (with quantization)
+:::
+
+## Implementation Details
+
+### RGB Range Weights (LUT)
+```cpp
+// Precompute all possible RGB distances
+std::vector<double> range_lut(MAX_RGB_DIST_SQ + 1);
+for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) {
+    range_lut[i] = gaussian(std::sqrt(i), sigma_range);
+}
+
+// Later, during filtering:
+const int dr = r_neighbor - r_center;
+const int dg = g_neighbor - g_center;
+const int db = b_neighbor - b_center;
+const int dist_sq = dr*dr + dg*dg + db*db;
+double w_range = range_lut[dist_sq]; // O(1) lookup
+```
+
+### CIELAB Range Weights (On-the-fly)
+```cpp
+// Precompute full-image RGB → LAB conversion
+std::vector<double> cie_image(width * height * 4);
+for (each pixel) {
+    rgb_to_lab(r, g, b, L, A, B);
+    cie_image[idx] = L; cie_image[idx+1] = A; cie_image[idx+2] = B;
+}
+
+// Later, during filtering:
+double dL = L_neighbor - L_center;
+double dA = A_neighbor - A_center;
+double dB = B_neighbor - B_center;
+double dist = std::sqrt(dL*dL + dA*dA + dB*dB);
+double w_range = gaussian(dist, sigma_range); // Computed on-the-fly
+```
+
+## Recommendations
+
+### Default Choice
+For most applications, **RGB** is the recommended default:
+- ✅ Faster processing
+- ✅ Good results for general images
+- ✅ Predictable behavior
+
+### When to Switch to CIELAB
+Consider CIELAB when you observe:
+- Inconsistent smoothing across different hues
+- Need for perceptually uniform processing
+- Working with images where color accuracy is critical
+- Willing to accept 2-4× performance cost
+
+### Parameter Tuning
+
+**Starting values**:
+- `sigma_spatial = 3.0` (both color spaces)
+- `sigma_range = 50.0` (CIELAB) or `sigma_range = 200.0` (RGB for similar visual effect)
+
+**Adjustment guidelines**:
+- Increase `sigma_range` → more blur, less edge preservation
+- Decrease `sigma_range` → sharper edges, less smoothing
+- Test with your specific images — optimal values vary by content
+
+## See Also
+
+- [Implementation Details](./implementation.md#range-weights) — Deep dive into LUT vs on-the-fly computation
+- [API Reference](./api.md) — `color_space` parameter documentation
+- [Keywords](./keywords.md) — Understanding range and spatial components

docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md

@@ -1,7 +1,7 @@
 ---
 id: explained
 title: Implementation Explained
-sidebar_position: 6
+sidebar_position: 5
 ---
 
 # Bilateral Filter — Implementation Explained
@@ -19,24 +19,37 @@ This prevents the "blurring" from crossing strong edges, where the color differe
 ## How It Works
 
 For each pixel in the image, we look at a local window (kernel) around it. The new pixel value is a weighted average of its neighbors:
-
 $$
 I_{new}(x) = \frac{1}{W_p} \sum_{x_i \in \Omega} I(x_i) \cdot w_{spatial}(\|x_i - x\|) \cdot w_{range}(|I(x_i) - I(x)|)
 $$
 
-Where:
-- $w_{spatial}$ is a Gaussian function of the distance.
-- $w_{range}$ is a Gaussian function of the intensity difference.
-- $W_p$ is the normalization factor (sum of all weights).
-
-
-
+Where each component means:
+
+- $x$: The coordinates of the **center pixel** being filtered.
+- $\Omega$: The set of **neighboring pixels** in the local kernel around $x$ (from `-radius` to `+radius`).
+- $I(x_i)$: The **color or intensity** of a neighbor pixel $x_i$.
+- $C(x_i)$: The **color vector** of pixel $x_i$.
+  - RGB: `[R, G, B]`
+  - CIELAB: `[L*, a*, b*]`
+- $w_{spatial}(|x_i - x|)$: A **Gaussian weight** based on the **spatial distance** between the neighbor and the center.
+  - Pixels closer to the center have **larger weights**.
+  - Formula: $\exp\Big(-\frac{\text{distance}^2}{2\sigma_s^2}\Big)$
+- $w_\text{range}(|C(x_i) - C(x)|)$: A **Gaussian weight** based on the **color difference** between neighbor and center.
+  - Pixels with **similar colors** have higher weights, preserving edges.
+  - Formula: $\exp\Big(-\frac{|C(x_i) - C(x)|^2}{2\sigma_r^2}\Big)$
+  - RGB: Precomputed via **LUT**
+  - CIELAB: Computed **on the fly**
+- $W_p = \sum_{x_i \in \Omega} w_{spatial} \cdot w_\text{range}$: **Normalization factor** to ensure the weighted average sums to a valid color.
+- **Result $I_\text{new}(x)$**: The **filtered color** of the center pixel after combining spatial and color-based weighting.
+- $|I(x_i) - I(x)|$: The Euclidean norm.
+  - **RGB**: $|I(x_i) - I(x)| = \sqrt{ \Delta R² + \Delta G² + \Delta B² }$
+  - **CIELAB**: $|I(x_i) - I(x)| = \sqrt{ \Delta L² + \Delta a² + \Delta b² }$
 
 :::info
 In this implementation, both weighting terms are **Gaussian kernels**:
 
 $$
-w_{\text{spatial}}(d) = \exp!\left(-\frac{d^2}{2\sigma_s^2}\right),
+w_{spatial}(d) = \exp!\left(-\frac{d^2}{2\sigma_s^2}\right),
 \quad
 w_{\text{range}}(d) = \exp!\left(-\frac{d^2}{2\sigma_r^2}\right)
 $$
@@ -45,9 +58,9 @@ where $ \sigma_s$ controls spatial smoothing and $\sigma_r$ controls edge sensit
 :::
 ## Implementation Details
 
-Our implementation uses a **naive sliding window** approach with **Look-Up Table (LUT) optimizations** to improve performance in WebAssembly.
+Our implementation uses a **naive sliding window** approach with **Look-Up Table (LUT) optimizations / "On the Fly" computations** to improve performance.
 
-### 1. Precomputed Look-Up Tables
+### 1. Precomputed Look-Up Tables (RGB color space)
 
 Calculating `std::exp()` inside the inner loop is expensive. We precompute the two Gaussian functions:
 - **Spatial Weights**: A 2D grid of weights based on the kernel radius. Since the spatial distance between a neighbor and the center never changes, this is calculated once per filter application.
@@ -60,13 +73,35 @@ for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) {
 }
 ```
 
-### 2. The Loop
+### 2. On-the-fly Range weights (CIE-LAB color space)
+
+When deriving range weights in the CIELAB color space, the LUT approach does not work
+(see the [Range Weights section in the implementation docs](../implementation/#range-weights) to understand why).
+Instead range weights are computed on the fly using the `gaussian` function.
+
+Since the RGB to CIELAB conversion is expensive, redundant computations are minimized by initially converting the full RGB image to CIELAB image.
+
+In the convolution step LAB distance is computed by reading those values from the CIELAB image buffer, and the gaussian is then evaluated.
+```
+dL = cie_image[neighbor_idx] - L0;
+dA = cie_image[neighbor_idx + 1] - A0;
+dB = cie_image[neighbor_idx + 2] - B0;
+
+dist = std::sqrt(dL * dL + dA * dA + dB * dB);
+w_range = gaussian(dist, sigma_range);
+```
+
+:::note
+`gaussian` itself is expensive to run. Future optimizations will include polynomial approximations of `exp(-x^2)` via Taylor expansion or Horner's method.
+:::
+
+### 3. The Loop
 
 We iterate over every pixel `(y, x)` and then over every neighbor `(ky, kx)` within the kernel radius:
 
 1.  **Load Neighbor**: Get RGB values of the neighbor.
 2.  **Spatial Weight**: Look up precomputed $G_{\sigma_{spatial}}$.
-3.  **Range Weight**: Calculate squared color distance $\|C_p - C_q\|^2$ and look up precomputed $G_{\sigma_{range}}$.
+3.  **Range Weight**: Calculate squared color distance $\|C_p - C_q\|^2$ and look up precomputed $G_{\sigma_{range}}$ if using RGB, or compute on the fly if using CIELAB.
 4.  **Accumulate**: `pixel_acc += neighbor_rgb * (spatial_w * range_w)`.
 5.  **Normalize**: Divide by probability sum.