scry/lib/recognition/perspectiveCorrection.ts

/**
 * Perspective correction to transform a quadrilateral region into a rectangle.
 * Uses homography transformation with bilinear sampling.
 */

import { Point, sampleBilinear } from "./imageUtils";

/** Standard output size for corrected card images (63:88 aspect ratio). */
export const OUTPUT_WIDTH = 480;
export const OUTPUT_HEIGHT = 670;

/**
 * Apply perspective correction to extract and normalize a card from an image.
 * 
 * @param pixels - Source RGBA pixel data
 * @param width - Source image width
 * @param height - Source image height
 * @param corners - Four corners in order: top-left, top-right, bottom-right, bottom-left
 * @param outputWidth - Width of output image
 * @param outputHeight - Height of output image
 * @returns Perspective-corrected RGBA pixel data
 */
export function warpPerspective(
  pixels: Uint8Array | Uint8ClampedArray,
  width: number,
  height: number,
  corners: Point[],
  outputWidth: number = OUTPUT_WIDTH,
  outputHeight: number = OUTPUT_HEIGHT
): { pixels: Uint8Array; width: number; height: number } {
  if (corners.length !== 4) {
    throw new Error("Exactly 4 corners required");
  }
  
  // Determine if card is landscape (rotated 90°)
  const width1 = distance(corners[0], corners[1]);
  const height1 = distance(corners[1], corners[2]);
  
  let orderedCorners: Point[];
  
  if (width1 > height1) {
    // Card is landscape - rotate corners to portrait
    orderedCorners = [corners[1], corners[2], corners[3], corners[0]];
  } else {
    orderedCorners = corners;
  }
  
  // Compute the perspective transform matrix
  const matrix = computePerspectiveTransform(orderedCorners, outputWidth, outputHeight);
  const inverse = invertMatrix3x3(matrix);
  
  // Apply the transform
  const result = new Uint8Array(outputWidth * outputHeight * 4);
  
  for (let y = 0; y < outputHeight; y++) {
    for (let x = 0; x < outputWidth; x++) {
      // Apply inverse transform to find source coordinates
      const srcPoint = applyTransform(inverse, x, y);
      
      // Bilinear interpolation for smooth sampling
      const [r, g, b, a] = sampleBilinear(pixels, width, height, srcPoint.x, srcPoint.y);
      
      const idx = (y * outputWidth + x) * 4;
      result[idx] = r;
      result[idx + 1] = g;
      result[idx + 2] = b;
      result[idx + 3] = a;
    }
  }
  
  return { pixels: result, width: outputWidth, height: outputHeight };
}

/**
 * Compute a perspective transform matrix from quad corners to rectangle.
 * Uses the Direct Linear Transform (DLT) algorithm.
 */
function computePerspectiveTransform(
  src: Point[],
  dstWidth: number,
  dstHeight: number
): number[] {
  // Destination corners (rectangle)
  const dst: Point[] = [
    { x: 0, y: 0 },
    { x: dstWidth - 1, y: 0 },
    { x: dstWidth - 1, y: dstHeight - 1 },
    { x: 0, y: dstHeight - 1 },
  ];
  
  // Build the 8x8 matrix for solving the homography
  const A: number[][] = [];
  const b: number[] = [];
  
  for (let i = 0; i < 4; i++) {
    const sx = src[i].x;
    const sy = src[i].y;
    const dx = dst[i].x;
    const dy = dst[i].y;
    
    A.push([sx, sy, 1, 0, 0, 0, -dx * sx, -dx * sy]);
    b.push(dx);
    
    A.push([0, 0, 0, sx, sy, 1, -dy * sx, -dy * sy]);
    b.push(dy);
  }
  
  // Solve using Gaussian elimination
  const h = solveLinearSystem(A, b);
  
  // Return 3x3 matrix as flat array [h11, h12, h13, h21, h22, h23, h31, h32, h33]
  return [h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], 1];
}

/**
 * Solve a linear system Ax = b using Gaussian elimination with partial pivoting.
 */
function solveLinearSystem(A: number[][], b: number[]): number[] {
  const n = b.length;
  
  // Create augmented matrix
  const augmented: number[][] = A.map((row, i) => [...row, b[i]]);
  
  // Forward elimination with partial pivoting
  for (let col = 0; col < n; col++) {
    // Find pivot
    let maxRow = col;
    for (let row = col + 1; row < n; row++) {
      if (Math.abs(augmented[row][col]) > Math.abs(augmented[maxRow][col])) {
        maxRow = row;
      }
    }
    
    // Swap rows
    [augmented[col], augmented[maxRow]] = [augmented[maxRow], augmented[col]];
    
    // Eliminate
    for (let row = col + 1; row < n; row++) {
      if (Math.abs(augmented[col][col]) < 1e-10) continue;
      
      const factor = augmented[row][col] / augmented[col][col];
      for (let j = col; j <= n; j++) {
        augmented[row][j] -= factor * augmented[col][j];
      }
    }
  }
  
  // Back substitution
  const x = new Array(n).fill(0);
  for (let i = n - 1; i >= 0; i--) {
    x[i] = augmented[i][n];
    for (let j = i + 1; j < n; j++) {
      x[i] -= augmented[i][j] * x[j];
    }
    if (Math.abs(augmented[i][i]) > 1e-10) {
      x[i] /= augmented[i][i];
    }
  }
  
  return x;
}

/**
 * Apply a 3x3 transform matrix to a point.
 */
function applyTransform(H: number[], x: number, y: number): Point {
  let w = H[6] * x + H[7] * y + H[8];
  if (Math.abs(w) < 1e-10) w = 1e-10;
  
  return {
    x: (H[0] * x + H[1] * y + H[2]) / w,
    y: (H[3] * x + H[4] * y + H[5]) / w,
  };
}

/**
 * Invert a 3x3 matrix.
 */
function invertMatrix3x3(m: number[]): number[] {
  const det =
    m[0] * (m[4] * m[8] - m[5] * m[7]) -
    m[1] * (m[3] * m[8] - m[5] * m[6]) +
    m[2] * (m[3] * m[7] - m[4] * m[6]);
  
  const invDet = Math.abs(det) < 1e-10 ? 1e10 : 1 / det;
  
  return [
    (m[4] * m[8] - m[5] * m[7]) * invDet,
    (m[2] * m[7] - m[1] * m[8]) * invDet,
    (m[1] * m[5] - m[2] * m[4]) * invDet,
    (m[5] * m[6] - m[3] * m[8]) * invDet,
    (m[0] * m[8] - m[2] * m[6]) * invDet,
    (m[2] * m[3] - m[0] * m[5]) * invDet,
    (m[3] * m[7] - m[4] * m[6]) * invDet,
    (m[1] * m[6] - m[0] * m[7]) * invDet,
    (m[0] * m[4] - m[1] * m[3]) * invDet,
  ];
}

/**
 * Calculate distance between two points.
 */
function distance(a: Point, b: Point): number {
  const dx = b.x - a.x;
  const dy = b.y - a.y;
  return Math.sqrt(dx * dx + dy * dy);
}