///|
/// Represents an 8-bit grayscale 2D image.
pub(all) struct Image {
  width : Int
  height : Int
  data : Array[Int]
} derive(Debug, Eq)

///|
/// Creates a new Image with width and height, filled with zeros.
pub fn Image::new(width : Int, height : Int) -> Image {
  let safe_width = if width > 0 { width } else { 0 }
  let safe_height = if height > 0 { height } else { 0 }
  let data = Array::make(safe_width * safe_height, 0)
  { width: safe_width, height: safe_height, data }
}

///|
/// Creates an Image from a flat array of pixels.
pub fn Image::from_array(width : Int, height : Int, data : Array[Int]) -> Image {
  let size = width * height
  let img_data = if data.length() == size {
    data.copy()
  } else {
    let d = Array::make(size, 0)
    let min_len = if data.length() < size { data.length() } else { size }
    for i in 0.. Int {
  if x < 0 || x >= self.width || y < 0 || y >= self.height {
    0
  } else {
    self.data[y * self.width + x]
  }
}

///|
/// Sets pixel value at (x, y).
pub fn Image::set(self : Image, x : Int, y : Int, val : Int) -> Unit {
  if x >= 0 && x < self.width && y >= 0 && y < self.height {
    let v = if val < 0 { 0 } else if val > 255 { 255 } else { val }
    self.data[y * self.width + x] = v
  }
}

///|
/// Fills the entire image with a constant pixel value.
pub fn Image::fill(self : Image, val : Int) -> Unit {
  let v = if val < 0 { 0 } else if val > 255 { 255 } else { val }
  let len = self.data.length()
  for i in 0.. Image {
  { width: self.width, height: self.height, data: self.data.copy() }
}

///|
/// Crops a subregion from the image.
pub fn Image::crop(
  self : Image,
  x : Int,
  y : Int,
  width : Int,
  height : Int,
) -> Image {
  let w = if width > 0 { width } else { 0 }
  let h = if height > 0 { height } else { 0 }
  let res = Image::new(w, h)
  for cy in 0.. Unit {
  for sy in 0..= 0 && dx < self.width && dy >= 0 && dy < self.height {
        self.set(dx, dy, src.get(sx, sy))
      }
    }
  }
}

///|
/// Resizes the image using bilinear interpolation.
pub fn Image::resize(self : Image, new_width : Int, new_height : Int) -> Image {
  if new_width <= 0 || new_height <= 0 {
    return Image::new(0, 0)
  }
  let res = Image::new(new_width, new_height)
  let sx = self.width.to_double() / new_width.to_double()
  let sy = self.height.to_double() / new_height.to_double()

  for y in 0.. Double {
  let x = pt.x
  let y = pt.y

  let x0 = x.floor().to_int()
  let y0 = y.floor().to_int()
  let x1 = x0 + 1
  let y1 = y0 + 1

  let dx = x - x0.to_double()
  let dy = y - y0.to_double()

  let clamp_x = fn(px : Int) -> Int {
    if px < 0 {
      0
    } else if px >= self.width {
      self.width - 1
    } else {
      px
    }
  }
  let clamp_y = fn(py : Int) -> Int {
    if py < 0 {
      0
    } else if py >= self.height {
      self.height - 1
    } else {
      py
    }
  }

  let p00 = self.get(clamp_x(x0), clamp_y(y0)).to_double()
  let p10 = self.get(clamp_x(x1), clamp_y(y0)).to_double()
  let p01 = self.get(clamp_x(x0), clamp_y(y1)).to_double()
  let p11 = self.get(clamp_x(x1), clamp_y(y1)).to_double()

  let top = p00 + dx * (p10 - p00)
  let bottom = p01 + dx * (p11 - p01)
  top + dy * (bottom - top)
}

///|
/// Bicubic interpolation sampling at sub-pixel coordinate pt with Catmull-Rom spline.
pub fn Image::sample_bicubic(self : Image, pt : Point2D) -> Double {
  let x = pt.x
  let y = pt.y
  let x0 = x.floor().to_int()
  let y0 = y.floor().to_int()
  let dx = x - x0.to_double()
  let dy = y - y0.to_double()

  let cubic_interp = fn(
    p0 : Double,
    p1 : Double,
    p2 : Double,
    p3 : Double,
    t : Double,
  ) -> Double {
    let a = -0.5 * p0 + 1.5 * p1 - 1.5 * p2 + 0.5 * p3
    let b = p0 - 2.5 * p1 + 2.0 * p2 - 0.5 * p3
    let c = -0.5 * p0 + 0.5 * p2
    let d = p1
    a * t * t * t + b * t * t + c * t + d
  }

  let clamp_x = fn(px : Int) -> Int {
    if px < 0 {
      0
    } else if px >= self.width {
      self.width - 1
    } else {
      px
    }
  }
  let clamp_y = fn(py : Int) -> Int {
    if py < 0 {
      0
    } else if py >= self.height {
      self.height - 1
    } else {
      py
    }
  }

  let row_vals = Array::make(4, 0.0)
  for m in 0..<4 {
    let cy = clamp_y(y0 - 1 + m)
    let p0 = self.get(clamp_x(x0 - 1), cy).to_double()
    let p1 = self.get(clamp_x(x0), cy).to_double()
    let p2 = self.get(clamp_x(x0 + 1), cy).to_double()
    let p3 = self.get(clamp_x(x0 + 2), cy).to_double()
    row_vals[m] = cubic_interp(p0, p1, p2, p3, dx)
  }

  let val = cubic_interp(row_vals[0], row_vals[1], row_vals[2], row_vals[3], dy)
  if val < 0.0 {
    0.0
  } else if val > 255.0 {
    255.0
  } else {
    val
  }
}

///|
/// Draws a line using Bresenham's algorithm.
pub fn Image::draw_line(
  self : Image,
  p0 : Point2D,
  p1 : Point2D,
  color : Int,
) -> Unit {
  let mut x0 = p0.x.round().to_int()
  let mut y0 = p0.y.round().to_int()
  let x1 = p1.x.round().to_int()
  let y1 = p1.y.round().to_int()

  let dx = (x1 - x0).abs()
  let dy = -(y1 - y0).abs()
  let sx = if x0 < x1 { 1 } else { -1 }
  let sy = if y0 < y1 { 1 } else { -1 }
  let mut err = dx + dy

  for ;; {
    self.set(x0, y0, color)
    if x0 == x1 && y0 == y1 {
      break
    }
    let e2 = 2 * err
    if e2 >= dy {
      err += dy
      x0 += sx
    }
    if e2 <= dx {
      err += dx
      y0 += sy
    }
  }
}

///|
/// Draws the 4 edges of a quadrilateral.
pub fn Image::draw_quad(self : Image, quad : Quad, color : Int) -> Unit {
  if quad.corners.length() != 4 {
    return
  }
  let p0 = quad.corners[0]
  let p1 = quad.corners[1]
  let p2 = quad.corners[2]
  let p3 = quad.corners[3]
  self.draw_line(p0, p1, color)
  self.draw_line(p1, p2, color)
  self.draw_line(p2, p3, color)
  self.draw_line(p3, p0, color)
}

///|
/// Draws a filled rectangle on the image.
pub fn Image::draw_rect_filled(
  self : Image,
  x : Int,
  y : Int,
  w : Int,
  h : Int,
  color : Int,
) -> Unit {
  for cy in 0.. Unit {
  let p0 = Point2D::new(x.to_double(), y.to_double())
  let p1 = Point2D::new((x + w).to_double(), y.to_double())
  let p2 = Point2D::new((x + w).to_double(), (y + h).to_double())
  let p3 = Point2D::new(x.to_double(), (y + h).to_double())
  self.draw_line(p0, p1, color)
  self.draw_line(p1, p2, color)
  self.draw_line(p2, p3, color)
  self.draw_line(p3, p0, color)
}

///|
/// Draws a circle using midpoint circle algorithm.
pub fn Image::draw_circle(
  self : Image,
  center : Point2D,
  radius : Int,
  color : Int,
) -> Unit {
  let cx = center.x.round().to_int()
  let cy = center.y.round().to_int()
  let mut x = radius
  let mut y = 0
  let mut err = 0

  while x >= y {
    self.set(cx + x, cy + y, color)
    self.set(cx + y, cy + x, color)
    self.set(cx - y, cy + x, color)
    self.set(cx - x, cy + y, color)
    self.set(cx - x, cy - y, color)
    self.set(cx - y, cy - x, color)
    self.set(cx + y, cy - x, color)
    self.set(cx + x, cy - y, color)

    if err <= 0 {
      y += 1
      err += 2 * y + 1
    }
    if err > 0 {
      x -= 1
      err -= 2 * x + 1
    }
  }
}

///|
/// Computes 2D Integral Image for fast O(1) box sum queries.
pub fn compute_integral_image(img : Image) -> Array[Int] {
  let w = img.width
  let h = img.height
  let size = (w + 1) * (h + 1)
  let integral = Array::make(size, 0)

  for y in 0.. Array[Double] {
  let w = img.width
  let h = img.height
  let size = (w + 1) * (h + 1)
  let integral_sq = Array::make(size, 0.0)

  for y in 0.. Int {
  let stride = w + 1
  let a = integral[y0 * stride + x0]
  let b = integral[y0 * stride + (x1 + 1)]
  let c = integral[(y1 + 1) * stride + x0]
  let d = integral[(y1 + 1) * stride + (x1 + 1)]
  d - b - c + a
}

///|
/// Exports image to NetPBM PGM (ASCII P2) format string.
pub fn Image::to_pgm(self : Image) -> String {
  let buf = StringBuilder()
  buf.write_string("P2\n")
  buf.write_string(self.width.to_string())
  buf.write_string(" ")
  buf.write_string(self.height.to_string())
  buf.write_string("\n255\n")
  for y in 0..