/*********************************************************/
/*   program 22: first-order polynomial warp             */
/*********************************************************/

#define  XSIZE  640
#define  YSIZE  480

main()
  {
  unsigned char *ir,*ig,*ib;
  unsigned char *or,*og,*ob;
  long  int  ix,iy;
  long  int  ox,oy;
  long  int  x1,x2,x3,x4;
  long  int  y1,y2,y3,y4;
  double     xa0,xa1,xa2,xa3;
  double     ya0,ya1,ya2,ya3;
  double     IX,IY;
  double     OX,OY;
  double     DX,DY;

  /* set corner warp control points */
  x1 = 140; y1 = 10;
  x2 = 535; y2 = 105;
  x3 = 610; y3 = 360;
  x4 = -30; y4 = 400;

  /* calculate X coefficients */
  xa0 = (double)x1;
  xa1 = (double)(x4 - x1) / (double)YSIZE;
  xa2 = (double)(x2 - x1) / (double)XSIZE;
  xa3 = (double)(x1-x2+x3-x4) / (double)XSIZE / (double)YSIZE;

  /* calculate Y coefficients */
  ya0 = (double)y1;
  ya1 = (double)(y4 - y1) / YSIZE;
  ya2 = (double)(y2 - y1) / XSIZE;
  ya3 = (double)(y1-y2+y3-y4) / XSIZE / YSIZE;

  /* allocate input and output memory buffers */
  ir = (unsigned char *) malloc (XSIZE*YSIZE);
  ig = (unsigned char *) malloc (XSIZE*YSIZE);
  ib = (unsigned char *) malloc (XSIZE*YSIZE);
  or = (unsigned char *) malloc (XSIZE*YSIZE);
  og = (unsigned char *) malloc (XSIZE*YSIZE);
  ob = (unsigned char *) malloc (XSIZE*YSIZE);

  /* read the input image */
  read_image("input",ir,ig,ib);

  /* process all lines of image */
  /* note: since this is a forward transform,
       compute location every 1/2 pixel */
  for (IY = 0; IY < YSIZE; IY = IY + 0.5)
    {
    for (IX = 0; IX < XSIZE; IX = IX + 0.5)
      {
      /* compute output pixel location */
      OX = xa0 + (xa1 * IY) + (xa2 * IX) + (xa3 * IX * IY);
      OY = ya0 + (ya1 * IY) + (ya2 * IX) + (ya3 * IX * IY);

      /* reduce input and output pixel locations to integer */
      ix = (long)IX;
      iy = (long)IY;
      ox = (long)(OX + 0.5);
      oy = (long)(OY + 0.5);

      /* select nearest neighbor */
      /* note: bilinear interpolation could also be performed */
      if ((ox >= 0) && (ox < XSIZE) && (oy >= 0) && (oy < YSIZE))
        {
        or[(oy*XSIZE)+ox] = ir[(iy*XSIZE)+ix];
        og[(oy*XSIZE)+ox] = ig[(iy*XSIZE)+ix];
        ob[(oy*XSIZE)+ox] = ib[(iy*XSIZE)+ix];
        }
      }
    }

  /* write the output image */
  write_image("output",or,og,ob);

  /* free memory buffers */
  free(ir);
  free(ig);
  free(ib);
  free(or);
  free(og);
  free(ob);
}