TV 正则化：ADMM 与 HQS 图像去噪

function tvAdmm(noisy, width, height, lambda, iterations) {
  const rho = 1.5;
  let x = noisy.slice();
  let zx = zerosLike(x), zy = zerosLike(x);
  let ux = zerosLike(x), uy = zerosLike(x);

  for (let k = 0; k < iterations; k += 1) {
    // x-update: (I + rho D^T D)x = b + rho D^T(z - u)
    const rhs = add(noisy, scale(adjoint(zx, zy, ux, uy), rho));
    x = conjugateGradient(x, rhs, rho);

    // z-update: isotropic vector soft-thresholding
    const [gx, gy] = gradient(x, width, height);
    [zx, zy] = shrink2D(gx, gy, ux, uy, lambda / rho);

    // scaled dual update
    ux = add(ux, subtract(gx, zx));
    uy = add(uy, subtract(gy, zy));
  }

  return x;
}

function tvHqs(noisy, width, height, lambda, stages) {
  let x = noisy.slice();
  let beta = lambda / 2;
  const B = fft2(noisy);
  const { dx, dy, dtd } = periodicDifferenceSpectrum(width, height);

  for (let k = 0; k < stages; k += 1) {
    const X = fft2(x);
    let gx = real(ifft2(X * dx));
    let gy = real(ifft2(X * dy));

    // Component-wise soft-thresholding from tv_reg.m.
    gx = soft(gx, lambda / beta);
    gy = soft(gy, lambda / beta);

    const numerator =
      B + beta * (fft2(gx) * conj(dx) + fft2(gy) * conj(dy));
    x = real(ifft2(numerator / (1 + beta * dtd)));
    beta *= 2;
  }

  return x;
}