///|
/// Computes Otsu's optimal global binarization threshold based on maximum inter-class variance.
pub fn threshold_otsu(img : Image) -> Int {
  let hist = Array::make(256, 0)
  let total_pixels = img.width * img.height
  if total_pixels == 0 {
    return 128
  }

  // Build histogram
  for y in 0.. 255 { 255 } else { p }
      hist[idx] = hist[idx] + 1
    }
  }

  // Compute total mean
  let mut sum_total = 0.0
  for i in 0..<256 {
    sum_total += i.to_double() * hist[i].to_double()
  }

  let mut sum_b = 0.0
  let mut w_b = 0
  let mut max_variance = 0.0
  let mut best_threshold = 128

  for t in 0..<256 {
    w_b += hist[t]
    if w_b == 0 {
      continue
    }
    let w_f = total_pixels - w_b
    if w_f == 0 {
      break
    }

    sum_b += t.to_double() * hist[t].to_double()
    let m_b = sum_b / w_b.to_double()
    let m_f = (sum_total - sum_b) / w_f.to_double()

    let diff = m_b - m_f
    let variance = w_b.to_double() * w_f.to_double() * diff * diff

    if variance > max_variance {
      max_variance = variance
      best_threshold = t
    }
  }

  best_threshold
}

///|
/// Applies a global binary threshold to the image (pixels >= threshold become 255, else 0).
pub fn threshold_global(img : Image, threshold : Int) -> Image {
  let res = Image::new(img.width, img.height)
  for y in 0..= threshold { 255 } else { 0 }
      res.set(x, y, bin_val)
    }
  }
  res
}

///|
/// Fast Bradley adaptive thresholding based on 2D Integral Image in O(1) time per pixel.
pub fn threshold_adaptive(
  img : Image,
  window_size : Int,
  constant_c : Int,
) -> Image {
  let w = img.width
  let h = img.height
  let res = Image::new(w, h)
  if w == 0 || h == 0 {
    return res
  }

  let integral = compute_integral_image(img)
  let s = if window_size % 2 == 0 { window_size + 1 } else { window_size }
  let r = s / 2

  for y in 0..= w { w - 1 } else { x + r }
      let y1 = if y + r >= h { h - 1 } else { y + r }

      let count = (x1 - x0 + 1) * (y1 - y0 + 1)
      let sum_val = query_integral_sum(integral, w, h, x0, y0, x1, y1)

      // Threshold: pixel * count >= sum * (100 - T)/100 or sum - constant * count
      let thresh = sum_val - constant_c * count
      let pixel = img.get(x, y)
      let bin_val = if pixel * count >= thresh { 255 } else { 0 }
      res.set(x, y, bin_val)
    }
  }

  res
}

///|
/// Sauvola local adaptive thresholding for degraded documents and non-uniform illumination.
/// T(x, y) = m(x, y) * (1 + k * (s(x, y) / R - 1))
pub fn threshold_sauvola(
  img : Image,
  window_size : Int,
  k_param : Double,
  dynamic_range_r : Double,
) -> Image {
  let w = img.width
  let h = img.height
  let res = Image::new(w, h)
  if w == 0 || h == 0 {
    return res
  }

  let integral = compute_integral_image(img)
  let integral_sq = compute_integral_image_sq(img)
  let s = if window_size % 2 == 0 { window_size + 1 } else { window_size }
  let r = s / 2
  let stride = w + 1

  for y in 0..= w { w - 1 } else { x + r }
      let y1 = if y + r >= h { h - 1 } else { y + r }

      let count = (x1 - x0 + 1) * (y1 - y0 + 1)
      let sum_val = query_integral_sum(integral, w, h, x0, y0, x1, y1)

      let a = integral_sq[y0 * stride + x0]
      let b = integral_sq[y0 * stride + (x1 + 1)]
      let c = integral_sq[(y1 + 1) * stride + x0]
      let d = integral_sq[(y1 + 1) * stride + (x1 + 1)]
      let sum_sq_val = d - b - c + a

      let mean = sum_val.to_double() / count.to_double()
      let variance = sum_sq_val / count.to_double() - mean * mean
      let std_dev = if variance > 0.0 { variance.sqrt() } else { 0.0 }

      let thresh = mean * (1.0 + k_param * (std_dev / dynamic_range_r - 1.0))
      let pixel = img.get(x, y).to_double()
      let bin_val = if pixel >= thresh { 255 } else { 0 }
      res.set(x, y, bin_val)
    }
  }

  res
}

///|
/// Niblack local adaptive thresholding: T(x, y) = m(x, y) + k * s(x, y).
pub fn threshold_niblack(
  img : Image,
  window_size : Int,
  k_param : Double,
) -> Image {
  let w = img.width
  let h = img.height
  let res = Image::new(w, h)
  if w == 0 || h == 0 {
    return res
  }

  let integral = compute_integral_image(img)
  let integral_sq = compute_integral_image_sq(img)
  let s = if window_size % 2 == 0 { window_size + 1 } else { window_size }
  let r = s / 2
  let stride = w + 1

  for y in 0..= w { w - 1 } else { x + r }
      let y1 = if y + r >= h { h - 1 } else { y + r }

      let count = (x1 - x0 + 1) * (y1 - y0 + 1)
      let sum_val = query_integral_sum(integral, w, h, x0, y0, x1, y1)

      let a = integral_sq[y0 * stride + x0]
      let b = integral_sq[y0 * stride + (x1 + 1)]
      let c = integral_sq[(y1 + 1) * stride + x0]
      let d = integral_sq[(y1 + 1) * stride + (x1 + 1)]
      let sum_sq_val = d - b - c + a

      let mean = sum_val.to_double() / count.to_double()
      let variance = sum_sq_val / count.to_double() - mean * mean
      let std_dev = if variance > 0.0 { variance.sqrt() } else { 0.0 }

      let thresh = mean + k_param * std_dev
      let pixel = img.get(x, y).to_double()
      let bin_val = if pixel >= thresh { 255 } else { 0 }
      res.set(x, y, bin_val)
    }
  }

  res
}