/****************************************************************************
 * Func: fant_resample                                                      *
 *                                                                          *
 * Desc: performs resample of image row or column using Fant's resampling   *
 *       algorithm presented in IEEE Computer Graphics & Applications Jan 86*
 *                                                                          *
 * Params: source- pointer to source image                                  *
 *         len- length of row or column being processed                     *
 *         pixel_offset- number of pixels to next pixe in row or col        *
 *         out_position- output coordinates for each input pixel            *
 *         dest- pointer to destination image                               *
 ****************************************************************************/

void fant_resample(image_ptr source, int len, int pixel_offset,
		   float *out_position, image_ptr dest)
    {
    int i, j;            /* loop variables */
    int current_pixel;   /* pixel being processed */
    int next_pixel;      /* one pixel beyond current_pixel */
    double intensity;    /* linearly interpolated pixel value */
    double sizfac;       /* size factor */
    double inseg;        /* part of current pixel to contribute */
    double outseg;       /* how much required to complete next pixel */
    double accum;        /* accumulator */
    double *in_position; /* normalized scaling data */


    in_position = malloc((len+1) * sizeof(double));
    /* normalize scaling factors to integer pixel boundaries */
    i=0;
    for(j=0; j<len; j++)
	{
	while(out_position[i+1] < j)
	    i++;
	in_position[j] = i + (double) (j - out_position[i]) /
				(out_position[i + 1] - out_position[i]);
	}
    in_position[len]=len;      /* initialize last element */

    /* initialize variables */
    inseg  = 1.0;
    outseg = in_position[1];
    sizfac = outseg;
    accum = 0.0;
    current_pixel = *source;
    source += pixel_offset;
    next_pixel = *source;
    source += pixel_offset;
    i=1;

    while(i<len)
	{
	intensity = inseg * current_pixel + (1.0 - inseg) * next_pixel;
	if(inseg < outseg) /* there is not sufficient current pixel */
	    {		   /* available to compute output pixel */
	    accum += (intensity * inseg);      /* add inseg contribution */
	    outseg -= inseg;             /* inseg portion now satisfied */
	    inseg = 1.0;                 /* reinitialize to only 1 pixel */
	    current_pixel = next_pixel;
	    next_pixel = *source;
	    source += pixel_offset;       /* increment to next pixel */
	    }
	else                /* there is sufficient input pixel left to */
	    {               /* compute a new output pixel */
	    accum += (intensity * outseg);   /* add outseg contribution */
	    accum /= sizfac;                 /* normalize pixel */
	    *dest = (int) MIN(accum, 0xff);  /* insure no overflow */
	    dest += pixel_offset;            /* increment dest pointer */
	    accum = 0.0;                     /* reset accumulator */
	    inseg -= outseg;                 /* outseg portion has been used */
	    outseg = in_position[i+1] - in_position[i];
	    sizfac = outseg;                 /* set outseg */
	    i++;
	    }
	}
    free(in_position);
    }