use crate::error::{NdimageError, NdimageResult};
use crate::filters::{
convolve, gaussian_filter_f32, gradient_magnitude, prewitt, scharr, sobel, BorderMode,
};
use scirs2_core::ndarray::{Array, Array2, ArrayD, Ix2};
use std::f32::consts::PI;
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GradientMethod {
Sobel,
Prewitt,
Scharr,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum EdgeDetectionAlgorithm {
Canny,
LoG,
Gradient,
}
#[derive(Debug, Clone)]
pub struct EdgeDetectionConfig {
pub algorithm: EdgeDetectionAlgorithm,
pub gradient_method: GradientMethod,
pub sigma: f32,
pub low_threshold: f32,
pub high_threshold: f32,
pub border_mode: BorderMode,
pub return_magnitude: bool,
}
impl Default for EdgeDetectionConfig {
fn default() -> Self {
Self {
algorithm: EdgeDetectionAlgorithm::Canny,
gradient_method: GradientMethod::Sobel,
sigma: 1.0,
low_threshold: 0.1,
high_threshold: 0.2,
border_mode: BorderMode::Reflect,
return_magnitude: false,
}
}
}
#[allow(dead_code)]
pub fn edge_detector(
image: &Array<f32, Ix2>,
config: EdgeDetectionConfig,
) -> NdimageResult<Array<f32, Ix2>> {
match config.algorithm {
EdgeDetectionAlgorithm::Canny => {
Ok(canny_impl(
image,
config.sigma,
config.low_threshold,
config.high_threshold,
config.gradient_method,
config.border_mode,
))
}
EdgeDetectionAlgorithm::LoG => {
let edges = laplacian_edges_impl(
image,
config.sigma,
config.low_threshold,
config.border_mode,
)?;
if !config.return_magnitude {
Ok(edges
.mapv(|v| v.abs() > config.low_threshold)
.mapv(|v| if v { 1.0 } else { 0.0 }))
} else {
Ok(edges)
}
}
EdgeDetectionAlgorithm::Gradient => {
let edges = gradient_edges_impl(
image,
config.gradient_method,
config.sigma,
config.border_mode,
)?;
if !config.return_magnitude {
Ok(edges
.mapv(|v| v > config.low_threshold)
.mapv(|v| if v { 1.0 } else { 0.0 }))
} else {
Ok(edges)
}
}
}
}
#[allow(dead_code)]
pub fn canny(
image: &Array<f32, Ix2>,
sigma: f32,
low_threshold: f32,
high_threshold: f32,
method: Option<GradientMethod>,
) -> Array<f32, Ix2> {
let method = method.unwrap_or(GradientMethod::Sobel);
let border_mode = BorderMode::Reflect;
canny_impl(
image,
sigma,
low_threshold,
high_threshold,
method,
border_mode,
)
}
#[allow(dead_code)]
fn canny_impl(
image: &Array<f32, Ix2>,
sigma: f32,
low_threshold: f32,
high_threshold: f32,
method: GradientMethod,
mode: BorderMode,
) -> Array<f32, Ix2> {
let image_dim = image.raw_dim();
let image_d = image.clone().into_dyn();
let smoothed =
gaussian_filter_f32(&image_d, sigma, Some(mode), None).unwrap_or_else(|_| image_d.clone());
let gradients = calculate_gradient(&smoothed, method, mode);
let gradient_x = &gradients.0;
let gradient_y = &gradients.1;
let (magnitude, direction) =
calculate_magnitude_and_direction(gradient_x, gradient_y, image_dim);
let suppressed = non_maximum_suppression(&magnitude, &direction);
hysteresis_thresholding(&suppressed, low_threshold, high_threshold)
}
#[allow(dead_code)]
fn calculate_gradient(
image: &ArrayD<f32>,
method: GradientMethod,
mode: BorderMode,
) -> (ArrayD<f32>, ArrayD<f32>) {
match method {
GradientMethod::Sobel => {
let gy = sobel(image, 0, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
let gx = sobel(image, 1, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
(gx, gy)
}
GradientMethod::Prewitt => {
let gy =
prewitt(image, 0, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
let gx =
prewitt(image, 1, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
(gx, gy)
}
GradientMethod::Scharr => {
let gy =
scharr(image, 0, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
let gx =
scharr(image, 1, Some(mode)).unwrap_or_else(|_| ArrayD::zeros(image.raw_dim()));
(gx, gy)
}
}
}
#[allow(dead_code)]
fn calculate_magnitude_and_direction(
gradient_x: &ArrayD<f32>,
gradient_y: &ArrayD<f32>,
shape: Ix2,
) -> (Array<f32, Ix2>, Array<f32, Ix2>) {
let magnitude = Array::<f32, Ix2>::zeros(shape);
let mut direction = Array::<f32, Ix2>::zeros(shape);
let mut mag_copy = Array::<f32, Ix2>::zeros(shape);
for (pos, _) in magnitude.indexed_iter() {
let idx_d = [pos.0, pos.1];
let gx = gradient_x[idx_d.as_ref()];
let gy = gradient_y[idx_d.as_ref()];
mag_copy[pos] = (gx * gx + gy * gy).sqrt();
let angle = gy.atan2(gx) * 180.0 / PI;
let angle = if angle < 0.0 { angle + 180.0 } else { angle };
direction[pos] = if !(22.5..157.5).contains(&angle) {
0.0 } else if (22.5..67.5).contains(&angle) {
45.0 } else if (67.5..112.5).contains(&angle) {
90.0 } else {
135.0 };
}
(mag_copy, direction)
}
#[allow(dead_code)]
fn non_maximum_suppression(
magnitude: &Array<f32, Ix2>,
direction: &Array<f32, Ix2>,
) -> Array<f32, Ix2> {
let shape = magnitude.dim();
let mut suppressed = Array::zeros(shape);
for row in 1..(shape.0 - 1) {
for col in 1..(shape.1 - 1) {
let dir = direction[(row, col)];
let mag = magnitude[(row, col)];
if mag == 0.0 {
continue;
}
let (neighbor1, neighbor2) = get_gradient_neighbors(row, col, dir, magnitude);
if mag >= neighbor1 && mag >= neighbor2 {
suppressed[(row, col)] = mag;
}
}
}
suppressed
}
#[allow(dead_code)]
fn hysteresis_thresholding(
suppressed: &Array<f32, Ix2>,
low_threshold: f32,
high_threshold: f32,
) -> Array<f32, Ix2> {
let shape = suppressed.dim();
let mut result = Array::from_elem(shape, 0.0);
let mut candidates = Vec::new();
for ((row, col), &val) in suppressed.indexed_iter() {
if val >= high_threshold {
result[(row, col)] = 1.0;
candidates.push((row, col));
} else if val >= low_threshold {
candidates.push((row, col));
}
}
let mut processed = Array::from_elem(shape, false);
for (row, col) in candidates.iter() {
if result[(*row, *col)] > 0.0 || processed[(*row, *col)] {
continue;
}
processed[(*row, *col)] = true;
if is_connected_to_strong_edge(*row, *col, &result) {
result[(*row, *col)] = 1.0;
let mut queue = Vec::with_capacity(8); queue.push((*row, *col));
while let Some((r, c)) = queue.pop() {
for nr in (r.saturating_sub(1))..=(r + 1).min(shape.0 - 1) {
for nc in (c.saturating_sub(1))..=(c + 1).min(shape.1 - 1) {
if processed[(nr, nc)] || result[(nr, nc)] > 0.0 {
continue;
}
if suppressed[(nr, nc)] >= low_threshold {
result[(nr, nc)] = 1.0;
processed[(nr, nc)] = true;
queue.push((nr, nc));
}
}
}
}
}
}
result
}
#[allow(dead_code)]
fn get_gradient_neighbors(
row: usize,
col: usize,
direction: f32,
magnitude: &Array<f32, Ix2>,
) -> (f32, f32) {
if direction == 0.0 {
(magnitude[(row, col - 1)], magnitude[(row, col + 1)])
}
else if direction == 45.0 {
(magnitude[(row - 1, col + 1)], magnitude[(row + 1, col - 1)])
}
else if direction == 90.0 {
(magnitude[(row - 1, col)], magnitude[(row + 1, col)])
}
else {
(magnitude[(row - 1, col - 1)], magnitude[(row + 1, col + 1)])
}
}
#[allow(dead_code)]
fn is_connected_to_strong_edge(row: usize, col: usize, edges: &Array<f32, Ix2>) -> bool {
let shape = edges.dim();
for i in (row.saturating_sub(1))..=(row + 1).min(shape.0 - 1) {
for j in (col.saturating_sub(1))..=(col + 1).min(shape.1 - 1) {
if !(i == row && j == col) && edges[(i, j)] > 0.0 {
return true;
}
}
}
false
}
#[allow(dead_code)]
pub fn laplacian_edges(
image: &Array<f32, Ix2>,
sigma: f32,
threshold: Option<f32>,
mode: Option<BorderMode>,
) -> NdimageResult<Array<f32, Ix2>> {
let mode = mode.unwrap_or(BorderMode::Reflect);
Ok(laplacian_edges_impl(
image,
sigma,
threshold.unwrap_or(0.0),
mode,
)?)
}
#[allow(dead_code)]
fn laplacian_edges_impl(
image: &Array<f32, Ix2>,
sigma: f32,
threshold: f32,
mode: BorderMode,
) -> NdimageResult<Array<f32, Ix2>> {
let image_d = image.clone().into_dyn();
let smoothed =
gaussian_filter_f32(&image_d, sigma, Some(mode), None).unwrap_or_else(|_| image_d.clone());
let ksize = ((2.0 * (3.0 * sigma).ceil() + 1.0).max(3.0)) as usize;
let mut kernel = ArrayD::zeros(vec![ksize; image.ndim()]);
let center = ksize / 2;
let center_value = -2.0 * image.ndim() as f32;
let center_idx = vec![center; image.ndim()];
kernel[center_idx.as_slice()] = center_value;
for dim in 0..image.ndim() {
let mut idx = vec![center; image.ndim()];
if center > 0 {
idx[dim] = center - 1;
kernel[idx.as_slice()] = 1.0;
}
idx[dim] = center + 1;
if idx[dim] < ksize {
kernel[idx.as_slice()] = 1.0;
}
}
let laplacian = convolve(&smoothed, &kernel, Some(mode)).map_err(|e| {
NdimageError::ComputationError(format!("Laplacian convolution failed: {}", e))
})?;
let mut result_copy = Array::zeros(image.dim());
for i in 0..image.dim().0 {
for j in 0..image.dim().1 {
result_copy[(i, j)] = laplacian[[i, j]];
}
}
if threshold > 0.0 {
return Ok(result_copy.mapv(|v| if v.abs() > threshold { v } else { 0.0 }));
}
Ok(result_copy)
}
#[allow(dead_code)]
pub fn gradient_edges(
image: &Array<f32, Ix2>,
method: Option<GradientMethod>,
sigma: Option<f32>,
mode: Option<BorderMode>,
) -> NdimageResult<Array<f32, Ix2>> {
let method = method.unwrap_or(GradientMethod::Sobel);
let mode = mode.unwrap_or(BorderMode::Reflect);
Ok(gradient_edges_impl(
image,
method,
sigma.unwrap_or(0.0),
mode,
)?)
}
#[allow(dead_code)]
fn gradient_edges_impl(
image: &Array<f32, Ix2>,
method: GradientMethod,
sigma: f32,
mode: BorderMode,
) -> NdimageResult<Array<f32, Ix2>> {
let image_d = image.clone().into_dyn();
let processed = if sigma > 0.0 {
gaussian_filter_f32(&image_d, sigma, Some(mode), None).map_err(|e| {
NdimageError::ComputationError(format!("Gaussian smoothing failed: {}", e))
})?
} else {
image_d
};
let method_str = match method {
GradientMethod::Sobel => "sobel",
GradientMethod::Prewitt => "prewitt",
GradientMethod::Scharr => "scharr",
};
let magnitude = gradient_magnitude(&processed, Some(mode), Some(method_str)).map_err(|e| {
NdimageError::ComputationError(format!("Gradient magnitude calculation failed: {}", e))
})?;
let mut result_copy = Array::zeros(image.dim());
for i in 0..image.dim().0 {
for j in 0..image.dim().1 {
result_copy[(i, j)] = magnitude[[i, j]];
}
}
Ok(result_copy)
}
#[allow(dead_code)]
pub fn sobel_edges(image: &ArrayD<f32>) -> NdimageResult<ArrayD<f32>> {
edge_detector_simple(image, Some(GradientMethod::Sobel), None)
}
#[allow(dead_code)]
pub fn edge_detector_simple(
image: &ArrayD<f32>,
method: Option<GradientMethod>,
mode: Option<BorderMode>,
) -> NdimageResult<ArrayD<f32>> {
let border_mode = mode.unwrap_or(BorderMode::Reflect);
let method_str = match method.unwrap_or(GradientMethod::Sobel) {
GradientMethod::Sobel => "sobel",
GradientMethod::Prewitt => "prewitt",
GradientMethod::Scharr => "scharr",
};
gradient_magnitude(image, Some(border_mode), Some(method_str))
}
#[cfg(test)]
mod tests {
use super::*;
use scirs2_core::ndarray::array;
#[test]
fn test_canny_edge_detector() {
let image = array![
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
];
let edges_sobel = canny(&image, 1.0, 0.1, 0.2, None);
assert!(
edges_sobel.fold(false, |acc, &x| acc || x > 0.0),
"No edges detected with Sobel"
);
let edges_scharr = canny(&image, 1.0, 0.1, 0.2, Some(GradientMethod::Scharr));
assert!(
edges_scharr.fold(false, |acc, &x| acc || x > 0.0),
"No edges detected with Scharr"
);
}
#[test]
fn test_unified_edge_detector() {
let image = array![
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
];
let edges_default = edge_detector(&image, EdgeDetectionConfig::default())
.expect("edge_detector should succeed for test with default config");
assert!(
edges_default.fold(false, |acc, &x| acc || (x > 0.0)),
"Edges should be detected with default config"
);
let gradient_config = EdgeDetectionConfig {
algorithm: EdgeDetectionAlgorithm::Gradient,
gradient_method: GradientMethod::Sobel,
sigma: 1.0,
low_threshold: 0.1,
high_threshold: 0.2,
border_mode: BorderMode::Reflect,
return_magnitude: true,
};
let gradient_edges = edge_detector(&image, gradient_config)
.expect("edge_detector should succeed for test with gradient config");
let max_magnitude = gradient_edges.fold(0.0, |acc, &x| if x > acc { x } else { acc });
assert!(
max_magnitude > 0.1,
"Gradient magnitudes should be above threshold"
);
let log_config = EdgeDetectionConfig {
algorithm: EdgeDetectionAlgorithm::LoG,
sigma: 1.0,
low_threshold: 0.05,
..EdgeDetectionConfig::default()
};
let log_edges = edge_detector(&image, log_config)
.expect("edge_detector should succeed for test with LoG config");
assert!(
log_edges.fold(false, |acc, &x| acc || (x > 0.0)),
"LoG should detect edges"
);
}
#[test]
fn test_gradient_edges() {
let image = array![
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 1.0, 1.0, 1.0],
];
let edges_sobel = gradient_edges(&image, Some(GradientMethod::Sobel), Some(1.0), None)
.expect("gradient_edges with Sobel should succeed for test");
let edges_prewitt = gradient_edges(&image, Some(GradientMethod::Prewitt), Some(1.0), None)
.expect("gradient_edges with Prewitt should succeed for test");
let edges_scharr = gradient_edges(&image, Some(GradientMethod::Scharr), Some(1.0), None)
.expect("gradient_edges with Scharr should succeed for test");
assert!(
edges_sobel.iter().any(|&x| x > 0.1),
"Sobel should detect edges"
);
assert!(
edges_prewitt.iter().any(|&x| x > 0.1),
"Prewitt should detect edges"
);
assert!(
edges_scharr.iter().any(|&x| x > 0.1),
"Scharr should detect edges"
);
let diagonalimage = array![
[1.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0],
];
let diag_sobel = gradient_edges(&diagonalimage, Some(GradientMethod::Sobel), None, None)
.expect("gradient_edges with Sobel should succeed for diagonal test");
let diag_scharr = gradient_edges(&diagonalimage, Some(GradientMethod::Scharr), None, None)
.expect("gradient_edges with Scharr should succeed for diagonal test");
let max_sobel = diag_sobel
.iter()
.fold(0.0, |acc, &x| if x > acc { x } else { acc });
let max_scharr = diag_scharr
.iter()
.fold(0.0, |acc, &x| if x > acc { x } else { acc });
assert!(
max_scharr > max_sobel,
"Scharr should have stronger diagonal response"
);
}
#[test]
fn test_laplacian_edges() {
let mut image = Array2::<f32>::zeros((5, 5));
image[[2, 2]] = 1.0;
let edges = laplacian_edges(&image, 1.0, None, None)
.expect("laplacian_edges should succeed for test");
assert!(edges.iter().any(|&x| x != 0.0), "LoG should detect edges");
assert!(
edges[[2, 2]] < 0.0,
"Center of point should have negative LoG value"
);
let edges_threshold = laplacian_edges(&image, 1.0, Some(0.1), None)
.expect("laplacian_edges with threshold should succeed for test");
let count_before = edges.iter().filter(|&&x| x != 0.0).count();
let count_after = edges_threshold.iter().filter(|&&x| x != 0.0).count();
assert!(
count_after <= count_before,
"Thresholding should reduce edge points"
);
}
}