dcraw.c

Go to the documentation of this file.

/*
   dcraw.c -- Dave Coffin's raw photo decoder
   Copyright 1997-2004 by Dave Coffin, dcoffin a cybercom o net

   This is a portable ANSI C program to convert raw image files from
   any digital camera into PPM format.  TIFF and CIFF parsing are
   based upon public specifications, but no such documentation is
   available for the raw sensor data, so writing this program has
   been an immense effort.

   This code is freely licensed for all uses, commercial and
   otherwise.  Comments, questions, and encouragement are welcome.

   $Revision: 1.251 $
   $Date: 2004/10/29 01:52:38 $
 */

#define _GNU_SOURCE
#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <setjmp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
/*
   By defining NO_JPEG, you lose only the ability to
   decode compressed .KDC files from the Kodak DC120.
 */
#ifndef NO_JPEG
#include <jpeglib.h>
#endif

#ifdef WIN32
#include <winsock2.h>
#pragma comment(lib, "ws2_32.lib")
#define strcasecmp stricmp
typedef __int64 INT64;
#else
#include <unistd.h>
#include <netinet/in.h>
typedef long long INT64;
#endif

#ifdef LJPEG_DECODE
#error Please compile dcraw.c by itself.
#error Do not link it with ljpeg_decode.
#endif

#ifndef LONG_BIT
#define LONG_BIT (8 * sizeof (long))
#endif

#define ushort UshORt
typedef unsigned char uchar;
typedef unsigned short ushort;

#define ABS_MAX ((ushort) -1)
#define RGB_MAX ((ushort) -1)

/*
   All global variables are defined here, and all functions that
   access them are prefixed with "CLASS".  Note that a thread-safe
   C++ class cannot have non-const static local variables.
 */



// Cinelerra
char dcraw_info[1024];
float **dcraw_data;
int dcraw_alpha;


FILE *ifp;
short order;
char *ifname, make[64], model[64], model2[64];
time_t timestamp;
int data_offset, curve_offset, curve_length;
int tiff_data_compression, kodak_data_compression;
int raw_height, raw_width, top_margin, left_margin;
int height, width, colors, black, rgb_max;
int iheight, iwidth, shrink;
int is_canon, is_cmy, is_foveon, use_coeff, trim, flip, xmag, ymag;
int zero_after_ff;
unsigned filters;
ushort (*image)[4], white[8][8];
void (*load_raw)();
float gamma_val=0.6, bright=1.0, red_scale=1.0, blue_scale=1.0;
int four_color_rgb=0, document_mode=0, quick_interpolate=0;
int verbose=0, use_auto_wb=0, use_camera_wb=0, use_secondary=0;
float camera_red, camera_blue;
float pre_mul[4], coeff[3][4];
float k1=1.5, k2=0.5, juice=0.0;
int histogram[0x2000];
void write_ppm(FILE *);


void (*write_fun)(FILE *) = write_ppm;


jmp_buf failure;

float green_scale=1.0;
float saturation = 1.0;
float contrast = 1.0;
int autoexposure=0, use_pivot=0, use_neutral_wb=0;
int alternate_scale = 0;
int center_weight = 0;
int use_camera_black = 1;
int user_black = -1;
float pivot_value = 0.75, exposure_compensation=0.0;
unsigned pivot_point[4], pivot_base[4];
float white_point_fraction = 0.99;

struct decode {
  struct decode *branch[2];
  int leaf;
} first_decode[2048], *second_decode, *free_decode;

#define CLASS

/*
   In order to inline this calculation, I make the risky
   assumption that all filter patterns can be described
   by a repeating pattern of eight rows and two columns

   Return values are either 0/1/2/3 = G/M/C/Y or 0/1/2/3 = R/G1/B/G2
 */
#define FC(row,col) \
        (filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3)

#define BAYER(row,col) \
        image[((row) >> shrink)*iwidth + ((col) >> shrink)][FC(row,col)]

/*
   PowerShot 600 uses 0xe1e4e1e4:

          0 1 2 3 4 5
        0 G M G M G M
        1 C Y C Y C Y
        2 M G M G M G
        3 C Y C Y C Y

   PowerShot A5 uses 0x1e4e1e4e:

          0 1 2 3 4 5
        0 C Y C Y C Y
        1 G M G M G M
        2 C Y C Y C Y
        3 M G M G M G

   PowerShot A50 uses 0x1b4e4b1e:

          0 1 2 3 4 5
        0 C Y C Y C Y
        1 M G M G M G
        2 Y C Y C Y C
        3 G M G M G M
        4 C Y C Y C Y
        5 G M G M G M
        6 Y C Y C Y C
        7 M G M G M G

   PowerShot Pro70 uses 0x1e4b4e1b:

          0 1 2 3 4 5
        0 Y C Y C Y C
        1 M G M G M G
        2 C Y C Y C Y
        3 G M G M G M
        4 Y C Y C Y C
        5 G M G M G M
        6 C Y C Y C Y
        7 M G M G M G

   PowerShots Pro90 and G1 use 0xb4b4b4b4:

          0 1 2 3 4 5
        0 G M G M G M
        1 Y C Y C Y C

   All RGB cameras use one of these Bayer grids:

        0x16161616:     0x61616161:     0x49494949:     0x94949494:

          0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5     0 1 2 3 4 5
        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
        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
        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
        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

 */

#ifndef __GLIBC__
char *memmem (char *haystack, size_t haystacklen,
              char *needle, size_t needlelen)
{
  char *c;
  for (c = haystack; c <= haystack + haystacklen - needlelen; c++)
    if (!memcmp (c, needle, needlelen))
      return c;
  return NULL;
}
#endif

void CLASS merror (void *ptr, char *where)
{
  if (ptr) return;
  fprintf (stderr, "%s: Out of memory in %s\n", ifname, where);
  longjmp (failure, 1);
}

/*
   Get a 2-byte integer, making no assumptions about CPU byte order.
   Nor should we assume that the compiler evaluates left-to-right.
 */
ushort CLASS fget2 (FILE *f)
{
  uchar a, b;

  a = fgetc(f);
  b = fgetc(f);
  if (order == 0x4949)          /* "II" means little-endian */
    return a + (b << 8);
  else                          /* "MM" means big-endian */
    return (a << 8) + b;
}

/*
   Same for a 4-byte integer.
 */
int CLASS fget4 (FILE *f)
{
  uchar a, b, c, d;

  a = fgetc(f);
  b = fgetc(f);
  c = fgetc(f);
  d = fgetc(f);
  if (order == 0x4949)
    return a + (b << 8) + (c << 16) + (d << 24);
  else
    return (a << 24) + (b << 16) + (c << 8) + d;
}

void CLASS canon_600_load_raw()
{
  uchar  data[1120], *dp;
  ushort pixel[896], *pix;
  int irow, orow, col;

  for (irow=orow=0; irow < height; irow++)
  {
    fread (data, 1120, 1, ifp);
    for (dp=data, pix=pixel; dp < data+1120; dp+=10, pix+=8)
    {
      pix[0] = (dp[0] << 2) + (dp[1] >> 6    );
      pix[1] = (dp[2] << 2) + (dp[1] >> 4 & 3);
      pix[2] = (dp[3] << 2) + (dp[1] >> 2 & 3);
      pix[3] = (dp[4] << 2) + (dp[1]      & 3);
      pix[4] = (dp[5] << 2) + (dp[9]      & 3);
      pix[5] = (dp[6] << 2) + (dp[9] >> 2 & 3);
      pix[6] = (dp[7] << 2) + (dp[9] >> 4 & 3);
      pix[7] = (dp[8] << 2) + (dp[9] >> 6    );
    }
    for (col=0; col < width; col++)
      BAYER(orow,col) = pixel[col] << 4;
    for (col=width; col < 896; col++)
      black += pixel[col];
    if ((orow+=2) > height)
      orow = 1;
  }
  black = ((INT64) black << 4) / ((896 - width) * height);
}

void CLASS canon_a5_load_raw()
{
  uchar  data[1940], *dp;
  ushort pixel[1552], *pix;
  int row, col;

  for (row=0; row < height; row++) {
    fread (data, raw_width * 10 / 8, 1, ifp);
    for (dp=data, pix=pixel; pix < pixel+raw_width; dp+=10, pix+=8)
    {
      pix[0] = (dp[1] << 2) + (dp[0] >> 6);
      pix[1] = (dp[0] << 4) + (dp[3] >> 4);
      pix[2] = (dp[3] << 6) + (dp[2] >> 2);
      pix[3] = (dp[2] << 8) + (dp[5]     );
      pix[4] = (dp[4] << 2) + (dp[7] >> 6);
      pix[5] = (dp[7] << 4) + (dp[6] >> 4);
      pix[6] = (dp[6] << 6) + (dp[9] >> 2);
      pix[7] = (dp[9] << 8) + (dp[8]     );
    }
    for (col=0; col < width; col++)
      BAYER(row,col) = (pixel[col] & 0x3ff) << 4;
    for (col=width; col < raw_width; col++)
      black += pixel[col] & 0x3ff;
  }
  if (raw_width > width)
    black = ((INT64) black << 4) / ((raw_width - width) * height);
}

/*
   getbits(-1) initializes the buffer
   getbits(n) where 0 <= n <= 25 returns an n-bit integer
 */
unsigned CLASS getbits (int nbits)
{
  static unsigned long bitbuf=0;
  static int vbits=0;
  unsigned c, ret;

  if (nbits == 0) return 0;
  if (nbits == -1)
    ret = bitbuf = vbits = 0;
  else {
    ret = bitbuf << (LONG_BIT - vbits) >> (LONG_BIT - nbits);
    vbits -= nbits;
  }
  while (vbits < LONG_BIT - 7) {
    c = fgetc(ifp);
    bitbuf = (bitbuf << 8) + c;
    if (c == 0xff && zero_after_ff)
      fgetc(ifp);
    vbits += 8;
  }
  return ret;
}

void CLASS init_decoder ()
{
  memset (first_decode, 0, sizeof first_decode);
  free_decode = first_decode;
}

/*
   Construct a decode tree according the specification in *source.
   The first 16 bytes specify how many codes should be 1-bit, 2-bit
   3-bit, etc.  Bytes after that are the leaf values.

   For example, if the source is

    { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,
      0x04,0x03,0x05,0x06,0x02,0x07,0x01,0x08,0x09,0x00,0x0a,0x0b,0xff  },

   then the code is

        00              0x04
        010             0x03
        011             0x05
        100             0x06
        101             0x02
        1100            0x07
        1101            0x01
        11100           0x08
        11101           0x09
        11110           0x00
        111110          0x0a
        1111110         0x0b
        1111111         0xff
 */
uchar * CLASS make_decoder (const uchar *source, int level)
{
  struct decode *cur;
  static int leaf;
  int i, next;

  if (level==0) leaf=0;
  cur = free_decode++;
  if (free_decode > first_decode+2048) {
    fprintf (stderr, "%s: decoder table overflow\n", ifname);
    longjmp (failure, 2);
  }
  for (i=next=0; i <= leaf && next < 16; )
    i += source[next++];
  if (i > leaf) {
    if (level < next) {
      cur->branch[0] = free_decode;
      make_decoder (source, level+1);
      cur->branch[1] = free_decode;
      make_decoder (source, level+1);
    } else
      cur->leaf = source[16 + leaf++];
  }
  return (uchar *) source + 16 + leaf;
}

void CLASS crw_init_tables (unsigned table)
{
  static const uchar first_tree[3][29] = {
    { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,
      0x04,0x03,0x05,0x06,0x02,0x07,0x01,0x08,0x09,0x00,0x0a,0x0b,0xff  },

    { 0,2,2,3,1,1,1,1,2,0,0,0,0,0,0,0,
      0x03,0x02,0x04,0x01,0x05,0x00,0x06,0x07,0x09,0x08,0x0a,0x0b,0xff  },

    { 0,0,6,3,1,1,2,0,0,0,0,0,0,0,0,0,
      0x06,0x05,0x07,0x04,0x08,0x03,0x09,0x02,0x00,0x0a,0x01,0x0b,0xff  },
  };

  static const uchar second_tree[3][180] = {
    { 0,2,2,2,1,4,2,1,2,5,1,1,0,0,0,139,
      0x03,0x04,0x02,0x05,0x01,0x06,0x07,0x08,
      0x12,0x13,0x11,0x14,0x09,0x15,0x22,0x00,0x21,0x16,0x0a,0xf0,
      0x23,0x17,0x24,0x31,0x32,0x18,0x19,0x33,0x25,0x41,0x34,0x42,
      0x35,0x51,0x36,0x37,0x38,0x29,0x79,0x26,0x1a,0x39,0x56,0x57,
      0x28,0x27,0x52,0x55,0x58,0x43,0x76,0x59,0x77,0x54,0x61,0xf9,
      0x71,0x78,0x75,0x96,0x97,0x49,0xb7,0x53,0xd7,0x74,0xb6,0x98,
      0x47,0x48,0x95,0x69,0x99,0x91,0xfa,0xb8,0x68,0xb5,0xb9,0xd6,
      0xf7,0xd8,0x67,0x46,0x45,0x94,0x89,0xf8,0x81,0xd5,0xf6,0xb4,
      0x88,0xb1,0x2a,0x44,0x72,0xd9,0x87,0x66,0xd4,0xf5,0x3a,0xa7,
      0x73,0xa9,0xa8,0x86,0x62,0xc7,0x65,0xc8,0xc9,0xa1,0xf4,0xd1,
      0xe9,0x5a,0x92,0x85,0xa6,0xe7,0x93,0xe8,0xc1,0xc6,0x7a,0x64,
      0xe1,0x4a,0x6a,0xe6,0xb3,0xf1,0xd3,0xa5,0x8a,0xb2,0x9a,0xba,
      0x84,0xa4,0x63,0xe5,0xc5,0xf3,0xd2,0xc4,0x82,0xaa,0xda,0xe4,
      0xf2,0xca,0x83,0xa3,0xa2,0xc3,0xea,0xc2,0xe2,0xe3,0xff,0xff  },

    { 0,2,2,1,4,1,4,1,3,3,1,0,0,0,0,140,
      0x02,0x03,0x01,0x04,0x05,0x12,0x11,0x06,
      0x13,0x07,0x08,0x14,0x22,0x09,0x21,0x00,0x23,0x15,0x31,0x32,
      0x0a,0x16,0xf0,0x24,0x33,0x41,0x42,0x19,0x17,0x25,0x18,0x51,
      0x34,0x43,0x52,0x29,0x35,0x61,0x39,0x71,0x62,0x36,0x53,0x26,
      0x38,0x1a,0x37,0x81,0x27,0x91,0x79,0x55,0x45,0x28,0x72,0x59,
      0xa1,0xb1,0x44,0x69,0x54,0x58,0xd1,0xfa,0x57,0xe1,0xf1,0xb9,
      0x49,0x47,0x63,0x6a,0xf9,0x56,0x46,0xa8,0x2a,0x4a,0x78,0x99,
      0x3a,0x75,0x74,0x86,0x65,0xc1,0x76,0xb6,0x96,0xd6,0x89,0x85,
      0xc9,0xf5,0x95,0xb4,0xc7,0xf7,0x8a,0x97,0xb8,0x73,0xb7,0xd8,
      0xd9,0x87,0xa7,0x7a,0x48,0x82,0x84,0xea,0xf4,0xa6,0xc5,0x5a,
      0x94,0xa4,0xc6,0x92,0xc3,0x68,0xb5,0xc8,0xe4,0xe5,0xe6,0xe9,
      0xa2,0xa3,0xe3,0xc2,0x66,0x67,0x93,0xaa,0xd4,0xd5,0xe7,0xf8,
      0x88,0x9a,0xd7,0x77,0xc4,0x64,0xe2,0x98,0xa5,0xca,0xda,0xe8,
      0xf3,0xf6,0xa9,0xb2,0xb3,0xf2,0xd2,0x83,0xba,0xd3,0xff,0xff  },

    { 0,0,6,2,1,3,3,2,5,1,2,2,8,10,0,117,
      0x04,0x05,0x03,0x06,0x02,0x07,0x01,0x08,
      0x09,0x12,0x13,0x14,0x11,0x15,0x0a,0x16,0x17,0xf0,0x00,0x22,
      0x21,0x18,0x23,0x19,0x24,0x32,0x31,0x25,0x33,0x38,0x37,0x34,
      0x35,0x36,0x39,0x79,0x57,0x58,0x59,0x28,0x56,0x78,0x27,0x41,
      0x29,0x77,0x26,0x42,0x76,0x99,0x1a,0x55,0x98,0x97,0xf9,0x48,
      0x54,0x96,0x89,0x47,0xb7,0x49,0xfa,0x75,0x68,0xb6,0x67,0x69,
      0xb9,0xb8,0xd8,0x52,0xd7,0x88,0xb5,0x74,0x51,0x46,0xd9,0xf8,
      0x3a,0xd6,0x87,0x45,0x7a,0x95,0xd5,0xf6,0x86,0xb4,0xa9,0x94,
      0x53,0x2a,0xa8,0x43,0xf5,0xf7,0xd4,0x66,0xa7,0x5a,0x44,0x8a,
      0xc9,0xe8,0xc8,0xe7,0x9a,0x6a,0x73,0x4a,0x61,0xc7,0xf4,0xc6,
      0x65,0xe9,0x72,0xe6,0x71,0x91,0x93,0xa6,0xda,0x92,0x85,0x62,
      0xf3,0xc5,0xb2,0xa4,0x84,0xba,0x64,0xa5,0xb3,0xd2,0x81,0xe5,
      0xd3,0xaa,0xc4,0xca,0xf2,0xb1,0xe4,0xd1,0x83,0x63,0xea,0xc3,
      0xe2,0x82,0xf1,0xa3,0xc2,0xa1,0xc1,0xe3,0xa2,0xe1,0xff,0xff  }
  };

  if (table > 2) table = 2;
  init_decoder();
  make_decoder ( first_tree[table], 0);
  second_decode = free_decode;
  make_decoder (second_tree[table], 0);
}

/*
   Return 0 if the image starts with compressed data,
   1 if it starts with uncompressed low-order bits.

   In Canon compressed data, 0xff is always followed by 0x00.
 */
int CLASS canon_has_lowbits()
{
  uchar test[0x4000];
  int ret=1, i;

  fseek (ifp, 0, SEEK_SET);
  fread (test, 1, sizeof test, ifp);
  for (i=540; i < sizeof test - 1; i++)
    if (test[i] == 0xff) {
      if (test[i+1]) return 1;
      ret=0;
    }
  return ret;
}

void CLASS canon_compressed_load_raw()
{
  ushort *pixel, *prow;
  int lowbits, shift, i, row, r, col, save, val;
  unsigned irow, icol;
  struct decode *decode, *dindex;
  int block, diffbuf[64], leaf, len, diff, carry=0, pnum=0, base[2];
  uchar c;
  INT64 bblack=0;

  pixel = calloc (raw_width*8, sizeof *pixel);
  merror (pixel, "canon_compressed_load_raw()");
  lowbits = canon_has_lowbits();
  shift = 4 - lowbits*2;
  fseek (ifp, 540 + lowbits*raw_height*raw_width/4, SEEK_SET);
  zero_after_ff = 1;
  getbits(-1);

  for (row = 0; row < raw_height; row += 8) {
    for (block=0; block < raw_width >> 3; block++) {
      memset (diffbuf, 0, sizeof diffbuf);
      decode = first_decode;
      for (i=0; i < 64; i++ ) {
        for (dindex=decode; dindex->branch[0]; )
          dindex = dindex->branch[getbits(1)];
        leaf = dindex->leaf;
        decode = second_decode;
        if (leaf == 0 && i) break;
        if (leaf == 0xff) continue;
        i  += leaf >> 4;
        len = leaf & 15;
        if (len == 0) continue;
        diff = getbits(len);
        if ((diff & (1 << (len-1))) == 0)
          diff -= (1 << len) - 1;
        if (i < 64) diffbuf[i] = diff;
      }
      diffbuf[0] += carry;
      carry = diffbuf[0];
      for (i=0; i < 64; i++ ) {
        if (pnum++ % raw_width == 0)
          base[0] = base[1] = 512;
        pixel[(block << 6) + i] = ( base[i & 1] += diffbuf[i] );
      }
    }
    if (lowbits) {
      save = ftell(ifp);                        /* Don't lose our place */
      fseek (ifp, 26 + row*raw_width/4, SEEK_SET);
      for (prow=pixel, i=0; i < raw_width*2; i++) {
        c = fgetc(ifp);
        for (r=0; r < 8; r+=2, prow++) {
          val = (*prow << 2) + ((c >> r) & 3);
          if (raw_width == 2672 && val < 512) val += 2;
          *prow = val;
        }
      }
      fseek (ifp, save, SEEK_SET);
    }
    for (r=0; r < 8; r++) {
      irow = row - top_margin + r;
      if (irow >= height) continue;
      for (col = 0; col < raw_width; col++) {
        icol = col - left_margin;
        if (icol < width)
          BAYER(irow,icol) = pixel[r*raw_width+col] << shift;
        else
          bblack += pixel[r*raw_width+col];
      }
    }
  }
  free(pixel);
  if (raw_width > width)
    black = (bblack << shift) / ((raw_width - width) * height);
}

void CLASS kodak_curve (ushort *curve)
{
  int i, entries, tag, type, len, val;

  for (i=0; i < 0x1000; i++)
    curve[i] = i;
  if (strcasecmp(make,"KODAK")) return;
  if (!curve_offset) {
    fseek (ifp, 12, SEEK_SET);
    entries = fget2(ifp);
    while (entries--) {
      tag  = fget2(ifp);
      type = fget2(ifp);
      len  = fget4(ifp);
      val  = fget4(ifp);
      if (tag == 0x90d) {
        curve_offset = val;
        curve_length = len;
      }
    }
  }
  if (curve_offset) {
    fseek (ifp, curve_offset, SEEK_SET);
    for (i=0; i < curve_length; i++)
      curve[i] = fget2(ifp);
    for ( ; i < 0x1000; i++)
      curve[i] = curve[i-1];
    rgb_max = curve[i-1] << 2;
  }
  fseek (ifp, data_offset, SEEK_SET);
}

/*
   Not a full implementation of Lossless JPEG,
   just enough to decode Canon and Kodak images.
 */
void CLASS lossless_jpeg_load_raw()
{
  int tag, len, jhigh=0, jwide=0, jrow, jcol, jidx, diff, i, row, col;
  uchar data[256], *dp;
  int vpred[2] = { 0x800, 0x800 }, hpred[2];
  struct decode *dstart[2], *dindex;
  ushort curve[0x10000];
  INT64 bblack=0;
  int min=INT_MAX;

  kodak_curve(curve);
  order = 0x4d4d;
  if (fget2(ifp) != 0xffd8) return;

  do {
    tag = fget2(ifp);
    len = fget2(ifp) - 2;
    if (tag <= 0xff00 || len > 255) return;
    fread (data, 1, len, ifp);
    switch (tag) {
      case 0xffc3:
        jhigh = (data[1] << 8) + data[2];
        jwide =((data[3] << 8) + data[4])*2;
        break;
      case 0xffc4:
        init_decoder();
        dstart[0] = dstart[1] = free_decode;
        for (dp = data; dp < data+len && *dp < 2; ) {
          dstart[*dp] = free_decode;
          dp = make_decoder (++dp, 0);
        }
    }
  } while (tag != 0xffda);

  zero_after_ff = 1;
  getbits(-1);
  for (jrow=0; jrow < jhigh; jrow++)
  {
 
 
 
 
    for (jcol=0; jcol < jwide; jcol++)
    {
      for (dindex = dstart[jcol & 1]; dindex->branch[0]; )
        dindex = dindex->branch[getbits(1)];
      len = dindex->leaf;
      diff = getbits(len);
      if ((diff & (1 << (len-1))) == 0)
        diff -= (1 << len) - 1;
      if (jcol < 2) {
        vpred[jcol] += diff;
        hpred[jcol] = vpred[jcol];
      } else
        hpred[jcol & 1] += diff;
      diff = hpred[jcol & 1];
      if (diff < 0) diff = 0;
      jidx = jrow*jwide + jcol;
      if (raw_width == 5108) {
        i = jidx / (1680*jhigh);
        if (i < 2) {
          row = jidx / 1680 % jhigh;
          col = jidx % 1680 + i*1680;
        } else {
          jidx -= 2*1680*jhigh;
          row = jidx / 1748;
          col = jidx % 1748 + 2*1680;
        }
      } else {
        row = jidx / raw_width;
        col = jidx % raw_width;
      }
      if ((unsigned) (row-top_margin) >= height)
        continue;
      if ((unsigned) (col-left_margin) < width) {
        BAYER(row-top_margin,col-left_margin) = curve[diff] << 2;
        if (min > curve[diff])
            min = curve[diff];
      } else
        bblack += curve[diff];
    }
  }











  if (raw_width > width)
    black = (bblack << 2) / ((raw_width - width) * height);
  if (!strcasecmp(make,"KODAK"))
    black = min << 2;
}

void CLASS nikon_compressed_load_raw()
{
  static const uchar nikon_tree[] = {
    0,1,5,1,1,1,1,1,1,2,0,0,0,0,0,0,
    5,4,3,6,2,7,1,0,8,9,11,10,12
  };
  int vpred[4], hpred[2], csize, row, col, i, len, diff;
  ushort *curve;
  struct decode *dindex;

  init_decoder();
  make_decoder (nikon_tree, 0);

  fseek (ifp, curve_offset, SEEK_SET);
  for (i=0; i < 4; i++)
    vpred[i] = fget2(ifp);
  csize = fget2(ifp);
  curve = calloc (csize, sizeof *curve);
  merror (curve, "nikon_compressed_load_raw()");
  for (i=0; i < csize; i++)
    curve[i] = fget2(ifp);

  fseek (ifp, data_offset, SEEK_SET);
  getbits(-1);

  for (row=0; row < height; row++)
    for (col=0; col < raw_width; col++)
    {
      for (dindex=first_decode; dindex->branch[0]; )
        dindex = dindex->branch[getbits(1)];
      len = dindex->leaf;
      diff = getbits(len);
      if ((diff & (1 << (len-1))) == 0)
        diff -= (1 << len) - 1;
      if (col < 2) {
        i = 2*(row & 1) + (col & 1);
        vpred[i] += diff;
        hpred[col] = vpred[i];
      } else
        hpred[col & 1] += diff;
      if ((unsigned) (col-left_margin) >= width) continue;
      diff = hpred[col & 1];
      if (diff < 0) diff = 0;
      if (diff >= csize) diff = csize-1;
      BAYER(row,col-left_margin) = curve[diff] << 2;
    }
  free(curve);
}

void CLASS nikon_load_raw()
{
  int irow, row, col, i;

  getbits(-1);
  for (irow=0; irow < height; irow++) {
    row = irow;
    if (model[0] == 'E') {
      row = irow * 2 % height + irow / (height/2);
      if (row == 1 && atoi(model+1) < 5000) {
        fseek (ifp, 0, SEEK_END);
        fseek (ifp, ftell(ifp)/2, SEEK_SET);
        getbits(-1);
      }
    }
    for (col=0; col < raw_width; col++) {
      i = getbits(12);
      if ((unsigned) (col-left_margin) < width)
        BAYER(row,col-left_margin) = i << 2;
      if (tiff_data_compression == 34713 && (col % 10) == 9)
        getbits(8);
    }
  }
}

/*
   Figure out if a NEF file is compressed.  These fancy heuristics
   are only needed for the D100, thanks to a bug in some cameras
   that tags all images as "compressed".
 */
int CLASS nikon_is_compressed()
{
  uchar test[256];
  int i;

  if (tiff_data_compression != 34713)
    return 0;
  if (strcmp(model,"D100"))
    return 1;
  fseek (ifp, data_offset, SEEK_SET);
  fread (test, 1, 256, ifp);
  for (i=15; i < 256; i+=16)
    if (test[i]) return 1;
  return 0;
}

/*
   Returns 1 for a Coolpix 990, 0 for a Coolpix 995.
 */
int CLASS nikon_e990()
{
  int i, histo[256];
  const uchar often[] = { 0x00, 0x55, 0xaa, 0xff };

  memset (histo, 0, sizeof histo);
  fseek (ifp, 2064*1540*3/4, SEEK_SET);
  for (i=0; i < 2000; i++)
    histo[fgetc(ifp)]++;
  for (i=0; i < 4; i++)
    if (histo[often[i]] > 400)
      return 1;
  return 0;
}

/*
   Returns 1 for a Coolpix 2100, 0 for anything else.
 */
int CLASS nikon_e2100()
{
  uchar t[12];
  int i;

  fseek (ifp, 0, SEEK_SET);
  for (i=0; i < 1024; i++) {
    fread (t, 1, 12, ifp);
    if (((t[2] & t[4] & t[7] & t[9]) >> 4
        & t[1] & t[6] & t[8] & t[11] & 3) != 3)
      return 0;
  }
  return 1;
}

/*
   Separates a Minolta DiMAGE Z2 from a Nikon E4300.
 */
int CLASS minolta_z2()
{
  int i;
  char tail[424];

  fseek (ifp, -sizeof tail, SEEK_END);
  fread (tail, 1, sizeof tail, ifp);
  for (i=0; i < sizeof tail; i++)
    if (tail[i]) return 1;
  return 0;
}

void CLASS nikon_e2100_load_raw()
{
  uchar   data[3432], *dp;
  ushort pixel[2288], *pix;
  int row, col;

  for (row=0; row <= height; row+=2) {
    if (row == height) {
      fseek (ifp, width==1616 ? 8792:424, SEEK_CUR);
      row = 1;
    }
    fread (data, 1, width*3/2, ifp);
    for (dp=data, pix=pixel; pix < pixel+width; dp+=12, pix+=8) {
      pix[0] = (dp[2] >> 4) + (dp[ 3] << 4);
      pix[1] = (dp[2] << 8) +  dp[ 1];
      pix[2] = (dp[7] >> 4) + (dp[ 0] << 4);
      pix[3] = (dp[7] << 8) +  dp[ 6];
      pix[4] = (dp[4] >> 4) + (dp[ 5] << 4);
      pix[5] = (dp[4] << 8) +  dp[11];
      pix[6] = (dp[9] >> 4) + (dp[10] << 4);
      pix[7] = (dp[9] << 8) +  dp[ 8];
    }
    for (col=0; col < width; col++)
      BAYER(row,col) = (pixel[col] & 0xfff) << 2;
  }
}

void CLASS nikon_e950_load_raw()
{
  int irow, row, col;

  getbits(-1);
  for (irow=0; irow < height; irow++) {
    row = irow * 2 % height;
    for (col=0; col < width; col++)
      BAYER(row,col) = getbits(10) << 4;
    for (col=28; col--; )
      getbits(8);
  }
}

/*
   The Fuji Super CCD is just a Bayer grid rotated 45 degrees.
 */
void CLASS fuji_s2_load_raw()
{
  ushort pixel[2944];
  int row, col, r, c;

  fseek (ifp, (2944*24+32)*2, SEEK_CUR);
  for (row=0; row < 2144; row++) {
    fread (pixel, 2, 2944, ifp);
    for (col=0; col < 2880; col++) {
      r = row + ((col+1) >> 1);
      c = 2143 - row + (col >> 1);
      BAYER(r,c) = ntohs(pixel[col]) << 2;
    }
  }
}

void CLASS fuji_common_load_raw (int ncol, int icol, int nrow)
{
  ushort pixel[2048];
  int row, col, r, c;

  for (row=0; row < nrow; row++) {
    fread (pixel, 2, ncol, ifp);
    if (ntohs(0xaa55) == 0xaa55)        /* data is little-endian */
      swab (pixel, pixel, ncol*2);
    for (col=0; col <= icol; col++) {
      r = icol - col + (row >> 1);
      c = col + ((row+1) >> 1);
      BAYER(r,c) = pixel[col] << 2;
    }
  }
}

void CLASS fuji_s5000_load_raw()
{
  fseek (ifp, (1472*4+24)*2, SEEK_CUR);
  fuji_common_load_raw (1472, 1423, 2152);
}

void CLASS fuji_s7000_load_raw()
{
  fuji_common_load_raw (2048, 2047, 3080);
}

/*
   The Fuji Super CCD SR has two photodiodes for each pixel.
   The secondary has about 1/16 the sensitivity of the primary,
   but this ratio may vary.
 */
void CLASS fuji_f700_load_raw()
{
  ushort pixel[2944];
  int row, col, r, c, val;

  for (row=0; row < 2168; row++) {
    fread (pixel, 2, 2944, ifp);
    if (ntohs(0xaa55) == 0xaa55)        /* data is little-endian */
      swab (pixel, pixel, 2944*2);
    for (col=0; col < 1440; col++) {
      r = 1439 - col + (row >> 1);
      c = col + ((row+1) >> 1);
      val = pixel[col+16 + use_secondary*1472];
      BAYER(r,c) = val;
    }
  }
}

void CLASS rollei_load_raw()
{
  uchar pixel[10];
  unsigned iten=0, isix, i, buffer=0, row, col, todo[16];

  isix = raw_width * raw_height * 5 / 8;
  while (fread (pixel, 1, 10, ifp) == 10) {
    for (i=0; i < 10; i+=2) {
      todo[i]   = iten++;
      todo[i+1] = pixel[i] << 8 | pixel[i+1];
      buffer    = pixel[i] >> 2 | buffer << 6;
    }
    for (   ; i < 16; i+=2) {
      todo[i]   = isix++;
      todo[i+1] = buffer >> (14-i)*5;
    }
    for (i=0; i < 16; i+=2) {
      row = todo[i] / raw_width - top_margin;
      col = todo[i] % raw_width - left_margin;
      if (row < height && col < width)
        BAYER(row,col) = (todo[i+1] & 0x3ff) << 4;
    }
  }
}

void CLASS phase_one_load_raw()
{
  int row, col, a, b;
  ushort pixel[4134], akey, bkey;

  fseek (ifp, 8, SEEK_CUR);
  fseek (ifp, fget4(ifp) + 296, SEEK_CUR);
  akey = fget2(ifp);
  bkey = fget2(ifp);
  fseek (ifp, data_offset + 12 + top_margin*raw_width*2, SEEK_SET);
  for (row=0; row < height; row++) {
    fread (pixel, 2, raw_width, ifp);
    for (col=0; col < raw_width; col+=2) {
      a = ntohs(pixel[col+0]) ^ akey;
      b = ntohs(pixel[col+1]) ^ bkey;
      pixel[col+0] = (b & 0xaaaa) | (a & 0x5555);
      pixel[col+1] = (a & 0xaaaa) | (b & 0x5555);
    }
    for (col=0; col < width; col++)
      BAYER(row,col) = pixel[col+left_margin];
  }
}

void CLASS ixpress_load_raw()
{
  ushort pixel[4090];
  int row, col;

  fseek (ifp, 304 + 6*2*4090, SEEK_SET);
  for (row=height; --row >= 0; ) {
    fread (pixel, 2, 4090, ifp);
    if (ntohs(0xaa55) == 0xaa55)        /* data is little-endian */
      swab (pixel, pixel, 4090*2);
    for (col=0; col < width; col++)
      BAYER(row,col) = pixel[width-1-col];
  }
}

/* For this function only, raw_width is in bytes, not pixels! */
void CLASS packed_12_load_raw()
{
  int row, col;

  getbits(-1);
  for (row=0; row < height; row++) {
    for (col=0; col < width; col++)
      BAYER(row,col) = getbits(12) << 2;
    for (col = width*3/2; col < raw_width; col++)
      getbits(8);
  }
}

void CLASS unpacked_load_raw (int order, int rsh)
{
  ushort *pixel;
  int row, col;

  pixel = calloc (raw_width, sizeof *pixel);
  merror (pixel, "unpacked_load_raw()");
  for (row=0; row < height; row++) {
    fread (pixel, 2, raw_width, ifp);
    if (order != ntohs(0x55aa))
      swab (pixel, pixel, width*2);
    for (col=0; col < width; col++)
      BAYER(row,col) = pixel[col] << 8 >> (8+rsh);
  }
  free(pixel);
}

void CLASS be_16_load_raw()             /* "be" = "big-endian" */
{
  unpacked_load_raw (0x55aa, 0);
}

void CLASS be_high_12_load_raw()
{
  unpacked_load_raw (0x55aa, 2);
}

void CLASS be_low_12_load_raw()
{
  unpacked_load_raw (0x55aa,-2);
}

void CLASS be_low_10_load_raw()
{
  unpacked_load_raw (0x55aa,-4);
}

void CLASS le_high_12_load_raw()        /* "le" = "little-endian" */
{
  unpacked_load_raw (0xaa55, 2);
}

void CLASS olympus_cseries_load_raw()
{
  int irow, row, col;

  for (irow=0; irow < height; irow++) {
    row = irow * 2 % height + irow / (height/2);
    if (row < 2) {
      fseek (ifp, data_offset - row*(-width*height*3/4 & -2048), SEEK_SET);
      getbits(-1);
    }
    for (col=0; col < width; col++)
      BAYER(row,col) = getbits(12) << 2;
  }
}

void CLASS eight_bit_load_raw()
{
  uchar *pixel;
  int row, col;

  pixel = calloc (raw_width, sizeof *pixel);
  merror (pixel, "eight_bit_load_raw()");
  for (row=0; row < height; row++) {
    fread (pixel, 1, raw_width, ifp);
    for (col=0; col < width; col++)
      BAYER(row,col) = pixel[col] << 6;
  }
  free (pixel);
}

void CLASS casio_qv5700_load_raw()
{
  uchar  data[3232],  *dp;
  ushort pixel[2576], *pix;
  int row, col;

  for (row=0; row < height; row++) {
    fread (data, 1, 3232, ifp);
    for (dp=data, pix=pixel; dp < data+3220; dp+=5, pix+=4) {
      pix[0] = (dp[0] << 2) + (dp[1] >> 6);
      pix[1] = (dp[1] << 4) + (dp[2] >> 4);
      pix[2] = (dp[2] << 6) + (dp[3] >> 2);
      pix[3] = (dp[3] << 8) + (dp[4]     );
    }
    for (col=0; col < width; col++)
      BAYER(row,col) = (pixel[col] & 0x3ff) << 4;
  }
}

void CLASS nucore_load_raw()
{
  uchar *data, *dp;
  int irow, row, col;

  data = calloc (width, 2);
  merror (data, "nucore_load_raw()");
  for (irow=0; irow < height; irow++) {
    fread (data, 2, width, ifp);
    if (model[0] == 'B' && width == 2598)
      row = height - 1 - irow/2 - height/2 * (irow & 1);
    else
      row = irow;
    for (dp=data, col=0; col < width; col++, dp+=2)
      BAYER(row,col) = (dp[0] << 2) + (dp[1] << 10);
  }
  free(data);
}

const int * CLASS make_decoder_int (const int *source, int level)
{
  struct decode *cur;

  cur = free_decode++;
  if (level < source[0]) {
    cur->branch[0] = free_decode;
    source = make_decoder_int (source, level+1);
    cur->branch[1] = free_decode;
    source = make_decoder_int (source, level+1);
  } else {
    cur->leaf = source[1];
    source += 2;
  }
  return source;
}

int CLASS radc_token (int tree)
{
  int t;
  static struct decode *dstart[18], *dindex;
  static const int *s, source[] = {
    1,1, 2,3, 3,4, 4,2, 5,7, 6,5, 7,6, 7,8,
    1,0, 2,1, 3,3, 4,4, 5,2, 6,7, 7,6, 8,5, 8,8,
    2,1, 2,3, 3,0, 3,2, 3,4, 4,6, 5,5, 6,7, 6,8,
    2,0, 2,1, 2,3, 3,2, 4,4, 5,6, 6,7, 7,5, 7,8,
    2,1, 2,4, 3,0, 3,2, 3,3, 4,7, 5,5, 6,6, 6,8,
    2,3, 3,1, 3,2, 3,4, 3,5, 3,6, 4,7, 5,0, 5,8,
    2,3, 2,6, 3,0, 3,1, 4,4, 4,5, 4,7, 5,2, 5,8,
    2,4, 2,7, 3,3, 3,6, 4,1, 4,2, 4,5, 5,0, 5,8,
    2,6, 3,1, 3,3, 3,5, 3,7, 3,8, 4,0, 5,2, 5,4,
    2,0, 2,1, 3,2, 3,3, 4,4, 4,5, 5,6, 5,7, 4,8,
    1,0, 2,2, 2,-2,
    1,-3, 1,3,
    2,-17, 2,-5, 2,5, 2,17,
    2,-7, 2,2, 2,9, 2,18,
    2,-18, 2,-9, 2,-2, 2,7,
    2,-28, 2,28, 3,-49, 3,-9, 3,9, 4,49, 5,-79, 5,79,
    2,-1, 2,13, 2,26, 3,39, 4,-16, 5,55, 6,-37, 6,76,
    2,-26, 2,-13, 2,1, 3,-39, 4,16, 5,-55, 6,-76, 6,37
  };

  if (free_decode == first_decode)
    for (s=source, t=0; t < 18; t++) {
      dstart[t] = free_decode;
      s = make_decoder_int (s, 0);
    }
  if (tree == 18) {
    if (model[2] == '4')
      return (getbits(5) << 3) + 4;     /* DC40 */
    else
      return (getbits(6) << 2) + 2;     /* DC50 */
  }
  for (dindex = dstart[tree]; dindex->branch[0]; )
    dindex = dindex->branch[getbits(1)];
  return dindex->leaf;
}

#define FORYX for (y=1; y < 3; y++) for (x=col+1; x >= col; x--)

#define PREDICTOR (c ? (buf[c][y-1][x] + buf[c][y][x+1]) / 2 \
: (buf[c][y-1][x+1] + 2*buf[c][y-1][x] + buf[c][y][x+1]) / 4)

void CLASS kodak_radc_load_raw()
{
  int row, col, tree, nreps, rep, step, i, c, s, r, x, y, val;
  short last[3] = { 16,16,16 }, mul[3], buf[3][3][386];

  init_decoder();
  getbits(-1);
  for (i=0; i < sizeof(buf)/sizeof(short); i++)
    buf[0][0][i] = 2048;
  for (row=0; row < height; row+=4) {
    for (i=0; i < 3; i++)
      mul[i] = getbits(6);
    for (c=0; c < 3; c++) {
      val = ((0x1000000/last[c] + 0x7ff) >> 12) * mul[c];
      s = val > 65564 ? 10:12;
      x = ~(-1 << (s-1));
      val <<= 12-s;
      for (i=0; i < sizeof(buf[0])/sizeof(short); i++)
        buf[c][0][i] = (buf[c][0][i] * val + x) >> s;
      last[c] = mul[c];
      for (r=0; r <= !c; r++) {
        buf[c][1][width/2] = buf[c][2][width/2] = mul[c] << 7;
        for (tree=1, col=width/2; col > 0; ) {
          if ((tree = radc_token(tree))) {
            col -= 2;
            if (tree == 8)
              FORYX buf[c][y][x] = radc_token(tree+10) * mul[c];
            else
              FORYX buf[c][y][x] = radc_token(tree+10) * 16 + PREDICTOR;
          } else
            do {
              nreps = (col > 2) ? radc_token(9) + 1 : 1;
              for (rep=0; rep < 8 && rep < nreps && col > 0; rep++) {
                col -= 2;
                FORYX buf[c][y][x] = PREDICTOR;
                if (rep & 1) {
                  step = radc_token(10) << 4;
                  FORYX buf[c][y][x] += step;
                }
              }
            } while (nreps == 9);
        }
        for (y=0; y < 2; y++)
          for (x=0; x < width/2; x++) {
            val = (buf[c][y+1][x] << 4) / mul[c];
            if (val < 0) val = 0;
            if (c)
              BAYER(row+y*2+c-1,x*2+2-c) = val;
            else
              BAYER(row+r*2+y,x*2+y) = val;
          }
        memcpy (buf[c][0]+!c, buf[c][2], sizeof buf[c][0]-2*!c);
      }
    }
    for (y=row; y < row+4; y++)
      for (x=0; x < width; x++)
        if ((x+y) & 1) {
          val = (BAYER(y,x)-2048)*2 + (BAYER(y,x-1)+BAYER(y,x+1))/2;
          if (val < 0) val = 0;
          BAYER(y,x) = val;
        }
  }
}

#undef FORYX
#undef PREDICTOR

#ifdef NO_JPEG
void CLASS kodak_jpeg_load_raw() {}
#else

METHODDEF(boolean)
fill_input_buffer (j_decompress_ptr cinfo)
{
  static char jpeg_buffer[4096];
  size_t nbytes;

  nbytes = fread (jpeg_buffer, 1, 4096, ifp);
  swab (jpeg_buffer, jpeg_buffer, nbytes);
  cinfo->src->next_input_byte = jpeg_buffer;
  cinfo->src->bytes_in_buffer = nbytes;
  return TRUE;
}

void CLASS kodak_jpeg_load_raw()
{
  struct jpeg_decompress_struct cinfo;
  struct jpeg_error_mgr jerr;
  JSAMPARRAY buf;
  JSAMPLE (*pixel)[3];
  int row, col;

  cinfo.err = jpeg_std_error (&jerr);
  jpeg_create_decompress (&cinfo);
  jpeg_stdio_src (&cinfo, ifp);
  cinfo.src->fill_input_buffer = fill_input_buffer;
  jpeg_read_header (&cinfo, TRUE);
  jpeg_start_decompress (&cinfo);
  if ((cinfo.output_width      != width  ) ||
      (cinfo.output_height*2   != height ) ||
      (cinfo.output_components != 3      )) {
    fprintf (stderr, "%s: incorrect JPEG dimensions\n", ifname);
    jpeg_destroy_decompress (&cinfo);
    longjmp (failure, 3);
  }
  buf = (*cinfo.mem->alloc_sarray)
                ((j_common_ptr) &cinfo, JPOOL_IMAGE, width*3, 1);

  while (cinfo.output_scanline < cinfo.output_height) {
    row = cinfo.output_scanline * 2;
    jpeg_read_scanlines (&cinfo, buf, 1);
    pixel = (void *) buf[0];
    for (col=0; col < width; col+=2) {
      BAYER(row+0,col+0) = pixel[col+0][1] << 6;
      BAYER(row+1,col+1) = pixel[col+1][1] << 6;
      BAYER(row+0,col+1) = (pixel[col][0] + pixel[col+1][0]) << 5;
      BAYER(ro