/*
** Copyright 2005 Huxtable.com. All rights reserved.
*/

package com.jhlabs.image;

/**
 * A class containing static math methods useful for image processing.
 */
public class ImageMath {

	public final static float PI = (float)Math.PI;
	public final static float HALF_PI = (float)Math.PI/2.0f;
	public final static float QUARTER_PI = (float)Math.PI/4.0f;
	public final static float TWO_PI = (float)Math.PI*2.0f;

	/**
	 * Apply a bias to a number in the unit interval, moving numbers towards 0 or 1
	 * according to the bias parameter.
	 * @param a the number to bias
	 * @param b the bias parameter. 0.5 means no change, smaller values bias towards 0, larger towards 1.
	 * @return the output value
	 */
	public static float bias(float a, float b) {
//		return (float)Math.pow(a, Math.log(b) / Math.log(0.5));
		return a/((1.0f/b-2)*(1.0f-a)+1);
	}

	/**
	 * A variant of the gamma function.
	 * @param a the number to apply gain to
	 * @param b the gain parameter. 0.5 means no change, smaller values reduce gain, larger values increase gain.
	 * @return the output value
	 */
	public static float gain(float a, float b) {
/*
		float p = (float)Math.log(1.0 - b) / (float)Math.log(0.5);

		if (a < .001)
			return 0.0f;
		else if (a > .999)
			return 1.0f;
		if (a < 0.5)
			return (float)Math.pow(2 * a, p) / 2;
		else
			return 1.0f - (float)Math.pow(2 * (1. - a), p) / 2;
*/
		float c = (1.0f/b-2.0f) * (1.0f-2.0f*a);
		if (a < 0.5)
			return a/(c+1.0f);
		else
			return (c-a)/(c-1.0f);
	}

	/**
	 * The step function. Returns 0 below a threshold, 1 above.
	 * @param a the threshold position
	 * @param x the input parameter
	 * @return the output value - 0 or 1
	 */
	public static float step(float a, float x) {
		return (x < a) ? 0.0f : 1.0f;
	}

	/**
	 * The pulse function. Returns 1 between two thresholds, 0 outside.
	 * @param a the lower threshold position
	 * @param b the upper threshold position
	 * @param x the input parameter
	 * @return the output value - 0 or 1
	 */
	public static float pulse(float a, float b, float x) {
		return (x < a || x >= b) ? 0.0f : 1.0f;
	}

	/**
	 * A smoothed pulse function. A cubic function is used to smooth the step between two thresholds.
	 * @param a1 the lower threshold position for the start of the pulse
	 * @param a2 the upper threshold position for the start of the pulse
	 * @param b1 the lower threshold position for the end of the pulse
	 * @param b2 the upper threshold position for the end of the pulse
	 * @param x the input parameter
	 * @return the output value
	 */
	public static float smoothPulse(float a1, float a2, float b1, float b2, float x) {
		if (x < a1 || x >= b2)
			return 0;
		if (x >= a2) {
			if (x < b1)
				return 1.0f;
			x = (x - b1) / (b2 - b1);
			return 1.0f - (x*x * (3.0f - 2.0f*x));
		}
		x = (x - a1) / (a2 - a1);
		return x*x * (3.0f - 2.0f*x);
	}

	/**
	 * A smoothed step function. A cubic function is used to smooth the step between two thresholds.
	 * @param a the lower threshold position
	 * @param b the upper threshold position
	 * @param x the input parameter
	 * @return the output value
	 */
	public static float smoothStep(float a, float b, float x) {
		if (x < a)
			return 0;
		if (x >= b)
			return 1;
		x = (x - a) / (b - a);
		return x*x * (3 - 2*x);
	}

	/**
	 * A "circle up" function. Returns y on a unit circle given 1-x. Useful for forming bevels.
	 * @param x the input parameter in the range 0..1
	 * @return the output value
	 */
	public static float circleUp(float x) {
		x = 1-x;
		return (float)Math.sqrt(1-x*x);
	}

	/**
	 * A "circle down" function. Returns 1-y on a unit circle given x. Useful for forming bevels.
	 * @param x the input parameter in the range 0..1
	 * @return the output value
	 */
	public static float circleDown(float x) {
		return 1.0f-(float)Math.sqrt(1-x*x);
	}

	/**
	 * Clamp a value to an interval.
	 * @param a the lower clamp threshold
	 * @param b the upper clamp threshold
	 * @param x the input parameter
	 * @return the clamped value
	 */
	public static float clamp(float x, float a, float b) {
		return (x < a) ? a : (x > b) ? b : x;
	}

	/**
	 * Clamp a value to an interval.
	 * @param a the lower clamp threshold
	 * @param b the upper clamp threshold
	 * @param x the input parameter
	 * @return the clamped value
	 */
	public static int clamp(int x, int a, int b) {
		return (x < a) ? a : (x > b) ? b : x;
	}

	/**
	 * Return a mod b. This differs from the % operator with respect to negative numbers.
	 * @param a the dividend
	 * @param b the divisor
	 * @return a mod b
	 */
	public static double mod(double a, double b) {
		int n = (int)(a/b);
		
		a -= n*b;
		if (a < 0)
			return a + b;
		return a;
	}

	/**
	 * Return a mod b. This differs from the % operator with respect to negative numbers.
	 * @param a the dividend
	 * @param b the divisor
	 * @return a mod b
	 */
	public static float mod(float a, float b) {
		int n = (int)(a/b);
		
		a -= n*b;
		if (a < 0)
			return a + b;
		return a;
	}

	/**
	 * Return a mod b. This differs from the % operator with respect to negative numbers.
	 * @param a the dividend
	 * @param b the divisor
	 * @return a mod b
	 */
	public static int mod(int a, int b) {
		int n = a/b;
		
		a -= n*b;
		if (a < 0)
			return a + b;
		return a;
	}

	/**
	 * The triangle function. Returns a repeating triangle shape in the range 0..1 with wavelength 1.0
	 * @param x the input parameter
	 * @return the output value
	 */
	public static float triangle(float x) {
		float r = mod(x, 1.0f);
		return 2.0f*(r < 0.5 ? r : 1-r);
	}

	/**
	 * Linear interpolation.
	 * @param t the interpolation parameter
	 * @param a the lower interpolation range
	 * @param b the upper interpolation range
	 * @return the interpolated value
	 */
	public static float lerp(float t, float a, float b) {
		return a + t * (b - a);
	}
	
	/**
	 * Linear interpolation.
	 * @param t the interpolation parameter
	 * @param a the lower interpolation range
	 * @param b the upper interpolation range
	 * @return the interpolated value
	 */
	public static int lerp(float t, int a, int b) {
		return (int)(a + t * (b - a));
	}

	/**
	 * Linear interpolation of ARGB values.
	 * @param t the interpolation parameter
	 * @param rgb1 the lower interpolation range
	 * @param rgb2 the upper interpolation range
	 * @return the interpolated value
	 */
	public static int mixColors(float t, int rgb1, int rgb2) {
		int a1 = (rgb1 >> 24) & 0xff;
		int r1 = (rgb1 >> 16) & 0xff;
		int g1 = (rgb1 >> 8) & 0xff;
		int b1 = rgb1 & 0xff;
		int a2 = (rgb2 >> 24) & 0xff;
		int r2 = (rgb2 >> 16) & 0xff;
		int g2 = (rgb2 >> 8) & 0xff;
		int b2 = rgb2 & 0xff;
		a1 = lerp(t, a1, a2);
		r1 = lerp(t, r1, r2);
		g1 = lerp(t, g1, g2);
		b1 = lerp(t, b1, b2);
		return (a1 << 24) | (r1 << 16) | (g1 << 8) | b1;
	}

	/**
	 * Bilinear interpolation of ARGB values.
	 * @param x the X interpolation parameter 0..1
	 * @param y the y interpolation parameter 0..1
	 * @param rgb array of four ARGB values in the order NW, NE, SW, SE
	 * @return the interpolated value
	 */
	public static int bilinearInterpolate(float x, float y, int[] p) {
		float m0, m1;
		int a0 = (p[0] >> 24) & 0xff;
		int r0 = (p[0] >> 16) & 0xff;
		int g0 = (p[0] >> 8) & 0xff;
		int b0 = p[0] & 0xff;
		int a1 = (p[1] >> 24) & 0xff;
		int r1 = (p[1] >> 16) & 0xff;
		int g1 = (p[1] >> 8) & 0xff;
		int b1 = p[1] & 0xff;
		int a2 = (p[2] >> 24) & 0xff;
		int r2 = (p[2] >> 16) & 0xff;
		int g2 = (p[2] >> 8) & 0xff;
		int b2 = p[2] & 0xff;
		int a3 = (p[3] >> 24) & 0xff;
		int r3 = (p[3] >> 16) & 0xff;
		int g3 = (p[3] >> 8) & 0xff;
		int b3 = p[3] & 0xff;

		float cx = 1.0f-x;
		float cy = 1.0f-y;

		m0 = cx * a0 + x * a1;
		m1 = cx * a2 + x * a3;
		int a = (int)(cy * m0 + y * m1);

		m0 = cx * r0 + x * r1;
		m1 = cx * r2 + x * r3;
		int r = (int)(cy * m0 + y * m1);

		m0 = cx * g0 + x * g1;
		m1 = cx * g2 + x * g3;
		int g = (int)(cy * m0 + y * m1);

		m0 = cx * b0 + x * b1;
		m1 = cx * b2 + x * b3;
		int b = (int)(cy * m0 + y * m1);

		return (a << 24) | (r << 16) | (g << 8) | b;
	}

	/**
	 * Return the NTSC gray level of an RGB value.
	 * @param rgb1 the input pixel
	 * @return the gray level (0-255)
	 */
	public static int brightnessNTSC(int rgb) {
		int r = (rgb >> 16) & 0xff;
		int g = (rgb >> 8) & 0xff;
		int b = rgb & 0xff;
		return (int)(r*0.299f + g*0.587f + b*0.114f);
	}
	
	// Catmull-Rom splines
	private final static float m00 = -0.5f;
	private final static float m01 =  1.5f;
	private final static float m02 = -1.5f;
	private final static float m03 =  0.5f;
	private final static float m10 =  1.0f;
	private final static float m11 = -2.5f;
	private final static float m12 =  2.0f;
	private final static float m13 = -0.5f;
	private final static float m20 = -0.5f;
	private final static float m21 =  0.0f;
	private final static float m22 =  0.5f;
	private final static float m23 =  0.0f;
	private final static float m30 =  0.0f;
	private final static float m31 =  1.0f;
	private final static float m32 =  0.0f;
	private final static float m33 =  0.0f;

	/**
	 * Compute a Catmull-Rom spline.
	 * @param x the input parameter
	 * @param numKnots the number of knots in the spline
	 * @param knots the array of knots
	 * @return the spline value
	 */
	public static float spline(float x, int numKnots, float[] knots) {
		int span;
		int numSpans = numKnots - 3;
		float k0, k1, k2, k3;
		float c0, c1, c2, c3;
		
		if (numSpans < 1)
			throw new IllegalArgumentException("Too few knots in spline");
		
		x = clamp(x, 0, 1) * numSpans;
		span = (int)x;
		if (span > numKnots-4)
			span = numKnots-4;
		x -= span;

		k0 = knots[span];
		k1 = knots[span+1];
		k2 = knots[span+2];
		k3 = knots[span+3];
		
		c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
		c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
		c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
		c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
		
		return ((c3*x + c2)*x + c1)*x + c0;
	}
	
	/**
	 * Compute a Catmull-Rom spline, but with variable knot spacing.
	 * @param x the input parameter
	 * @param numKnots the number of knots in the spline
	 * @param xknots the array of knot x values
	 * @param yknots the array of knot y values
	 * @return the spline value
	 */
	public static float spline(float x, int numKnots, int[] xknots, int[] yknots) {
		int span;
		int numSpans = numKnots - 3;
		float k0, k1, k2, k3;
		float c0, c1, c2, c3;
		
		if (numSpans < 1)
			throw new IllegalArgumentException("Too few knots in spline");
		
		for (span = 0; span < numSpans; span++)
			if (xknots[span+1] > x)
				break;
		if (span > numKnots-3)
			span = numKnots-3;
		float t = (float)(x-xknots[span]) / (xknots[span+1]-xknots[span]);
		span--;
		if (span < 0) {
			span = 0;
			t = 0;
		}

		k0 = yknots[span];
		k1 = yknots[span+1];
		k2 = yknots[span+2];
		k3 = yknots[span+3];
		
		c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
		c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
		c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
		c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
		
		return ((c3*t + c2)*t + c1)*t + c0;
	}

	/**
	 * Compute a Catmull-Rom spline for RGB values.
	 * @param x the input parameter
	 * @param numKnots the number of knots in the spline
	 * @param knots the array of knots
	 * @return the spline value
	 */
	public static int colorSpline(float x, int numKnots, int[] knots) {
		int span;
		int numSpans = numKnots - 3;
		float k0, k1, k2, k3;
		float c0, c1, c2, c3;
		
		if (numSpans < 1)
			throw new IllegalArgumentException("Too few knots in spline");
		
		x = clamp(x, 0, 1) * numSpans;
		span = (int)x;
		if (span > numKnots-4)
			span = numKnots-4;
		x -= span;

		int v = 0;
		for (int i = 0; i < 4; i++) {
			int shift = i * 8;
			
			k0 = (knots[span] >> shift) & 0xff;
			k1 = (knots[span+1] >> shift) & 0xff;
			k2 = (knots[span+2] >> shift) & 0xff;
			k3 = (knots[span+3] >> shift) & 0xff;
			
			c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
			c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
			c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
			c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
			int n = (int)(((c3*x + c2)*x + c1)*x + c0);
			if (n < 0)
				n = 0;
			else if (n > 255)
				n = 255;
			v |= n << shift;
		}
		
		return v;
	}

	/**
	 * Compute a Catmull-Rom spline for RGB values, but with variable knot spacing.
	 * @param x the input parameter
	 * @param numKnots the number of knots in the spline
	 * @param xknots the array of knot x values
	 * @param yknots the array of knot y values
	 * @return the spline value
	 */
	public static int colorSpline(int x, int numKnots, int[] xknots, int[] yknots) {
		int span;
		int numSpans = numKnots - 3;
		float k0, k1, k2, k3;
		float c0, c1, c2, c3;
		
		if (numSpans < 1)
			throw new IllegalArgumentException("Too few knots in spline");
		
		for (span = 0; span < numSpans; span++)
			if (xknots[span+1] > x)
				break;
		if (span > numKnots-3)
			span = numKnots-3;
		float t = (float)(x-xknots[span]) / (xknots[span+1]-xknots[span]);
		span--;
		if (span < 0) {
			span = 0;
			t = 0;
		}

		int v = 0;
		for (int i = 0; i < 4; i++) {
			int shift = i * 8;
			
			k0 = (yknots[span] >> shift) & 0xff;
			k1 = (yknots[span+1] >> shift) & 0xff;
			k2 = (yknots[span+2] >> shift) & 0xff;
			k3 = (yknots[span+3] >> shift) & 0xff;
			
			c3 = m00*k0 + m01*k1 + m02*k2 + m03*k3;
			c2 = m10*k0 + m11*k1 + m12*k2 + m13*k3;
			c1 = m20*k0 + m21*k1 + m22*k2 + m23*k3;
			c0 = m30*k0 + m31*k1 + m32*k2 + m33*k3;
			int n = (int)(((c3*t + c2)*t + c1)*t + c0);
			if (n < 0)
				n = 0;
			else if (n > 255)
				n = 255;
			v |= n << shift;
		}
		
		return v;
	}

	/**
	 * An implementation of Fant's resampling algorithm.
	 * @param source the source pixels
	 * @param dest the destination pixels
	 * @param length the length of the scanline to resample
	 * @param offset the start offset into the arrays
	 * @param stride the offset between pixels in consecutive rows
	 * @param out an array of output positions for each pixel
	 */
	public static void resample(int[] source, int[] dest, int length, int offset, int stride, float[] out) {
		int i, j;
		float intensity;
		float sizfac;
		float inSegment;
		float outSegment;
		int a, r, g, b, nextA, nextR, nextG, nextB;
		float aSum, rSum, gSum, bSum;
		float[] in;
		int srcIndex = offset;
		int destIndex = offset;
		int lastIndex = source.length;
		int rgb;

		in = new float[length+1];
		i = 0;
		for (j = 0; j < length; j++) {
			while (out[i+1] < j)
				i++;
			in[j] = i + (float) (j - out[i]) / (out[i + 1] - out[i]);
		}
		in[length] = length;

		inSegment  = 1.0f;
		outSegment = in[1];
		sizfac = outSegment;
		aSum = rSum = gSum = bSum = 0.0f;
		rgb = source[srcIndex];
		a = (rgb >> 24) & 0xff;
		r = (rgb >> 16) & 0xff;
		g = (rgb >> 8) & 0xff;
		b = rgb & 0xff;
		srcIndex += stride;
		rgb = source[srcIndex];
		nextA = (rgb >> 24) & 0xff;
		nextR = (rgb >> 16) & 0xff;
		nextG = (rgb >> 8) & 0xff;
		nextB = rgb & 0xff;
		srcIndex += stride;
		i = 1;

		while (i < length) {
			float aIntensity = inSegment * a + (1.0f - inSegment) * nextA;
			float rIntensity = inSegment * r + (1.0f - inSegment) * nextR;
			float gIntensity = inSegment * g + (1.0f - inSegment) * nextG;
			float bIntensity = inSegment * b + (1.0f - inSegment) * nextB;
			if (inSegment < outSegment) {
				aSum += (aIntensity * inSegment);
				rSum += (rIntensity * inSegment);
				gSum += (gIntensity * inSegment);
				bSum += (bIntensity * inSegment);
				outSegment -= inSegment;
				inSegment = 1.0f;
				a = nextA;
				r = nextR;
				g = nextG;
				b = nextB;
				if (srcIndex < lastIndex)
					rgb = source[srcIndex];
				nextA = (rgb >> 24) & 0xff;
				nextR = (rgb >> 16) & 0xff;
				nextG = (rgb >> 8) & 0xff;
				nextB = rgb & 0xff;
				srcIndex += stride;
			} else {
				aSum += (aIntensity * outSegment);
				rSum += (rIntensity * outSegment);
				gSum += (gIntensity * outSegment);
				bSum += (bIntensity * outSegment);
				dest[destIndex] = 
					((int)Math.min(aSum/sizfac, 255) << 24) |
					((int)Math.min(rSum/sizfac, 255) << 16) |
					((int)Math.min(gSum/sizfac, 255) << 8) |
					(int)Math.min(bSum/sizfac, 255);
				destIndex += stride;
				rSum = gSum = bSum = 0.0f;
				inSegment -= outSegment;
				outSegment = in[i+1] - in[i];
				sizfac = outSegment;
				i++;
			}
		}
	}

}