scry/lib/recognition/clahe.ts

/**
 * CLAHE (Contrast Limited Adaptive Histogram Equalization) implementation.
 * 
 * This is a pure TypeScript implementation for when OpenCV is not available.
 * For better performance, use OpenCV's createCLAHE when available.
 */

interface CLAHEOptions {
  clipLimit?: number;
  tileGridSize?: { width: number; height: number };
}

const DEFAULT_OPTIONS: Required<CLAHEOptions> = {
  clipLimit: 2.0,
  tileGridSize: { width: 8, height: 8 },
};

/**
 * Convert RGB to LAB color space.
 * Returns L channel (0-100 range, scaled to 0-255 for processing).
 */
function rgbToLab(r: number, g: number, b: number): { l: number; a: number; b: number } {
  // Normalize RGB to 0-1
  let rNorm = r / 255;
  let gNorm = g / 255;
  let bNorm = b / 255;

  // Apply gamma correction
  rNorm = rNorm > 0.04045 ? Math.pow((rNorm + 0.055) / 1.055, 2.4) : rNorm / 12.92;
  gNorm = gNorm > 0.04045 ? Math.pow((gNorm + 0.055) / 1.055, 2.4) : gNorm / 12.92;
  bNorm = bNorm > 0.04045 ? Math.pow((bNorm + 0.055) / 1.055, 2.4) : bNorm / 12.92;

  // Convert to XYZ
  const x = (rNorm * 0.4124564 + gNorm * 0.3575761 + bNorm * 0.1804375) / 0.95047;
  const y = rNorm * 0.2126729 + gNorm * 0.7151522 + bNorm * 0.0721750;
  const z = (rNorm * 0.0193339 + gNorm * 0.1191920 + bNorm * 0.9503041) / 1.08883;

  // Convert to LAB
  const xNorm = x > 0.008856 ? Math.pow(x, 1 / 3) : 7.787 * x + 16 / 116;
  const yNorm = y > 0.008856 ? Math.pow(y, 1 / 3) : 7.787 * y + 16 / 116;
  const zNorm = z > 0.008856 ? Math.pow(z, 1 / 3) : 7.787 * z + 16 / 116;

  return {
    l: Math.max(0, 116 * yNorm - 16), // 0-100
    a: 500 * (xNorm - yNorm),
    b: 200 * (yNorm - zNorm),
  };
}

/**
 * Convert LAB to RGB color space.
 */
function labToRgb(l: number, a: number, b: number): { r: number; g: number; b: number } {
  // Convert LAB to XYZ
  const yNorm = (l + 16) / 116;
  const xNorm = a / 500 + yNorm;
  const zNorm = yNorm - b / 200;

  const x3 = Math.pow(xNorm, 3);
  const y3 = Math.pow(yNorm, 3);
  const z3 = Math.pow(zNorm, 3);

  const x = (x3 > 0.008856 ? x3 : (xNorm - 16 / 116) / 7.787) * 0.95047;
  const y = y3 > 0.008856 ? y3 : (yNorm - 16 / 116) / 7.787;
  const z = (z3 > 0.008856 ? z3 : (zNorm - 16 / 116) / 7.787) * 1.08883;

  // Convert XYZ to RGB
  let rNorm = x * 3.2404542 + y * -1.5371385 + z * -0.4985314;
  let gNorm = x * -0.9692660 + y * 1.8760108 + z * 0.0415560;
  let bNorm = x * 0.0556434 + y * -0.2040259 + z * 1.0572252;

  // Apply gamma correction
  rNorm = rNorm > 0.0031308 ? 1.055 * Math.pow(rNorm, 1 / 2.4) - 0.055 : 12.92 * rNorm;
  gNorm = gNorm > 0.0031308 ? 1.055 * Math.pow(gNorm, 1 / 2.4) - 0.055 : 12.92 * gNorm;
  bNorm = bNorm > 0.0031308 ? 1.055 * Math.pow(bNorm, 1 / 2.4) - 0.055 : 12.92 * bNorm;

  return {
    r: Math.round(Math.max(0, Math.min(255, rNorm * 255))),
    g: Math.round(Math.max(0, Math.min(255, gNorm * 255))),
    b: Math.round(Math.max(0, Math.min(255, bNorm * 255))),
  };
}

/**
 * Compute histogram for a region.
 */
function computeHistogram(
  data: Uint8Array,
  startX: number,
  startY: number,
  width: number,
  height: number,
  stride: number
): number[] {
  const histogram = new Array(256).fill(0);
  
  for (let y = startY; y < startY + height; y++) {
    for (let x = startX; x < startX + width; x++) {
      const value = data[y * stride + x];
      histogram[value]++;
    }
  }
  
  return histogram;
}

/**
 * Clip histogram and redistribute excess.
 */
function clipHistogram(histogram: number[], clipLimit: number, numPixels: number): number[] {
  const limit = Math.floor(clipLimit * numPixels / 256);
  const clipped = [...histogram];
  
  let excess = 0;
  for (let i = 0; i < 256; i++) {
    if (clipped[i] > limit) {
      excess += clipped[i] - limit;
      clipped[i] = limit;
    }
  }
  
  // Redistribute excess
  const increment = Math.floor(excess / 256);
  const remainder = excess % 256;
  
  for (let i = 0; i < 256; i++) {
    clipped[i] += increment;
    if (i < remainder) {
      clipped[i]++;
    }
  }
  
  return clipped;
}

/**
 * Create lookup table from clipped histogram.
 */
function createLUT(histogram: number[], numPixels: number): Uint8Array {
  const lut = new Uint8Array(256);
  let cumSum = 0;
  
  for (let i = 0; i < 256; i++) {
    cumSum += histogram[i];
    lut[i] = Math.round((cumSum / numPixels) * 255);
  }
  
  return lut;
}

/**
 * Bilinear interpolation between four values.
 */
function bilinearInterpolate(
  topLeft: number,
  topRight: number,
  bottomLeft: number,
  bottomRight: number,
  xFrac: number,
  yFrac: number
): number {
  const top = topLeft + (topRight - topLeft) * xFrac;
  const bottom = bottomLeft + (bottomRight - bottomLeft) * xFrac;
  return Math.round(top + (bottom - top) * yFrac);
}

/**
 * Apply CLAHE to RGBA pixel data.
 * Processes only the L channel in LAB color space.
 * 
 * @param pixels RGBA pixel data
 * @param width Image width
 * @param height Image height
 * @param options CLAHE options
 * @returns Processed RGBA pixel data
 */
export function applyCLAHE(
  pixels: Uint8Array | Uint8ClampedArray,
  width: number,
  height: number,
  options: CLAHEOptions = {}
): Uint8Array {
  const opts = { ...DEFAULT_OPTIONS, ...options };
  const { clipLimit, tileGridSize } = opts;
  
  const tileWidth = Math.floor(width / tileGridSize.width);
  const tileHeight = Math.floor(height / tileGridSize.height);
  
  // Extract L channel from LAB
  const lChannel = new Uint8Array(width * height);
  const aChannel = new Float32Array(width * height);
  const bChannel = new Float32Array(width * height);
  
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const idx = (y * width + x) * 4;
      const lab = rgbToLab(pixels[idx], pixels[idx + 1], pixels[idx + 2]);
      const pixelIdx = y * width + x;
      lChannel[pixelIdx] = Math.round((lab.l / 100) * 255);
      aChannel[pixelIdx] = lab.a;
      bChannel[pixelIdx] = lab.b;
    }
  }
  
  // Compute LUTs for each tile
  const luts: Uint8Array[][] = [];
  const numPixelsPerTile = tileWidth * tileHeight;
  
  for (let ty = 0; ty < tileGridSize.height; ty++) {
    luts[ty] = [];
    for (let tx = 0; tx < tileGridSize.width; tx++) {
      const startX = tx * tileWidth;
      const startY = ty * tileHeight;
      
      const histogram = computeHistogram(
        lChannel,
        startX,
        startY,
        tileWidth,
        tileHeight,
        width
      );
      
      const clipped = clipHistogram(histogram, clipLimit, numPixelsPerTile);
      luts[ty][tx] = createLUT(clipped, numPixelsPerTile);
    }
  }
  
  // Apply CLAHE with bilinear interpolation
  const result = new Uint8Array(pixels.length);
  
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      const pixelIdx = y * width + x;
      const rgbaIdx = pixelIdx * 4;
      
      // Determine tile position
      const tileX = Math.min(x / tileWidth, tileGridSize.width - 1);
      const tileY = Math.min(y / tileHeight, tileGridSize.height - 1);
      
      const tx = Math.floor(tileX);
      const ty = Math.floor(tileY);
      
      const xFrac = tileX - tx;
      const yFrac = tileY - ty;
      
      // Get surrounding tiles (handle edges)
      const tx1 = Math.min(tx + 1, tileGridSize.width - 1);
      const ty1 = Math.min(ty + 1, tileGridSize.height - 1);
      
      const lValue = lChannel[pixelIdx];
      
      // Interpolate LUT values
      const newL = bilinearInterpolate(
        luts[ty][tx][lValue],
        luts[ty][tx1][lValue],
        luts[ty1][tx][lValue],
        luts[ty1][tx1][lValue],
        xFrac,
        yFrac
      );
      
      // Convert back to RGB
      const rgb = labToRgb(
        (newL / 255) * 100,
        aChannel[pixelIdx],
        bChannel[pixelIdx]
      );
      
      result[rgbaIdx] = rgb.r;
      result[rgbaIdx + 1] = rgb.g;
      result[rgbaIdx + 2] = rgb.b;
      result[rgbaIdx + 3] = pixels[rgbaIdx + 3]; // Preserve alpha
    }
  }
  
  return result;
}