diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json new file mode 100644 index 000000000..2ad8ea622 --- /dev/null +++ b/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" + } +} diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md new file mode 100644 index 000000000..37485ec28 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -0,0 +1,53 @@ +--- +id: api +title: Bilateral Filter — API & Reference +sidebar_label: API / Usage +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 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 | +| :-------------- | :--------- | :--------------------------------------------------------------------------------------------------------------------------------------------- | +| `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). | +| `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++) +- **Export**: Exposed to WASM via `extern "C"` wrapper as `bilateral_filter`. + ::: + +:::tip Color Space Discrepancies +As noted on the +[Color Space Selection page](../color-spaces/#why-the-scaling-factor-exists-and-why-418-works), +the bilateral filter will produce **different results** depending on the selected +`color_space`, even with identical parameters. + +To achieve **visually equivalent filtering behavior** between CIELAB and RGB, +treat CIELAB as the reference space and scale `sigma_range` for RGB: + +$$ +\sigma_{\text{range, RGB}} \approx 4.18 \times \sigma_{\text{range, CIELAB}} +$$ + +> The factor **4.18** is empirically derived for natural images and equalizes bilateral +> range weights across color spaces. Any value in the range **[4.1, 4.3]** will typically +> produce comparable results. This is a recommended default, not a universal constant. +> ::: diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md new file mode 100644 index 000000000..ec4ded405 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md @@ -0,0 +1,332 @@ +--- +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'; + + + +**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` | + +## Why the Scaling Factor Exists (and Why ~4.18 Works) + +RGB and CIELAB do **not** measure color differences on the same numeric scale. As a result, identical `sigma_range` values will generally not produce equivalent range weights or visual results. + +### What “equivalent behavior” means mathematically + +The bilateral filter’s range weight is defined as: + +$$ +w_{\text{range}} = \exp!\left(-\frac{d^2}{2\sigma_{\text{range}}^2}\right) +$$ + +For RGB and CIELAB to behave equivalently, they must produce **the same range weight** for corresponding color differences: + +$$ +\exp\left(-\frac{d_{\text{RGB}}^2}{2\sigma_{\text{RGB}}^2}\right) +\approx +\exp\left(-\frac{d_{\text{LAB}}^2}{2\sigma_{\text{LAB}}^2}\right) +$$ + +Taking the logarithm and simplifying yields: + +$$ +\frac{d_{\text{RGB}}}{\sigma_{\text{RGB}}} +\approx +\frac{d_{\text{LAB}}}{\sigma_{\text{LAB}}} +$$ + +This implies the required relationship: + +$$ +\sigma_{RGB} \approx +\frac{d_{\text{RGB}}}{d_{\text{LAB}}} +\sigma_{LAB} +$$ + +So the scaling factor is **not arbitrary** — it is the **ratio of typical RGB distances to LAB distances** for the same pixel differences. + +### Where the value ~4.18 comes from + +For natural images (photographic content, sRGB, D65): + +1. Sample many _local_ pixel pairs (neighbors). +2. Measure: + - $$d_{\text{RGB}} = \sqrt{\Delta R^2 + \Delta G^2 + \Delta B^2}$$ + - $$d_{\text{LAB}} = \sqrt{\Delta L^2 + \Delta a^2 + \Delta b^2}$$ + +3. Compute the ratio $\frac{d_{RGB}}{d_{LAB}}$. +4. Aggregate (mean or median). + +Across a wide range of natural images, this ratio consistently clusters around: + +$$ +\boxed{4.1 \text{ to } 4.3} +$$ + +The value **4.18** lies near the center of this empirical range and provides a strong default for matching bilateral range behavior between RGB and CIELAB. + +### Why this ratio is stable (but not universal) + +The factor remains stable for natural images because: + +- **LAB compresses perceptual differences** + Equal perceived color changes produce smaller numeric deltas than in RGB. +- **RGB channels are highly correlated** + Euclidean RGB distance accumulates redundant energy across channels. +- **Bilateral filters operate locally** + In the small-delta regime, the RGB→LAB transform is locally quasi-linear. + +However, the factor may vary if: + +- Images are synthetic or heavily quantized +- A different RGB color space or white point is used +- LAB components are re-weighted or normalized differently + +### Practical guidance + +- **Recommended default** + + For visually comparable smoothing on natural images, use: + + $$ + \sigma_{range_{RGB}} \approx 4.18 \times \sigma_{range_{CIELAB}} + $$ + +- **Advanced usage** + For strict equivalence, compute the ratio + $$ + k = \frac{\mathbb{E}[d_{RGB}]}{\mathbb{E}[d_{LAB}]} + $$ + on your image set and scale `sigma_range` accordingly. + +## 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 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 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 diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md new file mode 100644 index 000000000..642d456ef --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -0,0 +1,118 @@ +--- +id: explained +title: Implementation Explained +sidebar_position: 5 +--- + +# 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 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_{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 / "On the Fly" computations** to improve performance. + +### 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. +- **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 range_lut(MAX_RGB_DIST_SQ + 1); +for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) { + range_lut[i] = std::exp(-static_cast(i) / two_sigma_range_sq); +} +``` + +### 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. + +```cpp +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}}$ 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. + +### 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. diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md new file mode 100644 index 000000000..674ceeb9f --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md @@ -0,0 +1,219 @@ +--- +id: implementation +title: Bilateral Filter — Implementation details +sidebar_label: Implementation +sidebar_position: 2 +--- + +# Bilateral Filter — Implementation details + +This page maps the conceptual steps of the Bilateral Filter to the concrete implementation. + +:::important Opacity +The opacity of individual pixels in this implementation are ignored. +::: + +## 1. Parameters & Kernel Size + +The filter first calculates the kernel size based on the spatial standard deviation ($\sigma_{spatial}$). + +```cpp +const int radius = static_cast(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 + the center pixel (width 19x19). + +
+ + {/* background grid */} + {[...Array(19 * 19)].map((_, i) => { + const x = i % 19 + const y = Math.floor(i / 19) + return ( + + ) + })} + + {/* radius labels */} + {[...Array(9)].map((_, r) => ( + + {r + 1} + + ))} + + {/* spatial support: r=9 + 0.5 to account for center */} + + + {/* center pixel */} + + + +
+ +:::important Kernel Dimensions +The filter kernel is always square; width = height = 2 \* radius + 1. +::: + +## 2. Computing Weights + +To avoid computing `std::exp` millions of times per frame, we precompute the spatial and range weights in the RGB color space +and only the spatial weights in the CIELAB color space. + +:::note CIELAB Color Space Range Weight Computation +For the CIELAB color space, the range weights are computed "on the fly" to reduce processing time. +::: + +:::info +To calculate the weights, we use `gaussian`, a simple Gaussian function that performs the calculation: +$\exp\left(-\frac{x^2}{2\sigma_{spatial}^2}\right)$ + +```cpp +double gaussian(double x, double sigma) { + return std::exp(-(x * x) / (2.0 * sigma * sigma)); +} +``` + +::: + +### 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_diameter + (kx + radius)] = gaussian(dist, sigma_spatial); +``` + +:::note Similarity Between Color Spaces +Since the color spaces (CIELAB and RGB) represent the **range component** of the image and not the spatial component, +the logic here is the same regardless of the color space (x & y are the spatial component present in every image). +::: + +### Range Weights + +#### Precomputed RGB Range Weights (LUT) + +In the RGB color space, $I(x_i) - I(x)$ is simple (each channel is $[0,255]$), so +we calculate the `similarity score` for every possible color difference ahead of time and store the results in an LUT. +This saves computation time by allowing us to measure the color difference and look up the precomputed weight in the table during the main body of the filter. + +```cpp +range_lut[i] = gaussian(static_cast(std::sqrt(i)), sigma_range); +``` + +:::important LUT Usage +We precompute a lookup table for all possible differences (0–255 for each channel, or squared Euclidean differences 0–195075) in +the RGB color space because it is small enough to store. + +See the corresponding information block for CIELAB to see why this differs between the color spaces. +::: + +#### "On the Fly" CIELAB Range Weights + +In the CIELAB color space, the pixels are not bounded $[0,255]$ per channel like RGB. +Thus, we calculate the `similarity score` for every possible color difference as we need them during the main body of the bilateral filter ("on the fly"). + +```cpp +w_range = gaussian(dist, sigma_range); +``` + +:::important "On the fly" vs. LUT +In **CIELAB**, the pixels are not bounded 0–255 per channel in the same way: + +- L: $[0,100]$ +- a: roughly $[−128,127]$ +- b: roughly $[−128,127]$ + +But more importantly: + +1. **Continuous values:** After conversion from RGB, the values are floating-point. + The differences ($|L^*a^*b^* - L^*a^*b^*|^2$) are continuous, not integers. + So the LUT would need to store **all possible floating-point differences**, which is essentially impossible. +2. **Large dynamic range:** The squared Euclidean distance in Lab can be **much larger than in 8-bit RGB**, especially when using floating-point precision. + Precomputing a LUT with sufficient precision would be huge. +3. **Precision matters:** Small errors in range weights in Lab are more noticeable because the filter is very sensitive to perceptual color distances. + A coarse LUT could lead to visible artifacts. + ::: + +
+ +

RGB LUT vs CIELAB On-the-Fly Weights

+
+ +Bilateral filtering computes a **range weight** for each pixel in the kernel based on the color difference between the center pixel and its neighbor. + +import RgbVsLabRangeKernel from '@site/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel'; + + + +##### RGB LUT (Red curve and shaded area) + +In the RGB color space, the maximum possible color difference is limited (0–255 per channel). +This allows us to **precompute all possible weights** in a **Lookup Table (LUT)**. +During filtering, we simply **look up the weight** instead of recomputing it with +$\exp\left(-\frac{x^2}{2\sigma_{spatial}^2}\right)$ +for each neighbor. +The discrete nature of the LUT is represented by the shaded area and the curve shows how weight decays with increasing ΔRGB. + +##### CIELAB (Blue curve and shaded area) + +In the CIELAB color space, the number of possible differences is much larger and continuous. +Precomputing a LUT would require enormous memory, so weights are **computed on-the-fly**. +The curve represents the weight for a given color difference ΔLAB, and the shaded area illustrates the range of influence. + +:::note +This visual shows why RGB weights can be precomputed while CIELAB weights must be computed dynamically. +The **height of the curve/area corresponds to the weight** given by the Gaussian function: higher means more influence in the filtered pixel. +::: + +
+ +## 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` (from left to right inside the effective circular range and domain of the kernel). +2. **Fetch** neighbor RGB values. +3. **Calculate** color difference using squared Euclidean distance. +4. **Lookup weights:** + - **Spatial weights:** Precomputed at the start of the bilateral filter. + - **Range weights:** + - _RGB_: From LUT (precomputed at the start of the bilateral filter). + - _CIELAB_: Calculate "on the fly". +5. **Accumulate** the weighted sum and the sum of weights. + +## 4. Normalization + +Finally, we normalize the accumulated color values by the total weight (clamped to valid RGB values: $[0,255]$) to get the filtered pixel value: + +```cpp +result[center_idx] = static_cast(std::clamp(r_acc / weight_acc, 0.0, 255.0)); +``` + +This ensures the pixel brightness remains consistent with the local area. diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md new file mode 100644 index 000000000..a57010075 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md @@ -0,0 +1,27 @@ +--- +id: keywords +title: Keywords +--- + +- **spatial component**: The part of an image related to the **pixel positions** (x and y coordinates). + In the bilateral filter, this determines how **neighboring pixels are weighted based on distance** from the center pixel. + +- **range component**: The part of an image related to **pixel values**, such as color, brightness, or intensity. + In the bilateral filter, this determines how **neighboring pixels are weighted based on similarity in color or intensity**. + +- **kernel / window / bounding box**: A local subset of pixels around the center pixel. + - This is the region over which the bilateral filter computes weighted averages. + - Gaussian functions define the **weights for each pixel in the kernel**, considering both spatial and range components. + +- **standard deviation ($\sigma$)**: A measure of how spread out values are from their mean. + - In the bilateral filter, $\sigma$ controls the **width of the Gaussian weighting**. + - **$\sigma_{spatial}$**: Controls the influence of **distance** — larger values allow more distant pixels to contribute. + - **$\sigma_{range}$**: Controls the influence of **color/intensity differences** — larger values make edges less sharp. + +- **LUT (Look-Up Table)**: A precomputed array mapping input values to output values to **avoid repeated computation**. + - In the bilateral filter, RGB range weights are often stored in a LUT for **fast access**, while CIELAB weights are computed on the fly. + +- **weighted average**: A sum of values multiplied by their corresponding weights, then normalized by the total weight. + - The bilateral filter uses this to combine neighbor pixels into the **filtered center pixel value**. + +- **edge preservation**: The ability of the filter to **smooth flat regions while maintaining sharp transitions** at boundaries between different colors or intensities. diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md new file mode 100644 index 000000000..0bb808967 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md @@ -0,0 +1,42 @@ +--- +id: overview +title: Bilateral Filter +sidebar_label: Overview +sidebar_position: 1 +--- + +# 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. + +:::important Why Img2Num Uses Bilateral Filters +In Img2Num, the bilateral filter is used to **reduce noise while preserving edges**, which is critical for accurate image +segmentation (via methods like K-Means clustering), contour extraction and vectorization. + +Similarly to Gaussian blurs, it acts as a _low-pass filter_ that reduces noise. +Conversely, it is _less aggressive than Gaussian blurs, since it takes spatial position (x & y coordinates) into account_ - +allowing it to preserve sharp edges. +::: + +## At a glance + +- **Algorithm:** Bilateral Filter (Non-linear, edge-preserving). +- **Input/Output image data types:** `uint8_t` (8-bit unsigned integer channels). +- **Color spaces:** RGB & CIELAB can be chosen (see `color_space` in the [**API / Usage** section](../api/)). +- **Key steps:** + 1. For each pixel, inspect neighbors in radius $R$. + 2. Weight neighbors by **spatial distance** (Gaussian). + 3. Weight neighbors by **intensity/color difference** (Gaussian). + 4. Normalize and average. + +## Keywords + +The [keywords section](../keywords/) will help you in case the terminology confuses you. + +:::tip Gaussian functions +Bilateral filters rely on Gaussian weighting, so understanding Gaussian functions will help when reading the implementation. +::: diff --git a/docs/docs/reference/wasm/modules/image/cielab/_category_.json b/docs/docs/reference/wasm/modules/image/cielab/_category_.json new file mode 100644 index 000000000..13cce7a36 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/cielab/_category_.json @@ -0,0 +1,10 @@ +{ + "label": "cielab.h", + "position": 5, + "link": { + "type": "generated-index", + "title": "CIELAB Utilities", + "description": "Documentation for the CIELAB Utilities in the Image WebAssembly (WASM) module in Img2Num.", + "slug": "/reference/wasm/modules/image/cielab" + } +} diff --git a/docs/docs/reference/wasm/modules/image/cielab/api.md b/docs/docs/reference/wasm/modules/image/cielab/api.md new file mode 100644 index 000000000..b8b4af65d --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/cielab/api.md @@ -0,0 +1,218 @@ +--- +id: cielab-api +title: CIELAB Color Space API +sidebar_label: API Reference +sidebar_position: 2 +--- + +# CIELAB Color Space Conversion API + +## Overview + +The CIELAB module provides bidirectional conversion between 8-bit sRGB and CIELAB (CIE L\*a\*b\*) color spaces. CIELAB is a perceptually uniform color space designed to approximate human vision, where equal Euclidean distances correspond to roughly equal perceived color differences. + +## Functions + +### `rgb_to_lab` + +Convert 8-bit sRGB to CIELAB color space. + +```cpp +void rgb_to_lab( + const uint8_t r_u8, + const uint8_t g_u8, + const uint8_t b_u8, + double& out_l, + double& out_a, + double& out_b +); +``` + +#### Parameters + +| Parameter | Type | Range | Description | +| :-------- | :-------- | :----------- | :------------------------------ | +| `r_u8` | `uint8_t` | [0, 255] | Red channel (input) | +| `g_u8` | `uint8_t` | [0, 255] | Green channel (input) | +| `b_u8` | `uint8_t` | [0, 255] | Blue channel (input) | +| `out_l` | `double&` | [0, 100] | L\* lightness (output, clamped) | +| `out_a` | `double&` | ~[-128, 127] | a\* green-red axis (output) | +| `out_b` | `double&` | ~[-128, 127] | b\* blue-yellow axis (output) | + +#### Transformation Pipeline + +1. **sRGB → Linear RGB**: Inverse gamma correction (gamma expansion) + - Applies IEC 61966-2-1:1999 sRGB transfer function + - Converts [0, 255] → [0, 1] → linear [0, 1] + +2. **Linear RGB → XYZ**: Matrix multiplication + - Uses D65 illuminant (standard daylight, 6500K) + - Applies ITU-R BT.709 color primaries + +3. **XYZ → CIELAB**: Normalization and nonlinear transform + - Normalizes by D65 reference white point + - Applies CIE-defined piecewise function (cube root or linear near zero) + +#### Example + +```cpp +#include "cielab.h" + +// Convert bright red to LAB +uint8_t r = 255, g = 0, b = 0; +double L, a, b_lab; +rgb_to_lab(r, g, b, L, a, b_lab); +// Result: L ≈ 53.2, a ≈ 80.1, b ≈ 67.2 +``` + +--- + +### `lab_to_rgb` + +Convert CIELAB to 8-bit sRGB color space. + +```cpp +void lab_to_rgb( + const double L, + const double A, + const double B, + uint8_t& r_u8, + uint8_t& g_u8, + uint8_t& b_u8 +); +``` + +#### Parameters + +| Parameter | Type | Range | Description | +| :-------- | :--------- | :----------- | :------------------------------ | +| `L` | `double` | [0, 100] | L\* lightness (input) | +| `A` | `double` | ~[-128, 127] | a\* green-red axis (input) | +| `B` | `double` | ~[-128, 127] | b\* blue-yellow axis (input) | +| `r_u8` | `uint8_t&` | [0, 255] | Red channel (output, clamped) | +| `g_u8` | `uint8_t&` | [0, 255] | Green channel (output, clamped) | +| `b_u8` | `uint8_t&` | [0, 255] | Blue channel (output, clamped) | + +#### Transformation Pipeline + +1. **CIELAB → XYZ**: Inverse nonlinear transform and denormalization + - Applies inverse piecewise function (cube or linear) + - Denormalizes by D65 white point + +2. **XYZ → Linear RGB**: Inverse matrix multiplication + - May produce out-of-gamut values (negative or >1.0) + +3. **Linear RGB → sRGB**: Gamma correction + - Applies gamma compression using sRGB transfer function + - Clamps to [0, 1], rounds, and converts to [0, 255] + +#### Out-of-Gamut Handling + +:::warning Out-of-Gamut Colors +Not all LAB colors are representable in sRGB. Colors outside the sRGB gamut are clamped to the nearest valid RGB value, which may result in color shifts or loss of hue. +::: + +#### Example + +```cpp +#include "cielab.h" + +// Convert LAB back to RGB +double L = 53.2, a = 80.1, b = 67.2; +uint8_t r, g, b_rgb; +lab_to_rgb(L, a, b, r, g, b_rgb); +// Result: r ≈ 255, g ≈ 0, b ≈ 0 (bright red) +``` + +--- + +## Technical Specifications + +### Color Space Standards + +| Property | Value | +| :---------------- | :---------------------------- | +| **Color space** | sRGB (IEC 61966-2-1:1999) | +| **Illuminant** | D65 (6500K daylight) | +| **Observer** | CIE 1931 2° Standard Observer | +| **Gamma** | 2.4 (sRGB standard) | +| **RGB primaries** | ITU-R BT.709 | + +### CIELAB Coordinate System + +- **L\* (Lightness)**: Perceptual lightness + - 0 = black + - 100 = white + - 50 ≈ mid-gray +- **a\* (Green-Red)**: Color opponent dimension + - Negative values = green + - Positive values = red + - 0 = neutral (gray axis) +- **b\* (Blue-Yellow)**: Color opponent dimension + - Negative values = blue + - Positive values = yellow + - 0 = neutral (gray axis) + +### Distance Metric + +Euclidean distance in CIELAB space approximates perceptual color difference: + +$$ +\Delta E = \sqrt{(\Delta L^*)^2 + (\Delta a^*)^2 + (\Delta b^*)^2} +$$ + +**Interpretation**: + +- ΔE < 1: Imperceptible difference +- ΔE < 2: Perceptible with close observation +- ΔE < 10: Noticeable at a glance +- ΔE > 10: Significant color difference + +--- + +## Usage in Bilateral Filter + +The CIELAB color space is used in the bilateral filter to compute perceptually uniform range weights: + +```cpp +// In bilateral_filter.cpp (CIELAB mode) +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); +``` + +This produces more perceptually consistent smoothing compared to RGB Euclidean distance. + +--- + +## Performance Considerations + +### Computational Cost + +| Operation | Complexity | +| :----------- | :------------------------------- | +| `rgb_to_lab` | ~20-30 floating-point operations | +| `lab_to_rgb` | ~25-35 floating-point operations | + +**Key operations**: + +- Gamma correction: piecewise with `pow()` for nonlinear segment +- Matrix multiplication: 3×3 matrix +- XYZ transform: `cbrt()` or linear approximation + +### Optimization Notes + +- **Batch conversion**: When processing entire images, consider vectorization (SIMD) +- **LUT for gamma**: Can be precomputed for all 256 uint8_t values +- **Fast approximations**: Polynomial approximations for `pow()` and `cbrt()` can improve speed at slight accuracy cost + +--- + +## See Also + +- [Bilateral Filter Color Spaces](../../bilateral_filter/color-spaces) — How CIELAB is used in filtering +- [Bilateral Filter Implementation](../../bilateral_filter/implementation) — Implementation details +- [CIE 1976 L\*a\*b\* Color Space (Wikipedia)](https://en.wikipedia.org/wiki/CIELAB_color_space) +- [sRGB Specification (IEC 61966-2-1:1999)](https://www.color.org/chardata/rgb/srgb.xalter) diff --git a/docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md b/docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md new file mode 100644 index 000000000..de5588d1e --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md @@ -0,0 +1,187 @@ +--- +id: explained +title: Implementation Explained +sidebar_position: 1 +--- + +# RGB ↔ CIELAB Conversion Guide + +This explains the full mathematical conversion pipeline between **sRGB** and **CIELAB (Lab)** color spaces. + +# 1. Conversion Pipeline Overview + +## RGB → CIELAB + +```mermaid +flowchart LR + A[sRGB] --> B[Linear RGB] + B --> C["XYZ (D65)"] + C --> D[CIELAB] +``` + +1. sRGB → Linear RGB +2. Linear RGB → XYZ +3. XYZ → CIELAB + +## CIELAB → RGB + +```mermaid +flowchart LR + A[CIELAB] --> B["XYZ (D65)"] + B --> C[Linear RGB] + C --> D[sRGB] +``` + +1. CIELAB → XYZ +2. XYZ → Linear RGB +3. Linear RGB → sRGB + +# 2. sRGB to Linear RGB + +sRGB values are gamma‑compressed. Convert them to linear light: + +```math + +C_\text{lin} = +\begin{cases} +\frac{C_\text{srgb}}{12.92}, & C_\text{srgb} \le 0.04045 \\ +\left(\frac{C_\text{srgb} + 0.055}{1.055}\right)^{2.4}, & C_\text{srgb} > 0.04045 +\end{cases} + +``` + +This is applied independently to \(R\), \(G\), and \(B\). + +# 3. Linear RGB to XYZ + +Using the sRGB color space matrix with a D65 white point: + +```math + +\begin{bmatrix} +X \\ Y \\ Z +\end{bmatrix} += +\begin{bmatrix} +0.4124564 & 0.3575761 & 0.1804375 \\ +0.2126729 & 0.7151522 & 0.0721750 \\ +0.0193339 & 0.1191920 & 0.9503041 +\end{bmatrix} +\begin{bmatrix} +R_\text{lin} \\ G_\text{lin} \\ B_\text{lin} +\end{bmatrix} + +``` + +# 4. XYZ to CIELAB + +Normalize XYZ by the D65 reference white: + +```math +X_n = 0.95047,\quad Y_n = 1.00000,\quad Z_n = 1.08883 +``` + +```math +x = \frac{X}{X_n},\quad y = \frac{Y}{Y_n},\quad z = \frac{Z}{Z_n} +``` + +Define the nonlinear function: + +```math +f(t) = +\begin{cases} +t^{1/3}, & t > \left(\frac{6}{29}\right)^3 \\ +\frac{t}{3\left(\frac{6}{29}\right)^2} + \frac{4}{29}, & t \le \left(\frac{6}{29}\right)^3 +\end{cases} +``` + +Then compute Lab: + +```math +L^* = 116 f(y) - 16 +``` + +```math +a^* = 500 \left[f(x) - f(y)\right] +``` + +```math +b^* = 200 \left[f(y) - f(z)\right] +``` + +# 5. CIELAB to XYZ + +The inverse of \(f(t)\): + +```math +f^{-1}(t) = +\begin{cases} +t^3, & t > \frac{6}{29} \\ +3\left(\frac{6}{29}\right)^2 \left(t - \frac{4}{29}\right), & t \le \frac{6}{29} +\end{cases} +``` + +Compute: + +```math +f_y = \frac{L + 16}{116}, \quad +f_x = f_y + \frac{a}{500}, \quad +f_z = f_y - \frac{b}{200} +``` + +```math +X = X_n f^{-1}(f_x),\quad +Y = Y_n f^{-1}(f_y),\quad +Z = Z_n f^{-1}(f_z) +``` + +# 6. XYZ to Linear RGB + +```math +\begin{bmatrix} +R_\text{lin} \\ G_\text{lin} \\ B_\text{lin} +\end{bmatrix} += +\begin{bmatrix} + 3.240970 & -1.537383 & -0.498611 \\ +-0.969244 & 1.875968 & 0.041555 \\ + 0.055630 & -0.203977 & 1.056972 +\end{bmatrix} +\begin{bmatrix} +X \\ Y \\ Z +\end{bmatrix} +``` + +# 7. Linear RGB to sRGB + +```math +C_\text{srgb} = +\begin{cases} +12.92\, C_\text{lin}, & C_\text{lin} \le 0.0031308 \\ +1.055\, C_\text{lin}^{1/2.4} - 0.055, & C_\text{lin} > 0.0031308 +\end{cases} +``` + +Clamp results to \([0,1]\) and scaled by 255 before converting to 8‑bit. + +# 8. Summary + +## RGB → Lab + +- Remove gamma (sRGB → linear) +- Convert to XYZ +- Normalize by D65 +- Apply nonlinear transform +- Produce L\*, a\*, b\* + +## Lab → RGB + +- Convert Lab → XYZ via inverse nonlinear transform +- XYZ → linear RGB +- Linear RGB → sRGB (gamma) +- Clamp to valid output + +# 9. References + +- CIE 1976 L\*a\*b\* Specification +- IEC 61966‑2‑1 sRGB Standard diff --git a/docs/docs/reference/wasm/modules/image/overview.md b/docs/docs/reference/wasm/modules/image/overview.md index 929ad47b5..e227a11dd 100644 --- a/docs/docs/reference/wasm/modules/image/overview.md +++ b/docs/docs/reference/wasm/modules/image/overview.md @@ -26,6 +26,7 @@ src/wasm/modules/image/ │   ├── kmeans.h │   └── mergeSmallRegionsInPlace.h └── src + ├── bilateral_filter.cpp ├── fft_iterative.cpp ├── image_utils.cpp ├── kmeans.cpp @@ -41,6 +42,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. diff --git a/docs/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel.jsx b/docs/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel.jsx new file mode 100644 index 000000000..cbe12f6fb --- /dev/null +++ b/docs/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel.jsx @@ -0,0 +1,78 @@ +const RgbVsLabRangeKernel = () => ( +
+ + {/* Background */} + + + {/* RGB LUT area */} + { + const points = Array.from({ length: 101 }, (_, i) => { + const dx = 10 + i * 3; + const weight = Math.exp(-(i * i) / (2 * 20 * 20)); // sigma_range = 20 + const dy = 60 - weight * 50; + return `${dx},${dy}`; + }); + return ( + `M${points[0]} ` + + points + .slice(1) + .map((p) => `L${p}`) + .join(' ') + + ` L310,60 L10,60 Z` + ); + })()} + fill="rgba(255,107,107,0.3)" + stroke="#ff6b6b" + strokeWidth="2" + /> + + {/* CIELAB on-the-fly area */} + { + const points = Array.from({ length: 101 }, (_, i) => { + const dx = 10 + i * 3; + const weight = Math.exp(-(i * i) / (2 * 15 * 15)); // sigma_range = 15 + const dy = 120 - weight * 50; + return `${dx},${dy}`; + }); + return ( + `M${points[0]} ` + + points + .slice(1) + .map((p) => `L${p}`) + .join(' ') + + ` L310,120 L10,120 Z` + ); + })()} + fill="rgba(77,171,247,0.2)" + stroke="#4dabf7" + strokeWidth="2" + /> + + {/* Labels */} + + RGB LUT + + + CIELAB (on-the-fly) + + + {/* Axes */} + + + + Color difference ΔRGB + + + Color difference ΔLAB + + + + + Above, σr = 20 for RGB and σr = 15 for CIELAB. + +
+); + +export default RgbVsLabRangeKernel; diff --git a/src/components/WasmImageProcessor.jsx b/src/components/WasmImageProcessor.jsx index bb0448010..9616e84ab 100644 --- a/src/components/WasmImageProcessor.jsx +++ b/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); @@ -78,12 +78,17 @@ 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 imgBilateralFiltered = await bilateralFilter({ + pixels: fileData.pixels, + width, + height, + }); step(45); const thresholded = await blackThreshold({ ...fileData, - pixels: blurred, + pixels: imgBilateralFiltered, num_colors: 8, }); @@ -133,7 +138,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( diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 690f2aff5..bb48c0a26 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,6 +35,18 @@ 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, + color_space = 0, + }) => { + return ( + await call('bilateral_filter', { pixels, width, height, sigma_spatial, sigma_range, color_space }, ['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 +58,5 @@ export function useWasmWorker() { .output.pixels; }; - return { call, gaussianBlur, blackThreshold, kmeans, mergeSmallRegionsInPlace }; + return { call, gaussianBlur, bilateralFilter, blackThreshold, kmeans, mergeSmallRegionsInPlace }; } diff --git a/src/hooks/useWasmWorker.test.js b/src/hooks/useWasmWorker.test.js index ab3d00f3e..54da58936 100644 --- a/src/hooks/useWasmWorker.test.js +++ b/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()); diff --git a/src/wasm/modules/image/include/bilateral_filter.h b/src/wasm/modules/image/include/bilateral_filter.h new file mode 100644 index 000000000..2b5b938fe --- /dev/null +++ b/src/wasm/modules/image/include/bilateral_filter.h @@ -0,0 +1,24 @@ +#ifndef BILATERAL_FILTER_H +#define BILATERAL_FILTER_H + +#include // for size_t +#include // 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, + uint8_t color_space); + +} // namespace bilateral + +#endif // BILATERAL_FILTER_H diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h new file mode 100644 index 000000000..31625f9b4 --- /dev/null +++ b/src/wasm/modules/image/include/cielab.h @@ -0,0 +1,11 @@ +#ifndef CIELAB_H +#define CIELAB_H + +#include + +void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, + double &out_l, double &out_a, double &out_b); + +void lab_to_rgb(const double L, const double A, const double B, uint8_t &r_u8, + uint8_t &g_u8, uint8_t &b_u8); +#endif // CIELAB_H diff --git a/src/wasm/modules/image/include/kmeans.h b/src/wasm/modules/image/include/kmeans.h index 7dc18891f..20bbef4f3 100644 --- a/src/wasm/modules/image/include/kmeans.h +++ b/src/wasm/modules/image/include/kmeans.h @@ -31,6 +31,6 @@ EXPORTED void kmeans_clustering(uint8_t *data, int width, int height, int k, int max_iter); EXPORTED void kmeans_clustering_spatial(uint8_t *data, int width, int height, int k, int max_iter, - float spatial_weight); + float spatial_weight = 1.0); #endif diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp new file mode 100644 index 000000000..dfe807f7f --- /dev/null +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -0,0 +1,226 @@ +#include "bilateral_filter.h" +#include "cielab.h" +#include "exported.h" + +#include +#include +#include +#include +#include + +namespace bilateral { + +static constexpr double SIGMA_RADIUS_FACTOR{3.0}; // 3 standard deviations +static constexpr int MAX_KERNEL_RADIUS{50}; +// Max possible squared Euclidean distance in a 3-channel 8-bit image: 255^2 * 3 +// = 195075 Means max delta between images (imageA - imageB) in RGB channels +// (255^2 * 3) +static constexpr int MAX_RGB_DIST_SQ{255 * 255 * 3}; +static constexpr uint8_t COLOR_SPACE_OPTION_CIELAB{0}; +static constexpr uint8_t COLOR_SPACE_OPTION_RGB{1}; + +inline double gaussian(double x, double sigma) { + return std::exp(-(x * x) / (2.0 * sigma * sigma)); +} + +/* +The Bilateral Filter applies a composite weight based on both spatial distance +and radiometric difference (intensity) to return an image that is smoothed while +preserving edges. It reduces noise in flat regions while preserving edges by +assigning near-zero weight to pixels across high-contrast boundaries. + +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) +- color_space: Color space selector + ├── 0: CIELAB + └── 1: RGB +*/ +void bilateral_filter(uint8_t *image, size_t width, size_t height, + double sigma_spatial, double sigma_range, + uint8_t color_space) { + // bad data -> return + if (sigma_spatial <= 0.0 || sigma_range <= 0.0 || width <= 0 || height <= 0) + return; + if (color_space != COLOR_SPACE_OPTION_CIELAB && + color_space != COLOR_SPACE_OPTION_RGB) + return; + + const int raw_radius{ + static_cast(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatial))}; + const int radius{std::min(raw_radius, MAX_KERNEL_RADIUS)}; + const int kernel_diameter{2 * radius + 1}; + + std::vector result(width * height * 4); + + std::vector spatial_weights(kernel_diameter * kernel_diameter); + + // Precompute Spatial Weights (Gaussian Kernel) + for (int ky{-radius}; ky <= radius; ++ky) { + for (int kx{-radius}; kx <= radius; ++kx) { + const double dist{static_cast(std::sqrt(kx * kx + ky * ky))}; + spatial_weights[(ky + radius) * kernel_diameter + (kx + radius)] = + gaussian(dist, sigma_spatial); + } + } + + // ========= RGB-only section start ========= + // Precompute Range Weights + std::vector range_lut; + if (color_space == COLOR_SPACE_OPTION_RGB) { + range_lut.resize(MAX_RGB_DIST_SQ + 1); + for (int i{0}; i <= MAX_RGB_DIST_SQ; ++i) { + range_lut[i] = gaussian(static_cast(std::sqrt(i)), sigma_range); + } + } + // ========= RGB-only section end ========= + + // ========= CIELAB section start ========= + // Compute full image RGB - CIELAB conversion + std::vector cie_image; + if (color_space == COLOR_SPACE_OPTION_CIELAB) { + cie_image.resize(width * height * 4); + + for (int y{0}; y < height; y++) { + for (int x{0}; x < width; x++) { + int center_idx{(y * static_cast(width) + x) * 4}; + uint8_t r0{image[center_idx]}; + uint8_t g0{image[center_idx + 1]}; + uint8_t b0{image[center_idx + 2]}; + uint8_t a0{image[center_idx + 3]}; + double L0, A0, B0; + rgb_to_lab(r0, g0, b0, L0, A0, B0); + + cie_image[center_idx] = L0; + cie_image[center_idx + 1] = A0; + cie_image[center_idx + 2] = B0; + cie_image[center_idx + 3] = + 0.0; // unused but keep for indexing purposes + } + } + } + // ========= CIELAB section end ========= + + int h{static_cast(height)}; + int w{static_cast(width)}; + for (int y{0}; y < h; ++y) { + for (int x{0}; x < w; ++x) { + size_t center_idx{(y * width + x) * 4}; + + uint8_t r0{image[center_idx]}; + uint8_t g0{image[center_idx + 1]}; + uint8_t b0{image[center_idx + 2]}; + uint8_t a0{image[center_idx + 3]}; + + // ========= CIELAB-only section start ========= + double L0, A0, B0; + if (color_space == COLOR_SPACE_OPTION_CIELAB) { + L0 = cie_image[center_idx]; + A0 = cie_image[center_idx + 1]; + B0 = cie_image[center_idx + 2]; + } + // ========= CIELAB-only section end ========= + + // in RGB mode represents r,g,b accumulators + // in CIELAB mode represents L,A,B accumulators + double weight_acc_channel_0{0.0}, weight_acc_channel_1{0.0}, + weight_acc_channel_2{0.0}; + + double weight_acc{0.0}; + double w_space, w_range; + double dL, dA, dB, dist; + + for (int ky{-radius}; ky <= radius; ++ky) { + int ny{std::clamp(y + ky, 0, h - 1)}; + + for (int kx{-radius}; kx <= radius; ++kx) { + int nx{std::clamp(x + kx, 0, w - 1)}; + + size_t neighbor_idx{(ny * width + nx) * 4}; + + uint8_t r{image[neighbor_idx]}; + uint8_t g{image[neighbor_idx + 1]}; + uint8_t b{image[neighbor_idx + 2]}; + + double L{cie_image[neighbor_idx]}; + double A{cie_image[neighbor_idx + 1]}; + double B{cie_image[neighbor_idx + 2]}; + + w_space = + spatial_weights[(ky + radius) * kernel_diameter + (kx + radius)]; + + switch (color_space) { + case COLOR_SPACE_OPTION_RGB: { + const int dr{static_cast(r) - r0}; + const int dg{static_cast(g) - g0}; + const int db{static_cast(b) - b0}; + const int dist_sq{dr * dr + dg * dg + db * db}; + w_range = range_lut[dist_sq]; + + weight_acc_channel_0 += r * w_space * w_range; + weight_acc_channel_1 += g * w_space * w_range; + weight_acc_channel_2 += b * w_space * w_range; + break; + } + case COLOR_SPACE_OPTION_CIELAB: { + dL = L - L0; + dA = A - A0; + dB = B - B0; + + dist = std::sqrt(dL * dL + dA * dA + dB * dB); + w_range = gaussian(dist, sigma_range); + + weight_acc_channel_0 += L * w_space * w_range; + weight_acc_channel_1 += A * w_space * w_range; + weight_acc_channel_2 += B * w_space * w_range; + break; + } + } + + weight_acc += w_space * w_range; + } + } + + switch (color_space) { + case COLOR_SPACE_OPTION_RGB: { + result[center_idx] = static_cast( + std::clamp(weight_acc_channel_0 / weight_acc, 0.0, 255.0)); + result[center_idx + 1] = static_cast( + std::clamp(weight_acc_channel_1 / weight_acc, 0.0, 255.0)); + result[center_idx + 2] = static_cast( + std::clamp(weight_acc_channel_2 / weight_acc, 0.0, 255.0)); + result[center_idx + 3] = a0; + break; + } + case COLOR_SPACE_OPTION_CIELAB: { + double L{weight_acc_channel_0 / weight_acc}; + double A{weight_acc_channel_1 / weight_acc}; + double B{weight_acc_channel_2 / weight_acc}; + uint8_t r, g, b; + lab_to_rgb(L, A, B, r, g, b); + result[center_idx] = r; + result[center_idx + 1] = g; + result[center_idx + 2] = b; + result[center_idx + 3] = a0; + break; + } + } + } + } + + std::memcpy(image, result.data(), result.size()); +} + +} // namespace bilateral + +// Global wrapper for WASM export +EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, + double sigma_spatial, double sigma_range, + uint8_t color_space) { + bilateral::bilateral_filter(image, width, height, sigma_spatial, sigma_range, + color_space); +} diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp new file mode 100644 index 000000000..073594f4e --- /dev/null +++ b/src/wasm/modules/image/src/cielab.cpp @@ -0,0 +1,167 @@ +#include "cielab.h" +#include +#include + +// ====== Used in xyz_to_lab ======= +constexpr double DELTA{6.0 / 29.0}; // 0.2068966 +constexpr double DELTA_CUBED{DELTA * DELTA * DELTA}; // 0.008856 +constexpr double KAPPA{1.0 / (3.0 * DELTA * DELTA)}; // 7.787 +constexpr double EPSILON{16.0 / 116.0}; // 0.137931 + +// ====== Used in srgb_to_linear ====== +constexpr double SRGB_LINEAR_THRESHOLD{0.04045}; // linear segment boundary +constexpr double SRGB_LINEAR_FACTOR{12.92}; // scale factor for linear segment +constexpr double SRGB_GAMMA_OFFSET{0.055}; // offset for nonlinear segment +constexpr double SRGB_GAMMA{2.4}; // gamma exponent for nonlinear segment +constexpr double SRGB_GAMMA_INV{ + 1.0 / SRGB_GAMMA}; // gamma exponent for nonlinear segment + +/* + * ====== Used in rgb_to_lab ====== + * + * sRGB Khronos/W3C Transformation Matrices (D65 illuminant) + * + * sRGB → CIE XYZ: + * ┌ ┐ ┌ ┐ ┌ ┐ + * │ X │ │ 0.4124564 0.3575761 0.1804375 │ │ R │ + * │ Y │ = │ 0.2126729 0.7151522 0.0721750 │ │ G │ + * │ Z │ │ 0.0193339 0.1191920 0.9503041 │ │ B │ + * └ ┘ └ ┘ └ ┘ + * + * Reference: ITU-R BT.709 / sRGB standard (IEC 61966-2-1:1999) + */ +constexpr double SRGB_R_TO_X{0.4124564}; +constexpr double SRGB_G_TO_X{0.3575761}; +constexpr double SRGB_B_TO_X{0.1804375}; +constexpr double SRGB_R_TO_Y{0.2126729}; +constexpr double SRGB_G_TO_Y{0.7151522}; +constexpr double SRGB_B_TO_Y{0.0721750}; +constexpr double SRGB_R_TO_Z{0.0193339}; +constexpr double SRGB_G_TO_Z{0.1191920}; +constexpr double SRGB_B_TO_Z{0.9503041}; + +/* + * ====== Used in lab_to_rgb ====== + * + * sRGB Khronos/W3C Transformation Matrices (D65 illuminant) + * Inverse transformation (XYZ → sRGB) per Khronos/W3C, D65 white, slightly + * different from original BT.709 inverse. + * + * CIE XYZ → sRGB (inverse): + * ┌ ┐ ┌ ┐ ┌ ┐ + * │ R │ │ 3.240970 -1.537383 -0.498611 │ │ X │ + * │ G │ = │ -0.969244 1.875968 0.041555 │ │ Y │ + * │ B │ │ 0.055630 -0.203977 1.056972 │ │ Z │ + * └ ┘ └ ┘ └ ┘ + * + * Reference: ITU-R BT.709 / sRGB standard (IEC 61966-2-1:1999) + */ +constexpr double SRGB_X_TO_R{3.240970}; +constexpr double SRGB_Y_TO_R{-1.537383}; +constexpr double SRGB_Z_TO_R{-0.498611}; +constexpr double SRGB_X_TO_G{-0.969244}; +constexpr double SRGB_Y_TO_G{1.875968}; +constexpr double SRGB_Z_TO_G{0.041555}; +constexpr double SRGB_X_TO_B{0.055630}; +constexpr double SRGB_Y_TO_B{-0.203977}; +constexpr double SRGB_Z_TO_B{1.056972}; + +// Reference white point for D65 illuminant +constexpr double D65_Xn{0.95047}; +constexpr double D65_Yn{1.0}; +constexpr double D65_Zn{1.08883}; + +constexpr double LAB_L_FACTOR{116.0}; +constexpr double LAB_L_OFFSET{16.0}; +constexpr double LAB_A_FACTOR{500.0}; +constexpr double LAB_B_FACTOR{200.0}; + +// Function for the non-linear XYZ to Lab transformation +inline double xyz_to_lab(const double t) { + // prevent negative due to tiny floating errors + const double safe_t{std::max(0.0, t)}; + return (safe_t > DELTA_CUBED) ? std::cbrt(safe_t) + : (KAPPA * safe_t) + EPSILON; +} + +// Function for the non-linear sRGB to linear RGB transformation (inverse gamma +// correction) +inline double srgb_to_linear(const double c) { + const double safe_c{std::clamp(c, 0.0, 1.0)}; + return (safe_c <= SRGB_LINEAR_THRESHOLD) + ? safe_c / SRGB_LINEAR_FACTOR + : std::pow((safe_c + SRGB_GAMMA_OFFSET) / + (1.0 + SRGB_GAMMA_OFFSET), + SRGB_GAMMA); +} + +void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, + double &out_l, double &out_a, double &out_b) { + // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] + double r{srgb_to_linear(r_u8 / 255.0)}; + double g{srgb_to_linear(g_u8 / 255.0)}; + double b{srgb_to_linear(b_u8 / 255.0)}; + + // 2. Convert linear RGB to CIE XYZ (using D65 white point reference) + // The matrix below is for sRGB to XYZ (D65) + const double x{SRGB_R_TO_X * r + SRGB_G_TO_X * g + SRGB_B_TO_X * b}; + const double y{SRGB_R_TO_Y * r + SRGB_G_TO_Y * g + SRGB_B_TO_Y * b}; + const double z{SRGB_R_TO_Z * r + SRGB_G_TO_Z * g + SRGB_B_TO_Z * b}; + + // Normalize XYZ values by the white point + const double Xr{x / D65_Xn}; + const double Yr{y / D65_Yn}; + const double Zr{z / D65_Zn}; + + // 3. Convert CIE XYZ to CIE L*a*b* + const double fx{xyz_to_lab(Xr)}; + const double fy{xyz_to_lab(Yr)}; + const double fz{xyz_to_lab(Zr)}; + + // 4. Output values + out_l = LAB_L_FACTOR * fy - LAB_L_OFFSET; + out_a = LAB_A_FACTOR * (fx - fy); + out_b = LAB_B_FACTOR * (fy - fz); + + out_l = std::clamp(out_l, 0.0, 100.0); +} + +constexpr double inverse_xyz_to_lab(double t) { + return (t > DELTA) ? (t * t * t) : (3 * DELTA * DELTA * (t - EPSILON)); +} + +inline double gamma_encode(double u) { + // Guard against negative values from out-of-gamut colors + u = std::max(0.0, u); + return (u <= SRGB_LINEAR_THRESHOLD / SRGB_LINEAR_FACTOR) + ? SRGB_LINEAR_FACTOR * u + : (1.0 + SRGB_GAMMA_OFFSET) * std::pow(u, SRGB_GAMMA_INV) - + SRGB_GAMMA_OFFSET; +} + +void lab_to_rgb(const double L, const double A, const double B, + uint8_t &out_r_u8, uint8_t &out_g_u8, uint8_t &out_b_u8) { + // --- Lab → XYZ (D65 white point) + const double fy{(L + LAB_L_OFFSET) / LAB_L_FACTOR}; + const double fx{fy + A / LAB_A_FACTOR}; + const double fz{fy - B / LAB_B_FACTOR}; + + const double X{D65_Xn * inverse_xyz_to_lab(fx)}; + const double Y{D65_Yn * inverse_xyz_to_lab(fy)}; + const double Z{D65_Zn * inverse_xyz_to_lab(fz)}; + + // --- XYZ → linear RGB (sRGB) + double r{SRGB_X_TO_R * X + SRGB_Y_TO_R * Y + SRGB_Z_TO_R * Z}; + double g{SRGB_X_TO_G * X + SRGB_Y_TO_G * Y + SRGB_Z_TO_G * Z}; + double b{SRGB_X_TO_B * X + SRGB_Y_TO_B * Y + SRGB_Z_TO_B * Z}; + + // --- linear RGB → sRGB (gamma correction) + r = gamma_encode(std::clamp(r, 0.0, 1.0)); + g = gamma_encode(std::clamp(g, 0.0, 1.0)); + b = gamma_encode(std::clamp(b, 0.0, 1.0)); + + // --- Clamp and convert to 8-bit + out_r_u8 = static_cast(std::round(255.0 * std::clamp(r, 0.0, 1.0))); + out_g_u8 = static_cast(std::round(255.0 * std::clamp(g, 0.0, 1.0))); + out_b_u8 = static_cast(std::round(255.0 * std::clamp(b, 0.0, 1.0))); +} diff --git a/src/wasm/modules/image/src/kmeans.cpp b/src/wasm/modules/image/src/kmeans.cpp index d1951ac66..85ab85962 100644 --- a/src/wasm/modules/image/src/kmeans.cpp +++ b/src/wasm/modules/image/src/kmeans.cpp @@ -125,13 +125,12 @@ void kmeans_clustering_spatial(uint8_t *data, int width, int height, int k, for (int i = 0; i < height; ++i) { for (int j = 0; j < width; ++j) { int idx = i * width + j; - pixels[idx] = { - static_cast(data[idx * 3 + 0]), - static_cast(data[idx * 3 + 1]), - static_cast(data[idx * 3 + 2]), - static_cast(j), // x - static_cast(i) // y - }; + pixels[idx] = RGBXY{.r = static_cast(data[idx * 4 + 0]) / + 255, // normalize 0 -1 + .g = static_cast(data[idx * 4 + 1]) / 255, + .b = static_cast(data[idx * 4 + 2]) / 255, + .x = static_cast(j) / width, // normalize 0 - 1 + .y = static_cast(i) / height}; } } @@ -191,12 +190,11 @@ void kmeans_clustering_spatial(uint8_t *data, int width, int height, int k, } } } - - // Assign clustered colors back to data + // Assign clustered colors back to data (rescale pixel values 0 - 255) for (int i = 0; i < num_pixels; ++i) { int cluster = labels[i]; - data[i * 3 + 0] = static_cast(centroids[cluster].r); - data[i * 3 + 1] = static_cast(centroids[cluster].g); - data[i * 3 + 2] = static_cast(centroids[cluster].b); + data[i * 4 + 0] = static_cast(centroids[cluster].r * 255); + data[i * 4 + 1] = static_cast(centroids[cluster].g * 255); + data[i * 4 + 2] = static_cast(centroids[cluster].b * 255); } }