/*
 * kinostill.c
 *
 * version 0.1  2007-01-28
 *
 * copyright 2007  Mark J. Stock  mstock@umich.edu
 *
 *
 * Kinostill is a program for making animations from still images. It is
 * kind like stills2dv and slcreator, but it has no gui and does not do audio.
 *   http://www.deniscarl.com/stills2dv/
 *   http://slcreator.sourceforge.net/
 *
 * It does not assemble multiple videos together, try kino.
 *   http://www.kinodv.org/
 *
 * It does not do transitions, try any of the above.
 *
 * It does, however, have better interpolation algorithms so your
 * resulting video is not jumpy or blocky. So it does one thing, and I'd
 * like to think that it does it very well.
 *
 * It can support making images sequences with up to 100000 frames and
 * up to 32000 pixels, although these are only arbitrary limitations
 * designed to frustrate users from the year 2020.
 *
 * It can support smooth pans and zooms (Ken Burns-style), but only linear
 * ones for now. I am sure splines and acceleration will appear in a future
 * release if one person asks for them. I'm a pretty nice guy. Motion blur
 * is a possibility, too.
 *
 * It uses bilinear interpolation for both zooming in and out. If you 
 * have a better algorithm (and you probably do if you're reading this)
 * please implement it. And if you're feeling generous, please send it
 * back to me so that I can use it too.
 *
 * You can build the binary on Linux with:
 *
 * gcc -O2 -ffast-math -fomit-frame-pointer -funroll-loops -o kinostill kinostill.c -lpnm
 *
 * Or, if you received a Makefile with the code, type "make"
 *
 * Despite all of the optimizations, the program is still kind of a dog.
 *
 * To view the help file, run "kinostill -h", I'll repeat it here:
 *
 * usage:
 *   kinostill [options]
 *
 * where [-options] are zero or more of the following:
 *
 *    [-n num]    number of frames to output; default=100
 *
 *    [-x xres]   resolution in x-direction; default=720
 *
 *    [-y yres]   resolution in y-direction; default=480
 *
 *    [-xy xres [yres]]  resolution in x- and y-directions, or both;
 *                default=720 480
 *
 *    [-zoom start [end]]  magnification at start and end of movement;
 *                1.0 means that the input image is as large as is possible
 *                to completely maintain inside of the output image when
 *                the center is set to 0.5 0.5; if only one argument, that
 *                will be the zoom setting for the whole movement; default=1
 *
 *    [-sc x [y]]  set the center of the output images for the start of the
 *                movement to the given coordinates in the input image;
 *                coordinates begin at the bottom left and are normalized by
 *                the image size (0..1 in both directions); default is
 *                centered=0.5 0.5
 *
 *    [-ec x [y]]  set the center of the output images for the end of the
 *                movement to the given coordinates in the input image;
 *                coordinates are the same as above; default is start center
 *
 *                Example of movement definition: to make a series of 720x480
 *                images starting with the exact upper left quadrant of a
 *                1440x960 image and ending on the lower right quadrant, run
 *
 *                kinostill -zoom 0.5 -sc 0.25 0.75 -ec 0.75 0.25 < in.ppm
 *
 *    [-prefix name]  prefix for all output files; default=out_
 *
 *    [-help]     returns this help information
 *
 * Options may be abbreviated to an unambiguous length.
 * Output is to a series of PNM files.
 *
 *
 * copyright 2007  Mark J. Stock  mstock@umich.edu
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

#define TRUE 1
#define FALSE 0
#define MAXFRAMES 100000
#include <stdio.h>
#include <math.h>
#include <pam.h>

void makeImageSeries (int, float*, float, float*, float, struct pam*, struct pam*, tuple **indata, char*);

int main (int argc, char **argv) {

  // this struct defined in pam.h
  //   struct pam {
  //     int size
  //     int len
  //     FILE * file
  //     int format
  //     int plainformat
  //     int height
  //     int width
  //     int depth
  //     sample maxval
  //     int bytes_per_sample
  //     char tuple_type[256]; }
  struct pam inpam, outpam;
  // this holds all of the image data
  tuple **indata;
  // and this is the stuff that I use
  int numFrames;
  int xres, yres;
  int setStart,setEnd;
  float startCenter[2],endCenter[2];
  float startScale,endScale;
  char prefix[80];
  char progname[100];
  // if we didn't have to declare i, we might as well use F77
  int i;

  // define default parameters
  numFrames = 100;
  xres = 720;
  yres = 480;
  // define default start and end positions and scales
  setStart = FALSE;
  startCenter[0] = 0.5;
  startCenter[1] = 0.5;
  startScale = 1.0;
  setEnd = FALSE;
  endCenter[0] = 0.5;
  endCenter[1] = 0.5;
  endScale = 1.0;
  sprintf(prefix,"out_");

  // process command-line parameters

  // first, this is suggested by "man 3 libpnm"
  pnm_init(&argc, argv);

  // now we check them
  (void) strcpy(progname,argv[0]);
  if (argc < 1) (void) Usage(progname,0);
  for (i=1; i<argc; i++) {
    if (strncmp(argv[i], "-xy", 3) == 0) {
      xres = atoi(argv[++i]);
      if (argv[i+1][0] != '-') {
        yres = atoi(argv[++i]);
      } else {
        yres = xres;
      }
    } else if (strncmp(argv[i], "-x", 2) == 0) {
      xres = atoi(argv[++i]);
    } else if (strncmp(argv[i], "-y", 2) == 0) {
      yres = atoi(argv[++i]);
    } else if (strncmp(argv[i], "-zoom", 2) == 0) {
      startScale = atof(argv[++i]);
      if (argv[i+1][0] != '-') {
        endScale = atof(argv[++i]);
      } else {
        endScale = startScale;
      }
    } else if (strncmp(argv[i], "-n", 2) == 0) {
      numFrames = atoi(argv[++i]);
    } else if (strncmp(argv[i], "-sc", 2) == 0) {
      setStart = TRUE;
      startCenter[0] = atof(argv[++i]);
      if (argv[i+1][0] != '-') {
        startCenter[1] = atof(argv[++i]);
      } else {
        startCenter[1] = startCenter[0];
      }
    } else if (strncmp(argv[i], "-ec", 2) == 0) {
      setEnd = TRUE;
      endCenter[0] = atof(argv[++i]);
      if (argv[i+1][0] != '-') {
        endCenter[1] = atof(argv[++i]);
      } else {
        endCenter[1] = endCenter[0];
      }
    } else if (strncmp(argv[i], "-prefix", 2) == 0) {
      strcpy(prefix,argv[++i]);
    } else
      (void) Usage(progname,0);
  }

  // check the command-line parameters (after pnm_init is done with them)
  if (numFrames < 1) {
    fprintf(stderr,"Too few frames (%d); resetting to 100\n",numFrames);
    numFrames = 100;
  }
  if (numFrames > MAXFRAMES) {
    fprintf(stderr,"Too many frames (%d); are you sure?\n",numFrames);
    exit(1);
  }
  if (xres < 1 || xres > 32000) {
    fprintf(stderr,"X-resolution (%d) too weird.\n",xres);
    exit(1);
  }
  if (yres < 1 || yres > 32000) {
    fprintf(stderr,"Y-resolution (%d) too weird.\n",yres);
    exit(1);
  }
  if (startScale < 0.) {
    fprintf(stderr,"Zoom setting (%d) negative. I can't do that, Dave.\n",startScale);
    exit(1);
  }
  if (endScale < 0.) {
    fprintf(stderr,"Zoom setting (%d) negative. What are you thinking?\n",endScale);
    exit(1);
  }
  if (setStart && !setEnd) {
    endCenter[0] = startCenter[0];
    endCenter[1] = startCenter[1];
  }
  if (!setStart && setEnd) {
    startCenter[0] = endCenter[0];
    startCenter[1] = endCenter[1];
  }

  // read in the still image, return a pam struct and fill the data
  indata = pnm_readpam(stdin, &inpam, sizeof(inpam));

  // set the output image characteristics
  outpam = inpam;
  outpam.width = xres;
  outpam.height = yres;

  // process the data and write the images
  makeImageSeries (numFrames, startCenter, startScale,
    endCenter, endScale, &inpam, &outpam, indata, prefix);

  // help the user with some suggestions
  //if (outFile != stdout) {
    fprintf(stderr,"\nIf you'd like to encode these into a dv stream, try\n");
    fprintf(stderr,"   encodedv -e %d %s%%05d.ppm > output.dv\n",numFrames,prefix);
    fprintf(stderr,"\nIf you'd like to encode these into an mpeg2 file, try\n");
    fprintf(stderr,"   cat %s*.ppm | ppmtoy4m -S 420mpeg2 | mpeg2enc -f 8 -o output.mpg\n\n",prefix);
  //}

  // let's assume everything went well....?
  exit;
}


/*
 * Make the series
 */
void makeImageSeries ( int numFrames,
    float *startCenter, float startScale, float *endCenter, float endScale,
    struct pam *inpam, struct pam *outpam, tuple **indata,
    char *prefix) {

  int iFrame;
  unsigned int row, col, inrow, incol, plane;
  int srow,erow,scol,ecol;
  float thisCenter[2],thisScale,pixelScale;
  float tempdata[outpam->height][outpam->width][outpam->depth];
  float tempwt[outpam->height][outpam->width];
  float inY,inX,outX,outY;
  float rowConst1,rowConst2,colConst1,colConst2;
  float weight0,weight1,weight00,weight01,weight10,weight11;
  int inrow0,incol0,inrow1,incol1;
  int outrow0,outcol0,outrow1,outcol1;
  tuple **outdata;
  char outFileName[100];
  FILE *outFile;

  // allocate the data for the output
  outdata = pnm_allocpamarray (outpam);

  for (iFrame = 0; iFrame < numFrames; iFrame++) {

    // set center and scale for this image
    thisCenter[0] = startCenter[0] +
                    iFrame * (endCenter[0]-startCenter[0]) / (numFrames-1);
    thisCenter[1] = startCenter[1] +
                    iFrame * (endCenter[1]-startCenter[1]) / (numFrames-1);
    thisScale = startScale + iFrame * (endScale-startScale) / (numFrames-1);
    //fprintf(stderr,"frame %d, center %g %g, scale %g\n",iFrame,
    //  thisCenter[0],thisCenter[1],thisScale);

    // are the output pixels larger than the input pixels?
    if (outpam->width/(float)outpam->height >
        inpam->width/(float)inpam->height) {
      // output image is "wider" than input image
      // use width to denote pixel scale
      pixelScale = outpam->width/(inpam->width*thisScale);
      fprintf(stderr,"  frame %d, pixel scale is %g, using widths\n",
        iFrame,pixelScale);
    } else {
      // output image is "narrower" than input image
      // use height to denote pixel scale
      pixelScale = outpam->height/(inpam->height*thisScale);
      fprintf(stderr,"  frame %d, pixel scale is %g, using heights\n",
        iFrame,pixelScale);
    }

    // if output pixels are smaller than input pixels, loop over
    //   output pixels and interpolate from input
    if (pixelScale > 1.) {

      // set some useful parameters
      rowConst1 = inpam->height*(1.-thisCenter[1]);
      rowConst2 = 0.5*outpam->height;
      colConst1 = inpam->width*thisCenter[0];
      colConst2 = 0.5*outpam->width;
      //fprintf(stderr,"  const %g %g %g %g\n",rowConst1,rowConst2,colConst1,colConst2);

      // for each pixel in the output image...
      for (row = 0; row < outpam->height; row++) {
        for (col = 0; col < outpam->width; col++) {

          // find the coordinates of the center of the output pixel
          //   in the input image frame
          inY = rowConst1+(row+0.5-rowConst2)/pixelScale;
          inX = colConst1+(col+0.5-colConst2)/pixelScale;
          //fprintf(stderr,"  inYX %g %g\n",inY,inX);

          // find out which pixel from the input image goes here
          //inrow = (int)(row/pixelScale)%inpam->height;
          //incol = (int)(col/pixelScale)%inpam->width;
          inrow0 = (int)(inY)%inpam->height;
          incol0 = (int)(inX)%inpam->width;
          inrow1 = (int)(inY+1)%inpam->height;
          incol1 = (int)(inX+1)%inpam->width;
          //fprintf(stderr,"  index %d %d %d %d\n",inrow0,incol0,inrow1,incol1);

          weight0 = inY-(int)inY;
          weight1 = inX-(int)inX;
          weight00 = (1.-weight0)*(1.-weight1);
          weight01 = (1.-weight0)*weight1;
          weight10 = weight0*(1.-weight1);
          weight11 = weight0*weight1;
          //fprintf(stderr,"  weights %g %g %g %g\n",weight00,weight01,weight10,weight11);
          //fflush(stderr);

          // copy values for each channel from that pixel
          //for (plane = 0; plane < outpam->depth; plane++) {
          //  outdata[row][col][plane] = indata[inrow][incol][plane];
          //}

          // bilinearly interpolate from the original image
          for (plane = 0; plane < outpam->depth; plane++) {
            outdata[row][col][plane] =
              weight00 * indata[inrow0][incol0][plane] +
              weight01 * indata[inrow0][incol1][plane] +
              weight10 * indata[inrow1][incol0][plane] +
              weight11 * indata[inrow1][incol1][plane];
          }
          
        }
      }

    // if output pixels are larger than input pixels, loop over 
    //   input pixels and sum onto output
    } else {  // if (pixelScale < 1.)

      // zero each pixel in the output image
      for (row = 0; row < outpam->height; row++) {
        for (col = 0; col < outpam->width; col++) {
          for (plane = 0; plane < outpam->depth; plane++) {
            tempdata[row][col][plane] = 0.;
          }
          tempwt[row][col] = 0.;
        }
      }

      // set some useful parameters
      rowConst2 = inpam->height*(1.-thisCenter[1]);
      rowConst1 = 0.5*outpam->height;
      colConst2 = inpam->width*thisCenter[0];
      colConst1 = 0.5*outpam->width;
      //fprintf(stderr,"  const %g %g %g %g\n",rowConst1,rowConst2,colConst1,colConst2);

      // compute the row and column range from the input image
      srow = rowConst2 - 0.5 - (2+rowConst1)/pixelScale;
      erow = rowConst2 - 0.5 - (rowConst1-2-outpam->height)/pixelScale;
      scol = colConst2 - 0.5 - (2+colConst1)/pixelScale;
      ecol = colConst2 - 0.5 - (colConst1-2-outpam->width)/pixelScale;
      //fprintf(stderr,"rows %d %d  cols %d %d\n",srow,erow,scol,ecol);
      if (srow < 0) srow = 0;
      if (erow > inpam->height) erow = inpam->height;
      if (scol < 0) scol = 0;
      if (ecol > inpam->width) ecol = inpam->width;

      // for each pixel in the input image...
      //for (inrow = 0; inrow < inpam->height; inrow++) {
      //  for (incol = 0; incol < inpam->width; incol++) {
      // for each effective pixel
      for (inrow = srow; inrow < erow; inrow++) {
        for (incol = scol; incol < ecol; incol++) {

          // find the coordinates of the center of the input pixel
          //   in the output image frame (note * instead of /)
          outY = rowConst1+(inrow+0.5-rowConst2)*pixelScale;
          outX = colConst1+(incol+0.5-colConst2)*pixelScale;

          //if (outY > -2 && outY < outpam->height + 2) {
          //if (outX > -2 && outX < outpam->width + 2) {

          //fprintf(stderr,"outYX %g %g\n",outY,outX);

          // copy values for each channel from that pixel
          // find out which pixel from the input image goes here
          //row = (int)(inrow*pixelScale)%outpam->height;
          //col = (int)(incol*pixelScale)%outpam->width;
          outrow0 = (int)(outY+5) - 5;
          outcol0 = (int)(outX+5) - 5;
          //outrow1 = (int)(outY+1)%outpam->height;
          //outcol1 = (int)(outX+1)%outpam->width;
          outrow1 = outrow0 + 1;
          outcol1 = outcol0 + 1;
          //fprintf(stderr,"  index %d %d %d %d\n",outrow0,outcol0,outrow1,outcol1);

          //weight0 = outY-(int)outY;
          //weight1 = outX-(int)outX;
          weight0 = outY-outrow0;
          weight1 = outX-outcol0;
          weight00 = (1.-weight0)*(1.-weight1);
          weight01 = (1.-weight0)*weight1;
          weight10 = weight0*(1.-weight1);
          weight11 = weight0*weight1;

          // bilinearly interpolate from the original image
          if (outpam->depth == 1) {
            if (outrow0 > -1 && outrow0 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                tempdata[outrow0][outcol0][0] +=
                  weight00 * (float)indata[inrow][incol][0];
                tempwt[outrow0][outcol0] += weight00;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                tempdata[outrow0][outcol1][0] +=
                  weight01 * (float)indata[inrow][incol][0];
                tempwt[outrow0][outcol1] += weight01;
              }
            }
            if (outrow1 > -1 && outrow1 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                tempdata[outrow1][outcol0][0] +=
                  weight10 * (float)indata[inrow][incol][0];
                tempwt[outrow1][outcol0] += weight10;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                tempdata[outrow1][outcol1][0] +=
                  weight11 * (float)indata[inrow][incol][0];
                tempwt[outrow1][outcol1] += weight11;
              }
            }

          } else if (outpam->depth == 3) {
            if (outrow0 > -1 && outrow0 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                for (plane = 0; plane < 3; plane++) {
                  tempdata[outrow0][outcol0][plane] +=
                    weight00 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow0][outcol0] += weight00;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                for (plane = 0; plane < 3; plane++) {
                  tempdata[outrow0][outcol1][plane] +=
                    weight01 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow0][outcol1] += weight01;
              }
            }
            if (outrow1 > -1 && outrow1 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                for (plane = 0; plane < 3; plane++) {
                  tempdata[outrow1][outcol0][plane] +=
                    weight10 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow1][outcol0] += weight10;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                for (plane = 0; plane < 3; plane++) {
                  tempdata[outrow1][outcol1][plane] +=
                    weight11 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow1][outcol1] += weight11;
              }
            }

          } else {
            if (outrow0 > -1 && outrow0 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                for (plane = 0; plane < outpam->depth; plane++) {
                  tempdata[outrow0][outcol0][plane] +=
                    weight00 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow0][outcol0] += weight00;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                for (plane = 0; plane < outpam->depth; plane++) {
                  tempdata[outrow0][outcol1][plane] +=
                    weight01 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow0][outcol1] += weight01;
              }
            }
            if (outrow1 > -1 && outrow1 < outpam->height) {
              if (outcol0 > -1 && outcol0 < outpam->width) {
                for (plane = 0; plane < outpam->depth; plane++) {
                  tempdata[outrow1][outcol0][plane] +=
                    weight10 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow1][outcol0] += weight10;
              }
              if (outcol1 > -1 && outcol1 < outpam->width) {
                for (plane = 0; plane < outpam->depth; plane++) {
                  tempdata[outrow1][outcol1][plane] +=
                    weight11 * (float)indata[inrow][incol][plane];
                }
                tempwt[outrow1][outcol1] += weight11;
              }
            }
          }


          //}
          //}
        }
      }

      // then, scale the tempdata to the output pixels
      //scale = 1.*(pixelScale*pixelScale);
      if (outpam->bytes_per_sample == 1) {
        for (row = 0; row < outpam->height; row++) {
        for (col = 0; col < outpam->width; col++) {
        for (plane = 0; plane < outpam->depth; plane++) {
          //outdata[row][col][plane] = (unsigned char)(scale * tempdata[row][col][plane]);
          outdata[row][col][plane] = (unsigned char)(tempdata[row][col][plane]/tempwt[row][col]);
        }
        }
        }
      } else if (outpam->bytes_per_sample == 2) {
        for (row = 0; row < outpam->height; row++) {
        for (col = 0; col < outpam->width; col++) {
        for (plane = 0; plane < outpam->depth; plane++) {
          //outdata[row][col][plane] = (unsigned char)(scale * tempdata[row][col][plane]);
          outdata[row][col][plane] = (unsigned short int)(tempdata[row][col][plane]/tempwt[row][col]);
        }
        }
        }
      } else {
        fprintf(stderr,"ERROR: input image has more than 2 bytes per sample\n");
        fprintf(stderr,"Quitting.\n");
        exit(0);
      }
    }

    // write the image
    sprintf(outFileName,"%s%05d.ppm",prefix,iFrame);
    outFile = fopen(outFileName,"wb");
    outpam->file = outFile;
    pnm_writepam (outpam, outdata);
    fclose(outFile);
  }

  return;
}


/*
 * This function writes basic usage information to stderr,
 * and then quits. Too bad.
 */
int Usage(char progname[80],int status) {

   static char **cpp, *help_message[] = {
   "where [-options] are zero or more of the following:                        ",
   "                                                                           ",
   "   [-n num]    number of frames to output; default=100                     ",
   "                                                                           ",
   "   [-x xres]   resolution in x-direction; default=720                      ",
   "                                                                           ",
   "   [-y yres]   resolution in y-direction; default=480                      ",
   "                                                                           ",
   "   [-xy xres [yres]]  resolution in x- and y-directions, or both;          ",
   "               default=720 480                                             ",
   "                                                                           ",
   "   [-zoom start [end]]  magnification at start and end of movement;        ",
   "               1.0 means that the input image is as large as is possible   ",
   "               to completely maintain inside of the output image when      ",
   "               the center is set to 0.5 0.5; if only one argument, that    ",
   "               will be the zoom setting for the whole movement; default=1  ",
   "                                                                           ",
   "   [-sc x [y]]  set the center of the output images for the start of the   ",
   "               movement to the given coordinates in the input image;       ",
   "               coordinates begin at the bottom left and are normalized by  ",
   "               the image size (0..1 in both directions); default is        ",
   "               centered=0.5 0.5                                            ",
   "                                                                           ",
   "   [-ec x [y]]  set the center of the output images for the end of the     ",
   "               movement to the given coordinates in the input image;       ",
   "               coordinates are the same as above; default is start center  ",
   "                                                                           ",
   "               Example of movement definition: to make a series of 720x480 ",
   "               images starting with the exact upper left quadrant of a     ",
   "               1440x960 image and ending on the lower right quadrant, run  ",
   "                                                                           ",
   "               kinostill -zoom 0.5 -sc 0.25 0.75 -ec 0.75 0.25 < in.ppm   ",
   "                                                                           ",
   "   [-prefix name]  prefix for all output files; default=out_               ",
   "                                                                           ",
   "   [-help]     returns this help information                               ",
   " ",
   "Options may be abbreviated to an unambiguous length.",
   "Output is to a series of PNM files.",
   NULL
   };

   fprintf(stderr, "usage:\n  %s [options]\n\n", progname);
   for (cpp = help_message; *cpp; cpp++) fprintf(stderr, "%s\n", *cpp);
   fflush(stderr);
   exit(status);
   return(0);
}