/***************************************************************************
 * Func: perspective                                                       *
 *                                                                         *
 * Desc: performs a perspective transformation and writes out image file   *
 *                                                                         *
 * Params: buffer - pointer to image in memory                             *
 *         cols - number of columns in image                               *
 *         rows - number of rows in image                                  *
 *         fileout - name of output file                                   *
 *         x0, y0, x1, y1, x2, y2, x3, y3 - destination control points     *
 ***************************************************************************/

void perspective(void *buffer, int cols, int rows, int type, char *fileout,
			double x0, double x1, double x2, double x3,
			double y0, double y1, double y2, double y3)
    {
    double delx1, dely1, delx2, dely2, delx3, dely3; /* deltas */
    double a11, a12, a13;             /* coefficients */
    double a21, a22, a23;             /* coefficients */
    double a31, a32;                  /* coefficients */
    int u,v;                          /* integer coordinates to source image */
    double floatu, floatv;            /* floating point coordinates to source */
    int i,j;                          /* loop variables */
    double x,y;                       /* coordinates of destination image */
    unsigned long sourcebase;         /* one-dimensional index to image */
    image_ptr line_buff;              /* buffer for one line of gray image */
    pixel_ptr pixel_buff;             /* buffer for one line of pixels */
    image_ptr source;                 /* pointer to gray image */
    pixel_ptr pixel_source;           /* pointer to color image */
    double A, B, C, D, E, F, G, H, I; /* intermediate coefficients */
    FILE *fp;                         /* file pointer to output file */

    /* image specific (gray vs. color) customizations */
    switch(type)
	{
	case 5: source = (image_ptr) buffer;
		line_buff = (image_ptr) malloc(cols);
		break;
	case 6: pixel_source = (pixel_ptr) buffer;
		pixel_buff = (pixel_ptr) malloc(cols * sizeof(pixel_ptr));
		break;
	default: printf("Bad type value\n");
		 exit(1);
		 break;
	}

    /* open output file */
    if((fp=fopen(fileout, "wb")) == NULL)
	{
	printf("Unable to open %s for output\n", fileout);
	exit(1);
	}

    /* write out header */
    fprintf(fp, "P%d\n%d %d\n255\n", type, cols, rows);

    delx1=x1-x2; dely1=y1-y2;
    delx2=x3-x2; dely2=y3-y2;
    delx3=x0-x1+x2-x3; dely3=y0-y1+y2-y3;

    if((delx3==0.0) && (dely3==0.0))   /* is mapping affine? */
	{
	printf("Quadrilateral Case \n");
	a11=x1-x0; a21=x2-x1; a31=x0;
	a12=y1-y0; a22=y2-y1; a32=y0;
	a13=0.0; a23=0.0;
	}
    else
	{
	a13=(delx3*dely2-delx2*dely3)/(delx1*dely2-dely1*delx2);
	a23=(delx1*dely3-dely1*delx3)/(delx1*dely2-dely1*delx2);
	a11=x1-x0+a13*x1;
	a21=x3-x0+a23*x3;
	a31=x0;
	a12=y1-y0+a13*y1;
	a22=y3-y0+a23*y3;
	a32=y0;
	}

    A = a22-a32*a23;
    B = a31*a23 - a21;
    C = a21*a32 - a31*a22;
    D = a32*a13 - a12;
    E = a11 - a31*a13;
    F = a31*a12 - a11*a32;
    G = a12*a23 - a22*a13;
    H = a21*a13 - a11*a23;
    I = a11*a22 - a21*a12;

   for(i=0; i<rows; i++)
	{
	y = (double) (i) / (double) (rows);
	for(j=0; j<cols; j++)
	    {
	    x = (double) (j)/cols;
	    floatu = ((A*x + B*y +C)/(G*x + H*y + I)) * cols;
	    floatv = ((D*x + E*y +F)/(G*x + H*y + I)) * rows;
	    u = (int) (floatu + 0.5);     /* nearest neighbor interpolation */
	    v = (int) (floatv + 0.5);     /* nearest neighbor interpolation */
	    if((u>=cols) || (u<0) || (v>=rows) || (v<0))
		sourcebase = 0;
	    else
		sourcebase = (unsigned long) v * cols + u;
	    if(type == 6)
		{
		pixel_buff[j].r = pixel_source[sourcebase].r;
		pixel_buff[j].g = pixel_source[sourcebase].g;
		pixel_buff[j].b = pixel_source[sourcebase].b;
		}
	    else
		line_buff[j] = source[sourcebase];
	    }
	if(type == 6)
	    fwrite(pixel_buff, sizeof(pixel), cols, fp);
	else
	    fwrite(line_buff, 1, cols, fp);
	}
    fclose(fp);
    }