/***************************************************************************
 * Func: bilinear                                                          *
 *                                                                         *
 * Desc: performs a bilinear transformation and writes out image file      *
 *       This routine performs a the more complex map to transform a       *
 *       rectangle to a quadrilateral                                      *
 *                                                                         *
 * 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 bilinear(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 a11, a21, a31, a41;        /* coefficients */
    double a12, a22, a32, a42;        /* more coefficients */
    double first, second;             /* terms used in quadratic equation */
    double denom;                     /* denominator of quadratic equation */
    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;                   /* intermediate coefficients */
    FILE *fp;                         /* file pointer to output file */

    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);


    /* compute coefficients */
    a11=x0-x1-x3+x2;
    a21=x1-x0;
    a31=x3-x0;
    a41=x0;
    a12=y0-y1-y3+y2;
    a22=y1-y0;
    a32=y3-y0;
    a42=y0;

    A=a11*a22-a21*a12;

    if(A==0.0)          /* insure that A never equals 0 */
	A=0.00001;

    for(i=0; i<rows; i++)
	{
	y = (double) (i)/ (double) (rows);
	for(j=0; j<cols; j++)
	    {
	    x = (double) (j)  / (double) cols;
	    B=a12*x-a11*y+a11*a42-a41*a12+a31*a22-a21*a32;
	    C=a32*x-a31*y+a31*a42-a41*a32;
	    first = B*B;
	    second = 4.0*A*C;
	    if(first < second)
		floatu = B/(-2.0*A);     /* ignore 4AC term */
	    else
		floatu = (-1.0*B + sqrt(first-second))/(2.0*A);
	    if((floatu < 0.0) || (floatu > 1.0))   /* check if out of bounds */
		if(first <second)
		    floatu = B/-2.0*A;   /* ignore 4AC term */
		else
		    floatu = (-1.0*B - sqrt(first-second))/(2.0*A);
	    denom = (a11 * floatu + a31);
	    if(denom == 0.0)             /* avoid divide by 0 */
		denom = 0.000001;
	    floatv = (x-a21*floatu-a41)/ denom;
	    u = (int) (floatu * cols + 0.5);
	    v = (int) (floatv * rows + 0.5);
	    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);
    }