///|
fn pixel_with_border(
  image : GrayImage,
  x : Int,
  y : Int,
  policy : BorderPolicy,
) -> Double {
  if image.inside(x, y) {
    image.at(x, y)
  } else {
    match policy {
      Constant(value~) => value
      Clamp =>
        image.at(
          x.clamp(min=0, max=image.width - 1),
          y.clamp(min=0, max=image.height - 1),
        )
      Reflect => {
        let reflect_index = fn(index : Int, size : Int) {
          if size <= 1 {
            0
          } else {
            let period = size * 2 - 2
            let wrapped = index % period
            if wrapped < 0 {
              wrapped + period
            } else {
              wrapped
            }
          }
        }
        image.at(reflect_index(x, image.width), reflect_index(y, image.height))
      }
    }
  }
}

///|
pub fn normalize_minmax(image : GrayImage) -> GrayImage {
  let result = GrayImage(width=image.width, height=image.height)
  if image.pixels.is_empty() {
    result
  } else {
    let low = minimum_value(image.pixels)
    let high = maximum_value(image.pixels)
    let span = high - low
    for y in 0.. GrayImage {
  let result = GrayImage(width=image.width, height=image.height)
  let average = mean_value(image.pixels)
  let deviation = variance_value(image.pixels).sqrt()
  for y in 0.. DoubleBounds {
  { min: minimum_value(image.pixels), max: maximum_value(image.pixels) }
}

///|
pub fn image_mean(image : GrayImage) -> Double {
  mean_value(image.pixels)
}

///|
pub fn image_variance(image : GrayImage) -> Double {
  variance_value(image.pixels)
}

///|
pub fn crop_image(image : GrayImage, bounds : Bounds) -> GrayImage {
  let width = (bounds.max_x - bounds.min_x + 1).max(0)
  let height = (bounds.max_y - bounds.min_y + 1).max(0)
  let result = GrayImage(width~, height~)
  for y in 0.. GrayImage {
  let margin = padding.max(0)
  let result = GrayImage(
    width=image.width + margin * 2,
    height=image.height + margin * 2,
    fill=value,
  )
  for y in 0.. GrayImage {
  match try_convolve3(image, kernel, policy~) {
    Some(result) => result
    None => image.copy()
  }
}

///|
pub fn try_convolve3(
  image : GrayImage,
  kernel : Array[Double],
  policy? : BorderPolicy = Clamp,
) -> GrayImage? {
  if kernel.length() != 9 {
    None
  } else {
    let result = GrayImage(width=image.width, height=image.height)
    for y in 0.. GrayImage {
  match try_convolve_separable(image, kernel, policy~) {
    Some(result) => result
    None => image.copy()
  }
}

///|
pub fn try_convolve_separable(
  image : GrayImage,
  kernel : Array[Double],
  policy? : BorderPolicy = Clamp,
) -> GrayImage? {
  if kernel.is_empty() || kernel.length() % 2 == 0 {
    None
  } else {
    let radius = kernel.length() / 2
    let horizontal = GrayImage(width=image.width, height=image.height)
    for y in 0.. GrayImage {
  let width = radius.max(0) * 2 + 1
  let weight = 1.0 / Double::from_int(width)
  convolve_separable(image, Array::make(width, weight))
}

///|
pub fn gaussian_blur(
  image : GrayImage,
  radius : Int,
  sigma : Double,
) -> GrayImage {
  let safe_radius = radius.max(0)
  let safe_sigma = sigma.max(0.0001)
  let kernel : Array[Double] = []
  for i in 0..<(safe_radius * 2 + 1) {
    let offset = i - safe_radius
    let value = @math.exp(
      -Double::from_int(offset * offset) / (2.0 * safe_sigma * safe_sigma),
    )
    kernel.push(value)
  }
  convolve_separable(image, normalize_sum(kernel))
}

///|
fn neighborhood_values(
  image : GrayImage,
  x : Int,
  y : Int,
  radius : Int,
) -> Array[Double] {
  let values : Array[Double] = []
  let safe_radius = radius.max(0)
  for dy in -safe_radius..<=safe_radius {
    for dx in -safe_radius..<=safe_radius {
      values.push(image.at(x + dx, y + dy))
    }
  }
  values
}

///|
pub fn median_filter(image : GrayImage, radius~ : Int) -> GrayImage {
  let result = GrayImage(width=image.width, height=image.height)
  for y in 0.. GrayImage {
  convolve3(image, [-1.0, 0.0, 1.0, -2.0, 0.0, 2.0, -1.0, 0.0, 1.0])
}

///|
pub fn sobel_y(image : GrayImage) -> GrayImage {
  convolve3(image, [-1.0, -2.0, -1.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0])
}

///|
pub fn gradient_magnitude(image : GrayImage) -> GrayImage {
  let gx = sobel_x(image)
  let gy = sobel_y(image)
  let result = GrayImage(width=image.width, height=image.height)
  for y in 0.. GrayImage {
  let gx = sobel_x(image)
  let gy = sobel_y(image)
  let result = GrayImage(width=image.width, height=image.height)
  for y in 0.. GrayImage {
  let result = GrayImage(width=image.width, height=image.height)
  for y in 0.. GrayImage {
  let result = GrayImage(width=magnitude.width, height=magnitude.height)
  let safe_radius = radius.max(1)
  for y in 0.. center {
            largest = false
          }
        }
      }
      if largest {
        result.set(x, y, center)
      }
    }
  }
  result
}

///|
pub fn threshold_image(image : GrayImage, threshold : Double) -> BinaryMask {
  let mask = BinaryMask(width=image.width, height=image.height)
  for y in 0..= threshold)
    }
  }
  mask
}

///|
pub fn threshold_by_policy(
  image : GrayImage,
  threshold : PixelThreshold,
) -> BinaryMask {
  let mask = BinaryMask(width=image.width, height=image.height)
  for y in 0.. sample < value
        Above(value~) => sample > value
        Between(lower~, upper~) => sample >= lower && sample <= upper
      }
      mask.set(x, y, selected)
    }
  }
  mask
}

///|
pub fn histogram(
  image : GrayImage,
  bins : Int,
  bounds? : DoubleBounds = { min: 0.0, max: 1.0 },
) -> Histogram {
  let count = bins.max(1)
  let result = Array::make(count, 0)
  let width = bounds.max - bounds.min
  for value in image.pixels {
    let index = if width <= 0.0 {
      0
    } else {
      ((value - bounds.min) / width * Double::from_int(count))
      .floor()
      .to_int()
      .clamp(min=0, max=count - 1)
    }
    result[index] += 1
  }
  { minimum: bounds.min, maximum: bounds.max, bins: result }
}

///|
pub fn otsu_threshold(image : GrayImage) -> Double {
  let data = insertion_sorted(image.pixels)
  if data.is_empty() {
    0.0
  } else {
    let mut best = data[0]
    let mut best_score = -1.0
    for candidate in data {
      let lower = image.pixels.filter(fn(value) { value <= candidate })
      let upper = image.pixels.filter(fn(value) { value > candidate })
      if !lower.is_empty() && !upper.is_empty() {
        let weight = Double::from_int(lower.length() * upper.length())
        let difference = mean_value(lower) - mean_value(upper)
        let score = weight * difference * difference
        if score > best_score {
          best_score = score
          best = candidate
        }
      }
    }
    best
  }
}

///|
pub fn hysteresis_edges(
  image : GrayImage,
  low~ : Double,
  high~ : Double,
) -> BinaryMask {
  let magnitude = gradient_magnitude(image)
  let result = BinaryMask(width=image.width, height=image.height)
  let weak = BinaryMask(width=image.width, height=image.height)
  let strong = BinaryMask(width=image.width, height=image.height)
  for y in 0..= high {
        strong.set(x, y, true)
        result.set(x, y, true)
      } else if value >= low {
        weak.set(x, y, true)
      }
    }
  }
  let mut changed = true
  while changed {
    changed = false
    for y in 0.. Double {
  let mut checksum = 0.0
  for i in 0.. GrayImage {
  let result = GrayImage(width=image.width, height=image.height)
  for i in 0..