/*
 * Program to convert .dat files from the Mitsubishi DJ-1000 digital
 * camera to PBM format.
 *
 * Interprets data in the .dat file (the CCD data), compensates the
 * color sensitivity, enhances color saturation, normalizes, sharpens,
 * adjusts gamma, and resamples the image to a 4x3 aspect ratio.
 *
 * Run with option -? to get a list of command line options.
 *
 * We assume that the pixels in the CCD array (apparently a Sanyo
 * LC9997M) represent this information in patterns of 2x2:
 *
 *   A C
 *   D E
 *
 * A = Cyan + White, C = Yellow + Green, D = Cyan + Green, E = Yellow + White
 *
 * Compute: Red = E - D, Blue = A - C, Green = (C + D)/2
 * And gains: GainRed = Green/Red, GainBlue = Green/Blue
 *
 * Some useful links:
 *   http://www.globaldialog.com/~biggers/html/body_dj-1000.html
 *   http://www.itojun.org/diary/19970730/
 *   http://www.thok.org/intranet/djcam/foekmail.txt
 *   ftp://ftp.ado.co.jp/pub/dj1000/
 *   http://www.phonk.net/Gedachten/dj-1000
 *
 * Author: Bert Bos <bert@w3.org>
 * Created: January 2000
 * Version: $Id: dj1000toppm.c,v 1.7 2018/12/11 13:25:35 bbos Exp $ */

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <math.h>
#include <errno.h>
#include <assert.h>

#define DEBUG 0			/* Write intermediate files tmpNNN.ppm */
#define ADJUST_SAT 1		/* Reduce saturation in shadows */

#define VERSION "$Revision: 1.7 $"

#define LINES 248				/* Useful lines in DAT file */
#define COLUMNS 512				/* Columns */
#define TOP 2					/* Lines to discard */
#define BOTTOM 4				/* Lines to discard */
#define LEFT 2					/* Columns to discard */
#define RIGHT 9					/* Columns to discard */
#define SCALEDLINES ((LINES - TOP - BOTTOM) * 3 / 2) /* After cropping and */
#define SCALEDCOLUMNS (COLUMNS - LEFT - RIGHT)	     /*   scaling */
enum {Red, Green, Blue} _RGB;			/* Index into image array */

typedef unsigned char byte;			/* [0,255] */

#define GAMMABITS 16
#define GAMMASIZE (1 << GAMMABITS)		/* Steps in gamma function */

#define HISTBITS 12
#define HISTSIZE (1 << HISTBITS)		/* # steps in histogram */

#define EPS 0.0001

static char *progname;				/* Name of program */
static int verbose_flag = 0;			/* Whether to print info */
static float white[3] = {2.3, 0.43, 3.1};	/* White balance daylight */
static float yellow[3] = {0.64, 1.50, 0.86};	/* White balance lightbulb */
static float saturation = 1.0;			/* Adjust color saturation */
static float sharpval = 0.5;			/* Adjust sharpness */
static float gammaval = 0.45;			/* Adjust gamma */
static float normalization = 3;			/* 3% normalization */
static float low = -1.0, high = -1.0;		/* By default not used */
static int npasses = 2;				/* Median filter (smoothing) */


/* syserr -- print s + system error message and exit */
static void syserr(const char *s) {perror(s); exit(1);}


/* abs -- return absolute value */
#define abs(a) ((a) < 0 ? -(a) : (a))


/* min3 -- return minimum of 3 values */
#define min3(a, b, c) ((a)<(b) ? ((a)<(c) ? (a) : (c)) : ((b)<(c) ? (b) : (c)))


/* max3 -- return maximum of 3 values */
#define max3(a, b, c) ((a)>(b) ? ((a)>(c) ? (a) : (c)) : ((b)>(c) ? (b) : (c)))


/* min -- return minimum of 2 values */
#define min(a, b) ((a) < (b) ? (a) : (b))


/* max -- return maximum of 2 values */
#define max(a, b) ((a) > (b) ? (a) : (b))


/* rgb2hsl -- convert RGB to HSL */
static void rgb2hsl(float r, float g, float b, float *hr, float *sr, float *lr)
{
  double h, s, l, max, min, del, rc, gc, bc;

  assert(0 <= r && r <= 1);
  assert(0 <= g && g <= 1);
  assert(0 <= b && b <= 1);

  max = max3(r, g, b);
  min = min3(r, g, b);

  del = max - min;
  l = (min + max)/2;
  s = l == 1 ? 0 : del / (1 - abs(2 * l - 1));

  h = -1;					/* No hue */
  if (s != 0.0) {
    rc = (max - r) / del;
    gc = (max - g) / del;
    bc = (max - b) / del;

    if (r == max) h = bc - gc;
    else if (g == max) h = 2 + rc - bc;
    else /* if (b == max) */ h = 4 + gc - rc;

    h = h * 60;
    if (h < 0) h += 360;
  }

  *hr = h;  *sr = s;  *lr = l;
}


/* hsl2rgb -- convert HSL to RGB */
static void hsl2rgb(float h, float s, float l, float *r, float *g, float *b)
{
  float chroma, x, m;

  if (h >= 360.0) h -= 360.0;
  h = h / 60.0;
  chroma = (1 - abs(2 * l - 1)) * s;
  x = chroma * (1 - abs(fmodf(h, 2) - 1));
  if (h <= 1)      {*r = chroma; *g = x; *b = 0;} /* Between red and yellow */
  else if (h <= 2) {*r = x; *g = chroma; *b = 0;} /* Between yellow and green */
  else if (h <= 3) {*r = 0; *g = chroma; *b = x;} /* Between green and cyan */
  else if (h <= 4) {*r = 0; *g = x; *b = chroma;} /* Between cyan and blue */
  else if (h <= 5) {*r = x; *g = 0; *b = chroma;} /* Between blue and magenta */
  else             {*r = chroma; *g = 0; *b = x;} /* Between magenta and red */
  m = l - chroma/2;
  *r += m; *g += m; *b += m;
}


#if DEBUG
/* dump_gray -- dump an intermediate gray file (for debugging) */
static void dump_gray(float gray[LINES][COLUMNS], int seqno)
{
  int i, j;
  char s[20];
  FILE *tmpfile;

  sprintf(s, "tmp%03d.pgm", seqno);
  tmpfile = fopen(s, "w");
  fprintf(tmpfile, "P5 %d %d 255\n", COLUMNS, LINES);
  for (i = 0; i < LINES; i++)
    for (j = 0; j < COLUMNS; j++)
      putc(gray[i][j] * 255 + 0.5, tmpfile);
  fclose(tmpfile);
}


/* dump -- dump an intermediate image file (for debugging) */
static void dump(float image[LINES][COLUMNS][3], int seqno)
{
  int i, j, k;
  char s[20];
  FILE *tmpfile;

  sprintf(s, "tmp%03d.ppm", seqno);
  tmpfile = fopen(s, "w");
  fprintf(tmpfile, "P6 %d %d 255\n", COLUMNS, LINES);
  for (i = 0; i < LINES; i++)
    for (j = 0; j < COLUMNS; j++) {
      for (k = 0; k < 3; k++) {
	int c = image[i][j][k] * 255 + 0.5;
	putc(c < 0 ? 0 : c > 255 ? 255 : c, tmpfile);
      }
    }
  fclose(tmpfile);
}


/* dumphsl -- dump an intermediate HSl image file (for debugging) */
static void dumphsl(float hsl[LINES][COLUMNS][3], int seqno)
{
  int i, j, c;
  char s[20];
  FILE *tmpfile;
  float r, g, b;

  sprintf(s, "tmp%03d.ppm", seqno);
  tmpfile = fopen(s, "w");
  fprintf(tmpfile, "P6 %d %d 255\n", COLUMNS, LINES);
  for (i = 0; i < LINES; i++)
    for (j = 0; j < COLUMNS; j++) {
      hsl2rgb(hsl[i][j][0], hsl[i][j][1], hsl[i][j][2], &r, &g, &b);
      c = r * 255 + 0.5;
      putc(c < 0 ? 0 : c > 255 ? 255 : c, tmpfile);
      c = g * 255 + 0.5;
      putc(c < 0 ? 0 : c > 255 ? 255 : c, tmpfile);
      c = b * 255 + 0.5;
      putc(c < 0 ? 0 : c > 255 ? 255 : c, tmpfile);
    }
  fclose(tmpfile);
}
#endif /* DEBUG */


/* fill_gamma_table -- fill a table with precomputed gamma values */
static void fill_gamma_table(float table[GAMMASIZE], float gamma)
{
  unsigned int i;

  for (i = 0; i < GAMMASIZE; i++)
    table[i] = pow((double)i / GAMMASIZE, gamma);
}


/* interpolate_luminosity -- compute luminosity by adding horizontal pairs */
static void interpolate_luminosity(byte ccd[LINES][COLUMNS],
				   float gray[LINES][COLUMNS])
{
  /*
   * To get the luminosity of a pixel, take the average over a
   * complete ACDE quadruple, as follows. The luminosity of the center
   * pixel in this arrangement
   *
   *   ... A C A ...
   *   ....D E D ...
   *   ... A C A ...
   *
   * is the average of: the E itelf, the average of the two D's, the average of
   * the two C's and the average of the four A's.
   */
  unsigned int i, j;

  for (i = 1; i < LINES - 1; i++)
    for (j = 1; j < COLUMNS - 1; j++)
      gray[i][j] = (ccd[i][j] +			      /* Center */
		    (ccd[i][j-1] + ccd[i][j+1])/2.0 + /* Sides */
		    (ccd[i-1][j] + ccd[i+1][j])/2.0 + /* Above and below */
		    (ccd[i-1][j-1] + ccd[i-1][j+1] +
		     ccd[i+1][j-1] + ccd[i+1][j+1])/4.0) / 255.0 / 4;

  /* First and last columns set to zeros */
  for (i = 0; i < LINES; i++) gray[i][COLUMNS-1] = gray[i][0] = 0;

  /* First and last rows set to zeros */
  for (j = 0; j < COLUMNS; j++) gray[0][j] = gray[LINES-1][j] = 0;

#if DEBUG
  dump_gray(gray, 0);
#endif
}


/* set_pixel - set one pixel in the image array */
static void set_pixel(float pixel[3], float r, float g, float b)
{
  float gain_r = g > EPS ? r/g : g < -EPS ? r/g : 1;
  float gain_b = g > EPS ? b/g : g < -EPS ? r/g : 1;
  r *= gain_r;
  b *= gain_b;
  pixel[Red] = min(1.0, max(0.0, r));
  pixel[Green] = min(1.0, max(0.0, g));
  pixel[Blue] = min(1.0, max(0.0, b));
}


/* interpolate_color -- get colors from the differences between pixels */
static void interpolate_color(byte ccd[LINES][COLUMNS],
			      float image[LINES][COLUMNS][3])
{
  /*
   * We assume that the pixels in the CCD array represent this information
   * in patterns of 2x2:
   *
   * A C
   * D E
   *
   * A = Cyan + White, C = Yellow + Green, D = Cyan + Green, E = Yellow + White
   * A = R + 2G + 2B,  C = R + 2G,         D = 2G + B,       E = 2R + 2G + B
   *
   * Compute: Red = E - D, Blue = A - C, Green = (C + D)/2
   * And gains: GainRed = Green/Red, GainBlue = Green/Blue
   */

  unsigned int i, j;

  for (i = 1; i < LINES - 2; i++)
    for (j = 1; j < COLUMNS - 2; j++) {
      float A, C, D, E;
      /*
       *   od  ev
       * D  E  D  E
       * A  C  A  C  odd
       * D  E  D  E  even
       * A  C  A  C
       */
      if ((i & 1) == 0) {	/* Even row */
	if ((j & 1) == 0) {	/* Even column */
	  D = ccd[i][j];
	  E = (ccd[i][j-1] + ccd[i][j+1])/2.0;
	  A = (ccd[i-1][j] + ccd[i+1][j])/2.0;
	  C = (ccd[i-1][j-1] + ccd[i-1][j+1] + ccd[i+1][j-1] + ccd[i+1][j+1])/4;
	} else {		/* Odd column */
	  D = (ccd[i][j-1] + ccd[i][j+1])/2.0;
	  E = ccd[i][j];
	  A = (ccd[i-1][j-1] + ccd[i-1][j+1] + ccd[i+1][j-1] + ccd[i+1][j+1])/4;
	  C = (ccd[i-1][j] + ccd[i+1][j])/2.0;
	}
      } else {			/* Odd row */
	if ((j & 1) == 0) {	/* Even column */
	  D = (ccd[i-1][j] + ccd[i+1][j])/2.0;
	  E = (ccd[i-1][j-1] + ccd[i-1][j+1] + ccd[i+1][j-1] + ccd[i+1][j+1])/4;
	  A = ccd[i][j];
	  C = (ccd[i][j-1] + ccd[i][j+1])/2.0;
	} else {		/* Odd column */
	  D = (ccd[i-1][j-1] + ccd[i-1][j+1] + ccd[i+1][j-1] + ccd[i+1][j+1])/4;
	  E = (ccd[i-1][j] + ccd[i+1][j])/2.0;
	  A = (ccd[i][j-1] + ccd[i][j+1])/2.0;
	  C = ccd[i][j];
	}
      }
      D /= 255.0; E /= 255.0; A /= 255.0; C /= 255.0;
      set_pixel(image[i][j], E - D, (C + D)/2, A - C);
    }
  for (j = 0; j < COLUMNS; j++) /* Top row */
    set_pixel(image[0][j], 0, 0, 0);
  for (j = 0; j < COLUMNS; j++) /* Bottom row */
    set_pixel(image[LINES-1][j], 0, 0, 0);
  for (i = 0; i < LINES; i++)	/* Left column */
    set_pixel(image[i][0], 0, 0, 0);
  for (i = 0; i < LINES; i++)	/* Right column */
    set_pixel(image[i][COLUMNS-1], 0, 0, 0);

#if DEBUG
  dump(image, 1);
#endif /* DEBUG */
}


/* smooth -- apply a 3x3 median filter */
static void smooth(float image[LINES][COLUMNS][3], int npasses)
{
  float image2[LINES][COLUMNS][3];
  int i, j, k, n;

  for (n = 0; n < npasses; n++) {

    /* Smooth into image2 */
    for (i = 1; i < LINES - 1; i++)
      for (j = 1; j < COLUMNS - 1; j++)
	for (k = 0; k < 3; k++)
	  image2[i][j][k] =
	    (image[i-1][j-1][k] + image[i-1][j][k] + image[i-1][j+1][k] +
	     image[i][j-1][k]   + image[i][j][k]   + image[i][j+1][k]   +
	     image[i+1][j-1][k] + image[i+1][j][k] + image[i+1][j+1][k])/9;

    /* Copy back to image */
    for (i = 1; i < LINES - 1; i++)
      for (j = 1; j < COLUMNS - 1; j++)
	for (k = 0; k < 3; k++)
	  image[i][j][k] = image2[i][j][k];
  }

  if (verbose_flag)
    fprintf(stderr, "Median passes:           %-d\n", npasses);

#if DEBUG
  dump(image, 2);
#endif /* DEBUG */
}


/* sum_color -- sum the red, green or blue of all pixels */
static float sum_color(int color, float image[LINES][COLUMNS][3],
		       int row1, int row2, int col1, int col2)
{
  if (col1 == col2) {
    assert(row1 == row2);
    return image[row1][col1][color];
  } else if (row1 == row2) {
    int m = (col1 + col2)/2;
    return sum_color(color, image, row1, row1, col1, m) +
      sum_color(color, image, row1, row1, m + 1, col2); /* Recurse */
  } else {
    int m = (row1 + row2)/2;
    return sum_color(color, image, row1, m, col1, col2) +
      sum_color(color, image, m + 1, row2, col1, col2); /* Recurse */
  }
}

/* balance_white -- do some magic to balance the colors */
static void balance_white(float image[LINES][COLUMNS][3], float factor[3])
{
  unsigned int i, j, k;
  float sum[3] = {0, 0, 0};
  float m, avg;

  /* See if we should use heuristics to determine the white balance */
  if (factor[0] < 0 || factor[1] < 0 || factor[2] < 0) {

    /* Compute energy of each color */
    sum[Red] = sum_color(Red, image, TOP, LINES - BOTTOM - 1,
			 LEFT, COLUMNS - RIGHT - 1);
    sum[Green] = sum_color(Green, image, TOP, LINES - BOTTOM - 1,
			   LEFT, COLUMNS - RIGHT - 1);
    sum[Blue] = sum_color(Blue, image, TOP, LINES - BOTTOM - 1,
			  LEFT, COLUMNS - RIGHT - 1);
    avg = (sum[Red] + sum[Green] + sum[Blue])/3;

    /* Calculate factors so that energy in the three colors is the same */
    for (k = 0; k < 3; k++) factor[k] = avg / sum[k];

  } else {

    /* Normalize so that factor[0] + factor[1] + factor[2] = 3.0 */
    m = factor[0] + factor[1] + factor[2];
    for (k = 0; k < 3; k++) factor[k] *= 3.0 / m;
  }

  /* Multiply all values by the given factor */
  for (i = 1; i < LINES; i++)
    for (j = 1; j < COLUMNS; j++)
      for (k = 0; k < 3; k++)
	image[i][j][k] = min(factor[k] * image[i][j][k], 1.0);

  if (verbose_flag)
    fprintf(stderr,
	    "Color adjustment:        %-5.3g   %5.3g   %5.3g\n",
	    factor[0], factor[1], factor[2]);

#if DEBUG
  dump(image, 3);
#endif /* DEBUG */
}


/* convert2hsl -- convert RGB image to HSL image */
static void convert2hsl(float rgb[LINES][COLUMNS][3],
			float gray[LINES][COLUMNS],
			float hsl[LINES][COLUMNS][3])
{
  int i, j;

  for (i = 1; i < LINES; i++)
    for (j = 1; j < COLUMNS; j++) {
      rgb2hsl(rgb[i][j][0], rgb[i][j][1], rgb[i][j][2],
	      &hsl[i][j][0], &hsl[i][j][1], &hsl[i][j][2]);

      /* Replace the lightness by the previously computed gray
	 level. This makes the result sharper makes highlights white
	 instead of green. */
      hsl[i][j][2] = gray[i][j];
}

#if DEBUG
  dumphsl(hsl, 4);
#endif
}


/* adjust_saturation -- multiply each pixel's saturation */
static void adjust_saturation(float hsl[LINES][COLUMNS][3], float saturation)
{
  unsigned int i, j;

  if (saturation != 1.0) {
    for (i = 1; i < LINES; i++)
      for (j = 1; j < COLUMNS; j++) {
	hsl[i][j][1] *= saturation;
	if (hsl[i][j][1] > 1.0) hsl[i][j][1] = 1.0;
      }
  }

  if (verbose_flag)
    fprintf(stderr, "Saturation adjustment:   %-5.3g\n", saturation);

#if DEBUG
  dumphsl(hsl, 5);
#endif
}


/* adjust_range -- ignore lowest and highest norm percent of pixel values */
static void adjust_range(float hsl[LINES][COLUMNS][3],
			 float norm, float lo, float hi)
{
  unsigned int histogram[HISTSIZE];
  unsigned int i, j;
  unsigned int s, ignored;
  int h;

  /* If n is given, but not both of lo and hi, derive missing lo/hi from n. */
  if ((lo < 0 || hi < 0) && norm >= 0) {

   ignored = norm * (LINES - TOP - BOTTOM) * (COLUMNS - LEFT - RIGHT)/100 + 0.5;

   /* Make a histogram of values. */
   for (i = 0; i < HISTSIZE; i++) histogram[i] = 0;

    for (i = TOP; i < LINES - BOTTOM; i++)
      for (j = LEFT; j < COLUMNS - RIGHT; j++) {
	h = hsl[i][j][2] * (HISTSIZE - 1);
	histogram[h]++;
      }
  
    if (lo < 0) {
      /* Find value lo such that norm percent of the pixels is below it */
      for (s = 0, i = 0; i < HISTSIZE && s <= ignored; i++) s += histogram[i];
      lo = (float)i / HISTSIZE;
    }

    if (hi < 0) {
      /* Find value hi such that norm percent of the pixels is above it */
      for (s = 0, i = HISTSIZE-1; i > 0 && s <= ignored; i--) s+= histogram[i];
      hi = (float)i / HISTSIZE;
    }
  }

  /* If we have a lo and a hi, use them to normalize the values. */
  if (lo >= 0 && hi >= 0) {

    if (lo > hi) { float h = lo; lo = hi; hi = h; }

    /* Redistribute all pixels' values such that lo maps to 0 and hi to 1 */
    for (i = 0; i < LINES; i++) {
      for (j = 0; j < COLUMNS; j++) {
	if (hsl[i][j][2] <= lo) hsl[i][j][2] = 0.0;
	else if (hsl[i][j][2] >= hi) hsl[i][j][2] = 1.0;
	else hsl[i][j][2] = (hsl[i][j][2] - lo)/(hi - lo);
#if ADJUST_SAT
	/* Reduce saturation in shadows */
	hsl[i][j][1] *= sqrt(sqrt(hsl[i][j][2]));
#endif
      }
    }
  }

  if (verbose_flag && lo >= 0 && hi >= 0)
    fprintf(stderr, "Normalize range:         %-5.3g - %-5.3g\n", lo, hi);
  else if (verbose_flag)
    fprintf(stderr, "Normalize range:         none\n");

#if DEBUG
  dumphsl(hsl, 6);
#endif
}


/* sharpen -- sharpen the image by a factor between 0.0 and 1.0 */
static void sharpen(float hsl[LINES][COLUMNS][3], float sharpval)
{
  unsigned int i, j;
  float val[LINES][COLUMNS], avg;

  if (sharpval > 0) {

    /* Make new values in val. (Leave outer pixels unchanged.) */
    for (i = 1; i < LINES - 1; i++)
      for (j = 1; j < COLUMNS - 1; j++) {
	avg = (hsl[i-1][j-1][2] + hsl[i-1][j][2] + hsl[i-1][j+1][2]
	       + hsl[i][j-1][2] + hsl[i][j][2] + hsl[i][j+1][2]
	       + hsl[i+1][j-1][2] + hsl[i+1][j][2] + hsl[i+1][j+1][2])/9;
	val[i][j] = (hsl[i][j][2] - sharpval * avg)/(1.0 - sharpval);
	if (val[i][j] < 0.0) val[i][j] = 0.0;
	else if (val[i][j] > 1.0) val[i][j] = 1.0;
      }

    /* Copy back to hsl */
    for (i = 1; i < LINES - 1; i++)
      for (j = 1; j < COLUMNS - 1; j++)
	hsl[i][j][2] = val[i][j];
  }

  if (verbose_flag)
    fprintf(stderr, "Sharpened:               %-5.3g\n", sharpval);

#if DEBUG
  dumphsl(hsl, 7);
#endif
}


/* adjust_gamma -- apply gamma function */
static void adjust_gamma(float hsl[LINES][COLUMNS][3], float gammaval)
{
  unsigned int i, j;
  float gamma[GAMMASIZE];

  if (gammaval != 1.0) {
    fill_gamma_table(gamma, gammaval);
    for (i = 0; i < LINES; i++)
      for (j = 0; j < COLUMNS; j++)
	hsl[i][j][2] = gamma[(int)(hsl[i][j][2] * (GAMMASIZE - 1) + 0.5)];
  }

  if (verbose_flag)
    fprintf(stderr, "Gamma adjustment:        %-5.3g\n", gammaval);

#if DEBUG
  dumphsl(hsl, 8);
#endif
}


/* convert2rgb -- convert HSL image to RGB image */
static void convert2rgb(float hsl[LINES][COLUMNS][3],
			float rgb[LINES][COLUMNS][3])
{
  int i, j;

  for (i = 1; i < LINES; i++)
    for (j = 1; j < COLUMNS; j++)
      hsl2rgb(hsl[i][j][0], hsl[i][j][1], hsl[i][j][2],
	      &rgb[i][j][0], &rgb[i][j][1], &rgb[i][j][2]);
}


/* adjust_size -- crop edges and adjust for the 3x2 aspect ratio of the pixels */
static void adjust_size(float image[LINES][COLUMNS][3],
			float image2[SCALEDLINES][SCALEDCOLUMNS][3])
{
  unsigned int i, j, k;

  for (i = 0; i < SCALEDLINES; i++)
    for (j = 0; j < SCALEDCOLUMNS; j++) {
      float q = (i + 0.5) * 2.0 / 3.0;
      int h = q;
      float p = q - h;		/* 0 <= p < 1 */
      for (k = 0; k < 3; k++)
	image2[i][j][k] = (1.0 - p) * image[h+TOP][j+LEFT][k] +
	  p * image[h+TOP+1][j+LEFT][k];
    }
}


/* quantize -- reduce to 3 x 8 bits/pixel */
static void quantize(float image[SCALEDLINES][SCALEDCOLUMNS][3],
		     byte ppm[SCALEDLINES][SCALEDCOLUMNS][3])
{
  unsigned int i, j, k;

  for (i = 0; i < SCALEDLINES; i++)
    for (j = 0; j < SCALEDCOLUMNS; j++)
      for (k = 0; k < 3; k++)
	ppm[i][j][k] = image[i][j][k] * 255 + 0.5;
}


/* dj1000toppm -- read data in DJ-1000 format and write it in PPM format */
static void dj1000toppm(char *inname, int in, char *outname, int out)
{
  byte ccd[LINES][COLUMNS];			/* 8 bits/pixel */
  float gray[LINES][COLUMNS];			/* Horizontal interpolation */
  float image[LINES][COLUMNS][3];		/* Working copy */
  float hsl[LINES][COLUMNS][3];			/* Working copy */
  float image2[SCALEDLINES][SCALEDCOLUMNS][3];	/* Cropped, with square pixels */
  byte ppm[SCALEDLINES][SCALEDCOLUMNS][3], *p;	/* 3 x 8 bits/pixel */
  ssize_t result;
  char s[30], *q;
  size_t n;

  /* Read the image data; no checks, just assume it is DJ-1000 file... */
  if (read(in, ccd, LINES * COLUMNS) < 0) syserr(inname);

  /* Extract colors, adjust gamma, saturation, etc. */
  interpolate_luminosity(ccd, gray);
  interpolate_color(ccd, image);
  smooth(image, npasses);
  balance_white(image, white);
  convert2hsl(image, gray, hsl);
  adjust_saturation(hsl, saturation);
  adjust_range(hsl, normalization, low, high);
  sharpen(hsl, sharpval);
  adjust_gamma(hsl, gammaval);
  convert2rgb(hsl, image);
  adjust_size(image, image2);
  quantize(image2, ppm);

  /* Write as a raw PPM file */
  n = sprintf(s, "P6 %d %d 255\n", SCALEDCOLUMNS, SCALEDLINES);
  q = s;
  while (n != 0) {
    result = write(out, q, n);
    if (result < n && errno != 0) syserr(outname);
    q += result;
    n -= result;
  }
  n = 3 * SCALEDLINES * SCALEDCOLUMNS;
  p = &ppm[0][0][0];
  while (n != 0) {
    result = write(out, p, n);
    if (result < n && errno != 0) syserr(outname);
    p += result;
    n -= result;
  }
}


/* usage -- print usage message and exit with an error */
static void usage(void)
{
  fprintf(stderr, "Usage: %s [options] [dj1000file [ppmfile]]\n", progname);
  fprintf(stderr, "Options:\n");
  fprintf(stderr, "  -V         print \"%s\" and exit\n", VERSION);
  fprintf(stderr, "  -v         verbose, print info about the calculations\n");
  fprintf(stderr, "  -?         this help text\n");
  fprintf(stderr, "  -i         indoor (artificial light, equiv. to -r %g -g %g -b %g)\n", yellow[Red], yellow[Green], yellow[Blue]);
  fprintf(stderr, "  -r number  relative multiplier for red, <0 = use heuristic [%g]\n", white[Red]);
  fprintf(stderr, "  -g number  relative multiplier for green, <0 = use heuristic [%g]\n", white[Green]);
  fprintf(stderr, "  -b number  relative multiplier for blue, <0 = use heuristic [%g]\n", white[Blue]);
  fprintf(stderr, "  -s number  saturation [%f]\n", saturation);
  fprintf(stderr, "  -S number  sharpen (0.0-1.0) [%g]\n", sharpval);
  fprintf(stderr, "  -G number  gamma correction [%g]\n", gammaval);
  fprintf(stderr, "  -n number  normalize, discard lightest & darkest percent [%g]\n", normalization);
  fprintf(stderr, "  -l number  discard values lower than this, <0 = not used\n");
  fprintf(stderr, "  -h number  discard values higher than this, <0 = not used\n");
  fprintf(stderr, "  -m number  number of passes of median filter [%d]\n", npasses);
  exit(2);
}


int main(int argc, char *argv[])
{
  int k, c, infile, outfile;
  char *inname, *outname;

  progname = argv[0];
  while ((c = getopt(argc, argv, "Vvir:g:b:s:S:G:n:l:h:m:")) != -1) {
    switch (c) {
    case 'V': printf("%s\n", VERSION); exit(0); break;
    case 'v': verbose_flag = 1; break;
    case 'i': for (k = 0; k < 3; k++) white[k] = yellow[k]; break;
    case 'r': white[Red] = atof(optarg); break;
    case 'g': white[Green] = atof(optarg); break;
    case 'b': white[Blue] = atof(optarg); break;
    case 's': saturation = atof(optarg); break;
    case 'S': sharpval = atof(optarg); break;
    case 'G': gammaval = atof(optarg); break;
    case 'n': normalization = atof(optarg); break;
    case 'l': low = atof(optarg); break;
    case 'h': high = atof(optarg); break;
    case 'm': npasses = atoi(optarg); break;
    default: usage();
    }
  }

  /* Check for some errors. */
  if (sharpval < 0 || sharpval >= 1) {
    fprintf(stderr, "Error: sharpening (-S) must be >= 0 and < 1\n");
    exit(1);
  }
  if (gammaval <= 0) {
    fprintf(stderr, "Error: gamma (-G) must be > 0\n");
    exit(1);
  }

  if (optind + 2 < argc) usage();
  if (optind >= argc) {
    inname = "<stdin>";
    infile = 0;
  } else {
    inname = argv[optind];
    if ((infile = open(inname, O_RDONLY)) < 0) syserr(inname);
  }
  optind++;
  if (optind >= argc) {
    outname = "<stdout>";
    outfile = 1;
  } else {
    outname = argv[optind];
    if ((outfile = creat(outname, 0666)) < 0) syserr(outname);
  }

  dj1000toppm(inname, infile, outname, outfile);

  if (close(infile) < 0) syserr(inname);
  if (close(outfile) < 0) syserr(outname);
  return 0;
}