/****************************************************************************
 * Func: fields_warp                                                        *
 *                                                                          *
 * Desc: warps an image base on the fields warping algorithm presented  in  *
 *       Siggraph 1992 proceedings                                          *
 *                                                                          *
 * Params: source- pointer to source image                                  *
 *         ptr- pointer to destination image buffer                        *
 *         source_lines- array of source image control lines                *
 *         dest_lines- array of destination image control lines             *
 *         num_lines- number of control lines                               *
 *         num_frames- number of frames in warp sequence                    *
 *         outfile- base name of output file                                *
 *         rows- number of rows in image                                    *
 *         cols- number of cols in image                                    *
 *         type- image type                                                 *
 ****************************************************************************/

void fields_warp(void *source_ptr, void *dest_ptr, LINE *source_lines, LINE *dest_lines,
	    int num_lines, int num_frames, char *outfile, int rows, int cols,
	    int type)
    {
    int frame;                  /* current frame number */
    int line;                   /* current line number */
    int x, y;                   /* destination coordinate indices */
    int source_x, source_y;     /* source coordinate indices */
    float u, v;                 /* intermediate source coordinate indices */
    char suffix[3];             /* terminating output filename string */
    char file_name[100];        /* name of output file */
    unsigned long source_index; /* index into linear addressed source image */
    unsigned long dest_index;   /* index into linear addressed destination */
    LINE warp_lines[100];       /* array of intermediate control lines */
    float fweight;              /* frame weight used to interpolate lines */
    int last_row, last_col;     /* last valid row and column */
    float fraction, distance;   /* fraction along and distance from control line */
    float fdist;                /* intermediate variable used to compute distance */
    int dx, dy;                 /* delta X and delta Y */
    float weight_sum;           /* accumulated weight */
    double numerator;           /* numerator used to compute line weight */
    double denominator;         /* denominator used to compute line weight */
    float sum_x, sum_y;         /* delta sum X and Y */
    float weight;               /* line weight */
    float a=0.001;              /* small number used to avoid divide by 0 */
    float b=2.0;                /* exponent used in line weight calculation */
    float p=0.75;               /* exponent used in line weight calculation */
    image_ptr sourcei;          /* source pointer for gray-scale image */
    pixel_ptr sourcep;          /* source pointer for color image */
    pixel_ptr destp;            /* destination pointer for color image */
    image_ptr desti;            /* destination pointer for gray-scale */

    last_row = rows-1;          /* initialize last row */
    last_col = cols-1;          /* initialize last column */


    /* setup for a color or gray scale image */
    if(type == 5)
	{
	sourcei = (image_ptr) source_ptr;
	desti = (image_ptr) dest_ptr;
	}
    else
	{
	sourcep = (pixel_ptr) source_ptr;
	destp = (pixel_ptr) dest_ptr;
	}

    for(frame=1; frame<=num_frames; frame++)
	{
	fweight = (float)(frame) / num_frames;

	/* interpolate intermediate control lines */
	for(line=0; line<num_lines; line++)
	    {
	    warp_lines[line].P.x=source_lines[line].P.x +
	      ((dest_lines[line].P.x - source_lines[line].P.x)*fweight);
	    warp_lines[line].P.y=source_lines[line].P.y +
	      ((dest_lines[line].P.y - source_lines[line].P.y)*fweight);
	    warp_lines[line].Q.x=source_lines[line].Q.x +
	      ((dest_lines[line].Q.x - source_lines[line].Q.x)*fweight);
	    warp_lines[line].Q.y=source_lines[line].Q.y +
	      ((dest_lines[line].Q.y - source_lines[line].Q.y)*fweight);
	    warp_lines[line].dx = warp_lines[line].Q.x - warp_lines[line].P.x;
	    warp_lines[line].dy = warp_lines[line].Q.y - warp_lines[line].P.y;
	    warp_lines[line].length_squared =
		(long) (warp_lines[line].dx) * warp_lines[line].dx +
		(long) (warp_lines[line].dy) * warp_lines[line].dy;
	    warp_lines[line].length = sqrt(warp_lines[line].length_squared);
	    }

	for(y=0; y<rows; y++)
	    {
	    printf("processing frame %d: line %d of %d\n", frame, y, rows);
	    for(x=0; x<cols; x++)
		{
		weight_sum = 0.0;
		sum_x = 0.0;
		sum_y = 0.0;
		for (line = 0; line < num_lines; line++)
		    {
		    dx = x - warp_lines[line].P.x;
		    dy = y - warp_lines[line].P.y;
		    fraction = (dx * (long) warp_lines[line].dx +
				dy * (long) warp_lines[line].dy) /
				(float) warp_lines[line].length_squared;
		    fdist = (dy * (long) warp_lines[line].dx -
			     dx * (long) warp_lines[line].dy) /
			     warp_lines[line].length;
		    if (fraction <= 0 )
			distance = sqrt(dx * (long) dx + dy * (long) dy);
		    else if (fraction >= 1)
			{
			dx = x - warp_lines[line].Q.x;
			dy = y - warp_lines[line].Q.y;
			distance = sqrt(dx * (long) dx + dy * (long) dy);
			}
		    else if (fdist >= 0)
			distance =  fdist;
		    else
			distance = -1.0 * fdist;
		    u = source_lines[line].P.x +
			fraction * source_lines[line].dx -
			fdist * source_lines[line].dy /
			source_lines[line].length;
		    v = source_lines[line].P.y +
			fraction * source_lines[line].dy +
			fdist * source_lines[line].dx /
			source_lines[line].length;

		    /* calculate line weight */
		    numerator = pow((double)(warp_lines[line].length),p);
		    denominator = a + distance;
		    weight = pow((numerator / denominator) , b);

		    sum_x += (u - x) * weight;
		    sum_y += (v - y) * weight;
		    weight_sum += weight;
		    }

		source_x = x + sum_x / weight_sum + 0.5;
		source_y = y + sum_y / weight_sum + 0.5;
		CLIP(source_x, 0, last_col);
		CLIP(source_y, 0, last_row);
		source_index = source_x + (unsigned long) cols * source_y;
		dest_index = x + (unsigned long) cols * y;
		if(type == 5)
		    desti[dest_index] = sourcei[source_index];
		else
		    {
		    destp[dest_index].r = sourcep[source_index].r;
		    destp[dest_index].g = sourcep[source_index].g;
		    destp[dest_index].b = sourcep[source_index].b;
		    }
		}
	    }
	suffix[0]='.';
	itoa(frame, &suffix[1], 10);
	suffix[2]=0;
	strcpy(file_name, outfile);
	strcat(file_name, suffix);
	printf("Writing out frame %d\n",frame);
	if(type == 5)
	    write_pnm(desti, file_name, rows, cols, type);
	else
	    write_pnm((image_ptr) destp, file_name, rows, cols, type);
	}
    }