Skip to content

Bilateral speedup #798

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CattleProdigy wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Bilateral speedup #798

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ef82400
Add LUT-based color-distance for bilateral filter
CattleProdigy Jul 10, 2026
f6fad12
Rework bilateral filter to minimizing clamping compute
CattleProdigy Jul 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
247 changes: 208 additions & 39 deletions src/filter/bilateral.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
//! Bilateral Filter and associated items.

use image::{GenericImage, Pixel};
use image::{GenericImage, ImageBuffer, Luma, Pixel, Rgb, Rgba};
use num::cast::AsPrimitive;

use crate::definitions::Image;
Expand Down Expand Up @@ -56,6 +56,158 @@ where
}
}

/// A gaussian function of the euclidean distance between two pixel's colors, implemented using
/// a look up table. This gives substantial speed up for pixel's with a u8 subpixel.
///
/// This implements [`ColorDistance`].
pub struct LutGaussianEuclideanColorDistance {

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

needs testing

lut: Vec<f32>,
}

impl LutGaussianEuclideanColorDistance {
/// Creates a new [`GaussianEuclideanColorDistance`] using a given sigma value, which

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[LutGaussianEuclideanColorDistance]

/// must be positive.
///
/// Internally, this is stored as sigma squared for performance.
///
/// # Panics
///
/// 1. If `sigma <= 0`
pub fn new(sigma: f32) -> Self {
assert!(
sigma > 0.0,
"GaussianEuclideanColorDistance sigma must be positive"

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LutGaussianEuclideanColorDistance

);
let mut lut = Vec::with_capacity(512);

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe lets allocate [f32; 512] on the stack?

for diff in 0..512 {
let diff = (diff - 255) as f32;
let weight = fast_exp_negative(-0.5 * diff.powi(2) / sigma.powi(2));
lut.push(weight);
}
Self { lut }
}
}

impl ColorDistance<Luma<u8>> for LutGaussianEuclideanColorDistance {
#[inline(always)]
fn color_distance(&self, pixel1: &Luma<u8>, pixel2: &Luma<u8>) -> f32 {
let diff = (255 + pixel1.0[0] as usize) - pixel2.0[0] as usize;
self.lut[diff]
}
}

impl ColorDistance<Rgb<u8>> for LutGaussianEuclideanColorDistance {
#[inline(always)]
fn color_distance(&self, pixel1: &Rgb<u8>, pixel2: &Rgb<u8>) -> f32 {
let diffr = (255 + pixel1.0[0] as usize) - pixel2.0[0] as usize;
let diffg = (255 + pixel1.0[1] as usize) - pixel2.0[1] as usize;
let diffb = (255 + pixel1.0[2] as usize) - pixel2.0[2] as usize;
let lutr = self.lut[diffr];
let lutg = self.lut[diffg];
let lutb = self.lut[diffb];
lutr * lutg * lutb
}
}

impl ColorDistance<Rgba<u8>> for LutGaussianEuclideanColorDistance {
#[inline(always)]
fn color_distance(&self, pixel1: &Rgba<u8>, pixel2: &Rgba<u8>) -> f32 {
let diffr = (255 + pixel1.0[0] as usize) - pixel2.0[0] as usize;
let diffg = (255 + pixel1.0[1] as usize) - pixel2.0[1] as usize;
let diffb = (255 + pixel1.0[2] as usize) - pixel2.0[2] as usize;
let diffa = (255 + pixel1.0[3] as usize) - pixel2.0[3] as usize;
let lutr = self.lut[diffr];
let lutg = self.lut[diffg];
let lutb = self.lut[diffb];
let luta = self.lut[diffa];
lutr * lutg * lutb * luta
}
}

/// Loop-invariant context for a single `bilateral_filter` run.
///
/// The per-pixel work is identical for interior and border pixels except for how the
/// window coordinates are computed, so it is expressed once in [`BilateralCtx::filter`]
/// and specialized at compile time via the `CLAMP` const generic.
struct BilateralCtx<'a, I, C> {
image: &'a I,
color_distance: &'a C,
spatial_distance_lookup: &'a [f32],
radius: u32,
radius_range: u32,
width: u32,
height: u32,
}

impl<I, P, C> BilateralCtx<'_, I, C>
where
I: GenericImage<Pixel = P>,
P: Pixel,
C: ColorDistance<P>,
<P as image::Pixel>::Subpixel: 'static,
f32: From<P::Subpixel> + AsPrimitive<P::Subpixel>,
{
/// Filter a single output pixel.
///
/// When `CLAMP` is `true` the window coordinates are clamped into the image bounds
/// (used for border pixels); when `false` the window is assumed fully in-bounds
/// (used for the interior). Because `CLAMP` is a const generic, each variant
/// monomorphizes to straight-line code with the unused branch eliminated.
#[inline(always)]
fn filter<const CLAMP: bool>(&self, x: u32, y: u32) -> P {

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should be unsafe

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and needs testing

const MAX_CHANNELS: usize = 4;

debug_assert!(self.image.in_bounds(x, y));
// Safety: callers only pass (x, y) within the image bounds.
let center_pixel = unsafe { self.image.unsafe_get_pixel(x, y) };

let mut channel_sums = [0f32; MAX_CHANNELS];
let mut weight_sum = 0f32;

// Both branches walk the window in the same (w_y, w_x) order that
// `spatial_distance_lookup` was built in, so a sequential iterator gives the
// right weight without the per-pixel `window_len * w_y + w_x` index math and
// bounds check.
let mut spatial = self.spatial_distance_lookup.iter();
for w_y in 0..self.radius_range {
for w_x in 0..self.radius_range {
let (window_x, window_y) = if CLAMP {
(
(x + w_x).saturating_sub(self.radius).min(self.width - 1),
(y + w_y).saturating_sub(self.radius).min(self.height - 1),
)
} else {
((x + w_x) - self.radius, (y + w_y) - self.radius)
};

debug_assert!(self.image.in_bounds(window_x, window_y));
// Safety: for CLAMP=true the coords are clamped in-bounds; for CLAMP=false
// the interior loop ranges guarantee the whole window is in-bounds.
let window_pixel = unsafe { self.image.unsafe_get_pixel(window_x, window_y) };

let spatial_weight = spatial.next().unwrap();
let color_weight = self
.color_distance
.color_distance(&center_pixel, &window_pixel);
let weight = spatial_weight * color_weight;

weight_sum += weight;
for (i, c) in window_pixel.channels().iter().enumerate() {
channel_sums[i] += weight * f32::from(*c);
}
}
}

let mut out_pixel = center_pixel;
let num_channels = P::CHANNEL_COUNT as usize;
let out_channels = out_pixel.channels_mut();
for i in 0..num_channels {
out_channels[i] = (channel_sums[i] / weight_sum).as_();

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably faster with zip because out_channels is a slice, plus will not go out of bounds

}
out_pixel
}
}

/// Denoise an 8-bit image while preserving edges using bilateral filtering.
///
/// # Arguments
Expand Down Expand Up @@ -124,7 +276,9 @@ where
assert_ne!(image.height(), 0);
assert!(spatial_sigma > 0.0, "spatial_sigma must be positive");

let radius = i16::from(radius);
let (width, height) = image.dimensions();
let radius = (radius as u32).min(width).min(height);

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should document that radius is clamped to be in bounds

let radius = radius as i16;

let spatial_sigma_squared = spatial_sigma.powi(2);
let mut spatial_distance_lookup =
Expand All @@ -138,54 +292,69 @@ where
}
}

let (width, height) = image.dimensions();
let window_len = 2 * radius + 1;

let bilateral_pixel_filter = |x, y| {
debug_assert!(image.in_bounds(x, y));
// Safety: `Image::from_fn` yields `col` in [0, width) and `row` in [0, height).
let center_pixel = unsafe { image.unsafe_get_pixel(x, y) };
let radius = radius as u32;
let radius_range = 2 * radius + 1;

let mut channel_sums = [0f32; MAX_CHANNELS];
let mut weight_sum = 0f32;
let mut out_image = ImageBuffer::new(width, height);

for w_y in -radius..=radius {
for w_x in -radius..=radius {
// these casts will always be correct due to asserts made at the beginning of the
// function about the image width/height
//
// the subtraction will also never overflow due to the `is_empty()` assert
let window_y = (i32::from(w_y) + (y as i32)).clamp(0, (height as i32) - 1);
let window_x = (i32::from(w_x) + (x as i32)).clamp(0, (width as i32) - 1);
let ctx = BilateralCtx {
image,
color_distance: &color_distance,
spatial_distance_lookup: &spatial_distance_lookup,
radius,
radius_range,
width,
height,
};

let (window_y, window_x) = (window_y as u32, window_x as u32);
// Top edge
for y in 0..radius {
for x in 0..width {
let val = ctx.filter::<true>(x, y);
unsafe {
out_image.unsafe_put_pixel(x, y, val);
}
}
}

debug_assert!(image.in_bounds(window_x, window_y));
// Safety: we clamped `window_x` and `window_y` to be in bounds.
let window_pixel = unsafe { image.unsafe_get_pixel(window_x, window_y) };
// Middle Rows
for y in radius..(height - radius) {
// Left Edge
for x in 0..radius {
let val = ctx.filter::<true>(x, y);
unsafe {
out_image.unsafe_put_pixel(x, y, val);
}
}

let spatial_weight = spatial_distance_lookup
[(window_len * (w_y + radius) + (w_x + radius)) as usize];
let color_weight = color_distance.color_distance(&center_pixel, &window_pixel);
let weight = spatial_weight * color_weight;
// Middle, unclamped section
for x in radius..(width - radius) {
let val = ctx.filter::<false>(x, y);
unsafe {
out_image.unsafe_put_pixel(x, y, val);
}
}

weight_sum += weight;
for (i, c) in window_pixel.channels().iter().enumerate() {
channel_sums[i] += weight * f32::from(*c);
}
// Right edge
for x in (width - radius)..width {

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

overlaps similar to comment below

let val = ctx.filter::<true>(x, y);
unsafe {
out_image.unsafe_put_pixel(x, y, val);
}
}
}

let mut out_pixel = center_pixel;
let num_channels = P::CHANNEL_COUNT as usize;
let out_channels = out_pixel.channels_mut();
for i in 0..num_channels {
out_channels[i] = (channel_sums[i] / weight_sum).as_();
// Bottom edge
for y in (height - radius)..height {

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this overlaps with the previous loop of y in 0..radius, for example:

radius = 7
height = 10

0..radius <=> 0..7
(height - radius)..height <=> (10 - 7)..10 <=> 3..10

@cospectrum cospectrum Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will do

debug_assert!(radius <= height)
let interior_end = (height - radius).max(radius);

for y in 0..radius {}

for y in radius..interior_end {}

for y in interior_end..height {}

for x in 0..width {
let val = ctx.filter::<true>(x, y);
unsafe {
out_image.unsafe_put_pixel(x, y, val);
}
}
out_pixel
};
}

Image::from_fn(width, height, bilateral_pixel_filter)
out_image
}

/// Un-normalized Gaussian Weight
Expand Down
Loading