PageRenderTime 798ms CodeModel.GetById 7ms app.highlight 699ms RepoModel.GetById 1ms app.codeStats 6ms

/src/FreeImage/Source/LibRawLite/dcraw/dcraw.c

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