M19: Scientific imaging — ML segmentation, 3D ops, histology, registration, time series (506 tests)
1. ML Segmentation (new ri-ml crate): multi-scale feature extraction (Gaussian/ LoG/Hessian/Sobel/DoG at 5 scales), random forest classifier from scratch (decision trees with Gini impurity, bagging, random feature subsets), train_from_rois, predict probabilities, save/load model 2. 3D Operations (new ri-3d-ops crate): gaussian_blur_3d (separable), median_filter_3d, dilate_3d/erode_3d, distance_transform_3d, connected_components_3d (6/26 connectivity), measure_3d_objects (volume/surface/centroid/sphericity), resample_3d (trilinear) 3. Color Deconvolution: Beer-Lambert unmixing, presets (H&E, H-DAB, Masson Trichrome, Alcian Blue), auto stain vector estimation (PCA/NMF) 4. Advanced Registration: landmark-based (rigid/similarity/affine via least squares), B-spline non-rigid, Thirion's demons deformable 5. Time Series: bleach correction (ratio/exponential/histogram matching), temporal median/difference, delta_F/F0 for calcium imaging, reslice, orthogonal views 6. OME-TIFF: XML metadata parser, multi-series support, tiled TIFF, bio_formats_metadata common struct (pixel_size, z_step, channels) 7. Frequency Domain: notch filter, homomorphic filter, spectral analysis, autocorrelation, cross-correlation, phase symmetry detector 8. Advanced Measurements: fractal dimension (box-counting), lacunarity, orientation analysis (structure tensor), granulometry, radial distribution, 3D intensity profiles
This commit is contained in:
osobh
@@ -0,0 +1,666 @@
|
||||
//! Advanced image registration: landmark-based, B-spline, and demons.
|
||||
|
||||
use ri_core::TypedBuffer;
|
||||
|
||||
/// Transform types for landmark registration.
|
||||
#[derive(Clone, Debug)]
|
||||
pub enum TransformType {
|
||||
Rigid, // rotation + translation
|
||||
Similarity, // rotation + translation + scale
|
||||
Affine, // full 6 parameters
|
||||
}
|
||||
|
||||
/// A 2D affine transform: [a, b, tx; c, d, ty]
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct AffineTransform {
|
||||
pub a: f64,
|
||||
pub b: f64,
|
||||
pub tx: f64,
|
||||
pub c: f64,
|
||||
pub d: f64,
|
||||
pub ty: f64,
|
||||
}
|
||||
|
||||
impl AffineTransform {
|
||||
pub fn identity() -> Self {
|
||||
Self { a: 1.0, b: 0.0, tx: 0.0, c: 0.0, d: 1.0, ty: 0.0 }
|
||||
}
|
||||
|
||||
pub fn transform_point(&self, x: f64, y: f64) -> (f64, f64) {
|
||||
(self.a * x + self.b * y + self.tx,
|
||||
self.c * x + self.d * y + self.ty)
|
||||
}
|
||||
}
|
||||
|
||||
/// Compute a transform from source→target landmark correspondences.
|
||||
pub fn landmark_registration(
|
||||
source_points: &[(f64, f64)],
|
||||
target_points: &[(f64, f64)],
|
||||
transform_type: TransformType,
|
||||
) -> AffineTransform {
|
||||
let n = source_points.len().min(target_points.len());
|
||||
if n < 2 {
|
||||
return AffineTransform::identity();
|
||||
}
|
||||
|
||||
match transform_type {
|
||||
TransformType::Rigid => fit_rigid(source_points, target_points, n),
|
||||
TransformType::Similarity => fit_similarity(source_points, target_points, n),
|
||||
TransformType::Affine => fit_affine(source_points, target_points, n),
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply an affine transform to an image with bilinear interpolation.
|
||||
pub fn apply_transform(
|
||||
image: &TypedBuffer,
|
||||
width: u32,
|
||||
height: u32,
|
||||
transform: &AffineTransform,
|
||||
) -> TypedBuffer {
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
let n = w * h;
|
||||
let src: Vec<f64> = (0..n).map(|i| image.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
|
||||
// Invert the transform to map output → input
|
||||
let inv = invert_affine(transform);
|
||||
let mut out = vec![0.0f32; n];
|
||||
|
||||
for oy in 0..h {
|
||||
for ox in 0..w {
|
||||
let (sx, sy) = inv.transform_point(ox as f64, oy as f64);
|
||||
out[oy * w + ox] = bilinear_sample(&src, w, h, sx, sy) as f32;
|
||||
}
|
||||
}
|
||||
TypedBuffer::F32(out)
|
||||
}
|
||||
|
||||
/// B-spline non-rigid registration with coarse-to-fine refinement.
|
||||
/// Uses SSD similarity metric.
|
||||
pub fn bspline_registration(
|
||||
source: &TypedBuffer,
|
||||
target: &TypedBuffer,
|
||||
width: u32,
|
||||
height: u32,
|
||||
grid_spacing: u32,
|
||||
iterations: usize,
|
||||
) -> TypedBuffer {
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
let n = w * h;
|
||||
|
||||
let src: Vec<f64> = (0..n).map(|i| source.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let tgt: Vec<f64> = (0..n).map(|i| target.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
|
||||
// Initialize displacement field
|
||||
let mut dx_field = vec![0.0f64; n];
|
||||
let mut dy_field = vec![0.0f64; n];
|
||||
|
||||
// Coarse-to-fine: start with large grid, refine
|
||||
let spacings = [grid_spacing * 4, grid_spacing * 2, grid_spacing];
|
||||
|
||||
for &spacing in &spacings {
|
||||
let sp = spacing.max(2) as usize;
|
||||
let gw = (w + sp - 1) / sp + 1;
|
||||
let gh = (h + sp - 1) / sp + 1;
|
||||
let mut ctrl_dx = vec![0.0f64; gw * gh];
|
||||
let mut ctrl_dy = vec![0.0f64; gw * gh];
|
||||
|
||||
let step_size = 0.5;
|
||||
|
||||
for _iter in 0..iterations {
|
||||
// Compute gradient for each control point
|
||||
let mut grad_dx = vec![0.0f64; gw * gh];
|
||||
let mut grad_dy = vec![0.0f64; gw * gh];
|
||||
|
||||
for gy in 0..gh {
|
||||
for gx in 0..gw {
|
||||
let cx = (gx * sp) as f64;
|
||||
let cy = (gy * sp) as f64;
|
||||
|
||||
// Sample region around control point
|
||||
let mut grad_x = 0.0;
|
||||
let mut grad_y = 0.0;
|
||||
let mut count = 0.0;
|
||||
|
||||
let r = sp as i32;
|
||||
for dy in -r..=r {
|
||||
for dxx in -r..=r {
|
||||
let px = cx as i32 + dxx;
|
||||
let py = cy as i32 + dy;
|
||||
if px < 0 || px >= w as i32 || py < 0 || py >= h as i32 {
|
||||
continue;
|
||||
}
|
||||
let pi = py as usize * w + px as usize;
|
||||
let sx = px as f64 + dx_field[pi] + ctrl_dx[gy * gw + gx];
|
||||
let sy = py as f64 + dy_field[pi] + ctrl_dy[gy * gw + gx];
|
||||
let warped = bilinear_sample(&src, w, h, sx, sy);
|
||||
let diff = warped - tgt[pi];
|
||||
|
||||
// Image gradient at warped position
|
||||
let ix = bilinear_sample(&src, w, h, sx + 0.5, sy)
|
||||
- bilinear_sample(&src, w, h, sx - 0.5, sy);
|
||||
let iy = bilinear_sample(&src, w, h, sx, sy + 0.5)
|
||||
- bilinear_sample(&src, w, h, sx, sy - 0.5);
|
||||
|
||||
grad_x += diff * ix;
|
||||
grad_y += diff * iy;
|
||||
count += 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
if count > 0.0 {
|
||||
grad_dx[gy * gw + gx] = grad_x / count;
|
||||
grad_dy[gy * gw + gx] = grad_y / count;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Update control points
|
||||
for i in 0..gw * gh {
|
||||
ctrl_dx[i] -= step_size * grad_dx[i];
|
||||
ctrl_dy[i] -= step_size * grad_dy[i];
|
||||
}
|
||||
}
|
||||
|
||||
// Interpolate control points to displacement field
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let gx = x as f64 / sp as f64;
|
||||
let gy = y as f64 / sp as f64;
|
||||
dx_field[y * w + x] += bspline_interpolate(&ctrl_dx, gw, gh, gx, gy);
|
||||
dy_field[y * w + x] += bspline_interpolate(&ctrl_dy, gw, gh, gx, gy);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Apply displacement field
|
||||
let mut result = vec![0.0f32; n];
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = y * w + x;
|
||||
let sx = x as f64 + dx_field[i];
|
||||
let sy = y as f64 + dy_field[i];
|
||||
result[i] = bilinear_sample(&src, w, h, sx, sy) as f32;
|
||||
}
|
||||
}
|
||||
TypedBuffer::F32(result)
|
||||
}
|
||||
|
||||
/// Thirion's demons algorithm for deformable registration.
|
||||
pub fn demons_registration(
|
||||
source: &TypedBuffer,
|
||||
target: &TypedBuffer,
|
||||
width: u32,
|
||||
height: u32,
|
||||
iterations: usize,
|
||||
sigma: f64,
|
||||
) -> TypedBuffer {
|
||||
let w = width as usize;
|
||||
let h = height as usize;
|
||||
let n = w * h;
|
||||
|
||||
let src: Vec<f64> = (0..n).map(|i| source.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let tgt: Vec<f64> = (0..n).map(|i| target.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
|
||||
let mut dx = vec![0.0f64; n];
|
||||
let mut dy = vec![0.0f64; n];
|
||||
|
||||
for _iter in 0..iterations {
|
||||
let mut update_dx = vec![0.0f64; n];
|
||||
let mut update_dy = vec![0.0f64; n];
|
||||
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = y * w + x;
|
||||
let sx = x as f64 + dx[i];
|
||||
let sy = y as f64 + dy[i];
|
||||
let warped = bilinear_sample(&src, w, h, sx, sy);
|
||||
let diff = warped - tgt[i];
|
||||
|
||||
// Compute gradient of warped image
|
||||
let gx = bilinear_sample(&src, w, h, sx + 0.5, sy)
|
||||
- bilinear_sample(&src, w, h, sx - 0.5, sy);
|
||||
let gy = bilinear_sample(&src, w, h, sx, sy + 0.5)
|
||||
- bilinear_sample(&src, w, h, sx, sy - 0.5);
|
||||
|
||||
let grad_sq = gx * gx + gy * gy;
|
||||
let denom = grad_sq + diff * diff;
|
||||
|
||||
if denom > 1e-10 {
|
||||
update_dx[i] = -diff * gx / denom;
|
||||
update_dy[i] = -diff * gy / denom;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Smooth the update field with Gaussian
|
||||
let smooth_dx = gaussian_smooth_1d(&update_dx, w, h, sigma);
|
||||
let smooth_dy = gaussian_smooth_1d(&update_dy, w, h, sigma);
|
||||
|
||||
// Compose
|
||||
for i in 0..n {
|
||||
dx[i] += smooth_dx[i];
|
||||
dy[i] += smooth_dy[i];
|
||||
}
|
||||
|
||||
// Regularize the total displacement field
|
||||
dx = gaussian_smooth_1d(&dx, w, h, sigma);
|
||||
dy = gaussian_smooth_1d(&dy, w, h, sigma);
|
||||
}
|
||||
|
||||
// Apply final displacement
|
||||
let mut result = vec![0.0f32; n];
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = y * w + x;
|
||||
result[i] = bilinear_sample(&src, w, h, x as f64 + dx[i], y as f64 + dy[i]) as f32;
|
||||
}
|
||||
}
|
||||
TypedBuffer::F32(result)
|
||||
}
|
||||
|
||||
// ---- Helper functions ----
|
||||
|
||||
fn fit_rigid(src: &[(f64, f64)], tgt: &[(f64, f64)], n: usize) -> AffineTransform {
|
||||
// Compute centroids
|
||||
let (mut sx, mut sy, mut tx, mut ty) = (0.0, 0.0, 0.0, 0.0);
|
||||
for i in 0..n {
|
||||
sx += src[i].0; sy += src[i].1;
|
||||
tx += tgt[i].0; ty += tgt[i].1;
|
||||
}
|
||||
let nf = n as f64;
|
||||
let (scx, scy) = (sx / nf, sy / nf);
|
||||
let (tcx, tcy) = (tx / nf, ty / nf);
|
||||
|
||||
// Compute rotation using SVD-like approach
|
||||
let mut num = 0.0;
|
||||
let mut den = 0.0;
|
||||
for i in 0..n {
|
||||
let dx_s = src[i].0 - scx;
|
||||
let dy_s = src[i].1 - scy;
|
||||
let dx_t = tgt[i].0 - tcx;
|
||||
let dy_t = tgt[i].1 - tcy;
|
||||
num += dx_s * dy_t - dy_s * dx_t;
|
||||
den += dx_s * dx_t + dy_s * dy_t;
|
||||
}
|
||||
|
||||
let theta = num.atan2(den);
|
||||
let cos_t = theta.cos();
|
||||
let sin_t = theta.sin();
|
||||
|
||||
AffineTransform {
|
||||
a: cos_t,
|
||||
b: -sin_t,
|
||||
tx: tcx - cos_t * scx + sin_t * scy,
|
||||
c: sin_t,
|
||||
d: cos_t,
|
||||
ty: tcy - sin_t * scx - cos_t * scy,
|
||||
}
|
||||
}
|
||||
|
||||
fn fit_similarity(src: &[(f64, f64)], tgt: &[(f64, f64)], n: usize) -> AffineTransform {
|
||||
let nf = n as f64;
|
||||
let (mut scx, mut scy, mut tcx, mut tcy) = (0.0, 0.0, 0.0, 0.0);
|
||||
for i in 0..n {
|
||||
scx += src[i].0; scy += src[i].1;
|
||||
tcx += tgt[i].0; tcy += tgt[i].1;
|
||||
}
|
||||
scx /= nf; scy /= nf;
|
||||
tcx /= nf; tcy /= nf;
|
||||
|
||||
let mut num = 0.0;
|
||||
let mut den = 0.0;
|
||||
let mut src_var = 0.0;
|
||||
for i in 0..n {
|
||||
let dx_s = src[i].0 - scx;
|
||||
let dy_s = src[i].1 - scy;
|
||||
let dx_t = tgt[i].0 - tcx;
|
||||
let dy_t = tgt[i].1 - tcy;
|
||||
num += dx_s * dy_t - dy_s * dx_t;
|
||||
den += dx_s * dx_t + dy_s * dy_t;
|
||||
src_var += dx_s * dx_s + dy_s * dy_s;
|
||||
}
|
||||
|
||||
let theta = num.atan2(den);
|
||||
let scale = if src_var > 1e-10 {
|
||||
(den * theta.cos() + num * theta.sin()) / src_var
|
||||
} else {
|
||||
1.0
|
||||
};
|
||||
|
||||
let cos_t = scale * theta.cos();
|
||||
let sin_t = scale * theta.sin();
|
||||
|
||||
AffineTransform {
|
||||
a: cos_t,
|
||||
b: -sin_t,
|
||||
tx: tcx - cos_t * scx + sin_t * scy,
|
||||
c: sin_t,
|
||||
d: cos_t,
|
||||
ty: tcy - sin_t * scx - cos_t * scy,
|
||||
}
|
||||
}
|
||||
|
||||
fn fit_affine(src: &[(f64, f64)], tgt: &[(f64, f64)], n: usize) -> AffineTransform {
|
||||
if n < 3 {
|
||||
return fit_similarity(src, tgt, n);
|
||||
}
|
||||
|
||||
// Least squares: solve for [a, b, tx; c, d, ty]
|
||||
// A * [a b tx]^T = [tx_1..tx_n]^T
|
||||
// A = [x_i, y_i, 1]
|
||||
let mut ata = [[0.0f64; 3]; 3];
|
||||
let mut atb_x = [0.0f64; 3];
|
||||
let mut atb_y = [0.0f64; 3];
|
||||
|
||||
for i in 0..n {
|
||||
let row = [src[i].0, src[i].1, 1.0];
|
||||
for j in 0..3 {
|
||||
for k in 0..3 {
|
||||
ata[j][k] += row[j] * row[k];
|
||||
}
|
||||
atb_x[j] += row[j] * tgt[i].0;
|
||||
atb_y[j] += row[j] * tgt[i].1;
|
||||
}
|
||||
}
|
||||
|
||||
let params_x = solve_3x3(&ata, &atb_x);
|
||||
let params_y = solve_3x3(&ata, &atb_y);
|
||||
|
||||
AffineTransform {
|
||||
a: params_x[0],
|
||||
b: params_x[1],
|
||||
tx: params_x[2],
|
||||
c: params_y[0],
|
||||
d: params_y[1],
|
||||
ty: params_y[2],
|
||||
}
|
||||
}
|
||||
|
||||
fn solve_3x3(a: &[[f64; 3]; 3], b: &[f64; 3]) -> [f64; 3] {
|
||||
// Gaussian elimination with partial pivoting
|
||||
let mut m = [
|
||||
[a[0][0], a[0][1], a[0][2], b[0]],
|
||||
[a[1][0], a[1][1], a[1][2], b[1]],
|
||||
[a[2][0], a[2][1], a[2][2], b[2]],
|
||||
];
|
||||
|
||||
for col in 0..3 {
|
||||
// Pivot
|
||||
let mut max_row = col;
|
||||
for row in col + 1..3 {
|
||||
if m[row][col].abs() > m[max_row][col].abs() {
|
||||
max_row = row;
|
||||
}
|
||||
}
|
||||
m.swap(col, max_row);
|
||||
|
||||
let pivot = m[col][col];
|
||||
if pivot.abs() < 1e-15 {
|
||||
continue;
|
||||
}
|
||||
|
||||
for row in col + 1..3 {
|
||||
let factor = m[row][col] / pivot;
|
||||
for j in col..4 {
|
||||
m[row][j] -= factor * m[col][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Back substitution
|
||||
let mut x = [0.0; 3];
|
||||
for i in (0..3).rev() {
|
||||
if m[i][i].abs() < 1e-15 {
|
||||
continue;
|
||||
}
|
||||
x[i] = m[i][3];
|
||||
for j in i + 1..3 {
|
||||
x[i] -= m[i][j] * x[j];
|
||||
}
|
||||
x[i] /= m[i][i];
|
||||
}
|
||||
x
|
||||
}
|
||||
|
||||
fn invert_affine(t: &AffineTransform) -> AffineTransform {
|
||||
let det = t.a * t.d - t.b * t.c;
|
||||
if det.abs() < 1e-15 {
|
||||
return AffineTransform::identity();
|
||||
}
|
||||
let inv_det = 1.0 / det;
|
||||
AffineTransform {
|
||||
a: t.d * inv_det,
|
||||
b: -t.b * inv_det,
|
||||
tx: (t.b * t.ty - t.d * t.tx) * inv_det,
|
||||
c: -t.c * inv_det,
|
||||
d: t.a * inv_det,
|
||||
ty: (t.c * t.tx - t.a * t.ty) * inv_det,
|
||||
}
|
||||
}
|
||||
|
||||
fn bilinear_sample(data: &[f64], w: usize, h: usize, x: f64, y: f64) -> f64 {
|
||||
if x < 0.0 || y < 0.0 || x >= (w - 1) as f64 || y >= (h - 1) as f64 {
|
||||
// Clamp to edge
|
||||
let cx = x.clamp(0.0, (w - 1) as f64);
|
||||
let cy = y.clamp(0.0, (h - 1) as f64);
|
||||
return data[cy.round() as usize * w + cx.round() as usize];
|
||||
}
|
||||
|
||||
let x0 = x.floor() as usize;
|
||||
let y0 = y.floor() as usize;
|
||||
let x1 = x0 + 1;
|
||||
let y1 = y0 + 1;
|
||||
let fx = x - x0 as f64;
|
||||
let fy = y - y0 as f64;
|
||||
|
||||
data[y0 * w + x0] * (1.0 - fx) * (1.0 - fy)
|
||||
+ data[y0 * w + x1] * fx * (1.0 - fy)
|
||||
+ data[y1 * w + x0] * (1.0 - fx) * fy
|
||||
+ data[y1 * w + x1] * fx * fy
|
||||
}
|
||||
|
||||
fn bspline_interpolate(ctrl: &[f64], gw: usize, gh: usize, gx: f64, gy: f64) -> f64 {
|
||||
let ix = gx.floor() as i32;
|
||||
let iy = gy.floor() as i32;
|
||||
let fx = gx - ix as f64;
|
||||
let fy = gy - iy as f64;
|
||||
|
||||
// Bilinear interpolation on control grid
|
||||
let get = |x: i32, y: i32| -> f64 {
|
||||
let cx = x.clamp(0, gw as i32 - 1) as usize;
|
||||
let cy = y.clamp(0, gh as i32 - 1) as usize;
|
||||
ctrl[cy * gw + cx]
|
||||
};
|
||||
|
||||
get(ix, iy) * (1.0 - fx) * (1.0 - fy)
|
||||
+ get(ix + 1, iy) * fx * (1.0 - fy)
|
||||
+ get(ix, iy + 1) * (1.0 - fx) * fy
|
||||
+ get(ix + 1, iy + 1) * fx * fy
|
||||
}
|
||||
|
||||
fn gaussian_smooth_1d(data: &[f64], w: usize, h: usize, sigma: f64) -> Vec<f64> {
|
||||
if sigma < 0.01 {
|
||||
return data.to_vec();
|
||||
}
|
||||
let radius = (sigma * 3.0).ceil() as usize;
|
||||
let size = 2 * radius + 1;
|
||||
let mut kernel = vec![0.0; size];
|
||||
let mut sum = 0.0;
|
||||
for i in 0..size {
|
||||
let x = i as f64 - radius as f64;
|
||||
let v = (-x * x / (2.0 * sigma * sigma)).exp();
|
||||
kernel[i] = v;
|
||||
sum += v;
|
||||
}
|
||||
for v in &mut kernel {
|
||||
*v /= sum;
|
||||
}
|
||||
|
||||
// Horizontal
|
||||
let mut temp = vec![0.0; w * h];
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let mut s = 0.0;
|
||||
for (ki, &kv) in kernel.iter().enumerate() {
|
||||
let sx = (x as isize + ki as isize - radius as isize).clamp(0, w as isize - 1) as usize;
|
||||
s += data[y * w + sx] * kv;
|
||||
}
|
||||
temp[y * w + x] = s;
|
||||
}
|
||||
}
|
||||
|
||||
// Vertical
|
||||
let mut output = vec![0.0; w * h];
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let mut s = 0.0;
|
||||
for (ki, &kv) in kernel.iter().enumerate() {
|
||||
let sy = (y as isize + ki as isize - radius as isize).clamp(0, h as isize - 1) as usize;
|
||||
s += temp[sy * w + x] * kv;
|
||||
}
|
||||
output[y * w + x] = s;
|
||||
}
|
||||
}
|
||||
output
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn make_test_image(w: usize, h: usize, cx: f64, cy: f64) -> TypedBuffer {
|
||||
let mut data = vec![0.0f32; w * h];
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let dx = x as f64 - cx;
|
||||
let dy = y as f64 - cy;
|
||||
data[y * w + x] = (-0.05 * (dx * dx + dy * dy)).exp() as f32;
|
||||
}
|
||||
}
|
||||
TypedBuffer::F32(data)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn landmark_rigid_identity() {
|
||||
let points = vec![(10.0, 10.0), (20.0, 10.0), (10.0, 20.0)];
|
||||
let t = landmark_registration(&points, &points, TransformType::Rigid);
|
||||
assert!((t.a - 1.0).abs() < 1e-6);
|
||||
assert!((t.d - 1.0).abs() < 1e-6);
|
||||
assert!(t.tx.abs() < 1e-6);
|
||||
assert!(t.ty.abs() < 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn landmark_rigid_translation() {
|
||||
let src = vec![(0.0, 0.0), (10.0, 0.0), (0.0, 10.0)];
|
||||
let tgt = vec![(5.0, 3.0), (15.0, 3.0), (5.0, 13.0)];
|
||||
let t = landmark_registration(&src, &tgt, TransformType::Rigid);
|
||||
assert!((t.tx - 5.0).abs() < 1e-3, "tx={}", t.tx);
|
||||
assert!((t.ty - 3.0).abs() < 1e-3, "ty={}", t.ty);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn landmark_similarity_with_scale() {
|
||||
let src = vec![(0.0, 0.0), (10.0, 0.0), (0.0, 10.0)];
|
||||
let tgt = vec![(0.0, 0.0), (20.0, 0.0), (0.0, 20.0)];
|
||||
let t = landmark_registration(&src, &tgt, TransformType::Similarity);
|
||||
// Scale should be ~2
|
||||
let scale = (t.a * t.a + t.c * t.c).sqrt();
|
||||
assert!((scale - 2.0).abs() < 0.1, "scale={}", scale);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn landmark_affine_exact() {
|
||||
let src = vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)];
|
||||
let tgt = vec![(1.0, 2.0), (3.0, 2.0), (1.0, 5.0)];
|
||||
// Expected: x' = 2x + 0y + 1, y' = 0x + 3y + 2
|
||||
let t = landmark_registration(&src, &tgt, TransformType::Affine);
|
||||
assert!((t.a - 2.0).abs() < 1e-3, "a={}", t.a);
|
||||
assert!((t.tx - 1.0).abs() < 1e-3, "tx={}", t.tx);
|
||||
assert!((t.d - 3.0).abs() < 1e-3, "d={}", t.d);
|
||||
assert!((t.ty - 2.0).abs() < 1e-3, "ty={}", t.ty);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn apply_transform_identity() {
|
||||
let img = make_test_image(16, 16, 8.0, 8.0);
|
||||
let t = AffineTransform::identity();
|
||||
let result = apply_transform(&img, 16, 16, &t);
|
||||
let src = img.as_f32_slice().unwrap();
|
||||
let dst = result.as_f32_slice().unwrap();
|
||||
for i in 0..src.len() {
|
||||
assert!((src[i] - dst[i]).abs() < 1e-3);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn apply_transform_translation() {
|
||||
let w = 32u32;
|
||||
let h = 32u32;
|
||||
let img = make_test_image(w as usize, h as usize, 16.0, 16.0);
|
||||
let t = AffineTransform {
|
||||
a: 1.0, b: 0.0, tx: 3.0,
|
||||
c: 0.0, d: 1.0, ty: 2.0,
|
||||
};
|
||||
let result = apply_transform(&img, w, h, &t);
|
||||
let dst = result.as_f32_slice().unwrap();
|
||||
let src = img.as_f32_slice().unwrap();
|
||||
// Center of result should be at (19, 18)
|
||||
let idx_src = 16 * w as usize + 16;
|
||||
let idx_dst = 18 * w as usize + 19;
|
||||
assert!((src[idx_src] - dst[idx_dst]).abs() < 0.1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn bspline_registration_reduces_error() {
|
||||
let w = 32u32;
|
||||
let h = 32u32;
|
||||
let source = make_test_image(w as usize, h as usize, 14.0, 14.0);
|
||||
let target = make_test_image(w as usize, h as usize, 16.0, 16.0);
|
||||
|
||||
let result = bspline_registration(&source, &target, w, h, 8, 5);
|
||||
|
||||
// Compute SSD before and after
|
||||
let n = (w * h) as usize;
|
||||
let src: Vec<f64> = (0..n).map(|i| source.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let tgt: Vec<f64> = (0..n).map(|i| target.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let res: Vec<f64> = (0..n).map(|i| result.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
|
||||
let ssd_before: f64 = src.iter().zip(&tgt).map(|(s, t)| (s - t).powi(2)).sum();
|
||||
let ssd_after: f64 = res.iter().zip(&tgt).map(|(r, t)| (r - t).powi(2)).sum();
|
||||
|
||||
assert!(ssd_after < ssd_before, "bspline should reduce error: before={}, after={}", ssd_before, ssd_after);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn demons_registration_reduces_error() {
|
||||
let w = 32u32;
|
||||
let h = 32u32;
|
||||
let source = make_test_image(w as usize, h as usize, 14.0, 14.0);
|
||||
let target = make_test_image(w as usize, h as usize, 16.0, 16.0);
|
||||
|
||||
let result = demons_registration(&source, &target, w, h, 20, 1.5);
|
||||
|
||||
let n = (w * h) as usize;
|
||||
let src: Vec<f64> = (0..n).map(|i| source.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let tgt: Vec<f64> = (0..n).map(|i| target.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
let res: Vec<f64> = (0..n).map(|i| result.get_as_f64(i).unwrap_or(0.0)).collect();
|
||||
|
||||
let ssd_before: f64 = src.iter().zip(&tgt).map(|(s, t)| (s - t).powi(2)).sum();
|
||||
let ssd_after: f64 = res.iter().zip(&tgt).map(|(r, t)| (r - t).powi(2)).sum();
|
||||
|
||||
assert!(ssd_after < ssd_before, "demons should reduce error: before={}, after={}", ssd_before, ssd_after);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn invert_affine_roundtrip() {
|
||||
let t = AffineTransform { a: 0.866, b: -0.5, tx: 10.0, c: 0.5, d: 0.866, ty: 5.0 };
|
||||
let inv = invert_affine(&t);
|
||||
let (x, y) = t.transform_point(3.0, 7.0);
|
||||
let (rx, ry) = inv.transform_point(x, y);
|
||||
assert!((rx - 3.0).abs() < 1e-6);
|
||||
assert!((ry - 7.0).abs() < 1e-6);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user