From 57863c431996ac3197bd39d00bd225e9ab0c8b4c Mon Sep 17 00:00:00 2001 From: Krasner Date: Wed, 31 Dec 2025 16:18:47 +0000 Subject: [PATCH 01/53] bilateral filter with CIELAB distance --- src/components/WasmImageProcessor.jsx | 7 +- src/hooks/useWasmWorker.js | 5 +- src/hooks/useWasmWorker.test.js | 1 + src/wasm/modules/image/include/cielab.h | 60 ++++++++++++ src/wasm/modules/image/include/image_utils.h | 3 + src/wasm/modules/image/src/image_utils.cpp | 99 ++++++++++++++++++++ 6 files changed, 171 insertions(+), 4 deletions(-) create mode 100644 src/wasm/modules/image/include/cielab.h diff --git a/src/components/WasmImageProcessor.jsx b/src/components/WasmImageProcessor.jsx index 715456a1d..5ce3f750e 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 { gaussianBlur, bilateralFilter, blackThreshold, kmeans, mergeSmallRegionsInPlace } = useWasmWorker(); const [originalSrc, setOriginalSrc] = useState(null); const [fileData, setFileData] = useState(null); @@ -81,7 +81,8 @@ const WasmImageProcessor = () => { const { width, height } = fileData; step(20); - const blurred = await gaussianBlur(fileData); + // const blurred = await gaussianBlur(fileData); + const blurred = await bilateralFilter(fileData); step(45); const thresholded = await blackThreshold({ @@ -139,7 +140,7 @@ const WasmImageProcessor = () => { step(0); }, 800); } - }, [fileData, gaussianBlur, blackThreshold, kmeans, mergeSmallRegionsInPlace, navigate, step]); + }, [fileData, gaussianBlur, 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..00b2058b3 100644 --- a/src/hooks/useWasmWorker.js +++ b/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_pixels = width * 0.005, sigma_range = 5.0 }) => { + return (await call('bilateral_filter', { pixels, width, height, sigma_pixels, 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 }; } diff --git a/src/hooks/useWasmWorker.test.js b/src/hooks/useWasmWorker.test.js index ab3d00f3e..cc742c194 100644 --- a/src/hooks/useWasmWorker.test.js +++ b/src/hooks/useWasmWorker.test.js @@ -65,6 +65,7 @@ describe('useWasmWorker', () => { 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'); diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h new file mode 100644 index 000000000..c2ef4ea4a --- /dev/null +++ b/src/wasm/modules/image/include/cielab.h @@ -0,0 +1,60 @@ +#include +#include +#include + +// Function for the non-linear XYZ to Lab transformation +double f_xyz(double t) { + if (t > 0.008856) { + return std::pow(t, 1.0/3.0); + } else { + return (7.787 * t) + (16.0 / 116.0); + } +} + +// Function for the non-linear sRGB to linear RGB transformation (inverse gamma correction) +double inverse_gamma(double c) { + if (c > 0.04045) { + return std::pow((c + 0.055) / 1.055, 2.4); + } else { + return c / 12.92; + } +} + +void rgb_to_lab(unsigned char r_u8, unsigned char g_u8, unsigned char b_u8, double L, double A, double B) { + // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] + double r = r_u8 / 255.0; + double g = g_u8 / 255.0; + double b = b_u8 / 255.0; + + r = inverse_gamma(r); + g = inverse_gamma(g); + b = inverse_gamma(b); + + // 2. Convert linear RGB to CIE XYZ (using D65 white point reference) + // The matrix below is for sRGB to XYZ (D65) + double x = 0.4124564 * r + 0.3575761 * g + 0.1804375 * b; + double y = 0.2126729 * r + 0.7151522 * g + 0.0721750 * b; + double z = 0.0193339 * r + 0.1191920 * g + 0.9503041 * b; + + // Reference white point for D65 illuminant + const double Xn = 0.95047; + const double Yn = 1.00000; + const double Zn = 1.08883; + + // Normalize XYZ values by the white point + double Xr = x / Xn; + double Yr = y / Yn; + double Zr = z / Zn; + + // 3. Convert CIE XYZ to CIE L*a*b* + double fx = f_xyz(Xr); + double fy = f_xyz(Yr); + double fz = f_xyz(Zr); + + L = 116.0 * fy - 16.0; + A = 500.0 * (fx - fy); + B = 200.0 * (fy - fz); + + // Clamp L channel to standard range [0, 100] + L = std::max(0.0, std::min(100.0, L)); +} \ No newline at end of file diff --git a/src/wasm/modules/image/include/image_utils.h b/src/wasm/modules/image/include/image_utils.h index da86ab524..c16df0c75 100644 --- a/src/wasm/modules/image/include/image_utils.h +++ b/src/wasm/modules/image/include/image_utils.h @@ -14,6 +14,9 @@ uint8_t quantize(uint8_t value, uint8_t region_size); EXPORTED void gaussian_blur_fft(uint8_t *image, size_t width, size_t height, double sigma); + +EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range); + EXPORTED void invert_image(uint8_t *ptr, int width, int height); EXPORTED void threshold_image(uint8_t *ptr, const int width, const int height, diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index 9e06b8f72..d7138dc4e 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -1,5 +1,6 @@ #include "image_utils.h" #include "fft_iterative.h" +#include "cielab.h" #include #include @@ -7,6 +8,104 @@ #include #include +double gaussian(float x, double sigma) { + return exp(-(pow(x, 2))/(2 * pow(sigma, 2))) / (2 * M_PI * pow(sigma, 2)); +} + +void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range) +{ + if (!image || width == 0 || height == 0 || sigma_pixels <= 0 || sigma_range <= 0) + return; + + const int radius = static_cast(1.5 * sigma_range); + const size_t diameter = radius * 2 + 1; + + // precompute + double range_filter[diameter * diameter]; + for (int i = 0; i < diameter; i++){ + for (int j = 0; j < diameter; j++){ + float dist = static_cast(sqrt(pow(i - radius, 2) + pow(j - radius, 2))); + range_filter[i*diameter + j] = gaussian(dist, sigma_range); + } + } + + uint8_t result[4 * height * width]; + + for (int i = 2; i < height - 2; i++){ + for (int j = 2; j < width - 2; j++) { + int center_index = 4 * (i * width + j); + uint8_t r0 = image[center_index]; + uint8_t g0 = image[center_index + 1]; + uint8_t b0 = image[center_index + 2]; + uint8_t a0 = image[center_index + 3]; + + double L0, A0, B0; + rgb_to_lab(r0, g0, b0, L0, A0, B0); + + double rf = 0.0; + double gf = 0.0; + double bf = 0.0; + double rW = 0.0; + double gW = 0.0; + double bW = 0.0; + + for (int ki = -radius; ki < radius; ki++){ + for (int kj = -radius; kj < radius; kj ++){ + int _i = i + ki; + int _j = j + kj; + if (_i < 0) + _i = 0; + if (_i > height - 1) + _i = height - 1; + if (_j < 0) + _j = 0; + if (_j < width - 1) + _j = width - 1; + int index = (_i * width) + _j; // std::clamp(i + ki, 0, height-1) * width + std::clamp(j + kj, 0, width-1); + + uint8_t r = image[index]; + uint8_t g = image[index + 1]; + uint8_t b = image[index + 2]; + + // independent weighting per channel + //double wr = gaussian(static_cast(r-r0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wg = gaussian(static_cast(g-g0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wb = gaussian(static_cast(b-b0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + + // euclidean color distance + //float dist = sqrt(pow(static_cast(r-r0), 2) + pow(static_cast(g-g0), 2) + pow(static_cast(b-b0), 2)); + + double L, A, B; + rgb_to_lab(r, g, b, L, A, B); + float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); + + double w_euc = gaussian(dist, sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + double wr = w_euc; + double wg = w_euc; + double wb = w_euc; + + rf += r * wr; + rW += wr; + + gf += g * wg; + gW += wg; + + bf += b * wb; + bW += wb; + } + } + + result[center_index] = static_cast(rf / rW); + result[center_index + 1] = static_cast(gf / gW); + result[center_index + 2] = static_cast(bf / bW); + result[center_index + 3] = a0; + } + } + + image = result; +} + + // image: pointer to RGBA data // width, height: dimensions // sigma: standard deviation of Gaussian blur From d32d1cf7d0710e859952eb40010c9a2320929cd6 Mon Sep 17 00:00:00 2001 From: Krasner Date: Wed, 31 Dec 2025 20:11:55 +0000 Subject: [PATCH 02/53] fix nomenclature --- src/hooks/useWasmWorker.js | 2 +- src/wasm/modules/image/src/image_utils.cpp | 13 +++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 00b2058b3..986ea420f 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,7 +35,7 @@ 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_pixels = width * 0.005, sigma_range = 5.0 }) => { + const bilateralFilter = async ({ pixels, width, height, sigma_pixels = width * 0.005, sigma_range = 30.0 }) => { return (await call('bilateral_filter', { pixels, width, height, sigma_pixels, sigma_range }, ['pixels'])).output.pixels; }; const blackThreshold = async ({ pixels, width, height, num_colors }) => { diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index d7138dc4e..3486458c6 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -14,10 +14,11 @@ double gaussian(float x, double sigma) { void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range) { + // sigma_pixel = spatial kernel if (!image || width == 0 || height == 0 || sigma_pixels <= 0 || sigma_range <= 0) return; - const int radius = static_cast(1.5 * sigma_range); + const int radius = static_cast(1.5 * sigma_pixels); const size_t diameter = radius * 2 + 1; // precompute @@ -25,7 +26,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ for (int i = 0; i < diameter; i++){ for (int j = 0; j < diameter; j++){ float dist = static_cast(sqrt(pow(i - radius, 2) + pow(j - radius, 2))); - range_filter[i*diameter + j] = gaussian(dist, sigma_range); + range_filter[i*diameter + j] = gaussian(dist, sigma_pixels); } } @@ -68,9 +69,9 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ uint8_t b = image[index + 2]; // independent weighting per channel - //double wr = gaussian(static_cast(r-r0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wg = gaussian(static_cast(g-g0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wb = gaussian(static_cast(b-b0), sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wr = gaussian(static_cast(r-r0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wg = gaussian(static_cast(g-g0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wb = gaussian(static_cast(b-b0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; // euclidean color distance //float dist = sqrt(pow(static_cast(r-r0), 2) + pow(static_cast(g-g0), 2) + pow(static_cast(b-b0), 2)); @@ -79,7 +80,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ rgb_to_lab(r, g, b, L, A, B); float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); - double w_euc = gaussian(dist, sigma_pixels) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + double w_euc = gaussian(dist, sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; double wr = w_euc; double wg = w_euc; double wb = w_euc; From a6ee113df202591e2f47c71bfd1d4b475ffa3cf7 Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 06:11:18 +0000 Subject: [PATCH 03/53] fix bug --- src/wasm/modules/image/src/image_utils.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index 3486458c6..13c11678c 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -32,8 +32,8 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ uint8_t result[4 * height * width]; - for (int i = 2; i < height - 2; i++){ - for (int j = 2; j < width - 2; j++) { + for (int i = 0; i < height; i++){ + for (int j = 0; j < width; j++) { int center_index = 4 * (i * width + j); uint8_t r0 = image[center_index]; uint8_t g0 = image[center_index + 1]; @@ -60,7 +60,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ _i = height - 1; if (_j < 0) _j = 0; - if (_j < width - 1) + if (_j > width - 1) _j = width - 1; int index = (_i * width) + _j; // std::clamp(i + ki, 0, height-1) * width + std::clamp(j + kj, 0, width-1); @@ -78,6 +78,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ double L, A, B; rgb_to_lab(r, g, b, L, A, B); + // needs sqrt? float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); double w_euc = gaussian(dist, sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; From 60893f59362f5180eb12ee15f92b3ec42f713d82 Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 06:14:22 +0000 Subject: [PATCH 04/53] bug fix --- src/wasm/modules/image/src/image_utils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index 13c11678c..a5c7f9f9b 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -62,7 +62,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ _j = 0; if (_j > width - 1) _j = width - 1; - int index = (_i * width) + _j; // std::clamp(i + ki, 0, height-1) * width + std::clamp(j + kj, 0, width-1); + int index = 4 * (_i * width + _j); uint8_t r = image[index]; uint8_t g = image[index + 1]; From 04a7ff91c5734729cd582679e86b6e71d1ee13d3 Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 06:16:45 +0000 Subject: [PATCH 05/53] variable name fix --- src/wasm/modules/image/src/image_utils.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index a5c7f9f9b..758ef7350 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -22,11 +22,11 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ const size_t diameter = radius * 2 + 1; // precompute - double range_filter[diameter * diameter]; + double spatial_filter[diameter * diameter]; for (int i = 0; i < diameter; i++){ for (int j = 0; j < diameter; j++){ float dist = static_cast(sqrt(pow(i - radius, 2) + pow(j - radius, 2))); - range_filter[i*diameter + j] = gaussian(dist, sigma_pixels); + spatial_filter[i*diameter + j] = gaussian(dist, sigma_pixels); } } @@ -69,9 +69,9 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ uint8_t b = image[index + 2]; // independent weighting per channel - //double wr = gaussian(static_cast(r-r0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wg = gaussian(static_cast(g-g0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wb = gaussian(static_cast(b-b0), sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wr = gaussian(static_cast(r-r0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wg = gaussian(static_cast(g-g0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; + //double wb = gaussian(static_cast(b-b0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; // euclidean color distance //float dist = sqrt(pow(static_cast(r-r0), 2) + pow(static_cast(g-g0), 2) + pow(static_cast(b-b0), 2)); @@ -81,7 +81,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ // needs sqrt? float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); - double w_euc = gaussian(dist, sigma_range) * range_filter[ (ki + radius) * diameter + (kj + radius) ]; + double w_euc = gaussian(dist, sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; double wr = w_euc; double wg = w_euc; double wb = w_euc; From 0021a76d841dcf48e424761347952b41ff976701 Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 06:19:56 +0000 Subject: [PATCH 06/53] incorporate git actions advice --- src/wasm/modules/image/include/cielab.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index c2ef4ea4a..baa8a1705 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -20,7 +20,7 @@ double inverse_gamma(double c) { } } -void rgb_to_lab(unsigned char r_u8, unsigned char g_u8, unsigned char b_u8, double L, double A, double B) { +void rgb_to_lab(unsigned char r_u8, unsigned char g_u8, unsigned char b_u8, double& L, double& A, double& B) { // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] double r = r_u8 / 255.0; double g = g_u8 / 255.0; From e97becb939974e0528cf80f876e04fcce6b5305f Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 21:51:15 +0000 Subject: [PATCH 07/53] cleanup --- src/hooks/useWasmWorker.js | 2 +- src/wasm/modules/image/include/cielab.h | 2 +- src/wasm/modules/image/src/image_utils.cpp | 22 +++++++++------------- 3 files changed, 11 insertions(+), 15 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 986ea420f..5b85e4600 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,7 +35,7 @@ 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_pixels = width * 0.005, sigma_range = 30.0 }) => { + const bilateralFilter = async ({ pixels, width, height, sigma_pixels = width * 0.005, sigma_range = 50.0 }) => { return (await call('bilateral_filter', { pixels, width, height, sigma_pixels, sigma_range }, ['pixels'])).output.pixels; }; const blackThreshold = async ({ pixels, width, height, num_colors }) => { diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index baa8a1705..fe9af6425 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -20,7 +20,7 @@ double inverse_gamma(double c) { } } -void rgb_to_lab(unsigned char r_u8, unsigned char g_u8, unsigned char b_u8, double& L, double& A, double& B) { +void rgb_to_lab(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B) { // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] double r = r_u8 / 255.0; double g = g_u8 / 255.0; diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index 758ef7350..a2d99c97d 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -8,7 +8,8 @@ #include #include -double gaussian(float x, double sigma) { +double evaluate_gaussian(float x, double sigma) { + // evaluates 1d gaussian function desribed by sigma at x return exp(-(pow(x, 2))/(2 * pow(sigma, 2))) / (2 * M_PI * pow(sigma, 2)); } @@ -26,7 +27,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ for (int i = 0; i < diameter; i++){ for (int j = 0; j < diameter; j++){ float dist = static_cast(sqrt(pow(i - radius, 2) + pow(j - radius, 2))); - spatial_filter[i*diameter + j] = gaussian(dist, sigma_pixels); + spatial_filter[i*diameter + j] = evaluate_gaussian(dist, sigma_pixels); } } @@ -51,7 +52,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ double bW = 0.0; for (int ki = -radius; ki < radius; ki++){ - for (int kj = -radius; kj < radius; kj ++){ + for (int kj = -radius; kj < radius; kj++){ int _i = i + ki; int _j = j + kj; if (_i < 0) @@ -68,20 +69,15 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ uint8_t g = image[index + 1]; uint8_t b = image[index + 2]; - // independent weighting per channel - //double wr = gaussian(static_cast(r-r0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wg = gaussian(static_cast(g-g0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; - //double wb = gaussian(static_cast(b-b0), sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; - - // euclidean color distance - //float dist = sqrt(pow(static_cast(r-r0), 2) + pow(static_cast(g-g0), 2) + pow(static_cast(b-b0), 2)); - + /* + as described in https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Tomasi98.pdf + use euclidean distance in LAB color space for less artifacts + */ double L, A, B; rgb_to_lab(r, g, b, L, A, B); - // needs sqrt? float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); - double w_euc = gaussian(dist, sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; + double w_euc = evaluate_gaussian(dist, sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; double wr = w_euc; double wg = w_euc; double wb = w_euc; From d7065e623636c310ce044f05b7876f3f17da4cdb Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 22:16:34 +0000 Subject: [PATCH 08/53] Fix indexing problem in kmeans_clustering_spatial. Verify that it works but don't call explicitly in useWasmWorker.json - commented out for future --- src/hooks/useWasmWorker.js | 4 ++++ src/wasm/modules/image/include/kmeans.h | 2 +- src/wasm/modules/image/src/kmeans.cpp | 21 ++++++++++----------- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 5b85e4600..84f90769c 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -44,6 +44,10 @@ export function useWasmWorker() { const kmeans = async ({ pixels, width, height, num_colors, max_iter = 100 }) => { return (await call('kmeans_clustering', { pixels, width, height, num_colors, max_iter }, ['pixels'])).output.pixels; }; + + /*const kmeans = async ({ pixels, width, height, num_colors, max_iter = 100, spatial_weight=0.1 }) => { + return (await call('kmeans_clustering_spatial', { pixels, width, height, num_colors, max_iter, spatial_weight }, ['pixels'])).output.pixels; + };*/ const mergeSmallRegionsInPlace = async ({ pixels, width, height, minArea, minWidth, minHeight }) => { return (await call('mergeSmallRegionsInPlace', { pixels, width, height, minArea, minWidth, minHeight }, ['pixels'])) .output.pixels; 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/kmeans.cpp b/src/wasm/modules/image/src/kmeans.cpp index d1951ac66..0a5f91f04 100644 --- a/src/wasm/modules/image/src/kmeans.cpp +++ b/src/wasm/modules/image/src/kmeans.cpp @@ -125,12 +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 +191,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); } } From ba234a8a55e4691761bb0077b7695e42684dbc09 Mon Sep 17 00:00:00 2001 From: Krasner Date: Thu, 1 Jan 2026 22:40:46 +0000 Subject: [PATCH 09/53] kernel range fix --- src/wasm/modules/image/src/image_utils.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index a2d99c97d..ca7426ac4 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -51,8 +51,8 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ double gW = 0.0; double bW = 0.0; - for (int ki = -radius; ki < radius; ki++){ - for (int kj = -radius; kj < radius; kj++){ + for (int ki = -radius; ki <= radius; ki++){ + for (int kj = -radius; kj <= radius; kj++){ int _i = i + ki; int _j = j + kj; if (_i < 0) From 0d4200582d4c6ff4f2e8af98c47679a496746661 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Tue, 30 Dec 2025 20:31:06 -0800 Subject: [PATCH 10/53] feat(WASM): add bilateral filter --- src/components/WasmImageProcessor.jsx | 15 ++- src/hooks/useWasmWorker.js | 5 +- src/wasm/modules/image/include/image_utils.h | 3 + .../modules/image/src/bilateral_filter.cpp | 103 ++++++++++++++++++ 4 files changed, 121 insertions(+), 5 deletions(-) create mode 100644 src/wasm/modules/image/src/bilateral_filter.cpp diff --git a/src/components/WasmImageProcessor.jsx b/src/components/WasmImageProcessor.jsx index 715456a1d..3fffba2c3 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); @@ -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({ + image: 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( diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 690f2aff5..252e8b793 100644 --- a/src/hooks/useWasmWorker.js +++ b/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 ({ image, width, height, sigma_spatial = 3.0, sigma_range = 50.0 }) => { + return (await call('bilateral_filter', { image, width, height, sigma_spatial, sigma_range }, ['image'])).output.image; + }; 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 }; } diff --git a/src/wasm/modules/image/include/image_utils.h b/src/wasm/modules/image/include/image_utils.h index da86ab524..9c549712b 100644 --- a/src/wasm/modules/image/include/image_utils.h +++ b/src/wasm/modules/image/include/image_utils.h @@ -21,4 +21,7 @@ EXPORTED void threshold_image(uint8_t *ptr, const int width, const int height, EXPORTED void black_threshold_image(uint8_t *ptr, const int width, const int height, const int num_thresholds); +EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, + double sigma_spatial, double sigma_range); + #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..acd05a9a7 --- /dev/null +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -0,0 +1,103 @@ +#include "image_utils.h" +#include +#include +#include +#include + +static constexpr double SIGMA_RADIUS_FACTOR = 3.0; +static constexpr int MAX_PIXEL_VAL = 255; +// 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; + +/* +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) +*/ +void bilateral_filter(uint8_t *image, size_t width, size_t height, + double sigma_spatial, double sigma_range) { + if (sigma_spatial <= 0.0 || sigma_range <= 0.0) return; + + const int radius = static_cast(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatial)); + const int kernel_width = 2 * radius + 1; + const size_t stride = width * 4; + std::vector result(width * height * 4); + + // NOTE: precompute Spatial Weights (Gaussian Kernel) + std::vector spatial_weights(kernel_width * kernel_width); + double two_sigma_space_sq = 2 * sigma_spatial * sigma_spatial; + + for (int ky = -radius; ky <= radius; ++ky) { + for (int kx = -radius; kx <= radius; ++kx) { + double dist2 = static_cast(kx * kx + ky * ky); + spatial_weights[(ky + radius) * kernel_width + (kx + radius)] = + std::exp(-dist2 / two_sigma_space_sq); + } + } + + // NOTE: precompute Range Weights + std::vector range_lut(MAX_RGB_DIST_SQ + 1); + double two_sigma_range_sq = 2 * sigma_range * sigma_range; + + for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) { + range_lut[i] = std::exp(-static_cast(i) / two_sigma_range_sq); + } + + 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]; + + double r_acc = 0.0, g_acc = 0.0, b_acc = 0.0, weight_acc = 0.0; + + 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 w_space = spatial_weights[(ky + radius) * kernel_width + (kx + radius)]; + + int dr = static_cast(r) - r0; + int dg = static_cast(g) - g0; + int db = static_cast(b) - b0; + int dist_sq = dr*dr + dg*dg + db*db; + + 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; + } + } + + result[center_idx] = static_cast(std::clamp(r_acc / weight_acc, 0.0, 255.0)); + result[center_idx + 1] = static_cast(std::clamp(g_acc / weight_acc, 0.0, 255.0)); + result[center_idx + 2] = static_cast(std::clamp(b_acc / weight_acc, 0.0, 255.0)); + result[center_idx + 3] = a0; + } + } + + std::memcpy(image, result.data(), result.size()); +} From cd4da10388bca78e73e25142ee0998564d49db62 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Tue, 30 Dec 2025 20:44:16 -0800 Subject: [PATCH 11/53] test: add call method test, and custom parameters test --- src/hooks/useWasmWorker.test.js | 53 +++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/hooks/useWasmWorker.test.js b/src/hooks/useWasmWorker.test.js index ab3d00f3e..266170747 100644 --- a/src/hooks/useWasmWorker.test.js +++ b/src/hooks/useWasmWorker.test.js @@ -202,6 +202,59 @@ describe('useWasmWorker', () => { }); }); + describe('bilateralFilter', () => { + it('should call worker with bilateral_filter function', async () => { + const { result } = renderHook(() => useWasmWorker()); + + const image = new Uint8ClampedArray([255, 0, 0, 255]); + const width = 1; + const height = 1; + + act(() => { + result.current.bilateralFilter({ image, width, height }); + }); + + expect(mockWorkerInstance.postMessage).toHaveBeenCalledWith( + expect.objectContaining({ + funcName: 'bilateral_filter', + args: expect.objectContaining({ + image, + width, + height, + sigma_spatial: 3.0, + sigma_range: 50.0, + }), + bufferKeys: ['image'], + }) + ); + }); + + it('should use custom sigma_spatial and sigma_range when provided', async () => { + const { result } = renderHook(() => useWasmWorker()); + + const image = new Uint8ClampedArray([255, 0, 0, 255]); + + act(() => { + result.current.bilateralFilter({ + image, + 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()); From 8d296f126c7985a89b79b8b01c0213e20703e632 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Tue, 30 Dec 2025 20:56:53 -0800 Subject: [PATCH 12/53] refactor: rename image parameter to pixels to match lib conventions --- src/components/WasmImageProcessor.jsx | 2 +- src/hooks/useWasmWorker.js | 4 ++-- src/hooks/useWasmWorker.test.js | 12 ++++++------ 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/components/WasmImageProcessor.jsx b/src/components/WasmImageProcessor.jsx index 3fffba2c3..c347c3ded 100644 --- a/src/components/WasmImageProcessor.jsx +++ b/src/components/WasmImageProcessor.jsx @@ -85,7 +85,7 @@ const WasmImageProcessor = () => { // const blurred = await gaussianBlur(fileData); const imgBilateralFiltered = await bilateralFilter({ - image: fileData.pixels, + pixels: fileData.pixels, width, height, }); diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 252e8b793..0a0e480d0 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,8 +35,8 @@ 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 ({ image, width, height, sigma_spatial = 3.0, sigma_range = 50.0 }) => { - return (await call('bilateral_filter', { image, width, height, sigma_spatial, sigma_range }, ['image'])).output.image; + 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; diff --git a/src/hooks/useWasmWorker.test.js b/src/hooks/useWasmWorker.test.js index 266170747..ba3adffbb 100644 --- a/src/hooks/useWasmWorker.test.js +++ b/src/hooks/useWasmWorker.test.js @@ -206,25 +206,25 @@ describe('useWasmWorker', () => { it('should call worker with bilateral_filter function', async () => { const { result } = renderHook(() => useWasmWorker()); - const image = new Uint8ClampedArray([255, 0, 0, 255]); + const pixels = new Uint8ClampedArray([255, 0, 0, 255]); const width = 1; const height = 1; act(() => { - result.current.bilateralFilter({ image, width, height }); + result.current.bilateralFilter({ pixels, width, height }); }); expect(mockWorkerInstance.postMessage).toHaveBeenCalledWith( expect.objectContaining({ funcName: 'bilateral_filter', args: expect.objectContaining({ - image, + pixels, width, height, sigma_spatial: 3.0, sigma_range: 50.0, }), - bufferKeys: ['image'], + bufferKeys: ['pixels'], }) ); }); @@ -232,11 +232,11 @@ describe('useWasmWorker', () => { it('should use custom sigma_spatial and sigma_range when provided', async () => { const { result } = renderHook(() => useWasmWorker()); - const image = new Uint8ClampedArray([255, 0, 0, 255]); + const pixels = new Uint8ClampedArray([255, 0, 0, 255]); act(() => { result.current.bilateralFilter({ - image, + pixels, width: 100, height: 100, sigma_spatial: 5.0, From c21ef2156391e5b51563a2a3bc971f50a3e49197 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Tue, 30 Dec 2025 21:01:00 -0800 Subject: [PATCH 13/53] test: add bilateralFilter as part of helper methods return test --- src/hooks/useWasmWorker.test.js | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/hooks/useWasmWorker.test.js b/src/hooks/useWasmWorker.test.js index ba3adffbb..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'); From a39db1627c1cf607d086826d1ff457aeb9c0ab02 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Tue, 30 Dec 2025 21:01:37 -0800 Subject: [PATCH 14/53] refactor: remove unused variable MAX_PIXEL_VAL --- src/wasm/modules/image/src/bilateral_filter.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index acd05a9a7..e82da0c7e 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -5,7 +5,6 @@ #include static constexpr double SIGMA_RADIUS_FACTOR = 3.0; -static constexpr int MAX_PIXEL_VAL = 255; // 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; From 6317a8192623addc723736aae4849e40dab2442a Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Wed, 31 Dec 2025 08:56:41 -0800 Subject: [PATCH 15/53] feat: add headers to be used by WASM (best practice) --- .../modules/image/include/bilateral_filter.h | 20 +++++++++++++++++++ .../modules/image/src/bilateral_filter.cpp | 16 +++++++++++++-- 2 files changed, 34 insertions(+), 2 deletions(-) create mode 100644 src/wasm/modules/image/include/bilateral_filter.h 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..b05b824b4 --- /dev/null +++ b/src/wasm/modules/image/include/bilateral_filter.h @@ -0,0 +1,20 @@ +#ifndef BILATERAL_FILTER_H +#define BILATERAL_FILTER_H + +#include // for size_t +#include // for uint8_t + +namespace bilateral { + +// Apply bilateral filter to an image. +// 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 diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index e82da0c7e..89c79321b 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -1,8 +1,13 @@ -#include "image_utils.h" +#include "bilateral_filter.h" +#include "exported.h" + #include #include #include #include +#include + +namespace bilateral { static constexpr double SIGMA_RADIUS_FACTOR = 3.0; // Max possible squared Euclidean distance in a 3-channel 8-bit image: 255^2 * 3 = 195075 @@ -25,7 +30,6 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, const int radius = static_cast(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatial)); const int kernel_width = 2 * radius + 1; - const size_t stride = width * 4; std::vector result(width * height * 4); // NOTE: precompute Spatial Weights (Gaussian Kernel) @@ -100,3 +104,11 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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) { + bilateral::bilateral_filter(image, width, height, sigma_spatial, sigma_range); +} From 47199957fa5b25ec2ae23349e3a789ee8752408e Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Wed, 31 Dec 2025 09:28:01 -0800 Subject: [PATCH 16/53] docs(WASM): add bilateral_filter documentation (overview, explained, implementation, and api) --- .../image/bilateral_filter/_category_.json | 10 +++ .../modules/image/bilateral_filter/api.md | 29 ++++++++ .../image/bilateral_filter/explained.md | 65 +++++++++++++++++ .../image/bilateral_filter/implementation.md | 73 +++++++++++++++++++ .../image/bilateral_filter/overview.md | 31 ++++++++ .../reference/wasm/modules/image/overview.md | 2 + 6 files changed, 210 insertions(+) create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/api.md create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md 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..2517950d2 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json @@ -0,0 +1,10 @@ +{ + "label": "bilateral_filter.h", + "position": 2, + "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" + } +} \ No newline at end of file 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..3451cc049 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -0,0 +1,29 @@ +--- +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. + +| Function | Signature | Purpose | +| :--- | :--- | :--- | +| `bilateral_filter` | `void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, double sigma_range)` | Applies a bilateral filter to an RGBA image. | + +## 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`. +::: 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..d98038924 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -0,0 +1,65 @@ +--- +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). + +## 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 +// 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. 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_s}$. +3. **Range Weight**: Calculate squared color distance $\|C_p - C_q\|^2$ and look up precomputed $G_{\sigma_r}$. +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..ab5c6b9fb --- /dev/null +++ b/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(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(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(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/overview.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md new file mode 100644 index 000000000..91705a4c7 --- /dev/null +++ b/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/) diff --git a/docs/docs/reference/wasm/modules/image/overview.md b/docs/docs/reference/wasm/modules/image/overview.md index c55f6a01a..001b834f2 100644 --- a/docs/docs/reference/wasm/modules/image/overview.md +++ b/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 image denoising. - `fft_iterative` — Fast Fourier Transform utilities. - Used by **Gaussian Blur** inside image_utils.h. - `kmeans` — K-means clustering used for quantization. From 9e71927046c6a67c05038def4e6d18b578ce3755 Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Wed, 31 Dec 2025 09:40:38 -0800 Subject: [PATCH 17/53] feat(bilateral_filter): add upper bound validation for sigma_spatial --- src/wasm/modules/image/src/bilateral_filter.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 89c79321b..7c7a19e35 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -9,7 +9,8 @@ namespace bilateral { -static constexpr double SIGMA_RADIUS_FACTOR = 3.0; +static constexpr double SIGMA_RADIUS_FACTOR = 3.0; +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; @@ -27,8 +28,10 @@ It reduces noise in flat regions while preserving edges by assigning near-zero w void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, double sigma_range) { if (sigma_spatial <= 0.0 || sigma_range <= 0.0) return; + if (width <= 0 || height <= 0) return; - const int radius = static_cast(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatial)); + 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_width = 2 * radius + 1; std::vector result(width * height * 4); From fc1af428eeef8fc1591a0dd41a5aa218659bb2ea Mon Sep 17 00:00:00 2001 From: Francisco Sanchez Date: Wed, 31 Dec 2025 09:41:26 -0800 Subject: [PATCH 18/53] docs: improve bilateral headers documentation --- src/wasm/modules/image/include/bilateral_filter.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/wasm/modules/image/include/bilateral_filter.h b/src/wasm/modules/image/include/bilateral_filter.h index b05b824b4..d6065cb49 100644 --- a/src/wasm/modules/image/include/bilateral_filter.h +++ b/src/wasm/modules/image/include/bilateral_filter.h @@ -7,6 +7,7 @@ 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) From 7e2c9426963431482f879b5d8e7ec27fc5e1cd5a Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:24:15 +0200 Subject: [PATCH 19/53] docs(bilateral filter): explain use of Gaussian kernels inside formula - Simple info admonition that explains the link to Gaussian functions --- .../modules/image/bilateral_filter/explained.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index d98038924..e54106593 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -29,6 +29,20 @@ Where: - $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. From bf2fb0aa3b54c02ff0018959bd34d1014a5b7653 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:24:45 +0200 Subject: [PATCH 20/53] docs(bilateral filter): better styling --- .../reference/wasm/modules/image/bilateral_filter/explained.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index e54106593..8e83f4808 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -53,8 +53,7 @@ Calculating `std::exp()` inside the inner loop is expensive. We precompute the t - **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 -// Precomputing Range Weights +```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); From 457bc287279f807842ea643ea8fd8e386ac719f0 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:25:04 +0200 Subject: [PATCH 21/53] docs(bilateral filter): better styling --- .../wasm/modules/image/bilateral_filter/explained.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index 8e83f4808..07be50365 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -65,8 +65,8 @@ for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) { 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_s}$. -3. **Range Weight**: Calculate squared color distance $\|C_p - C_q\|^2$ and look up precomputed $G_{\sigma_r}$. +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. From 7df317326f1073596de3bf50723d5e90f18684cc Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:25:31 +0200 Subject: [PATCH 22/53] docs(bilateral filter): correct sidebar_position --- .../wasm/modules/image/bilateral_filter/_category_.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json index 2517950d2..dea27276f 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json @@ -1,6 +1,6 @@ { "label": "bilateral_filter.h", - "position": 2, + "position": 4, "link": { "type": "generated-index", "title": "Bilateral Filter", From 273f207ab843553498a0fe00f1f2541e7d3c6567 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:25:52 +0200 Subject: [PATCH 23/53] docs(bilateral filter): mobile accessibility --- .../reference/wasm/modules/image/bilateral_filter/api.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index 3451cc049..0a23beb6f 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -9,9 +9,11 @@ sidebar_position: 5 Quick reference for the function implemented in the header. -| Function | Signature | Purpose | -| :--- | :--- | :--- | -| `bilateral_filter` | `void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, double sigma_range)` | Applies a bilateral filter to an RGBA image. | +```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 From 6dad78d3148642d3a4484a5881e9bac2471db676 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 17:26:28 +0200 Subject: [PATCH 24/53] docs(bilateral filter): explicit description in module overview --- docs/docs/reference/wasm/modules/image/overview.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/docs/reference/wasm/modules/image/overview.md b/docs/docs/reference/wasm/modules/image/overview.md index 001b834f2..3620261f8 100644 --- a/docs/docs/reference/wasm/modules/image/overview.md +++ b/docs/docs/reference/wasm/modules/image/overview.md @@ -39,7 +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 image denoising. +- `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. From a4c38497a20c7112e697b35b279f295e053d7c83 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 17:40:21 +0200 Subject: [PATCH 25/53] === end of commits from #176 === From c0e9d0761b57e915972bea30eaa7e5975676525f Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 19:44:41 +0200 Subject: [PATCH 26/53] fix(duplicate symbol: bilateral_filter): temp rename cielab -> bilateral_filter_cielab --- src/wasm/modules/image/include/image_utils.h | 2 +- src/wasm/modules/image/src/image_utils.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/wasm/modules/image/include/image_utils.h b/src/wasm/modules/image/include/image_utils.h index 3a373490c..3383c33ea 100644 --- a/src/wasm/modules/image/include/image_utils.h +++ b/src/wasm/modules/image/include/image_utils.h @@ -15,7 +15,7 @@ uint8_t quantize(uint8_t value, uint8_t region_size); EXPORTED void gaussian_blur_fft(uint8_t *image, size_t width, size_t height, double sigma); -EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range); +EXPORTED void bilateral_filter_cielab(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range); EXPORTED void invert_image(uint8_t *ptr, int width, int height); diff --git a/src/wasm/modules/image/src/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index ca7426ac4..b3c5d217b 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -13,12 +13,12 @@ double evaluate_gaussian(float x, double sigma) { return exp(-(pow(x, 2))/(2 * pow(sigma, 2))) / (2 * M_PI * pow(sigma, 2)); } -void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range) +void bilateral_filter_cielab(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range) { // sigma_pixel = spatial kernel if (!image || width == 0 || height == 0 || sigma_pixels <= 0 || sigma_range <= 0) return; - + const int radius = static_cast(1.5 * sigma_pixels); const size_t diameter = radius * 2 + 1; @@ -69,14 +69,14 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_ uint8_t g = image[index + 1]; uint8_t b = image[index + 2]; - /* + /* as described in https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Tomasi98.pdf use euclidean distance in LAB color space for less artifacts */ double L, A, B; rgb_to_lab(r, g, b, L, A, B); float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); - + double w_euc = evaluate_gaussian(dist, sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; double wr = w_euc; double wg = w_euc; From 7604897b2701db9d4075db7c773f7b1ab251ada2 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 20:16:44 +0200 Subject: [PATCH 27/53] docs(bilateral filter): fix api.md styling --- docs/docs/reference/wasm/modules/image/bilateral_filter/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index 0a23beb6f..d5a5276f2 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -14,6 +14,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, double sigma_range) +``` ## Parameters From 31887a72514751a4dc3082a0d96cc4ff556e043a Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 22:01:03 +0200 Subject: [PATCH 28/53] feat(bilateral filter): combine CIELAB & RGB implementations into single function --- src/hooks/useWasmWorker.js | 8 +- .../modules/image/include/bilateral_filter.h | 3 +- src/wasm/modules/image/include/image_utils.h | 6 - .../modules/image/src/bilateral_filter.cpp | 144 +++++++++++------- src/wasm/modules/image/src/image_utils.cpp | 97 ------------ 5 files changed, 93 insertions(+), 165 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index e6b060186..f1c502d22 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,8 +35,8 @@ 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 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; @@ -44,10 +44,6 @@ export function useWasmWorker() { const kmeans = async ({ pixels, width, height, num_colors, max_iter = 100 }) => { return (await call('kmeans_clustering', { pixels, width, height, num_colors, max_iter }, ['pixels'])).output.pixels; }; - - /*const kmeans = async ({ pixels, width, height, num_colors, max_iter = 100, spatial_weight=0.1 }) => { - return (await call('kmeans_clustering_spatial', { pixels, width, height, num_colors, max_iter, spatial_weight }, ['pixels'])).output.pixels; - };*/ const mergeSmallRegionsInPlace = async ({ pixels, width, height, minArea, minWidth, minHeight }) => { return (await call('mergeSmallRegionsInPlace', { pixels, width, height, minArea, minWidth, minHeight }, ['pixels'])) .output.pixels; diff --git a/src/wasm/modules/image/include/bilateral_filter.h b/src/wasm/modules/image/include/bilateral_filter.h index d6065cb49..b5d78d1e7 100644 --- a/src/wasm/modules/image/include/bilateral_filter.h +++ b/src/wasm/modules/image/include/bilateral_filter.h @@ -14,7 +14,8 @@ namespace bilateral { // - 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); + double sigma_spatial, double sigma_range, + uint8_t color_space); } // namespace bilateral diff --git a/src/wasm/modules/image/include/image_utils.h b/src/wasm/modules/image/include/image_utils.h index 3383c33ea..da86ab524 100644 --- a/src/wasm/modules/image/include/image_utils.h +++ b/src/wasm/modules/image/include/image_utils.h @@ -14,9 +14,6 @@ uint8_t quantize(uint8_t value, uint8_t region_size); EXPORTED void gaussian_blur_fft(uint8_t *image, size_t width, size_t height, double sigma); - -EXPORTED void bilateral_filter_cielab(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range); - EXPORTED void invert_image(uint8_t *ptr, int width, int height); EXPORTED void threshold_image(uint8_t *ptr, const int width, const int height, @@ -24,7 +21,4 @@ EXPORTED void threshold_image(uint8_t *ptr, const int width, const int height, EXPORTED void black_threshold_image(uint8_t *ptr, const int width, const int height, const int num_thresholds); -EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, - double sigma_spatial, double sigma_range); - #endif diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 7c7a19e35..bcd0e537b 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -1,5 +1,6 @@ #include "bilateral_filter.h" #include "exported.h" +#include "cielab.h" #include #include @@ -9,14 +10,21 @@ namespace bilateral { -static constexpr double SIGMA_RADIUS_FACTOR = 3.0; -static constexpr int MAX_KERNEL_RADIUS = 50; +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 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. +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: @@ -24,72 +32,97 @@ It reduces noise in flat regions while preserving edges by assigning near-zero w - 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) { - if (sigma_spatial <= 0.0 || sigma_range <= 0.0) return; - if (width <= 0 || height <= 0) 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_width = 2 * radius + 1; + 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; + + 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); - // NOTE: precompute Spatial Weights (Gaussian Kernel) - std::vector spatial_weights(kernel_width * kernel_width); - double two_sigma_space_sq = 2 * sigma_spatial * sigma_spatial; + std::vector spatial_weights(kernel_diameter * kernel_diameter); - for (int ky = -radius; ky <= radius; ++ky) { - for (int kx = -radius; kx <= radius; ++kx) { - double dist2 = static_cast(kx * kx + ky * ky); - spatial_weights[(ky + radius) * kernel_width + (kx + radius)] = - std::exp(-dist2 / two_sigma_space_sq); + // 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); } } - // NOTE: precompute Range Weights - std::vector range_lut(MAX_RGB_DIST_SQ + 1); - double two_sigma_range_sq = 2 * sigma_range * sigma_range; + // ========= 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] = std::exp(-static_cast(i) / two_sigma_range_sq); + 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 ========= + + 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) { + rgb_to_lab(r0, g0, b0, L0, A0, B0); + } + // ========= CIELAB-only 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]; - - double r_acc = 0.0, g_acc = 0.0, b_acc = 0.0, weight_acc = 0.0; + double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}, weight_acc{0.0}; - for (int ky = -radius; ky <= radius; ++ky) { - int ny = std::clamp(y + ky, 0, h - 1); + 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); + for (int kx{-radius}; kx <= radius; ++kx) { + int nx{std::clamp(x + kx, 0, w - 1)}; - size_t neighbor_idx = (ny * width + nx) * 4; + 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]; + uint8_t r{image[neighbor_idx]}; + uint8_t g{image[neighbor_idx + 1]}; + uint8_t b{image[neighbor_idx + 2]}; - double w_space = spatial_weights[(ky + radius) * kernel_width + (kx + radius)]; + double w_space{spatial_weights[(ky + radius) * kernel_diameter + (kx + radius)]}; - int dr = static_cast(r) - r0; - int dg = static_cast(g) - g0; - int db = static_cast(b) - b0; - int dist_sq = dr*dr + dg*dg + db*db; + double w_range; + 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]; + break; + } + case COLOR_SPACE_OPTION_CIELAB: { + double L, A, B; + rgb_to_lab(r, g, b, L, A, B); + const double dist = std::sqrt((L-L0)*(L-L0) + (A-A0)*(A-A0) + (B-B0)*(B-B0)); + w_range = gaussian(dist, sigma_range); + break; + } + } - double w_range = range_lut[dist_sq]; - double w = w_space * w_range; + double w{w_space * w_range}; r_acc += r * w; g_acc += g * w; @@ -103,7 +136,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, result[center_idx + 2] = static_cast(std::clamp(b_acc / weight_acc, 0.0, 255.0)); result[center_idx + 3] = a0; } - } + } std::memcpy(image, result.data(), result.size()); } @@ -112,6 +145,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, // Global wrapper for WASM export EXPORTED void bilateral_filter(uint8_t *image, size_t width, size_t height, - double sigma_spatial, double sigma_range) { - bilateral::bilateral_filter(image, width, height, sigma_spatial, sigma_range); + 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/image_utils.cpp b/src/wasm/modules/image/src/image_utils.cpp index b3c5d217b..9e06b8f72 100644 --- a/src/wasm/modules/image/src/image_utils.cpp +++ b/src/wasm/modules/image/src/image_utils.cpp @@ -1,6 +1,5 @@ #include "image_utils.h" #include "fft_iterative.h" -#include "cielab.h" #include #include @@ -8,102 +7,6 @@ #include #include -double evaluate_gaussian(float x, double sigma) { - // evaluates 1d gaussian function desribed by sigma at x - return exp(-(pow(x, 2))/(2 * pow(sigma, 2))) / (2 * M_PI * pow(sigma, 2)); -} - -void bilateral_filter_cielab(uint8_t *image, size_t width, size_t height, double sigma_pixels, double sigma_range) -{ - // sigma_pixel = spatial kernel - if (!image || width == 0 || height == 0 || sigma_pixels <= 0 || sigma_range <= 0) - return; - - const int radius = static_cast(1.5 * sigma_pixels); - const size_t diameter = radius * 2 + 1; - - // precompute - double spatial_filter[diameter * diameter]; - for (int i = 0; i < diameter; i++){ - for (int j = 0; j < diameter; j++){ - float dist = static_cast(sqrt(pow(i - radius, 2) + pow(j - radius, 2))); - spatial_filter[i*diameter + j] = evaluate_gaussian(dist, sigma_pixels); - } - } - - uint8_t result[4 * height * width]; - - for (int i = 0; i < height; i++){ - for (int j = 0; j < width; j++) { - int center_index = 4 * (i * width + j); - uint8_t r0 = image[center_index]; - uint8_t g0 = image[center_index + 1]; - uint8_t b0 = image[center_index + 2]; - uint8_t a0 = image[center_index + 3]; - - double L0, A0, B0; - rgb_to_lab(r0, g0, b0, L0, A0, B0); - - double rf = 0.0; - double gf = 0.0; - double bf = 0.0; - double rW = 0.0; - double gW = 0.0; - double bW = 0.0; - - for (int ki = -radius; ki <= radius; ki++){ - for (int kj = -radius; kj <= radius; kj++){ - int _i = i + ki; - int _j = j + kj; - if (_i < 0) - _i = 0; - if (_i > height - 1) - _i = height - 1; - if (_j < 0) - _j = 0; - if (_j > width - 1) - _j = width - 1; - int index = 4 * (_i * width + _j); - - uint8_t r = image[index]; - uint8_t g = image[index + 1]; - uint8_t b = image[index + 2]; - - /* - as described in https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Tomasi98.pdf - use euclidean distance in LAB color space for less artifacts - */ - double L, A, B; - rgb_to_lab(r, g, b, L, A, B); - float dist = sqrt(pow(static_cast(L-L0), 2) + pow(static_cast(A-A0), 2) + pow(static_cast(B-B0), 2)); - - double w_euc = evaluate_gaussian(dist, sigma_range) * spatial_filter[ (ki + radius) * diameter + (kj + radius) ]; - double wr = w_euc; - double wg = w_euc; - double wb = w_euc; - - rf += r * wr; - rW += wr; - - gf += g * wg; - gW += wg; - - bf += b * wb; - bW += wb; - } - } - - result[center_index] = static_cast(rf / rW); - result[center_index + 1] = static_cast(gf / gW); - result[center_index + 2] = static_cast(bf / bW); - result[center_index + 3] = a0; - } - } - - image = result; -} - - // image: pointer to RGBA data // width, height: dimensions // sigma: standard deviation of Gaussian blur From e3299c56eda290778cb74299695b1cdd7e87614c Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 22:06:18 +0200 Subject: [PATCH 29/53] refactor(cielab.h): split into .h & .cpp, include guard --- src/wasm/modules/image/include/cielab.h | 62 ++----------------- src/wasm/modules/image/src/cielab.cpp | 81 +++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 56 deletions(-) create mode 100644 src/wasm/modules/image/src/cielab.cpp diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index fe9af6425..320be1621 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -1,60 +1,10 @@ +#ifndef CIELAB_H +#define CIELAB_H + #include #include #include -// Function for the non-linear XYZ to Lab transformation -double f_xyz(double t) { - if (t > 0.008856) { - return std::pow(t, 1.0/3.0); - } else { - return (7.787 * t) + (16.0 / 116.0); - } -} - -// Function for the non-linear sRGB to linear RGB transformation (inverse gamma correction) -double inverse_gamma(double c) { - if (c > 0.04045) { - return std::pow((c + 0.055) / 1.055, 2.4); - } else { - return c / 12.92; - } -} - -void rgb_to_lab(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B) { - // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] - double r = r_u8 / 255.0; - double g = g_u8 / 255.0; - double b = b_u8 / 255.0; - - r = inverse_gamma(r); - g = inverse_gamma(g); - b = inverse_gamma(b); - - // 2. Convert linear RGB to CIE XYZ (using D65 white point reference) - // The matrix below is for sRGB to XYZ (D65) - double x = 0.4124564 * r + 0.3575761 * g + 0.1804375 * b; - double y = 0.2126729 * r + 0.7151522 * g + 0.0721750 * b; - double z = 0.0193339 * r + 0.1191920 * g + 0.9503041 * b; - - // Reference white point for D65 illuminant - const double Xn = 0.95047; - const double Yn = 1.00000; - const double Zn = 1.08883; - - // Normalize XYZ values by the white point - double Xr = x / Xn; - double Yr = y / Yn; - double Zr = z / Zn; - - // 3. Convert CIE XYZ to CIE L*a*b* - double fx = f_xyz(Xr); - double fy = f_xyz(Yr); - double fz = f_xyz(Zr); - - L = 116.0 * fy - 16.0; - A = 500.0 * (fx - fy); - B = 200.0 * (fy - fz); - - // Clamp L channel to standard range [0, 100] - L = std::max(0.0, std::min(100.0, L)); -} \ No newline at end of file +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); +#endif // CIELAB_H diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp new file mode 100644 index 000000000..14c350863 --- /dev/null +++ b/src/wasm/modules/image/src/cielab.cpp @@ -0,0 +1,81 @@ +#include "cielab.h" + +// ====== 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 + +// ====== Used in rgb_to_lab ====== +// Multipliers for RGB to XYZ +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}; + +// 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); +} From a482a16e354a6f0581575a0df54232a6f0f61574 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Sat, 3 Jan 2026 23:01:38 +0200 Subject: [PATCH 30/53] fix(bilateral filter): guard against unknown color_space param Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> --- src/wasm/modules/image/src/bilateral_filter.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index bcd0e537b..263bb3eb5 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -41,6 +41,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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)}; From 4425cce777611b737d7f54b72823de762d3deeb0 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Sat, 3 Jan 2026 23:04:18 +0200 Subject: [PATCH 31/53] fix(cielab.{h,cpp} includes): properly structured --- src/wasm/modules/image/include/cielab.h | 4 +--- src/wasm/modules/image/src/cielab.cpp | 4 ++++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 320be1621..15350359b 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -1,9 +1,7 @@ #ifndef CIELAB_H #define CIELAB_H -#include -#include -#include +#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); diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 14c350863..eb00bd9ed 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -1,5 +1,9 @@ #include "cielab.h" +#include +#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 From 0ee78a6a9bf99d7685a52f457455f80e71350d41 Mon Sep 17 00:00:00 2001 From: Krasner Date: Sat, 3 Jan 2026 22:39:02 +0000 Subject: [PATCH 32/53] debugging memory and run time... --- src/hooks/useWasmWorker.js | 2 +- src/wasm/modules/image/CMakeLists.txt | 2 +- src/wasm/modules/image/include/cielab.h | 1 + .../modules/image/src/bilateral_filter.cpp | 29 +++++---- src/wasm/modules/image/src/cielab.cpp | 61 +++++++++++++++++++ 5 files changed, 78 insertions(+), 17 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index f1c502d22..acc90cbd7 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,7 +35,7 @@ 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 }) => { + const bilateralFilter = async ({ pixels, width, height, sigma_spatial = width * 0.005, 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 }) => { diff --git a/src/wasm/modules/image/CMakeLists.txt b/src/wasm/modules/image/CMakeLists.txt index 3391f7aa5..ef66eb4bf 100644 --- a/src/wasm/modules/image/CMakeLists.txt +++ b/src/wasm/modules/image/CMakeLists.txt @@ -49,7 +49,7 @@ target_link_options(${MODULE_NAME}_wasm PRIVATE ${COMMON_FLAGS}) # Build-type specific flags if(CMAKE_BUILD_TYPE STREQUAL "Debug") - target_compile_options(${MODULE_NAME}_wasm PRIVATE -O0 -g4) + target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -g4) target_link_options(${MODULE_NAME}_wasm PRIVATE "SHELL:-s ASSERTIONS=2" -g4 diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 15350359b..ed94f7457 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -5,4 +5,5 @@ 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 rgb_to_lab2(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B); #endif // CIELAB_H diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 263bb3eb5..d186e3f51 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -7,6 +7,7 @@ #include #include #include +#include namespace bilateral { @@ -18,7 +19,7 @@ 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) { +double gaussian(double x, double sigma) { return std::exp(-(x * x) / (2.0 * sigma * sigma)); } @@ -63,7 +64,6 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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); } @@ -84,11 +84,14 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, // ========= CIELAB-only section start ========= double L0, A0, B0; if (color_space == COLOR_SPACE_OPTION_CIELAB) { - rgb_to_lab(r0, g0, b0, L0, A0, B0); + rgb_to_lab2(r0, g0, b0, L0, A0, B0); + // std::cout << "cielab0 done" << std::endl; } // ========= CIELAB-only section end ========= double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}, weight_acc{0.0}; + double w_space, w_range; + double L, A, B, dist; for (int ky{-radius}; ky <= radius; ++ky) { int ny{std::clamp(y + ky, 0, h - 1)}; @@ -102,9 +105,8 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, uint8_t g{image[neighbor_idx + 1]}; uint8_t b{image[neighbor_idx + 2]}; - double w_space{spatial_weights[(ky + radius) * kernel_diameter + (kx + radius)]}; - - double w_range; + 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; @@ -115,20 +117,17 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, break; } case COLOR_SPACE_OPTION_CIELAB: { - double L, A, B; - rgb_to_lab(r, g, b, L, A, B); - const double dist = std::sqrt((L-L0)*(L-L0) + (A-A0)*(A-A0) + (B-B0)*(B-B0)); + rgb_to_lab2(r, g, b, L, A, B); + dist = std::sqrt((L-L0)*(L-L0) + (A-A0)*(A-A0) + (B-B0)*(B-B0)); w_range = gaussian(dist, sigma_range); break; } } - double w{w_space * w_range}; - - r_acc += r * w; - g_acc += g * w; - b_acc += b * w; - weight_acc += w; + r_acc += r * w_space * w_range; + g_acc += g * w_space * w_range; + b_acc += b * w_space * w_range; + weight_acc += w_space * w_range; } } diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index eb00bd9ed..772ac0624 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -83,3 +83,64 @@ void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, out_l = std::clamp(out_l, 0.0, 100.0); } + +#include +#include +#include + +// Function for the non-linear XYZ to Lab transformation +double f_xyz(double t) { + if (t > 0.008856) { + return std::pow(t, 1.0/3.0); + } else { + return (7.787 * t) + (16.0 / 116.0); + } +} + +// Function for the non-linear sRGB to linear RGB transformation (inverse gamma correction) +double inverse_gamma(double c) { + if (c > 0.04045) { + return std::pow((c + 0.055) / 1.055, 2.4); + } else { + return c / 12.92; + } +} + +void rgb_to_lab2(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B) { + // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] + double r = r_u8 / 255.0; + double g = g_u8 / 255.0; + double b = b_u8 / 255.0; + + r = inverse_gamma(r); + g = inverse_gamma(g); + b = inverse_gamma(b); + + // 2. Convert linear RGB to CIE XYZ (using D65 white point reference) + // The matrix below is for sRGB to XYZ (D65) + double x = 0.4124564 * r + 0.3575761 * g + 0.1804375 * b; + double y = 0.2126729 * r + 0.7151522 * g + 0.0721750 * b; + double z = 0.0193339 * r + 0.1191920 * g + 0.9503041 * b; + + // Reference white point for D65 illuminant + const double Xn = 0.95047; + const double Yn = 1.00000; + const double Zn = 1.08883; + + // Normalize XYZ values by the white point + double Xr = x / Xn; + double Yr = y / Yn; + double Zr = z / Zn; + + // 3. Convert CIE XYZ to CIE L*a*b* + double fx = f_xyz(Xr); + double fy = f_xyz(Yr); + double fz = f_xyz(Zr); + + L = 116.0 * fy - 16.0; + A = 500.0 * (fx - fy); + B = 200.0 * (fy - fz); + + // Clamp L channel to standard range [0, 100] + L = std::max(0.0, std::min(100.0, L)); +} From dcb4cb76992b3ab6fee51cf203b42cf297f0b93f Mon Sep 17 00:00:00 2001 From: Krasner Date: Sun, 4 Jan 2026 15:05:41 +0000 Subject: [PATCH 33/53] Convert full RGB image to CIELAB then look up during convolution step --- src/hooks/useWasmWorker.js | 2 +- src/wasm/modules/image/CMakeLists.txt | 4 +- src/wasm/modules/image/include/cielab.h | 7 +- .../modules/image/src/bilateral_filter.cpp | 37 ++++++++-- src/wasm/modules/image/src/cielab.cpp | 67 +------------------ 5 files changed, 40 insertions(+), 77 deletions(-) diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index acc90cbd7..f1c502d22 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,7 +35,7 @@ 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 = width * 0.005, sigma_range = 50.0, color_space = 0 }) => { + 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 }) => { diff --git a/src/wasm/modules/image/CMakeLists.txt b/src/wasm/modules/image/CMakeLists.txt index ef66eb4bf..4d80a5e36 100644 --- a/src/wasm/modules/image/CMakeLists.txt +++ b/src/wasm/modules/image/CMakeLists.txt @@ -49,13 +49,13 @@ target_link_options(${MODULE_NAME}_wasm PRIVATE ${COMMON_FLAGS}) # Build-type specific flags if(CMAKE_BUILD_TYPE STREQUAL "Debug") - target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -g4) + target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -g4 -ffast-math) target_link_options(${MODULE_NAME}_wasm PRIVATE "SHELL:-s ASSERTIONS=2" -g4 ) else() - target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3) + target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -ffast-math) target_link_options(${MODULE_NAME}_wasm PRIVATE "SHELL:-s SINGLE_FILE=0" ) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index ed94f7457..9447cb86a 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -2,8 +2,11 @@ #define CIELAB_H #include +#include +#include +#include +#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 rgb_to_lab2(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B); + double& out_l, double& out_a, double& out_b); #endif // CIELAB_H diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index d186e3f51..6f275d594 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -70,6 +70,27 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, } // ========= RGB-only section end ========= + // ========= CIELAB section start ========= + // Compute full image RGB - CIELAB conversion + std::vector cie_image(width * height * 4); + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + int 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]; + 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) { @@ -84,14 +105,15 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, // ========= CIELAB-only section start ========= double L0, A0, B0; if (color_space == COLOR_SPACE_OPTION_CIELAB) { - rgb_to_lab2(r0, g0, b0, L0, A0, B0); - // std::cout << "cielab0 done" << std::endl; + L0 = cie_image[center_idx]; + A0 = cie_image[center_idx + 1]; + B0 = cie_image[center_idx + 2]; } // ========= CIELAB-only section end ========= double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}, weight_acc{0.0}; double w_space, w_range; - double L, A, B, dist; + double dL, dA, dB, dist; for (int ky{-radius}; ky <= radius; ++ky) { int ny{std::clamp(y + ky, 0, h - 1)}; @@ -117,9 +139,12 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, break; } case COLOR_SPACE_OPTION_CIELAB: { - rgb_to_lab2(r, g, b, L, A, B); - dist = std::sqrt((L-L0)*(L-L0) + (A-A0)*(A-A0) + (B-B0)*(B-B0)); - w_range = gaussian(dist, sigma_range); + 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); // this is fast break; } } diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 772ac0624..d810ff6c8 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -1,9 +1,5 @@ #include "cielab.h" -#include -#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 @@ -82,65 +78,4 @@ void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, out_b = LAB_B_FACTOR * (fy - fz); out_l = std::clamp(out_l, 0.0, 100.0); -} - -#include -#include -#include - -// Function for the non-linear XYZ to Lab transformation -double f_xyz(double t) { - if (t > 0.008856) { - return std::pow(t, 1.0/3.0); - } else { - return (7.787 * t) + (16.0 / 116.0); - } -} - -// Function for the non-linear sRGB to linear RGB transformation (inverse gamma correction) -double inverse_gamma(double c) { - if (c > 0.04045) { - return std::pow((c + 0.055) / 1.055, 2.4); - } else { - return c / 12.92; - } -} - -void rgb_to_lab2(uint8_t r_u8, uint8_t g_u8, uint8_t b_u8, double& L, double& A, double& B) { - // 1. Convert 8-bit RGB [0, 255] to linear RGB [0.0, 1.0] - double r = r_u8 / 255.0; - double g = g_u8 / 255.0; - double b = b_u8 / 255.0; - - r = inverse_gamma(r); - g = inverse_gamma(g); - b = inverse_gamma(b); - - // 2. Convert linear RGB to CIE XYZ (using D65 white point reference) - // The matrix below is for sRGB to XYZ (D65) - double x = 0.4124564 * r + 0.3575761 * g + 0.1804375 * b; - double y = 0.2126729 * r + 0.7151522 * g + 0.0721750 * b; - double z = 0.0193339 * r + 0.1191920 * g + 0.9503041 * b; - - // Reference white point for D65 illuminant - const double Xn = 0.95047; - const double Yn = 1.00000; - const double Zn = 1.08883; - - // Normalize XYZ values by the white point - double Xr = x / Xn; - double Yr = y / Yn; - double Zr = z / Zn; - - // 3. Convert CIE XYZ to CIE L*a*b* - double fx = f_xyz(Xr); - double fy = f_xyz(Yr); - double fz = f_xyz(Zr); - - L = 116.0 * fy - 16.0; - A = 500.0 * (fx - fy); - B = 200.0 * (fy - fz); - - // Clamp L channel to standard range [0, 100] - L = std::max(0.0, std::min(100.0, L)); -} +} \ No newline at end of file From 8d267d71b16eee4b4a27fb526108c865d144de42 Mon Sep 17 00:00:00 2001 From: Krasner Date: Sun, 4 Jan 2026 15:13:25 +0000 Subject: [PATCH 34/53] put in missing if --- .../modules/image/src/bilateral_filter.cpp | 34 +++++++++++-------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 6f275d594..ba099bb69 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -72,21 +72,25 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, // ========= CIELAB section start ========= // Compute full image RGB - CIELAB conversion - std::vector cie_image(width * height * 4); - for (int y = 0; y < height; y++) { - for (int x = 0; x < width; x++) { - int 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]; - 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 + 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 * 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 ========= From 60a6e30ad37d67b8ad2e2f24cd0e41143e74875d Mon Sep 17 00:00:00 2001 From: Krasner Date: Sun, 4 Jan 2026 21:39:49 +0000 Subject: [PATCH 35/53] updates --- src/wasm/modules/image/include/cielab.h | 4 ---- src/wasm/modules/image/src/bilateral_filter.cpp | 1 - src/wasm/modules/image/src/cielab.cpp | 3 +++ 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 9447cb86a..97f1bfa9e 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -2,10 +2,6 @@ #define CIELAB_H #include -#include -#include -#include -#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); diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index ba099bb69..83dd35f39 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -7,7 +7,6 @@ #include #include #include -#include namespace bilateral { diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index d810ff6c8..6e8e4b204 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -1,4 +1,7 @@ #include "cielab.h" +#include +#include +#include // ====== Used in xyz_to_lab ======= constexpr double DELTA{6.0 / 29.0}; // 0.2068966 From 583b3247730a6d3ce085942b8494e73fbd56de4f Mon Sep 17 00:00:00 2001 From: Krasner Date: Sun, 4 Jan 2026 21:56:21 +0000 Subject: [PATCH 36/53] update bilateral filter docs to include cielab color space --- .../modules/image/bilateral_filter/api.md | 4 +++- .../image/bilateral_filter/explained.md | 24 ++++++++++++++++--- .../modules/image/src/bilateral_filter.cpp | 2 +- 3 files changed, 25 insertions(+), 5 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index d5a5276f2..e053c09c8 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -13,7 +13,8 @@ Quick reference for the function implemented in the header. void bilateral_filter(uint8_t *image, size_t width, size_t height, double sigma_spatial, - double sigma_range) + double sigma_range, + uint8_t color_space) ``` ## Parameters @@ -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++) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index 07be50365..8a8351bca 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -47,7 +47,7 @@ where $ \sigma_s$ controls spatial smoothing and $\sigma_r$ controls edge sensit Our implementation uses a **naive sliding window** approach with **Look-Up Table (LUT) optimizations** to improve performance in WebAssembly. -### 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 +60,31 @@ 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. 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 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. diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 83dd35f39..f07dd2aa5 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -147,7 +147,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, dB = cie_image[neighbor_idx + 2] - B0; dist = std::sqrt(dL * dL + dA * dA + dB * dB); - w_range = gaussian(dist, sigma_range); // this is fast + w_range = gaussian(dist, sigma_range); break; } } From 56bd7928ede1a05c7ae83e07727c58d31e09bdb0 Mon Sep 17 00:00:00 2001 From: Krasner Date: Mon, 5 Jan 2026 03:33:46 +0000 Subject: [PATCH 37/53] for cielab apply bilateral filter weights on LAB components then convert to RGB --- src/wasm/modules/image/include/cielab.h | 5 +- .../modules/image/src/bilateral_filter.cpp | 63 +++++++++++++++---- src/wasm/modules/image/src/cielab.cpp | 59 +++++++++++++++++ 3 files changed, 114 insertions(+), 13 deletions(-) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 97f1bfa9e..3eadd57df 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -4,5 +4,8 @@ #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); + 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/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index f07dd2aa5..32c991d27 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -46,6 +46,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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); @@ -114,7 +115,13 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, } // ========= CIELAB-only section end ========= - double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}, weight_acc{0.0}; + // double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}; + + // in RGB mode represents r,g,b accumulators + // in CIELAB mode represents L,A,B accumulators + double acc0{0.0}, acc1{0.0}, acc2{0.0}; + + double weight_acc{0.0}; double w_space, w_range; double dL, dA, dB, dist; @@ -130,6 +137,10 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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) { @@ -142,9 +153,9 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, break; } case COLOR_SPACE_OPTION_CIELAB: { - dL = cie_image[neighbor_idx] - L0; - dA = cie_image[neighbor_idx + 1] - A0; - dB = cie_image[neighbor_idx + 2] - B0; + dL = L - L0; + dA = A - A0; + dB = B - B0; dist = std::sqrt(dL * dL + dA * dA + dB * dB); w_range = gaussian(dist, sigma_range); @@ -152,20 +163,48 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, } } - r_acc += r * w_space * w_range; - g_acc += g * w_space * w_range; - b_acc += b * w_space * w_range; + switch (color_space) { + case COLOR_SPACE_OPTION_RGB: { + acc0 += r * w_space * w_range; + acc1 += g * w_space * w_range; + acc2 += b * w_space * w_range; + break; + } + case COLOR_SPACE_OPTION_CIELAB: { + acc0 += L * w_space * w_range; + acc1 += A * w_space * w_range; + acc2 += B * w_space * w_range; + break; + } + } weight_acc += w_space * w_range; } } - result[center_idx] = static_cast(std::clamp(r_acc / weight_acc, 0.0, 255.0)); - result[center_idx + 1] = static_cast(std::clamp(g_acc / weight_acc, 0.0, 255.0)); - result[center_idx + 2] = static_cast(std::clamp(b_acc / weight_acc, 0.0, 255.0)); - result[center_idx + 3] = a0; + switch (color_space) { + case COLOR_SPACE_OPTION_RGB: { + result[center_idx] = static_cast(std::clamp(acc0 / weight_acc, 0.0, 255.0)); + result[center_idx + 1] = static_cast(std::clamp(acc1 / weight_acc, 0.0, 255.0)); + result[center_idx + 2] = static_cast(std::clamp(acc2 / weight_acc, 0.0, 255.0)); + result[center_idx + 3] = a0; + break; + } + case COLOR_SPACE_OPTION_CIELAB: { + double L = acc0 / weight_acc; + double A = acc1 / weight_acc; + double B = acc2 / 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()); } diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 6e8e4b204..9c40f749b 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -14,6 +14,7 @@ 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 / 2.4}; // gamma exponent for nonlinear segment // ====== Used in rgb_to_lab ====== // Multipliers for RGB to XYZ @@ -27,6 +28,16 @@ constexpr double SRGB_R_TO_Z{0.0193339}; constexpr double SRGB_G_TO_Z{0.1191920}; constexpr double SRGB_B_TO_Z{0.9503041}; +constexpr double SRGB_X_TO_R{3.2406}; +constexpr double SRGB_Y_TO_R{-1.5372}; +constexpr double SRGB_Z_TO_R{-0.4986}; +constexpr double SRGB_X_TO_G{-0.9689}; +constexpr double SRGB_Y_TO_G{1.8758}; +constexpr double SRGB_Z_TO_G{0.0415}; +constexpr double SRGB_X_TO_B{0.0557}; +constexpr double SRGB_Y_TO_B{0.2040}; +constexpr double SRGB_Z_TO_B{1.0570}; + // Reference white point for D65 illuminant constexpr double D65_Xn = 0.95047; constexpr double D65_Yn = 1.0; @@ -81,4 +92,52 @@ void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, out_b = LAB_B_FACTOR * (fy - fz); out_l = std::clamp(out_l, 0.0, 100.0); +} + +inline double finv(double t) +{ + if (t > DELTA) + return t * t * t; + else + return 3 * DELTA * DELTA * (t - EPSILON); +} + +inline double clamp01(double v) +{ + return std::min(1.0, std::max(0.0, v)); +} + +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) +{ + // --- 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; + + double X = D65_Xn * finv(fx); + double Y = D65_Yn * finv(fy); + double Z = D65_Zn * finv(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) + auto gamma_encode = [](double u) -> double { + if (u <= SRGB_LINEAR_THRESHOLD / SRGB_LINEAR_FACTOR) + return SRGB_LINEAR_FACTOR * u; + else + return (1.0 + SRGB_GAMMA_OFFSET) * std::pow(u, SRGB_GAMMA_INV) - SRGB_GAMMA_OFFSET; + }; + + r = gamma_encode(r); + g = gamma_encode(g); + b = gamma_encode(b); + + // --- Clamp and convert to 8-bit + r_u8 = static_cast(std::round(255.0 * clamp01(r))); + g_u8 = static_cast(std::round(255.0 * clamp01(g))); + b_u8 = static_cast(std::round(255.0 * clamp01(b))); } \ No newline at end of file From 1d7a5a45fdfa1802ef382145aa3fc0bc07c17836 Mon Sep 17 00:00:00 2001 From: Krasner Date: Mon, 5 Jan 2026 04:23:50 +0000 Subject: [PATCH 38/53] fix bug --- src/wasm/modules/image/src/cielab.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 9c40f749b..144932ad6 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -35,7 +35,7 @@ constexpr double SRGB_X_TO_G{-0.9689}; constexpr double SRGB_Y_TO_G{1.8758}; constexpr double SRGB_Z_TO_G{0.0415}; constexpr double SRGB_X_TO_B{0.0557}; -constexpr double SRGB_Y_TO_B{0.2040}; +constexpr double SRGB_Y_TO_B{-0.2040}; constexpr double SRGB_Z_TO_B{1.0570}; // Reference white point for D65 illuminant From 9a7580a00046d39e7568055af2868724ee745d41 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Mon, 5 Jan 2026 22:07:40 +0200 Subject: [PATCH 39/53] docs(bilateral filter): update for better clarity --- .../modules/image/bilateral_filter/api.md | 14 +- .../image/bilateral_filter/explained.md | 43 ++- .../image/bilateral_filter/implementation.md | 250 ++++++++++++++++-- .../image/bilateral_filter/keywords.md | 28 ++ .../image/bilateral_filter/overview.md | 26 +- 5 files changed, 304 insertions(+), 57 deletions(-) create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index e053c09c8..39317d478 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -2,21 +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, - uint8_t color_space) +```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 | diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index 8a8351bca..49ae79358 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/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,7 +58,7 @@ 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 (RGB color space) @@ -62,7 +75,9 @@ for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) { ### 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. Instead range weights are computed on the fly using the `gaussian` function. +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. @@ -76,7 +91,9 @@ dist = std::sqrt(dL * dL + dA * dA + dB * dB); w_range = gaussian(dist, sigma_range); ``` -NOTE: `gaussian` itself is expensive to run. Future optimizations include polynomial approximations of `exp(-x^2)` via Taylor expansion or Horner's method. +:::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 diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md index ab5c6b9fb..645fff30e 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md @@ -2,14 +2,18 @@ id: implementation title: Bilateral Filter — Implementation details sidebar_label: Implementation -sidebar_position: 4 +sidebar_position: 2 --- # Bilateral Filter — Implementation details -This page maps the conceptual steps of the Bilateral Filter to the concrete functions and loops in the implementation. +This page maps the conceptual steps of the Bilateral Filter to the concrete implementation. -## 1. Parameters & Window Size +:::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}$). @@ -18,53 +22,241 @@ const int radius = static_cast(std::ceil(SIGMA_RADIUS_FACTOR * sigma_spatia 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). +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 */} + + +
-## 2. Precomputing Weights (Optimization) +:::important Kernel Dimensions +The filter kernel is always square; width = height = 2 * radius + 1. +::: -To avoid computing `std::exp` millions of times per frame, we precalculate the weights. +## 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_width + (kx + radius)] = - std::exp(-dist2 / two_sigma_space_sq); +spatial_weights[(ky + radius) * kernel_diameter + (kx + radius)] = gaussian(dist, sigma_spatial); ``` -### 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. +:::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 -for (int i = 0; i <= MAX_RGB_DIST_SQ; ++i) { - range_lut[i] = std::exp(-static_cast(i) / two_sigma_range_sq); -} +range_lut[i] = gaussian(static_cast(std::sqrt(i)), sigma_range); ``` -## 3. Sliding Window Loop +:::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. -The core processing happens in a nested loop over every pixel $(y, x)$. For each pixel, we: +See the corresponding information block for CIELAB to see why this differs between the color spaces. +::: + +#### "On the Fly" CIELAB Range Weights -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. +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 -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; +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. + +
+ + {/* 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, $\sigma_{r} = 20$ for RGB and $\sigma_{r} = 15$ for CIELAB. +
+ +##### 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 to get the filtered pixel value: +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)); 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..844ba0c65 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md @@ -0,0 +1,28 @@ +--- +id: keywords +title: Keywords +sidebar_position: 6 +--- + +- **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 index 91705a4c7..3dfd6731b 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md @@ -2,7 +2,7 @@ id: overview title: Bilateral Filter sidebar_label: Overview -sidebar_position: 2 +sidebar_position: 1 --- # Bilateral Filter @@ -13,19 +13,29 @@ This section introduces the **bilateral filter** used in the Img2Num project 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). -- **Data type:** `uint8_t` (8-bit unsigned integer channels). +- **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 difference** (Gaussian). + 3. Weight neighbors by **intensity/color difference** (Gaussian). 4. Normalize and average. -## Pages in this mini-guide +## Keywords -* **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. +The [keywords section](../keywords/) will help you in case the terminology confuses you. -Jump to implementation: [Implementation details](../implementation/) +:::tip Gaussian functions +Bilateral filters rely on Gaussian weighting, so understanding Gaussian functions will help when reading the implementation. +::: From 85fd528d39e2871aa32fd4c2a3165ef9f3a0c02e Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Mon, 5 Jan 2026 22:37:09 +0200 Subject: [PATCH 40/53] refactor(cielab): use constexpr functions and improve naming --- src/wasm/modules/image/include/cielab.h | 6 +- src/wasm/modules/image/src/cielab.cpp | 74 +++++++++++-------------- 2 files changed, 35 insertions(+), 45 deletions(-) diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 3eadd57df..2dc7e9616 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -4,8 +4,8 @@ #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); + 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); +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/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 144932ad6..3820b340b 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -1,7 +1,6 @@ #include "cielab.h" #include #include -#include // ====== Used in xyz_to_lab ======= constexpr double DELTA{6.0 / 29.0}; // 0.2068966 @@ -10,11 +9,11 @@ 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 / 2.4}; // gamma exponent for nonlinear segment +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 ====== // Multipliers for RGB to XYZ @@ -39,26 +38,26 @@ constexpr double SRGB_Y_TO_B{-0.2040}; constexpr double SRGB_Z_TO_B{1.0570}; // 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 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; +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; + 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 + 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); } @@ -94,30 +93,28 @@ void rgb_to_lab(const uint8_t r_u8, const uint8_t g_u8, const uint8_t b_u8, out_l = std::clamp(out_l, 0.0, 100.0); } -inline double finv(double t) +constexpr double inverse_xyz_to_lab(double t) { - if (t > DELTA) - return t * t * t; - else - return 3 * DELTA * DELTA * (t - EPSILON); + return (t > DELTA) ? (t * t * t) : (3 * DELTA * DELTA * (t - EPSILON)); } -inline double clamp01(double v) -{ - return std::min(1.0, std::max(0.0, v)); -} +constexpr double gamma_encode(double 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& r_u8, uint8_t& g_u8, uint8_t& b_u8) + 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 fy{(L + LAB_L_OFFSET) / LAB_L_FACTOR}; + const double fx{fy + A / LAB_A_FACTOR}; + const double fz{fy - B / LAB_B_FACTOR}; - double X = D65_Xn * finv(fx); - double Y = D65_Yn * finv(fy); - double Z = D65_Zn * finv(fz); + 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}; @@ -125,19 +122,12 @@ void lab_to_rgb(const double L, const double A, const double B, double b{SRGB_X_TO_B * X + SRGB_Y_TO_B * Y + SRGB_Z_TO_B * Z}; // --- linear RGB → sRGB (gamma correction) - auto gamma_encode = [](double u) -> double { - if (u <= SRGB_LINEAR_THRESHOLD / SRGB_LINEAR_FACTOR) - return SRGB_LINEAR_FACTOR * u; - else - return (1.0 + SRGB_GAMMA_OFFSET) * std::pow(u, SRGB_GAMMA_INV) - SRGB_GAMMA_OFFSET; - }; - r = gamma_encode(r); g = gamma_encode(g); b = gamma_encode(b); // --- Clamp and convert to 8-bit - r_u8 = static_cast(std::round(255.0 * clamp01(r))); - g_u8 = static_cast(std::round(255.0 * clamp01(g))); - b_u8 = static_cast(std::round(255.0 * clamp01(b))); -} \ No newline at end of file + 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))); +} From ef96ab47231d8723c76d3d16cfec7e8b583fdf0b Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Mon, 5 Jan 2026 22:47:59 +0200 Subject: [PATCH 41/53] refactor(bilateral_filter): improve readability and use brace initialization - Replace acc0/acc1/acc2 with descriptive weight_acc_channel_0/1/2 - Use brace-initialization for ints and doubles - Minor spacing and formatting cleanup --- .../modules/image/src/bilateral_filter.cpp | 58 +++++++++---------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 32c991d27..621ada57b 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -46,7 +46,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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); @@ -75,14 +75,14 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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 * 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]; + + 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); @@ -115,11 +115,9 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, } // ========= CIELAB-only section end ========= - // double r_acc{0.0}, g_acc{0.0}, b_acc{0.0}; - // in RGB mode represents r,g,b accumulators // in CIELAB mode represents L,A,B accumulators - double acc0{0.0}, acc1{0.0}, acc2{0.0}; + 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; @@ -142,13 +140,13 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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; + 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]; break; } @@ -165,15 +163,15 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, switch (color_space) { case COLOR_SPACE_OPTION_RGB: { - acc0 += r * w_space * w_range; - acc1 += g * w_space * w_range; - acc2 += b * w_space * w_range; + 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: { - acc0 += L * w_space * w_range; - acc1 += A * w_space * w_range; - acc2 += B * w_space * w_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; } } @@ -183,16 +181,16 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, switch (color_space) { case COLOR_SPACE_OPTION_RGB: { - result[center_idx] = static_cast(std::clamp(acc0 / weight_acc, 0.0, 255.0)); - result[center_idx + 1] = static_cast(std::clamp(acc1 / weight_acc, 0.0, 255.0)); - result[center_idx + 2] = static_cast(std::clamp(acc2 / weight_acc, 0.0, 255.0)); + 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 = acc0 / weight_acc; - double A = acc1 / weight_acc; - double B = acc2 / weight_acc; + 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; @@ -204,7 +202,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, } } } - + std::memcpy(image, result.data(), result.size()); } From fc770b00e1317754f10f3a69e992b382938e8bf3 Mon Sep 17 00:00:00 2001 From: Krasner Date: Mon, 5 Jan 2026 20:59:37 +0000 Subject: [PATCH 42/53] condense 2 switch cases --- .../modules/image/src/bilateral_filter.cpp | 24 +++++++------------ 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 621ada57b..45bdeb23c 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -148,6 +148,10 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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: { @@ -157,24 +161,14 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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; } } - - switch (color_space) { - case COLOR_SPACE_OPTION_RGB: { - 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: { - 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; } } From 65d4b542a1bf5d9fcdc86e37d05c4a58259729ef Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Mon, 5 Jan 2026 23:30:23 +0200 Subject: [PATCH 43/53] refactor(bilateral filter): inline gaussian function --- src/wasm/modules/image/src/bilateral_filter.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wasm/modules/image/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 45bdeb23c..71faf9c51 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -18,7 +18,7 @@ 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}; -double gaussian(double x, double sigma) { +inline double gaussian(double x, double sigma) { return std::exp(-(x * x) / (2.0 * sigma * sigma)); } @@ -168,7 +168,7 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, break; } } - + weight_acc += w_space * w_range; } } From 16c00974e26401641749029bfe851447ea6e2987 Mon Sep 17 00:00:00 2001 From: Ryan Millard <142347829+Ryan-Millard@users.noreply.github.com> Date: Tue, 6 Jan 2026 00:07:25 +0200 Subject: [PATCH 44/53] fix(cpp: cielab): gaussian function inlined now Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> --- src/wasm/modules/image/src/cielab.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index 3820b340b..ca8da11e2 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -98,11 +98,11 @@ constexpr double inverse_xyz_to_lab(double t) return (t > DELTA) ? (t * t * t) : (3 * DELTA * DELTA * (t - EPSILON)); } -constexpr double gamma_encode(double u) { +inline double gamma_encode(double 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) From 6e4cd68866fed736ccf33ffa7492c2890315da02 Mon Sep 17 00:00:00 2001 From: Krasner Date: Tue, 6 Jan 2026 02:21:24 +0000 Subject: [PATCH 45/53] RGB-CIELAB conversion write up --- .../wasm/modules/image/cielab/explained.md | 183 ++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100644 docs/docs/reference/wasm/modules/image/cielab/explained.md diff --git a/docs/docs/reference/wasm/modules/image/cielab/explained.md b/docs/docs/reference/wasm/modules/image/cielab/explained.md new file mode 100644 index 000000000..21e5e5b1a --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/cielab/explained.md @@ -0,0 +1,183 @@ +--- +id: explained +title: Implementation Explained +sidebar_position: 5 +--- + + +# RGB ↔ CIELAB Conversion Guide + +This explains the full mathematical conversion pipeline between **sRGB** and **CIELAB (Lab)** color spaces. + +--- + +# 1. Conversion Pipeline Overview + +## RGB → CIELAB +1. sRGB → Linear RGB +2. Linear RGB → XYZ +3. XYZ → CIELAB + +## CIELAB → RGB +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.2406 & -1.5372 & -0.4986 \\ +-0.9689 & 1.8758 & 0.0415 \\ + 0.0557 & -0.2040 & 1.0570 +\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 From 7baacc886a43d30131caf03584afcde3727834f6 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 05:13:45 +0200 Subject: [PATCH 46/53] docs(bilateral filter): add color space docs --- .../image/bilateral_filter/color-spaces.md | 219 ++++++++++++++++++ .../image/bilateral_filter/implementation.md | 55 +---- .../image/bilateral_filter/keywords.md | 1 - .../wasm/modules/image/cielab/_category_.json | 10 + .../wasm/modules/image/cielab/api.md | 218 +++++++++++++++++ .../image/cielab/{explained.md => index.md} | 3 +- .../bilateral_filter/RgbVsLabRangeKernel.jsx | 60 +++++ 7 files changed, 510 insertions(+), 56 deletions(-) create mode 100644 docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md create mode 100644 docs/docs/reference/wasm/modules/image/cielab/_category_.json create mode 100644 docs/docs/reference/wasm/modules/image/cielab/api.md rename docs/docs/reference/wasm/modules/image/cielab/{explained.md => index.md} (99%) create mode 100644 docs/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel.jsx 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..778a1fd22 --- /dev/null +++ b/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'; + + + +**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 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/implementation.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md index 645fff30e..428731ad3 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md @@ -164,60 +164,9 @@ A coarse LUT could lead to visible artifacts. Bilateral filtering computes a **range weight** for each pixel in the kernel based on the color difference between the center pixel and its neighbor. -
- - {/* 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) +import RgbVsLabRangeKernel from '@site/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel'; - {/* Axes */} - - - Color difference ΔRGB - Color difference ΔLAB - - - Above, $\sigma_{r} = 20$ for RGB and $\sigma_{r} = 15$ for CIELAB. -
+ ##### RGB LUT (Red curve and shaded 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 index 844ba0c65..44a0aeea9 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md @@ -1,7 +1,6 @@ --- id: keywords title: Keywords -sidebar_position: 6 --- - **spatial component**: The part of an image related to the **pixel positions** (x and y coordinates). 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..399432a17 --- /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: 1 +--- + +# 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/explained.md b/docs/docs/reference/wasm/modules/image/cielab/index.md similarity index 99% rename from docs/docs/reference/wasm/modules/image/cielab/explained.md rename to docs/docs/reference/wasm/modules/image/cielab/index.md index 21e5e5b1a..5d244753c 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/explained.md +++ b/docs/docs/reference/wasm/modules/image/cielab/index.md @@ -1,10 +1,9 @@ --- id: explained title: Implementation Explained -sidebar_position: 5 +sidebar_position: 2 --- - # RGB ↔ CIELAB Conversion Guide This explains the full mathematical conversion pipeline between **sRGB** and **CIELAB (Lab)** color spaces. 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..f95444357 --- /dev/null +++ b/docs/src/components/docs/reference/wasm/modules/image/bilateral_filter/RgbVsLabRangeKernel.jsx @@ -0,0 +1,60 @@ +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; From fd3cff06b41ce6e81d7d02781e7c66d1b8de0950 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 05:22:55 +0200 Subject: [PATCH 47/53] style(formatting): fix formatting issues on all files --- .../image/bilateral_filter/_category_.json | 2 +- .../modules/image/bilateral_filter/api.md | 19 +- .../image/bilateral_filter/color-spaces.md | 55 ++- .../image/bilateral_filter/explained.md | 5 + .../image/bilateral_filter/implementation.md | 25 +- .../image/bilateral_filter/keywords.md | 14 +- .../image/bilateral_filter/overview.md | 5 +- .../wasm/modules/image/cielab/api.md | 54 +-- .../wasm/modules/image/cielab/index.md | 45 ++- .../bilateral_filter/RgbVsLabRangeKernel.jsx | 44 ++- src/hooks/useWasmWorker.js | 13 +- .../modules/image/include/bilateral_filter.h | 6 +- src/wasm/modules/image/include/cielab.h | 6 +- .../modules/image/src/bilateral_filter.cpp | 340 +++++++++--------- src/wasm/modules/image/src/cielab.cpp | 150 ++++---- src/wasm/modules/image/src/kmeans.cpp | 13 +- 16 files changed, 443 insertions(+), 353 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json index dea27276f..2ad8ea622 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/_category_.json @@ -7,4 +7,4 @@ "description": "Documentation for the Bilateral Filter in the Image WebAssembly (WASM) module in Img2Num.", "slug": "/reference/wasm/modules/image/bilateral_filter" } -} \ No newline at end of file +} diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index 39317d478..c85a7fd4c 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -19,16 +19,17 @@ The alpha channel, `image[i + 3]`, is left untouched - it is not part of the bil ## 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. | +| 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`. -::: + ::: 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 index 778a1fd22..757b9ac8a 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md @@ -11,22 +11,24 @@ The bilateral filter in Img2Num supports two color spaces for computing range (c ## 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 | +| 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 @@ -39,6 +41,7 @@ The bilateral filter in Img2Num supports two color spaces for computing range (c 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: $$ @@ -46,12 +49,15 @@ $$ $$ 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) @@ -61,15 +67,17 @@ $$ $$ 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 @@ -93,6 +101,7 @@ import RgbVsLabRangeKernel from '@site/src/components/docs/reference/wasm/module **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** @@ -106,38 +115,41 @@ 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; + // Scale RGB sigma_range to match CIELAB perceptually + sigma_range_actual = sigma_range_base * 4.18; } else { - sigma_range_actual = sigma_range_base; + 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` | +| 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 @@ -146,14 +158,16 @@ Using the same `sigma_range = 50`: :::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); @@ -170,6 +184,7 @@ 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); @@ -189,13 +204,17 @@ 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 @@ -204,10 +223,12 @@ Consider CIELAB when you observe: ### 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 diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index 49ae79358..0f9475980 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -11,6 +11,7 @@ This section explains the inner workings of the **bilateral filter** implementat ## 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. @@ -19,6 +20,7 @@ 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)|) $$ @@ -56,6 +58,7 @@ $$ 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. @@ -63,6 +66,7 @@ Our implementation uses a **naive sliding window** approach with **Look-Up Table ### 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. @@ -82,6 +86,7 @@ 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; diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md index 428731ad3..de1a5c523 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md @@ -73,11 +73,12 @@ We primarily use $\sigma_{spatial} \approx 3.0$, which results in a kernel radiu {/* center pixel */} + :::important Kernel Dimensions -The filter kernel is always square; width = height = 2 * radius + 1. +The filter kernel is always square; width = height = 2 \* radius + 1. ::: ## 2. Computing Weights @@ -98,6 +99,7 @@ double gaussian(double x, double sigma) { return std::exp(-(x * x) / (2.0 * sigma * sigma)); } ``` + ::: ### Spatial Weights (constant per kernel) @@ -143,19 +145,21 @@ 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. + 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. + 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. -::: + A coarse LUT could lead to visible artifacts. + :::
@@ -178,6 +182,7 @@ 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. @@ -197,10 +202,10 @@ The core processing happens in a nested loop over every pixel $(y, x)$. For each 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". + - **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 diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md index 44a0aeea9..a57010075 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/keywords.md @@ -9,19 +9,19 @@ title: Keywords - **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. +- **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. +- **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**. +- **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. +- **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 index 3dfd6731b..0bb808967 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/overview.md @@ -17,12 +17,13 @@ and where the corresponding code lives so you can jump straight into the impleme 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* - +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/)). diff --git a/docs/docs/reference/wasm/modules/image/cielab/api.md b/docs/docs/reference/wasm/modules/image/cielab/api.md index 399432a17..39970c829 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/api.md +++ b/docs/docs/reference/wasm/modules/image/cielab/api.md @@ -30,14 +30,14 @@ void rgb_to_lab( #### 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) | +| 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 @@ -84,14 +84,14 @@ void lab_to_rgb( #### 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) | +| 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 @@ -130,13 +130,13 @@ lab_to_rgb(L, a, b, r, g, b_rgb); ### 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 | +| 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 @@ -144,12 +144,10 @@ lab_to_rgb(L, a, b, r, g, b_rgb); - 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 @@ -164,6 +162,7 @@ $$ $$ **Interpretation**: + - ΔE < 1: Imperceptible difference - ΔE < 2: Perceptible with close observation - ΔE < 10: Noticeable at a glance @@ -192,12 +191,13 @@ This produces more perceptually consistent smoothing compared to RGB Euclidean d ### Computational Cost -| Operation | Complexity | -|:---|:---| +| 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 diff --git a/docs/docs/reference/wasm/modules/image/cielab/index.md b/docs/docs/reference/wasm/modules/image/cielab/index.md index 5d244753c..5268f1612 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/index.md +++ b/docs/docs/reference/wasm/modules/image/cielab/index.md @@ -6,27 +6,30 @@ sidebar_position: 2 # RGB ↔ CIELAB Conversion Guide -This explains the full mathematical conversion pipeline between **sRGB** and **CIELAB (Lab)** color spaces. +This explains the full mathematical conversion pipeline between **sRGB** and **CIELAB (Lab)** color spaces. --- # 1. Conversion Pipeline Overview ## RGB → CIELAB -1. sRGB → Linear RGB -2. Linear RGB → XYZ -3. XYZ → CIELAB + +1. sRGB → Linear RGB +2. Linear RGB → XYZ +3. XYZ → CIELAB ## CIELAB → RGB -1. CIELAB → XYZ -2. XYZ → Linear RGB -3. Linear RGB → 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} = @@ -44,6 +47,7 @@ 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} @@ -60,6 +64,7 @@ R_\text{lin} \\ G_\text{lin} \\ B_\text{lin} \end{bmatrix} ``` + --- # 4. XYZ to CIELAB @@ -69,6 +74,7 @@ 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} ``` @@ -163,20 +169,23 @@ 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\* + +- 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 + +- 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 + +- CIE 1976 L\*a\*b\* Specification +- IEC 61966‑2‑1 sRGB Standard 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 index f95444357..cbe12f6fb 100644 --- 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 @@ -1,10 +1,6 @@ const RgbVsLabRangeKernel = () => (
- + {/* Background */} @@ -13,11 +9,18 @@ const RgbVsLabRangeKernel = () => ( d={(() => { 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 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`; + 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" @@ -29,11 +32,18 @@ const RgbVsLabRangeKernel = () => ( d={(() => { 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 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`; + 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" @@ -41,14 +51,22 @@ const RgbVsLabRangeKernel = () => ( /> {/* Labels */} - RGB LUT - CIELAB (on-the-fly) + + RGB LUT + + + CIELAB (on-the-fly) + {/* Axes */} - Color difference ΔRGB - Color difference ΔLAB + + Color difference ΔRGB + + + Color difference ΔLAB + diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index f1c502d22..bb48c0a26 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -35,8 +35,17 @@ 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 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; diff --git a/src/wasm/modules/image/include/bilateral_filter.h b/src/wasm/modules/image/include/bilateral_filter.h index b5d78d1e7..2b5b938fe 100644 --- a/src/wasm/modules/image/include/bilateral_filter.h +++ b/src/wasm/modules/image/include/bilateral_filter.h @@ -11,8 +11,10 @@ namespace bilateral { // 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) +// - 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); diff --git a/src/wasm/modules/image/include/cielab.h b/src/wasm/modules/image/include/cielab.h index 2dc7e9616..31625f9b4 100644 --- a/src/wasm/modules/image/include/cielab.h +++ b/src/wasm/modules/image/include/cielab.h @@ -4,8 +4,8 @@ #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); + 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); +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/src/bilateral_filter.cpp b/src/wasm/modules/image/src/bilateral_filter.cpp index 71faf9c51..dfe807f7f 100644 --- a/src/wasm/modules/image/src/bilateral_filter.cpp +++ b/src/wasm/modules/image/src/bilateral_filter.cpp @@ -1,37 +1,41 @@ #include "bilateral_filter.h" -#include "exported.h" #include "cielab.h" +#include "exported.h" -#include -#include #include -#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) +// 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)); + 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. +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) +- 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 @@ -39,165 +43,176 @@ It reduces noise in flat regions while preserving edges by assigning near-zero w 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; + // 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 ========= - 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}; + // 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}; - std::vector result(width * height * 4); + double weight_acc{0.0}; + double w_space, w_range; + double dL, dA, dB, dist; - std::vector spatial_weights(kernel_diameter * kernel_diameter); + for (int ky{-radius}; ky <= radius; ++ky) { + int ny{std::clamp(y + ky, 0, h - 1)}; - // 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); + 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; } - } + } - // ========= 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); + 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; + } } } - // ========= 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()); + std::memcpy(image, result.data(), result.size()); } } // namespace bilateral @@ -206,5 +221,6 @@ void bilateral_filter(uint8_t *image, size_t width, size_t height, 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); + 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 index ca8da11e2..cb0b98470 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -1,19 +1,20 @@ #include "cielab.h" -#include #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 +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 +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 ====== // Multipliers for RGB to XYZ @@ -51,83 +52,86 @@ constexpr double LAB_B_FACTOR{200.0}; 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; + 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) +// 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); + 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); + 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) -{ +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) { 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; + ? 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(r); - g = gamma_encode(g); - b = gamma_encode(b); - - // --- 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))); + 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(r); + g = gamma_encode(g); + b = gamma_encode(b); + + // --- 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 0a5f91f04..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] = 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 - }; + 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}; } } From 9538f6e67a8b2f32cdab938102ee9fa9d174f4d2 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 06:10:42 +0200 Subject: [PATCH 48/53] style(docs files): fix display styles --- .../wasm/modules/image/bilateral_filter/explained.md | 6 +++--- .../wasm/modules/image/bilateral_filter/implementation.md | 4 ++-- src/wasm/modules/image/CMakeLists.txt | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md index 0f9475980..642d456ef 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/explained.md @@ -51,9 +51,9 @@ Where each component means: In this implementation, both weighting terms are **Gaussian kernels**: $$ -w_{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) +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. @@ -87,7 +87,7 @@ Since the RGB to CIELAB conversion is expensive, redundant computations are mini 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; diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md index de1a5c523..674ceeb9f 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/implementation.md @@ -92,7 +92,7 @@ For the CIELAB color space, the range weights are computed "on the fly" to reduc :::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)$ +$\exp\left(-\frac{x^2}{2\sigma_{spatial}^2}\right)$ ```cpp double gaussian(double x, double sigma) { @@ -177,7 +177,7 @@ import RgbVsLabRangeKernel from '@site/src/components/docs/reference/wasm/module 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)$ +$\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. diff --git a/src/wasm/modules/image/CMakeLists.txt b/src/wasm/modules/image/CMakeLists.txt index 4d80a5e36..3391f7aa5 100644 --- a/src/wasm/modules/image/CMakeLists.txt +++ b/src/wasm/modules/image/CMakeLists.txt @@ -49,13 +49,13 @@ target_link_options(${MODULE_NAME}_wasm PRIVATE ${COMMON_FLAGS}) # Build-type specific flags if(CMAKE_BUILD_TYPE STREQUAL "Debug") - target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -g4 -ffast-math) + target_compile_options(${MODULE_NAME}_wasm PRIVATE -O0 -g4) target_link_options(${MODULE_NAME}_wasm PRIVATE "SHELL:-s ASSERTIONS=2" -g4 ) else() - target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3 -ffast-math) + target_compile_options(${MODULE_NAME}_wasm PRIVATE -O3) target_link_options(${MODULE_NAME}_wasm PRIVATE "SHELL:-s SINGLE_FILE=0" ) From c89d3ca499700f967db6efbb7296191fe3521f87 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 07:46:28 +0200 Subject: [PATCH 49/53] feat(cielab): improve RGB to Lab documentation and matrix precision MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add Mermaid diagrams for RGB → Lab and Lab → RGB - Update RGB↔XYZ and XYZ↔RGB matrices to high-precision Khronos/W3C values - Add explanatory comments for matrices in cielab.cpp --- .../wasm/modules/image/cielab/index.md | 38 +++++++------- src/wasm/modules/image/src/cielab.cpp | 49 ++++++++++++++----- 2 files changed, 55 insertions(+), 32 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/cielab/index.md b/docs/docs/reference/wasm/modules/image/cielab/index.md index 5268f1612..b1063d0f8 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/index.md +++ b/docs/docs/reference/wasm/modules/image/cielab/index.md @@ -8,24 +8,34 @@ sidebar_position: 2 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: @@ -42,8 +52,6 @@ C_\text{lin} = 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: @@ -65,8 +73,6 @@ R_\text{lin} \\ G_\text{lin} \\ B_\text{lin} ``` ---- - # 4. XYZ to CIELAB Normalize XYZ by the D65 reference white: @@ -103,8 +109,6 @@ a^* = 500 \left[f(x) - f(y)\right] b^* = 200 \left[f(y) - f(z)\right] ``` ---- - # 5. CIELAB to XYZ The inverse of \(f(t)\): @@ -131,8 +135,6 @@ Y = Y_n f^{-1}(f_y),\quad Z = Z_n f^{-1}(f_z) ``` ---- - # 6. XYZ to Linear RGB ```math @@ -141,17 +143,15 @@ R_\text{lin} \\ G_\text{lin} \\ B_\text{lin} \end{bmatrix} = \begin{bmatrix} - 3.2406 & -1.5372 & -0.4986 \\ --0.9689 & 1.8758 & 0.0415 \\ - 0.0557 & -0.2040 & 1.0570 + 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 @@ -164,8 +164,6 @@ C_\text{srgb} = Clamp results to \([0,1]\) and scaled by 255 before converting to 8‑bit. ---- - # 8. Summary ## RGB → Lab @@ -183,8 +181,6 @@ Clamp results to \([0,1]\) and scaled by 255 before converting to 8‑bit. - Linear RGB → sRGB (gamma) - Clamp to valid output ---- - # 9. References - CIE 1976 L\*a\*b\* Specification diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index cb0b98470..aab4732d5 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -16,8 +16,20 @@ 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 ====== -// Multipliers for RGB to XYZ +/* + * ====== 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}; @@ -28,15 +40,30 @@ constexpr double SRGB_R_TO_Z{0.0193339}; constexpr double SRGB_G_TO_Z{0.1191920}; constexpr double SRGB_B_TO_Z{0.9503041}; -constexpr double SRGB_X_TO_R{3.2406}; -constexpr double SRGB_Y_TO_R{-1.5372}; -constexpr double SRGB_Z_TO_R{-0.4986}; -constexpr double SRGB_X_TO_G{-0.9689}; -constexpr double SRGB_Y_TO_G{1.8758}; -constexpr double SRGB_Z_TO_G{0.0415}; -constexpr double SRGB_X_TO_B{0.0557}; -constexpr double SRGB_Y_TO_B{-0.2040}; -constexpr double SRGB_Z_TO_B{1.0570}; +/* + * ====== 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}; From 32f483e4c50609dc80d5451e77ee2eb2130f12f4 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 22:25:47 +0200 Subject: [PATCH 50/53] style(docs: cielab): fix comment formatting - line 47 --- src/wasm/modules/image/src/cielab.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index aab4732d5..bcdc40c94 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -44,7 +44,8 @@ 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. + * Inverse transformation (XYZ → sRGB) per Khronos/W3C, D65 white, slightly + * different from original BT.709 inverse. * * CIE XYZ → sRGB (inverse): * ┌ ┐ ┌ ┐ ┌ ┐ From 14a3aa62c54c2e44a7342180e708e7b0aa70bd1b Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 23:04:04 +0200 Subject: [PATCH 51/53] docs(cpp: bilateral filter): document the difference between sigma_range in the color spaces and recommend a value --- .../modules/image/bilateral_filter/api.md | 18 ++++ .../image/bilateral_filter/color-spaces.md | 94 ++++++++++++++++++- .../wasm/modules/image/cielab/api.md | 2 +- .../{index.md => implementation-explained.md} | 2 +- 4 files changed, 113 insertions(+), 3 deletions(-) rename docs/docs/reference/wasm/modules/image/cielab/{index.md => implementation-explained.md} (99%) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index c85a7fd4c..d706e518c 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -33,3 +33,21 @@ The alpha channel, `image[i + 3]`, is left untouched - it is not part of the bil - **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 index 757b9ac8a..cabf1286b 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md @@ -128,7 +128,7 @@ The scaling factor of **~4.18** is empirically derived and works well for natura - 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 @@ -140,6 +140,98 @@ Using the same `sigma_range = 50`: | **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 diff --git a/docs/docs/reference/wasm/modules/image/cielab/api.md b/docs/docs/reference/wasm/modules/image/cielab/api.md index 39970c829..b8b4af65d 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/api.md +++ b/docs/docs/reference/wasm/modules/image/cielab/api.md @@ -2,7 +2,7 @@ id: cielab-api title: CIELAB Color Space API sidebar_label: API Reference -sidebar_position: 1 +sidebar_position: 2 --- # CIELAB Color Space Conversion API diff --git a/docs/docs/reference/wasm/modules/image/cielab/index.md b/docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md similarity index 99% rename from docs/docs/reference/wasm/modules/image/cielab/index.md rename to docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md index b1063d0f8..de5588d1e 100644 --- a/docs/docs/reference/wasm/modules/image/cielab/index.md +++ b/docs/docs/reference/wasm/modules/image/cielab/implementation-explained.md @@ -1,7 +1,7 @@ --- id: explained title: Implementation Explained -sidebar_position: 2 +sidebar_position: 1 --- # RGB ↔ CIELAB Conversion Guide From 45f62748c2ca13c1e442e83c7b067d914997c03d Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 23:06:24 +0200 Subject: [PATCH 52/53] style(docs: bilateral filter): fix formatting --- .../modules/image/bilateral_filter/api.md | 6 ++--- .../image/bilateral_filter/color-spaces.md | 26 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md index d706e518c..37485ec28 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/api.md @@ -48,6 +48,6 @@ $$ $$ > 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. -::: +> 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 index cabf1286b..ec4ded405 100644 --- a/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md +++ b/docs/docs/reference/wasm/modules/image/bilateral_filter/color-spaces.md @@ -128,7 +128,7 @@ The scaling factor of **~4.18** is empirically derived and works well for natura - 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 @@ -140,7 +140,6 @@ Using the same `sigma_range = 50`: | **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. @@ -183,11 +182,11 @@ So the scaling factor is **not arbitrary** — it is the **ratio of typical RGB For natural images (photographic content, sRGB, D65): -1. Sample many *local* pixel pairs (neighbors). +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}$$ - * $$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). @@ -203,29 +202,30 @@ The value **4.18** lies near the center of this empirical range and provides a s The factor remains stable for natural images because: -* **LAB compresses perceptual differences** +- **LAB compresses perceptual differences** Equal perceived color changes produce smaller numeric deltas than in RGB. -* **RGB channels are highly correlated** +- **RGB channels are highly correlated** Euclidean RGB distance accumulates redundant energy across channels. -* **Bilateral filters operate locally** +- **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 +- 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** +- **Recommended default** For visually comparable smoothing on natural images, use: + $$ \sigma_{range_{RGB}} \approx 4.18 \times \sigma_{range_{CIELAB}} $$ -* **Advanced usage** +- **Advanced usage** For strict equivalence, compute the ratio $$ k = \frac{\mathbb{E}[d_{RGB}]}{\mathbb{E}[d_{LAB}]} From abb0402d28d1f4019ac2dbf2dccffebbfceba444 Mon Sep 17 00:00:00 2001 From: Ryan-Millard Date: Tue, 6 Jan 2026 23:30:27 +0200 Subject: [PATCH 53/53] fix(cpp: cielab.h): add clamp values to protect against bad data --- src/wasm/modules/image/src/cielab.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/wasm/modules/image/src/cielab.cpp b/src/wasm/modules/image/src/cielab.cpp index bcdc40c94..073594f4e 100644 --- a/src/wasm/modules/image/src/cielab.cpp +++ b/src/wasm/modules/image/src/cielab.cpp @@ -131,6 +131,8 @@ constexpr double inverse_xyz_to_lab(double t) { } 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) - @@ -154,9 +156,9 @@ void lab_to_rgb(const double L, const double A, const double B, 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(r); - g = gamma_encode(g); - b = gamma_encode(b); + 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)));