A set of helper functions for manipulating 3D bitmask in sparse (run-based) format

bitmaskHelpers3d.ts

import SparseBitmask from "./SparseBitmask";
import { bitmaskFromWeightedRuns } from "./bitmaskHelpers";
import { OrderedMap } from "./OrderedCollections";
import { range } from "./base";

/**
 * [Claude generated]
 *
 * `cutAndPasteBitmask3d` copies a 3-D sub-volume (the "cut" box) out of `from`
 * and XOR-pastes it into `to` at an offset position. Both the cut and paste
 * boxes are clipped to their respective canvas bounds before use, and their
 * intersection in source-space is the only region actually processed.
 *
 * It works entirely in run-length space: `getWindowedRunLengths3d` streams the
 * source region as a sequence of signed run lengths, and the consumer loop
 * replays those runs into the destination by calling `flip` on the two boundary
 * indices of each 1-run — never touching individual bits. The three traversal
 * modes (row / yMode / zMode) mirror the sentinel protocol of
 * `getWindowedRunLengths3d`: a sentinel value of `windowSizeX` or `windowSizeXY`
 * signals that the next value is a full-row or full-slice count rather than a
 * single-row run, letting the loop skip entire rows or slices in one step.
 *
 * Convenience wrappers:
 *   `cutBitmask3d`    — paste at origin (crop / resize the canvas)
 *   `pasteBitmask3d`  — cut from origin (stamp the full source onto `to`)
 *   `offsetBitmask3d` — same canvas size, shift content by (dx, dy, dz)
 */
export function cutAndPasteBitmask3d(
  from: SparseBitmask,
  fromWidth: number,
  fromHeight: number,
  fromDepth: number,
  cutX: number,
  cutY: number,
  cutZ: number,
  cutWidth: number,
  cutHeight: number,
  cutDepth: number,
  pasteX: number,
  pasteY: number,
  pasteZ: number,
  toWidth: number,
  toHeight: number,
  toDepth: number,
  to = new SparseBitmask(toWidth * toHeight * toDepth, from.isRunEnds)
): SparseBitmask {
  if (from.getIndicesCount() === 0) return to;

  const appliedTo = to.isRunEnds ? to : new SparseBitmask(to.length, true);

  // compute the part of the input mask that is within the cut window
  const cutMinX = Math.max(cutX, 0);
  const cutMinY = Math.max(cutY, 0);
  const cutMinZ = Math.max(cutZ, 0);
  const cutMaxX = Math.min(cutX + cutWidth, fromWidth);
  const cutMaxY = Math.min(cutY + cutHeight, fromHeight);
  const cutMaxZ = Math.min(cutZ + cutDepth, fromDepth);

  // compute the part of cut window that is within the output canvas translated to input frame
  const pasteMinX = Math.max(-pasteX, 0) + cutX;
  const pasteMinY = Math.max(-pasteY, 0) + cutY;
  const pasteMinZ = Math.max(-pasteZ, 0) + cutZ;
  const pasteMaxX = Math.min(-pasteX + toWidth, cutWidth) + cutX;
  const pasteMaxY = Math.min(-pasteY + toHeight, cutHeight) + cutY;
  const pasteMaxZ = Math.min(-pasteZ + toDepth, cutDepth) + cutZ;

  // compute the intersection
  const minX = Math.max(cutMinX, pasteMinX);
  const minY = Math.max(cutMinY, pasteMinY);
  const minZ = Math.max(cutMinZ, pasteMinZ);
  const maxX = Math.min(cutMaxX, pasteMaxX);
  const maxY = Math.min(cutMaxY, pasteMaxY);
  const maxZ = Math.min(cutMaxZ, pasteMaxZ);

  if (minX >= maxX || minY >= maxY || minZ >= maxZ) return to;

  const translateX = pasteX - cutX;
  const translateY = pasteY - cutY;
  const translateZ = pasteZ - cutZ;
  // convert to run lengths windowed by row for subsequent interlacing
  const windowedRunLengths = getWindowedRunLengths3d(
    from,
    fromWidth,
    fromHeight,
    fromDepth,
    minX,
    minY,
    minZ,
    maxX,
    maxY,
    maxZ
  );

  const toWH = toWidth * toHeight;
  const windowSizeX = maxX - minX;
  const windowSizeY = maxY - minY;
  const windowSizeXY = windowSizeX * windowSizeY;
  const gapY = toWidth - windowSizeX;
  const gapZ = (toHeight - windowSizeY) * toWidth;

  let offset =
    (minZ + translateZ) * toWH +
    (minY + translateY) * toWidth +
    minX +
    translateX;
  let rowLength = 0;
  let sliceLength = 0;

  let yMode = false;
  let zMode = false;
  for (const runLength of windowedRunLengths) {
    if (rowLength >= windowSizeX) {
      offset += gapY;
      rowLength = 0;
    }
    if (sliceLength >= windowSizeXY) {
      offset += gapZ;
      sliceLength = 0;
    }

    if (zMode) {
      if (runLength < 0) {
        if (gapY > 0) {
          for (const _ of range(runLength, 0)) {
            for (const __ of range(windowSizeY)) {
              appliedTo.flip(offset);
              appliedTo.flip(offset + windowSizeX);
              offset += toWidth;
            }
            offset += gapZ;
          }
        } else if (gapZ > 0) {
          for (const _ of range(runLength, 0)) {
            appliedTo.flip(offset);
            appliedTo.flip(offset + windowSizeXY);
            offset += toWH;
          }
        } else {
          appliedTo.flip(offset);
          appliedTo.flip(offset - runLength * toWH);
          offset -= runLength * toWH;
        }
      } else {
        offset += runLength * toWH;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        if (gapY > 0) {
          for (const _ of range(runLength, 0)) {
            appliedTo.flip(offset);
            appliedTo.flip(offset + windowSizeX);
            offset += toWidth;
          }
        } else {
          appliedTo.flip(offset);
          appliedTo.flip(offset - runLength * toWidth);
          offset -= runLength * toWidth;
        }
        sliceLength -= runLength * windowSizeX;
      } else {
        offset += runLength * toWidth;
        sliceLength += runLength * windowSizeX;
      }
      yMode = false;
    } else if (runLength === windowSizeXY) {
      zMode = true;
    } else if (runLength === windowSizeX) {
      yMode = true;
    } else if (runLength < 0) {
      appliedTo.flip(offset);
      appliedTo.flip(offset - runLength);
      offset -= runLength;
      rowLength -= runLength;
      sliceLength -= runLength;
    } else {
      offset += runLength;
      rowLength += runLength;
      sliceLength += runLength;
    }
  }

  if (!to.isRunEnds) {
    for (const index of appliedTo.getIndices(false)) {
      to.set(index, 1);
    }
  }

  return to;
}

export function cutBitmask3d(
  from: SparseBitmask,
  fromWidth: number,
  fromHeight: number,
  fromDepth: number,
  x: number,
  y: number,
  z: number,
  toWidth: number,
  toHeight: number,
  toDepth: number
): SparseBitmask {
  if (
    x === 0 &&
    y === 0 &&
    z === 0 &&
    fromWidth === toWidth &&
    fromHeight === toHeight &&
    fromDepth === toDepth
  ) {
    return from;
  }
  return cutAndPasteBitmask3d(
    from,
    fromWidth,
    fromHeight,
    fromDepth,
    x,
    y,
    z,
    toWidth,
    toHeight,
    toDepth,
    0,
    0,
    0,
    toWidth,
    toHeight,
    toDepth
  );
}

export function pasteBitmask3d(
  from: SparseBitmask,
  fromWidth: number,
  fromHeight: number,
  fromDepth: number,
  x: number,
  y: number,
  z: number,
  toWidth: number,
  toHeight: number,
  toDepth: number,
  to = new SparseBitmask(toWidth * toHeight * toDepth, from.isRunEnds)
): SparseBitmask {
  return cutAndPasteBitmask3d(
    from,
    fromWidth,
    fromHeight,
    fromDepth,
    0,
    0,
    0,
    fromWidth,
    fromHeight,
    fromDepth,
    x,
    y,
    z,
    toWidth,
    toHeight,
    toDepth,
    to
  );
}

export function offsetBitmask3d(
  mask: SparseBitmask,
  offsetX: number,
  offsetY: number,
  offsetZ: number,
  width: number,
  height: number,
  depth: number
): SparseBitmask {
  return cutAndPasteBitmask3d(
    mask,
    width,
    height,
    depth,
    0,
    0,
    0,
    width,
    height,
    depth,
    offsetX,
    offsetY,
    offsetZ,
    width,
    height,
    depth
  );
}

/**
 * [Claude generated]
 *
 * The core problem: run-end indices are flat 1-D positions, so a single run
 * can straddle a row or slice boundary inside a windowed sub-volume. That makes
 * them awkward for 2-D/3-D operations that need to reason row-by-row or
 * slice-by-slice.
 *
 * Run lengths solve this cleanly — a run that crosses a boundary is simply split
 * into two shorter runs at that boundary. Because run lengths don't have to
 * alternate sign the way run-end indices do, two adjacent runs of the same sign
 * are perfectly valid and arise naturally at every split point.
 *
 * `sliceRunEnds` clips the raw run-end stream to [startIndex, endIndex] first,
 * injecting synthetic boundary markers so the main loop only ever sees indices
 * within the window.
 */
/**
 * Decompose run-end indices into 1 or 0 run lengths that are windowed by each row
 * Stream a sequence of run lengths
 * Negative value means run of 1s, positive value means run of 0s
 * Runs that crosses window boundaries will be decomposed into shorter runs so every run length stays within a single row
 *
 * The exception is full row runs or full slice runs (y traversal or z traversal mode):
 * - First a "special" value that equates to window size is yielded
 * - The user should then take the next yield value to represent the number of full rows or full slices this run spans
 * - Sign logic applies as above
 * - Sequence will resume in row traversal mode unless another special value is yielded
 *
 * Unlike run end indices that always indicate a change of sign, adjacent run lengths can be of the same sign
 * This is the base function for many 2D and 3D bitmask operations
 */
export function* getWindowedRunLengths3d(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth = 1,
  minX = 0,
  minY = 0,
  minZ = 0,
  maxX = width,
  maxY = height,
  maxZ = depth
): Iterable<number> {
  if (minX >= maxX || minY >= maxY || minZ >= maxZ) return;
  const wh = width * height;
  const windowSizeX = maxX - minX;
  const windowSizeY = maxY - minY;
  const windowSizeXY = windowSizeX * windowSizeY;
  let sign = 1;
  const startIndex = minZ * wh + minY * width + minX;
  const endIndex = (maxZ - 1) * wh + (maxY - 1) * width + maxX;
  let currRowStart = startIndex;
  let currSliceStart = currRowStart;
  let currRowEnd = currRowStart + windowSizeX;
  let currSliceEnd = currRowStart + (windowSizeY - 1) * width + windowSizeX;
  let lastIndex = currRowStart;

  for (const index of sliceRunEnds(
    mask.getIndices(true),
    startIndex,
    endIndex
  )) {
    while (index >= currRowEnd) {
      if (lastIndex > currRowStart) {
        // not full row, yield run until the end of the row
        if (currRowEnd > lastIndex) {
          yield sign * (currRowEnd - lastIndex);
        }
        if (currRowEnd >= currSliceEnd) {
          currSliceStart += wh;
          currSliceEnd += wh;
          currRowStart = currSliceStart;
          currRowEnd = currRowStart + windowSizeX;
        } else {
          currRowStart += width;
          currRowEnd += width;
        }
        lastIndex = currRowStart;
      } else if (index >= currSliceEnd) {
        // handle slice boundary
        while (index >= currSliceEnd) {
          if (lastIndex > currSliceStart) {
            // not full slice, yield run until the end of the slice
            if (currSliceEnd > lastIndex) {
              const k = Math.floor((currSliceEnd - currRowEnd) / width) + 1;
              yield windowSizeX;
              yield sign * k;
            }
            currSliceStart += wh;
            currSliceEnd += wh;
          } else {
            // full slice run
            const k = Math.floor((index - currSliceEnd) / wh) + 1;
            yield windowSizeXY;
            yield sign * k;
            currSliceStart += k * wh;
            currSliceEnd += k * wh;
          }
          currRowStart = currSliceStart;
          currRowEnd = currRowStart + windowSizeX;
          lastIndex = currRowStart;
        }
      } else {
        // full row run
        const k = Math.floor((index - currRowEnd) / width) + 1;
        yield windowSizeX;
        yield sign * k;
        currRowStart += k * width;
        currRowEnd += k * width;
        lastIndex = currRowStart;
      }
    }

    if (index > currRowStart) {
      yield sign * (index - lastIndex);
      lastIndex = index;
    }

    sign = -sign;
  }
}

function* sliceRunEnds(
  runEnds: Iterable<number>,
  startIndex: number,
  endIndex: number
): Iterable<number> {
  let curr = 0;
  let started = false;
  let ended = false;
  for (const index of runEnds) {
    if (!started && index > startIndex) {
      if (curr) {
        yield startIndex;
      }
      started = true;
    }
    if (index >= endIndex) {
      yield endIndex;
      ended = true;
      break;
    }
    if (started) yield index;
    curr = 1 - curr;
  }
  if (!started && curr) {
    yield startIndex;
  }
  if (!ended) {
    yield endIndex;
  }
}

/**
 * [Claude generated]
 *
 * Transposes the X and Y axes (x↔y, z unchanged).
 *
 * Rather than iterating every set bit, it exploits run structure: subtracting
 * the mask shifted by ±1 in Y isolates the first and last rows of each Y-run.
 * Only those boundary rows need to be mapped into the transposed coordinate
 * space (y + x*height + z*wh), keeping the operation proportional to the
 * number of run boundaries rather than the number of set bits.
 */
export function switchPrimaryAxis(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): SparseBitmask {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  const flipped = new SparseBitmask(mask.length, true);
  const wh = width * height;
  const yMinusOne = offsetBitmask3d(mask, 0, -1, 0, width, height, depth);
  const yPlusOne = offsetBitmask3d(mask, 0, 1, 0, width, height, depth);
  for (const index of mask.subtraction(yPlusOne).getIndices(false)) {
    const x = index % width;
    const y = Math.trunc((index % wh) / width);
    const z = Math.trunc(index / wh);
    flipped.flip(z * wh + x * height + y);
  }
  for (const index of mask.subtraction(yMinusOne).getIndices(false)) {
    const x = index % width;
    const y = Math.trunc((index % wh) / width);
    const z = Math.trunc(index / wh);
    flipped.flip(z * wh + x * height + y + 1);
  }
  return mask.isRunEnds ? flipped : flipped.asOneBitIndices();
}

/**
 * [Claude generated]
 *
 * Transposes the Y and Z axes (y↔z, x unchanged).
 *
 * Uses getWindowedRunLengths3d to stream run boundaries and maps each directly
 * to the transposed coordinate space (x + z*width + y*wd). The three traversal
 * modes handle row, y, and z runs separately since each maps to a different
 * target stride in the output layout.
 */
export function switchSecondaryAxis(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): SparseBitmask {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  const flipped = new SparseBitmask(mask.length, true);
  const wh = width * height;
  const wd = width * depth;
  let xStart = 0;
  let yStart = 0;
  let zStart = 0;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, width, height, depth)) {
    if (xStart >= width) {
      xStart = 0;
      yStart += 1;
    }
    if (yStart >= height) {
      yStart = 0;
      zStart += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const zEnd = zStart - runLength;
        for (const y of range(height)) {
          flipped.flip(y * wd + zStart * width);
          flipped.flip(y * wd + zEnd * width);
        }
        zStart = zEnd;
      } else {
        zStart += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        const yEnd = yStart - runLength;
        for (const y of range(yStart, yEnd)) {
          flipped.flip(y * wd + zStart * width);
          flipped.flip(y * wd + (zStart + 1) * width);
        }
        yStart = yEnd;
      } else {
        yStart += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      const xEnd = xStart - runLength;
      flipped.flip(yStart * wd + zStart * width + xStart);
      flipped.flip(yStart * wd + zStart * width + xEnd);
      xStart = xEnd;
    } else {
      xStart += runLength;
    }
  }
  return mask.isRunEnds ? flipped : flipped.asOneBitIndices();
}

/**
 * [Claude generated]
 *
 * Axis order convention:
 *   0, 1, 2   → x, y, z (identity, no flip)
 *  -1, -2, -3 → x, y, z flipped
 * The position in the array is the output axis; the value is which input axis
 * (and whether it is flipped) maps there. A trailing [0,1,2] is filled in for
 * any omitted axes.
 *
 * Returns the inverse permutation: given axisOrder mapping input→output,
 * produces the axisOrder that maps output→input, preserving flip signs.
 */
export function getInverseAxisOrder(axisOrder: number[]): number[] {
  const normAxisOrder = [...axisOrder, ...[0, 1, 2].slice(axisOrder.length)];
  const inverseAxisOrder: number[] = [];
  normAxisOrder.forEach((axis, index) => {
    inverseAxisOrder[axis < 0 ? -axis - 1 : axis] =
      axis < 0 ? -index - 1 : index;
  });
  return inverseAxisOrder;
}

/**
 * [Claude generated] Composes two axis orders: equivalent to applying `to` first, then `apply`.
 * Flip signs are combined — flipping a flipped axis yields an unflipped one.
 */
export function chainAxisOrder(apply: number[], to: number[]): number[] {
  const normTo = [...to, ...[0, 1, 2].slice(to.length)];
  const normApply = [...apply, ...[0, 1, 2].slice(apply.length)];
  return normApply.map(i => {
    const v = normTo[i < 0 ? -i - 1 : i];
    return i < 0 ? -v - 1 : v;
  });
}

/**
 * [Claude generated] Returns the [width, height, depth] of the output after applying axisOrder.
 * Flip signs are ignored — only the axis mapping affects dimension sizes.
 */
export function getRotatedDimensions(
  axisOrder: number[],
  width: number,
  height: number,
  depth = 1
): [number, number, number] {
  const normAxisOrder = [...axisOrder, ...[0, 1, 2].slice(axisOrder.length)];
  const dims = [width, height, depth];
  return [
    dims[normAxisOrder[0] < 0 ? -normAxisOrder[0] - 1 : normAxisOrder[0]],
    dims[normAxisOrder[1] < 0 ? -normAxisOrder[1] - 1 : normAxisOrder[1]],
    dims[normAxisOrder[2] < 0 ? -normAxisOrder[2] - 1 : normAxisOrder[2]],
  ];
}

/**
 * [Claude generated]
 *
 * Applies an axis permutation (ignoring flip signs) by decomposing it into at
 * most two switchPrimaryAxis / switchSecondaryAxis calls. All 6 permutations of
 * xyz are covered; flip signs in the axis order are handled separately by
 * flipAxes / permuteAndFlipAxes.
 */
export function permuteAxes(
  mask: SparseBitmask,
  axisOrder: number[],
  width: number,
  height: number,
  depth = 1
): SparseBitmask {
  const normAxisOrder = [...axisOrder, ...[0, 1, 2].slice(axisOrder.length)];
  const permutation = normAxisOrder.map(i => (i < 0 ? -i - 1 : i)).join("");

  if (permutation === "012") {
    return mask;
  }
  if (permutation === "120") {
    const yxz = switchPrimaryAxis(mask, width, height, depth);
    return switchSecondaryAxis(yxz, height, width, depth);
  }
  if (permutation === "201") {
    const xzy = switchSecondaryAxis(mask, width, height, depth);
    return switchPrimaryAxis(xzy, width, depth, height);
  }
  if (permutation === "210") {
    const xzy = switchSecondaryAxis(mask, width, height, depth);
    const zxy = switchPrimaryAxis(xzy, width, depth, height);
    return switchSecondaryAxis(zxy, depth, width, height);
  }
  if (permutation === "102") {
    return switchPrimaryAxis(mask, width, height, depth);
  }
  if (permutation === "021") {
    return switchSecondaryAxis(mask, width, height, depth);
  }

  throw new Error(`unsupported axis order: ${permutation}`);
}

/**
 * [Claude generated]
 *
 * Mirrors the mask along any combination of axes. Uses getWindowedRunLengths3d
 * to stream runs and reflects each boundary index across the relevant axis
 * (e.g. flipped x: xStart → width - xEnd). Permutation is handled separately
 * by permuteAxes; this only mirrors, never reorders axes.
 */
export function flipAxes(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number,
  flipX = false,
  flipY = false,
  flipZ = false
): SparseBitmask {
  if (!flipX && !flipY && !flipZ) return mask;
  const flipped = new SparseBitmask(mask.length, true);

  const wh = width * height;
  let x = 0;
  let y = 0;
  let z = 0;
  let zMode = false;
  let yMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, width, height, depth)) {
    if (x >= width) {
      x = 0;
      y += 1;
    }
    if (y >= height) {
      y = 0;
      z += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const zStart = flipZ ? depth - (z - runLength) : z;
        const zEnd = flipZ ? depth - z : z - runLength;
        flipped.flip(zStart * wh);
        flipped.flip(zEnd * wh);
        z -= runLength;
      } else {
        z += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        const yStart = flipY ? height - (y - runLength) : y;
        const yEnd = flipY ? height - y : y - runLength;
        const zStart = flipZ ? depth - (z + 1) : z;
        flipped.flip(zStart * wh + yStart * width);
        flipped.flip(zStart * wh + yEnd * width);
        y -= runLength;
      } else {
        y += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      const xStart = flipX ? width - (x - runLength) : x;
      const xEnd = flipX ? width - x : x - runLength;
      const yStart = flipY ? height - (y + 1) : y;
      const zStart = flipZ ? depth - (z + 1) : z;
      flipped.flip(zStart * wh + yStart * width + xStart);
      flipped.flip(zStart * wh + yStart * width + xEnd);
      x -= runLength;
    } else {
      x += runLength;
    }
  }

  return mask.isRunEnds ? flipped : flipped.asOneBitIndices();
}

/**
 * [Claude generated] Applies a full axis order (permutation + flips) in two steps: permuteAxes
 * first (ignoring signs), then flipAxes on the rotated dimensions using the
 * sign bits to determine which axes to mirror.
 */
export function permuteAndFlipAxes(
  mask: SparseBitmask,
  axisOrder: number[],
  width: number,
  height: number,
  depth: number
): SparseBitmask {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  const normAxisOrder = [...axisOrder, ...[0, 1, 2].slice(axisOrder.length)];
  const permuted = permuteAxes(mask, normAxisOrder, width, height, depth);
  const [rotatedWidth, rotatedHeight, rotatedDepth] = getRotatedDimensions(
    axisOrder,
    width,
    height,
    depth
  );
  const flipped = flipAxes(
    permuted,
    rotatedWidth,
    rotatedHeight,
    rotatedDepth,
    normAxisOrder[0] < 0,
    normAxisOrder[1] < 0,
    normAxisOrder[2] < 0
  );
  return mask.isRunEnds ? flipped : flipped.asOneBitIndices();
}

/**
 * [Claude generated]
 *
 * MPR (multi-planar reconstruction) support.
 *
 * `slice3d` decomposes a 3D mask into its full set of 2D planar views — one
 * XY mask per z-slice, one XZ mask per y-column, one YZ mask per x-row — in
 * two passes: the original layout covers XY and XZ in one sweep; a second
 * sweep over the y↔x transposed mask covers YZ.
 *
 * The `propagate*` functions are the inverse: given a 2D annotation on one
 * plane, they project it back into all three planar views plus a thin 3D mask
 * (xyz) representing the annotated slice in 3D space. This lets callers
 * display the effect of a 2D edit across all MPR views at once.
 */

export interface Sliced3dBitmask {
  xyz: SparseBitmask[];
  xy: SparseBitmask[];
  xz: SparseBitmask[];
  yz: SparseBitmask[];
}

// slice a 3D bitmask into three sets of 2D bitmasks on the XY, XZ, and YZ planes
export function slice3d(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): Sliced3dBitmask {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }
  const xy = Array.from(
    { length: depth },
    () => new SparseBitmask(width * height, true)
  );
  const xz = Array.from(
    { length: height },
    () => new SparseBitmask(width * depth, true)
  );
  const yz = Array.from(
    { length: width },
    () => new SparseBitmask(height * depth, true)
  );

  const wh = width * height;
  let xStart = 0;
  let yStart = 0;
  let zStart = 0;
  let xMode = false;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, width, height, depth)) {
    if (xStart >= width) {
      xStart = 0;
      yStart += 1;
    }
    if (yStart >= height) {
      yStart = 0;
      zStart += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const zEnd = zStart - runLength;
        for (const z of range(zStart, zEnd)) {
          xy[z].flip(0);
        }
        for (const y of range(height)) {
          xz[y].flip(zStart * width);
          xz[y].flip(zEnd * width);
        }
        zStart = zEnd;
      } else {
        zStart += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        const yEnd = yStart - runLength;
        xy[zStart].flip(yStart * width);
        xy[zStart].flip(yEnd * width);
        for (const y of range(yStart, yEnd)) {
          xz[y].flip(zStart * width);
          xz[y].flip((zStart + 1) * width);
        }
        yStart = yEnd;
      } else {
        yStart += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      const xEnd = xStart - runLength;
      xy[zStart].flip(yStart * width + xStart);
      xy[zStart].flip(yStart * width + xEnd);
      xz[yStart].flip(zStart * width + xStart);
      xz[yStart].flip(zStart * width + xEnd);
      xStart = xEnd;
    } else {
      xStart += runLength;
    }
  }

  const yxz = switchPrimaryAxis(mask, width, height, depth);
  // reset
  xStart = 0;
  yStart = 0;
  zStart = 0;
  for (const runLength of getWindowedRunLengths3d(yxz, height, width, depth)) {
    if (yStart >= height) {
      yStart = 0;
      xStart += 1;
    }
    if (xStart >= width) {
      xStart = 0;
      zStart += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const zEnd = zStart - runLength;
        for (const x of range(width)) {
          yz[x].flip(zStart * height);
          yz[x].flip(zEnd * height);
        }
        zStart = zEnd;
      } else {
        zStart += runLength;
      }
      zMode = false;
    } else if (xMode) {
      if (runLength < 0) {
        const xEnd = xStart - runLength;
        for (const x of range(xStart, xEnd)) {
          yz[x].flip(zStart * height);
          yz[x].flip((zStart + 1) * height);
        }
        xStart = xEnd;
      } else {
        xStart += runLength;
      }
      xMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === height) {
      xMode = true;
    } else if (runLength < 0) {
      const yEnd = yStart - runLength;
      yz[xStart].flip(zStart * height + yStart);
      yz[xStart].flip(zStart * height + yEnd);
      yStart = yEnd;
    } else {
      yStart += runLength;
    }
  }

  return { xyz: [mask], xy, xz, yz };
}

/** [Claude generated] Projects a 2D XY mask at z=zIndex into all three MPR views and a 3D mask. */
export function propagateXY(
  mask: SparseBitmask,
  zIndex: number,
  width: number,
  height: number,
  depth: number
): Sliced3dBitmask {
  if (mask.length !== width * height) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }
  const xy = Array.from(
    { length: width },
    () => new SparseBitmask(width * height, true)
  );
  const xz = Array.from(
    { length: height },
    () => new SparseBitmask(width * depth, true)
  );
  const yz = Array.from(
    { length: width },
    () => new SparseBitmask(height * depth, true)
  );
  const xyz = new SparseBitmask(width * height * depth, true);

  xy[zIndex] = mask;

  const wh = width * height;
  let xStart = 0;
  let yStart = 0;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, width, height)) {
    if (xStart >= width) {
      xStart = 0;
      yStart += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        xyz.flip(zIndex * wh);
        xyz.flip((zIndex + 1) * wh);
        for (const y of range(height)) {
          xz[y].flip(zIndex * width);
          xz[y].flip((zIndex + 1) * width);
        }
      }
      break;
    } else if (yMode) {
      if (runLength < 0) {
        const yEnd = yStart - runLength;
        xyz.flip(zIndex * wh + yStart * width);
        xyz.flip(zIndex * wh + yEnd * width);
        for (const y of range(yStart, yEnd)) {
          xz[y].flip(zIndex * width);
          xz[y].flip((zIndex + 1) * width);
        }
        yStart = yEnd;
      } else {
        yStart += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      const xEnd = xStart - runLength;
      xyz.flip(zIndex * wh + yStart * width + xStart);
      xyz.flip(zIndex * wh + yStart * width + xEnd);
      xz[yStart].flip(zIndex * width + xStart);
      xz[yStart].flip(zIndex * width + xEnd);
      xStart = xEnd;
    } else {
      xStart += runLength;
    }
  }

  const yx = switchPrimaryAxis(mask, width, height, 1);

  // reset
  xStart = 0;
  yStart = 0;
  let xMode = false;
  zMode = false;
  for (const runLength of getWindowedRunLengths3d(yx, height, width)) {
    if (yStart >= height) {
      yStart = 0;
      xStart += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        for (const x of range(width)) {
          yz[x].flip(zIndex * height);
          yz[x].flip((zIndex + 1) * height);
        }
      }
      break;
    } else if (xMode) {
      if (runLength < 0) {
        const xEnd = xStart - runLength;
        for (const x of range(xStart, xEnd)) {
          yz[x].flip(zIndex * height);
          yz[x].flip((zIndex + 1) * height);
        }
        xStart = xEnd;
      } else {
        xStart += runLength;
      }
      xMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === height) {
      xMode = true;
    } else if (runLength < 0) {
      const yEnd = yStart - runLength;
      yz[xStart].flip(zIndex * height + yStart);
      yz[xStart].flip(zIndex * height + yEnd);
      yStart = yEnd;
    } else {
      yStart += runLength;
    }
  }

  return { xyz: [xyz], xy, xz, yz };
}

/** [Claude generated] Projects a 2D XZ mask at y=yIndex into all three MPR views and a 3D mask. */
export function propagateXZ(
  mask: SparseBitmask,
  yIndex: number,
  width: number,
  height: number,
  depth: number
): Sliced3dBitmask {
  if (mask.length !== width * depth) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }
  const xy = Array.from(
    { length: depth },
    () => new SparseBitmask(width * height, true)
  );
  const xz = Array.from(
    { length: height },
    () => new SparseBitmask(width * depth, true)
  );
  const yz = Array.from(
    { length: width },
    () => new SparseBitmask(height * depth, true)
  );
  const xyz = new SparseBitmask(width * height * depth, true);

  xz[yIndex] = mask;

  const wh = width * height;
  const wd = width * depth;
  let xStart = 0;
  let zStart = 0;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, width, depth)) {
    if (xStart >= width) {
      xStart = 0;
      zStart += 1;
    }

    if (yMode) {
      if (runLength < 0) {
        for (const z of range(depth)) {
          xyz.flip(z * wh + yIndex * width);
          xyz.flip(z * wh + (yIndex + 1) * width);
          xy[z].flip(yIndex * width);
          xy[z].flip((yIndex + 1) * width);
          for (const x of range(width)) {
            yz[x].flip(z * height + yIndex);
            yz[x].flip(z * height + yIndex + 1);
          }
        }
      }
      break;
    } else if (zMode) {
      if (runLength < 0) {
        const zEnd = zStart - runLength;
        for (const z of range(zStart, zEnd)) {
          xyz.flip(z * wh + yIndex * width);
          xyz.flip(z * wh + (yIndex + 1) * width);
          xy[z].flip(yIndex * width);
          xy[z].flip((yIndex + 1) * width);
          for (const x of range(width)) {
            yz[x].flip(z * height + yIndex);
            yz[x].flip(z * height + yIndex + 1);
          }
        }
        zStart = zEnd;
      } else {
        zStart += runLength;
      }
      zMode = false;
    } else if (runLength === wd) {
      yMode = true;
    } else if (runLength === width) {
      zMode = true;
    } else if (runLength < 0) {
      const xEnd = xStart - runLength;
      xyz.flip(zStart * wh + yIndex * width + xStart);
      xyz.flip(zStart * wh + yIndex * width + xEnd);
      xy[zStart].flip(yIndex * width + xStart);
      xy[zStart].flip(yIndex * width + xEnd);
      for (const x of range(xStart, xEnd)) {
        yz[x].flip(zStart * height + yIndex);
        yz[x].flip(zStart * height + yIndex + 1);
      }
      xStart = xEnd;
    } else {
      xStart += runLength;
    }
  }

  return { xyz: [xyz], xy, xz, yz };
}

/** [Claude generated] Projects a 2D YZ mask at x=xIndex into all three MPR views and a 3D mask. */
export function propagateYZ(
  mask: SparseBitmask,
  xIndex: number,
  width: number,
  height: number,
  depth: number
): Sliced3dBitmask {
  if (mask.length !== height * depth) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }
  const xy = Array.from(
    { length: depth },
    () => new SparseBitmask(width * height, true)
  );
  const xz = Array.from(
    { length: height },
    () => new SparseBitmask(width * depth, true)
  );
  const yz = Array.from(
    { length: width },
    () => new SparseBitmask(height * depth, true)
  );
  const xyz = new SparseBitmask(width * height * depth, true);

  yz[xIndex] = mask;

  const wh = width * height;
  const hd = height * depth;
  let yStart = 0;
  let zStart = 0;
  let xMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(mask, height, depth)) {
    if (yStart >= height) {
      yStart = 0;
      zStart += 1;
    }

    if (xMode) {
      if (runLength < 0) {
        for (const z of range(depth)) {
          for (const y of range(height)) {
            xyz.flip(z * wh + y * width + xIndex);
            xyz.flip(z * wh + y * width + xIndex + 1);
            xy[z].flip(y * width + xIndex);
            xy[z].flip(y * width + xIndex + 1);
            xz[y].flip(z * width + xIndex);
            xz[y].flip(z * width + xIndex + 1);
          }
        }
      }
      break;
    } else if (zMode) {
      if (runLength < 0) {
        const zEnd = zStart - runLength;
        for (const z of range(zStart, zEnd)) {
          for (const y of range(height)) {
            xyz.flip(z * wh + y * width + xIndex);
            xyz.flip(z * wh + y * width + xIndex + 1);
            xy[z].flip(y * width + xIndex);
            xy[z].flip(y * width + xIndex + 1);
            xz[y].flip(z * width + xIndex);
            xz[y].flip(z * width + xIndex + 1);
          }
        }
        zStart = zEnd;
      } else {
        zStart += runLength;
      }
      zMode = false;
    } else if (runLength === hd) {
      xMode = true;
    } else if (runLength === height) {
      zMode = true;
    } else if (runLength < 0) {
      const yEnd = yStart - runLength;
      for (const y of range(yStart, yEnd)) {
        xyz.flip(zStart * wh + y * width + xIndex);
        xyz.flip(zStart * wh + y * width + xIndex + 1);
        xy[zStart].flip(y * width + xIndex);
        xy[zStart].flip(y * width + xIndex + 1);
        xz[y].flip(zStart * width + xIndex);
        xz[y].flip(zStart * width + xIndex + 1);
      }
      yStart = yEnd;
    } else {
      yStart += runLength;
    }
  }

  return { xyz: [xyz], xy, xz, yz };
}

/** [Claude generated] Like slice3d but returns only the XY slices (one per z-layer). */
export function slice3dXY(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): SparseBitmask[] {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }
  const sliced = Array.from(
    { length: depth },
    () => new SparseBitmask(width * height, true)
  );
  const wh = width * height;
  const whd = wh * depth;
  let xy = 0;
  let z = 0;
  let zMode = false;
  let xyzMode = false;
  for (const runLength of getWindowedRunLengths3d(
    mask,
    width * height,
    depth
  )) {
    if (xy >= wh) {
      xy = 0;
      z += 1;
    }

    if (xyzMode) {
      if (runLength < 0) {
        for (const _ of range(depth)) {
          sliced[z].flip(0);
          z += 1;
        }
        break;
      }
    } else if (zMode) {
      if (runLength < 0) {
        for (const _ of range(runLength, 0)) {
          sliced[z].flip(0);
          z += 1;
        }
      } else {
        z += runLength;
      }
      zMode = false;
    } else if (runLength === whd) {
      xyzMode = true;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength < 0) {
      sliced[z].flip(xy);
      sliced[z].flip(xy - runLength);
      xy -= runLength;
    } else {
      xy += runLength;
    }
  }
  return sliced;
}

/**
 * [Claude generated]
 *
 * Resamples a 3-D bitmask to new dimensions using area-weighted interpolation.
 *
 * Each source 1-run is mapped to fractional destination coordinates via
 * (scaleX, scaleY, scaleZ). The overlap of that mapped run with each destination
 * voxel is accumulated as a weight in an OrderedMap keyed by run-end index.
 * Partial voxels at the start/end of a run contribute their fractional overlap
 * (wIStart, wIEnd, etc.); fully covered voxels contribute 1. The weights are
 * stored as deltas at transition points — the same run-end principle as the
 * bitmask itself — so the final pass in `bitmaskFromWeightedRuns` just sweeps
 * through them to recover a running sum, thresholding at 0.5 to produce the
 * output mask.
 */
export function resizeBitmask3d(
  mask: SparseBitmask,
  fromWidth: number,
  fromHeight: number,
  fromDepth: number,
  toWidth: number,
  toHeight: number,
  toDepth: number
): SparseBitmask {
  if (fromWidth * fromHeight * fromDepth !== mask.length) {
    throw new Error("mismatch parameters");
  }
  if (!mask.isRunEnds) {
    throw new Error("mask must be in run end format");
  }

  const runWeights = new OrderedMap<number, number>((_, k) => k);
  const scaleX = toWidth / fromWidth;
  const scaleY = toHeight / fromHeight;
  const scaleZ = toDepth / fromDepth;

  const fromWH = fromWidth * fromHeight;
  const toWH = toWidth * toHeight;
  let x = 0;
  let y = 0;
  let z = 0;
  let yMode = false;
  let zMode = false;

  for (const runLength of getWindowedRunLengths3d(
    mask,
    fromWidth,
    fromHeight,
    fromDepth
  )) {
    if (x >= fromWidth) {
      x = 0;
      y += 1;
    }
    if (y >= fromHeight) {
      y = 0;
      z += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const kStart = z * scaleZ;
        const kEnd = (z - runLength) * scaleZ;
        const kStartInt = Math.trunc(kStart);
        const kEndInt = Math.trunc(kEnd);
        const wKStart = kStart - kStartInt;
        const wKEnd = kEnd - kEndInt;

        if (wKEnd > 0) {
          runWeights.update(kEndInt * toWH, (weight = 0) => weight + wKEnd);
          runWeights.update(
            (kEndInt + 1) * toWH,
            (weight = 0) => weight - wKEnd
          );
        }
        if (kEndInt > kStartInt) {
          runWeights.update(kStartInt * toWH, (weight = 0) => weight + 1);
          runWeights.update(kEndInt * toWH, (weight = 0) => weight - 1);
        }
        if (wKStart > 0) {
          runWeights.update(kStartInt * toWH, (weight = 0) => weight - wKStart);
          runWeights.update(
            (kStartInt + 1) * toWH,
            (weight = 0) => weight + wKStart
          );
        }

        z -= runLength;
      } else {
        z += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        const kStart = z * scaleZ;
        const kEnd = (z + 1) * scaleZ;
        const kStartInt = Math.trunc(kStart);
        const kEndInt = Math.trunc(kEnd);
        const wKStart = kStart - kStartInt;
        const wKEnd = kEnd - kEndInt;
        const jStart = y * scaleY;
        const jEnd = (y - runLength) * scaleY;
        const jStartInt = Math.trunc(jStart);
        const jEndInt = Math.trunc(jEnd);
        const wJStart = jStart - jStartInt;
        const wJEnd = jEnd - jEndInt;

        const processSlice = (k: number, w: number) => {
          if (wJEnd > 0) {
            runWeights.update(
              k * toWH + jEndInt * toWidth,
              (weight = 0) => weight + wJEnd * w
            );
            runWeights.update(
              k * toWH + (jEndInt + 1) * toWidth,
              (weight = 0) => weight - wJEnd * w
            );
          }
          if (jEndInt > jStartInt) {
            runWeights.update(
              k * toWH + jStartInt * toWidth,
              (weight = 0) => weight + w
            );
            runWeights.update(
              k * toWH + jEndInt * toWidth,
              (weight = 0) => weight - w
            );
          }
          if (wJStart > 0) {
            runWeights.update(
              k * toWH + jStartInt * toWidth,
              (weight = 0) => weight - wJStart * w
            );
            runWeights.update(
              k * toWH + (jStartInt + 1) * toWidth,
              (weight = 0) => weight + wJStart * w
            );
          }
        };

        if (wKEnd > 0) processSlice(kEndInt, wKEnd);
        for (const k of range(kStartInt, kEndInt)) {
          processSlice(k, 1);
        }
        if (wKStart > 0) processSlice(kStartInt, -wKStart);

        y -= runLength;
      } else {
        y += runLength;
      }
      yMode = false;
    } else if (runLength === fromWH) {
      zMode = true;
    } else if (runLength === fromWidth) {
      yMode = true;
    } else if (runLength < 0) {
      const kStart = z * scaleZ;
      const kEnd = (z + 1) * scaleZ;
      const kStartInt = Math.trunc(kStart);
      const kEndInt = Math.trunc(kEnd);
      const wKStart = kStart - kStartInt;
      const wKEnd = kEnd - kEndInt;
      const jStart = y * scaleY;
      const jEnd = (y + 1) * scaleY;
      const jStartInt = Math.trunc(jStart);
      const jEndInt = Math.trunc(jEnd);
      const wJStart = jStart - jStartInt;
      const wJEnd = jEnd - jEndInt;
      const iStart = x * scaleX;
      const iEnd = (x - runLength) * scaleX;
      const iStartInt = Math.trunc(iStart);
      const iEndInt = Math.trunc(iEnd);
      const wIStart = iStart - iStartInt;
      const wIEnd = iEnd - iEndInt;

      const processRow = (j: number, k: number, w: number) => {
        if (wIEnd > 0) {
          runWeights.update(
            k * toWH + j * toWidth + iEndInt,
            (weight = 0) => weight + wIEnd * w
          );
          runWeights.update(
            k * toWH + j * toWidth + iEndInt + 1,
            (weight = 0) => weight - wIEnd * w
          );
        }
        if (iEndInt > iStartInt) {
          runWeights.update(
            k * toWH + j * toWidth + iStartInt,
            (weight = 0) => weight + w
          );
          runWeights.update(
            k * toWH + j * toWidth + iEndInt,
            (weight = 0) => weight - w
          );
        }
        if (wIStart > 0) {
          runWeights.update(
            k * toWH + j * toWidth + iStartInt,
            (weight = 0) => weight - wIStart * w
          );
          runWeights.update(
            k * toWH + j * toWidth + iStartInt + 1,
            (weight = 0) => weight + wIStart * w
          );
        }
      };

      const processSlice = (k: number, w: number) => {
        if (wJEnd > 0) processRow(jEndInt, k, wJEnd * w);
        for (const j of range(jStartInt, jEndInt)) {
          processRow(j, k, w);
        }
        if (wJStart > 0) processRow(jStartInt, k, -wJStart * w);
      };

      if (wKEnd > 0) processSlice(kEndInt, wKEnd);
      for (const k of range(kStartInt, kEndInt)) {
        processSlice(k, 1);
      }
      if (wKStart > 0) processSlice(kStartInt, -wKStart);

      x -= runLength;
    } else {
      x += runLength;
    }
  }
  return bitmaskFromWeightedRuns(runWeights, toWidth * toHeight * toDepth, 0.5);
}

/**
 * [Claude generated] Returns the centroid (centre of mass) of all set bits as [x, y, z].
 * Accumulates weighted sums over runs via getWindowedRunLengths3d so it never
 * expands individual bits. Falls back to geometric centre if the mask is empty.
 */
export function getCentroid3d(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): [number, number, number] {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }

  if (mask.getIndicesCount() === 0) return [width / 2, height / 2, depth / 2];

  let sumX = 0;
  let sumY = 0;
  let sumZ = 0;
  let weight = 0;

  const wh = width * height;
  let x = 0;
  let y = 0;
  let z = 0;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(
    mask.asRunEndIndices(),
    width,
    height,
    depth
  )) {
    if (x >= width) {
      x = 0;
      y += 1;
    }
    if (y >= height) {
      y = 0;
      z += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        const dz = -runLength;
        const dWeight = wh * dz;
        sumX += 0.5 * width * dWeight;
        sumY += 0.5 * height * dWeight;
        sumZ = (z + dz * 0.5) * dWeight;
        weight += dWeight;
        z -= runLength;
      } else {
        z += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        const dy = -runLength;
        const dWeight = width * dy;
        sumX += 0.5 * width * dWeight;
        sumY += (y + dy * 0.5) * dWeight;
        sumZ += (z + 0.5) * dWeight;
        weight += dWeight;
        y -= runLength;
      } else {
        y += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      const dx = -runLength;
      const dWeight = dx;
      sumX += (x + dx * 0.5) * dWeight;
      sumY += (y + 0.5) * dWeight;
      sumZ += (z + 0.5) * dWeight;
      weight += dWeight;
      x -= runLength;
    } else {
      x += runLength;
    }
  }
  return [sumX / weight, sumY / weight, sumZ / weight];
}

/**
 * [Claude generated] Returns the axis-aligned bounding box of all set bits as [minX, minY, minZ, maxX, maxY, maxZ].
 * If the mask is empty, returns an inverted box (min > max) as a sentinel.
 */
export function getBoundingBox3d(
  mask: SparseBitmask,
  width: number,
  height: number,
  depth: number
): [number, number, number, number, number, number] {
  if (mask.length !== width * height * depth) {
    throw new Error("mismatch parameters");
  }
  let minX = width;
  let minY = height;
  let minZ = depth;
  let maxX = 0;
  let maxY = 0;
  let maxZ = 0;

  if (mask.getIndicesCount() === 0) return [minX, minY, minZ, maxX, maxY, maxZ];

  const wh = width * height;
  let x = 0;
  let y = 0;
  let z = 0;
  let yMode = false;
  let zMode = false;
  for (const runLength of getWindowedRunLengths3d(
    mask.asRunEndIndices(),
    width,
    height,
    depth
  )) {
    if (x >= width) {
      x = 0;
      y += 1;
    }
    if (y >= height) {
      y = 0;
      z += 1;
    }

    if (zMode) {
      if (runLength < 0) {
        minX = Math.min(minX, 0);
        minY = Math.min(minY, 0);
        minZ = Math.min(minZ, z);
        maxX = Math.max(maxX, width);
        maxY = Math.max(maxY, height);
        maxZ = Math.max(maxZ, z - runLength);
        z -= runLength;
        zMode = false;
      } else {
        z += runLength;
      }
      zMode = false;
    } else if (yMode) {
      if (runLength < 0) {
        minX = Math.min(minX, 0);
        minY = Math.min(minY, y);
        minZ = Math.min(minZ, z);
        maxX = Math.max(maxX, width);
        maxY = Math.max(maxY, y + 1);
        maxZ = Math.max(maxZ, z + 1);
        y -= runLength;
      } else {
        y += runLength;
      }
      yMode = false;
    } else if (runLength === wh) {
      zMode = true;
    } else if (runLength === width) {
      yMode = true;
    } else if (runLength < 0) {
      minX = Math.min(minX, x);
      minY = Math.min(minY, y);
      minZ = Math.min(minZ, z);
      maxX = Math.max(maxX, x - runLength);
      maxY = Math.max(maxY, y + 1);
      maxZ = Math.max(maxZ, z + 1);
      x -= runLength;
    } else {
      x += runLength;
    }
  }
  return [minX, minY, minZ, maxX, maxY, maxZ];
}

/**
 * Using a digital differential analyzer to compute the sphere edge efficiently
 * without trigonometry, or Pythagoras.
 *
 * Returns run-end indices (sorted) for a hollow sphere of diameter n centred
 * at the origin. Exploits 8-fold octant symmetry — each DDA step emits up to
 * 8 reflected runs at once. Coincident boundary indices (where two runs meet
 * exactly) cancel out via XOR-merge: if the same index appears twice it is
 * popped rather than duplicated, which correctly collapses adjacent runs of the
 * same value in the output.
 */
export function drawSphere(n: number): number[] {
  const n2 = n * n;
  const runs = new Map<number, number>();
  const offset = (n - 1) / 2;
  const indexOffset = offset * (n2 + n + 1);
  let x0 = -offset;
  const y0 = n & 1 ? 0 : -0.5;
  let z = n & 1 ? 0 : -0.5;

  let d0 = x0 * x0 + y0 * y0 + z * z - (n * n) / 4;
  let deltaOut = 1 - 2 * z;
  const deltaUp0 = 1 - 2 * y0;
  let deltaRight0 = 2 * x0 + 1;

  const innerLoop = (
    x: number,
    y: number,
    d: number,
    deltaUp: number,
    deltaRight: number
  ) => {
    while (y > x) {
      runs.set(z * n2 + y * n + x, z * n2 + y * n - x + 1);
      runs.set(z * n2 + -y * n + x, z * n2 + -y * n - x + 1);
      runs.set(-z * n2 + y * n + x, -z * n2 + y * n - x + 1);
      runs.set(-z * n2 + -y * n + x, -z * n2 + -y * n - x + 1);
      runs.set(y * n2 + z * n + x, y * n2 + z * n - x + 1);
      runs.set(y * n2 + -z * n + x, y * n2 + -z * n - x + 1);
      runs.set(-y * n2 + z * n + x, -y * n2 + z * n - x + 1);
      runs.set(-y * n2 + -z * n + x, -y * n2 + -z * n - x + 1);
      y -= 1;
      d += deltaUp;
      deltaUp += 2;

      while (d > 0) {
        runs.set(z * n2 + x * n + y + 1, z * n2 + x * n - y);
        runs.set(z * n2 + -x * n + y + 1, z * n2 + -x * n - y);
        runs.set(-z * n2 + x * n + y + 1, -z * n2 + x * n - y);
        runs.set(-z * n2 + -x * n + y + 1, -z * n2 + -x * n - y);
        runs.set(x * n2 + z * n + y + 1, x * n2 + z * n - y);
        runs.set(-x * n2 + z * n + y + 1, -x * n2 + z * n - y);
        runs.set(x * n2 + -z * n + y + 1, x * n2 + -z * n - y);
        runs.set(-x * n2 + -z * n + y + 1, -x * n2 + -z * n - y);
        x += 1;
        d += deltaRight;
        deltaRight += 2;
      }
    }

    if (y === x) {
      runs.set(z * n2 + y * n + x, z * n2 + y * n - x + 1);
      runs.set(z * n2 + -y * n + x, z * n2 + -y * n - x + 1);
      runs.set(-z * n2 + y * n + x, -z * n2 + y * n - x + 1);
      runs.set(-z * n2 + -y * n + x, -z * n2 + -y * n - x + 1);
      runs.set(y * n2 + z * n + x, y * n2 + z * n - x + 1);
      runs.set(y * n2 + -z * n + x, y * n2 + -z * n - x + 1);
      runs.set(-y * n2 + z * n + x, -y * n2 + z * n - x + 1);
      runs.set(-y * n2 + -z * n + x, -y * n2 + -z * n - x + 1);
    }
  };

  while (z > x0) {
    innerLoop(x0, y0, d0, deltaUp0, deltaRight0);
    z -= 1;
    d0 += deltaOut;
    deltaOut += 2;

    while (d0 > 0) {
      x0 += 1;
      d0 += deltaRight0;
      deltaRight0 += 2;
    }
  }

  if (z === x0) {
    innerLoop(x0, y0, d0, deltaUp0, deltaRight0);
  }

  const indices: number[] = [];
  runs.forEach((value, key) => {
    indices.push(key, value);
  });
  let lastIndex = -1;
  const sortedIndices: number[] = [];
  indices
    .sort((a, b) => a - b)
    .forEach(index => {
      const properIndex = index + indexOffset;
      if (properIndex === lastIndex) {
        sortedIndices.pop();
        lastIndex = -1;
      } else {
        sortedIndices.push(properIndex);
        lastIndex = properIndex;
      }
    });
  return sortedIndices;
}