#include <math.h>
#include <string.h>

#include <exec/memory.h>
#include <exec/types.h>

#include <clib/exec_protos.h>

#include <pragmas/exec_pragmas.h>

#include "scale.h"

/* Function protoytypes */
void scale_region( struct PixelRegion *srcPR, struct PixelRegion *destPR, BOOL cubic_interpolation );
double cubic( double dx, int jm1, int j, int jp1, int jp2 );
void pixel_region_get_row( struct PixelRegion *srcPR, int src_row, int orig_width, UBYTE *src, int num_rows );
void pixel_region_set_row( struct PixelRegion *destPR, int dest_row, int width, UBYTE *dest );

void scale_region( struct PixelRegion *srcPR, struct PixelRegion *destPR, BOOL cubic_interpolation )
{
	unsigned char *src_m1, *src, *src_p1, *src_p2;
	unsigned char *s_m1, *s, *s_p1, *s_p2;
	unsigned char *dest, *d;
	double *row, *r;
	int src_row, src_col;
	int bytes, b;
	int width, height;
	int orig_width, orig_height;
	double x_rat, y_rat;
	double x_cum, y_cum;
	double x_last, y_last;
	double *x_frac, y_frac, tot_frac;
	float dx, dy;
	int i, j;
	int frac;
	int advance_dest_x, advance_dest_y;
	int minus_x, plus_x, plus2_x;
	ScaleType scale_type;

	orig_width = srcPR->w;
	orig_height = srcPR->h;

	width = destPR->w;
	height = destPR->h;

	/* Some calculations */
	bytes = destPR->bytes;

	/* The data pointers */
	src_m1 = AllocVec( orig_width * bytes, MEMF_PUBLIC | MEMF_CLEAR );
	src = AllocVec( orig_width * bytes, MEMF_PUBLIC | MEMF_CLEAR );
	src_p1 = AllocVec( orig_width * bytes, MEMF_PUBLIC | MEMF_CLEAR );
	src_p2 = AllocVec( orig_width * bytes, MEMF_PUBLIC | MEMF_CLEAR );
	dest = AllocVec( width * bytes, MEMF_PUBLIC | MEMF_CLEAR );

	/* Find the ratios of old x to new x and old y to new y */
	x_rat = (double)orig_width / (double)width;
	y_rat = (double)orig_height / (double)height;

	/* Determine the scale type */
	if ( x_rat < 1.0 && y_rat < 1.0 )
	{
		scale_type = MagnifyX_MagnifyY;
	}
	else if ( x_rat < 1.0 && y_rat >= 1.0 )
	{
		scale_type = MagnifyX_MinifyY;
	}
	else if ( x_rat >= 1.0 && y_rat < 1.0 )
	{
		scale_type = MinifyX_MagnifyY;
	}
	else
	{
		scale_type = MinifyX_MinifyY;
	}

	/* Allocate an array to help with the calculations */
	row = (double *)AllocVec( sizeof( double ) * width * bytes, MEMF_PUBLIC | MEMF_CLEAR );
	x_frac = (double *)AllocVec( sizeof( double ) * ( width + orig_width ), MEMF_PUBLIC | MEMF_CLEAR );

	/* Initialize the pre-calculated pixel fraction array */
	src_col = 0;
	x_cum = (double)src_col;
	x_last = x_cum;

	for ( i = 0; i < width + orig_width; i++ )
	{
		if ( x_cum + x_rat <= ( src_col + 1 + EPSILON ) )
		{
			x_cum += x_rat;
			x_frac[i] = x_cum - x_last;
		}
		else
		{
			src_col ++;
			x_frac[i] = src_col - x_last;
		}

		x_last += x_frac[i];
	}

	/* Clear the "row" array */
	memset( row, 0, sizeof( double ) * width * bytes );

	/* Counters */
	src_row = 0;
	y_cum = (double)src_row;
	y_last = y_cum;

	/* Get the first src row */
	pixel_region_get_row( srcPR, src_row, orig_width, src, 1 );

	/* Get the next two if possible */
	if ( src_row < ( orig_height - 1 ) )
	{
		pixel_region_get_row( srcPR, ( src_row + 1 ), orig_width, src_p1, 1 );
	}

	if ( ( src_row + 1 ) < ( orig_height - 1 ) )
	{
		pixel_region_get_row( srcPR, ( src_row + 2 ), orig_width, src_p2, 1 );
	}

	/* Scale the selected region */
	i = height;

	while ( i )
	{
		src_col = 0;
		x_cum = (double)src_col;

		/* Determine the fraction of the src pixel we are using for y */
		if ( y_cum + y_rat <= ( src_row + 1 + EPSILON ) )
		{
			y_cum += y_rat;
			dy = y_cum - src_row;
			y_frac = y_cum - y_last;

			advance_dest_y = TRUE;
		}
		else
		{
			y_frac = ( src_row + 1 ) - y_last;
			dy = 1.0;

			advance_dest_y = FALSE;
		}

		y_last += y_frac;

		s = src;
		s_m1 = ( src_row > 0 ) ? src_m1 : src;
		s_p1 = ( src_row < ( orig_height - 1 ) ) ? src_p1 : src;
		s_p2 = ( ( src_row + 1 ) < ( orig_height - 1 ) ) ? src_p2 : s_p1;

		r = row;

		frac = 0;
		j = width;

		while ( j )
		{
			if ( x_cum + x_rat <= ( src_col + 1 + EPSILON ) )
			{
				x_cum += x_rat;
				dx = x_cum - src_col;

				advance_dest_x = TRUE;
			}
			else
			{
				dx = 1.0;

				advance_dest_x = FALSE;
			}

			tot_frac = x_frac[ frac++ ] * y_frac;

			minus_x = ( src_col > 0 ) ? -bytes : 0;
			plus_x = ( src_col < ( orig_width - 1 ) ) ? bytes : 0;
			plus2_x = ( ( src_col + 1 ) < ( orig_width - 1 ) ) ? bytes * 2 : plus_x;

			if ( cubic_interpolation )
			{
				switch ( scale_type )
				{
					case MagnifyX_MagnifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += cubic( dy, (int)cubic( dx, s_m1[ b + minus_x ], s_m1[b], s_m1[ b + plus_x ], s_m1[ b + plus2_x ] ),
								(int)cubic( dx, s[ b + minus_x ], s[b],	s[ b + plus_x ], s[ b + plus2_x ] ),
								(int)cubic( dx, s_p1[ b + minus_x ], s_p1[b], s_p1[ b + plus_x ], s_p1[ b + plus2_x ] ),
								(int)cubic( dx, s_p2[ b + minus_x ], s_p2[b], s_p2[ b + plus_x ], s_p2[ b + plus2_x ] ) ) * tot_frac;
						}
					break;

					case MagnifyX_MinifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += cubic( dx, s[ b + minus_x ], s[b], s[ b + plus_x ], s[ b + plus2_x ] ) * tot_frac;
						}
					break;

					case MinifyX_MagnifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += cubic( dy, s_m1[b], s[b], s_p1[b], s_p2[b] ) * tot_frac;
						}
					break;

					case MinifyX_MinifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += s[b] * tot_frac;
						}
					break;
				}
			}
			else
			{
				switch ( scale_type )
				{
					case MagnifyX_MagnifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += ( ( 1 - dy ) * ( ( 1 - dx ) * s[b] + dx * s[ b + plus_x ] ) +
							dy  * ( ( 1 - dx ) * s_p1[b] + dx * s_p1[ b + plus_x ] ) ) * tot_frac;
						}
					break;

					case MagnifyX_MinifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += ( s[b] * ( 1 - dx ) + s[ b + plus_x ] * dx ) * tot_frac;
						}
					break;

					case MinifyX_MagnifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += ( s[b] * ( 1 - dy ) + s_p1[b] * dy ) * tot_frac;
						}
					break;

					case MinifyX_MinifyY:
						for ( b = 0; b < bytes; b++ )
						{
							r[b] += s[b] * tot_frac;
						}
					break;
				}
			}

			if ( advance_dest_x )
			{
				r += bytes;
				j--;
			}
			else
			{
				s_m1 += bytes;
				s += bytes;
				s_p1 += bytes;
				s_p2 += bytes;
				src_col++;
			}
		}

		if ( advance_dest_y )
		{
			tot_frac = 1.0 / ( x_rat * y_rat );

			/* Copy "row" to "dest" */
			d = dest;
			r = row;
			j = width;

			while ( j-- )
			{
				b = bytes;

				while (b--)
				{
					*d++ = (unsigned char)( *r++ * tot_frac );
				}
			}

			/* Set the pixel region span */
			pixel_region_set_row( destPR, ( height - i ), width, dest );

			/* Clear the "row" array */
			memset( row, 0, sizeof( double ) * width * bytes );

			i--;
		}
		else
		{
			/* Shuffle pointers */
			s = src_m1;
			src_m1 = src;
			src = src_p1;
			src_p1 = src_p2;
			src_p2 = s;
			src_row++;

			if ( ( src_row + 1 ) < ( orig_height - 1 ) )
			{
				pixel_region_get_row( srcPR, ( src_row + 2 ), orig_width, src_p2, 1 );
			}
		}
	}

	/* Free up temporary arrays */
  FreeVec( row );
  FreeVec( x_frac );
  FreeVec( src_m1 );
  FreeVec( src );
  FreeVec( src_p1 );
  FreeVec( src_p2 );
  FreeVec( dest );
}

double cubic( double dx, int jm1, int j, int jp1, int jp2 )
{
	double dx1, dx2, dx3;
	double h1, h2, h3, h4;
	double result;

	/* Constraint parameter = -1 */
	dx1 = fabs( dx );
	dx2 = dx1 * dx1;
	dx3 = dx2 * dx1;
	h1 = dx3 - 2 * dx2 + 1;
	result = h1 * j;

	dx1 = fabs( dx - 1.0 );
	dx2 = dx1 * dx1;
	dx3 = dx2 * dx1;
	h2 = dx3 - 2 * dx2 + 1;
	result += h2 * jp1;

	dx1 = fabs( dx - 2.0 );
	dx2 = dx1 * dx1;
	dx3 = dx2 * dx1;
	h3 = -dx3 + 5 * dx2 - 8 * dx1 + 4;
	result += h3 * jp2;

	dx1 = fabs( dx + 1.0 );
	dx2 = dx1 * dx1;
	dx3 = dx2 * dx1;
	h4 = -dx3 + 5 * dx2 - 8 * dx1 + 4;
	result += h4 * jm1;

	if ( result < 0.0 ) result = 0.0;
  if ( result > 255.0 ) result = 255.0;

  return result;
}

void pixel_region_get_row( struct PixelRegion *srcPR, int src_row, int orig_width, UBYTE *src, int num_rows )
{
	CopyMem( &srcPR->buf[ src_row * orig_width * srcPR->bytes ], src, orig_width * srcPR->bytes );
}

void pixel_region_set_row( struct PixelRegion *destPR, int dest_row, int width, UBYTE *dest )
{
	CopyMem( dest, &destPR->buf[ dest_row * width * destPR->bytes ], width * destPR->bytes );
}