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

docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json

@@ -0,0 +1,10 @@
+{
+  "label": "bilateral_filter.h",
+  "position": 4,
+  "link": {
+    "type": "generated-index",
+    "title": "Bilateral Filter",
+    "description": "Documentation for the Bilateral Filter in the Image WebAssembly (WASM) module in Img2Num.",
+    "slug": "/reference/wasm/modules/image/bilateral_filter"
+  }
+}

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

@@ -0,0 +1,31 @@
+---
+id: api
+title: Bilateral Filter — API & Reference
+sidebar_label: API / Usage
+sidebar_position: 5
+---
+
+# 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)
+
+## Parameters
+
+| Parameter | Type | Description |
+| :--- | :--- | :--- |
+| `image` | `uint8_t*` | Pointer to the RGBA image data (4 bytes per pixel). Modified in-place. |
+| `width` | `size_t` | Width of the image in pixels. |
+| `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). |
+
+:::info Implementation Details
+- **Namespace**: `bilateral` (C++)
+- **Export**: Exposed to WASM via `extern "C"` wrapper as `bilateral_filter`.
+:::

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

@@ -0,0 +1,78 @@
+---
+id: explained
+title: Implementation Explained
+sidebar_position: 6
+---
+
+# Bilateral Filter — Implementation Explained
+
+This section explains the inner workings of the **bilateral filter** implementation.
+
+## Overview
+
+The bilateral filter smoothes an image while **preserving edges**. It achieves this by weighting neighboring pixels based on two criteria:
+1.  **Spatial Distance**: Pixels closer to the center have higher weight.
+2.  **Range (Color) Difference**: Pixels with similar colors to the center have higher weight.
+
+This prevents the "blurring" from crossing strong edges, where the color difference is large.
+
+## 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).
+
+
+
+
+:::info
+In this implementation, both weighting terms are **Gaussian kernels**:
+
+$$
+w_{\text{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)
+$$
+
+where $ \sigma_s$ controls spatial smoothing and $\sigma_r$ controls edge sensitivity.
+:::
+## Implementation Details
+
+Our implementation uses a **naive sliding window** approach with **Look-Up Table (LUT) optimizations** to improve performance in WebAssembly.
+
+### 1. Precomputed Look-Up Tables
+
+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.
+- **Range Weights**: A 1D array mapping squared color distance ($0$ to $255^2 \times 3$) to a weight. This allows O(1) lookups for the "edge preservation" factor.
+
+```cpp title="Precomputing Range Weights"
+std::vector<double> range_lut(MAX_RGB_DIST_SQ + 1);
+for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) {
+    range_lut[i] = std::exp(-static_cast<double>(i) / two_sigma_range_sq);
+}
+```
+
+### 2. 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}}$.
+4.  **Accumulate**: `pixel_acc += neighbor_rgb * (spatial_w * range_w)`.
+5.  **Normalize**: Divide by probability sum.
+
+### Complexity
+
+- **Time Complexity**: $O(W \cdot H \cdot R^2)$, where $R$ is the kernel radius.
+- **Space Complexity**: $O(W \cdot H)$ for the output buffer.
+
+This complexity is why the filter can be slow for large radii ($\sigma_{spatial} > 5.0$), but we currently parameterize the radius to be small ($\sigma_{spatial} \leq 3.0$) so it is not a problem.

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

@@ -0,0 +1,73 @@
+---
+id: implementation
+title: Bilateral Filter — Implementation details
+sidebar_label: Implementation
+sidebar_position: 4
+---
+
+# Bilateral Filter — Implementation details
+
+This page maps the conceptual steps of the Bilateral Filter to the concrete functions and loops in the implementation.
+
+## 1. Parameters & Window Size
+
+The filter first calculates the kernel size based on the spatial standard deviation ($\sigma_{spatial}$).
+
+```cpp
+const int radius = static_cast<int>(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatial));
+const int kernel_width = 2 * radius + 1;
+```
+
+We primarily use $\sigma_{spatial} \approx 3.0$, which results in a kernel radius of 9 (width 19x19).
+
+## 2. Precomputing Weights (Optimization)
+
+To avoid computing `std::exp` millions of times per frame, we precalculate the weights.
+
+### Spatial Weights (constant per kernel)
+The distance pattern is the same for every pixel, so we calculate the distance-based weights once at the start.
+
+```cpp
+spatial_weights[(ky + radius) * kernel_width + (kx + radius)] =
+    std::exp(-dist2 / two_sigma_space_sq);
+```
+
+### Range Weights (LUT)
+We calculate the `similarity score` for every possible color difference ahead of time. We just measure the color difference and look up the precomputed weight in the table.
+
+```cpp
+for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) {
+    range_lut[i] = std::exp(-static_cast<double>(i) / two_sigma_range_sq);
+}
+```
+
+## 3. Sliding Window Loop
+
+The core processing happens in a nested loop over every pixel $(y, x)$. For each pixel, we:
+
+1.  **Iterate** over the window (from $-radius$ to $+radius$).
+2.  **Fetch** neighbor RGB values.
+3.  **Calculate** color difference (squared Euclidean distance).
+4.  **Lookup** spatial weight (from array) and range weight (from LUT).
+5.  **Accumulate** the weighted sum and the sum of weights.
+
+```cpp
+double w_space = spatial_weights[...];
+double w_range = range_lut[dist_sq];
+double w = w_space * w_range;
+
+r_acc += r * w;
+g_acc += g * w;
+b_acc += b * w;
+weight_acc += w;
+```
+
+## 4. Normalization
+
+Finally, we normalize the accumulated color values by the total weight to get the filtered pixel value:
+
+```cpp
+result[center_idx] = static_cast<uint8_t>(std::clamp(r_acc / weight_acc, 0.0, 255.0));
+```
+
+This ensures the pixel brightness remains consistent with the local area.

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

@@ -0,0 +1,31 @@
+---
+id: overview
+title: Bilateral Filter
+sidebar_label: Overview
+sidebar_position: 2
+---
+
+# Bilateral Filter
+
+This section introduces the **bilateral filter** used in the Img2Num project
+(see [`bilateral_filter.h`](https://github.com/Ryan-Millard/Img2Num/blob/main/src/wasm/modules/image/include/bilateral_filter.h)
+& [`bilateral_filter.cpp`](https://github.com/Ryan-Millard/Img2Num/blob/main/src/wasm/modules/image/src/bilateral_filter.cpp)).
+It focuses on how the algorithm is implemented, why each step is necessary,
+and where the corresponding code lives so you can jump straight into the implementation.
+
+## At a glance
+- **Algorithm:** Bilateral Filter (Non-linear, edge-preserving).
+- **Data type:** `uint8_t` (8-bit unsigned integer channels).
+- **Key steps:**
+  1. For each pixel, inspect neighbors in radius $R$.
+  2. Weight neighbors by **spatial distance** (Gaussian).
+  3. Weight neighbors by **intensity difference** (Gaussian).
+  4. Normalize and average.
+
+## Pages in this mini-guide
+
+* **Overview** (this page)
+* **Implementation details** — step-by-step mapping between theory and the actual C++ code.
+* **API & reference** — brief function signatures and purpose for quick lookup.
+
+Jump to implementation: [Implementation details](../implementation/)

docs/docs/reference/wasm/modules/image/overview.md

@@ -25,6 +25,7 @@ src/wasm/modules/image/
 │   ├── kmeans.h
 │   └── mergeSmallRegionsInPlace.h
 └── src
+    ├── bilateral_filter.cpp
     ├── fft_iterative.cpp
     ├── image_utils.cpp
     ├── kmeans.cpp
@@ -38,6 +39,7 @@ Each header corresponds to a major subsystem:
 - `Image.h` — Core image class.
     - Internally uses a **Pixel type**.
 - `PixelConverters` — Functions for converting between pixel formats.
+- `bilateral_filter` — Bilateral filter for edge-conserving image denoising.
 - `fft_iterative` — Fast Fourier Transform utilities.
     - Used by **Gaussian Blur** inside image_utils.h.
 - `kmeans` — K-means clustering used for quantization.

src/components/WasmImageProcessor.jsx

@@ -13,7 +13,7 @@ const WasmImageProcessor = () => {
   const inputId = useId();
   const inputRef = useRef(null);
 
-  const { gaussianBlur, blackThreshold, kmeans, mergeSmallRegionsInPlace } = useWasmWorker();
+  const { bilateralFilter, blackThreshold, kmeans, mergeSmallRegionsInPlace } = useWasmWorker();
 
   const [originalSrc, setOriginalSrc] = useState(null);
   const [fileData, setFileData] = useState(null);
@@ -81,12 +81,19 @@ const WasmImageProcessor = () => {
       const { width, height } = fileData;
 
       step(20);
-      const blurred = await gaussianBlur(fileData);
+      // NOTE: Gaussian blur destroys the sharp outlines first, preventing the Bilateral filter from detecting and preserving them
+      // const blurred = await gaussianBlur(fileData);
+
+      const imgBilateralFiltered = await bilateralFilter({
+        pixels: fileData.pixels,
+        width,
+        height,
+      });
 
       step(45);
       const thresholded = await blackThreshold({
         ...fileData,
-        pixels: blurred,
+        pixels: imgBilateralFiltered,
         num_colors: 8,
       });
 
@@ -139,7 +146,7 @@ const WasmImageProcessor = () => {
         step(0);
       }, 800);
     }
-  }, [fileData, gaussianBlur, blackThreshold, kmeans, mergeSmallRegionsInPlace, navigate, step]);
+  }, [fileData, bilateralFilter, blackThreshold, kmeans, mergeSmallRegionsInPlace, navigate, step]);
 
   /* Memo'd UI fragments */
   const EmptyState = useMemo(

src/hooks/useWasmWorker.js

@@ -35,6 +35,9 @@ export function useWasmWorker() {
   const gaussianBlur = async ({ pixels, width, height, sigma_pixels = width * 0.005 }) => {
     return (await call('gaussian_blur_fft', { pixels, width, height, sigma_pixels }, ['pixels'])).output.pixels;
   };
+  const bilateralFilter = async ({ pixels, width, height, sigma_spatial = 3.0, sigma_range = 50.0 }) => {
+    return (await call('bilateral_filter', { pixels, width, height, sigma_spatial, sigma_range }, ['pixels'])).output.pixels;
+  };
   const blackThreshold = async ({ pixels, width, height, num_colors }) => {
     return (await call('black_threshold_image', { pixels, width, height, num_colors }, ['pixels'])).output.pixels;
   };
@@ -46,5 +49,5 @@ export function useWasmWorker() {
       .output.pixels;
   };
 
-  return { call, gaussianBlur, blackThreshold, kmeans, mergeSmallRegionsInPlace };
+  return { call, gaussianBlur, bilateralFilter, blackThreshold, kmeans, mergeSmallRegionsInPlace };
 }

src/hooks/useWasmWorker.test.js

@@ -59,12 +59,14 @@ describe('useWasmWorker', () => {
 
     expect(result.current).toHaveProperty('call');
     expect(result.current).toHaveProperty('gaussianBlur');
+    expect(result.current).toHaveProperty('bilateralFilter');
     expect(result.current).toHaveProperty('blackThreshold');
     expect(result.current).toHaveProperty('kmeans');
     expect(result.current).toHaveProperty('mergeSmallRegionsInPlace');
 
     expect(typeof result.current.call).toBe('function');
     expect(typeof result.current.gaussianBlur).toBe('function');
+    expect(typeof result.current.bilateralFilter).toBe('function');
     expect(typeof result.current.blackThreshold).toBe('function');
     expect(typeof result.current.kmeans).toBe('function');
     expect(typeof result.current.mergeSmallRegionsInPlace).toBe('function');
@@ -202,6 +204,59 @@ describe('useWasmWorker', () => {
     });
   });
 
+  describe('bilateralFilter', () => {
+    it('should call worker with bilateral_filter function', async () => {
+      const { result } = renderHook(() => useWasmWorker());
+
+      const pixels = new Uint8ClampedArray([255, 0, 0, 255]);
+      const width = 1;
+      const height = 1;
+
+      act(() => {
+        result.current.bilateralFilter({ pixels, width, height });
+      });
+
+      expect(mockWorkerInstance.postMessage).toHaveBeenCalledWith(
+        expect.objectContaining({
+          funcName: 'bilateral_filter',
+          args: expect.objectContaining({
+            pixels,
+            width,
+            height,
+            sigma_spatial: 3.0,
+            sigma_range: 50.0,
+          }),
+          bufferKeys: ['pixels'],
+        })
+      );
+    });
+
+    it('should use custom sigma_spatial and sigma_range when provided', async () => {
+      const { result } = renderHook(() => useWasmWorker());
+
+      const pixels = new Uint8ClampedArray([255, 0, 0, 255]);
+
+      act(() => {
+        result.current.bilateralFilter({
+          pixels,
+          width: 100,
+          height: 100,
+          sigma_spatial: 5.0,
+          sigma_range: 25.0,
+        });
+      });
+
+      expect(mockWorkerInstance.postMessage).toHaveBeenCalledWith(
+        expect.objectContaining({
+          args: expect.objectContaining({
+            sigma_spatial: 5.0,
+            sigma_range: 25.0,
+          }),
+        })
+      );
+    });
+  });
+
   describe('blackThreshold', () => {
     it('should call worker with black_threshold_image function', async () => {
       const { result } = renderHook(() => useWasmWorker());

src/wasm/modules/image/include/bilateral_filter.h

@@ -0,0 +1,21 @@
+#ifndef BILATERAL_FILTER_H
+#define BILATERAL_FILTER_H
+
+#include <cstddef> // for size_t
+#include <cstdint> // for uint8_t
+
+namespace bilateral {
+
+// Apply bilateral filter to an image.
+// The filter modifies the image buffer in-place.
+// Parameters:
+//  - image: Pointer to RGBA pixel buffer
+//  - width, height: Image dimensions (px)
+//  - sigma_spatial: Gaussian standard deviation for spatial proximity (spatial decay)
+//  - sigma_range: Gaussian standard deviation for intensity difference (radiometric decay)
+void bilateral_filter(uint8_t *image, size_t width, size_t height,
+                      double sigma_spatial, double sigma_range);
+
+} // namespace bilateral
+
+#endif // BILATERAL_FILTER_H