libs/libkdcraw/libraw/internal/foveon.cpp

   1 /*
   2    GENERATED FILE, DO NOT EDIT
   3    Generated from dcraw/dcraw.c at Thu Jan  8 19:26:36 2009
   4    Look into original file (probably http://cybercom.net/~dcoffin/dcraw/dcraw.c)
   5    for copyright information.
   6 */
   7
   8 #define CLASS LibRaw::
   9 #include "libraw/libraw_types.h"
  10 #define LIBRAW_LIBRARY_BUILD
  11 #include "libraw/libraw.h"
  12 #include "internal/defines.h"
  13 #define SRC_USES_SHRINK
  14 #define SRC_USES_BLACK
  15 #define SRC_USES_CURVE
  16 #include "internal/var_defines.h"
  17 #define sget4(s) sget4((uchar *)s)
  18
  19 /* RESTRICTED code starts here */
  20
  21 void CLASS foveon_decoder (unsigned size, unsigned code)
  22 {
  23 #ifndef LIBRAW_NOTHREADS
  24 #define huff tls->foveon_decoder_huff
  25 #else
  26   static unsigned huff[1024];
  27 #endif
  28   struct decode *cur;
  29   int i, len;
  30
  31   if (!code) {
  32     for (i=0; i < size; i++)
  33       huff[i] = get4();
  34     init_decoder();
  35   }
  36   cur = free_decode++;
  37   if (free_decode > first_decode+2048) {
  38 #ifdef LIBRAW_LIBRARY_BUILD
  39       throw LIBRAW_EXCEPTION_DECODE_RAW;
  40 #else
  41     fprintf (stderr,_("%s: decoder table overflow\n"), ifname);
  42     longjmp (failure, 2);
  43 #endif
  44   }
  45   if (code)
  46     for (i=0; i < size; i++)
  47       if (huff[i] == code) {
  48         cur->leaf = i;
  49         return;
  50       }
  51   if ((len = code >> 27) > 26) return;
  52   code = (len+1) << 27 | (code & 0x3ffffff) << 1;
  53
  54   cur->branch[0] = free_decode;
  55   foveon_decoder (size, code);
  56   cur->branch[1] = free_decode;
  57   foveon_decoder (size, code+1);
  58 #ifndef LIBRAW_NOTHREADS
  59 #undef huff
  60 #endif
  61 }
  62
  63 void CLASS foveon_thumb (FILE *tfp)
  64 {
  65   unsigned bwide, row, col, bitbuf=0, bit=1, c, i;
  66   char *buf;
  67   struct decode *dindex;
  68   short pred[3];
  69
  70   bwide = get4();
  71   fprintf (tfp, "P6\n%d %d\n255\n", thumb_width, thumb_height);
  72   if (bwide > 0) {
  73     if (bwide < thumb_width*3) return;
  74     buf = (char *) malloc (bwide);
  75     merror (buf, "foveon_thumb()");
  76     for (row=0; row < thumb_height; row++) {
  77       fread  (buf, 1, bwide, ifp);
  78       fwrite (buf, 3, thumb_width, tfp);
  79     }
  80     free (buf);
  81     return;
  82   }
  83   foveon_decoder (256, 0);
  84
  85   for (row=0; row < thumb_height; row++) {
  86     memset (pred, 0, sizeof pred);
  87     if (!bit) get4();
  88     for (bit=col=0; col < thumb_width; col++)
  89       FORC3 {
  90         for (dindex=first_decode; dindex->branch[0]; ) {
  91           if ((bit = (bit-1) & 31) == 31)
  92             for (i=0; i < 4; i++)
  93               bitbuf = (bitbuf << 8) + fgetc(ifp);
  94           dindex = dindex->branch[bitbuf >> bit & 1];
  95         }
  96         pred[c] += dindex->leaf;
  97         fputc (pred[c], tfp);
  98       }
  99   }
 100 }
 101
 102 void CLASS foveon_load_camf()
 103 {
 104   unsigned key, i, val;
 105
 106   fseek (ifp, meta_offset, SEEK_SET);
 107   key = get4();
 108   fread (meta_data, 1, meta_length, ifp);
 109   for (i=0; i < meta_length; i++) {
 110     key = (key * 1597 + 51749) % 244944;
 111     val = key * (INT64) 301593171 >> 24;
 112     meta_data[i] ^= ((((key << 8) - val) >> 1) + val) >> 17;
 113   }
 114 }
 115
 116 void CLASS foveon_load_raw()
 117 {
 118   struct decode *dindex;
 119   short diff[1024];
 120   unsigned bitbuf=0;
 121   int pred[3], fixed, row, col, bit=-1, c, i;
 122
 123   fixed = get4();
 124   read_shorts ((ushort *) diff, 1024);
 125   if (!fixed) foveon_decoder (1024, 0);
 126
 127   for (row=0; row < height; row++) {
 128     memset (pred, 0, sizeof pred);
 129     if (!bit && !fixed && atoi(model+2) < 14) get4();
 130     for (col=bit=0; col < width; col++) {
 131       if (fixed) {
 132         bitbuf = get4();
 133         FORC3 pred[2-c] += diff[bitbuf >> c*10 & 0x3ff];
 134       }
 135       else FORC3 {
 136         for (dindex=first_decode; dindex->branch[0]; ) {
 137           if ((bit = (bit-1) & 31) == 31)
 138             for (i=0; i < 4; i++)
 139               bitbuf = (bitbuf << 8) + fgetc(ifp);
 140           dindex = dindex->branch[bitbuf >> bit & 1];
 141         }
 142         pred[c] += diff[dindex->leaf];
 143         if (pred[c] >> 16 && ~pred[c] >> 16) derror();
 144       }
 145       FORC3 image[row*width+col][c] = pred[c];
 146     }
 147   }
 148   if (document_mode)
 149     for (i=0; i < height*width*4; i++)
 150       if ((short) image[0][i] < 0) image[0][i] = 0;
 151   foveon_load_camf();
 152 }
 153
 154 const char * CLASS foveon_camf_param (const char *block, const char *param)
 155 {
 156   unsigned idx, num;
 157   char *pos, *cp, *dp;
 158
 159   for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
 160     pos = meta_data + idx;
 161     if (strncmp (pos, "CMb", 3)) break;
 162     if (pos[3] != 'P') continue;
 163     if (strcmp (block, pos+sget4(pos+12))) continue;
 164     cp = pos + sget4(pos+16);
 165     num = sget4(cp);
 166     dp = pos + sget4(cp+4);
 167     while (num--) {
 168       cp += 8;
 169       if (!strcmp (param, dp+sget4(cp)))
 170         return dp+sget4(cp+4);
 171     }
 172   }
 173   return 0;
 174 }
 175
 176 void * CLASS foveon_camf_matrix (unsigned dim[3], const char *name)
 177 {
 178   unsigned i, idx, type, ndim, size, *mat;
 179   char *pos, *cp, *dp;
 180   double dsize;
 181
 182   for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
 183     pos = meta_data + idx;
 184     if (strncmp (pos, "CMb", 3)) break;
 185     if (pos[3] != 'M') continue;
 186     if (strcmp (name, pos+sget4(pos+12))) continue;
 187     dim[0] = dim[1] = dim[2] = 1;
 188     cp = pos + sget4(pos+16);
 189     type = sget4(cp);
 190     if ((ndim = sget4(cp+4)) > 3) break;
 191     dp = pos + sget4(cp+8);
 192     for (i=ndim; i--; ) {
 193       cp += 12;
 194       dim[i] = sget4(cp);
 195     }
 196     if ((dsize = (double) dim[0]*dim[1]*dim[2]) > meta_length/4) break;
 197     mat = (unsigned *) malloc ((size = dsize) * 4);
 198     merror (mat, "foveon_camf_matrix()");
 199     for (i=0; i < size; i++)
 200       if (type && type != 6)
 201         mat[i] = sget4(dp + i*4);
 202       else
 203         mat[i] = sget4(dp + i*2) & 0xffff;
 204     return mat;
 205   }
 206 #ifdef LIBRAW_LIBRARY_BUILD
 207   imgdata.process_warnings |= LIBRAW_WARN_FOVEON_NOMATRIX;
 208 #endif
 209 #ifdef DCRAW_VERBOSE
 210   fprintf (stderr,_("%s: \"%s\" matrix not found!\n"), ifname, name);
 211 #endif
 212   return 0;
 213 }
 214
 215 int CLASS foveon_fixed (void *ptr, int size, const char *name)
 216 {
 217   void *dp;
 218   unsigned dim[3];
 219
 220   dp = foveon_camf_matrix (dim, name);
 221   if (!dp) return 0;
 222   memcpy (ptr, dp, size*4);
 223   free (dp);
 224   return 1;
 225 }
 226
 227 float CLASS foveon_avg (short *pix, int range[2], float cfilt)
 228 {
 229   int i;
 230   float val, min=FLT_MAX, max=-FLT_MAX, sum=0;
 231
 232   for (i=range[0]; i <= range[1]; i++) {
 233     sum += val = pix[i*4] + (pix[i*4]-pix[(i-1)*4]) * cfilt;
 234     if (min > val) min = val;
 235     if (max < val) max = val;
 236   }
 237   if (range[1] - range[0] == 1) return sum/2;
 238   return (sum - min - max) / (range[1] - range[0] - 1);
 239 }
 240
 241 short * CLASS foveon_make_curve (double max, double mul, double filt)
 242 {
 243   short *curve;
 244   unsigned i, size;
 245   double x;
 246
 247   if (!filt) filt = 0.8;
 248   size = 4*M_PI*max / filt;
 249   if (size == UINT_MAX) size--;
 250   curve = (short *) calloc (size+1, sizeof *curve);
 251   merror (curve, "foveon_make_curve()");
 252   curve[0] = size;
 253   for (i=0; i < size; i++) {
 254     x = i*filt/max/4;
 255     curve[i+1] = (cos(x)+1)/2 * tanh(i*filt/mul) * mul + 0.5;
 256   }
 257   return curve;
 258 }
 259
 260 void CLASS foveon_make_curves
 261         (short **curvep, float dq[3], float div[3], float filt)
 262 {
 263   double mul[3], max=0;
 264   int c;
 265
 266   FORC3 mul[c] = dq[c]/div[c];
 267   FORC3 if (max < mul[c]) max = mul[c];
 268   FORC3 curvep[c] = foveon_make_curve (max, mul[c], filt);
 269 }
 270
 271 int CLASS foveon_apply_curve (short *curve, int i)
 272 {
 273   if (abs(i) >= curve[0]) return 0;
 274   return i < 0 ? -curve[1-i] : curve[1+i];
 275 }
 276
 277 #ifdef image
 278 #undef image
 279 #endif
 280 #define image ((short(*)[4]) imgdata.image)
 281
 282 void CLASS foveon_interpolate()
 283 {
 284   static const short hood[] = { -1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1 };
 285   short *pix, prev[3], *curve[8], (*shrink)[3];
 286   float cfilt=0, ddft[3][3][2], ppm[3][3][3];
 287   float cam_xyz[3][3], correct[3][3], last[3][3], trans[3][3];
 288   float chroma_dq[3], color_dq[3], diag[3][3], div[3];
 289   float (*black)[3], (*sgain)[3], (*sgrow)[3];
 290   float fsum[3], val, frow, num;
 291   int row, col, c, i, j, diff, sgx, irow, sum, min, max, limit;
 292   int dscr[2][2], dstb[4], (*smrow[7])[3], total[4], ipix[3];
 293   int work[3][3], smlast, smred, smred_p=0, dev[3];
 294   int satlev[3], keep[4], active[4];
 295   unsigned dim[3], *badpix;
 296   double dsum=0, trsum[3];
 297   char str[128];
 298   const char* cp;
 299
 300
 301 #ifdef DCRAW_VERBOSE
 302   if (verbose)
 303     fprintf(stderr,_("Foveon interpolation...\n"));
 304 #endif
 305 #ifdef LIBRAW_LIBRARY_BUILD
 306   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,0,9);
 307 #endif
 308
 309   foveon_fixed (dscr, 4, "DarkShieldColRange");
 310   foveon_fixed (ppm[0][0], 27, "PostPolyMatrix");
 311   foveon_fixed (satlev, 3, "SaturationLevel");
 312   foveon_fixed (keep, 4, "KeepImageArea");
 313   foveon_fixed (active, 4, "ActiveImageArea");
 314   foveon_fixed (chroma_dq, 3, "ChromaDQ");
 315   foveon_fixed (color_dq, 3,
 316         foveon_camf_param ("IncludeBlocks", "ColorDQ") ?
 317                 "ColorDQ" : "ColorDQCamRGB");
 318   if (foveon_camf_param ("IncludeBlocks", "ColumnFilter"))
 319                  foveon_fixed (&cfilt, 1, "ColumnFilter");
 320
 321   memset (ddft, 0, sizeof ddft);
 322   if (!foveon_camf_param ("IncludeBlocks", "DarkDrift")
 323          || !foveon_fixed (ddft[1][0], 12, "DarkDrift"))
 324     for (i=0; i < 2; i++) {
 325       foveon_fixed (dstb, 4, i ? "DarkShieldBottom":"DarkShieldTop");
 326       for (row = dstb[1]; row <= dstb[3]; row++)
 327         for (col = dstb[0]; col <= dstb[2]; col++)
 328           FORC3 ddft[i+1][c][1] += (short) image[row*width+col][c];
 329       FORC3 ddft[i+1][c][1] /= (dstb[3]-dstb[1]+1) * (dstb[2]-dstb[0]+1);
 330     }
 331
 332   if (!(cp = foveon_camf_param ("WhiteBalanceIlluminants", model2)))
 333   {
 334 #ifdef DCRAW_VERBOSE
 335       fprintf (stderr,_("%s: Invalid white balance \"%s\"\n"), ifname, model2);
 336 #endif
 337 #ifdef LIBRAW_LIBRARY_BUILD
 338       imgdata.process_warnings |= LIBRAW_WARN_FOVEON_INVALIDWB;
 339 #endif
 340       return;
 341   }
 342   foveon_fixed (cam_xyz, 9, cp);
 343   foveon_fixed (correct, 9,
 344         foveon_camf_param ("WhiteBalanceCorrections", model2));
 345   memset (last, 0, sizeof last);
 346   for (i=0; i < 3; i++)
 347     for (j=0; j < 3; j++)
 348       FORC3 last[i][j] += correct[i][c] * cam_xyz[c][j];
 349
 350   #define LAST(x,y) last[(i+x)%3][(c+y)%3]
 351   for (i=0; i < 3; i++)
 352     FORC3 diag[c][i] = LAST(1,1)*LAST(2,2) - LAST(1,2)*LAST(2,1);
 353   #undef LAST
 354   FORC3 div[c] = diag[c][0]*0.3127 + diag[c][1]*0.329 + diag[c][2]*0.3583;
 355   sprintf (str, "%sRGBNeutral", model2);
 356   if (foveon_camf_param ("IncludeBlocks", str))
 357     foveon_fixed (div, 3, str);
 358   num = 0;
 359   FORC3 if (num < div[c]) num = div[c];
 360   FORC3 div[c] /= num;
 361
 362   memset (trans, 0, sizeof trans);
 363   for (i=0; i < 3; i++)
 364     for (j=0; j < 3; j++)
 365       FORC3 trans[i][j] += rgb_cam[i][c] * last[c][j] * div[j];
 366   FORC3 trsum[c] = trans[c][0] + trans[c][1] + trans[c][2];
 367   dsum = (6*trsum[0] + 11*trsum[1] + 3*trsum[2]) / 20;
 368   for (i=0; i < 3; i++)
 369     FORC3 last[i][c] = trans[i][c] * dsum / trsum[i];
 370   memset (trans, 0, sizeof trans);
 371   for (i=0; i < 3; i++)
 372     for (j=0; j < 3; j++)
 373       FORC3 trans[i][j] += (i==c ? 32 : -1) * last[c][j] / 30;
 374
 375   foveon_make_curves (curve, color_dq, div, cfilt);
 376   FORC3 chroma_dq[c] /= 3;
 377   foveon_make_curves (curve+3, chroma_dq, div, cfilt);
 378   FORC3 dsum += chroma_dq[c] / div[c];
 379   curve[6] = foveon_make_curve (dsum, dsum, cfilt);
 380   curve[7] = foveon_make_curve (dsum*2, dsum*2, cfilt);
 381
 382   sgain = (float (*)[3]) foveon_camf_matrix (dim, "SpatialGain");
 383   if (!sgain) return;
 384   sgrow = (float (*)[3]) calloc (dim[1], sizeof *sgrow);
 385   sgx = (width + dim[1]-2) / (dim[1]-1);
 386
 387   black = (float (*)[3]) calloc (height, sizeof *black);
 388 #ifdef LIBRAW_LIBRARY_BUILD
 389   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,1,9);
 390 #endif
 391   for (row=0; row < height; row++) {
 392     for (i=0; i < 6; i++)
 393       ddft[0][0][i] = ddft[1][0][i] +
 394         row / (height-1.0) * (ddft[2][0][i] - ddft[1][0][i]);
 395     FORC3 black[row][c] =
 396         ( foveon_avg (image[row*width]+c, dscr[0], cfilt) +
 397           foveon_avg (image[row*width]+c, dscr[1], cfilt) * 3
 398           - ddft[0][c][0] ) / 4 - ddft[0][c][1];
 399   }
 400   memcpy (black, black+8, sizeof *black*8);
 401   memcpy (black+height-11, black+height-22, 11*sizeof *black);
 402   memcpy (last, black, sizeof last);
 403
 404   for (row=1; row < height-1; row++) {
 405     FORC3 if (last[1][c] > last[0][c]) {
 406         if (last[1][c] > last[2][c])
 407           black[row][c] = (last[0][c] > last[2][c]) ? last[0][c]:last[2][c];
 408       } else
 409         if (last[1][c] < last[2][c])
 410           black[row][c] = (last[0][c] < last[2][c]) ? last[0][c]:last[2][c];
 411     memmove (last, last+1, 2*sizeof last[0]);
 412     memcpy (last[2], black[row+1], sizeof last[2]);
 413   }
 414   FORC3 black[row][c] = (last[0][c] + last[1][c])/2;
 415   FORC3 black[0][c] = (black[1][c] + black[3][c])/2;
 416
 417   val = 1 - exp(-1/24.0);
 418   memcpy (fsum, black, sizeof fsum);
 419   for (row=1; row < height; row++)
 420     FORC3 fsum[c] += black[row][c] =
 421         (black[row][c] - black[row-1][c])*val + black[row-1][c];
 422   memcpy (last[0], black[height-1], sizeof last[0]);
 423   FORC3 fsum[c] /= height;
 424   for (row = height; row--; )
 425     FORC3 last[0][c] = black[row][c] =
 426         (black[row][c] - fsum[c] - last[0][c])*val + last[0][c];
 427
 428   memset (total, 0, sizeof total);
 429   for (row=2; row < height; row+=4)
 430     for (col=2; col < width; col+=4) {
 431       FORC3 total[c] += (short) image[row*width+col][c];
 432       total[3]++;
 433     }
 434   for (row=0; row < height; row++)
 435     FORC3 black[row][c] += fsum[c]/2 + total[c]/(total[3]*100.0);
 436
 437 #ifdef LIBRAW_LIBRARY_BUILD
 438   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,2,9);
 439 #endif
 440   for (row=0; row < height; row++) {
 441     for (i=0; i < 6; i++)
 442       ddft[0][0][i] = ddft[1][0][i] +
 443         row / (height-1.0) * (ddft[2][0][i] - ddft[1][0][i]);
 444     pix = image[row*width];
 445     memcpy (prev, pix, sizeof prev);
 446     frow = row / (height-1.0) * (dim[2]-1);
 447     if ((irow = frow) == dim[2]-1) irow--;
 448     frow -= irow;
 449     for (i=0; i < dim[1]; i++)
 450       FORC3 sgrow[i][c] = sgain[ irow   *dim[1]+i][c] * (1-frow) +
 451                           sgain[(irow+1)*dim[1]+i][c] *    frow;
 452     for (col=0; col < width; col++) {
 453       FORC3 {
 454         diff = pix[c] - prev[c];
 455         prev[c] = pix[c];
 456         ipix[c] = pix[c] + floor ((diff + (diff*diff >> 14)) * cfilt
 457                 - ddft[0][c][1] - ddft[0][c][0] * ((float) col/width - 0.5)
 458                 - black[row][c] );
 459       }
 460       FORC3 {
 461         work[0][c] = ipix[c] * ipix[c] >> 14;
 462         work[2][c] = ipix[c] * work[0][c] >> 14;
 463         work[1][2-c] = ipix[(c+1) % 3] * ipix[(c+2) % 3] >> 14;
 464       }
 465       FORC3 {
 466         for (val=i=0; i < 3; i++)
 467           for (  j=0; j < 3; j++)
 468             val += ppm[c][i][j] * work[i][j];
 469         ipix[c] = floor ((ipix[c] + floor(val)) *
 470                 ( sgrow[col/sgx  ][c] * (sgx - col%sgx) +
 471                   sgrow[col/sgx+1][c] * (col%sgx) ) / sgx / div[c]);
 472         if (ipix[c] > 32000) ipix[c] = 32000;
 473         pix[c] = ipix[c];
 474       }
 475       pix += 4;
 476     }
 477   }
 478   free (black);
 479   free (sgrow);
 480   free (sgain);
 481
 482 #ifdef LIBRAW_LIBRARY_BUILD
 483   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,3,9);
 484 #endif
 485   if ((badpix = (unsigned int *) foveon_camf_matrix (dim, "BadPixels"))) {
 486     for (i=0; i < dim[0]; i++) {
 487       col = (badpix[i] >> 8 & 0xfff) - keep[0];
 488       row = (badpix[i] >> 20       ) - keep[1];
 489       if ((unsigned)(row-1) > height-3 || (unsigned)(col-1) > width-3)
 490         continue;
 491       memset (fsum, 0, sizeof fsum);
 492       for (sum=j=0; j < 8; j++)
 493         if (badpix[i] & (1 << j)) {
 494           FORC3 fsum[c] += (short)
 495                 image[(row+hood[j*2])*width+col+hood[j*2+1]][c];
 496           sum++;
 497         }
 498       if (sum) FORC3 image[row*width+col][c] = fsum[c]/sum;
 499     }
 500     free (badpix);
 501   }
 502
 503   /* Array for 5x5 Gaussian averaging of red values */
 504   smrow[6] = (int (*)[3]) calloc (width*5, sizeof **smrow);
 505   merror (smrow[6], "foveon_interpolate()");
 506   for (i=0; i < 5; i++)
 507     smrow[i] = smrow[6] + i*width;
 508
 509   /* Sharpen the reds against these Gaussian averages */
 510   for (smlast=-1, row=2; row < height-2; row++) {
 511     while (smlast < row+2) {
 512       for (i=0; i < 6; i++)
 513         smrow[(i+5) % 6] = smrow[i];
 514       pix = image[++smlast*width+2];
 515       for (col=2; col < width-2; col++) {
 516         smrow[4][col][0] =
 517           (pix[0]*6 + (pix[-4]+pix[4])*4 + pix[-8]+pix[8] + 8) >> 4;
 518         pix += 4;
 519       }
 520     }
 521     pix = image[row*width+2];
 522     for (col=2; col < width-2; col++) {
 523       smred = ( 6 *  smrow[2][col][0]
 524               + 4 * (smrow[1][col][0] + smrow[3][col][0])
 525               +      smrow[0][col][0] + smrow[4][col][0] + 8 ) >> 4;
 526       if (col == 2)
 527         smred_p = smred;
 528       i = pix[0] + ((pix[0] - ((smred*7 + smred_p) >> 3)) >> 3);
 529       if (i > 32000) i = 32000;
 530       pix[0] = i;
 531       smred_p = smred;
 532       pix += 4;
 533     }
 534   }
 535
 536 #ifdef LIBRAW_LIBRARY_BUILD
 537   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,4,9);
 538 #endif
 539   /* Adjust the brighter pixels for better linearity */
 540   min = 0xffff;
 541   FORC3 {
 542     i = satlev[c] / div[c];
 543     if (min > i) min = i;
 544   }
 545   limit = min * 9 >> 4;
 546   for (pix=image[0]; pix < image[height*width]; pix+=4) {
 547     if (pix[0] <= limit || pix[1] <= limit || pix[2] <= limit)
 548       continue;
 549     min = max = pix[0];
 550     for (c=1; c < 3; c++) {
 551       if (min > pix[c]) min = pix[c];
 552       if (max < pix[c]) max = pix[c];
 553     }
 554     if (min >= limit*2) {
 555       pix[0] = pix[1] = pix[2] = max;
 556     } else {
 557       i = 0x4000 - ((min - limit) << 14) / limit;
 558       i = 0x4000 - (i*i >> 14);
 559       i = i*i >> 14;
 560       FORC3 pix[c] += (max - pix[c]) * i >> 14;
 561     }
 562   }
 563 /*
 564    Because photons that miss one detector often hit another,
 565    the sum R+G+B is much less noisy than the individual colors.
 566    So smooth the hues without smoothing the total.
 567  */
 568 #ifdef LIBRAW_LIBRARY_BUILD
 569   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,5,9);
 570 #endif
 571   for (smlast=-1, row=2; row < height-2; row++) {
 572     while (smlast < row+2) {
 573       for (i=0; i < 6; i++)
 574         smrow[(i+5) % 6] = smrow[i];
 575       pix = image[++smlast*width+2];
 576       for (col=2; col < width-2; col++) {
 577         FORC3 smrow[4][col][c] = (pix[c-4]+2*pix[c]+pix[c+4]+2) >> 2;
 578         pix += 4;
 579       }
 580     }
 581     pix = image[row*width+2];
 582     for (col=2; col < width-2; col++) {
 583       FORC3 dev[c] = -foveon_apply_curve (curve[7], pix[c] -
 584         ((smrow[1][col][c] + 2*smrow[2][col][c] + smrow[3][col][c]) >> 2));
 585       sum = (dev[0] + dev[1] + dev[2]) >> 3;
 586       FORC3 pix[c] += dev[c] - sum;
 587       pix += 4;
 588     }
 589   }
 590   for (smlast=-1, row=2; row < height-2; row++) {
 591     while (smlast < row+2) {
 592       for (i=0; i < 6; i++)
 593         smrow[(i+5) % 6] = smrow[i];
 594       pix = image[++smlast*width+2];
 595       for (col=2; col < width-2; col++) {
 596         FORC3 smrow[4][col][c] =
 597                 (pix[c-8]+pix[c-4]+pix[c]+pix[c+4]+pix[c+8]+2) >> 2;
 598         pix += 4;
 599       }
 600     }
 601     pix = image[row*width+2];
 602     for (col=2; col < width-2; col++) {
 603       for (total[3]=375, sum=60, c=0; c < 3; c++) {
 604         for (total[c]=i=0; i < 5; i++)
 605           total[c] += smrow[i][col][c];
 606         total[3] += total[c];
 607         sum += pix[c];
 608       }
 609       if (sum < 0) sum = 0;
 610       j = total[3] > 375 ? (sum << 16) / total[3] : sum * 174;
 611       FORC3 pix[c] += foveon_apply_curve (curve[6],
 612                 ((j*total[c] + 0x8000) >> 16) - pix[c]);
 613       pix += 4;
 614     }
 615   }
 616
 617 #ifdef LIBRAW_LIBRARY_BUILD
 618   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,6,9);
 619 #endif
 620   /* Transform the image to a different colorspace */
 621   for (pix=image[0]; pix < image[height*width]; pix+=4) {
 622     FORC3 pix[c] -= foveon_apply_curve (curve[c], pix[c]);
 623     sum = (pix[0]+pix[1]+pix[1]+pix[2]) >> 2;
 624     FORC3 pix[c] -= foveon_apply_curve (curve[c], pix[c]-sum);
 625     FORC3 {
 626       for (dsum=i=0; i < 3; i++)
 627         dsum += trans[c][i] * pix[i];
 628       if (dsum < 0)  dsum = 0;
 629       if (dsum > 24000) dsum = 24000;
 630       ipix[c] = dsum + 0.5;
 631     }
 632     FORC3 pix[c] = ipix[c];
 633   }
 634
 635   /* Smooth the image bottom-to-top and save at 1/4 scale */
 636   shrink = (short (*)[3]) calloc ((width/4) * (height/4), sizeof *shrink);
 637   merror (shrink, "foveon_interpolate()");
 638   for (row = height/4; row--; )
 639     for (col=0; col < width/4; col++) {
 640       ipix[0] = ipix[1] = ipix[2] = 0;
 641       for (i=0; i < 4; i++)
 642         for (j=0; j < 4; j++)
 643           FORC3 ipix[c] += image[(row*4+i)*width+col*4+j][c];
 644       FORC3
 645         if (row+2 > height/4)
 646           shrink[row*(width/4)+col][c] = ipix[c] >> 4;
 647         else
 648           shrink[row*(width/4)+col][c] =
 649             (shrink[(row+1)*(width/4)+col][c]*1840 + ipix[c]*141 + 2048) >> 12;
 650     }
 651   /* From the 1/4-scale image, smooth right-to-left */
 652 #ifdef LIBRAW_LIBRARY_BUILD
 653   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,7,9);
 654 #endif
 655   for (row=0; row < (height & ~3); row++) {
 656     ipix[0] = ipix[1] = ipix[2] = 0;
 657     if ((row & 3) == 0)
 658       for (col = width & ~3 ; col--; )
 659         FORC3 smrow[0][col][c] = ipix[c] =
 660           (shrink[(row/4)*(width/4)+col/4][c]*1485 + ipix[c]*6707 + 4096) >> 13;
 661
 662   /* Then smooth left-to-right */
 663     ipix[0] = ipix[1] = ipix[2] = 0;
 664     for (col=0; col < (width & ~3); col++)
 665       FORC3 smrow[1][col][c] = ipix[c] =
 666         (smrow[0][col][c]*1485 + ipix[c]*6707 + 4096) >> 13;
 667
 668   /* Smooth top-to-bottom */
 669     if (row == 0)
 670       memcpy (smrow[2], smrow[1], sizeof **smrow * width);
 671     else
 672       for (col=0; col < (width & ~3); col++)
 673         FORC3 smrow[2][col][c] =
 674           (smrow[2][col][c]*6707 + smrow[1][col][c]*1485 + 4096) >> 13;
 675
 676   /* Adjust the chroma toward the smooth values */
 677     for (col=0; col < (width & ~3); col++) {
 678       for (i=j=30, c=0; c < 3; c++) {
 679         i += smrow[2][col][c];
 680         j += image[row*width+col][c];
 681       }
 682       j = (j << 16) / i;
 683       for (sum=c=0; c < 3; c++) {
 684         ipix[c] = foveon_apply_curve (curve[c+3],
 685           ((smrow[2][col][c] * j + 0x8000) >> 16) - image[row*width+col][c]);
 686         sum += ipix[c];
 687       }
 688       sum >>= 3;
 689       FORC3 {
 690         i = image[row*width+col][c] + ipix[c] - sum;
 691         if (i < 0) i = 0;
 692         image[row*width+col][c] = i;
 693       }
 694     }
 695   }
 696   free (shrink);
 697   free (smrow[6]);
 698   for (i=0; i < 8; i++)
 699     free (curve[i]);
 700 #ifdef LIBRAW_LIBRARY_BUILD
 701   RUN_CALLBACK(LIBRAW_PROGRESS_FOVEON_INTERPOLATE,8,9);
 702 #endif
 703
 704   /* Trim off the black border */
 705   active[1] -= keep[1];
 706   active[3] -= 2;
 707   i = active[2] - active[0];
 708   for (row=0; row < active[3]-active[1]; row++)
 709     memcpy (image[row*i], image[(row+active[1])*width+active[0]],
 710          i * sizeof *image);
 711   width = i;
 712   height = row;
 713 }
 714 #undef image
 715
 716 /* RESTRICTED code ends here */
 717 char * CLASS foveon_gets (int offset, char *str, int len)
 718 {
 719   int i;
 720   fseek (ifp, offset, SEEK_SET);
 721   for (i=0; i < len-1; i++)
 722     if ((str[i] = get2()) == 0) break;
 723   str[i] = 0;
 724   return str;
 725 }
 726
 727 void CLASS parse_foveon()
 728 {
 729   int entries, img=0, off, len, tag, save, i, wide, high, pent, poff[256][2];
 730   char name[64], value[64];
 731
 732   order = 0x4949;                       /* Little-endian */
 733   fseek (ifp, 36, SEEK_SET);
 734   flip = get4();
 735   fseek (ifp, -4, SEEK_END);
 736   fseek (ifp, get4(), SEEK_SET);
 737   if (get4() != 0x64434553) return;     /* SECd */
 738   entries = (get4(),get4());
 739   while (entries--) {
 740     off = get4();
 741     len = get4();
 742     tag = get4();
 743     save = ftell(ifp);
 744     fseek (ifp, off, SEEK_SET);
 745     if (get4() != (0x20434553 | (tag << 24))) return;
 746     switch (tag) {
 747       case 0x47414d49:                  /* IMAG */
 748       case 0x32414d49:                  /* IMA2 */
 749         fseek (ifp, 12, SEEK_CUR);
 750         wide = get4();
 751         high = get4();
 752         if (wide > raw_width && high > raw_height) {
 753           raw_width  = wide;
 754           raw_height = high;
 755           data_offset = off+24;
 756         }
 757         fseek (ifp, off+28, SEEK_SET);
 758         if (fgetc(ifp) == 0xff && fgetc(ifp) == 0xd8
 759                 && thumb_length < len-28) {
 760           thumb_offset = off+28;
 761           thumb_length = len-28;
 762           write_thumb = &CLASS jpeg_thumb;
 763         }
 764         if (++img == 2 && !thumb_length) {
 765           thumb_offset = off+24;
 766           thumb_width = wide;
 767           thumb_height = high;
 768           write_thumb = &CLASS foveon_thumb;
 769         }
 770         break;
 771       case 0x464d4143:                  /* CAMF */
 772         meta_offset = off+24;
 773         meta_length = len-28;
 774         if (meta_length > 0x20000)
 775             meta_length = 0x20000;
 776         break;
 777       case 0x504f5250:                  /* PROP */
 778         pent = (get4(),get4());
 779         fseek (ifp, 12, SEEK_CUR);
 780         off += pent*8 + 24;
 781         if ((unsigned) pent > 256) pent=256;
 782         for (i=0; i < pent*2; i++)
 783           poff[0][i] = off + get4()*2;
 784         for (i=0; i < pent; i++) {
 785           foveon_gets (poff[i][0], name, 64);
 786           foveon_gets (poff[i][1], value, 64);
 787           if (!strcmp (name, "ISO"))
 788             iso_speed = atoi(value);
 789           if (!strcmp (name, "CAMMANUF"))
 790             strcpy (make, value);
 791           if (!strcmp (name, "CAMMODEL"))
 792             strcpy (model, value);
 793           if (!strcmp (name, "WB_DESC"))
 794             strcpy (model2, value);
 795           if (!strcmp (name, "TIME"))
 796             timestamp = atoi(value);
 797           if (!strcmp (name, "EXPTIME"))
 798             shutter = atoi(value) / 1000000.0;
 799           if (!strcmp (name, "APERTURE"))
 800             aperture = atof(value);
 801           if (!strcmp (name, "FLENGTH"))
 802             focal_len = atof(value);
 803         }
 804 #ifdef LOCALTIME
 805         timestamp = mktime (gmtime (&timestamp));
 806 #endif
 807     }
 808     fseek (ifp, save, SEEK_SET);
 809   }
 810   is_foveon = 1;
 811 }