void CLASS foveon_interpolate()
{
  static const short hood[] = {-1, -1, -1, 0, -1, 1, 0, -1, 0, 1, 1, -1, 1, 0, 1, 1};
  short *pix, prev[3], *curve[8], (*shrink)[3];
  float cfilt = 0, ddft[3][3][2], ppm[3][3][3];
  float cam_xyz[3][3], correct[3][3], last[3][3], trans[3][3];
  float chroma_dq[3], color_dq[3], diag[3][3], div[3];
  float(*black)[3], (*sgain)[3], (*sgrow)[3];
  float fsum[3], val, frow, num;
  int row, col, c, i, j, diff, sgx, irow, sum, min, max, limit;
  int dscr[2][2], dstb[4], (*smrow[7])[3], total[4], ipix[3];
  int work[3][3], smlast, smred, smred_p = 0, dev[3];
  int satlev[3], keep[4], active[4];
  unsigned dim[3], *badpix;
  double dsum = 0, trsum[3];
  char str[128];
  const char *cp;

#ifdef DCRAW_VERBOSE
  if (verbose)
    fprintf(stderr, _("Foveon interpolation...\n"));
#endif

  foveon_load_camf();
  foveon_fixed(dscr, 4, "DarkShieldColRange");
  foveon_fixed(ppm[0][0], 27, "PostPolyMatrix");
  foveon_fixed(satlev, 3, "SaturationLevel");
  foveon_fixed(keep, 4, "KeepImageArea");
  foveon_fixed(active, 4, "ActiveImageArea");
  foveon_fixed(chroma_dq, 3, "ChromaDQ");
  foveon_fixed(color_dq, 3, foveon_camf_param("IncludeBlocks", "ColorDQ") ? "ColorDQ" : "ColorDQCamRGB");
  if (foveon_camf_param("IncludeBlocks", "ColumnFilter"))
    foveon_fixed(&cfilt, 1, "ColumnFilter");

  memset(ddft, 0, sizeof ddft);
  if (!foveon_camf_param("IncludeBlocks", "DarkDrift") || !foveon_fixed(ddft[1][0], 12, "DarkDrift"))
    for (i = 0; i < 2; i++)
    {
      foveon_fixed(dstb, 4, i ? "DarkShieldBottom" : "DarkShieldTop");
      for (row = dstb[1]; row <= dstb[3]; row++)
        for (col = dstb[0]; col <= dstb[2]; col++)
          FORC3 ddft[i + 1][c][1] += (short)image[row * width + col][c];
      FORC3 ddft[i + 1][c][1] /= (dstb[3] - dstb[1] + 1) * (dstb[2] - dstb[0] + 1);
    }

  if (!(cp = foveon_camf_param("WhiteBalanceIlluminants", model2)))
  {
#ifdef DCRAW_VERBOSE
    fprintf(stderr, _("%s: Invalid white balance \"%s\"\n"), ifname, model2);
#endif
    return;
  }
  foveon_fixed(cam_xyz, 9, cp);
  foveon_fixed(correct, 9, foveon_camf_param("WhiteBalanceCorrections", model2));
  memset(last, 0, sizeof last);
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      FORC3 last[i][j] += correct[i][c] * cam_xyz[c][j];

#define LAST(x, y) last[(i + x) % 3][(c + y) % 3]
  for (i = 0; i < 3; i++)
    FORC3 diag[c][i] = LAST(1, 1) * LAST(2, 2) - LAST(1, 2) * LAST(2, 1);
#undef LAST
  FORC3 div[c] = diag[c][0] * 0.3127 + diag[c][1] * 0.329 + diag[c][2] * 0.3583;
  sprintf(str, "%sRGBNeutral", model2);
  if (foveon_camf_param("IncludeBlocks", str))
    foveon_fixed(div, 3, str);
  num = 0;
  FORC3 if (num < div[c]) num = div[c];
  FORC3 div[c] /= num;

  memset(trans, 0, sizeof trans);
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      FORC3 trans[i][j] += rgb_cam[i][c] * last[c][j] * div[j];
  FORC3 trsum[c] = trans[c][0] + trans[c][1] + trans[c][2];
  dsum = (6 * trsum[0] + 11 * trsum[1] + 3 * trsum[2]) / 20;
  for (i = 0; i < 3; i++)
    FORC3 last[i][c] = trans[i][c] * dsum / trsum[i];
  memset(trans, 0, sizeof trans);
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      FORC3 trans[i][j] += (i == c ? 32 : -1) * last[c][j] / 30;

  foveon_make_curves(curve, color_dq, div, cfilt);
  FORC3 chroma_dq[c] /= 3;
  foveon_make_curves(curve + 3, chroma_dq, div, cfilt);
  FORC3 dsum += chroma_dq[c] / div[c];
  curve[6] = foveon_make_curve(dsum, dsum, cfilt);
  curve[7] = foveon_make_curve(dsum * 2, dsum * 2, cfilt);

  sgain = (float(*)[3])foveon_camf_matrix(dim, "SpatialGain");
  if (!sgain)
    return;
  sgrow = (float(*)[3])calloc(dim[1], sizeof *sgrow);
  sgx = (width + dim[1] - 2) / (dim[1] - 1);

  black = (float(*)[3])calloc(height, sizeof *black);
  for (row = 0; row < height; row++)
  {
    for (i = 0; i < 6; i++)
      ((float *)ddft[0])[i] =
          ((float *)ddft[1])[i] + row / (height - 1.0) * (((float *)ddft[2])[i] - ((float *)ddft[1])[i]);
    FORC3 black[row][c] = (foveon_avg(image[row * width] + c, dscr[0], cfilt) +
                           foveon_avg(image[row * width] + c, dscr[1], cfilt) * 3 - ddft[0][c][0]) /
                              4 -
                          ddft[0][c][1];
  }
  memcpy(black, black + 8, sizeof *black * 8);
  memcpy(black + height - 11, black + height - 22, 11 * sizeof *black);
  memcpy(last, black, sizeof last);

  for (row = 1; row < height - 1; row++)
  {
    FORC3 if (last[1][c] > last[0][c])
    {
      if (last[1][c] > last[2][c])
        black[row][c] = (last[0][c] > last[2][c]) ? last[0][c] : last[2][c];
    }
    else if (last[1][c] < last[2][c]) black[row][c] = (last[0][c] < last[2][c]) ? last[0][c] : last[2][c];
    memmove(last, last + 1, 2 * sizeof last[0]);
    memcpy(last[2], black[row + 1], sizeof last[2]);
  }
  FORC3 black[row][c] = (last[0][c] + last[1][c]) / 2;
  FORC3 black[0][c] = (black[1][c] + black[3][c]) / 2;

  val = 1 - exp(-1 / 24.0);
  memcpy(fsum, black, sizeof fsum);
  for (row = 1; row < height; row++)
    FORC3 fsum[c] += black[row][c] = (black[row][c] - black[row - 1][c]) * val + black[row - 1][c];
  memcpy(last[0], black[height - 1], sizeof last[0]);
  FORC3 fsum[c] /= height;
  for (row = height; row--;)
    FORC3 last[0][c] = black[row][c] = (black[row][c] - fsum[c] - last[0][c]) * val + last[0][c];

  memset(total, 0, sizeof total);
  for (row = 2; row < height; row += 4)
    for (col = 2; col < width; col += 4)
    {
      FORC3 total[c] += (short)image[row * width + col][c];
      total[3]++;
    }
  for (row = 0; row < height; row++)
    FORC3 black[row][c] += fsum[c] / 2 + total[c] / (total[3] * 100.0);

  for (row = 0; row < height; row++)
  {
    for (i = 0; i < 6; i++)
      ((float *)ddft[0])[i] =
          ((float *)ddft[1])[i] + row / (height - 1.0) * (((float *)ddft[2])[i] - ((float *)ddft[1])[i]);
    pix = image[row * width];
    memcpy(prev, pix, sizeof prev);
    frow = row / (height - 1.0) * (dim[2] - 1);
    if ((irow = frow) == dim[2] - 1)
      irow--;
    frow -= irow;
    for (i = 0; i < dim[1]; i++)
      FORC3 sgrow[i][c] = sgain[irow * dim[1] + i][c] * (1 - frow) + sgain[(irow + 1) * dim[1] + i][c] * frow;
    for (col = 0; col < width; col++)
    {
      FORC3
      {
        diff = pix[c] - prev[c];
        prev[c] = pix[c];
        ipix[c] = pix[c] + floor((diff + (diff * diff >> 14)) * cfilt - ddft[0][c][1] -
                                 ddft[0][c][0] * ((float)col / width - 0.5) - black[row][c]);
      }
      FORC3
      {
        work[0][c] = ipix[c] * ipix[c] >> 14;
        work[2][c] = ipix[c] * work[0][c] >> 14;
        work[1][2 - c] = ipix[(c + 1) % 3] * ipix[(c + 2) % 3] >> 14;
      }
      FORC3
      {
        for (val = i = 0; i < 3; i++)
          for (j = 0; j < 3; j++)
            val += ppm[c][i][j] * work[i][j];
        ipix[c] =
            floor((ipix[c] + floor(val)) *
                  (sgrow[col / sgx][c] * (sgx - col % sgx) + sgrow[col / sgx + 1][c] * (col % sgx)) / sgx / div[c]);
        if (ipix[c] > 32000)
          ipix[c] = 32000;
        pix[c] = ipix[c];
      }
      pix += 4;
    }
  }
  free(black);
  free(sgrow);
  free(sgain);

  if ((badpix = (unsigned *)foveon_camf_matrix(dim, "BadPixels")))
  {
    for (i = 0; i < dim[0]; i++)
    {
      col = (badpix[i] >> 8 & 0xfff) - keep[0];
      row = (badpix[i] >> 20) - keep[1];
      if ((unsigned)(row - 1) > height - 3 || (unsigned)(col - 1) > width - 3)
        continue;
      memset(fsum, 0, sizeof fsum);
      for (sum = j = 0; j < 8; j++)
        if (badpix[i] & (1 << j))
        {
          FORC3 fsum[c] += (short)image[(row + hood[j * 2]) * width + col + hood[j * 2 + 1]][c];
          sum++;
        }
      if (sum)
        FORC3 image[row * width + col][c] = fsum[c] / sum;
    }
    free(badpix);
  }

  /* Array for 5x5 Gaussian averaging of red values */
  smrow[6] = (int(*)[3])calloc(width * 5, sizeof **smrow);
  merror(smrow[6], "foveon_interpolate()");
  for (i = 0; i < 5; i++)
    smrow[i] = smrow[6] + i * width;

  /* Sharpen the reds against these Gaussian averages */
  for (smlast = -1, row = 2; row < height - 2; row++)
  {
    while (smlast < row + 2)
    {
      for (i = 0; i < 6; i++)
        smrow[(i + 5) % 6] = smrow[i];
      pix = image[++smlast * width + 2];
      for (col = 2; col < width - 2; col++)
      {
        smrow[4][col][0] = (pix[0] * 6 + (pix[-4] + pix[4]) * 4 + pix[-8] + pix[8] + 8) >> 4;
        pix += 4;
      }
    }
    pix = image[row * width + 2];
    for (col = 2; col < width - 2; col++)
    {
      smred = (6 * smrow[2][col][0] + 4 * (smrow[1][col][0] + smrow[3][col][0]) + smrow[0][col][0] + smrow[4][col][0] +
               8) >>
              4;
      if (col == 2)
        smred_p = smred;
      i = pix[0] + ((pix[0] - ((smred * 7 + smred_p) >> 3)) >> 3);
      if (i > 32000)
        i = 32000;
      pix[0] = i;
      smred_p = smred;
      pix += 4;
    }
  }

  /* Adjust the brighter pixels for better linearity */
  min = 0xffff;
  FORC3
  {
    i = satlev[c] / div[c];
    if (min > i)
      min = i;
  }
  limit = min * 9 >> 4;
  for (pix = image[0]; pix < image[height * width]; pix += 4)
  {
    if (pix[0] <= limit || pix[1] <= limit || pix[2] <= limit)
      continue;
    min = max = pix[0];
    for (c = 1; c < 3; c++)
    {
      if (min > pix[c])
        min = pix[c];
      if (max < pix[c])
        max = pix[c];
    }
    if (min >= limit * 2)
    {
      pix[0] = pix[1] = pix[2] = max;
    }
    else
    {
      i = 0x4000 - ((min - limit) << 14) / limit;
      i = 0x4000 - (i * i >> 14);
      i = i * i >> 14;
      FORC3 pix[c] += (max - pix[c]) * i >> 14;
    }
  }
  /*
     Because photons that miss one detector often hit another,
     the sum R+G+B is much less noisy than the individual colors.
     So smooth the hues without smoothing the total.
   */
  for (smlast = -1, row = 2; row < height - 2; row++)
  {
    while (smlast < row + 2)
    {
      for (i = 0; i < 6; i++)
        smrow[(i + 5) % 6] = smrow[i];
      pix = image[++smlast * width + 2];
      for (col = 2; col < width - 2; col++)
      {
        FORC3 smrow[4][col][c] = (pix[c - 4] + 2 * pix[c] + pix[c + 4] + 2) >> 2;
        pix += 4;
      }
    }
    pix = image[row * width + 2];
    for (col = 2; col < width - 2; col++)
    {
      FORC3 dev[c] =
          -foveon_apply_curve(curve[7], pix[c] - ((smrow[1][col][c] + 2 * smrow[2][col][c] + smrow[3][col][c]) >> 2));
      sum = (dev[0] + dev[1] + dev[2]) >> 3;
      FORC3 pix[c] += dev[c] - sum;
      pix += 4;
    }
  }
  for (smlast = -1, row = 2; row < height - 2; row++)
  {
    while (smlast < row + 2)
    {
      for (i = 0; i < 6; i++)
        smrow[(i + 5) % 6] = smrow[i];
      pix = image[++smlast * width + 2];
      for (col = 2; col < width - 2; col++)
      {
        FORC3 smrow[4][col][c] = (pix[c - 8] + pix[c - 4] + pix[c] + pix[c + 4] + pix[c + 8] + 2) >> 2;
        pix += 4;
      }
    }
    pix = image[row * width + 2];
    for (col = 2; col < width - 2; col++)
    {
      for (total[3] = 375, sum = 60, c = 0; c < 3; c++)
      {
        for (total[c] = i = 0; i < 5; i++)
          total[c] += smrow[i][col][c];
        total[3] += total[c];
        sum += pix[c];
      }
      if (sum < 0)
        sum = 0;
      j = total[3] > 375 ? (sum << 16) / total[3] : sum * 174;
      FORC3 pix[c] += foveon_apply_curve(curve[6], ((j * total[c] + 0x8000) >> 16) - pix[c]);
      pix += 4;
    }
  }

  /* Transform the image to a different colorspace */
  for (pix = image[0]; pix < image[height * width]; pix += 4)
  {
    FORC3 pix[c] -= foveon_apply_curve(curve[c], pix[c]);
    sum = (pix[0] + pix[1] + pix[1] + pix[2]) >> 2;
    FORC3 pix[c] -= foveon_apply_curve(curve[c], pix[c] - sum);
    FORC3
    {
      for (dsum = i = 0; i < 3; i++)
        dsum += trans[c][i] * pix[i];
      if (dsum < 0)
        dsum = 0;
      if (dsum > 24000)
        dsum = 24000;
      ipix[c] = dsum + 0.5;
    }
    FORC3 pix[c] = ipix[c];
  }

  /* Smooth the image bottom-to-top and save at 1/4 scale */
  shrink = (short(*)[3])calloc((height / 4), (width / 4) * sizeof *shrink);
  merror(shrink, "foveon_interpolate()");
  for (row = height / 4; row--;)
    for (col = 0; col < width / 4; col++)
    {
      ipix[0] = ipix[1] = ipix[2] = 0;
      for (i = 0; i < 4; i++)
        for (j = 0; j < 4; j++)
          FORC3 ipix[c] += image[(row * 4 + i) * width + col * 4 + j][c];
      FORC3
      if (row + 2 > height / 4)
        shrink[row * (width / 4) + col][c] = ipix[c] >> 4;
      else
        shrink[row * (width / 4) + col][c] =
            (shrink[(row + 1) * (width / 4) + col][c] * 1840 + ipix[c] * 141 + 2048) >> 12;
    }
  /* From the 1/4-scale image, smooth right-to-left */
  for (row = 0; row < (height & ~3); row++)
  {
    ipix[0] = ipix[1] = ipix[2] = 0;
    if ((row & 3) == 0)
      for (col = width & ~3; col--;)
        FORC3 smrow[0][col][c] = ipix[c] =
            (shrink[(row / 4) * (width / 4) + col / 4][c] * 1485 + ipix[c] * 6707 + 4096) >> 13;

    /* Then smooth left-to-right */
    ipix[0] = ipix[1] = ipix[2] = 0;
    for (col = 0; col < (width & ~3); col++)
      FORC3 smrow[1][col][c] = ipix[c] = (smrow[0][col][c] * 1485 + ipix[c] * 6707 + 4096) >> 13;

    /* Smooth top-to-bottom */
    if (row == 0)
      memcpy(smrow[2], smrow[1], sizeof **smrow * width);
    else
      for (col = 0; col < (width & ~3); col++)
        FORC3 smrow[2][col][c] = (smrow[2][col][c] * 6707 + smrow[1][col][c] * 1485 + 4096) >> 13;

    /* Adjust the chroma toward the smooth values */
    for (col = 0; col < (width & ~3); col++)
    {
      for (i = j = 30, c = 0; c < 3; c++)
      {
        i += smrow[2][col][c];
        j += image[row * width + col][c];
      }
      j = (j << 16) / i;
      for (sum = c = 0; c < 3; c++)
      {
        ipix[c] =
            foveon_apply_curve(curve[c + 3], ((smrow[2][col][c] * j + 0x8000) >> 16) - image[row * width + col][c]);
        sum += ipix[c];
      }
      sum >>= 3;
      FORC3
      {
        i = image[row * width + col][c] + ipix[c] - sum;
        if (i < 0)
          i = 0;
        image[row * width + col][c] = i;
      }
    }
  }
  free(shrink);
  free(smrow[6]);
  for (i = 0; i < 8; i++)
    free(curve[i]);

  /* Trim off the black border */
  active[1] -= keep[1];
  active[3] -= 2;
  i = active[2] - active[0];
  for (row = 0; row < active[3] - active[1]; row++)
    memcpy(image[row * i], image[(row + active[1]) * width + active[0]], i * sizeof *image);
  width = i;
  height = row;
}