const clampPercent = (p) => {
  if (!Number.isFinite(p)) return 0
  if (p < 0) return 0
  if (p > 100) return 100
  return p
}

const buildHistogram = (gray) => {
  const hist = new Uint32Array(256)
  for (let i = 0; i < gray.length; i++) hist[gray[i]] += 1
  return hist
}

const otsuThreshold = (gray) => {
  const hist = buildHistogram(gray)
  const total = gray.length || 1
  let sum = 0
  for (let t = 0; t < 256; t++) sum += t * hist[t]

  let sumB = 0
  let wB = 0
  let varMax = -1
  let threshold = 127

  for (let t = 0; t < 256; t++) {
    wB += hist[t]
    if (wB === 0) continue
    const wF = total - wB
    if (wF === 0) break

    sumB += t * hist[t]
    const mB = sumB / wB
    const mF = (sum - sumB) / wF
    const between = wB * wF * (mB - mF) * (mB - mF)
    if (between > varMax) {
      varMax = between
      threshold = t
    }
  }

  return threshold
}

const percentile = (gray, p) => {
  const pct = clampPercent(p)
  const hist = buildHistogram(gray)
  const total = gray.length
  if (!total) return 0
  const target = Math.max(0, Math.min(total - 1, Math.round((pct / 100) * (total - 1))))
  let acc = 0
  for (let value = 0; value < hist.length; value++) {
    acc += hist[value]
    if (acc > target) return value
  }
  return 255
}

const clipRescale = (gray, low, high) => {
  const out = new Uint8Array(gray.length)
  if (!gray.length) return out
  if (high <= low) {
    out.fill(0)
    return out
  }
  const range = high - low
  for (let i = 0; i < gray.length; i++) {
    const clamped = Math.min(high, Math.max(low, gray[i]))
    out[i] = Math.max(0, Math.min(255, Math.round(((clamped - low) * 255) / range)))
  }
  return out
}

const binarize = (gray, threshold, invert = false) => {
  const out = new Uint8Array(gray.length)
  if (!invert) {
    for (let i = 0; i < gray.length; i++) out[i] = gray[i] >= threshold ? 255 : 0
  } else {
    for (let i = 0; i < gray.length; i++) out[i] = gray[i] < threshold ? 255 : 0
  }
  return out
}

const morphDilate = (bin, width, height, radius = 1) => {
  const out = new Uint8Array(bin.length)
  const R = Math.max(0, Math.round(radius))
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      let any = 0
      for (let dy = -R; dy <= R && !any; dy++) {
        const ny = y + dy
        if (ny < 0 || ny >= height) continue
        for (let dx = -R; dx <= R; dx++) {
          const nx = x + dx
          if (nx < 0 || nx >= width) continue
          if (bin[ny * width + nx]) {
            any = 1
            break
          }
        }
      }
      out[y * width + x] = any ? 255 : 0
    }
  }
  return out
}

const morphErode = (bin, width, height, radius = 1) => {
  const out = new Uint8Array(bin.length)
  const R = Math.max(0, Math.round(radius))
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      let all = 1
      for (let dy = -R; dy <= R && all; dy++) {
        const ny = y + dy
        if (ny < 0 || ny >= height) {
          all = 0
          break
        }
        for (let dx = -R; dx <= R; dx++) {
          const nx = x + dx
          if (nx < 0 || nx >= width) {
            all = 0
            break
          }
          if (!bin[ny * width + nx]) {
            all = 0
            break
          }
        }
      }
      out[y * width + x] = all ? 255 : 0
    }
  }
  return out
}

const morphOpen = (bin, width, height, radius = 1) =>
  morphDilate(morphErode(bin, width, height, radius), width, height, radius)

const morphClose = (bin, width, height, radius = 1) =>
  morphErode(morphDilate(bin, width, height, radius), width, height, radius)

const removeSmallComponents = (bin, width, height, minArea = 0) => {
  const total = width * height
  if (!total || bin.length < total || minArea <= 0) return Uint8Array.from(bin)
  const out = Uint8Array.from(bin)
  const visited = new Uint8Array(total)
  const queue = new Int32Array(total)
  const component = []

  for (let i = 0; i < total; i++) {
    if (!out[i] || visited[i]) continue
    component.length = 0
    let head = 0
    let tail = 0
    queue[tail++] = i
    visited[i] = 1

    while (head < tail) {
      const idx = queue[head++]
      component.push(idx)
      const x = idx % width
      const y = (idx / width) | 0
      const tryPush = (nx, ny) => {
        if (nx < 0 || ny < 0 || nx >= width || ny >= height) return
        const nidx = ny * width + nx
        if (visited[nidx] || !out[nidx]) return
        visited[nidx] = 1
        queue[tail++] = nidx
      }
      tryPush(x - 1, y)
      tryPush(x + 1, y)
      tryPush(x, y - 1)
      tryPush(x, y + 1)
    }

    if (component.length < minArea) {
      for (const idx of component) out[idx] = 0
    }
  }
  return out
}

const fillSmallHoles = (bin, width, height, maxArea = 0) => {
  const total = width * height
  if (!total || bin.length < total || maxArea <= 0) return Uint8Array.from(bin)
  const out = Uint8Array.from(bin)
  const visited = new Uint8Array(total)
  const queue = new Int32Array(total)
  const component = []

  for (let i = 0; i < total; i++) {
    if (out[i] || visited[i]) continue
    component.length = 0
    let head = 0
    let tail = 0
    let touchesBorder = false
    queue[tail++] = i
    visited[i] = 1

    while (head < tail) {
      const idx = queue[head++]
      component.push(idx)
      const x = idx % width
      const y = (idx / width) | 0
      if (x === 0 || y === 0 || x === width - 1 || y === height - 1) touchesBorder = true
      const tryPush = (nx, ny) => {
        if (nx < 0 || ny < 0 || nx >= width || ny >= height) return
        const nidx = ny * width + nx
        if (visited[nidx] || out[nidx]) return
        visited[nidx] = 1
        queue[tail++] = nidx
      }
      tryPush(x - 1, y)
      tryPush(x + 1, y)
      tryPush(x, y - 1)
      tryPush(x, y + 1)
    }

    if (!touchesBorder && component.length <= maxArea) {
      for (const idx of component) out[idx] = 255
    }
  }
  return out
}

const clampByte = (v) => {
  if (!Number.isFinite(v)) return 0
  if (v < 0) return 0
  if (v > 255) return 255
  return v & 0xff
}

const boxBlur = (gray, width, height, radius = 1) => {
  const r = Math.max(1, Math.round(radius))
  const stride = width + 1
  const integral = new Float64Array((height + 1) * stride)

  const idx = (x, y) => y * stride + x

  for (let y = 1; y <= height; y++) {
    let rowSum = 0
    for (let x = 1; x <= width; x++) {
      const value = gray[(y - 1) * width + (x - 1)]
      rowSum += value
      integral[idx(x, y)] = integral[idx(x, y - 1)] + rowSum
    }
  }

  const out = new Float32Array(width * height)
  for (let y = 0; y < height; y++) {
    const y0 = Math.max(0, y - r)
    const y1 = Math.min(height - 1, y + r)
    for (let x = 0; x < width; x++) {
      const x0 = Math.max(0, x - r)
      const x1 = Math.min(width - 1, x + r)
      const A = idx(x0, y0)
      const B = idx(x1 + 1, y0)
      const C = idx(x0, y1 + 1)
      const D = idx(x1 + 1, y1 + 1)
      const area = (x1 - x0 + 1) * (y1 - y0 + 1)
      const sum = integral[D] - integral[B] - integral[C] + integral[A]
      out[y * width + x] = area > 0 ? sum / area : gray[y * width + x]
    }
  }
  return out
}

const illuminationNormalize = (gray, width, height, radius = 0, amount = 1) => {
  const r = Math.max(0, Math.round(radius))
  const amt = Math.max(0, Math.min(1, Number.isFinite(amount) ? amount : 1))
  if (!r || !amt) return Uint8Array.from(gray)
  const blurred = boxBlur(gray, width, height, r)
  const out = new Uint8Array(gray.length)
  for (let i = 0; i < gray.length; i++) {
    const base = gray[i]
    const refl = clampByte(base - blurred[i] + 128)
    const mixed = clampByte(base + amt * (refl - base))
    out[i] = mixed
  }
  return out
}

const sobel = (gray, width, height) => {
  const out = new Uint8Array(gray.length)
  const gxK = [-1, 0, 1, -2, 0, 2, -1, 0, 1]
  const gyK = [-1, -2, -1, 0, 0, 0, 1, 2, 1]
  for (let y = 1; y < height - 1; y++) {
    for (let x = 1; x < width - 1; x++) {
      let gx = 0
      let gy = 0
      let k = 0
      for (let dy = -1; dy <= 1; dy++) {
        const ny = y + dy
        for (let dx = -1; dx <= 1; dx++, k++) {
          const nx = x + dx
          const value = gray[ny * width + nx]
          gx += value * gxK[k]
          gy += value * gyK[k]
        }
      }
      const magnitude = Math.min(255, Math.round(Math.hypot(gx, gy)))
      out[y * width + x] = magnitude
    }
  }
  return out
}

const removeBorderTouching = (bin, width, height, minArea = 0) => {
  const total = bin.length
  const labels = new Int32Array(total)
  const areas = [0]
  const stack = []
  let id = 0

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const index = y * width + x
      if (!bin[index] || labels[index]) continue
      id += 1
      let area = 0
      let touchesBorder = false
      stack.push(index)
      labels[index] = id

      while (stack.length) {
        const current = stack.pop()
        area += 1
        const cx = current % width
        const cy = Math.floor(current / width)
        if (cx === 0 || cy === 0 || cx === width - 1 || cy === height - 1) touchesBorder = true
        const neighbours = [current - 1, current + 1, current - width, current + width]
        for (const nb of neighbours) {
          if (nb < 0 || nb >= total) continue
          if (!bin[nb] || labels[nb]) continue
          labels[nb] = id
          stack.push(nb)
        }
      }

      areas[id] = touchesBorder ? 0 : area
    }
  }

  const out = new Uint8Array(total)
  for (let i = 0; i < total; i++) {
    const area = areas[labels[i]] || 0
    if (area >= minArea) out[i] = 255
  }
  return out
}

const keepLargest = (bin, width, height, count = 1) => {
  const total = bin.length
  const labels = new Int32Array(total)
  const areas = [0]
  const stack = []
  let id = 0

  for (let index = 0; index < total; index++) {
    if (!bin[index] || labels[index]) continue
    id += 1
    let area = 0
    stack.push(index)
    labels[index] = id

    while (stack.length) {
      const current = stack.pop()
      area += 1
      const neighbours = [current - 1, current + 1, current - width, current + width]
      for (const nb of neighbours) {
        if (nb < 0 || nb >= total) continue
        if (!bin[nb] || labels[nb]) continue
        labels[nb] = id
        stack.push(nb)
      }
    }

    areas[id] = area
  }

  const ids = Array.from({ length: id }, (_, i) => i + 1)
    .sort((a, b) => (areas[b] || 0) - (areas[a] || 0))
    .slice(0, Math.max(1, Math.round(count)))
  const keep = new Uint8Array(id + 1)
  for (const entry of ids) keep[entry] = 1

  const out = new Uint8Array(total)
  for (let i = 0; i < total; i++) {
    if (keep[labels[i]]) out[i] = 255
  }
  return out
}

const invert255 = (mask) => {
  const out = new Uint8Array(mask.length)
  for (let i = 0; i < mask.length; i++) out[i] = mask[i] ? 0 : 255
  return out
}

const largestBorderConnected = (bin, width, height) => {
  const total = width * height
  if (total === 0) return new Uint8Array(0)
  const visited = new Uint8Array(total)
  const labels = new Int32Array(total)
  let label = 0
  let bestLabel = 0
  let bestArea = 0
  const stack = []

  const push = (index) => {
    if (index < 0 || index >= total) return
    if (visited[index] || !bin[index]) return
    visited[index] = 1
    stack.push(index)
  }

  const starts = []
  for (let x = 0; x < width; x++) {
    const top = x
    const bottom = (height - 1) * width + x
    if (bin[top]) starts.push(top)
    if (bin[bottom]) starts.push(bottom)
  }
  for (let y = 1; y < height - 1; y++) {
    const left = y * width
    const right = y * width + (width - 1)
    if (bin[left]) starts.push(left)
    if (bin[right]) starts.push(right)
  }

  for (const start of starts) {
    if (visited[start]) continue
    label += 1
    let area = 0
    push(start)
    while (stack.length) {
      const current = stack.pop()
      labels[current] = label
      area += 1
      const x = current % width
      const y = (current / width) | 0
      if (x > 0) push(current - 1)
      if (x + 1 < width) push(current + 1)
      if (y > 0) push(current - width)
      if (y + 1 < height) push(current + width)
    }
    if (area > bestArea) {
      bestArea = area
      bestLabel = label
    }
  }

  const out = new Uint8Array(total)
  if (!bestArea) return out
  for (let i = 0; i < total; i++) if (labels[i] === bestLabel) out[i] = 255
  return out
}

const dropSmall = (bin, width, height, minArea = 1) => {
  const total = bin.length
  if (total === 0) return new Uint8Array(0)
  const labels = new Int32Array(total)
  const areas = [0]
  const stack = []
  let id = 0

  for (let index = 0; index < total; index++) {
    if (!bin[index] || labels[index]) continue
    id += 1
    let area = 0
    labels[index] = id
    stack.push(index)
    while (stack.length) {
      const current = stack.pop()
      area += 1
      const x = current % width
      const y = (current / width) | 0
      if (x > 0) {
        const nb = current - 1
        if (!labels[nb] && bin[nb]) {
          labels[nb] = id
          stack.push(nb)
        }
      }
      if (x + 1 < width) {
        const nb = current + 1
        if (!labels[nb] && bin[nb]) {
          labels[nb] = id
          stack.push(nb)
        }
      }
      if (y > 0) {
        const nb = current - width
        if (!labels[nb] && bin[nb]) {
          labels[nb] = id
          stack.push(nb)
        }
      }
      if (y + 1 < height) {
        const nb = current + width
        if (!labels[nb] && bin[nb]) {
          labels[nb] = id
          stack.push(nb)
        }
      }
    }
    areas[id] = area
  }

  const out = new Uint8Array(total)
  for (let i = 0; i < total; i++) {
    const area = areas[labels[i]] || 0
    if (area >= minArea) out[i] = 255
  }
  return out
}

const subtractEdges = (bin, gray, width, height, opts = {}) => {
  const thr = opts.thr ?? 48
  const dilate = opts.dilate ?? 1
  let edges = sobel(gray, width, height)
  edges = binarize(edges, thr, false)
  if (dilate > 0) edges = morphDilate(edges, width, height, dilate)
  const out = new Uint8Array(bin.length)
  for (let i = 0; i < bin.length; i++) out[i] = bin[i] && !edges[i] ? 255 : 0
  return out
}

const maskOr = (a, b) => {
  const len = Math.min(a.length, b.length)
  const out = new Uint8Array(len)
  for (let i = 0; i < len; i++) out[i] = (a[i] || b[i]) ? 255 : 0
  return out
}

const maskAnd = (a, b) => {
  const len = Math.min(a.length, b.length)
  const out = new Uint8Array(len)
  for (let i = 0; i < len; i++) out[i] = (a[i] && b[i]) ? 255 : 0
  return out
}

const maskMinus = (a, b) => {
  const len = Math.min(a.length, b.length)
  const out = new Uint8Array(len)
  for (let i = 0; i < len; i++) out[i] = a[i] && !b[i] ? 255 : 0
  return out
}

const stripBorder = (mask, width, height, margin = 0) => {
  const m = Math.max(0, Math.round(margin))
  const out = Uint8Array.from(mask)
  if (!m) return out
  const w = Math.max(1, width)
  const h = Math.max(1, height)
  for (let y = 0; y < h; y++) {
    const row = y * w
    if (y < m || y >= h - m) {
      for (let x = 0; x < w; x++) out[row + x] = 0
      continue
    }
    for (let x = 0; x < m; x++) out[row + x] = 0
    for (let x = Math.max(m, w - m); x < w; x++) out[row + x] = 0
  }
  return out
}

const distanceTransformL1 = (mask, width, height) => {
  const total = width * height
  const dist = new Int32Array(total)
  const INF = 1e9
  for (let i = 0; i < total; i++) dist[i] = mask[i] ? 0 : INF

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const idx = y * width + x
      if (x > 0) dist[idx] = Math.min(dist[idx], dist[idx - 1] + 1)
      if (y > 0) dist[idx] = Math.min(dist[idx], dist[idx - width] + 1)
    }
  }

  for (let y = height - 1; y >= 0; y--) {
    for (let x = width - 1; x >= 0; x--) {
      const idx = y * width + x
      if (x + 1 < width) dist[idx] = Math.min(dist[idx], dist[idx + 1] + 1)
      if (y + 1 < height) dist[idx] = Math.min(dist[idx], dist[idx + width] + 1)
    }
  }

  return dist
}

const gateByDistance = (edges, seed, width, height, maxDist) => {
  const dist = distanceTransformL1(seed, width, height)
  const limit = Math.max(0, Math.round(maxDist))
  const out = new Uint8Array(edges.length)
  for (let i = 0; i < edges.length; i++) out[i] = edges[i] && dist[i] <= limit ? 255 : 0
  return out
}

const geodesicGrow = (seed, mask, width, height, iterations = 6) => {
  const iters = Math.max(0, Math.round(iterations))
  let current = new Uint8Array(seed)
  if (!iters) {
    for (let i = 0; i < current.length; i++) current[i] = current[i] && mask[i] ? 255 : 0
    return current
  }
  for (let k = 0; k < iters; k++) {
    current = morphDilate(current, width, height, 1)
    for (let i = 0; i < current.length; i++) current[i] = current[i] && mask[i] ? 255 : 0
  }
  return current
}

const suppressGridEdges = (edges, width, height, vFrac = 0.55, hFrac = 0.55, pad = 1) => {
  if (!edges.length) return new Uint8Array(0)
  const out = new Uint8Array(edges.length)
  const columnCounts = new Int32Array(width)
  const rowCounts = new Int32Array(height)
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      if (edges[y * width + x]) {
        columnCounts[x] += 1
        rowCounts[y] += 1
      }
    }
  }

  const vMask = new Uint8Array(width)
  const hMask = new Uint8Array(height)
  const colThreshold = Math.max(0, Math.min(height, Math.round(vFrac * height)))
  const rowThreshold = Math.max(0, Math.min(width, Math.round(hFrac * width)))
  const radius = Math.max(0, Math.round(pad))

  for (let x = 0; x < width; x++) {
    if (columnCounts[x] >= colThreshold) {
      for (let dx = -radius; dx <= radius; dx++) {
        const nx = x + dx
        if (nx >= 0 && nx < width) vMask[nx] = 1
      }
    }
  }

  for (let y = 0; y < height; y++) {
    if (rowCounts[y] >= rowThreshold) {
      for (let dy = -radius; dy <= radius; dy++) {
        const ny = y + dy
        if (ny >= 0 && ny < height) hMask[ny] = 1
      }
    }
  }

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const idx = y * width + x
      if (edges[idx] && !(vMask[x] || hMask[y])) out[idx] = 255
    }
  }

  return out
}

const suppressGridEdgesAdaptive = (edges, width, height, sigma = 2.5, pad = 1) => {
  if (!edges.length) return new Uint8Array(0)
  const out = new Uint8Array(edges.length)
  const columnCounts = new Float64Array(width)
  const rowCounts = new Float64Array(height)

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      if (edges[y * width + x]) {
        columnCounts[x] += 1
        rowCounts[y] += 1
      }
    }
  }

  const meanAndStd = (arr) => {
    let sum = 0
    for (let i = 0; i < arr.length; i++) sum += arr[i]
    const mean = arr.length ? sum / arr.length : 0
    let variance = 0
    for (let i = 0; i < arr.length; i++) {
      const diff = arr[i] - mean
      variance += diff * diff
    }
    const std = arr.length ? Math.sqrt(variance / arr.length) : 0
    return { mean, std }
  }

  const { mean: meanCol, std: stdCol } = meanAndStd(columnCounts)
  const { mean: meanRow, std: stdRow } = meanAndStd(rowCounts)
  const colThreshold = meanCol + sigma * stdCol
  const rowThreshold = meanRow + sigma * stdRow
  const radius = Math.max(0, Math.round(pad))
  const killCols = new Uint8Array(width)
  const killRows = new Uint8Array(height)

  for (let x = 0; x < width; x++) {
    if (columnCounts[x] >= colThreshold) {
      for (let dx = -radius; dx <= radius; dx++) {
        const nx = x + dx
        if (nx >= 0 && nx < width) killCols[nx] = 1
      }
    }
  }

  for (let y = 0; y < height; y++) {
    if (rowCounts[y] >= rowThreshold) {
      for (let dy = -radius; dy <= radius; dy++) {
        const ny = y + dy
        if (ny >= 0 && ny < height) killRows[ny] = 1
      }
    }
  }

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const idx = y * width + x
      if (edges[idx] && !(killCols[x] || killRows[y])) out[idx] = 255
    }
  }

  return out
}

const suppressThinLines = (mask, width, height, maxWidth = 3, minLengthFrac = 0.25) => {
  if (!mask.length) return new Uint8Array(0)
  const total = width * height
  if (total !== mask.length || width <= 0 || height <= 0) return Uint8Array.from(mask)

  const working = Uint8Array.from(mask)
  const keep = new Uint8Array(total)
  const radiusMin = Math.max(1, Math.ceil(Math.max(1, Math.round(maxWidth)) / 2))
  const dist = distanceTransformL1(working, width, height)
  for (let i = 0; i < total; i++) {
    if (working[i] && dist[i] >= radiusMin) keep[i] = 255
  }

  const minFrac = Math.max(0, Math.min(1, Number(minLengthFrac)))
  const minRun = Math.max(1, Math.round(minFrac * Math.max(width, height)))

  for (let y = 0; y < height; y++) {
    let run = 0
    let runStart = 0
    for (let x = 0; x <= width; x++) {
      const idx = y * width + x
      const slender = x < width && working[idx] && !keep[idx]
      if (slender) {
        if (run === 0) runStart = x
        run += 1
      } else if (run) {
        if (run >= minRun) {
          for (let xx = runStart; xx < runStart + run; xx++) working[y * width + xx] = 0
        }
        run = 0
      }
    }
  }

  for (let x = 0; x < width; x++) {
    let run = 0
    let runStart = 0
    for (let y = 0; y <= height; y++) {
      const idx = y * width + x
      const slender = y < height && working[idx] && !keep[idx]
      if (slender) {
        if (run === 0) runStart = y
        run += 1
      } else if (run) {
        if (run >= minRun) {
          for (let yy = runStart; yy < runStart + run; yy++) working[yy * width + x] = 0
        }
        run = 0
      }
    }
  }

  const out = new Uint8Array(total)
  for (let i = 0; i < total; i++) out[i] = keep[i] || working[i] ? 255 : 0
  return out
}

const localStd = (gray, width, height, window = 7) => {
  const win = Math.max(1, Math.round(window))
  const stride = width + 1
  const integral = new Float64Array(stride * (height + 1))
  const integralSq = new Float64Array(stride * (height + 1))

  const idx = (x, y) => y * stride + x

  for (let y = 1; y <= height; y++) {
    let rowSum = 0
    let rowSumSq = 0
    for (let x = 1; x <= width; x++) {
      const value = gray[(y - 1) * width + (x - 1)]
      rowSum += value
      rowSumSq += value * value
      integral[idx(x, y)] = integral[idx(x, y - 1)] + rowSum
      integralSq[idx(x, y)] = integralSq[idx(x, y - 1)] + rowSumSq
    }
  }

  const out = new Float32Array(width * height)
  for (let y = 0; y < height; y++) {
    const y0 = Math.max(0, y - win)
    const y1 = Math.min(height - 1, y + win)
    for (let x = 0; x < width; x++) {
      const x0 = Math.max(0, x - win)
      const x1 = Math.min(width - 1, x + win)
      const A = idx(x0, y0)
      const B = idx(x1 + 1, y0)
      const C = idx(x0, y1 + 1)
      const D = idx(x1 + 1, y1 + 1)
      const area = (x1 - x0 + 1) * (y1 - y0 + 1)
      if (area <= 0) {
        out[y * width + x] = 0
        continue
      }
      const sum = integral[D] - integral[B] - integral[C] + integral[A]
      const sumSq = integralSq[D] - integralSq[B] - integralSq[C] + integralSq[A]
      const mean = sum / area
      const variance = Math.max(0, sumSq / area - mean * mean)
      out[y * width + x] = Math.sqrt(variance)
    }
  }
  return out
}

const hysteresisEdges = (gradient, width, height, highThreshold, lowThreshold) => {
  const total = width * height
  const strong = 255
  const weak = 128
  const high = Math.max(0, Math.round(highThreshold))
  const low = Math.max(0, Math.min(high, Math.round(lowThreshold)))
  const labels = new Uint8Array(total)
  const output = new Uint8Array(total)

  for (let i = 0; i < total; i++) {
    const value = gradient[i]
    labels[i] = value >= high ? strong : value >= low ? weak : 0
  }

  const stack = []
  const offsets = [-width - 1, -width, -width + 1, -1, 1, width - 1, width, width + 1]

  for (let i = 0; i < total; i++) {
    if (labels[i] !== strong || output[i]) continue
    output[i] = strong
    stack.push(i)
    while (stack.length) {
      const current = stack.pop()
      const cy = Math.floor(current / width)
      const cx = current - cy * width
      for (const offset of offsets) {
        const neighbour = current + offset
        if (neighbour < 0 || neighbour >= total) continue
        const ny = Math.floor(neighbour / width)
        const nx = neighbour - ny * width
        if (Math.abs(nx - cx) > 1 || Math.abs(ny - cy) > 1) continue
        if (labels[neighbour] === weak && !output[neighbour]) {
          output[neighbour] = strong
          stack.push(neighbour)
        }
      }
    }
  }

  return output
}

const meanStdMasked = (gray, mask) => {
  let sum = 0
  let sumSq = 0
  let count = 0
  const len = Math.min(gray.length, mask.length)
  for (let i = 0; i < len; i++) {
    if (!mask[i]) continue
    const v = gray[i]
    sum += v
    sumSq += v * v
    count += 1
  }
  if (!count) return { mean: 128, std: 16 }
  const mean = sum / count
  const variance = Math.max(0, sumSq / count - mean * mean)
  const std = Math.sqrt(variance)
  return { mean, std }
}

module.exports = {
  otsuThreshold,
  percentile,
  clipRescale,
  binarize,
  morphDilate,
  morphErode,
  morphOpen,
  morphClose,
  removeSmallComponents,
  fillSmallHoles,
  boxBlur,
  illuminationNormalize,
  sobel,
  removeBorderTouching,
  keepLargest,
  invert255,
  largestBorderConnected,
  dropSmall,
  subtractEdges,
  maskOr,
  maskAnd,
  maskMinus,
  stripBorder,
  distanceTransformL1,
  gateByDistance,
  geodesicGrow,
  suppressGridEdges,
  suppressGridEdgesAdaptive,
  suppressThinLines,
  localStd,
  hysteresisEdges,
  meanStdMasked,
}