///|
/// Represents a detected feature corner point with sub-pixel position and response score.
pub(all) struct CornerPoint {
pt : Point2D
response : Double
} derive(Debug, Eq)
///|
/// Creates a new corner point.
pub fn CornerPoint::new(pt : Point2D, response : Double) -> CornerPoint {
{ pt, response }
}
///|
/// Harris Corner Detector: Computes corner response map and extracts local maxima.
pub fn detect_harris_corners(
img : Image,
k_param? : Double = 0.04,
threshold? : Double = 1000.0,
nms_radius? : Int = 3,
) -> Array[CornerPoint] {
let w = img.width
let h = img.height
if w < 5 || h < 5 {
return []
}
let (gx, gy) = filter_sobel(img)
let size = w * h
// Compute products of derivatives
let ixx = Array::make(size, 0.0)
let iyy = Array::make(size, 0.0)
let ixy = Array::make(size, 0.0)
for i in 0.. Array[Double] {
let temp = Array::make(size, 0.0)
for y in 0..= w { w - 1 } else { ix }
sum += buf[y * w + cx] * kernel[k]
}
temp[y * w + x] = sum
}
}
let res = Array::make(size, 0.0)
for y in 0..= h { h - 1 } else { iy }
sum += temp[cy * w + x] * kernel[k]
}
res[y * w + x] = sum
}
}
res
}
let s_ixx = smooth_buffer(ixx)
let s_iyy = smooth_buffer(iyy)
let s_ixy = smooth_buffer(ixy)
// Compute Harris response: R = det(M) - k * trace(M)^2
let response_map = Array::make(size, 0.0)
for i in 0.. threshold {
let mut is_max = true
for dy in -r..<=r {
for dx in -r..<=r {
if dx != 0 || dy != 0 {
if response_map[(y + dy) * w + (x + dx)] >= val {
is_max = false
break
}
}
}
if !is_max {
break
}
}
if is_max {
corners.push(
CornerPoint::new(Point2D::new(x.to_double(), y.to_double()), val),
)
}
}
}
}
corners
}
///|
/// Shi-Tomasi Corner Detector (Good Features to Track): Uses minimum eigenvalue of M.
pub fn detect_shi_tomasi_corners(
img : Image,
threshold? : Double = 50.0,
nms_radius? : Int = 3,
) -> Array[CornerPoint] {
let w = img.width
let h = img.height
if w < 5 || h < 5 {
return []
}
let (gx, gy) = filter_sobel(img)
let size = w * h
let ixx = Array::make(size, 0.0)
let iyy = Array::make(size, 0.0)
let ixy = Array::make(size, 0.0)
for i in 0.. 0.0 {
0.5 * (trace - disc.sqrt())
} else {
0.5 * trace
}
response_map[y * w + x] = min_eigen
}
}
// Non-maximum suppression
let corners : Array[CornerPoint] = []
let r = nms_radius
for y in r..<(h - r) {
for x in r..<(w - r) {
let val = response_map[y * w + x]
if val > threshold {
let mut is_max = true
for dy in -r..<=r {
for dx in -r..<=r {
if dx != 0 || dy != 0 {
if response_map[(y + dy) * w + (x + dx)] >= val {
is_max = false
break
}
}
}
if !is_max {
break
}
}
if is_max {
corners.push(
CornerPoint::new(Point2D::new(x.to_double(), y.to_double()), val),
)
}
}
}
}
corners
}
///|
/// FAST-9 Corner Detector (Features from Accelerated Segment Test).
pub fn detect_fast_corners(
img : Image,
threshold? : Int = 20,
nonmax_suppression? : Bool = true,
) -> Array[Point2D] {
let w = img.width
let h = img.height
if w < 7 || h < 7 {
return []
}
// 16-pixel Bresenham circle offsets relative to center (dx, dy)
let circle_offsets : Array[(Int, Int)] = [
(0, -3),
(1, -3),
(2, -2),
(3, -1),
(3, 0),
(3, 1),
(2, 2),
(1, 3),
(0, 3),
(-1, 3),
(-2, 2),
(-3, 1),
(-3, 0),
(-3, -1),
(-2, -2),
(-1, -3),
]
let scores = Array::make(w * h, 0)
let candidates : Array[(Int, Int)] = []
for y in 3..<(h - 3) {
for x in 3..<(w - 3) {
let ip = img.get(x, y)
let t = threshold
let high_t = ip + t
let low_t = ip - t
// Fast check on 4 cardinal pixels (1, 5, 9, 13)
let p1 = img.get(x, y - 3)
let p5 = img.get(x + 3, y)
let p9 = img.get(x, y + 3)
let p13 = img.get(x - 3, y)
let mut count_bright = 0
let mut count_dark = 0
if p1 > high_t {
count_bright += 1
} else if p1 < low_t {
count_dark += 1
}
if p5 > high_t {
count_bright += 1
} else if p5 < low_t {
count_dark += 1
}
if p9 > high_t {
count_bright += 1
} else if p9 < low_t {
count_dark += 1
}
if p13 > high_t {
count_bright += 1
} else if p13 < low_t {
count_dark += 1
}
if count_bright >= 3 || count_dark >= 3 {
// Full 16-point circle test
let vals = Array::make(16, 0)
for k in 0..<16 {
let off = circle_offsets[k]
vals[k] = img.get(x + off.0, y + off.1)
}
// Test for at least 9 contiguous pixels brighter or darker
let mut is_corner = false
let mut score = 0
// Check bright arc
for start in 0..<16 {
let mut contiguous = 0
for step in 0..<9 {
let idx = (start + step) % 16
if vals[idx] > high_t {
contiguous += 1
} else {
break
}
}
if contiguous >= 9 {
is_corner = true
break
}
}
// Check dark arc
if !is_corner {
for start in 0..<16 {
let mut contiguous = 0
for step in 0..<9 {
let idx = (start + step) % 16
if vals[idx] < low_t {
contiguous += 1
} else {
break
}
}
if contiguous >= 9 {
is_corner = true
break
}
}
}
if is_corner {
// Compute corner response score: sum of absolute differences
let mut sum_diff = 0
for k in 0..<16 {
let diff = (vals[k] - ip).abs()
if diff > t {
sum_diff += diff - t
}
}
score = sum_diff
scores[y * w + x] = score
candidates.push((x, y))
}
}
}
}
if !nonmax_suppression {
let res : Array[Point2D] = []
for c in candidates {
res.push(Point2D::new(c.0.to_double(), c.1.to_double()))
}
return res
}
// Non-maximum suppression on 3x3 neighborhood
let res : Array[Point2D] = []
for c in candidates {
let x = c.0
let y = c.1
let val = scores[y * w + x]
let mut is_max = true
for dy in -1..<=1 {
for dx in -1..<=1 {
if dx != 0 || dy != 0 {
if scores[(y + dy) * w + (x + dx)] > val {
is_max = false
break
}
}
}
if !is_max {
break
}
}
if is_max {
res.push(Point2D::new(x.to_double(), y.to_double()))
}
}
res
}
///|
/// Subpixel Corner Refinement using Inverted Hessian and Gradient Orthogonality.
/// Iteratively solves: sum (nabla I(p_i) * nabla I(p_i)^T) * (q - p_i) = 0 for refined corner q.
pub fn refine_corner_subpixel(
img : Image,
initial_corner : Point2D,
window_radius? : Int = 3,
max_iterations? : Int = 20,
epsilon? : Double = 0.01,
) -> Point2D {
let w = img.width
let h = img.height
let mut q = initial_corner
let r = window_radius
for _iter in 0..= w - 1 ||
qy_floor - r < 1 ||
qy_floor + r >= h - 1 {
break
}
for dy in -r..<=r {
for dx in -r..<=r {
let px = qx_floor + dx
let py = qy_floor + dy
let pt = Point2D::new(px.to_double(), py.to_double())
// Gaussian spatial weight
let dist_sq = (dx * dx + dy * dy).to_double()
let weight = @math.exp(
-dist_sq / (2.0 * (r.to_double() * 0.5) * (r.to_double() * 0.5)),
)
// Central difference gradient
let gx = (img.get(px + 1, py) - img.get(px - 1, py)).to_double() * 0.5
let gy = (img.get(px, py + 1) - img.get(px, py - 1)).to_double() * 0.5
let gxx = gx * gx * weight
let gyy = gy * gy * weight
let gxy = gx * gy * weight
sum_a += gxx
sum_b += gyy
sum_c += gxy
sum_bx += gxx * pt.x + gxy * pt.y
sum_by += gxy * pt.x + gyy * pt.y
}
}
let det = sum_a * sum_b - sum_c * sum_c
if det.abs() < 1.0e-12 {
break
}
let inv_det = 1.0 / det
let new_qx = (sum_b * sum_bx - sum_c * sum_by) * inv_det
let new_qy = (sum_a * sum_by - sum_c * sum_bx) * inv_det
let delta_x = new_qx - q.x
let delta_y = new_qy - q.y
q = Point2D::new(new_qx, new_qy)
if delta_x * delta_x + delta_y * delta_y < epsilon * epsilon {
break
}
}
q
}
///|
/// Refines all 4 corners of a quadrilateral to sub-pixel accuracy.
pub fn refine_quad_corners(
img : Image,
quad : Quad,
window_radius? : Int = 3,
) -> Quad {
if quad.corners.length() != 4 {
return quad
}
let c0 = refine_corner_subpixel(img, quad.corners[0], window_radius~)
let c1 = refine_corner_subpixel(img, quad.corners[1], window_radius~)
let c2 = refine_corner_subpixel(img, quad.corners[2], window_radius~)
let c3 = refine_corner_subpixel(img, quad.corners[3], window_radius~)
Quad::new(c0, c1, c2, c3)
}