12e36160e3f0bccd55ae405ee0dfbfa35d225e6f
[goodguy/cinelerra.git] / cinelerra-5.1 / cinelerra / dcraw.C
1 /*
2    dcraw.c -- Dave Coffin's raw photo decoder
3    Copyright 1997-2018 by Dave Coffin, dcoffin a cybercom o net
4
5    This is a command-line ANSI C program to convert raw photos from
6    any digital camera on any computer running any operating system.
7
8    No license is required to download and use dcraw.c.  However,
9    to lawfully redistribute dcraw, you must either (a) offer, at
10    no extra charge, full source code* for all executable files
11    containing RESTRICTED functions, (b) distribute this code under
12    the GPL Version 2 or later, (c) remove all RESTRICTED functions,
13    re-implement them, or copy them from an earlier, unrestricted
14    Revision of dcraw.c, or (d) purchase a license from the author.
15
16    The functions that process Foveon images have been RESTRICTED
17    since Revision 1.237.  All other code remains free for all uses.
18
19    *If you have not modified dcraw.c in any way, a link to my
20    homepage qualifies as "full source code".
21
22    $Revision: 1.478 $
23    $Date: 2018/06/01 20:36:25 $
24  */
25
26 #define DCRAW_VERSION "9.28"
27
28 #ifndef _GNU_SOURCE
29 #define _GNU_SOURCE
30 #endif
31 #define _USE_MATH_DEFINES
32 #include <ctype.h>
33 #include <errno.h>
34 #include <fcntl.h>
35 #include <float.h>
36 #include <limits.h>
37 #include <math.h>
38 #include <setjmp.h>
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <time.h>
43 #include <sys/types.h>
44
45 #if defined(DJGPP) || defined(__MINGW32__)
46 #define fseeko fseek
47 #define ftello ftell
48 #else
49 #define fgetc getc_unlocked
50 #endif
51 #ifdef __CYGWIN__
52 #include <io.h>
53 #endif
54 #ifdef WIN32
55 #include <sys/utime.h>
56 #include <winsock2.h>
57 #pragma comment(lib, "ws2_32.lib")
58 #define snprintf _snprintf
59 #define strcasecmp stricmp
60 #define strncasecmp strnicmp
61 typedef __int64 INT64;
62 typedef unsigned __int64 UINT64;
63 #else
64 #include <unistd.h>
65 #include <utime.h>
66 #include <netinet/in.h>
67 #include <stdint.h>
68 typedef int64_t INT64;
69 typedef uint64_t UINT64;
70 #endif
71
72 #ifdef NODEPS
73 #define NO_JASPER
74 #define NO_JPEG
75 #define NO_LCMS
76 #endif
77 #ifndef NO_JASPER
78 #include <jasper/jasper.h>      /* Decode Red camera movies */
79 #endif
80 #ifndef NO_JPEG
81 #include <jpeglib.h>            /* Decode compressed Kodak DC120 photos */
82 #endif                          /* and Adobe Lossy DNGs */
83 #ifndef NO_LCMS
84 #include <lcms2.h>              /* Support color profiles */
85 #endif
86 #ifdef LOCALEDIR
87 #include <libintl.h>
88 #define _(String) gettext(String)
89 #else
90 #define _(String) (String)
91 #endif
92
93 #include "dcraw.h"
94
95 //const static data
96 const double DCRaw_data::xyz_rgb[3][3] = {              /* XYZ from RGB */
97         { 0.412453, 0.357580, 0.180423 },
98         { 0.212671, 0.715160, 0.072169 },
99         { 0.019334, 0.119193, 0.950227 } };
100 const float DCRaw_data::d65_white[3] = {
101           0.950456, 1.000000, 1.088754 };
102
103 /*
104    All global variables are defined here, and all functions that
105    access them are prefixed with "CLASS".  For thread-safety, all
106    non-const static local variables except cbrt[] must be declared
107    "thread_local".
108  */
109
110 #define FORC(cnt) for (c=0; c < cnt; c++)
111 #define FORC3 FORC(3)
112 #define FORC4 FORC(4)
113 #define FORCC FORC(colors)
114
115 #define SQR(x) ((x)*(x))
116 #define ABS(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31))
117 #define MIN(a,b) ((a) < (b) ? (a) : (b))
118 #define MAX(a,b) ((a) > (b) ? (a) : (b))
119 #define LIM(x,min,max) MAX(min,MIN(x,max))
120 #define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y))
121 #define CLIP(x) LIM((int)(x),0,65535)
122 #define SWAP(a,b) { a=a+b; b=a-b; a=a-b; }
123 #define ZERO(var) memset(&var, 0, sizeof var)
124 /*
125    In order to inline this calculation, I make the risky
126    assumption that all filter patterns can be described
127    by a repeating pattern of eight rows and two columns
128
129    Do not use the FC or BAYER macros with the Leaf CatchLight,
130    because its pattern is 16x16, not 2x8.
131
132    Return values are either 0/1/2/3 = G/M/C/Y or 0/1/2/3 = R/G1/B/G2
133
134         PowerShot 600   PowerShot A50   PowerShot Pro70 Pro90 & G1
135         0xe1e4e1e4:     0x1b4e4b1e:     0x1e4b4e1b:     0xb4b4b4b4:
136
137           0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5
138         0 G M G M G M   0 C Y C Y C Y   0 Y C Y C Y C   0 G M G M G M
139         1 C Y C Y C Y   1 M G M G M G   1 M G M G M G   1 Y C Y C Y C
140         2 M G M G M G   2 Y C Y C Y C   2 C Y C Y C Y
141         3 C Y C Y C Y   3 G M G M G M   3 G M G M G M
142                         4 C Y C Y C Y   4 Y C Y C Y C
143         PowerShot A5    5 G M G M G M   5 G M G M G M
144         0x1e4e1e4e:     6 Y C Y C Y C   6 C Y C Y C Y
145                         7 M G M G M G   7 M G M G M G
146           0 1 2 3 4 5
147         0 C Y C Y C Y
148         1 G M G M G M
149         2 C Y C Y C Y
150         3 M G M G M G
151
152    All RGB cameras use one of these Bayer grids:
153
154         0x16161616:     0x61616161:     0x49494949:     0x94949494:
155
156           0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5
157         0 B G B G B G   0 G R G R G R   0 G B G B G B   0 R G R G R G
158         1 G R G R G R   1 B G B G B G   1 R G R G R G   1 G B G B G B
159         2 B G B G B G   2 G R G R G R   2 G B G B G B   2 R G R G R G
160         3 G R G R G R   3 B G B G B G   3 R G R G R G   3 G B G B G B
161  */
162
163 #define RAW(row,col) \
164         raw_image[(row)*raw_width+(col)]
165
166 #define FC(row,col) \
167         (filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3)
168
169 #define BAYER(row,col) \
170         image[((row) >> shrink)*iwidth + ((col) >> shrink)][FC(row,col)]
171
172 #define BAYER2(row,col) \
173         image[((row) >> shrink)*iwidth + ((col) >> shrink)][fcol(row,col)]
174
175 int CLASS fcol (int row, int col)
176 {
177   static const char filter[16][16] =
178   { { 2,1,1,3,2,3,2,0,3,2,3,0,1,2,1,0 },
179     { 0,3,0,2,0,1,3,1,0,1,1,2,0,3,3,2 },
180     { 2,3,3,2,3,1,1,3,3,1,2,1,2,0,0,3 },
181     { 0,1,0,1,0,2,0,2,2,0,3,0,1,3,2,1 },
182     { 3,1,1,2,0,1,0,2,1,3,1,3,0,1,3,0 },
183     { 2,0,0,3,3,2,3,1,2,0,2,0,3,2,2,1 },
184     { 2,3,3,1,2,1,2,1,2,1,1,2,3,0,0,1 },
185     { 1,0,0,2,3,0,0,3,0,3,0,3,2,1,2,3 },
186     { 2,3,3,1,1,2,1,0,3,2,3,0,2,3,1,3 },
187     { 1,0,2,0,3,0,3,2,0,1,1,2,0,1,0,2 },
188     { 0,1,1,3,3,2,2,1,1,3,3,0,2,1,3,2 },
189     { 2,3,2,0,0,1,3,0,2,0,1,2,3,0,1,0 },
190     { 1,3,1,2,3,2,3,2,0,2,0,1,1,0,3,0 },
191     { 0,2,0,3,1,0,0,1,1,3,3,2,3,2,2,1 },
192     { 2,1,3,2,3,1,2,1,0,3,0,2,0,2,0,2 },
193     { 0,3,1,0,0,2,0,3,2,1,3,1,1,3,1,3 } };
194
195   if (filters == 1) return filter[(row+top_margin)&15][(col+left_margin)&15];
196   if (filters == 9) return xtrans[(row+6) % 6][(col+6) % 6];
197   return FC(row,col);
198 }
199
200 void CLASS reset()
201 {
202 // zero data segment
203         DCRaw_data *data = (DCRaw_data *)this;
204         memset(data, 0, sizeof(*data));
205 // non-zero init data
206         aber[0] = aber[1] = aber[2] = aber[3] = 1;
207         gamm[0] = 0.45; gamm[1] = 4.5;
208         bright = 1;
209         use_camera_matrix = 1;
210         output_color = 1;
211         output_bps = 8;
212         greybox[2] = UINT_MAX;  greybox[3] = UINT_MAX;
213 }
214
215 #if 0
216 char *my_memmem (char *haystack, size_t haystacklen,
217               char *needle, size_t needlelen)
218 {
219   char *c;
220   for (c = haystack; c <= haystack + haystacklen - needlelen; c++)
221     if (!memcmp (c, needle, needlelen))
222       return c;
223   return 0;
224 }
225 #define memmem my_memmem
226 char *my_strcasestr (char *haystack, const char *needle)
227 {
228   char *c;
229   for (c = haystack; *c; c++)
230     if (!strncasecmp(c, needle, strlen(needle)))
231       return c;
232   return 0;
233 }
234 #define strcasestr my_strcasestr
235 #endif
236
237 void CLASS merror (void *ptr, const char *where)
238 {
239   if (ptr) return;
240   fprintf (stderr,_("%s: Out of memory in %s\n"), ifname, where);
241   longjmp (failure, 1);
242 }
243
244 void CLASS derror()
245 {
246   if (!data_error) {
247     fprintf (stderr, "%s: ", ifname);
248     if (feof(ifp))
249       fprintf (stderr,_("Unexpected end of file\n"));
250     else
251       fprintf (stderr,_("Corrupt data near 0x%jx\n"), (INT64) ftello(ifp));
252   }
253   data_error++;
254 }
255
256 ushort CLASS sget2 (uchar *s)
257 {
258   if (order == 0x4949)          /* "II" means little-endian */
259     return s[0] | s[1] << 8;
260   else                          /* "MM" means big-endian */
261     return s[0] << 8 | s[1];
262 }
263
264 ushort CLASS get2()
265 {
266   uchar str[2] = { 0xff,0xff };
267   fread (str, 1, 2, ifp);
268   return sget2(str);
269 }
270
271 unsigned CLASS sget4 (uchar *s)
272 {
273   if (order == 0x4949)
274     return s[0] | s[1] << 8 | s[2] << 16 | s[3] << 24;
275   else
276     return s[0] << 24 | s[1] << 16 | s[2] << 8 | s[3];
277 }
278 #define sget4(s) sget4((uchar *)s)
279
280 unsigned CLASS get4()
281 {
282   uchar str[4] = { 0xff,0xff,0xff,0xff };
283   fread (str, 1, 4, ifp);
284   return sget4(str);
285 }
286
287 unsigned CLASS getint (int type)
288 {
289   return type == 3 ? get2() : get4();
290 }
291
292 float CLASS int_to_float (int i)
293 {
294   union { int i; float f; } u;
295   u.i = i;
296   return u.f;
297 }
298
299 double CLASS getreal (int type)
300 {
301   union { char c[8]; double d; } u;
302   int i, rev;
303
304   switch (type) {
305     case 3: return (unsigned short) get2();
306     case 4: return (unsigned int) get4();
307     case 5:  u.d = (unsigned int) get4();
308       return u.d / (unsigned int) get4();
309     case 8: return (signed short) get2();
310     case 9: return (signed int) get4();
311     case 10: u.d = (signed int) get4();
312       return u.d / (signed int) get4();
313     case 11: return int_to_float (get4());
314     case 12:
315       rev = 7 * ((order == 0x4949) == (ntohs(0x1234) == 0x1234));
316       for (i=0; i < 8; i++)
317         u.c[i ^ rev] = fgetc(ifp);
318       return u.d;
319     default: return fgetc(ifp);
320   }
321 }
322
323 static void swabb(const void *bfrom, void *bto, ssize_t n)
324 {
325   const char *from = (const char *) bfrom;
326   char *to = (char *) bto;
327   n &= ~((ssize_t) 1);
328   while (n > 1) {
329     const char b0 = from[--n], b1 = from[--n];
330     to[n] = b0;  to[n + 1] = b1;
331   }
332 }
333
334 void CLASS read_shorts (ushort *pixel, int count)
335 {
336   if (fread (pixel, 2, count, ifp) < count) derror();
337   if ((order == 0x4949) == (ntohs(0x1234) == 0x1234))
338     swabb(pixel, pixel, count*2);
339 }
340
341 void CLASS cubic_spline (const int *x_, const int *y_, const int len)
342 {
343   float **A, *b, *c, *d, *x, *y;
344   int i, j;
345
346   A = (float **) calloc (((2*len + 4)*sizeof **A + sizeof *A), 2*len);
347   if (!A) return;
348   A[0] = (float *) (A + 2*len);
349   for (i = 1; i < 2*len; i++)
350     A[i] = A[0] + 2*len*i;
351   y = len + (x = i + (d = i + (c = i + (b = A[0] + i*i))));
352   for (i = 0; i < len; i++) {
353     x[i] = x_[i] / 65535.0;
354     y[i] = y_[i] / 65535.0;
355   }
356   for (i = len-1; i > 0; i--) {
357     b[i] = (y[i] - y[i-1]) / (x[i] - x[i-1]);
358     d[i-1] = x[i] - x[i-1];
359   }
360   for (i = 1; i < len-1; i++) {
361     A[i][i] = 2 * (d[i-1] + d[i]);
362     if (i > 1) {
363       A[i][i-1] = d[i-1];
364       A[i-1][i] = d[i-1];
365     }
366     A[i][len-1] = 6 * (b[i+1] - b[i]);
367   }
368   for(i = 1; i < len-2; i++) {
369     float v = A[i+1][i] / A[i][i];
370     for(j = 1; j <= len-1; j++)
371       A[i+1][j] -= v * A[i][j];
372   }
373   for(i = len-2; i > 0; i--) {
374     float acc = 0;
375     for(j = i; j <= len-2; j++)
376       acc += A[i][j]*c[j];
377     c[i] = (A[i][len-1] - acc) / A[i][i];
378   }
379   for (i = 0; i < 0x10000; i++) {
380     float x_out = (float)(i / 65535.0);
381     float y_out = 0;
382     for (j = 0; j < len-1; j++) {
383       if (x[j] <= x_out && x_out <= x[j+1]) {
384         float v = x_out - x[j];
385         y_out = y[j] +
386           ((y[j+1] - y[j]) / d[j] - (2 * d[j] * c[j] + c[j+1] * d[j])/6) * v
387            + (c[j] * 0.5) * v*v + ((c[j+1] - c[j]) / (6 * d[j])) * v*v*v;
388       }
389     }
390     curve[i] = y_out < 0.0 ? 0 : (y_out >= 1.0 ? 65535 :
391                 (ushort)(y_out * 65535.0 + 0.5));
392   }
393   free (A);
394 }
395
396 void CLASS canon_600_fixed_wb (int temp)
397 {
398   static const short mul[4][5] = {
399     {  667, 358,397,565,452 },
400     {  731, 390,367,499,517 },
401     { 1119, 396,348,448,537 },
402     { 1399, 485,431,508,688 } };
403   int lo, hi, i;
404   float frac=0;
405
406   for (lo=4; --lo; )
407     if (*mul[lo] <= temp) break;
408   for (hi=0; hi < 3; hi++)
409     if (*mul[hi] >= temp) break;
410   if (lo != hi)
411     frac = (float) (temp - *mul[lo]) / (*mul[hi] - *mul[lo]);
412   for (i=1; i < 5; i++)
413     pre_mul[i-1] = 1 / (frac * mul[hi][i] + (1-frac) * mul[lo][i]);
414 }
415
416 /* Return values:  0 = white  1 = near white  2 = not white */
417 int CLASS canon_600_color (int ratio[2], int mar)
418 {
419   int clipped=0, target, miss;
420
421   if (flash_used) {
422     if (ratio[1] < -104)
423       { ratio[1] = -104; clipped = 1; }
424     if (ratio[1] >   12)
425       { ratio[1] =   12; clipped = 1; }
426   } else {
427     if (ratio[1] < -264 || ratio[1] > 461) return 2;
428     if (ratio[1] < -50)
429       { ratio[1] = -50; clipped = 1; }
430     if (ratio[1] > 307)
431       { ratio[1] = 307; clipped = 1; }
432   }
433   target = flash_used || ratio[1] < 197
434         ? -38 - (398 * ratio[1] >> 10)
435         : -123 + (48 * ratio[1] >> 10);
436   if (target - mar <= ratio[0] &&
437       target + 20  >= ratio[0] && !clipped) return 0;
438   miss = target - ratio[0];
439   if (abs(miss) >= mar*4) return 2;
440   if (miss < -20) miss = -20;
441   if (miss > mar) miss = mar;
442   ratio[0] = target - miss;
443   return 1;
444 }
445
446 void CLASS canon_600_auto_wb()
447 {
448   int mar, row, col, i, j, st, count[] = { 0,0 };
449   int test[8], total[2][8], ratio[2][2], stat[2];
450
451   memset (&total, 0, sizeof total);
452   i = canon_ev + 0.5;
453   if      (i < 10) mar = 150;
454   else if (i > 12) mar = 20;
455   else mar = 280 - 20 * i;
456   if (flash_used) mar = 80;
457   for (row=14; row < height-14; row+=4)
458     for (col=10; col < width; col+=2) {
459       for (i=0; i < 8; i++)
460         test[(i & 4) + FC(row+(i >> 1),col+(i & 1))] =
461                     BAYER(row+(i >> 1),col+(i & 1));
462       for (i=0; i < 8; i++)
463         if (test[i] < 150 || test[i] > 1500) goto next;
464       for (i=0; i < 4; i++)
465         if (abs(test[i] - test[i+4]) > 50) goto next;
466       for (i=0; i < 2; i++) {
467         for (j=0; j < 4; j+=2)
468           ratio[i][j >> 1] = ((test[i*4+j+1]-test[i*4+j]) << 10) / test[i*4+j];
469         stat[i] = canon_600_color (ratio[i], mar);
470       }
471       if ((st = stat[0] | stat[1]) > 1) goto next;
472       for (i=0; i < 2; i++)
473         if (stat[i])
474           for (j=0; j < 2; j++)
475             test[i*4+j*2+1] = test[i*4+j*2] * (0x400 + ratio[i][j]) >> 10;
476       for (i=0; i < 8; i++)
477         total[st][i] += test[i];
478       count[st]++;
479 next: ;
480     }
481   if (count[0] | count[1]) {
482     st = count[0]*200 < count[1];
483     for (i=0; i < 4; i++)
484       pre_mul[i] = 1.0 / (total[st][i] + total[st][i+4]);
485   }
486 }
487
488 void CLASS canon_600_coeff()
489 {
490   static const short table[6][12] = {
491     { -190,702,-1878,2390,   1861,-1349,905,-393, -432,944,2617,-2105  },
492     { -1203,1715,-1136,1648, 1388,-876,267,245,  -1641,2153,3921,-3409 },
493     { -615,1127,-1563,2075,  1437,-925,509,3,     -756,1268,2519,-2007 },
494     { -190,702,-1886,2398,   2153,-1641,763,-251, -452,964,3040,-2528  },
495     { -190,702,-1878,2390,   1861,-1349,905,-393, -432,944,2617,-2105  },
496     { -807,1319,-1785,2297,  1388,-876,769,-257,  -230,742,2067,-1555  } };
497   int t=0, i, c;
498   float mc, yc;
499
500   mc = pre_mul[1] / pre_mul[2];
501   yc = pre_mul[3] / pre_mul[2];
502   if (mc > 1 && mc <= 1.28 && yc < 0.8789) t=1;
503   if (mc > 1.28 && mc <= 2) {
504     if  (yc < 0.8789) t=3;
505     else if (yc <= 2) t=4;
506   }
507   if (flash_used) t=5;
508   for (raw_color = i=0; i < 3; i++)
509     FORCC rgb_cam[i][c] = table[t][i*4 + c] / 1024.0;
510 }
511
512 void CLASS canon_600_load_raw()
513 {
514   uchar  data[1120], *dp;
515   ushort *pix;
516   int irow, row;
517
518   for (irow=row=0; irow < height; irow++) {
519     if (fread (data, 1, 1120, ifp) < 1120) derror();
520     pix = raw_image + row*raw_width;
521     for (dp=data; dp < data+1120;  dp+=10, pix+=8) {
522       pix[0] = (dp[0] << 2) + (dp[1] >> 6    );
523       pix[1] = (dp[2] << 2) + (dp[1] >> 4 & 3);
524       pix[2] = (dp[3] << 2) + (dp[1] >> 2 & 3);
525       pix[3] = (dp[4] << 2) + (dp[1]      & 3);
526       pix[4] = (dp[5] << 2) + (dp[9]      & 3);
527       pix[5] = (dp[6] << 2) + (dp[9] >> 2 & 3);
528       pix[6] = (dp[7] << 2) + (dp[9] >> 4 & 3);
529       pix[7] = (dp[8] << 2) + (dp[9] >> 6    );
530     }
531     if ((row+=2) > height) row = 1;
532   }
533 }
534
535 void CLASS canon_600_correct()
536 {
537   int row, col, val;
538   static const short mul[4][2] =
539   { { 1141,1145 }, { 1128,1109 }, { 1178,1149 }, { 1128,1109 } };
540
541   for (row=0; row < height; row++)
542     for (col=0; col < width; col++) {
543       if ((val = BAYER(row,col) - black) < 0) val = 0;
544       val = val * mul[row & 3][col & 1] >> 9;
545       BAYER(row,col) = val;
546     }
547   canon_600_fixed_wb(1311);
548   canon_600_auto_wb();
549   canon_600_coeff();
550   maximum = (0x3ff - black) * 1109 >> 9;
551   black = 0;
552 }
553
554 int CLASS canon_s2is()
555 {
556   unsigned row;
557
558   for (row=0; row < 100; row++) {
559     fseek (ifp, row*3340 + 3284, SEEK_SET);
560     if (getc(ifp) > 15) return 1;
561   }
562   return 0;
563 }
564
565 unsigned CLASS getbithuff (int nbits, ushort *huff)
566 {
567   unsigned c;
568   unsigned bitbuf;
569   int vbits, reset;
570
571   if (nbits > 25) return 0;
572   if (nbits < 0)
573     return gbh_bitbuf = gbh_vbits = gbh_reset = 0;
574   if (nbits == 0 || gbh_vbits < 0) return 0;
575   bitbuf = gbh_bitbuf;  vbits = gbh_vbits;  reset = gbh_reset;
576   while (!reset && vbits < nbits && (c = fgetc(ifp)) != EOF &&
577     !(reset = zero_after_ff && c == 0xff && fgetc(ifp))) {
578     bitbuf = (bitbuf << 8) + (uchar) c;
579     vbits += 8;
580   }
581   c = bitbuf << (32-vbits) >> (32-nbits);
582   if (huff) {
583     vbits -= huff[c] >> 8;
584     c = (uchar) huff[c];
585   } else
586     vbits -= nbits;
587   if (vbits < 0) derror();
588   gbh_bitbuf = bitbuf;  gbh_vbits = vbits;  gbh_reset = reset;
589   return c;
590 }
591
592 #define getbits(n) getbithuff(n,0)
593 #define gethuff(h) getbithuff(*h,h+1)
594
595 /*
596    Construct a decode tree according the specification in *source.
597    The first 16 bytes specify how many codes should be 1-bit, 2-bit
598    3-bit, etc.  Bytes after that are the leaf values.
599
600    For example, if the source is
601
602     { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,
603       0x04,0x03,0x05,0x06,0x02,0x07,0x01,0x08,0x09,0x00,0x0a,0x0b,0xff  },
604
605    then the code is
606
607         00              0x04
608         010             0x03
609         011             0x05
610         100             0x06
611         101             0x02
612         1100            0x07
613         1101            0x01
614         11100           0x08
615         11101           0x09
616         11110           0x00
617         111110          0x0a
618         1111110         0x0b
619         1111111         0xff
620  */
621 ushort * CLASS make_decoder_ref (const uchar **source)
622 {
623   int max, len, h, i, j;
624   const uchar *count;
625   ushort *huff;
626
627   count = (*source += 16) - 17;
628   for (max=16; max && !count[max]; max--);
629   huff = (ushort *) calloc (1 + (1 << max), sizeof *huff);
630   merror (huff, "make_decoder()");
631   huff[0] = max;
632   for (h=len=1; len <= max; len++)
633     for (i=0; i < count[len]; i++, ++*source)
634       for (j=0; j < 1 << (max-len); j++)
635         if (h <= 1 << max)
636           huff[h++] = len << 8 | **source;
637   return huff;
638 }
639
640 ushort * CLASS make_decoder (const uchar *source)
641 {
642   return make_decoder_ref (&source);
643 }
644
645 void CLASS crw_init_tables (unsigned table, ushort *huff[2])
646 {
647   static const uchar first_tree[3][29] = {
648     { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,
649       0x04,0x03,0x05,0x06,0x02,0x07,0x01,0x08,0x09,0x00,0x0a,0x0b,0xff  },
650     { 0,2,2,3,1,1,1,1,2,0,0,0,0,0,0,0,
651       0x03,0x02,0x04,0x01,0x05,0x00,0x06,0x07,0x09,0x08,0x0a,0x0b,0xff  },
652     { 0,0,6,3,1,1,2,0,0,0,0,0,0,0,0,0,
653       0x06,0x05,0x07,0x04,0x08,0x03,0x09,0x02,0x00,0x0a,0x01,0x0b,0xff  },
654   };
655   static const uchar second_tree[3][180] = {
656     { 0,2,2,2,1,4,2,1,2,5,1,1,0,0,0,139,
657       0x03,0x04,0x02,0x05,0x01,0x06,0x07,0x08,
658       0x12,0x13,0x11,0x14,0x09,0x15,0x22,0x00,0x21,0x16,0x0a,0xf0,
659       0x23,0x17,0x24,0x31,0x32,0x18,0x19,0x33,0x25,0x41,0x34,0x42,
660       0x35,0x51,0x36,0x37,0x38,0x29,0x79,0x26,0x1a,0x39,0x56,0x57,
661       0x28,0x27,0x52,0x55,0x58,0x43,0x76,0x59,0x77,0x54,0x61,0xf9,
662       0x71,0x78,0x75,0x96,0x97,0x49,0xb7,0x53,0xd7,0x74,0xb6,0x98,
663       0x47,0x48,0x95,0x69,0x99,0x91,0xfa,0xb8,0x68,0xb5,0xb9,0xd6,
664       0xf7,0xd8,0x67,0x46,0x45,0x94,0x89,0xf8,0x81,0xd5,0xf6,0xb4,
665       0x88,0xb1,0x2a,0x44,0x72,0xd9,0x87,0x66,0xd4,0xf5,0x3a,0xa7,
666       0x73,0xa9,0xa8,0x86,0x62,0xc7,0x65,0xc8,0xc9,0xa1,0xf4,0xd1,
667       0xe9,0x5a,0x92,0x85,0xa6,0xe7,0x93,0xe8,0xc1,0xc6,0x7a,0x64,
668       0xe1,0x4a,0x6a,0xe6,0xb3,0xf1,0xd3,0xa5,0x8a,0xb2,0x9a,0xba,
669       0x84,0xa4,0x63,0xe5,0xc5,0xf3,0xd2,0xc4,0x82,0xaa,0xda,0xe4,
670       0xf2,0xca,0x83,0xa3,0xa2,0xc3,0xea,0xc2,0xe2,0xe3,0xff,0xff  },
671     { 0,2,2,1,4,1,4,1,3,3,1,0,0,0,0,140,
672       0x02,0x03,0x01,0x04,0x05,0x12,0x11,0x06,
673       0x13,0x07,0x08,0x14,0x22,0x09,0x21,0x00,0x23,0x15,0x31,0x32,
674       0x0a,0x16,0xf0,0x24,0x33,0x41,0x42,0x19,0x17,0x25,0x18,0x51,
675       0x34,0x43,0x52,0x29,0x35,0x61,0x39,0x71,0x62,0x36,0x53,0x26,
676       0x38,0x1a,0x37,0x81,0x27,0x91,0x79,0x55,0x45,0x28,0x72,0x59,
677       0xa1,0xb1,0x44,0x69,0x54,0x58,0xd1,0xfa,0x57,0xe1,0xf1,0xb9,
678       0x49,0x47,0x63,0x6a,0xf9,0x56,0x46,0xa8,0x2a,0x4a,0x78,0x99,
679       0x3a,0x75,0x74,0x86,0x65,0xc1,0x76,0xb6,0x96,0xd6,0x89,0x85,
680       0xc9,0xf5,0x95,0xb4,0xc7,0xf7,0x8a,0x97,0xb8,0x73,0xb7,0xd8,
681       0xd9,0x87,0xa7,0x7a,0x48,0x82,0x84,0xea,0xf4,0xa6,0xc5,0x5a,
682       0x94,0xa4,0xc6,0x92,0xc3,0x68,0xb5,0xc8,0xe4,0xe5,0xe6,0xe9,
683       0xa2,0xa3,0xe3,0xc2,0x66,0x67,0x93,0xaa,0xd4,0xd5,0xe7,0xf8,
684       0x88,0x9a,0xd7,0x77,0xc4,0x64,0xe2,0x98,0xa5,0xca,0xda,0xe8,
685       0xf3,0xf6,0xa9,0xb2,0xb3,0xf2,0xd2,0x83,0xba,0xd3,0xff,0xff  },
686     { 0,0,6,2,1,3,3,2,5,1,2,2,8,10,0,117,
687       0x04,0x05,0x03,0x06,0x02,0x07,0x01,0x08,
688       0x09,0x12,0x13,0x14,0x11,0x15,0x0a,0x16,0x17,0xf0,0x00,0x22,
689       0x21,0x18,0x23,0x19,0x24,0x32,0x31,0x25,0x33,0x38,0x37,0x34,
690       0x35,0x36,0x39,0x79,0x57,0x58,0x59,0x28,0x56,0x78,0x27,0x41,
691       0x29,0x77,0x26,0x42,0x76,0x99,0x1a,0x55,0x98,0x97,0xf9,0x48,
692       0x54,0x96,0x89,0x47,0xb7,0x49,0xfa,0x75,0x68,0xb6,0x67,0x69,
693       0xb9,0xb8,0xd8,0x52,0xd7,0x88,0xb5,0x74,0x51,0x46,0xd9,0xf8,
694       0x3a,0xd6,0x87,0x45,0x7a,0x95,0xd5,0xf6,0x86,0xb4,0xa9,0x94,
695       0x53,0x2a,0xa8,0x43,0xf5,0xf7,0xd4,0x66,0xa7,0x5a,0x44,0x8a,
696       0xc9,0xe8,0xc8,0xe7,0x9a,0x6a,0x73,0x4a,0x61,0xc7,0xf4,0xc6,
697       0x65,0xe9,0x72,0xe6,0x71,0x91,0x93,0xa6,0xda,0x92,0x85,0x62,
698       0xf3,0xc5,0xb2,0xa4,0x84,0xba,0x64,0xa5,0xb3,0xd2,0x81,0xe5,
699       0xd3,0xaa,0xc4,0xca,0xf2,0xb1,0xe4,0xd1,0x83,0x63,0xea,0xc3,
700       0xe2,0x82,0xf1,0xa3,0xc2,0xa1,0xc1,0xe3,0xa2,0xe1,0xff,0xff  }
701   };
702   if (table > 2) table = 2;
703   huff[0] = make_decoder ( first_tree[table]);
704   huff[1] = make_decoder (second_tree[table]);
705 }
706
707 /*
708    Return 0 if the image starts with compressed data,
709    1 if it starts with uncompressed low-order bits.
710
711    In Canon compressed data, 0xff is always followed by 0x00.
712  */
713 int CLASS canon_has_lowbits()
714 {
715   uchar test[0x4000];
716   int ret=1, i;
717
718   fseek (ifp, 0, SEEK_SET);
719   fread (test, 1, sizeof test, ifp);
720   for (i=540; i < sizeof test - 1; i++)
721     if (test[i] == 0xff) {
722       if (test[i+1]) return 1;
723       ret=0;
724     }
725   return ret;
726 }
727
728 void CLASS canon_load_raw()
729 {
730   ushort *pixel, *prow, *huff[2];
731   int nblocks, lowbits, i, c, row, r, save, val;
732   int block, diffbuf[64], leaf, len, diff, carry=0, pnum=0, base[2];
733
734   crw_init_tables (tiff_compress, huff);
735   lowbits = canon_has_lowbits();
736   if (!lowbits) maximum = 0x3ff;
737   fseek (ifp, 540 + lowbits*raw_height*raw_width/4, SEEK_SET);
738   zero_after_ff = 1;
739   getbits(-1);
740   for (row=0; row < raw_height; row+=8) {
741     pixel = raw_image + row*raw_width;
742     nblocks = MIN (8, raw_height-row) * raw_width >> 6;
743     for (block=0; block < nblocks; block++) {
744       memset (diffbuf, 0, sizeof diffbuf);
745       for (i=0; i < 64; i++ ) {
746         leaf = gethuff(huff[i > 0]);
747         if (leaf == 0 && i) break;
748         if (leaf == 0xff) continue;
749         i  += leaf >> 4;
750         len = leaf & 15;
751         if (len == 0) continue;
752         diff = getbits(len);
753         if ((diff & (1 << (len-1))) == 0)
754           diff -= (1 << len) - 1;
755         if (i < 64) diffbuf[i] = diff;
756       }
757       diffbuf[0] += carry;
758       carry = diffbuf[0];
759       for (i=0; i < 64; i++ ) {
760         if (pnum++ % raw_width == 0)
761           base[0] = base[1] = 512;
762         if ((pixel[(block << 6) + i] = base[i & 1] += diffbuf[i]) >> 10)
763           derror();
764       }
765     }
766     if (lowbits) {
767       save = ftell(ifp);
768       fseek (ifp, 26 + row*raw_width/4, SEEK_SET);
769       for (prow=pixel, i=0; i < raw_width*2; i++) {
770         c = fgetc(ifp);
771         for (r=0; r < 8; r+=2, prow++) {
772           val = (*prow << 2) + ((c >> r) & 3);
773           if (raw_width == 2672 && val < 512) val += 2;
774           *prow = val;
775         }
776       }
777       fseek (ifp, save, SEEK_SET);
778     }
779   }
780   FORC(2) free (huff[c]);
781 }
782
783 struct jhead {
784   int algo, bits, high, wide, clrs, sraw, psv, restart, vpred[6];
785   ushort quant[64], idct[64], *huff[20], *free[20], *row;
786 };
787
788 int CLASS ljpeg_start (struct jhead *jh, int info_only)
789 {
790   ushort c, tag, len;
791   uchar data[0x10000];
792   const uchar *dp;
793
794   memset (jh, 0, sizeof *jh);
795   jh->restart = INT_MAX;
796   if ((fgetc(ifp),fgetc(ifp)) != 0xd8) return 0;
797   do {
798     if (!fread (data, 2, 2, ifp)) return 0;
799     tag =  data[0] << 8 | data[1];
800     len = (data[2] << 8 | data[3]) - 2;
801     if (tag <= 0xff00) return 0;
802     fread (data, 1, len, ifp);
803     switch (tag) {
804       case 0xffc3:
805         jh->sraw = ((data[7] >> 4) * (data[7] & 15) - 1) & 3;
806       case 0xffc1:
807       case 0xffc0:
808         jh->algo = tag & 0xff;
809         jh->bits = data[0];
810         jh->high = data[1] << 8 | data[2];
811         jh->wide = data[3] << 8 | data[4];
812         jh->clrs = data[5] + jh->sraw;
813         if (len == 9 && !dng_version) getc(ifp);
814         break;
815       case 0xffc4:
816         if (info_only) break;
817         for (dp = data; dp < data+len && !((c = *dp++) & -20); )
818           jh->free[c] = jh->huff[c] = make_decoder_ref (&dp);
819         break;
820       case 0xffda:
821         jh->psv = data[1+data[0]*2];
822         jh->bits -= data[3+data[0]*2] & 15;
823         break;
824       case 0xffdb:
825         FORC(64) jh->quant[c] = data[c*2+1] << 8 | data[c*2+2];
826         break;
827       case 0xffdd:
828         jh->restart = data[0] << 8 | data[1];
829     }
830   } while (tag != 0xffda);
831   if (jh->bits > 16 || jh->clrs > 6 ||
832      !jh->bits || !jh->high || !jh->wide || !jh->clrs) return 0;
833   if (info_only) return 1;
834   if (!jh->huff[0]) return 0;
835   FORC(19) if (!jh->huff[c+1]) jh->huff[c+1] = jh->huff[c];
836   if (jh->sraw) {
837     FORC(4)        jh->huff[2+c] = jh->huff[1];
838     FORC(jh->sraw) jh->huff[1+c] = jh->huff[0];
839   }
840   jh->row = (ushort *) calloc (jh->wide*jh->clrs, 4);
841   merror (jh->row, "ljpeg_start()");
842   return zero_after_ff = 1;
843 }
844
845 void CLASS ljpeg_end (struct jhead *jh)
846 {
847   int c;
848   FORC4 if (jh->free[c]) free (jh->free[c]);
849   free (jh->row);
850 }
851
852 int CLASS ljpeg_diff (ushort *huff)
853 {
854   int len, diff;
855
856   len = gethuff(huff);
857   if (len == 16 && (!dng_version || dng_version >= 0x1010000))
858     return -32768;
859   diff = getbits(len);
860   if ((diff & (1 << (len-1))) == 0)
861     diff -= (1 << len) - 1;
862   return diff;
863 }
864
865 ushort * CLASS ljpeg_row (int jrow, struct jhead *jh)
866 {
867   int col, c, diff, pred, spred=0;
868   ushort mark=0, *row[3];
869
870   if (jrow * jh->wide % jh->restart == 0) {
871     FORC(6) jh->vpred[c] = 1 << (jh->bits-1);
872     if (jrow) {
873       fseek (ifp, -2, SEEK_CUR);
874       do mark = (mark << 8) + (c = fgetc(ifp));
875       while (c != EOF && mark >> 4 != 0xffd);
876     }
877     getbits(-1);
878   }
879   FORC3 row[c] = jh->row + jh->wide*jh->clrs*((jrow+c) & 1);
880   for (col=0; col < jh->wide; col++)
881     FORC(jh->clrs) {
882       diff = ljpeg_diff (jh->huff[c]);
883       if (jh->sraw && c <= jh->sraw && (col | c))
884                     pred = spred;
885       else if (col) pred = row[0][-jh->clrs];
886       else          pred = (jh->vpred[c] += diff) - diff;
887       if (jrow && col) switch (jh->psv) {
888         case 1: break;
889         case 2: pred = row[1][0];                                       break;
890         case 3: pred = row[1][-jh->clrs];                               break;
891         case 4: pred = pred +   row[1][0] - row[1][-jh->clrs];          break;
892         case 5: pred = pred + ((row[1][0] - row[1][-jh->clrs]) >> 1);   break;
893         case 6: pred = row[1][0] + ((pred - row[1][-jh->clrs]) >> 1);   break;
894         case 7: pred = (pred + row[1][0]) >> 1;                         break;
895         default: pred = 0;
896       }
897       if ((**row = pred + diff) >> jh->bits) derror();
898       if (c <= jh->sraw) spred = **row;
899       row[0]++; row[1]++;
900     }
901   return row[2];
902 }
903
904 void CLASS lossless_jpeg_load_raw()
905 {
906   int jwide, jrow, jcol, val, jidx, i, j, row=0, col=0;
907   struct jhead jh;
908   ushort *rp;
909
910   if (!ljpeg_start (&jh, 0)) return;
911   jwide = jh.wide * jh.clrs;
912
913   for (jrow=0; jrow < jh.high; jrow++) {
914     rp = ljpeg_row (jrow, &jh);
915     if (load_flags & 1)
916       row = jrow & 1 ? height-1-jrow/2 : jrow/2;
917     for (jcol=0; jcol < jwide; jcol++) {
918       val = curve[*rp++];
919       if (cr2_slice[0]) {
920         jidx = jrow*jwide + jcol;
921         i = jidx / (cr2_slice[1]*raw_height);
922         if ((j = i >= cr2_slice[0]))
923                  i  = cr2_slice[0];
924         jidx -= i * (cr2_slice[1]*raw_height);
925         row = jidx / cr2_slice[1+j];
926         col = jidx % cr2_slice[1+j] + i*cr2_slice[1];
927       }
928       if (raw_width == 3984 && (col -= 2) < 0)
929         col += (row--,raw_width);
930       if ((unsigned) row < raw_height) RAW(row,col) = val;
931       if (++col >= raw_width)
932         col = (row++,0);
933     }
934   }
935   ljpeg_end (&jh);
936 }
937
938 void CLASS canon_sraw_load_raw()
939 {
940   struct jhead jh;
941   short *rp=0, (*ip)[4];
942   int jwide, slice, scol, ecol, row, col, jrow=0, jcol=0, pix[3], c;
943   int v[3]={0,0,0}, ver, hue;
944   char *cp;
945
946   if (!ljpeg_start (&jh, 0) || jh.clrs < 4) return;
947   jwide = (jh.wide >>= 1) * jh.clrs;
948
949   for (ecol=slice=0; slice <= cr2_slice[0]; slice++) {
950     scol = ecol;
951     ecol += cr2_slice[1] * 2 / jh.clrs;
952     if (!cr2_slice[0] || ecol > raw_width-1) ecol = raw_width & -2;
953     for (row=0; row < height; row += (jh.clrs >> 1) - 1) {
954       ip = (short (*)[4]) image + row*width;
955       for (col=scol; col < ecol; col+=2, jcol+=jh.clrs) {
956         if ((jcol %= jwide) == 0)
957           rp = (short *) ljpeg_row (jrow++, &jh);
958         if (col >= width) continue;
959         FORC (jh.clrs-2)
960           ip[col + (c >> 1)*width + (c & 1)][0] = rp[jcol+c];
961         ip[col][1] = rp[jcol+jh.clrs-2] - 16384;
962         ip[col][2] = rp[jcol+jh.clrs-1] - 16384;
963       }
964     }
965   }
966   for (cp=model2; *cp && !isdigit(*cp); cp++);
967   sscanf (cp, "%d.%d.%d", v, v+1, v+2);
968   ver = (v[0]*1000 + v[1])*1000 + v[2];
969   hue = (jh.sraw+1) << 2;
970   if (unique_id >= 0x80000281 || (unique_id == 0x80000218 && ver > 1000006))
971     hue = jh.sraw << 1;
972   ip = (short (*)[4]) image;
973   rp = ip[0];
974   for (row=0; row < height; row++, ip+=width) {
975     if (row & (jh.sraw >> 1)) {
976       for (col=0; col < width; col+=2)
977         for (c=1; c < 3; c++)
978           if (row == height-1)
979                ip[col][c] =  ip[col-width][c];
980           else ip[col][c] = (ip[col-width][c] + ip[col+width][c] + 1) >> 1;
981     }
982     for (col=1; col < width; col+=2)
983       for (c=1; c < 3; c++)
984         if (col == width-1)
985              ip[col][c] =  ip[col-1][c];
986         else ip[col][c] = (ip[col-1][c] + ip[col+1][c] + 1) >> 1;
987   }
988   for ( ; rp < ip[0]; rp+=4) {
989     if (unique_id == 0x80000218 ||
990         unique_id == 0x80000250 ||
991         unique_id == 0x80000261 ||
992         unique_id == 0x80000281 ||
993         unique_id == 0x80000287) {
994       rp[1] = (rp[1] << 2) + hue;
995       rp[2] = (rp[2] << 2) + hue;
996       pix[0] = rp[0] + ((   50*rp[1] + 22929*rp[2]) >> 14);
997       pix[1] = rp[0] + ((-5640*rp[1] - 11751*rp[2]) >> 14);
998       pix[2] = rp[0] + ((29040*rp[1] -   101*rp[2]) >> 14);
999     } else {
1000       if (unique_id < 0x80000218) rp[0] -= 512;
1001       pix[0] = rp[0] + rp[2];
1002       pix[2] = rp[0] + rp[1];
1003       pix[1] = rp[0] + ((-778*rp[1] - (rp[2] << 11)) >> 12);
1004     }
1005     FORC3 rp[c] = (int)CLIP(pix[c] * sraw_mul[c] >> 10);
1006   }
1007   ljpeg_end (&jh);
1008   maximum = 0x3fff;
1009 }
1010
1011 void CLASS adobe_copy_pixel (unsigned row, unsigned col, ushort **rp)
1012 {
1013   int c;
1014
1015   if (tiff_samples == 2 && shot_select) (*rp)++;
1016   if (raw_image) {
1017     if (row < raw_height && col < raw_width)
1018       RAW(row,col) = curve[**rp];
1019     *rp += tiff_samples;
1020   } else {
1021     if (row < height && col < width)
1022       FORC(tiff_samples)
1023         image[row*width+col][c] = curve[(*rp)[c]];
1024     *rp += tiff_samples;
1025   }
1026   if (tiff_samples == 2 && shot_select) (*rp)--;
1027 }
1028
1029 void CLASS ljpeg_idct (struct jhead *jh)
1030 {
1031   int c, i, j, len, skip, coef;
1032   float work[3][8][8], *cs = ljpeg_cs;
1033   static const uchar zigzag[80] =
1034   {  0, 1, 8,16, 9, 2, 3,10,17,24,32,25,18,11, 4, 5,12,19,26,33,
1035     40,48,41,34,27,20,13, 6, 7,14,21,28,35,42,49,56,57,50,43,36,
1036     29,22,15,23,30,37,44,51,58,59,52,45,38,31,39,46,53,60,61,54,
1037     47,55,62,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63 };
1038
1039   if (!cs[0])
1040     FORC(106) cs[c] = cos((c & 31)*M_PI/16)/2;
1041   memset (work, 0, sizeof work);
1042   work[0][0][0] = jh->vpred[0] += ljpeg_diff (jh->huff[0]) * jh->quant[0];
1043   for (i=1; i < 64; i++ ) {
1044     len = gethuff (jh->huff[16]);
1045     i += skip = len >> 4;
1046     if (!(len &= 15) && skip < 15) break;
1047     coef = getbits(len);
1048     if ((coef & (1 << (len-1))) == 0)
1049       coef -= (1 << len) - 1;
1050     ((float *)work)[zigzag[i]] = coef * jh->quant[i];
1051   }
1052   FORC(8) work[0][0][c] *= M_SQRT1_2;
1053   FORC(8) work[0][c][0] *= M_SQRT1_2;
1054   for (i=0; i < 8; i++)
1055     for (j=0; j < 8; j++)
1056       FORC(8) work[1][i][j] += work[0][i][c] * cs[(j*2+1)*c];
1057   for (i=0; i < 8; i++)
1058     for (j=0; j < 8; j++)
1059       FORC(8) work[2][i][j] += work[1][c][j] * cs[(i*2+1)*c];
1060
1061   FORC(64) jh->idct[c] = CLIP(((float *)work[2])[c]+0.5);
1062 }
1063
1064 void CLASS lossless_dng_load_raw()
1065 {
1066   unsigned save, trow=0, tcol=0, jwide, jrow, jcol, row, col, i, j;
1067   struct jhead jh;
1068   ushort *rp;
1069
1070   while (trow < raw_height) {
1071     save = ftell(ifp);
1072     if (tile_length < INT_MAX)
1073       fseek (ifp, get4(), SEEK_SET);
1074     if (!ljpeg_start (&jh, 0)) break;
1075     jwide = jh.wide;
1076     if (filters) jwide *= jh.clrs;
1077     jwide /= MIN (is_raw, tiff_samples);
1078     switch (jh.algo) {
1079       case 0xc1:
1080         jh.vpred[0] = 16384;
1081         getbits(-1);
1082         for (jrow=0; jrow+7 < jh.high; jrow += 8) {
1083           for (jcol=0; jcol+7 < jh.wide; jcol += 8) {
1084             ljpeg_idct (&jh);
1085             rp = jh.idct;
1086             row = trow + jcol/tile_width + jrow*2;
1087             col = tcol + jcol%tile_width;
1088             for (i=0; i < 16; i+=2)
1089               for (j=0; j < 8; j++)
1090                 adobe_copy_pixel (row+i, col+j, &rp);
1091           }
1092         }
1093         break;
1094       case 0xc3:
1095         for (row=col=jrow=0; jrow < jh.high; jrow++) {
1096           rp = ljpeg_row (jrow, &jh);
1097           for (jcol=0; jcol < jwide; jcol++) {
1098             adobe_copy_pixel (trow+row, tcol+col, &rp);
1099             if (++col >= tile_width || col >= raw_width)
1100               row += 1 + (col = 0);
1101           }
1102         }
1103     }
1104     fseek (ifp, save+4, SEEK_SET);
1105     if ((tcol += tile_width) >= raw_width)
1106       trow += tile_length + (tcol = 0);
1107     ljpeg_end (&jh);
1108   }
1109 }
1110
1111 void CLASS packed_dng_load_raw()
1112 {
1113   ushort *pixel, *rp;
1114   int row, col;
1115
1116   pixel = (ushort *) calloc (raw_width, tiff_samples*sizeof *pixel);
1117   merror (pixel, "packed_dng_load_raw()");
1118   for (row=0; row < raw_height; row++) {
1119     if (tiff_bps == 16)
1120       read_shorts (pixel, raw_width * tiff_samples);
1121     else {
1122       getbits(-1);
1123       for (col=0; col < raw_width * tiff_samples; col++)
1124         pixel[col] = getbits(tiff_bps);
1125     }
1126     for (rp=pixel, col=0; col < raw_width; col++)
1127       adobe_copy_pixel (row, col, &rp);
1128   }
1129   free (pixel);
1130 }
1131
1132 void CLASS pentax_load_raw()
1133 {
1134   ushort bit[2][15], huff[4097];
1135   int dep, row, col, diff, c, i;
1136   ushort vpred[2][2] = {{0,0},{0,0}}, hpred[2];
1137
1138   fseek (ifp, meta_offset, SEEK_SET);
1139   dep = (get2() + 12) & 15;
1140   fseek (ifp, 12, SEEK_CUR);
1141   FORC(dep) bit[0][c] = get2();
1142   FORC(dep) bit[1][c] = fgetc(ifp);
1143   FORC(dep)
1144     for (i=bit[0][c]; i <= ((bit[0][c]+(4096 >> bit[1][c])-1) & 4095); )
1145       huff[++i] = bit[1][c] << 8 | c;
1146   huff[0] = 12;
1147   fseek (ifp, data_offset, SEEK_SET);
1148   getbits(-1);
1149   for (row=0; row < raw_height; row++)
1150     for (col=0; col < raw_width; col++) {
1151       diff = ljpeg_diff (huff);
1152       if (col < 2) hpred[col] = vpred[row & 1][col] += diff;
1153       else         hpred[col & 1] += diff;
1154       RAW(row,col) = hpred[col & 1];
1155       if (hpred[col & 1] >> tiff_bps) derror();
1156     }
1157 }
1158
1159 void CLASS nikon_load_raw()
1160 {
1161   static const uchar nikon_tree[][32] = {
1162     { 0,1,5,1,1,1,1,1,1,2,0,0,0,0,0,0,  /* 12-bit lossy */
1163       5,4,3,6,2,7,1,0,8,9,11,10,12 },
1164     { 0,1,5,1,1,1,1,1,1,2,0,0,0,0,0,0,  /* 12-bit lossy after split */
1165       0x39,0x5a,0x38,0x27,0x16,5,4,3,2,1,0,11,12,12 },
1166     { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,  /* 12-bit lossless */
1167       5,4,6,3,7,2,8,1,9,0,10,11,12 },
1168     { 0,1,4,3,1,1,1,1,1,2,0,0,0,0,0,0,  /* 14-bit lossy */
1169       5,6,4,7,8,3,9,2,1,0,10,11,12,13,14 },
1170     { 0,1,5,1,1,1,1,1,1,1,2,0,0,0,0,0,  /* 14-bit lossy after split */
1171       8,0x5c,0x4b,0x3a,0x29,7,6,5,4,3,2,1,0,13,14 },
1172     { 0,1,4,2,2,3,1,2,0,0,0,0,0,0,0,0,  /* 14-bit lossless */
1173       7,6,8,5,9,4,10,3,11,12,2,0,1,13,14 } };
1174   ushort *huff, ver0, ver1, vpred[2][2], hpred[2], csize;
1175   int i, min, max, step=0, tree=0, split=0, row, col, len, shl, diff;
1176
1177   fseek (ifp, meta_offset, SEEK_SET);
1178   ver0 = fgetc(ifp);
1179   ver1 = fgetc(ifp);
1180   if (ver0 == 0x49 || ver1 == 0x58)
1181     fseek (ifp, 2110, SEEK_CUR);
1182   if (ver0 == 0x46) tree = 2;
1183   if (tiff_bps == 14) tree += 3;
1184   read_shorts (vpred[0], 4);
1185   max = 1 << tiff_bps & 0x7fff;
1186   if ((csize = get2()) > 1)
1187     step = max / (csize-1);
1188   if (ver0 == 0x44 && ver1 == 0x20 && step > 0) {
1189     for (i=0; i < csize; i++)
1190       curve[i*step] = get2();
1191     for (i=0; i < max; i++)
1192       curve[i] = ( curve[i-i%step]*(step-i%step) +
1193                    curve[i-i%step+step]*(i%step) ) / step;
1194     fseek (ifp, meta_offset+562, SEEK_SET);
1195     split = get2();
1196   } else if (ver0 != 0x46 && csize <= 0x4001)
1197     read_shorts (curve, max=csize);
1198   while (curve[max-2] == curve[max-1]) max--;
1199   huff = make_decoder (nikon_tree[tree]);
1200   fseek (ifp, data_offset, SEEK_SET);
1201   getbits(-1);
1202   for (min=row=0; row < height; row++) {
1203     if (split && row == split) {
1204       free (huff);
1205       huff = make_decoder (nikon_tree[tree+1]);
1206       max += (min = 16) << 1;
1207     }
1208     for (col=0; col < raw_width; col++) {
1209       i = gethuff(huff);
1210       len = i & 15;
1211       shl = i >> 4;
1212       diff = ((getbits(len-shl) << 1) + 1) << shl >> 1;
1213       if ((diff & (1 << (len-1))) == 0)
1214         diff -= (1 << len) - !shl;
1215       if (col < 2) hpred[col] = vpred[row & 1][col] += diff;
1216       else         hpred[col & 1] += diff;
1217       if ((ushort)(hpred[col & 1] + min) >= max) derror();
1218       RAW(row,col) = curve[LIM((short)hpred[col & 1],0,0x3fff)];
1219     }
1220   }
1221   free (huff);
1222 }
1223
1224 void CLASS nikon_yuv_load_raw()
1225 {
1226   int row, col, yuv[4], rgb[3], b, c;
1227   UINT64 bitbuf=0;
1228
1229   for (row=0; row < raw_height; row++)
1230     for (col=0; col < raw_width; col++) {
1231       if (!(b = col & 1)) {
1232         bitbuf = 0;
1233         FORC(6) bitbuf |= (UINT64) fgetc(ifp) << c*8;
1234         FORC(4) yuv[c] = (bitbuf >> c*12 & 0xfff) - (c >> 1 << 11);
1235       }
1236       rgb[0] = yuv[b] + 1.370705*yuv[3];
1237       rgb[1] = yuv[b] - 0.337633*yuv[2] - 0.698001*yuv[3];
1238       rgb[2] = yuv[b] + 1.732446*yuv[2];
1239       FORC3 image[row*width+col][c] = curve[LIM(rgb[c],0,0xfff)] / cam_mul[c];
1240     }
1241 }
1242
1243 /*
1244    Returns 1 for a Coolpix 995, 0 for anything else.
1245  */
1246 int CLASS nikon_e995()
1247 {
1248   int i, histo[256];
1249   const uchar often[] = { 0x00, 0x55, 0xaa, 0xff };
1250
1251   memset (histo, 0, sizeof histo);
1252   fseek (ifp, -2000, SEEK_END);
1253   for (i=0; i < 2000; i++)
1254     histo[fgetc(ifp)]++;
1255   for (i=0; i < 4; i++)
1256     if (histo[often[i]] < 200)
1257       return 0;
1258   return 1;
1259 }
1260
1261 /*
1262    Returns 1 for a Coolpix 2100, 0 for anything else.
1263  */
1264 int CLASS nikon_e2100()
1265 {
1266   uchar t[12];
1267   int i;
1268
1269   fseek (ifp, 0, SEEK_SET);
1270   for (i=0; i < 1024; i++) {
1271     fread (t, 1, 12, ifp);
1272     if (((t[2] & t[4] & t[7] & t[9]) >> 4
1273         & t[1] & t[6] & t[8] & t[11] & 3) != 3)
1274       return 0;
1275   }
1276   return 1;
1277 }
1278
1279 void CLASS nikon_3700()
1280 {
1281   int bits, i;
1282   uchar dp[24];
1283   static const struct {
1284     int bits;
1285     char make[12], model[15];
1286   } table[] = {
1287     { 0x00, "Pentax",  "Optio 33WR" },
1288     { 0x03, "Nikon",   "E3200" },
1289     { 0x32, "Nikon",   "E3700" },
1290     { 0x33, "Olympus", "C740UZ" } };
1291
1292   fseek (ifp, 3072, SEEK_SET);
1293   fread (dp, 1, 24, ifp);
1294   bits = (dp[8] & 3) << 4 | (dp[20] & 3);
1295   for (i=0; i < sizeof table / sizeof *table; i++)
1296     if (bits == table[i].bits) {
1297       strcpy (make,  table[i].make );
1298       strcpy (model, table[i].model);
1299     }
1300 }
1301
1302 /*
1303    Separates a Minolta DiMAGE Z2 from a Nikon E4300.
1304  */
1305 int CLASS minolta_z2()
1306 {
1307   int i, nz;
1308   char tail[424];
1309
1310   fseek (ifp, -sizeof tail, SEEK_END);
1311   fread (tail, 1, sizeof tail, ifp);
1312   for (nz=i=0; i < sizeof tail; i++)
1313     if (tail[i]) nz++;
1314   return nz > 20;
1315 }
1316
1317 void CLASS ppm_thumb()
1318 {
1319   char *thumb;
1320   thumb_length = thumb_width*thumb_height*3;
1321   thumb = (char *) malloc (thumb_length);
1322   merror (thumb, "ppm_thumb()");
1323   fprintf (ofp, "P6\n%d %d\n255\n", thumb_width, thumb_height);
1324   fread  (thumb, 1, thumb_length, ifp);
1325   fwrite (thumb, 1, thumb_length, ofp);
1326   free (thumb);
1327 }
1328
1329 void CLASS ppm16_thumb()
1330 {
1331   int i;
1332   char *thumb;
1333   thumb_length = thumb_width*thumb_height*3;
1334   thumb = (char *) calloc (thumb_length, 2);
1335   merror (thumb, "ppm16_thumb()");
1336   read_shorts ((ushort *) thumb, thumb_length);
1337   for (i=0; i < thumb_length; i++)
1338     thumb[i] = ((ushort *) thumb)[i] >> 8;
1339   fprintf (ofp, "P6\n%d %d\n255\n", thumb_width, thumb_height);
1340   fwrite (thumb, 1, thumb_length, ofp);
1341   free (thumb);
1342 }
1343
1344 void CLASS layer_thumb()
1345 {
1346   int i, c;
1347   char *thumb, map[][4] = { "012","102" };
1348
1349   colors = thumb_misc >> 5 & 7;
1350   thumb_length = thumb_width*thumb_height;
1351   thumb = (char *) calloc (colors, thumb_length);
1352   merror (thumb, "layer_thumb()");
1353   fprintf (ofp, "P%d\n%d %d\n255\n",
1354         5 + (colors >> 1), thumb_width, thumb_height);
1355   fread (thumb, thumb_length, colors, ifp);
1356   for (i=0; i < thumb_length; i++)
1357     FORCC putc (thumb[i+thumb_length*(map[thumb_misc >> 8][c]-'0')], ofp);
1358   free (thumb);
1359 }
1360
1361 void CLASS rollei_thumb()
1362 {
1363   unsigned i;
1364   ushort *thumb;
1365
1366   thumb_length = thumb_width * thumb_height;
1367   thumb = (ushort *) calloc (thumb_length, 2);
1368   merror (thumb, "rollei_thumb()");
1369   fprintf (ofp, "P6\n%d %d\n255\n", thumb_width, thumb_height);
1370   read_shorts (thumb, thumb_length);
1371   for (i=0; i < thumb_length; i++) {
1372     putc (thumb[i] << 3, ofp);
1373     putc (thumb[i] >> 5  << 2, ofp);
1374     putc (thumb[i] >> 11 << 3, ofp);
1375   }
1376   free (thumb);
1377 }
1378
1379 void CLASS rollei_load_raw()
1380 {
1381   uchar pixel[10];
1382   unsigned iten=0, isix, i, buffer=0, todo[16];
1383
1384   isix = raw_width * raw_height * 5 / 8;
1385   while (fread (pixel, 1, 10, ifp) == 10) {
1386     for (i=0; i < 10; i+=2) {
1387       todo[i]   = iten++;
1388       todo[i+1] = pixel[i] << 8 | pixel[i+1];
1389       buffer    = pixel[i] >> 2 | buffer << 6;
1390     }
1391     for (   ; i < 16; i+=2) {
1392       todo[i]   = isix++;
1393       todo[i+1] = buffer >> (14-i)*5;
1394     }
1395     for (i=0; i < 16; i+=2)
1396       raw_image[todo[i]] = (todo[i+1] & 0x3ff);
1397   }
1398   maximum = 0x3ff;
1399 }
1400
1401 int CLASS raw (unsigned row, unsigned col)
1402 {
1403   return (row < raw_height && col < raw_width) ? RAW(row,col) : 0;
1404 }
1405
1406 void CLASS phase_one_flat_field (int is_float, int nc)
1407 {
1408   ushort head[8];
1409   unsigned wide, high, y, x, c, rend, cend, row, col;
1410   float *mrow, num, mult[4];
1411
1412   read_shorts (head, 8);
1413   if (head[2] * head[3] * head[4] * head[5] == 0) return;
1414   wide = head[2] / head[4] + (head[2] % head[4] != 0);
1415   high = head[3] / head[5] + (head[3] % head[5] != 0);
1416   mrow = (float *) calloc (nc*wide, sizeof *mrow);
1417   merror (mrow, "phase_one_flat_field()");
1418   for (y=0; y < high; y++) {
1419     for (x=0; x < wide; x++)
1420       for (c=0; c < nc; c+=2) {
1421         num = is_float ? getreal(11) : get2()/32768.0;
1422         if (y==0) mrow[c*wide+x] = num;
1423         else mrow[(c+1)*wide+x] = (num - mrow[c*wide+x]) / head[5];
1424       }
1425     if (y==0) continue;
1426     rend = head[1] + y*head[5];
1427     for (row = rend-head[5];
1428          row < raw_height && row < rend &&
1429          row < head[1]+head[3]-head[5]; row++) {
1430       for (x=1; x < wide; x++) {
1431         for (c=0; c < nc; c+=2) {
1432           mult[c] = mrow[c*wide+x-1];
1433           mult[c+1] = (mrow[c*wide+x] - mult[c]) / head[4];
1434         }
1435         cend = head[0] + x*head[4];
1436         for (col = cend-head[4];
1437              col < raw_width &&
1438              col < cend && col < head[0]+head[2]-head[4]; col++) {
1439           c = nc > 2 ? FC(row-top_margin,col-left_margin) : 0;
1440           if (!(c & 1)) {
1441             c = RAW(row,col) * mult[c];
1442             RAW(row,col) = LIM(c,0,65535);
1443           }
1444           for (c=0; c < nc; c+=2)
1445             mult[c] += mult[c+1];
1446         }
1447       }
1448       for (x=0; x < wide; x++)
1449         for (c=0; c < nc; c+=2)
1450           mrow[c*wide+x] += mrow[(c+1)*wide+x];
1451     }
1452   }
1453   free (mrow);
1454 }
1455
1456 void CLASS phase_one_correct()
1457 {
1458   unsigned entries, tag, data, save, col, row, type;
1459   int len, i, j, k, cip, val[4], dev[4], sum, max;
1460   int head[9], diff, mindiff=INT_MAX, off_412=0;
1461   static const signed char dir[12][2] =
1462     { {-1,-1}, {-1,1}, {1,-1}, {1,1}, {-2,0}, {0,-2}, {0,2}, {2,0},
1463       {-2,-2}, {-2,2}, {2,-2}, {2,2} };
1464   float poly[8], num, cfrac, frac, mult[2], *yval[2];
1465   ushort *xval[2];
1466   int qmult_applied = 0, qlin_applied = 0;
1467
1468   if (half_size || !meta_length) return;
1469   if (verbose) fprintf (stderr,_("Phase One correction...\n"));
1470   fseek (ifp, meta_offset, SEEK_SET);
1471   order = get2();
1472   fseek (ifp, 6, SEEK_CUR);
1473   fseek (ifp, meta_offset+get4(), SEEK_SET);
1474   entries = get4();  get4();
1475   while (entries--) {
1476     tag  = get4();
1477     len  = get4();
1478     data = get4();
1479     save = ftell(ifp);
1480     fseek (ifp, meta_offset+data, SEEK_SET);
1481     if (tag == 0x419) {                         /* Polynomial curve */
1482       for (get4(), i=0; i < 8; i++)
1483         poly[i] = getreal(11);
1484       poly[3] += (ph1.tag_210 - poly[7]) * poly[6] + 1;
1485       for (i=0; i < 0x10000; i++) {
1486         num = (poly[5]*i + poly[3])*i + poly[1];
1487         curve[i] = LIM(num,0,65535);
1488       } goto apply;                             /* apply to right half */
1489     } else if (tag == 0x41a) {                  /* Polynomial curve */
1490       for (i=0; i < 4; i++)
1491         poly[i] = getreal(11);
1492       for (i=0; i < 0x10000; i++) {
1493         for (num=0, j=4; j--; )
1494           num = num * i + poly[j];
1495         curve[i] = LIM(num+i,0,65535);
1496       } apply:                                  /* apply to whole image */
1497       for (row=0; row < raw_height; row++)
1498         for (col = (tag & 1)*ph1.split_col; col < raw_width; col++)
1499           RAW(row,col) = curve[RAW(row,col)];
1500     } else if (tag == 0x400) {                  /* Sensor defects */
1501       while ((len -= 8) >= 0) {
1502         col  = get2();
1503         row  = get2();
1504         type = get2(); get2();
1505         if (col >= raw_width) continue;
1506         if (type == 131 || type == 137)         /* Bad column */
1507           for (row=0; row < raw_height; row++)
1508             if (FC(row-top_margin,col-left_margin) == 1) {
1509               for (sum=i=0; i < 4; i++)
1510                 sum += val[i] = raw (row+dir[i][0], col+dir[i][1]);
1511               for (max=i=0; i < 4; i++) {
1512                 dev[i] = abs((val[i] << 2) - sum);
1513                 if (dev[max] < dev[i]) max = i;
1514               }
1515               RAW(row,col) = (sum - val[max])/3.0 + 0.5;
1516             } else {
1517               for (sum=0, i=8; i < 12; i++)
1518                 sum += raw (row+dir[i][0], col+dir[i][1]);
1519               RAW(row,col) = 0.5 + sum * 0.0732233 +
1520                 (raw(row,col-2) + raw(row,col+2)) * 0.3535534;
1521             }
1522         else if (type == 129) {                 /* Bad pixel */
1523           if (row >= raw_height) continue;
1524           j = (FC(row-top_margin,col-left_margin) != 1) * 4;
1525           for (sum=0, i=j; i < j+8; i++)
1526             sum += raw (row+dir[i][0], col+dir[i][1]);
1527           RAW(row,col) = (sum + 4) >> 3;
1528         }
1529       }
1530     } else if (tag == 0x401) {                  /* All-color flat fields */
1531       phase_one_flat_field (1, 2);
1532     } else if (tag == 0x416 || tag == 0x410) {
1533       phase_one_flat_field (0, 2);
1534     } else if (tag == 0x40b) {                  /* Red+blue flat field */
1535       phase_one_flat_field (0, 4);
1536     } else if (tag == 0x412) {
1537       fseek (ifp, 36, SEEK_CUR);
1538       diff = abs (get2() - ph1.tag_21a);
1539       if (mindiff > diff) {
1540         mindiff = diff;
1541         off_412 = ftell(ifp) - 38;
1542       }
1543     } else if (tag == 0x41f && !qlin_applied) { /* Quadrant linearization */
1544       ushort lc[2][2][16], ref[16];
1545       int qr, qc;
1546       for (qr = 0; qr < 2; qr++)
1547         for (qc = 0; qc < 2; qc++)
1548           for (i = 0; i < 16; i++)
1549             lc[qr][qc][i] = get4();
1550       for (i = 0; i < 16; i++) {
1551         int v = 0;
1552         for (qr = 0; qr < 2; qr++)
1553           for (qc = 0; qc < 2; qc++)
1554             v += lc[qr][qc][i];
1555         ref[i] = (v + 2) >> 2;
1556       }
1557       for (qr = 0; qr < 2; qr++) {
1558         for (qc = 0; qc < 2; qc++) {
1559           int cx[19], cf[19];
1560           for (i = 0; i < 16; i++) {
1561             cx[1+i] = lc[qr][qc][i];
1562             cf[1+i] = ref[i];
1563           }
1564           cx[0] = cf[0] = 0;
1565           cx[17] = cf[17] = ((unsigned) ref[15] * 65535) / lc[qr][qc][15];
1566           cx[18] = cf[18] = 65535;
1567           cubic_spline(cx, cf, 19);
1568           for (row = (qr ? ph1.split_row : 0);
1569                row < (qr ? raw_height : ph1.split_row); row++)
1570             for (col = (qc ? ph1.split_col : 0);
1571                  col < (qc ? raw_width : ph1.split_col); col++)
1572               RAW(row,col) = curve[RAW(row,col)];
1573         }
1574       }
1575       qlin_applied = 1;
1576     } else if (tag == 0x41e && !qmult_applied) { /* Quadrant multipliers */
1577       float qmult[2][2] = { { 1, 1 }, { 1, 1 } };
1578       get4(); get4(); get4(); get4();
1579       qmult[0][0] = 1.0 + getreal(11);
1580       get4(); get4(); get4(); get4(); get4();
1581       qmult[0][1] = 1.0 + getreal(11);
1582       get4(); get4(); get4();
1583       qmult[1][0] = 1.0 + getreal(11);
1584       get4(); get4(); get4();
1585       qmult[1][1] = 1.0 + getreal(11);
1586       for (row=0; row < raw_height; row++)
1587         for (col=0; col < raw_width; col++) {
1588           i = qmult[row >= ph1.split_row][col >= ph1.split_col] * RAW(row,col);
1589           RAW(row,col) = LIM(i,0,65535);
1590         }
1591       qmult_applied = 1;
1592     } else if (tag == 0x431 && !qmult_applied) { /* Quadrant combined */
1593       ushort lc[2][2][7], ref[7];
1594       int qr, qc;
1595       for (i = 0; i < 7; i++)
1596         ref[i] = get4();
1597       for (qr = 0; qr < 2; qr++)
1598         for (qc = 0; qc < 2; qc++)
1599           for (i = 0; i < 7; i++)
1600             lc[qr][qc][i] = get4();
1601       for (qr = 0; qr < 2; qr++) {
1602         for (qc = 0; qc < 2; qc++) {
1603           int cx[9], cf[9];
1604           for (i = 0; i < 7; i++) {
1605             cx[1+i] = ref[i];
1606             cf[1+i] = ((unsigned) ref[i] * lc[qr][qc][i]) / 10000;
1607           }
1608           cx[0] = cf[0] = 0;
1609           cx[8] = cf[8] = 65535;
1610           cubic_spline(cx, cf, 9);
1611           for (row = (qr ? ph1.split_row : 0);
1612                row < (qr ? raw_height : ph1.split_row); row++)
1613             for (col = (qc ? ph1.split_col : 0);
1614                  col < (qc ? raw_width : ph1.split_col); col++)
1615               RAW(row,col) = curve[RAW(row,col)];
1616         }
1617       }
1618       qmult_applied = 1;
1619       qlin_applied = 1;
1620     }
1621     fseek (ifp, save, SEEK_SET);
1622   }
1623   if (off_412) {
1624     fseek (ifp, off_412, SEEK_SET);
1625     for (i=0; i < 9; i++) head[i] = get4() & 0x7fff;
1626     yval[0] = (float *) calloc (head[1]*head[3] + head[2]*head[4], 6);
1627     merror (yval[0], "phase_one_correct()");
1628     yval[1] = (float  *) (yval[0] + head[1]*head[3]);
1629     xval[0] = (ushort *) (yval[1] + head[2]*head[4]);
1630     xval[1] = (ushort *) (xval[0] + head[1]*head[3]);
1631     get2();
1632     for (i=0; i < 2; i++)
1633       for (j=0; j < head[i+1]*head[i+3]; j++)
1634         yval[i][j] = getreal(11);
1635     for (i=0; i < 2; i++)
1636       for (j=0; j < head[i+1]*head[i+3]; j++)
1637         xval[i][j] = get2();
1638     for (row=0; row < raw_height; row++)
1639       for (col=0; col < raw_width; col++) {
1640         cfrac = (float) col * head[3] / raw_width;
1641         cfrac -= cip = cfrac;
1642         num = RAW(row,col) * 0.5;
1643         for (i=cip; i < cip+2; i++) {
1644           for (k=j=0; j < head[1]; j++)
1645             if (num < xval[0][k = head[1]*i+j]) break;
1646           frac = (j == 0 || j == head[1]) ? 0 :
1647                 (xval[0][k] - num) / (xval[0][k] - xval[0][k-1]);
1648           mult[i-cip] = yval[0][k-1] * frac + yval[0][k] * (1-frac);
1649         }
1650         i = ((mult[0] * (1-cfrac) + mult[1] * cfrac) * row + num) * 2;
1651         RAW(row,col) = LIM(i,0,65535);
1652       }
1653     free (yval[0]);
1654   }
1655 }
1656
1657 void CLASS phase_one_load_raw()
1658 {
1659   int a, b, i;
1660   ushort akey, bkey, mask;
1661
1662   fseek (ifp, ph1.key_off, SEEK_SET);
1663   akey = get2();
1664   bkey = get2();
1665   mask = ph1.format == 1 ? 0x5555:0x1354;
1666   fseek (ifp, data_offset, SEEK_SET);
1667   read_shorts (raw_image, raw_width*raw_height);
1668   if (ph1.format)
1669     for (i=0; i < raw_width*raw_height; i+=2) {
1670       a = raw_image[i+0] ^ akey;
1671       b = raw_image[i+1] ^ bkey;
1672       raw_image[i+0] = (a & mask) | (b & ~mask);
1673       raw_image[i+1] = (b & mask) | (a & ~mask);
1674     }
1675 }
1676
1677 unsigned CLASS ph1_bithuff (int nbits, ushort *huff)
1678 {
1679   UINT64 bitbuf;
1680   int vbits;
1681   unsigned c;
1682
1683   if (nbits == -1)
1684     return ph1_bitbuf = ph1_vbits = 0;
1685   if (nbits == 0) return 0;
1686   bitbuf = ph1_bitbuf;  vbits = ph1_vbits;
1687   if (vbits < nbits) {
1688     bitbuf = bitbuf << 32 | get4();
1689     vbits += 32;
1690   }
1691   c = bitbuf << (64-vbits) >> (64-nbits);
1692   if (huff) {
1693     nbits = huff[c] >> 8;
1694     c = (uchar) huff[c];
1695   }
1696   vbits -= nbits;
1697   ph1_bitbuf = bitbuf;  ph1_vbits = vbits;
1698   return c;
1699 }
1700 #define ph1_bits(n) ph1_bithuff(n,0)
1701 #define ph1_huff(h) ph1_bithuff(*h,h+1)
1702
1703 void CLASS phase_one_load_raw_c()
1704 {
1705   static const int length[] = { 8,7,6,9,11,10,5,12,14,13 };
1706   int *offset, len[2], pred[2], row, col, i, j;
1707   ushort *pixel;
1708   short (*cblack)[2], (*rblack)[2];
1709
1710   pixel = (ushort *) calloc (raw_width*3 + raw_height*4, 2);
1711   merror (pixel, "phase_one_load_raw_c()");
1712   offset = (int *) (pixel + raw_width);
1713   fseek (ifp, strip_offset, SEEK_SET);
1714   for (row=0; row < raw_height; row++)
1715     offset[row] = get4();
1716   cblack = (short (*)[2]) (offset + raw_height);
1717   fseek (ifp, ph1.black_col, SEEK_SET);
1718   if (ph1.black_col)
1719     read_shorts ((ushort *) cblack[0], raw_height*2);
1720   rblack = cblack + raw_height;
1721   fseek (ifp, ph1.black_row, SEEK_SET);
1722   if (ph1.black_row)
1723     read_shorts ((ushort *) rblack[0], raw_width*2);
1724   for (i=0; i < 256; i++)
1725     curve[i] = i*i / 3.969 + 0.5;
1726   for (row=0; row < raw_height; row++) {
1727     fseek (ifp, data_offset + offset[row], SEEK_SET);
1728     ph1_bits(-1);
1729     pred[0] = pred[1] = 0;
1730     for (col=0; col < raw_width; col++) {
1731       if (col >= (raw_width & -8))
1732         len[0] = len[1] = 14;
1733       else if ((col & 7) == 0)
1734         for (i=0; i < 2; i++) {
1735           for (j=0; j < 5 && !ph1_bits(1); j++);
1736           if (j--) len[i] = length[j*2 + ph1_bits(1)];
1737         }
1738       if ((i = len[col & 1]) == 14)
1739         pixel[col] = pred[col & 1] = ph1_bits(16);
1740       else
1741         pixel[col] = pred[col & 1] += ph1_bits(i) + 1 - (1 << (i - 1));
1742       if (pred[col & 1] >> 16) derror();
1743       if (ph1.format == 5 && pixel[col] < 256)
1744         pixel[col] = curve[pixel[col]];
1745     }
1746     for (col=0; col < raw_width; col++) {
1747       i = (pixel[col] << 2*(ph1.format != 8)) - ph1.black
1748         + cblack[row][col >= ph1.split_col]
1749         + rblack[col][row >= ph1.split_row];
1750       if (i > 0) RAW(row,col) = i;
1751     }
1752   }
1753   free (pixel);
1754   maximum = 0xfffc - ph1.black;
1755 }
1756
1757 void CLASS hasselblad_load_raw()
1758 {
1759   struct jhead jh;
1760   int shot, row, col, *back[5], len[2], diff[12], pred, sh, f, s, c;
1761   unsigned upix, urow, ucol;
1762   ushort *ip;
1763
1764   if (!ljpeg_start (&jh, 0)) return;
1765   order = 0x4949;
1766   ph1_bits(-1);
1767   back[4] = (int *) calloc (raw_width, 3*sizeof **back);
1768   merror (back[4], "hasselblad_load_raw()");
1769   FORC3 back[c] = back[4] + c*raw_width;
1770   cblack[6] >>= sh = tiff_samples > 1;
1771   shot = LIM(shot_select, 1, tiff_samples) - 1;
1772   for (row=0; row < raw_height; row++) {
1773     FORC4 back[(c+3) & 3] = back[c];
1774     for (col=0; col < raw_width; col+=2) {
1775       for (s=0; s < tiff_samples*2; s+=2) {
1776         FORC(2) len[c] = ph1_huff(jh.huff[0]);
1777         FORC(2) {
1778           diff[s+c] = ph1_bits(len[c]);
1779           if ((diff[s+c] & (1 << (len[c]-1))) == 0)
1780             diff[s+c] -= (1 << len[c]) - 1;
1781           if (diff[s+c] == 65535) diff[s+c] = -32768;
1782         }
1783       }
1784       for (s=col; s < col+2; s++) {
1785         pred = 0x8000 + load_flags;
1786         if (col) pred = back[2][s-2];
1787         if (col && row > 1) switch (jh.psv) {
1788           case 11: pred += back[0][s]/2 - back[0][s-2]/2;  break;
1789         }
1790         f = (row & 1)*3 ^ ((col+s) & 1);
1791         FORC (tiff_samples) {
1792           pred += diff[(s & 1)*tiff_samples+c];
1793           upix = pred >> sh & 0xffff;
1794           if (raw_image && c == shot)
1795             RAW(row,s) = upix;
1796           if (image) {
1797             urow = row-top_margin  + (c & 1);
1798             ucol = col-left_margin - ((c >> 1) & 1);
1799             ip = &image[urow*width+ucol][f];
1800             if (urow < height && ucol < width)
1801               *ip = c < 4 ? upix : (*ip + upix) >> 1;
1802           }
1803         }
1804         back[2][s] = pred;
1805       }
1806     }
1807   }
1808   free (back[4]);
1809   ljpeg_end (&jh);
1810   if (image) mix_green = 1;
1811 }
1812
1813 void CLASS leaf_hdr_load_raw()
1814 {
1815   ushort *pixel=0;
1816   unsigned tile=0, r, c, row, col;
1817
1818   if (!filters) {
1819     pixel = (ushort *) calloc (raw_width, sizeof *pixel);
1820     merror (pixel, "leaf_hdr_load_raw()");
1821   }
1822   FORC(tiff_samples)
1823     for (r=0; r < raw_height; r++) {
1824       if (r % tile_length == 0) {
1825         fseek (ifp, data_offset + 4*tile++, SEEK_SET);
1826         fseek (ifp, get4(), SEEK_SET);
1827       }
1828       if (filters && c != shot_select) continue;
1829       if (filters) pixel = raw_image + r*raw_width;
1830       read_shorts (pixel, raw_width);
1831       if (!filters && (row = r - top_margin) < height)
1832         for (col=0; col < width; col++)
1833           image[row*width+col][c] = pixel[col+left_margin];
1834     }
1835   if (!filters) {
1836     maximum = 0xffff;
1837     raw_color = 1;
1838     free (pixel);
1839   }
1840 }
1841
1842 void CLASS unpacked_load_raw()
1843 {
1844   int row, col, bits=0;
1845
1846   while (1 << ++bits < maximum);
1847   read_shorts (raw_image, raw_width*raw_height);
1848   for (row=0; row < raw_height; row++)
1849     for (col=0; col < raw_width; col++)
1850       if ((RAW(row,col) >>= load_flags) >> bits
1851         && (unsigned) (row-top_margin) < height
1852         && (unsigned) (col-left_margin) < width) derror();
1853 }
1854
1855 void CLASS sinar_4shot_load_raw()
1856 {
1857   ushort *pixel;
1858   unsigned shot, row, col, r, c;
1859
1860   if (raw_image) {
1861     shot = LIM (shot_select, 1, 4) - 1;
1862     fseek (ifp, data_offset + shot*4, SEEK_SET);
1863     fseek (ifp, get4(), SEEK_SET);
1864     unpacked_load_raw();
1865     return;
1866   }
1867   pixel = (ushort *) calloc (raw_width, sizeof *pixel);
1868   merror (pixel, "sinar_4shot_load_raw()");
1869   for (shot=0; shot < 4; shot++) {
1870     fseek (ifp, data_offset + shot*4, SEEK_SET);
1871     fseek (ifp, get4(), SEEK_SET);
1872     for (row=0; row < raw_height; row++) {
1873       read_shorts (pixel, raw_width);
1874       if ((r = row-top_margin - (shot >> 1 & 1)) >= height) continue;
1875       for (col=0; col < raw_width; col++) {
1876         if ((c = col-left_margin - (shot & 1)) >= width) continue;
1877         image[r*width+c][(row & 1)*3 ^ (~col & 1)] = pixel[col];
1878       }
1879     }
1880   }
1881   free (pixel);
1882   mix_green = 1;
1883 }
1884
1885 void CLASS imacon_full_load_raw()
1886 {
1887   int row, col;
1888
1889   if (!image) return;
1890   for (row=0; row < height; row++)
1891     for (col=0; col < width; col++)
1892       read_shorts (image[row*width+col], 3);
1893 }
1894
1895 void CLASS packed_load_raw()
1896 {
1897   int vbits=0, bwide, rbits, bite, half, irow, row, col, val, i;
1898   UINT64 bitbuf=0;
1899
1900   bwide = raw_width * tiff_bps / 8;
1901   bwide += bwide & load_flags >> 9;
1902   rbits = bwide * 8 - raw_width * tiff_bps;
1903   if (load_flags & 1) bwide = bwide * 16 / 15;
1904   bite = 8 + (load_flags & 56);
1905   half = (raw_height+1) >> 1;
1906   for (irow=0; irow < raw_height; irow++) {
1907     row = irow;
1908     if (load_flags & 2 &&
1909         (row = irow % half * 2 + irow / half) == 1 &&
1910         load_flags & 4) {
1911       if (vbits=0, tiff_compress)
1912         fseek (ifp, data_offset - (-half*bwide & -2048), SEEK_SET);
1913       else {
1914         fseek (ifp, 0, SEEK_END);
1915         fseek (ifp, ftell(ifp) >> 3 << 2, SEEK_SET);
1916       }
1917     }
1918     for (col=0; col < raw_width; col++) {
1919       for (vbits -= tiff_bps; vbits < 0; vbits += bite) {
1920         bitbuf <<= bite;
1921         for (i=0; i < bite; i+=8)
1922           bitbuf |= ((UINT64) fgetc(ifp) << i);
1923       }
1924       val = bitbuf << (64-tiff_bps-vbits) >> (64-tiff_bps);
1925       RAW(row,col ^ (load_flags >> 6 & 3)) = val;
1926       if (load_flags & 1 && (col % 10) == 9 && fgetc(ifp) &&
1927         row < height+top_margin && col < width+left_margin) derror();
1928     }
1929     vbits -= rbits;
1930   }
1931 }
1932
1933 void CLASS nokia_load_raw()
1934 {
1935   uchar  *data,  *dp;
1936   int rev, dwide, row, col, c;
1937   double sum[]={0,0};
1938
1939   rev = 3 * (order == 0x4949);
1940   dwide = (raw_width * 5 + 1) / 4;
1941   data = (uchar *) malloc (dwide*2);
1942   merror (data, "nokia_load_raw()");
1943   for (row=0; row < raw_height; row++) {
1944     if (fread (data+dwide, 1, dwide, ifp) < dwide) derror();
1945     FORC(dwide) data[c] = data[dwide+(c ^ rev)];
1946     for (dp=data, col=0; col < raw_width; dp+=5, col+=4)
1947       FORC4 RAW(row,col+c) = (dp[c] << 2) | (dp[4] >> (c << 1) & 3);
1948   }
1949   free (data);
1950   maximum = 0x3ff;
1951   if (strcmp(make,"OmniVision")) return;
1952   row = raw_height/2;
1953   FORC(width-1) {
1954     sum[ c & 1] += SQR(RAW(row,c)-RAW(row+1,c+1));
1955     sum[~c & 1] += SQR(RAW(row+1,c)-RAW(row,c+1));
1956   }
1957   if (sum[1] > sum[0]) filters = 0x4b4b4b4b;
1958 }
1959
1960 void CLASS canon_rmf_load_raw()
1961 {
1962   int row, col, bits, orow, ocol, c;
1963
1964   for (row=0; row < raw_height; row++)
1965     for (col=0; col < raw_width-2; col+=3) {
1966       bits = get4();
1967       FORC3 {
1968         orow = row;
1969         if ((ocol = col+c-4) < 0) {
1970           ocol += raw_width;
1971           if ((orow -= 2) < 0)
1972             orow += raw_height;
1973         }
1974         RAW(orow,ocol) = curve[bits >> (10*c+2) & 0x3ff];
1975       }
1976     }
1977   maximum = curve[0x3ff];
1978 }
1979
1980 unsigned CLASS pana_bits (int nbits)
1981 {
1982   uchar *buf = pana_buf;
1983   int vbits = pana_vbits;
1984   int byte;
1985
1986   if (!nbits) return pana_vbits=0;
1987   if (!vbits) {
1988     fread (buf+load_flags, 1, 0x4000-load_flags, ifp);
1989     fread (buf, 1, load_flags, ifp);
1990   }
1991   vbits = (vbits - nbits) & 0x1ffff;
1992   byte = vbits >> 3 ^ 0x3ff0;
1993   pana_vbits = vbits;
1994   return (buf[byte] | buf[byte+1] << 8) >> (vbits & 7) & ~(-1 << nbits);
1995 }
1996
1997 void CLASS panasonic_load_raw()
1998 {
1999   int row, col, i, j, sh=0, pred[2], nonz[2];
2000
2001   pana_bits(0);
2002   for (row=0; row < height; row++)
2003     for (col=0; col < raw_width; col++) {
2004       if ((i = col % 14) == 0)
2005         pred[0] = pred[1] = nonz[0] = nonz[1] = 0;
2006       if (i % 3 == 2) sh = 4 >> (3 - pana_bits(2));
2007       if (nonz[i & 1]) {
2008         if ((j = pana_bits(8))) {
2009           if ((pred[i & 1] -= 0x80 << sh) < 0 || sh == 4)
2010                pred[i & 1] &= ~(-1 << sh);
2011           pred[i & 1] += j << sh;
2012         }
2013       } else if ((nonz[i & 1] = pana_bits(8)) || i > 11)
2014         pred[i & 1] = nonz[i & 1] << 4 | pana_bits(4);
2015       if ((RAW(row,col) = pred[col & 1]) > 4098 && col < width) derror();
2016     }
2017 }
2018
2019 void CLASS olympus_load_raw()
2020 {
2021   ushort huff[4096];
2022   int row, col, nbits, sign, low, high, i, c, w, n, nw;
2023   int acarry[2][3], *carry, pred, diff;
2024
2025   huff[n=0] = 0xc0c;
2026   for (i=12; i--; )
2027     FORC(2048 >> i) huff[++n] = (i+1) << 8 | i;
2028   fseek (ifp, 7, SEEK_CUR);
2029   getbits(-1);
2030   for (row=0; row < height; row++) {
2031     memset (acarry, 0, sizeof acarry);
2032     for (col=0; col < raw_width; col++) {
2033       carry = acarry[col & 1];
2034       i = 2 * (carry[2] < 3);
2035       for (nbits=2+i; (ushort) carry[0] >> (nbits+i); nbits++);
2036       low = (sign = getbits(3)) & 3;
2037       sign = sign << 29 >> 31;
2038       if ((high = getbithuff(12,huff)) == 12)
2039         high = getbits(16-nbits) >> 1;
2040       carry[0] = (high << nbits) | getbits(nbits);
2041       diff = (carry[0] ^ sign) + carry[1];
2042       carry[1] = (diff*3 + carry[1]) >> 5;
2043       carry[2] = carry[0] > 16 ? 0 : carry[2]+1;
2044       if (col >= width) continue;
2045       if (row < 2 && col < 2) pred = 0;
2046       else if (row < 2) pred = RAW(row,col-2);
2047       else if (col < 2) pred = RAW(row-2,col);
2048       else {
2049         w  = RAW(row,col-2);
2050         n  = RAW(row-2,col);
2051         nw = RAW(row-2,col-2);
2052         if ((w < nw && nw < n) || (n < nw && nw < w)) {
2053           if (ABS(w-nw) > 32 || ABS(n-nw) > 32)
2054             pred = w + n - nw;
2055           else pred = (w + n) >> 1;
2056         } else pred = ABS(w-nw) > ABS(n-nw) ? w : n;
2057       }
2058       if ((RAW(row,col) = pred + ((diff << 2) | low)) >> 12) derror();
2059     }
2060   }
2061 }
2062
2063 void CLASS canon_crx_load_raw()
2064 {
2065 }
2066
2067 void CLASS fuji_xtrans_load_raw()
2068 {
2069 }
2070
2071 void CLASS minolta_rd175_load_raw()
2072 {
2073   uchar pixel[768];
2074   unsigned irow, box, row, col;
2075
2076   for (irow=0; irow < 1481; irow++) {
2077     if (fread (pixel, 1, 768, ifp) < 768) derror();
2078     box = irow / 82;
2079     row = irow % 82 * 12 + ((box < 12) ? box | 1 : (box-12)*2);
2080     switch (irow) {
2081       case 1477: case 1479: continue;
2082       case 1476: row = 984; break;
2083       case 1480: row = 985; break;
2084       case 1478: row = 985; box = 1;
2085     }
2086     if ((box < 12) && (box & 1)) {
2087       for (col=0; col < 1533; col++, row ^= 1)
2088         if (col != 1) RAW(row,col) = (col+1) & 2 ?
2089                    pixel[col/2-1] + pixel[col/2+1] : pixel[col/2] << 1;
2090       RAW(row,1)    = pixel[1]   << 1;
2091       RAW(row,1533) = pixel[765] << 1;
2092     } else
2093       for (col=row & 1; col < 1534; col+=2)
2094         RAW(row,col) = pixel[col/2] << 1;
2095   }
2096   maximum = 0xff << 1;
2097 }
2098
2099 void CLASS quicktake_100_load_raw()
2100 {
2101   uchar pixel[484][644];
2102   static const short gstep[16] =
2103   { -89,-60,-44,-32,-22,-15,-8,-2,2,8,15,22,32,44,60,89 };
2104   static const short rstep[6][4] =
2105   { {  -3,-1,1,3  }, {  -5,-1,1,5  }, {  -8,-2,2,8  },
2106     { -13,-3,3,13 }, { -19,-4,4,19 }, { -28,-6,6,28 } };
2107   static const short curve[256] =
2108   { 0,1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,
2109     28,29,30,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,53,
2110     54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,74,75,76,77,78,
2111     79,80,81,82,83,84,86,88,90,92,94,97,99,101,103,105,107,110,112,114,116,
2112     118,120,123,125,127,129,131,134,136,138,140,142,144,147,149,151,153,155,
2113     158,160,162,164,166,168,171,173,175,177,179,181,184,186,188,190,192,195,
2114     197,199,201,203,205,208,210,212,214,216,218,221,223,226,230,235,239,244,
2115     248,252,257,261,265,270,274,278,283,287,291,296,300,305,309,313,318,322,
2116     326,331,335,339,344,348,352,357,361,365,370,374,379,383,387,392,396,400,
2117     405,409,413,418,422,426,431,435,440,444,448,453,457,461,466,470,474,479,
2118     483,487,492,496,500,508,519,531,542,553,564,575,587,598,609,620,631,643,
2119     654,665,676,687,698,710,721,732,743,754,766,777,788,799,810,822,833,844,
2120     855,866,878,889,900,911,922,933,945,956,967,978,989,1001,1012,1023 };
2121   int rb, row, col, sharp, val=0;
2122
2123   getbits(-1);
2124   memset (pixel, 0x80, sizeof pixel);
2125   for (row=2; row < height+2; row++) {
2126     for (col=2+(row & 1); col < width+2; col+=2) {
2127       val = ((pixel[row-1][col-1] + 2*pixel[row-1][col+1] +
2128                 pixel[row][col-2]) >> 2) + gstep[getbits(4)];
2129       pixel[row][col] = val = LIM(val,0,255);
2130       if (col < 4)
2131         pixel[row][col-2] = pixel[row+1][~row & 1] = val;
2132       if (row == 2)
2133         pixel[row-1][col+1] = pixel[row-1][col+3] = val;
2134     }
2135     pixel[row][col] = val;
2136   }
2137   for (rb=0; rb < 2; rb++)
2138     for (row=2+rb; row < height+2; row+=2)
2139       for (col=3-(row & 1); col < width+2; col+=2) {
2140         if (row < 4 || col < 4) sharp = 2;
2141         else {
2142           val = ABS(pixel[row-2][col] - pixel[row][col-2])
2143               + ABS(pixel[row-2][col] - pixel[row-2][col-2])
2144               + ABS(pixel[row][col-2] - pixel[row-2][col-2]);
2145           sharp = val <  4 ? 0 : val <  8 ? 1 : val < 16 ? 2 :
2146                   val < 32 ? 3 : val < 48 ? 4 : 5;
2147         }
2148         val = ((pixel[row-2][col] + pixel[row][col-2]) >> 1)
2149               + rstep[sharp][getbits(2)];
2150         pixel[row][col] = val = LIM(val,0,255);
2151         if (row < 4) pixel[row-2][col+2] = val;
2152         if (col < 4) pixel[row+2][col-2] = val;
2153       }
2154   for (row=2; row < height+2; row++)
2155     for (col=3-(row & 1); col < width+2; col+=2) {
2156       val = ((pixel[row][col-1] + (pixel[row][col] << 2) +
2157               pixel[row][col+1]) >> 1) - 0x100;
2158       pixel[row][col] = LIM(val,0,255);
2159     }
2160   for (row=0; row < height; row++)
2161     for (col=0; col < width; col++)
2162       RAW(row,col) = curve[pixel[row+2][col+2]];
2163   maximum = 0x3ff;
2164 }
2165
2166 #define radc_token(tree) ((signed char) getbithuff(8,huff[tree]))
2167
2168 #define FORYX for (y=1; y < 3; y++) for (x=col+1; x >= col; x--)
2169
2170 #define PREDICTOR (c ? (buf[c][y-1][x] + buf[c][y][x+1]) / 2 \
2171 : (buf[c][y-1][x+1] + 2*buf[c][y-1][x] + buf[c][y][x+1]) / 4)
2172
2173 void CLASS kodak_radc_load_raw()
2174 {
2175   static const char src[] = {
2176     1,1, 2,3, 3,4, 4,2, 5,7, 6,5, 7,6, 7,8,
2177     1,0, 2,1, 3,3, 4,4, 5,2, 6,7, 7,6, 8,5, 8,8,
2178     2,1, 2,3, 3,0, 3,2, 3,4, 4,6, 5,5, 6,7, 6,8,
2179     2,0, 2,1, 2,3, 3,2, 4,4, 5,6, 6,7, 7,5, 7,8,
2180     2,1, 2,4, 3,0, 3,2, 3,3, 4,7, 5,5, 6,6, 6,8,
2181     2,3, 3,1, 3,2, 3,4, 3,5, 3,6, 4,7, 5,0, 5,8,
2182     2,3, 2,6, 3,0, 3,1, 4,4, 4,5, 4,7, 5,2, 5,8,
2183     2,4, 2,7, 3,3, 3,6, 4,1, 4,2, 4,5, 5,0, 5,8,
2184     2,6, 3,1, 3,3, 3,5, 3,7, 3,8, 4,0, 5,2, 5,4,
2185     2,0, 2,1, 3,2, 3,3, 4,4, 4,5, 5,6, 5,7, 4,8,
2186     1,0, 2,2, 2,-2,
2187     1,-3, 1,3,
2188     2,-17, 2,-5, 2,5, 2,17,
2189     2,-7, 2,2, 2,9, 2,18,
2190     2,-18, 2,-9, 2,-2, 2,7,
2191     2,-28, 2,28, 3,-49, 3,-9, 3,9, 4,49, 5,-79, 5,79,
2192     2,-1, 2,13, 2,26, 3,39, 4,-16, 5,55, 6,-37, 6,76,
2193     2,-26, 2,-13, 2,1, 3,-39, 4,16, 5,-55, 6,-76, 6,37
2194   };
2195   ushort huff[19][256];
2196   int row, col, tree, nreps, rep, step, i, c, s, r, x, y, val;
2197   short last[3] = { 16,16,16 }, mul[3], buf[3][3][386];
2198   static const ushort pt[] =
2199     { 0,0, 1280,1344, 2320,3616, 3328,8000, 4095,16383, 65535,16383 };
2200
2201   for (i=2; i < 12; i+=2)
2202     for (c=pt[i-2]; c <= pt[i]; c++)
2203       curve[c] = (float)
2204         (c-pt[i-2]) / (pt[i]-pt[i-2]) * (pt[i+1]-pt[i-1]) + pt[i-1] + 0.5;
2205   for (s=i=0; i < sizeof src; i+=2)
2206     FORC(256 >> src[i])
2207       ((ushort *)huff)[s++] = src[i] << 8 | (uchar) src[i+1];
2208   s = kodak_cbpp == 243 ? 2 : 3;
2209   FORC(256) huff[18][c] = (8-s) << 8 | c >> s << s | 1 << (s-1);
2210   getbits(-1);
2211   for (i=0; i < sizeof(buf)/sizeof(short); i++)
2212     ((short *)buf)[i] = 2048;
2213   for (row=0; row < height; row+=4) {
2214     FORC3 mul[c] = getbits(6);
2215     FORC3 {
2216       val = ((0x1000000/last[c] + 0x7ff) >> 12) * mul[c];
2217       s = val > 65564 ? 10:12;
2218       x = ~(-1 << (s-1));
2219       val <<= 12-s;
2220       for (i=0; i < sizeof(buf[0])/sizeof(short); i++)
2221         ((short *)buf[c])[i] = (((short *)buf[c])[i] * val + x) >> s;
2222       last[c] = mul[c];
2223       for (r=0; r <= !c; r++) {
2224         buf[c][1][width/2] = buf[c][2][width/2] = mul[c] << 7;
2225         for (tree=1, col=width/2; col > 0; ) {
2226           if ((tree = radc_token(tree))) {
2227             col -= 2;
2228             if (tree == 8)
2229               FORYX buf[c][y][x] = (uchar) radc_token(18) * mul[c];
2230             else
2231               FORYX buf[c][y][x] = radc_token(tree+10) * 16 + PREDICTOR;
2232           } else
2233             do {
2234               nreps = (col > 2) ? radc_token(9) + 1 : 1;
2235               for (rep=0; rep < 8 && rep < nreps && col > 0; rep++) {
2236                 col -= 2;
2237                 FORYX buf[c][y][x] = PREDICTOR;
2238                 if (rep & 1) {
2239                   step = radc_token(10) << 4;
2240                   FORYX buf[c][y][x] += step;
2241                 }
2242               }
2243             } while (nreps == 9);
2244         }
2245         for (y=0; y < 2; y++)
2246           for (x=0; x < width/2; x++) {
2247             val = (buf[c][y+1][x] << 4) / mul[c];
2248             if (val < 0) val = 0;
2249             if (c) RAW(row+y*2+c-1,x*2+2-c) = val;
2250             else   RAW(row+r*2+y,x*2+y) = val;
2251           }
2252         memcpy (buf[c][0]+!c, buf[c][2], sizeof buf[c][0]-2*!c);
2253       }
2254     }
2255     for (y=row; y < row+4; y++)
2256       for (x=0; x < width; x++)
2257         if ((x+y) & 1) {
2258           r = x ? x-1 : x+1;
2259           s = x+1 < width ? x+1 : x-1;
2260           val = (RAW(y,x)-2048)*2 + (RAW(y,r)+RAW(y,s))/2;
2261           if (val < 0) val = 0;
2262           RAW(y,x) = val;
2263         }
2264   }
2265   for (i=0; i < height*width; i++)
2266     raw_image[i] = curve[raw_image[i]];
2267   maximum = 0x3fff;
2268 }
2269
2270 #undef FORYX
2271 #undef PREDICTOR
2272
2273 #ifdef NO_JPEG
2274 void CLASS kodak_jpeg_load_raw() {}
2275 void CLASS lossy_dng_load_raw() {}
2276 #else
2277
2278 METHODDEF(boolean)
2279 fill_input_buffer (j_decompress_ptr cinfo)
2280 {
2281   static uchar jpeg_buffer[4096];
2282   size_t nbytes;
2283
2284   nbytes = fread (jpeg_buffer, 1, 4096, ifp);
2285   swabb(jpeg_buffer, jpeg_buffer, nbytes);
2286   cinfo->src->next_input_byte = jpeg_buffer;
2287   cinfo->src->bytes_in_buffer = nbytes;
2288   return TRUE;
2289 }
2290
2291 void CLASS kodak_jpeg_load_raw()
2292 {
2293   struct jpeg_decompress_struct cinfo;
2294   struct jpeg_error_mgr jerr;
2295   JSAMPARRAY buf;
2296   JSAMPLE (*pixel)[3];
2297   int row, col;
2298
2299   cinfo.err = jpeg_std_error (&jerr);
2300   jpeg_create_decompress (&cinfo);
2301   jpeg_stdio_src (&cinfo, ifp);
2302   cinfo.src->fill_input_buffer = fill_input_buffer;
2303   jpeg_read_header (&cinfo, TRUE);
2304   jpeg_start_decompress (&cinfo);
2305   if ((cinfo.output_width      != width  ) ||
2306       (cinfo.output_height*2   != height ) ||
2307       (cinfo.output_components != 3      )) {
2308     fprintf (stderr,_("%s: incorrect JPEG dimensions\n"), ifname);
2309     jpeg_destroy_decompress (&cinfo);
2310     longjmp (failure, 3);
2311   }
2312   buf = (*cinfo.mem->alloc_sarray)
2313                 ((j_common_ptr) &cinfo, JPOOL_IMAGE, width*3, 1);
2314
2315   while (cinfo.output_scanline < cinfo.output_height) {
2316     row = cinfo.output_scanline * 2;
2317     jpeg_read_scanlines (&cinfo, buf, 1);
2318     pixel = (JSAMPLE (*)[3]) buf[0];
2319     for (col=0; col < width; col+=2) {
2320       RAW(row+0,col+0) = pixel[col+0][1] << 1;
2321       RAW(row+1,col+1) = pixel[col+1][1] << 1;
2322       RAW(row+0,col+1) = pixel[col][0] + pixel[col+1][0];
2323       RAW(row+1,col+0) = pixel[col][2] + pixel[col+1][2];
2324     }
2325   }
2326   jpeg_finish_decompress (&cinfo);
2327   jpeg_destroy_decompress (&cinfo);
2328   maximum = 0xff << 1;
2329 }
2330
2331 void CLASS gamma_curve (double pwr, double ts, int mode, int imax);
2332
2333 void CLASS lossy_dng_load_raw()
2334 {
2335   struct jpeg_decompress_struct cinfo;
2336   struct jpeg_error_mgr jerr;
2337   JSAMPARRAY buf;
2338   JSAMPLE (*pixel)[3];
2339   unsigned sorder=order, ntags, opcode, deg, i, j, c;
2340   unsigned save=data_offset-4, trow=0, tcol=0, row, col;
2341   ushort cur[3][256];
2342   double coeff[9], tot;
2343
2344   if (meta_offset) {
2345     fseek (ifp, meta_offset, SEEK_SET);
2346     order = 0x4d4d;
2347     ntags = get4();
2348     while (ntags--) {
2349       opcode = get4(); get4(); get4();
2350       if (opcode != 8)
2351       { fseek (ifp, get4(), SEEK_CUR); continue; }
2352       fseek (ifp, 20, SEEK_CUR);
2353       if ((c = get4()) > 2) break;
2354       fseek (ifp, 12, SEEK_CUR);
2355       if ((deg = get4()) > 8) break;
2356       for (i=0; i <= deg && i < 9; i++)
2357         coeff[i] = getreal(12);
2358       for (i=0; i < 256; i++) {
2359         for (tot=j=0; j <= deg; j++)
2360           tot += coeff[j] * pow(i/255.0, j);
2361         cur[c][i] = tot*0xffff;
2362       }
2363     }
2364     order = sorder;
2365   } else {
2366     gamma_curve (1/2.4, 12.92, 1, 255);
2367     FORC3 memcpy (cur[c], curve, sizeof cur[0]);
2368   }
2369   cinfo.err = jpeg_std_error (&jerr);
2370   jpeg_create_decompress (&cinfo);
2371   while (trow < raw_height) {
2372     fseek (ifp, save+=4, SEEK_SET);
2373     if (tile_length < INT_MAX)
2374       fseek (ifp, get4(), SEEK_SET);
2375     jpeg_stdio_src (&cinfo, ifp);
2376     jpeg_read_header (&cinfo, TRUE);
2377     jpeg_start_decompress (&cinfo);
2378     buf = (*cinfo.mem->alloc_sarray)
2379         ((j_common_ptr) &cinfo, JPOOL_IMAGE, cinfo.output_width*3, 1);
2380     while (cinfo.output_scanline < cinfo.output_height &&
2381         (row = trow + cinfo.output_scanline) < height) {
2382       jpeg_read_scanlines (&cinfo, buf, 1);
2383       pixel = (JSAMPLE (*)[3]) buf[0];
2384       for (col=0; col < cinfo.output_width && tcol+col < width; col++) {
2385         FORC3 image[row*width+tcol+col][c] = cur[c][pixel[col][c]];
2386       }
2387     }
2388     jpeg_abort_decompress (&cinfo);
2389     if ((tcol += tile_width) >= raw_width)
2390       trow += tile_length + (tcol = 0);
2391   }
2392   jpeg_destroy_decompress (&cinfo);
2393   maximum = 0xffff;
2394 }
2395 #endif
2396
2397 void CLASS kodak_dc120_load_raw()
2398 {
2399   static const int mul[4] = { 162, 192, 187,  92 };
2400   static const int add[4] = {   0, 636, 424, 212 };
2401   uchar pixel[848];
2402   int row, shift, col;
2403
2404   for (row=0; row < height; row++) {
2405     if (fread (pixel, 1, 848, ifp) < 848) derror();
2406     shift = row * mul[row & 3] + add[row & 3];
2407     for (col=0; col < width; col++)
2408       RAW(row,col) = (ushort) pixel[(col + shift) % 848];
2409   }
2410   maximum = 0xff;
2411 }
2412
2413 void CLASS eight_bit_load_raw()
2414 {
2415   uchar *pixel;
2416   unsigned row, col;
2417
2418   pixel = (uchar *) calloc (raw_width, sizeof *pixel);
2419   merror (pixel, "eight_bit_load_raw()");
2420   for (row=0; row < raw_height; row++) {
2421     if (fread (pixel, 1, raw_width, ifp) < raw_width) derror();
2422     for (col=0; col < raw_width; col++)
2423       RAW(row,col) = curve[pixel[col]];
2424   }
2425   free (pixel);
2426   maximum = curve[0xff];
2427 }
2428
2429 void CLASS kodak_c330_load_raw()
2430 {
2431   uchar *pixel;
2432   int row, col, y, cb, cr, rgb[3], c;
2433
2434   pixel = (uchar *) calloc (raw_width, 2*sizeof *pixel);
2435   merror (pixel, "kodak_c330_load_raw()");
2436   for (row=0; row < height; row++) {
2437     if (fread (pixel, raw_width, 2, ifp) < 2) derror();
2438     if (load_flags && (row & 31) == 31)
2439       fseek (ifp, raw_width*32, SEEK_CUR);
2440     for (col=0; col < width; col++) {
2441       y  = pixel[col*2];
2442       cb = pixel[(col*2 & -4) | 1] - 128;
2443       cr = pixel[(col*2 & -4) | 3] - 128;
2444       rgb[1] = y - ((cb + cr + 2) >> 2);
2445       rgb[2] = rgb[1] + cb;
2446       rgb[0] = rgb[1] + cr;
2447       FORC3 image[row*width+col][c] = curve[LIM(rgb[c],0,255)];
2448     }
2449   }
2450   free (pixel);
2451   maximum = curve[0xff];
2452 }
2453
2454 void CLASS kodak_c603_load_raw()
2455 {
2456   uchar *pixel;
2457   int row, col, y, cb, cr, rgb[3], c;
2458
2459   pixel = (uchar *) calloc (raw_width, 3*sizeof *pixel);
2460   merror (pixel, "kodak_c603_load_raw()");
2461   for (row=0; row < height; row++) {
2462     if (~row & 1)
2463       if (fread (pixel, raw_width, 3, ifp) < 3) derror();
2464     for (col=0; col < width; col++) {
2465       y  = pixel[width*2*(row & 1) + col];
2466       cb = pixel[width + (col & -2)]   - 128;
2467       cr = pixel[width + (col & -2)+1] - 128;
2468       rgb[1] = y - ((cb + cr + 2) >> 2);
2469       rgb[2] = rgb[1] + cb;
2470       rgb[0] = rgb[1] + cr;
2471       FORC3 image[row*width+col][c] = curve[LIM(rgb[c],0,255)];
2472     }
2473   }
2474   free (pixel);
2475   maximum = curve[0xff];
2476 }
2477
2478 void CLASS kodak_262_load_raw()
2479 {
2480   static const uchar kodak_tree[2][26] =
2481   { { 0,1,5,1,1,2,0,0,0,0,0,0,0,0,0,0, 0,1,2,3,4,5,6,7,8,9 },
2482     { 0,3,1,1,1,1,1,2,0,0,0,0,0,0,0,0, 0,1,2,3,4,5,6,7,8,9 } };
2483   ushort *huff[2];
2484   uchar *pixel;
2485   int *strip, ns, c, row, col, chess, pi=0, pi1, pi2, pred, val;
2486
2487   FORC(2) huff[c] = make_decoder (kodak_tree[c]);
2488   ns = (raw_height+63) >> 5;
2489   pixel = (uchar *) malloc (raw_width*32 + ns*4);
2490   merror (pixel, "kodak_262_load_raw()");
2491   strip = (int *) (pixel + raw_width*32);
2492   order = 0x4d4d;
2493   FORC(ns) strip[c] = get4();
2494   for (row=0; row < raw_height; row++) {
2495     if ((row & 31) == 0) {
2496       fseek (ifp, strip[row >> 5], SEEK_SET);
2497       getbits(-1);
2498       pi = 0;
2499     }
2500     for (col=0; col < raw_width; col++) {
2501       chess = (row + col) & 1;
2502       pi1 = chess ? pi-2           : pi-raw_width-1;
2503       pi2 = chess ? pi-2*raw_width : pi-raw_width+1;
2504       if (col <= chess) pi1 = -1;
2505       if (pi1 < 0) pi1 = pi2;
2506       if (pi2 < 0) pi2 = pi1;
2507       if (pi1 < 0 && col > 1) pi1 = pi2 = pi-2;
2508       pred = (pi1 < 0) ? 0 : (pixel[pi1] + pixel[pi2]) >> 1;
2509       pixel[pi] = val = pred + ljpeg_diff (huff[chess]);
2510       if (val >> 8) derror();
2511       val = curve[pixel[pi++]];
2512       RAW(row,col) = val;
2513     }
2514   }
2515   free (pixel);
2516   FORC(2) free (huff[c]);
2517 }
2518
2519 int CLASS kodak_65000_decode (short *out, int bsize)
2520 {
2521   uchar c, blen[768];
2522   ushort raw[6];
2523   INT64 bitbuf=0;
2524   int save, bits=0, i, j, len, diff;
2525
2526   save = ftell(ifp);
2527   bsize = (bsize + 3) & -4;
2528   for (i=0; i < bsize; i+=2) {
2529     c = fgetc(ifp);
2530     if ((blen[i  ] = c & 15) > 12 ||
2531         (blen[i+1] = c >> 4) > 12 ) {
2532       fseek (ifp, save, SEEK_SET);
2533       for (i=0; i < bsize; i+=8) {
2534         read_shorts (raw, 6);
2535         out[i  ] = raw[0] >> 12 << 8 | raw[2] >> 12 << 4 | raw[4] >> 12;
2536         out[i+1] = raw[1] >> 12 << 8 | raw[3] >> 12 << 4 | raw[5] >> 12;
2537         for (j=0; j < 6; j++)
2538           out[i+2+j] = raw[j] & 0xfff;
2539       }
2540       return 1;
2541     }
2542   }
2543   if ((bsize & 7) == 4) {
2544     bitbuf  = fgetc(ifp) << 8;
2545     bitbuf += fgetc(ifp);
2546     bits = 16;
2547   }
2548   for (i=0; i < bsize; i++) {
2549     len = blen[i];
2550     if (bits < len) {
2551       for (j=0; j < 32; j+=8)
2552         bitbuf += (INT64) fgetc(ifp) << (bits+(j^8));
2553       bits += 32;
2554     }
2555     diff = bitbuf & (0xffff >> (16-len));
2556     bitbuf >>= len;
2557     bits -= len;
2558     if ((diff & (1 << (len-1))) == 0)
2559       diff -= (1 << len) - 1;
2560     out[i] = diff;
2561   }
2562   return 0;
2563 }
2564
2565 void CLASS kodak_65000_load_raw()
2566 {
2567   short buf[256];
2568   int row, col, len, pred[2], ret, i;
2569
2570   for (row=0; row < height; row++)
2571     for (col=0; col < width; col+=256) {
2572       pred[0] = pred[1] = 0;
2573       len = MIN (256, width-col);
2574       ret = kodak_65000_decode (buf, len);
2575       for (i=0; i < len; i++)
2576         if ((RAW(row,col+i) =   curve[ret ? buf[i] :
2577                 (pred[i & 1] += buf[i])]) >> 12) derror();
2578     }
2579 }
2580
2581 void CLASS kodak_ycbcr_load_raw()
2582 {
2583   short buf[384], *bp;
2584   int row, col, len, c, i, j, k, y[2][2], cb, cr, rgb[3];
2585   ushort *ip;
2586
2587   if (!image) return;
2588   for (row=0; row < height; row+=2)
2589     for (col=0; col < width; col+=128) {
2590       len = MIN (128, width-col);
2591       kodak_65000_decode (buf, len*3);
2592       y[0][1] = y[1][1] = cb = cr = 0;
2593       for (bp=buf, i=0; i < len; i+=2, bp+=2) {
2594         cb += bp[4];
2595         cr += bp[5];
2596         rgb[1] = -((cb + cr + 2) >> 2);
2597         rgb[2] = rgb[1] + cb;
2598         rgb[0] = rgb[1] + cr;
2599         for (j=0; j < 2; j++)
2600           for (k=0; k < 2; k++) {
2601             if ((y[j][k] = y[j][k^1] + *bp++) >> 10) derror();
2602             ip = image[(row+j)*width + col+i+k];
2603             FORC3 ip[c] = curve[LIM(y[j][k]+rgb[c], 0, 0xfff)];
2604           }
2605       }
2606     }
2607 }
2608
2609 void CLASS kodak_rgb_load_raw()
2610 {
2611   short buf[768], *bp;
2612   int row, col, len, c, i, rgb[3];
2613   ushort *ip=image[0];
2614
2615   for (row=0; row < height; row++)
2616     for (col=0; col < width; col+=256) {
2617       len = MIN (256, width-col);
2618       kodak_65000_decode (buf, len*3);
2619       memset (rgb, 0, sizeof rgb);
2620       for (bp=buf, i=0; i < len; i++, ip+=4)
2621         FORC3 if ((ip[c] = rgb[c] += *bp++) >> 12) derror();
2622     }
2623 }
2624
2625 void CLASS kodak_thumb_load_raw()
2626 {
2627   int row, col;
2628   colors = thumb_misc >> 5;
2629   for (row=0; row < height; row++)
2630     for (col=0; col < width; col++)
2631       read_shorts (image[row*width+col], colors);
2632   maximum = (1 << (thumb_misc & 31)) - 1;
2633 }
2634
2635 void CLASS sony_decrypt (unsigned *data, int len, int start, int key)
2636 {
2637   unsigned p = sony_p, *pad = sony_pad;;
2638   if (start) {
2639     for (p=0; p < 4; p++)
2640       pad[p] = key = key * 48828125 + 1;
2641     pad[3] = pad[3] << 1 | (pad[0]^pad[2]) >> 31;
2642     for (p=4; p < 127; p++)
2643       pad[p] = (pad[p-4]^pad[p-2]) << 1 | (pad[p-3]^pad[p-1]) >> 31;
2644     for (p=0; p < 127; p++)
2645       pad[p] = htonl(pad[p]);
2646   }
2647   while (len-- && p++)
2648     *data++ ^= pad[(p-1) & 127] = pad[p & 127] ^ pad[(p+64) & 127];
2649    sony_p = p;
2650 }
2651
2652 void CLASS sony_load_raw()
2653 {
2654   uchar head[40];
2655   ushort *pixel;
2656   unsigned i, key, row, col;
2657
2658   fseek (ifp, 200896, SEEK_SET);
2659   fseek (ifp, (unsigned) fgetc(ifp)*4 - 1, SEEK_CUR);
2660   order = 0x4d4d;
2661   key = get4();
2662   fseek (ifp, 164600, SEEK_SET);
2663   fread (head, 1, 40, ifp);
2664   sony_decrypt ((unsigned *) head, 10, 1, key);
2665   for (i=26; i-- > 22; )
2666     key = key << 8 | head[i];
2667   fseek (ifp, data_offset, SEEK_SET);
2668   for (row=0; row < raw_height; row++) {
2669     pixel = raw_image + row*raw_width;
2670     if (fread (pixel, 2, raw_width, ifp) < raw_width) derror();
2671     sony_decrypt ((unsigned *) pixel, raw_width/2, !row, key);
2672     for (col=0; col < raw_width; col++)
2673       if ((pixel[col] = ntohs(pixel[col])) >> 14) derror();
2674   }
2675   maximum = 0x3ff0;
2676 }
2677
2678 void CLASS sony_arw_load_raw()
2679 {
2680   ushort huff[32770];
2681   static const ushort tab[18] =
2682   { 0xf11,0xf10,0xe0f,0xd0e,0xc0d,0xb0c,0xa0b,0x90a,0x809,
2683     0x708,0x607,0x506,0x405,0x304,0x303,0x300,0x202,0x201 };
2684   int i, c, n, col, row, sum=0;
2685
2686   huff[0] = 15;
2687   for (n=i=0; i < 18; i++)
2688     FORC(32768 >> (tab[i] >> 8)) huff[++n] = tab[i];
2689   getbits(-1);
2690   for (col = raw_width; col--; )
2691     for (row=0; row < raw_height+1; row+=2) {
2692       if (row == raw_height) row = 1;
2693       if ((sum += ljpeg_diff(huff)) >> 12) derror();
2694       if (row < height) RAW(row,col) = sum;
2695     }
2696 }
2697
2698 void CLASS sony_arw2_load_raw()
2699 {
2700   uchar *data, *dp;
2701   ushort pix[16];
2702   int row, col, val, max, min, imax, imin, sh, bit, i;
2703
2704   data = (uchar *) malloc (raw_width+1);
2705   merror (data, "sony_arw2_load_raw()");
2706   for (row=0; row < height; row++) {
2707     fread (data, 1, raw_width, ifp);
2708     for (dp=data, col=0; col < raw_width-30; dp+=16) {
2709       max = 0x7ff & (val = sget4(dp));
2710       min = 0x7ff & val >> 11;
2711       imax = 0x0f & val >> 22;
2712       imin = 0x0f & val >> 26;
2713       for (sh=0; sh < 4 && 0x80 << sh <= max-min; sh++);
2714       for (bit=30, i=0; i < 16; i++)
2715         if      (i == imax) pix[i] = max;
2716         else if (i == imin) pix[i] = min;
2717         else {
2718           pix[i] = ((sget2(dp+(bit >> 3)) >> (bit & 7) & 0x7f) << sh) + min;
2719           if (pix[i] > 0x7ff) pix[i] = 0x7ff;
2720           bit += 7;
2721         }
2722       for (i=0; i < 16; i++, col+=2)
2723         RAW(row,col) = curve[pix[i] << 1] >> 2;
2724       col -= col & 1 ? 1:31;
2725     }
2726   }
2727   free (data);
2728 }
2729
2730 void CLASS samsung_load_raw()
2731 {
2732   int row, col, c, i, dir, op[4], len[4];
2733
2734   order = 0x4949;
2735   for (row=0; row < raw_height; row++) {
2736     fseek (ifp, strip_offset+row*4, SEEK_SET);
2737     fseek (ifp, data_offset+get4(), SEEK_SET);
2738     ph1_bits(-1);
2739     FORC4 len[c] = row < 2 ? 7:4;
2740     for (col=0; col < raw_width; col+=16) {
2741       dir = ph1_bits(1);
2742       FORC4 op[c] = ph1_bits(2);
2743       FORC4 switch (op[c]) {
2744         case 3: len[c] = ph1_bits(4);   break;
2745         case 2: len[c]--;               break;
2746         case 1: len[c]++;
2747       }
2748       for (c=0; c < 16; c+=2) {
2749         i = len[((c & 1) << 1) | (c >> 3)];
2750         RAW(row,col+c) = ((signed) ph1_bits(i) << (32-i) >> (32-i)) +
2751           (dir ? RAW(row+(~c | -2),col+c) : col ? RAW(row,col+(c | -2)) : 128);
2752         if (c == 14) c = -1;
2753       }
2754     }
2755   }
2756   for (row=0; row < raw_height-1; row+=2)
2757     for (col=0; col < raw_width-1; col+=2)
2758       SWAP (RAW(row,col+1), RAW(row+1,col));
2759 }
2760
2761 void CLASS samsung2_load_raw()
2762 {
2763   static const ushort tab[14] =
2764   { 0x304,0x307,0x206,0x205,0x403,0x600,0x709,
2765     0x80a,0x90b,0xa0c,0xa0d,0x501,0x408,0x402 };
2766   ushort huff[1026], vpred[2][2] = {{0,0},{0,0}}, hpred[2];
2767   int i, c, n, row, col, diff;
2768
2769   huff[0] = 10;
2770   for (n=i=0; i < 14; i++)
2771     FORC(1024 >> (tab[i] >> 8)) huff[++n] = tab[i];
2772   getbits(-1);
2773   for (row=0; row < raw_height; row++)
2774     for (col=0; col < raw_width; col++) {
2775       diff = ljpeg_diff (huff);
2776       if (col < 2) hpred[col] = vpred[row & 1][col] += diff;
2777       else         hpred[col & 1] += diff;
2778       RAW(row,col) = hpred[col & 1];
2779       if (hpred[col & 1] >> tiff_bps) derror();
2780     }
2781 }
2782
2783 void CLASS samsung3_load_raw()
2784 {
2785   int opt, init, mag, pmode, row, tab, col, pred, diff, i, c;
2786   ushort lent[3][2], len[4], *prow[2];
2787
2788   order = 0x4949;
2789   fseek (ifp, 9, SEEK_CUR);
2790   opt = fgetc(ifp);
2791   init = (get2(),get2());
2792   for (row=0; row < raw_height; row++) {
2793     fseek (ifp, (data_offset-ftell(ifp)) & 15, SEEK_CUR);
2794     ph1_bits(-1);
2795     mag = 0; pmode = 7;
2796     FORC(6) ((ushort *)lent)[c] = row < 2 ? 7:4;
2797     prow[ row & 1] = &RAW(row-1,1-((row & 1) << 1));    // green
2798     prow[~row & 1] = &RAW(row-2,0);                     // red and blue
2799     for (tab=0; tab+15 < raw_width; tab+=16) {
2800       if (~opt & 4 && !(tab & 63)) {
2801         i = ph1_bits(2);
2802         mag = i < 3 ? mag-'2'+"204"[i] : ph1_bits(12);
2803       }
2804       if (opt & 2)
2805         pmode = 7 - 4*ph1_bits(1);
2806       else if (!ph1_bits(1))
2807         pmode = ph1_bits(3);
2808       if (opt & 1 || !ph1_bits(1)) {
2809         FORC4 len[c] = ph1_bits(2);
2810         FORC4 {
2811           i = ((row & 1) << 1 | (c & 1)) % 3;
2812           len[c] = len[c] < 3 ? lent[i][0]-'1'+"120"[len[c]] : ph1_bits(4);
2813           lent[i][0] = lent[i][1];
2814           lent[i][1] = len[c];
2815         }
2816       }
2817       FORC(16) {
2818         col = tab + (((c & 7) << 1)^(c >> 3)^(row & 1));
2819         pred = (pmode == 7 || row < 2)
2820              ? (tab ? RAW(row,tab-2+(col & 1)) : init)
2821              : (prow[col & 1][col-'4'+"0224468"[pmode]] +
2822                 prow[col & 1][col-'4'+"0244668"[pmode]] + 1) >> 1;
2823         diff = ph1_bits (i = len[c >> 2]);
2824         if (diff >> (i-1)) diff -= 1 << i;
2825         diff = diff * (mag*2+1) + mag;
2826         RAW(row,col) = pred + diff;
2827       }
2828     }
2829   }
2830 }
2831
2832 #define HOLE(row) ((holes >> (((row) - raw_height) & 7)) & 1)
2833
2834 /* Kudos to Rich Taylor for figuring out SMaL's compression algorithm. */
2835 void CLASS smal_decode_segment (unsigned seg[2][2], int holes)
2836 {
2837   uchar hist[3][18] = {
2838     { 7, 7, 0, 0, 63, 55, 47, 39, 31, 23, 15, 7, 0 },
2839     { 7, 7, 0, 0, 63, 55, 47, 39, 31, 23, 15, 7, 0 },
2840     { 3, 3, 0, 0, 63,     47,     31,     15,    0 } };
2841   int low, high=0xff, carry=0, nbits=8;
2842   int pix, s, count, bin, next, i, sym[3];
2843   uchar diff, pred[]={0,0};
2844   ushort data=0, range=0;
2845
2846   fseek (ifp, seg[0][1]+1, SEEK_SET);
2847   getbits(-1);
2848   if (seg[1][0] > raw_width*raw_height)
2849       seg[1][0] = raw_width*raw_height;
2850   for (pix=seg[0][0]; pix < seg[1][0]; pix++) {
2851     for (s=0; s < 3; s++) {
2852       data = data << nbits | getbits(nbits);
2853       if (carry < 0)
2854         carry = (nbits += carry+1) < 1 ? nbits-1 : 0;
2855       while (--nbits >= 0)
2856         if ((data >> nbits & 0xff) == 0xff) break;
2857       if (nbits > 0)
2858           data = ((data & ((1 << (nbits-1)) - 1)) << 1) |
2859         ((data + (((data & (1 << (nbits-1)))) << 1)) & (-1 << nbits));
2860       if (nbits >= 0) {
2861         data += getbits(1);
2862         carry = nbits - 8;
2863       }
2864       count = ((((data-range+1) & 0xffff) << 2) - 1) / (high >> 4);
2865       for (bin=0; hist[s][bin+5] > count; bin++);
2866                 low = hist[s][bin+5] * (high >> 4) >> 2;
2867       if (bin) high = hist[s][bin+4] * (high >> 4) >> 2;
2868       high -= low;
2869       for (nbits=0; high << nbits < 128; nbits++);
2870       range = (range+low) << nbits;
2871       high <<= nbits;
2872       next = hist[s][1];
2873       if (++hist[s][2] > hist[s][3]) {
2874         next = (next+1) & hist[s][0];
2875         hist[s][3] = (hist[s][next+4] - hist[s][next+5]) >> 2;
2876         hist[s][2] = 1;
2877       }
2878       if (hist[s][hist[s][1]+4] - hist[s][hist[s][1]+5] > 1) {
2879         if (bin < hist[s][1])
2880           for (i=bin; i < hist[s][1]; i++) hist[s][i+5]--;
2881         else if (next <= bin)
2882           for (i=hist[s][1]; i < bin; i++) hist[s][i+5]++;
2883       }
2884       hist[s][1] = next;
2885       sym[s] = bin;
2886     }
2887     diff = sym[2] << 5 | sym[1] << 2 | (sym[0] & 3);
2888     if (sym[0] & 4)
2889       diff = diff ? -diff : 0x80;
2890     if (ftell(ifp) + 12 >= seg[1][1])
2891       diff = 0;
2892     raw_image[pix] = pred[pix & 1] += diff;
2893     if (!(pix & 1) && HOLE(pix / raw_width)) pix += 2;
2894   }
2895   maximum = 0xff;
2896 }
2897
2898 void CLASS smal_v6_load_raw()
2899 {
2900   unsigned seg[2][2];
2901
2902   fseek (ifp, 16, SEEK_SET);
2903   seg[0][0] = 0;
2904   seg[0][1] = get2();
2905   seg[1][0] = raw_width * raw_height;
2906   seg[1][1] = INT_MAX;
2907   smal_decode_segment (seg, 0);
2908 }
2909
2910 int CLASS median4 (int *p)
2911 {
2912   int min, max, sum, i;
2913
2914   min = max = sum = p[0];
2915   for (i=1; i < 4; i++) {
2916     sum += p[i];
2917     if (min > p[i]) min = p[i];
2918     if (max < p[i]) max = p[i];
2919   }
2920   return (sum - min - max) >> 1;
2921 }
2922
2923 void CLASS fill_holes (int holes)
2924 {
2925   int row, col, val[4];
2926
2927   for (row=2; row < height-2; row++) {
2928     if (!HOLE(row)) continue;
2929     for (col=1; col < width-1; col+=4) {
2930       val[0] = RAW(row-1,col-1);
2931       val[1] = RAW(row-1,col+1);
2932       val[2] = RAW(row+1,col-1);
2933       val[3] = RAW(row+1,col+1);
2934       RAW(row,col) = median4(val);
2935     }
2936     for (col=2; col < width-2; col+=4)
2937       if (HOLE(row-2) || HOLE(row+2))
2938         RAW(row,col) = (RAW(row,col-2) + RAW(row,col+2)) >> 1;
2939       else {
2940         val[0] = RAW(row,col-2);
2941         val[1] = RAW(row,col+2);
2942         val[2] = RAW(row-2,col);
2943         val[3] = RAW(row+2,col);
2944         RAW(row,col) = median4(val);
2945       }
2946   }
2947 }
2948
2949 void CLASS smal_v9_load_raw()
2950 {
2951   unsigned seg[256][2], offset, nseg, holes, i;
2952
2953   fseek (ifp, 67, SEEK_SET);
2954   offset = get4();
2955   nseg = (uchar) fgetc(ifp);
2956   fseek (ifp, offset, SEEK_SET);
2957   for (i=0; i < nseg*2; i++)
2958     ((unsigned *)seg)[i] = get4() + data_offset*(i & 1);
2959   fseek (ifp, 78, SEEK_SET);
2960   holes = fgetc(ifp);
2961   fseek (ifp, 88, SEEK_SET);
2962   seg[nseg][0] = raw_height * raw_width;
2963   seg[nseg][1] = get4() + data_offset;
2964   for (i=0; i < nseg; i++)
2965     smal_decode_segment (seg+i, holes);
2966   if (holes) fill_holes (holes);
2967 }
2968
2969 void CLASS redcine_load_raw()
2970 {
2971 #ifndef NO_JASPER
2972   int c, row, col;
2973   jas_stream_t *in;
2974   jas_image_t *jimg;
2975   jas_matrix_t *jmat;
2976   jas_seqent_t *data;
2977   ushort *img, *pix;
2978
2979   jas_init();
2980   in = jas_stream_fopen (ifname, "rb");
2981   jas_stream_seek (in, data_offset+20, SEEK_SET);
2982   jimg = jas_image_decode (in, -1, 0);
2983   if (!jimg) longjmp (failure, 3);
2984   jmat = jas_matrix_create (height/2, width/2);
2985   merror (jmat, "redcine_load_raw()");
2986   img = (ushort *) calloc ((height+2), (width+2)*2);
2987   merror (img, "redcine_load_raw()");
2988   FORC4 {
2989     jas_image_readcmpt (jimg, c, 0, 0, width/2, height/2, jmat);
2990     data = jas_matrix_getref (jmat, 0, 0);
2991     for (row = c >> 1; row < height; row+=2)
2992       for (col = c & 1; col < width; col+=2)
2993         img[(row+1)*(width+2)+col+1] = data[(row/2)*(width/2)+col/2];
2994   }
2995   for (col=1; col <= width; col++) {
2996     img[col] = img[2*(width+2)+col];
2997     img[(height+1)*(width+2)+col] = img[(height-1)*(width+2)+col];
2998   }
2999   for (row=0; row < height+2; row++) {
3000     img[row*(width+2)] = img[row*(width+2)+2];
3001     img[(row+1)*(width+2)-1] = img[(row+1)*(width+2)-3];
3002   }
3003   for (row=1; row <= height; row++) {
3004     pix = img + row*(width+2) + (col = 1 + (FC(row,1) & 1));
3005     for (   ; col <= width; col+=2, pix+=2) {
3006       c = (((pix[0] - 0x800) << 3) +
3007         pix[-(width+2)] + pix[width+2] + pix[-1] + pix[1]) >> 2;
3008       pix[0] = LIM(c,0,4095);
3009     }
3010   }
3011   for (row=0; row < height; row++)
3012     for (col=0; col < width; col++)
3013       RAW(row,col) = curve[img[(row+1)*(width+2)+col+1]];
3014   free (img);
3015   jas_matrix_destroy (jmat);
3016   jas_image_destroy (jimg);
3017   jas_stream_close (in);
3018 #endif
3019 }
3020
3021 /* RESTRICTED code starts here */
3022
3023 void CLASS foveon_decoder (unsigned size, unsigned code)
3024 {
3025   unsigned *huff = fov_huff;
3026   struct decode *cur;
3027   int i, len;
3028
3029   if (!code) {
3030     for (i=0; i < size; i++)
3031       huff[i] = get4();
3032     memset (first_decode, 0, sizeof first_decode);
3033     free_decode = first_decode;
3034   }
3035   cur = free_decode++;
3036   if (free_decode > first_decode+2048) {
3037     fprintf (stderr,_("%s: decoder table overflow\n"), ifname);
3038     longjmp (failure, 2);
3039   }
3040   if (code)
3041     for (i=0; i < size; i++)
3042       if (huff[i] == code) {
3043         cur->leaf = i;
3044         return;
3045       }
3046   if ((len = code >> 27) > 26) return;
3047   code = (len+1) << 27 | (code & 0x3ffffff) << 1;
3048
3049   cur->branch[0] = free_decode;
3050   foveon_decoder (size, code);
3051   cur->branch[1] = free_decode;
3052   foveon_decoder (size, code+1);
3053 }
3054
3055 void CLASS foveon_thumb()
3056 {
3057   unsigned bwide, row, col, bitbuf=0, bit=1, c, i;
3058   char *buf;
3059   struct decode *dindex;
3060   short pred[3];
3061
3062   bwide = get4();
3063   fprintf (ofp, "P6\n%d %d\n255\n", thumb_width, thumb_height);
3064   if (bwide > 0) {
3065     if (bwide < thumb_width*3) return;
3066     buf = (char *) malloc (bwide);
3067     merror (buf, "foveon_thumb()");
3068     for (row=0; row < thumb_height; row++) {
3069       fread  (buf, 1, bwide, ifp);
3070       fwrite (buf, 3, thumb_width, ofp);
3071     }
3072     free (buf);
3073     return;
3074   }
3075   foveon_decoder (256, 0);
3076
3077   for (row=0; row < thumb_height; row++) {
3078     memset (pred, 0, sizeof pred);
3079     if (!bit) get4();
3080     for (bit=col=0; col < thumb_width; col++)
3081       FORC3 {
3082         for (dindex=first_decode; dindex->branch[0]; ) {
3083           if ((bit = (bit-1) & 31) == 31)
3084             for (i=0; i < 4; i++)
3085               bitbuf = (bitbuf << 8) + fgetc(ifp);
3086           dindex = dindex->branch[bitbuf >> bit & 1];
3087         }
3088         pred[c] += dindex->leaf;
3089         fputc (pred[c], ofp);
3090       }
3091   }
3092 }
3093
3094 void CLASS foveon_sd_load_raw()
3095 {
3096   struct decode *dindex;
3097   short diff[1024];
3098   unsigned bitbuf=0;
3099   int pred[3], row, col, bit=-1, c, i;
3100
3101   read_shorts ((ushort *) diff, 1024);
3102   if (!load_flags) foveon_decoder (1024, 0);
3103
3104   for (row=0; row < height; row++) {
3105     memset (pred, 0, sizeof pred);
3106     if (!bit && !load_flags && atoi(model+2) < 14) get4();
3107     for (col=bit=0; col < width; col++) {
3108       if (load_flags) {
3109         bitbuf = get4();
3110         FORC3 pred[2-c] += diff[bitbuf >> c*10 & 0x3ff];
3111       }
3112       else FORC3 {
3113         for (dindex=first_decode; dindex->branch[0]; ) {
3114           if ((bit = (bit-1) & 31) == 31)
3115             for (i=0; i < 4; i++)
3116               bitbuf = (bitbuf << 8) + fgetc(ifp);
3117           dindex = dindex->branch[bitbuf >> bit & 1];
3118         }
3119         pred[c] += diff[dindex->leaf];
3120         if (pred[c] >> 16 && ~pred[c] >> 16) derror();
3121       }
3122       FORC3 image[row*width+col][c] = pred[c];
3123     }
3124   }
3125 }
3126
3127 void CLASS foveon_huff (ushort *huff)
3128 {
3129   int i, j, clen, code;
3130
3131   huff[0] = 8;
3132   for (i=0; i < 13; i++) {
3133     clen = getc(ifp);
3134     code = getc(ifp);
3135     for (j=0; j < 256 >> clen; )
3136       huff[code+ ++j] = clen << 8 | i;
3137   }
3138   get2();
3139 }
3140
3141 void CLASS foveon_dp_load_raw()
3142 {
3143   unsigned c, roff[4], row, col, diff;
3144   ushort huff[512], vpred[2][2], hpred[2];
3145
3146   fseek (ifp, 8, SEEK_CUR);
3147   foveon_huff (huff);
3148   roff[0] = 48;
3149   FORC3 roff[c+1] = -(-(roff[c] + get4()) & -16);
3150   FORC3 {
3151     fseek (ifp, data_offset+roff[c], SEEK_SET);
3152     getbits(-1);
3153     vpred[0][0] = vpred[0][1] = vpred[1][0] = vpred[1][1] = 512;
3154     for (row=0; row < height; row++) {
3155       for (col=0; col < width; col++) {
3156         diff = ljpeg_diff(huff);
3157         if (col < 2) hpred[col] = vpred[row & 1][col] += diff;
3158         else hpred[col & 1] += diff;
3159         image[row*width+col][c] = hpred[col & 1];
3160       }
3161     }
3162   }
3163 }
3164
3165 void CLASS foveon_load_camf()
3166 {
3167   unsigned type, wide, high, i, j, row, col, diff;
3168   ushort huff[258], vpred[2][2] = {{512,512},{512,512}}, hpred[2];
3169
3170   fseek (ifp, meta_offset, SEEK_SET);
3171   type = get4();  get4();  get4();
3172   wide = get4();
3173   high = get4();
3174   if (type == 2) {
3175     fread (meta_data, 1, meta_length, ifp);
3176     for (i=0; i < meta_length; i++) {
3177       high = (high * 1597 + 51749) % 244944;
3178       wide = high * (INT64) 301593171 >> 24;
3179       meta_data[i] ^= ((((high << 8) - wide) >> 1) + wide) >> 17;
3180     }
3181   } else if (type == 4) {
3182     free (meta_data);
3183     meta_data = (char *) malloc (meta_length = wide*high*3/2);
3184     merror (meta_data, "foveon_load_camf()");
3185     foveon_huff (huff);
3186     get4();
3187     getbits(-1);
3188     for (j=row=0; row < high; row++) {
3189       for (col=0; col < wide; col++) {
3190         diff = ljpeg_diff(huff);
3191         if (col < 2) hpred[col] = vpred[row & 1][col] += diff;
3192         else         hpred[col & 1] += diff;
3193         if (col & 1) {
3194           meta_data[j++] = hpred[0] >> 4;
3195           meta_data[j++] = hpred[0] << 4 | hpred[1] >> 8;
3196           meta_data[j++] = hpred[1];
3197         }
3198       }
3199     }
3200   } else
3201     fprintf (stderr,_("%s has unknown CAMF type %d.\n"), ifname, type);
3202 }
3203
3204 const char * CLASS foveon_camf_param (const char *block, const char *param)
3205 {
3206   unsigned idx, num;
3207   char *pos, *cp, *dp;
3208
3209   for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
3210     pos = meta_data + idx;
3211     if (strncmp (pos, "CMb", 3)) break;
3212     if (pos[3] != 'P') continue;
3213     if (strcmp (block, pos+sget4(pos+12))) continue;
3214     cp = pos + sget4(pos+16);
3215     num = sget4(cp);
3216     dp = pos + sget4(cp+4);
3217     while (num--) {
3218       cp += 8;
3219       if (!strcmp (param, dp+sget4(cp)))
3220         return dp+sget4(cp+4);
3221     }
3222   }
3223   return 0;
3224 }
3225
3226 void * CLASS foveon_camf_matrix (unsigned dim[3], const char *name)
3227 {
3228   unsigned i, idx, type, ndim, size, *mat;
3229   char *pos, *cp, *dp;
3230   double dsize;
3231
3232   for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
3233     pos = meta_data + idx;
3234     if (strncmp (pos, "CMb", 3)) break;
3235     if (pos[3] != 'M') continue;
3236     if (strcmp (name, pos+sget4(pos+12))) continue;
3237     dim[0] = dim[1] = dim[2] = 1;
3238     cp = pos + sget4(pos+16);
3239     type = sget4(cp);
3240     if ((ndim = sget4(cp+4)) > 3) break;
3241     dp = pos + sget4(cp+8);
3242     for (i=ndim; i--; ) {
3243       cp += 12;
3244       dim[i] = sget4(cp);
3245     }
3246     if ((dsize = (double) dim[0]*dim[1]*dim[2]) > meta_length/4) break;
3247     mat = (unsigned *) malloc ((size = dsize) * 4);
3248     merror (mat, "foveon_camf_matrix()");
3249     for (i=0; i < size; i++)
3250       if (type && type != 6)
3251         mat[i] = sget4(dp + i*4);
3252       else
3253         mat[i] = sget4(dp + i*2) & 0xffff;
3254     return mat;
3255   }
3256   fprintf (stderr,_("%s: \"%s\" matrix not found!\n"), ifname, name);
3257   return 0;
3258 }
3259
3260 int CLASS foveon_fixed (void *ptr, int size, const char *name)
3261 {
3262   void *dp;
3263   unsigned dim[3];
3264
3265   if (!name) return 0;
3266   dp = foveon_camf_matrix (dim, name);
3267   if (!dp) return 0;
3268   memcpy (ptr, dp, size*4);
3269   free (dp);
3270   return 1;
3271 }
3272
3273 float CLASS foveon_avg (short *pix, int range[2], float cfilt)
3274 {
3275   int i;
3276   float val, min=FLT_MAX, max=-FLT_MAX, sum=0;
3277
3278   for (i=range[0]; i <= range[1]; i++) {
3279     sum += val = pix[i*4] + (pix[i*4]-pix[(i-1)*4]) * cfilt;
3280     if (min > val) min = val;
3281     if (max < val) max = val;
3282   }
3283   if (range[1] - range[0] == 1) return sum/2;
3284   return (sum - min - max) / (range[1] - range[0] - 1);
3285 }
3286
3287 short * CLASS foveon_make_curve (double max, double mul, double filt)
3288 {
3289   short *curve;
3290   unsigned i, size;
3291   double x;
3292
3293   if (!filt) filt = 0.8;
3294   size = 4*M_PI*max / filt;
3295   if (size > INT_MAX-1) size = INT_MAX-1;
3296   curve = (short *) calloc (size+1, sizeof *curve);
3297   merror (curve, "foveon_make_curve()");
3298   curve[0] = size;
3299   for (i=0; i < size; i++) {
3300     x = i*filt/max/4;
3301     curve[i+1] = (cos(x)+1)/2 * tanh(i*filt/mul) * mul + 0.5;
3302   }
3303   return curve;
3304 }
3305
3306 void CLASS foveon_make_curves
3307         (short **curvep, float dq[3], float div[3], float filt)
3308 {
3309   double mul[3], max=0;
3310   int c;
3311
3312   FORC3 mul[c] = dq[c]/div[c];
3313   FORC3 if (max < mul[c]) max = mul[c];
3314   FORC3 curvep[c] = foveon_make_curve (max, mul[c], filt);
3315 }
3316
3317 int CLASS foveon_apply_curve (short *curve, int i)
3318 {
3319   if (abs(i) >= curve[0]) return 0;
3320   return i < 0 ? -curve[1-i] : curve[1+i];
3321 }
3322
3323 #define image ((short (*)[4]) image)
3324
3325 void CLASS foveon_interpolate()
3326 {
3327   static const short hood[] = { -1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1 };
3328   short *pix, prev[3], *curve[8], (*shrink)[3];
3329   float cfilt, ddft[3][3][2], ppm[3][3][3];
3330   float cam_xyz[3][3], correct[3][3], last[3][3], trans[3][3];
3331   float chroma_dq[3], color_dq[3], diag[3][3], div[3];
3332   float (*black)[3], (*sgain)[3], (*sgrow)[3];
3333   float fsum[3], val, frow, num;
3334   int row, col, c, i, j, diff, sgx, irow, sum, min, max, limit;
3335   int dscr[2][2], dstb[4], (*smrow[7])[3], total[4], ipix[3];
3336   int work[3][3], smlast, smred, smred_p, dev[3];
3337   int satlev[3], keep[4], active[4];
3338   unsigned dim[3], *badpix;
3339   double dsum, trsum[3];
3340   char str[128];
3341   const char* cp;
3342 // clear local storage
3343   pix = 0; ZERO(prev); ZERO(curve); ZERO(shrink);
3344   cfilt = 0; ZERO(ddft); ZERO(ppm);
3345   ZERO(cam_xyz); ZERO(correct); ZERO(last); ZERO(trans);
3346   ZERO(chroma_dq); ZERO(color_dq); ZERO(diag); ZERO(div);
3347   ZERO(black); ZERO(sgain); ZERO(sgrow);
3348   ZERO(fsum); val = frow = num = 0;
3349   row = col = c = i = j = diff = sgx = irow = sum = min = max = limit = 0;
3350   ZERO(dscr); ZERO(dstb); ZERO(smrow); ZERO(total); ZERO(ipix);
3351   ZERO(work); ZERO(smlast); ZERO(smred); smred_p=0; ZERO(dev);
3352   ZERO(satlev); ZERO(keep); ZERO(active);
3353   ZERO(dim); badpix = 0;
3354   dsum = 0; ZERO(trsum);
3355   ZERO(str);
3356   cp = 0;
3357
3358   if (verbose)
3359     fprintf (stderr,_("Foveon interpolation...\n"));
3360
3361   foveon_load_camf();
3362   foveon_fixed (dscr, 4, "DarkShieldColRange");
3363   foveon_fixed (ppm[0][0], 27, "PostPolyMatrix");
3364   foveon_fixed (satlev, 3, "SaturationLevel");
3365   foveon_fixed (keep, 4, "KeepImageArea");
3366   foveon_fixed (active, 4, "ActiveImageArea");
3367   foveon_fixed (chroma_dq, 3, "ChromaDQ");
3368   foveon_fixed (color_dq, 3,
3369         foveon_camf_param ("IncludeBlocks", "ColorDQ") ?
3370                 "ColorDQ" : "ColorDQCamRGB");
3371   if (foveon_camf_param ("IncludeBlocks", "ColumnFilter"))
3372                  foveon_fixed (&cfilt, 1, "ColumnFilter");
3373
3374   memset (ddft, 0, sizeof ddft);
3375   if (!foveon_camf_param ("IncludeBlocks", "DarkDrift")
3376          || !foveon_fixed (ddft[1][0], 12, "DarkDrift"))
3377     for (i=0; i < 2; i++) {
3378       foveon_fixed (dstb, 4, i ? "DarkShieldBottom":"DarkShieldTop");
3379       for (row = dstb[1]; row <= dstb[3]; row++)
3380         for (col = dstb[0]; col <= dstb[2]; col++)
3381           FORC3 ddft[i+1][c][1] += (short) image[row*width+col][c];
3382       FORC3 ddft[i+1][c][1] /= (dstb[3]-dstb[1]+1) * (dstb[2]-dstb[0]+1);
3383     }
3384
3385   if (!(cp = foveon_camf_param ("WhiteBalanceIlluminants", model2)))
3386   { fprintf (stderr,_("%s: Invalid white balance \"%s\"\n"), ifname, model2);
3387     return; }
3388   foveon_fixed (cam_xyz, 9, cp);
3389   foveon_fixed (correct, 9,
3390         foveon_camf_param ("WhiteBalanceCorrections", model2));
3391   memset (last, 0, sizeof last);
3392   for (i=0; i < 3; i++)
3393     for (j=0; j < 3; j++)
3394       FORC3 last[i][j] += correct[i][c] * cam_xyz[c][j];
3395
3396   #define LAST(x,y) last[(i+x)%3][(c+y)%3]
3397   for (i=0; i < 3; i++)
3398     FORC3 diag[c][i] = LAST(1,1)*LAST(2,2) - LAST(1,2)*LAST(2,1);
3399   #undef LAST
3400   FORC3 div[c] = diag[c][0]*0.3127 + diag[c][1]*0.329 + diag[c][2]*0.3583;
3401   sprintf (str, "%sRGBNeutral", model2);
3402   if (foveon_camf_param ("IncludeBlocks", str))
3403     foveon_fixed (div, 3, str);
3404   num = 0;
3405   FORC3 if (num < div[c]) num = div[c];
3406   FORC3 div[c] /= num;
3407
3408   memset (trans, 0, sizeof trans);
3409   for (i=0; i < 3; i++)
3410     for (j=0; j < 3; j++)
3411       FORC3 trans[i][j] += rgb_cam[i][c] * last[c][j] * div[j];
3412   FORC3 trsum[c] = trans[c][0] + trans[c][1] + trans[c][2];
3413   dsum = (6*trsum[0] + 11*trsum[1] + 3*trsum[2]) / 20;
3414   for (i=0; i < 3; i++)
3415     FORC3 last[i][c] = trans[i][c] * dsum / trsum[i];
3416   memset (trans, 0, sizeof trans);
3417   for (i=0; i < 3; i++)
3418     for (j=0; j < 3; j++)
3419       FORC3 trans[i][j] += (i==c ? 32 : -1) * last[c][j] / 30;
3420
3421   foveon_make_curves (curve, color_dq, div, cfilt);
3422   FORC3 chroma_dq[c] /= 3;
3423   foveon_make_curves (curve+3, chroma_dq, div, cfilt);
3424   FORC3 dsum += chroma_dq[c] / div[c];
3425   curve[6] = foveon_make_curve (dsum, dsum, cfilt);
3426   curve[7] = foveon_make_curve (dsum*2, dsum*2, cfilt);
3427
3428   sgain = (float (*)[3]) foveon_camf_matrix (dim, "SpatialGain");
3429   if (!sgain) return;
3430   sgrow = (float (*)[3]) calloc (dim[1], sizeof *sgrow);
3431   sgx = (width + dim[1]-2) / (dim[1]-1);
3432
3433   black = (float (*)[3]) calloc (height, sizeof *black);
3434   for (row=0; row < height; row++) {
3435     for (i=0; i < 6; i++)
3436       ((float *)ddft[0])[i] = ((float *)ddft[1])[i] +
3437         row / (height-1.0) * (((float *)ddft[2])[i] - ((float *)ddft[1])[i]);
3438     FORC3 black[row][c] =
3439         ( foveon_avg (image[row*width]+c, dscr[0], cfilt) +
3440           foveon_avg (image[row*width]+c, dscr[1], cfilt) * 3
3441           - ddft[0][c][0] ) / 4 - ddft[0][c][1];
3442   }
3443   memcpy (black, black+8, sizeof *black*8);
3444   memcpy (black+height-11, black+height-22, 11*sizeof *black);
3445   memcpy (last, black, sizeof last);
3446
3447   for (row=1; row < height-1; row++) {
3448     FORC3 if (last[1][c] > last[0][c]) {
3449         if (last[1][c] > last[2][c])
3450           black[row][c] = (last[0][c] > last[2][c]) ? last[0][c]:last[2][c];
3451       } else
3452         if (last[1][c] < last[2][c])
3453           black[row][c] = (last[0][c] < last[2][c]) ? last[0][c]:last[2][c];
3454     memmove (last, last+1, 2*sizeof last[0]);
3455     memcpy (last[2], black[row+1], sizeof last[2]);
3456   }
3457   FORC3 black[row][c] = (last[0][c] + last[1][c])/2;
3458   FORC3 black[0][c] = (black[1][c] + black[3][c])/2;
3459
3460   val = 1 - exp(-1/24.0);
3461   memcpy (fsum, black, sizeof fsum);
3462   for (row=1; row < height; row++)
3463     FORC3 fsum[c] += black[row][c] =