M20: Deep learning, large images, GPU shaders, parallelism, interoperability (638 tests)
1. Neural Net Inference (new ri-nn crate): Conv2D, ReLU, Sigmoid, Softmax, MaxPool2D, BatchNorm, Dense, Upsample, U-Net architecture, StarDist-like detector — all from scratch, no ML deps 2. Virtual Stacks (new ri-virtual crate): VirtualStack with lazy loading + LRU cache, TiledImage for large 2D, multi-resolution pyramid 3. GPU Shaders: 7 new wgsl shaders (GLCM, optical flow, EDT, connected components, color deconvolution, separable convolution, 2D histogram) 4. Image Sequences: load/save numbered files, natural sorting, MJPEG AVI, Y4M export with ffmpeg command generation 5. ImageJ Interop: .roi binary format read/write, ROI zip sets, .lut format, Results Table CSV import/export 6. Parallel Processing: rayon for bilateral, median, EDT, connected components, particle analysis, ML features, random forest — set_parallelism config 7. Quality Metrics: PSNR, SSIM, MS-SSIM, NRMSE, UQI, VIF, noise estimation, SNR, contrast-to-noise ratio 8. Integration Tests (new ri-integration-tests crate): 15 end-to-end tests covering full pipelines across multiple crates
This commit is contained in:
osobh
@@ -0,0 +1,480 @@
|
||||
//! Image quality metrics for comparing reference and test images.
|
||||
//!
|
||||
//! Includes PSNR, SSIM, MS-SSIM, NRMSE, UQI, VIF, noise estimation, SNR, and CNR.
|
||||
|
||||
use ri_core::TypedBuffer;
|
||||
|
||||
fn to_f64_vec(buf: &TypedBuffer) -> Vec<f64> {
|
||||
match buf {
|
||||
TypedBuffer::U8(v) => v.iter().map(|&x| x as f64).collect(),
|
||||
TypedBuffer::U16(v) => v.iter().map(|&x| x as f64).collect(),
|
||||
TypedBuffer::F32(v) => v.iter().map(|&x| x as f64).collect(),
|
||||
}
|
||||
}
|
||||
|
||||
fn max_pixel_value(buf: &TypedBuffer) -> f64 {
|
||||
match buf {
|
||||
TypedBuffer::U8(_) => 255.0,
|
||||
TypedBuffer::U16(_) => 65535.0,
|
||||
TypedBuffer::F32(_) => 1.0,
|
||||
}
|
||||
}
|
||||
|
||||
/// Peak Signal-to-Noise Ratio (PSNR) in dB.
|
||||
pub fn psnr(reference: &TypedBuffer, test: &TypedBuffer) -> f64 {
|
||||
let ref_data = to_f64_vec(reference);
|
||||
let test_data = to_f64_vec(test);
|
||||
assert_eq!(ref_data.len(), test_data.len());
|
||||
|
||||
let mse: f64 = ref_data
|
||||
.iter()
|
||||
.zip(test_data.iter())
|
||||
.map(|(r, t)| (r - t) * (r - t))
|
||||
.sum::<f64>()
|
||||
/ ref_data.len() as f64;
|
||||
|
||||
if mse == 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
|
||||
let max_val = max_pixel_value(reference);
|
||||
10.0 * (max_val * max_val / mse).log10()
|
||||
}
|
||||
|
||||
/// Structural Similarity Index (SSIM).
|
||||
/// Returns mean SSIM over all local windows.
|
||||
pub fn ssim(reference: &TypedBuffer, test: &TypedBuffer, width: u32, height: u32) -> f64 {
|
||||
ssim_map(reference, test, width, height).0
|
||||
}
|
||||
|
||||
/// Compute SSIM returning (mean_ssim, ssim_map).
|
||||
fn ssim_map(reference: &TypedBuffer, test: &TypedBuffer, width: u32, height: u32) -> (f64, Vec<f64>) {
|
||||
let ref_data = to_f64_vec(reference);
|
||||
let test_data = to_f64_vec(test);
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
let max_val = max_pixel_value(reference);
|
||||
|
||||
let c1 = (0.01 * max_val) * (0.01 * max_val);
|
||||
let c2 = (0.03 * max_val) * (0.03 * max_val);
|
||||
let window_size = 7usize;
|
||||
let half = window_size / 2;
|
||||
|
||||
let mut ssim_values = Vec::new();
|
||||
|
||||
for y in half..h.saturating_sub(half) {
|
||||
for x in half..w.saturating_sub(half) {
|
||||
let mut sum_r = 0.0;
|
||||
let mut sum_t = 0.0;
|
||||
let mut sum_r2 = 0.0;
|
||||
let mut sum_t2 = 0.0;
|
||||
let mut sum_rt = 0.0;
|
||||
let mut count = 0.0;
|
||||
|
||||
for wy in 0..window_size {
|
||||
for wx in 0..window_size {
|
||||
let py = y - half + wy;
|
||||
let px = x - half + wx;
|
||||
if py < h && px < w {
|
||||
let r = ref_data[py * w + px];
|
||||
let t = test_data[py * w + px];
|
||||
sum_r += r;
|
||||
sum_t += t;
|
||||
sum_r2 += r * r;
|
||||
sum_t2 += t * t;
|
||||
sum_rt += r * t;
|
||||
count += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
let mu_r = sum_r / count;
|
||||
let mu_t = sum_t / count;
|
||||
let sigma_r2 = sum_r2 / count - mu_r * mu_r;
|
||||
let sigma_t2 = sum_t2 / count - mu_t * mu_t;
|
||||
let sigma_rt = sum_rt / count - mu_r * mu_t;
|
||||
|
||||
let numerator = (2.0 * mu_r * mu_t + c1) * (2.0 * sigma_rt + c2);
|
||||
let denominator = (mu_r * mu_r + mu_t * mu_t + c1) * (sigma_r2 + sigma_t2 + c2);
|
||||
|
||||
ssim_values.push(numerator / denominator);
|
||||
}
|
||||
}
|
||||
|
||||
let mean = if ssim_values.is_empty() {
|
||||
1.0
|
||||
} else {
|
||||
ssim_values.iter().sum::<f64>() / ssim_values.len() as f64
|
||||
};
|
||||
|
||||
(mean, ssim_values)
|
||||
}
|
||||
|
||||
/// Multi-Scale SSIM.
|
||||
pub fn ms_ssim(reference: &TypedBuffer, test: &TypedBuffer, width: u32, height: u32) -> f64 {
|
||||
let weights = [0.0448, 0.2856, 0.3001, 0.2363, 0.1333];
|
||||
let mut ref_data = to_f64_vec(reference);
|
||||
let mut test_data = to_f64_vec(test);
|
||||
let mut w = width as usize;
|
||||
let mut h = height as usize;
|
||||
|
||||
let mut result = 1.0;
|
||||
|
||||
for (level, &weight) in weights.iter().enumerate() {
|
||||
if w < 8 || h < 8 {
|
||||
break;
|
||||
}
|
||||
|
||||
let ref_buf = TypedBuffer::F32(ref_data.iter().map(|&v| v as f32).collect());
|
||||
let test_buf = TypedBuffer::F32(test_data.iter().map(|&v| v as f32).collect());
|
||||
let (s, _) = ssim_map(&ref_buf, &test_buf, w as u32, h as u32);
|
||||
|
||||
if level == weights.len() - 1 {
|
||||
result *= s.powf(weight);
|
||||
} else {
|
||||
// Just use SSIM contrast/structure for intermediate scales
|
||||
result *= s.max(0.0).powf(weight);
|
||||
}
|
||||
|
||||
// Downsample by 2
|
||||
let new_w = w / 2;
|
||||
let new_h = h / 2;
|
||||
if new_w == 0 || new_h == 0 { break; }
|
||||
|
||||
let mut new_ref = vec![0.0; new_w * new_h];
|
||||
let mut new_test = vec![0.0; new_w * new_h];
|
||||
for y in 0..new_h {
|
||||
for x in 0..new_w {
|
||||
let sy = y * 2;
|
||||
let sx = x * 2;
|
||||
new_ref[y * new_w + x] = (ref_data[sy * w + sx]
|
||||
+ ref_data[sy * w + sx + 1]
|
||||
+ ref_data[(sy + 1) * w + sx]
|
||||
+ ref_data[(sy + 1) * w + sx + 1]) / 4.0;
|
||||
new_test[y * new_w + x] = (test_data[sy * w + sx]
|
||||
+ test_data[sy * w + sx + 1]
|
||||
+ test_data[(sy + 1) * w + sx]
|
||||
+ test_data[(sy + 1) * w + sx + 1]) / 4.0;
|
||||
}
|
||||
}
|
||||
ref_data = new_ref;
|
||||
test_data = new_test;
|
||||
w = new_w;
|
||||
h = new_h;
|
||||
}
|
||||
|
||||
result
|
||||
}
|
||||
|
||||
/// Normalized Root Mean Square Error.
|
||||
pub fn nrmse(reference: &TypedBuffer, test: &TypedBuffer) -> f64 {
|
||||
let ref_data = to_f64_vec(reference);
|
||||
let test_data = to_f64_vec(test);
|
||||
assert_eq!(ref_data.len(), test_data.len());
|
||||
|
||||
let mse: f64 = ref_data
|
||||
.iter()
|
||||
.zip(test_data.iter())
|
||||
.map(|(r, t)| (r - t) * (r - t))
|
||||
.sum::<f64>()
|
||||
/ ref_data.len() as f64;
|
||||
let rmse = mse.sqrt();
|
||||
|
||||
let ref_min = ref_data.iter().cloned().fold(f64::MAX, f64::min);
|
||||
let ref_max = ref_data.iter().cloned().fold(f64::MIN, f64::max);
|
||||
let range = ref_max - ref_min;
|
||||
|
||||
if range == 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
rmse / range
|
||||
}
|
||||
|
||||
/// Universal Quality Index (UQI).
|
||||
pub fn uqi(reference: &TypedBuffer, test: &TypedBuffer) -> f64 {
|
||||
let ref_data = to_f64_vec(reference);
|
||||
let test_data = to_f64_vec(test);
|
||||
let n = ref_data.len() as f64;
|
||||
|
||||
let mu_r = ref_data.iter().sum::<f64>() / n;
|
||||
let mu_t = test_data.iter().sum::<f64>() / n;
|
||||
|
||||
let sigma_r2: f64 = ref_data.iter().map(|&r| (r - mu_r) * (r - mu_r)).sum::<f64>() / n;
|
||||
let sigma_t2: f64 = test_data.iter().map(|&t| (t - mu_t) * (t - mu_t)).sum::<f64>() / n;
|
||||
let sigma_rt: f64 = ref_data.iter().zip(test_data.iter())
|
||||
.map(|(&r, &t)| (r - mu_r) * (t - mu_t)).sum::<f64>() / n;
|
||||
|
||||
let numerator = 4.0 * sigma_rt * mu_r * mu_t;
|
||||
let denominator = (sigma_r2 + sigma_t2) * (mu_r * mu_r + mu_t * mu_t);
|
||||
|
||||
if denominator == 0.0 { 1.0 } else { numerator / denominator }
|
||||
}
|
||||
|
||||
/// Simplified Visual Information Fidelity (VIF).
|
||||
/// Uses a simple single-scale approach based on local statistics.
|
||||
pub fn vif(reference: &TypedBuffer, test: &TypedBuffer, width: u32, height: u32) -> f64 {
|
||||
let ref_data = to_f64_vec(reference);
|
||||
let test_data = to_f64_vec(test);
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
let block = 3usize;
|
||||
|
||||
let mut num = 0.0;
|
||||
let mut den = 0.0;
|
||||
let sigma_nsq = 2.0; // noise variance estimate
|
||||
|
||||
for y in (0..h - block).step_by(block) {
|
||||
for x in (0..w - block).step_by(block) {
|
||||
let mut sum_r = 0.0;
|
||||
let mut sum_t = 0.0;
|
||||
let mut sum_r2 = 0.0;
|
||||
let mut sum_t2 = 0.0;
|
||||
let mut sum_rt = 0.0;
|
||||
let n = (block * block) as f64;
|
||||
|
||||
for by in 0..block {
|
||||
for bx in 0..block {
|
||||
let r = ref_data[(y + by) * w + (x + bx)];
|
||||
let t = test_data[(y + by) * w + (x + bx)];
|
||||
sum_r += r;
|
||||
sum_t += t;
|
||||
sum_r2 += r * r;
|
||||
sum_t2 += t * t;
|
||||
sum_rt += r * t;
|
||||
}
|
||||
}
|
||||
|
||||
let mu_r = sum_r / n;
|
||||
let mu_t = sum_t / n;
|
||||
let sigma_r2 = (sum_r2 / n - mu_r * mu_r).max(0.0);
|
||||
let sigma_t2 = (sum_t2 / n - mu_t * mu_t).max(0.0);
|
||||
let sigma_rt = sum_rt / n - mu_r * mu_t;
|
||||
|
||||
let g = if sigma_r2 > 0.0 { sigma_rt / sigma_r2 } else { 0.0 };
|
||||
let sigma_v2 = (sigma_t2 - g * sigma_rt).max(0.0);
|
||||
|
||||
num += (1.0 + (g * g * sigma_r2) / (sigma_v2 + sigma_nsq)).ln();
|
||||
den += (1.0 + sigma_r2 / sigma_nsq).ln();
|
||||
}
|
||||
}
|
||||
|
||||
if den == 0.0 { 1.0 } else { num / den }
|
||||
}
|
||||
|
||||
/// Estimate noise level using MAD (Median Absolute Deviation) of the Laplacian.
|
||||
pub fn noise_estimation(buffer: &TypedBuffer, width: u32, height: u32) -> f64 {
|
||||
let data = to_f64_vec(buffer);
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
|
||||
// Apply Laplacian
|
||||
let mut laplacian_values = Vec::with_capacity((w - 2) * (h - 2));
|
||||
for y in 1..h - 1 {
|
||||
for x in 1..w - 1 {
|
||||
let lap = -data[(y - 1) * w + x]
|
||||
- data[y * w + (x - 1)]
|
||||
+ 4.0 * data[y * w + x]
|
||||
- data[y * w + (x + 1)]
|
||||
- data[(y + 1) * w + x];
|
||||
laplacian_values.push(lap.abs());
|
||||
}
|
||||
}
|
||||
|
||||
laplacian_values.sort_by(|a, b| a.partial_cmp(b).unwrap());
|
||||
let median = if laplacian_values.is_empty() {
|
||||
0.0
|
||||
} else {
|
||||
laplacian_values[laplacian_values.len() / 2]
|
||||
};
|
||||
|
||||
// MAD estimator: sigma = median / 0.6745 * sqrt(2) factor for Laplacian
|
||||
median * 1.4826 / (2.0f64.sqrt())
|
||||
}
|
||||
|
||||
/// Signal-to-Noise Ratio estimate (mean / std_dev).
|
||||
pub fn snr(buffer: &TypedBuffer) -> f64 {
|
||||
let data = to_f64_vec(buffer);
|
||||
let n = data.len() as f64;
|
||||
let mean = data.iter().sum::<f64>() / n;
|
||||
let variance = data.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>() / n;
|
||||
let std_dev = variance.sqrt();
|
||||
|
||||
if std_dev == 0.0 { f64::INFINITY } else { mean / std_dev }
|
||||
}
|
||||
|
||||
/// Contrast-to-Noise Ratio between a signal region and a background region.
|
||||
/// `signal_roi` and `background_roi` are lists of pixel indices.
|
||||
pub fn contrast_to_noise(
|
||||
buffer: &TypedBuffer,
|
||||
signal_pixels: &[usize],
|
||||
background_pixels: &[usize],
|
||||
) -> f64 {
|
||||
let data = to_f64_vec(buffer);
|
||||
|
||||
let signal_mean = signal_pixels.iter()
|
||||
.map(|&i| data[i])
|
||||
.sum::<f64>() / signal_pixels.len() as f64;
|
||||
|
||||
let bg_mean = background_pixels.iter()
|
||||
.map(|&i| data[i])
|
||||
.sum::<f64>() / background_pixels.len() as f64;
|
||||
|
||||
let bg_variance = background_pixels.iter()
|
||||
.map(|&i| (data[i] - bg_mean) * (data[i] - bg_mean))
|
||||
.sum::<f64>() / background_pixels.len() as f64;
|
||||
|
||||
let bg_std = bg_variance.sqrt();
|
||||
|
||||
if bg_std == 0.0 {
|
||||
f64::INFINITY
|
||||
} else {
|
||||
(signal_mean - bg_mean).abs() / bg_std
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn psnr_identical() {
|
||||
let buf = TypedBuffer::U8(vec![100; 64]);
|
||||
let p = psnr(&buf, &buf);
|
||||
assert!(p.is_infinite());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn psnr_different() {
|
||||
let ref_buf = TypedBuffer::U8(vec![100; 64]);
|
||||
let test_buf = TypedBuffer::U8(vec![110; 64]);
|
||||
let p = psnr(&ref_buf, &test_buf);
|
||||
assert!(p > 0.0);
|
||||
assert!(p < 100.0);
|
||||
// Expected: 10 * log10(255^2 / 100) ≈ 28.13
|
||||
assert!((p - 28.13).abs() < 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ssim_identical() {
|
||||
let buf = TypedBuffer::U8(vec![128; 64]);
|
||||
let s = ssim(&buf, &buf, 8, 8);
|
||||
assert!((s - 1.0).abs() < 0.01);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ssim_different() {
|
||||
let ref_buf = TypedBuffer::U8(vec![100; 64]);
|
||||
let test_buf = TypedBuffer::U8(vec![200; 64]);
|
||||
let s = ssim(&ref_buf, &test_buf, 8, 8);
|
||||
assert!(s < 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ms_ssim_identical() {
|
||||
let buf = TypedBuffer::U8(vec![128; 128 * 128]);
|
||||
let s = ms_ssim(&buf, &buf, 128, 128);
|
||||
assert!((s - 1.0).abs() < 0.01);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nrmse_identical() {
|
||||
let buf = TypedBuffer::U8(vec![100; 64]);
|
||||
let n = nrmse(&buf, &buf);
|
||||
assert!((n - 0.0).abs() < 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nrmse_different() {
|
||||
let ref_buf = TypedBuffer::U8((0..64).map(|i| i as u8 * 4).collect());
|
||||
let test_buf = TypedBuffer::U8((0..64).map(|i| (i as u8 * 4).wrapping_add(10)).collect());
|
||||
let n = nrmse(&ref_buf, &test_buf);
|
||||
assert!(n > 0.0);
|
||||
assert!(n < 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn uqi_identical() {
|
||||
let buf = TypedBuffer::U8(vec![100, 150, 200, 50]);
|
||||
let q = uqi(&buf, &buf);
|
||||
assert!((q - 1.0).abs() < 0.01);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn uqi_different() {
|
||||
let ref_buf = TypedBuffer::U8(vec![100, 150, 200, 50]);
|
||||
let test_buf = TypedBuffer::U8(vec![50, 200, 100, 150]);
|
||||
let q = uqi(&ref_buf, &test_buf);
|
||||
assert!(q < 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn vif_identical() {
|
||||
let buf = TypedBuffer::F32((0..81).map(|i| i as f32 / 81.0).collect());
|
||||
let v = vif(&buf, &buf, 9, 9);
|
||||
assert!(v > 0.9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noise_estimation_uniform() {
|
||||
let buf = TypedBuffer::U8(vec![128; 100]);
|
||||
let n = noise_estimation(&buf, 10, 10);
|
||||
assert!(n < 1.0); // Uniform image → near-zero noise
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noise_estimation_noisy() {
|
||||
// Create a noisy image
|
||||
let mut data: Vec<u8> = Vec::with_capacity(100);
|
||||
for i in 0..100 {
|
||||
data.push(if i % 2 == 0 { 100 } else { 200 });
|
||||
}
|
||||
let buf = TypedBuffer::U8(data);
|
||||
let n = noise_estimation(&buf, 10, 10);
|
||||
assert!(n > 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn snr_uniform() {
|
||||
let buf = TypedBuffer::U8(vec![128; 100]);
|
||||
let s = snr(&buf);
|
||||
assert!(s.is_infinite()); // No variance → infinite SNR
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn snr_variable() {
|
||||
let buf = TypedBuffer::F32(vec![10.0, 12.0, 8.0, 11.0, 9.0]);
|
||||
let s = snr(&buf);
|
||||
assert!(s > 0.0);
|
||||
assert!(s.is_finite());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cnr_basic() {
|
||||
let buf = TypedBuffer::F32(vec![
|
||||
100.0, 102.0, 98.0, 101.0, // Signal region (indices 0-3)
|
||||
10.0, 12.0, 8.0, 11.0, // Background region (indices 4-7)
|
||||
]);
|
||||
let signal = vec![0, 1, 2, 3];
|
||||
let background = vec![4, 5, 6, 7];
|
||||
let c = contrast_to_noise(&buf, &signal, &background);
|
||||
assert!(c > 10.0); // Large contrast relative to background noise
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn psnr_f32() {
|
||||
let ref_buf = TypedBuffer::F32(vec![0.5; 64]);
|
||||
let test_buf = TypedBuffer::F32(vec![0.6; 64]);
|
||||
let p = psnr(&ref_buf, &test_buf);
|
||||
assert!(p > 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn ssim_window_edge() {
|
||||
// Very small image where SSIM window can barely fit
|
||||
let buf = TypedBuffer::U8(vec![128; 49]);
|
||||
let s = ssim(&buf, &buf, 7, 7);
|
||||
// Should still produce a reasonable result
|
||||
assert!(s >= 0.0);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user