/* Copyright (c) 1998 Regents of the University of California */

/*
 *  pflare.c - program to blur, vignette, and tint pictures
 *
 *     3/1998 - 3/1999
 *     Paul Debevec
 *
 *  Templated off of psum.c by Greg Ward
 *
 */

#include  <stdio.h>
#include  <malloc.h>
#include  <math.h>
#include  <string.h>

#include  "color.h"

#include "resolu.h"

/* Function headers */

void xshift(int dx);
void yshift(int dy);
void dzoom();
void giblur0(int twosigmasquared);
void giblur0x(int twosigmasquared);
void giblur0xy(int twosigmasquared);
void giblur0yx(int twosigmasquared);
void giblur0y(int twosigmasquared);
void sepia(float cr,float cg,float cb);
void multiplyimage(float cr,float cg,float cb);
void scaleimage(float s);
void vignette(float r0, float r1);
void downsample();
void upsample();
void flip(COLOR **im);
void zeroimage(COLOR **imzero);
void copyimage(COLOR **imto, COLOR **imfrom);
void accumulateimage(COLOR **imto, COLOR **imfrom, float scale);
void undistort(double cx, double cy, double k1, double k2, double k3);
void falloff(float r1, float r2, float g1, float g2, float b1, float b2);
void sumthresh(int x0, int y0, int x1, int y1, float thresh);
void invalidate(char *ppmfilename);
void fixbug();
void hardcode();
void help();

int  xsiz, ysiz;

#define MAX(x,y) (((x)>(y))?(x):(y))
#define MIN(x,y) (((x)<(y))?(x):(y))

char  *progname;

char  *fname;			/* the file names */
FILE  *fptr;			/* the file pointers */
COLOR  scale;			/* scaling factors */
int  nfile;			/* number of files */

int DBM;

COLOR **imorig, **imin, **imout, **imfinal;
COLOR black;

tabputs(s)			/* print line preceded by a tab */
     char  *s;
{
  putc('\t', stdout);
  fputs(s, stdout);
}

main(argc, argv)
     int  argc;
     char  *argv[];
{
  double  d;
  int  xres, yres;
  int  an;
  int ai0,ai1,ai2,ai3,ai4,ai5;
  float af0,af1,af2,af3,af4,af5;
  double ad0,ad1,ad2;
  char * as0,as1,as2;
  int accumulate = 0;

  progname = argv[0];

  DBM = 0;

  black[0] = 0;   black[1] = 0;   black[2] = 0; 
  
  if (argc < 2) {
    fprintf(stderr, "%s: Please specify a file or '-' for stdin, or '-h' for help\n", progname);
    quit(1);
  }

  an = 1; // Set argument counter to 1

  if (strcmp("-h",argv[an]) == 0) {
    fprintf(stderr, "Help on using pflare...\n");
    help();
  }

  if (argv[an][0] == '-') {   /* Read from stdin if first arg is '-' */
    fname = "<stdin>";
    fptr = stdin;
  } else { /* Read from file */
    fname = argv[1];
    if ((fptr = fopen(argv[1], "r")) == NULL) {
      fprintf(stderr, "%s: can't open file: %s\n",
	      progname, argv[1]);
      quit(1);
    }
    an++;
  }

  /* get header */
  fputs(fname, stdout);
  fputs(":\n", stdout);
  getheader(fptr, tabputs);
  /* get picture size */
  if (fgetresolu(&xres, &yres, fptr) != (YMAJOR|YDECR)) {
    if (DBM) fprintf(stderr, "%s: bad picture size\n", progname);
    quit(1);
  } else if (nfile == 0) {
    xsiz = xres;
    ysiz = yres;
  } else if (xres != xsiz || yres != ysiz) {
    if (DBM) fprintf(stderr, "%s: pictures are different sizes\n",
	    progname);
    quit(1);
  }

  if (DBM) fprintf(stderr,"Opened the file...\n");

  /* add new header info. */
  printargs(argc, argv, stdout);
  putchar('\n');
  fputresolu(YMAJOR|YDECR, xsiz, ysiz, stdout);

  if (DBM) fprintf(stderr,"Wrote args...\n");

  allocateandread();

  if (DBM) fprintf(stderr,"Allocated and Read...\n");

  /************************************************/
  /* Parse the rest of the command line arguments */
  /************************************************/

  for (;an < argc; an++)
    if (strcmp("-v",argv[an]) == 0) {
      DBM = 1;
      if (DBM) if (DBM) fprintf(stderr, "Verbose on...\n");
    }
    else if (strcmp("-gb",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian Blurring %d...\n", ai0);
      giblur0(ai0);
    }
    else if (strcmp("-xs",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "XShift %d...\n", ai0);
      xshift(ai0);
    }
    else if (strcmp("-ys",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "YShift %d...\n", ai0);
      yshift(ai0);
    }
    else if (strcmp("-gbx",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian Blurring X %d...\n", ai0);
      giblur0x(ai0);
    }
    else if (strcmp("-gby",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian Blurring Y %d...\n", ai0);
      giblur0y(ai0);
    }
    else if (strcmp("-gbxy",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian Blurring X Y %d...\n", ai0);
      giblur0xy(ai0);
    }
    else if (strcmp("-gbyx",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian Blurring Y X %d...\n", ai0);
      giblur0yx(ai0);
    }
    else if (strcmp("-gb1",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Gaussian+ Blurring %d...\n", ai0);
      giblur1(ai0);
    }
    else if (strcmp("-bb",argv[an]) == 0) {
      an++;
      ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Box Blurring %d...\n", ai0);
      boxblur0(ai0);
    }
    else if (strcmp("-se",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      an++; af1 = atof(argv[an]);
      an++; af2 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Sepia %6.4f %6.4f %6.4f...\n", af0,af1,af2);
      sepia(af0,af1,af2);
    }
    else if (strcmp("-mu",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      an++; af1 = atof(argv[an]);
      an++; af2 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Multiply %6.4f %6.4f %6.4f...\n", af0,af1,af2);
      multiplyimage(af0,af1,af2);
    }
    else if (strcmp("-sc",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Scale %6.4f...\n", af0);
      scaleimage(af0);
    }
    else if (strcmp("-vi",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      an++; af1 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Vignette %6.4f %6.4f...\n", af0,af1);
      vignette(af0,af1);
    }
    else if (strcmp("-fa",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      an++; af1 = atof(argv[an]);
      an++; af2 = atof(argv[an]);
      an++; af3 = atof(argv[an]);
      an++; af4 = atof(argv[an]);
      an++; af5 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Falloff %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f...\n", af0,af1,af2,af3,af4,af5);
      falloff(af0,af1,af2,af3,af4,af5);
    }
    else if (strcmp("-ud",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      an++; af1 = atof(argv[an]);
      an++; af2 = atof(argv[an]);
      an++; af3 = atof(argv[an]);
      an++; af4 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Undistort %6.4f %6.4f %6.4f %6.4f %6.4f...\n", af0,af1,af2,af3,af4);
      undistort(af0,af1,af2,af3,af4);
    }
    else if (strcmp("-ds",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "Downsample\n");
      downsample();
    }
    else if (strcmp("-fb",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "FixBug\n");
      fixbug();
    }
    else if (strcmp("-in",argv[an]) == 0) {
      an++; as0 = argv[an];
      if (DBM) fprintf(stderr, "Invalidate %s...\n", as0);
      invalidate(as0);
    }
    else if (strcmp("-us",argv[an]) == 0) {
      /* an++; af0 = atof(argv[an]); */
      if (DBM) fprintf(stderr, "Upsample\n");
      upsample();
    }
    else if (strcmp("-dz",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "Digital Zoom...\n");
      dzoom();
    }
    else if (strcmp("-ac",argv[an]) == 0) {
      an++; af0 = atof(argv[an]);
      if (DBM) fprintf(stderr, "Accumulating %6.4f...\n", af0);
      accumulate = 1;
      accumulateimage(imfinal,imin,af0);
      copyimage(imin,imorig);
    }
    else if (strcmp("-fl",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "Flipping...\n", ai0);
      flip(imin);
    }
    else if (strcmp("-ya",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "Yanking...\n", af0);
      copyimage(imin,imfinal);
    }
    else if (strcmp("-ze",argv[an]) == 0) {
      if (DBM) fprintf(stderr, "Zeroing...\n", af0);
      zeroimage(imfinal);
    }
    else if (strcmp("-ha",argv[an]) == 0) {
      an++; ai0 = atoi(argv[an]);
      if (DBM) fprintf(stderr, "Hard-Coded Operation...\n", af0);
      hardcode();
    }
    else if (strcmp("-st",argv[an]) == 0) {
      an++; ai0 = atoi(argv[an]);
      an++; ai1 = atoi(argv[an]);
      an++; ai2 = atoi(argv[an]);
      an++; ai3 = atoi(argv[an]);
      an++; af0 = atof(argv[an]);
      if (DBM) fprintf(stderr, "SumThreshhold %d %d %d %d %6.4f...\n", ai0,ai1,ai2,ai3,af0);
      sumthresh(ai0,ai1,ai2,ai3,af0);
    }
    else {
      fprintf(stderr, "%s: Unknown command: %s\n", progname, argv[an]);
      quit(1);
    }

  if (!accumulate) copyimage(imfinal,imin);

  if (DBM) fprintf(stderr,"About to write file...\n");

  writeandfree();

  if (DBM) fprintf(stderr,"Quitting...\n");

  quit(0);
}

swapimages() {
  COLOR ** temp;

  temp = imin;
  imin = imout;
  imout = temp;
}

allocateandread() {
  int  y;

  imin  = (COLOR **)malloc(ysiz*sizeof(COLOR*));
  imout = (COLOR **)malloc(ysiz*sizeof(COLOR*));
  imorig = (COLOR **)malloc(ysiz*sizeof(COLOR*));
  imfinal = (COLOR **)malloc(ysiz*sizeof(COLOR*));

  if (imin == NULL || imout == NULL) {
    if (DBM) fprintf(stderr, "%s: out of memory\n", progname);
    quit(1);
  }

  for (y=0; y<ysiz; y++) {
    imorig[y] = (COLOR *)malloc(xsiz*sizeof(COLOR));
    if (freadscan(imorig[y], xsiz, fptr) < 0) {
      if (DBM) fprintf(stderr, "%s: read error on file: %s\n",
	      progname, fname[0]);
      quit(1);
    }
    imin[y] = (COLOR *)malloc(xsiz*sizeof(COLOR));
    imout[y] = (COLOR *)malloc(xsiz*sizeof(COLOR));
    imfinal[y] = (COLOR *)malloc(xsiz*sizeof(COLOR));
  }

  copyimage(imin,imorig);
  zeroimage(imfinal);
}

writeandfree() {
  int  y;

  for (y=0; y<ysiz; y++) {
    if (fwritescan(imfinal[y], xsiz, stdout) < 0) {
      if (DBM) fprintf(stderr, "%s: write error\n", progname);
      quit(1);
    }
    free(imin[y]);
    free(imout[y]);
    free(imfinal[y]);
  }
  
  free(imin);
  free(imout);
  free(imfinal);
}

quit(code)		/* exit gracefully */
     int  code;
{
  exit(code);
}

/*
int mkgaussian(float *gauss,float s2) {
  int x;
  int fsize = (int)(s2 * 4);

  for (x=0
  gauss[x]=exp(-x*x/(2*S2))/sqrt(2*M_PI*S2);
}
*/

void giblur0(int twosigsquared) {

  giblur0x(twosigsquared);
  giblur0y(twosigsquared);

}

void giblur0x(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;

  scalefac=1;
  for (fx = 0; fx < twosigsquared; fx++)
    scalefac *= 0.25;

  for (y = 0; y < ysiz; y++)
    
    /* Blur in X */
    for (fx = 0; fx < twosigsquared; fx++) {
      
      copycolor(imout[y][0], imin[y][0]);
      addcolor(imout[y][0], imin[y][0]);
      addcolor(imout[y][0], imin[y][1]);

      copycolor(imout[y][xsiz-1], imin[y][xsiz-1]);
      addcolor(imout[y][xsiz-1], imin[y][xsiz-1]);
      addcolor(imout[y][xsiz-1], imin[y][xsiz-2]);

      for (x = 1; x < xsiz-1; x++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x-1]);
	addcolor(imout[y][x], imin[y][x+1]);
      }
      for (x = 0; x < xsiz; x++)
	copycolor(imin[y][x], imout[y][x]);
      
    }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
}

void giblur0y(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;

  scalefac=1;
  for (fx = 0; fx < twosigsquared; fx++)
    scalefac *= 0.25;

  /* Blur in Y */
  for (x = 0; x < xsiz; x++)
    
    for (fx = 0; fx < twosigsquared; fx++) {
      
      copycolor(imout[0][x], imin[0][x]);
      addcolor(imout[0][x], imin[0][x]);
      addcolor(imout[0][x], imin[1][x]);

      copycolor(imout[ysiz-1][x], imin[ysiz-1][x]);
      addcolor(imout[ysiz-1][x], imin[ysiz-1][x]);
      addcolor(imout[ysiz-1][x], imin[ysiz-2][x]);

      for (y = 1; y < ysiz-1; y++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y-1][x]);
	addcolor(imout[y][x], imin[y+1][x]);
      }
      for (y = 0; y < ysiz; y++)
	copycolor(imin[y][x], imout[y][x]);
      
    }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
}

void giblur0xy(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;

  scalefac=1;
  for (fx = 0; fx < twosigsquared; fx++)
    scalefac *= 0.25;

  for (fx = 0; fx < twosigsquared; fx++) {
    for (y = 1; y < ysiz-1; y++) {
      for (x = 1; x < xsiz-1; x++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y-1][x-1]);
	addcolor(imout[y][x], imin[y+1][x+1]);
      }
    }
    swapimages();
  }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
}

void giblur0yx(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;

  scalefac=1;
  for (fx = 0; fx < twosigsquared; fx++)
    scalefac *= 0.25;

  for (fx = 0; fx < twosigsquared; fx++) {
    for (y = 1; y < ysiz-1; y++) {
      for (x = 1; x < xsiz-1; x++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y+1][x-1]);
	addcolor(imout[y][x], imin[y-1][x+1]);
      }
    }
    swapimages();
  }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
}

giblur1(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;
  COLOR tempc;

  scalefac=1;
  for (fx = 0; fx < twosigsquared; fx++)
    scalefac *= 0.25;

  for (y = 0; y < ysiz; y++)
    
    /* Blur in X */
    for (fx = 0; fx < twosigsquared; fx++) {
      
      copycolor(imout[y][0], imin[y][0]);
      addcolor(imout[y][0], imin[y][0]);
      addcolor(imout[y][0], imin[y][1]);

      copycolor(imout[y][xsiz-1], imin[y][xsiz-1]);
      addcolor(imout[y][xsiz-1], imin[y][xsiz-1]);
      addcolor(imout[y][xsiz-1], imin[y][xsiz-2]);

      for (x = 2; x < xsiz-2; x++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x-1]);
	addcolor(imout[y][x], imin[y][x+1]);
	copycolor(tempc, imin[y][x-2]);
	scalecolor(tempc, 0.1);
	addcolor(imout[y][x], tempc);
	copycolor(tempc, imin[y][x+2]);
	scalecolor(tempc, 0.1);
	addcolor(imout[y][x], tempc);
      }
      for (x = 0; x < xsiz; x++)
	copycolor(imin[y][x], imout[y][x]);
      
    }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
  
  /* Blur in Y */
  for (x = 0; x < xsiz; x++)
    
    for (fx = 0; fx < twosigsquared; fx++) {
      
      copycolor(imout[0][x], imin[0][x]);
      addcolor(imout[0][x], imin[0][x]);
      addcolor(imout[0][x], imin[1][x]);

      copycolor(imout[ysiz-1][x], imin[ysiz-1][x]);
      addcolor(imout[ysiz-1][x], imin[ysiz-1][x]);
      addcolor(imout[ysiz-1][x], imin[ysiz-2][x]);

      for (y = 2; y < ysiz-2; y++) {
	copycolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y][x]);
	addcolor(imout[y][x], imin[y-1][x]);
	addcolor(imout[y][x], imin[y+1][x]);
	copycolor(tempc, imin[y-2][x]);
	scalecolor(tempc, 0.1);
	addcolor(imout[y][x], tempc);
	copycolor(tempc, imin[y+2][x]);
	scalecolor(tempc, 0.1);
	addcolor(imout[y][x], tempc);
      }
      for (y = 0; y < ysiz; y++)
	copycolor(imin[y][x], imout[y][x]);
      
    }
  
  /* Renormalize */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++)
      scalecolor(imin[y][x],scalefac);
}

#define  subcolor(c1,c2)	((c1)[0]-=(c2)[0],(c1)[1]-=(c2)[1],(c1)[2]-=(c2)[2])

#define  lumcolor(c1)	(0.3*(c1)[0]+0.6*(c1)[1]+0.1*(c1)[2])

boxblur0(int twosigsquared) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac;
  COLOR total;

  scalefac=1 / (float)(2*twosigsquared+1);

  for (y = 0; y < ysiz; y++) {
    
    /* Blur in X */
    setcolor(total,0,0,0);
    
    for (x = 0; x < twosigsquared; x++)
      addcolor(total, imin[y][x]);

    for (x = 0; x < twosigsquared; x++) {
      addcolor(total, imin[y][x+twosigsquared]);
      copycolor(imout[y][x], total);
    }

    for (x = twosigsquared; x < xsiz-twosigsquared; x++) {
      addcolor(total, imin[y][x+twosigsquared]);
      subcolor(total, imin[y][x-twosigsquared]);
      copycolor(imout[y][x], total);
    }

    for (x = xsiz-twosigsquared; x < xsiz; x++) {
      subcolor(total, imin[y][x-twosigsquared]);
      copycolor(imout[y][x], total);
    }

  }

  /* Renormalize and copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      scalecolor(imout[y][x],scalefac);
      copycolor(imin[y][x],imout[y][x]);
    }

  /* Blur in Y */
  for (x = 0; x < xsiz; x++) {
    
    setcolor(total,0,0,0);
    
    for (y = 0; y < twosigsquared; y++)
      addcolor(total, imin[y][x]);

    for (y = 0; y < twosigsquared; y++) {
      addcolor(total, imin[y+twosigsquared][x]);
      copycolor(imout[y][x], total);
    }

    for (y = twosigsquared; y < ysiz-twosigsquared; y++) {
      addcolor(total, imin[y+twosigsquared][x]);
      subcolor(total, imin[y-twosigsquared][x]);
      copycolor(imout[y][x], total);
    }

    for (y = ysiz-twosigsquared; y < ysiz; y++) {
      subcolor(total, imin[y-twosigsquared][x]);
      copycolor(imout[y][x], total);
    }

  }

  /* Renormalize and copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      scalecolor(imout[y][x],scalefac);
      copycolor(imin[y][x],imout[y][x]);
    }
}

/* Flip an image top-to-bottom and left-to-right */
void flip(COLOR **im) {
  register int x, y;

  COLOR c;

  for (y = 0; y < ysiz/2; y++)
    for (x = 0; x < xsiz; x++) {
      copycolor(c, im[ysiz-y-1][xsiz-x-1]);
      copycolor(im[ysiz-y-1][xsiz-x-1],im[y][x]);
      copycolor(im[y][x], c);
    }

  if (ysiz % 1 == 1) {  /* If there's an odd number of lines, do the middle one */
    y = ysiz/2;
    for (x = 0; x < xsiz/2; x++) {
      copycolor(c, im[ysiz-y-1][xsiz-x-1]);
      copycolor(im[ysiz-y-1][xsiz-x-1],im[y][x]);
      copycolor(im[y][x], c);
    }
  }

}

void copyimage(COLOR **imto, COLOR **imfrom) {
  register int  x;
  register int  y;

  /* Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      copycolor(imto[y][x],imfrom[y][x]);
    }
}

void zeroimage(COLOR **imzero) {
  register int  x;
  register int  y;
  COLOR zero;

  setcolor(zero,0.0,0.0,0.0);

  /* Zero */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      copycolor(imzero[y][x],zero);
    }
}

void accumulateimage(COLOR **imto, COLOR **imfrom, float scale) {
  register int  x;
  register int  y;
  COLOR temp;

  /* Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      copycolor(temp,imfrom[y][x]);
      scalecolor(temp,scale);
      addcolor(imto[y][x],temp);
    }
}

void sepia(float cr,float cg,float cb) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  COLOR sepiac;

  /* Renormalize and Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      setcolor(sepiac,cr,cg,cb);

      lum = lumcolor(imin[y][x]);
      scalecolor(sepiac,lum);
      
      copycolor(imin[y][x],sepiac);
    }
}

void multiplyimage(float cr,float cg,float cb) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  COLOR multc;

  setcolor(multc,cr,cg,cb);

  /* Renormalize and Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      multcolor(imin[y][x],multc);
    }
}

void scaleimage(float s) {
  int  fx, fy, y, i;
  register int  x;

  /* Renormalize and Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      scalecolor(imin[y][x],s);
    }
}

/* move things one byte over to fix temporary bug in DKODAK */
void fixbug() {
  int  fx, fy, y, i;
  register int  x;
  COLOR c0;
  COLOR c1;
  COLOR cf;

  for (y = 0; y < ysiz; y++) {
    for (x = 0; x < xsiz-1; x++) {
      copycolor(c0,imin[y][x]);
      copycolor(c1,imin[y][x+1]);

      cf[0] = c0[1];
      cf[1] = c0[2];
      cf[2] = c1[0];

      copycolor(imout[y][x],cf);
    }

    cf[0]=c1[1]; cf[1]=c1[2]; cf[2]=c1[0];
    copycolor(imout[y][x],cf);
  }
  copyimage(imin,imout);
}

void sumthresh(int x0, int y0, int x1, int y1, float thresh) {
  int x,y;
  COLOR sum,maxb;
  int area = 0;

  setcolor(sum,0,0,0);
  setcolor(maxb,0,0,0);

  /* Calculate brightest pixel */
  for (y = y0; y < y1; y++) {
    for (x = x0; x < x1; x++) {
      if (imin[ysiz-y-1][x][1] > maxb[1]) {
	copycolor(maxb,imin[ysiz-y-1][x]);
      }
    }
  }

  /* Calculate average of all pixels at least thresh * the max brightness */
  for (y = y0; y < y1; y++) {
    for (x = x0; x < x1; x++) {
      if (imin[ysiz-y-1][x][1] > thresh*maxb[1]) {
	addcolor(sum,imin[ysiz-y-1][x]);
	area++;
      }
    }
  }

  fprintf(stderr,"Sum %6.3f %6.3f %6.3f Area %d Average %6.3f %6.3f %6.3f\n",sum[0],sum[1],sum[2],area,sum[0]/(float)area,sum[1]/(float)area,sum[2]/(float)area);

}

void falloff(float r1, float r2, float g1, float g2, float b1, float b2) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  COLOR boost;
  float rr;
  float a,t;

  float xr = (float)xsiz /2;
  float yr = (float)ysiz /2;

  for (y = 0; y < ysiz; y++) {
    for (x = 0; x < xsiz; x++) {

      /* Compute normalized squared radius */
      rr = ((x-xr)*(x-xr) + (y-ysiz/2)*(y-ysiz/2)) / (xr*xr);

      boost[0] = ((rr*r2 + r1) * rr + 1.0);
      boost[1] = ((rr*g2 + g1) * rr + 1.0);
      boost[2] = ((rr*b2 + b1) * rr + 1.0);

      multcolor(imin[y][x],boost);
    }
  }
}

void vignette(float r0, float r1) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  double radius;
  float reduce,rscale;
  float a,t;
  float bsiz;

  /* Renormalize and Copy */
  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {

      radius = sqrt((x-xsiz/2)*(x-xsiz/2) + (y-ysiz/2)*(y-ysiz/2));

      if (radius > r0) {
	if (radius < r1) {
	  t = 1.0 - (radius-r0)/(r1-r0);
	  a = t * 2 - 1;
	  reduce = (0.25 * M_PI + 0.5 * asin(a) + 0.5 * a * sqrt( 1 - a*a ))/(0.5 * M_PI);
	  scalecolor(imin[y][x],reduce);
	} else {
	  setcolor(imin[y][x],0.0,0.0,0.0);
	}
      }
    }
}

void undistort(double cx, double cy, double k1, double k2, double k3) {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  COLOR c0,c1,c2,c3,black;
  double radius;
  float reduce,rscale;
  float a,t;
  float bsiz;

  int funsize;
  double* distfun;
  double* idistfun;
  double r,rf,qxf,qyf;
  int ri,qxi,qyi;

  int width = xsiz;
  int height = ysiz;

  double px, py;
  double qx, qy;

  funsize = 2 * MAX(height,width);
  distfun = (double*)calloc(funsize,sizeof(double));
  idistfun = (double*)calloc(funsize,sizeof(double));

  // Sample forward function
  for (x=0;x<funsize;x++) {
    r = 0.001 * (double)x;
    distfun[x] = (double)x * (1 + r*r * (k1 + r*r*(k2 + r*r*k3)));
  }

  // Resample for backward function
  x = 0;
  y = 0;
  while (x<(funsize-1) && y<funsize) {
    while (y<funsize && distfun[x+1] >= (double)y) {
      idistfun[y] = (double)x + ((double)y - distfun[x]) / (distfun[x+1]-distfun[x]);
      if (y) idistfun[y] /= (double)y;
      else idistfun[y] = 1;
      y++;
    }
    x++;
  }
  free(distfun);

  /* Perform Undistortion */
  for (y = 0; y < ysiz; y++) {
    for (x = 0; x < xsiz; x++) {

      px = (double)x - cx;
      py = (double)y - cy;

      r = sqrt(px*px + py*py);
      ri = (int)r;     // Linearly Interpolate
      rf = r - (double)ri;

      if ((ri+1 < funsize)) {
	qx = cx + px * ((1-rf)*idistfun[ri] + rf*idistfun[ri+1]);
	qy = cy + py * ((1-rf)*idistfun[ri] + rf*idistfun[ri+1]);

/*        copycolor(imout[y][x], imin[(int)qy][(int)qx]);   - no bilerp*/

        /* Perform bilinear interpolation */
	qxi = (int) qx;
        qxf = qx - qxi;
	qyi = (int) qy;
        qyf = qy - qyi;

	/* Make sure point is not off other image boundary */
	if ( (qxi >= 0) && (qxi < xsiz-1) && (qyi >= 0) && (qyi < ysiz-1) ) {
	  copycolor(c0,imin[qyi][qxi]);
	  scalecolor(c0,1-qxf);
	  scalecolor(c0,1-qyf);

	  copycolor(c1,imin[qyi][qxi+1]);
	  scalecolor(c1,qxf);
	  scalecolor(c1,1-qyf);

	  copycolor(c2,imin[qyi+1][qxi]);
	  scalecolor(c2,1-qxf);
	  scalecolor(c2,qyf);

	  copycolor(c3,imin[qyi+1][qxi+1]);
          scalecolor(c3,qxf);
	  scalecolor(c3,qyf);

          copycolor(imout[y][x], c0);
          addcolor(imout[y][x], c1);
          addcolor(imout[y][x], c2);
          addcolor(imout[y][x], c3);
        }
        else copycolor(imout[y][x],black);
      }
    }
  }

  free(idistfun);
  copyimage(imin,imout);
}

void downsample() {
  int  fx, fy, y, i;
  register int  x;
  float scalefac,lum;
  double radius;
  float reduce,rscale;
  float a,t;
  float bsiz;

  /* Shrink */
  for (y = 0; y < ysiz/2; y++)
    for (x = 0; x < xsiz/2; x++) {
      copycolor(imout[y][x],imin[2*y][2*x]);
      addcolor(imout[y][x],imin[2*y][2*x+1]);
      addcolor(imout[y][x],imin[2*y+1][2*x+1]);
      addcolor(imout[y][x],imin[2*y+1][2*x]);
      scalecolor(imout[y][x],0.25);
    }

  swapimages();
  xsiz /= 2;
  ysiz /= 2;
}

void upsample() {
  register int  x, y;

  /* Expand */
  xsiz *= 2;
  ysiz *= 2;
  for (y = 0; y < ysiz/2; y++)
    for (x = 0; x < xsiz/2; x++) {
      copycolor(imout[2*y][2*x],imin[y][x]);
      copycolor(imout[2*y][2*x+1],imin[y][x]);
      copycolor(imout[2*y+1][2*x],imin[y][x]);
      copycolor(imout[2*y+1][2*x+1],imin[y][x]);
    }
  swapimages();

  giblur0(1);
}

void dzoom() {
  register int  x, y;

  for (y = 0; y < ysiz; y++)
    for (x = 0; x < xsiz; x++) {
      copycolor(imout[y][x],imin[(y+ysiz/2)/2][(x+xsiz/2)/2]);
    }
  swapimages();

  giblur0(1);
}

void xshift(int dx) {
  register int  x, y;

  for (y = 0; y < ysiz; y++) {
    for (x = 0; x < dx; x++) {
      copycolor(imout[y][x],black);
    }
    for (x = MAX(dx,0); x < MIN(xsiz+dx,xsiz); x++) {
      copycolor(imout[y][x],imin[y][x-dx]);
    }
    for (x = xsiz+dx; x < xsiz; x++) {
      copycolor(imout[y][x],black);
    }
  }
  swapimages();
}

void yshift(int dy) {
  register int  x, y;

  for (x = 0; x < xsiz; x++) {
    for (y = 0; y < dy; y++) {
      copycolor(imout[y][x],black);
    }
    for (y = MAX(dy,0); y < MIN(ysiz+dy,ysiz); y++) {
      copycolor(imout[y][x],imin[y-dy][x]);
    }
    for (y = ysiz+dy; y < ysiz; y++) {
      copycolor(imout[y][x],black);
    }
  }
  swapimages();
}

void hardcode() {
  downsample();
  downsample();
  downsample();
  giblur0(12);
  accumulateimage(imfinal,imin,0.1);
  giblur0(28);
  accumulateimage(imfinal,imin,0.05);
  copyimage(imin,imfinal);
  copyimage(imfinal,imin);
  
  xsiz *= 8;  ysiz *= 8;
  copyimage(imin,imorig);
  downsample();
  downsample();
  downsample();
  giblur0(3);
  accumulateimage(imfinal,imin,0.1);
  copyimage(imin,imfinal);
  upsample();
  upsample();
  copyimage(imfinal,imin);

  xsiz *= 2;  ysiz *= 2;
  copyimage(imin,imorig);
  downsample();
  giblur0(1);
  accumulateimage(imfinal,imin,0.75);
  copyimage(imin,imfinal);
  upsample();
}

void invalidate(char *ppmfilename) {
  
  int bpp;
  int i, x, y;
  char strbuf[100];
  unsigned char row[40960];
  int depth,width,height;
  FILE *image;
  int r,g,b;
  
  if( (image=fopen(ppmfilename,"r")) == NULL ) {
    if (DBM) fprintf(stderr,"Load_PPM: can't open input file %s\n",fname);
    return;
  }
  strcpy(strbuf,"#");while (strbuf[0] == '#') fgets(strbuf,100,image);
  if (strncmp(strbuf,"P6",2)==0)   depth = 3;
  else if (strncmp(strbuf,"Pf",2)==0)   depth = 4; // Floating point
  else depth = 1;
  strcpy(strbuf,"#");while (strbuf[0] == '#') fgets(strbuf,100,image);
  sscanf(strbuf,"%d %d",&width,&height);
  strcpy(strbuf,"#");while (strbuf[0] == '#') fgets(strbuf,100,image);

  // The last line contains the bits per channel info; we'll require 8
  sscanf(strbuf,"%d",&bpp);
  if (bpp != 255) {
    if (DBM) fprintf(stderr,"Read_PPM: Must use 8 bits per channel\n");
    return;
  }

  if (DBM) fprintf(stderr,"Size is %dx%d, %d channels.\n",width,height,depth);

  if (depth == 3) {
    for (y=0;y<height;y++) {
      fread(row,1,width * depth,image);
      for (x=0;x<width;x++) {
        unsigned char *pp = row + x * depth;
        /* SetRGBint(x,height - 1 - y,*pp,*(pp+1),*(pp+2)); */
	r = *pp; g = *(pp+1); b = *(pp+2);
	if (b == 255 && r == 0 && g == 0) {
	  setcolor(imin[y][x],0,0,0);
	}
      }
    }
  }
  fclose(image);
}

void help() {
  fprintf(stderr,"\npflare - blur, vignette, and tint RADIANCE pictures\n\n");
  fprintf(stderr,"0\n");
  fprintf(stderr,"-gb n                Gaussian Blur, n=2*sigma^2\n");
  fprintf(stderr,"-gbx n               Gaussian Blur (X only)\n");
  fprintf(stderr,"-gby n               Gaussian Blur (Y only)\n");
  fprintf(stderr,"-bb n                Box Blur\n");
  fprintf(stderr,"-se r g b            Tint to color (r,g,b)\n");
  fprintf(stderr,"-mu r g b            Tint to color (r,g,b)\n");
  fprintf(stderr,"-sc s                Scale image by s)\n");
  fprintf(stderr,"-vi r0 r1            Vignette with falloff radii r0 and r1\n");
  fprintf(stderr,"-ud cx cy k0 k1 k2   Undistort\n");
  fprintf(stderr,"-ds                  Downsample\n");
  fprintf(stderr,"-in fn               Invalidate from ppm\n");
  fprintf(stderr,"-ac x                Accumulate\n");
  fprintf(stderr,"-ya                  Yank (copy ifinal in imin)\n");
  fprintf(stderr,"-ze                  Zero final image\n");
  fprintf(stderr,"-ha n                Hard-Coded Operation #n\n");
  fprintf(stderr,"-ac x                Sum with Threshhold\n");
  
  fprintf(stderr,"\n");
  exit(1);
}