import java.util.Vector;
import Jama.Matrix;

/**
 * Filter Implementation.
 * Easy Filter V1.0, Copyright E.B/JavaZOOM - 1999.
 * Using the Jama Matrix Package.
 *
 * =========================================================================
 * Copyright (C) 1999 by E.B, JavaZOOM All rights reserved.
 *
 * Permission to use, copy, modify, and distribute this software is freely
 * granted by JavaZOOM, provided that the copyright notice
 * above and the following warranty disclaimer are preserved in human
 * readable form.
 *
 * Because this software is licenses free of charge, it is provided
 * "AS IS", with NO WARRANTY.  TO THE EXTENT PERMITTED BY LAW, JavaZOOM
 * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
 * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * JavaZOOM WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
 * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
 * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. 
 *
 * =========================================================================
 */
public class FilterImplementation extends GenericFilter
{
    /**
     * FilterImplementation Constructor.
     * Do not modify !
     */
    public FilterImplementation(Vector paramsValue,int[] imageArray,int imgWidth, int imgHeight)
    {
  		super(paramsValue,imageArray,imgWidth,imgHeight);
    }
	/**
	 * Gradient Computation based on Philippe Paillou article.
	 * IEEE Trans. Geo. & Remote Sensing VOL 35 No.1 Jan. 1997.
	 * It uses sigma,lamda,k (Canny's criteria) through c,alpha,w parameters.
 	 */
	public void computeFilter()
	{
		/*-- Images parameters --*/
		int[]	source = getInitialImageBuffer();
		int		samples = getInitialImageWidth();
		int		lines = getInitialImageHeight();
		
		/*-- Set Output Data --*/		
		setFinalImageWidth(samples);
		setFinalImageHeight(lines);
		int[]	target = new int[samples*lines];
		setFinalImageBuffer(target);		

		/*-- Get GUI Parameters --*/
		double alpha = getParamValue(0);
		double w = getParamValue(1);

		/*-- Filter Red Component only (coz greyscale image) --*/
		int[]	sourceR = new int[lines*samples];
		for (int loc=0;loc<lines*samples;loc++) sourceR[loc] = getRed(source[loc]);

		//int[]	sourceG = new int[lines*samples];
		//int[]	sourceB = new int[lines*samples];
		//for (int loc=0;loc<lines*samples;loc++) sourceG[loc] = getGreen(source[loc]);
		//for (int loc=0;loc<lines*samples;loc++) sourceB[loc] = getBlue(source[loc]);
		//int[] filteredG = reduce1D(rawFilter(sourceG,lines,samples,c,alpha,w),samples,lines);
		//int[] filteredB = reduce1D(rawFilter(sourceB,lines,samples,c,alpha,w),samples,lines);

		int[] filteredR = reduce1D(rawFilter(sourceR,lines,samples,alpha,w),samples,lines);
			
		/*-- Fill in target --*/
		for (int loc=0;loc<lines*samples;loc++) target[loc] = setPixel(filteredR[loc],filteredR[loc],filteredR[loc]);
	}

	/*------------------------------------------*/
	/*--           Start Computation          --*/
	/*------------------------------------------*/
	private double[][] rawFilter(int[] theSource,int N,int P,double alpha,double w)
	{
		/*-- Transform int Array in double Matrix representation --*/
		Matrix I = new Matrix(extend2D(theSource,P,N));	
		double q = w/alpha;
		double sigma = Math.sqrt(2/(alpha*(1-q*q)));
		double lamda = Math.sqrt(2*alpha);
		double k = Math.sqrt((1-q*q)/(5-q*q));
		/*------------------------------------------------------------------------*/
		/*--       Calculous of Ix and Iy by ARMA filtering operations          --*/
		/*------------------------------------------------------------------------*/
		double a1 = 2*Math.exp(-alpha)*cosh(w)-Math.exp(-2*alpha)-1;
		double b1 = -2*Math.exp(-alpha)*cosh(w);
		double b2 = Math.exp(-2*alpha);
		double d = (1-2*Math.exp(-alpha)*cosh(w)+Math.exp(-2*alpha))/(2*alpha*Math.exp(-alpha)*sinh(w)+w*(1-Math.exp(-2*alpha)));
		double c1 = alpha*d;
		double c2 = w*d;
		double a0p = c2;
		double a1p = (c1*sinh(w)-c2*cosh(w))*Math.exp(-alpha);
		double a1m = a1p-c2*b1;
		double a2m = -c2*b2;
		/*--------------------*/
		/*-- Ix Computation --*/
		/*--------------------*/
		double[] MAXp = {0.0,1.0};
		double[] ARXp = {1.0,b1,b2};
		double[] MAXm = {0.0,1.0};
		double[] ARXm = {1.0,b1,b2};
		Matrix Xp = new Matrix(N,P);
		Matrix Xm = new Matrix(N,P);
		Matrix Ix0 = new Matrix(N,P);
		for (int l=0;l<N;l++)
		{
			double[] ligne = (I.getMatrix(l,l,0,P-1)).getRowPackedCopy();
			Matrix locXp = (new Matrix(filterARMA(MAXp,ARXp,ligne),P)).transpose();
			Matrix locXm = (new Matrix(flip(filterARMA(MAXm,ARXm,flip(ligne))),P)).transpose();
			Xp.setMatrix(l,l,0,P-1,locXp);
			Xm.setMatrix(l,l,0,P-1,locXm);
		}
		Ix0 = (Xp.minus(Xm)).times(a1);
		double[] MAIxp = {a0p,a1p};
		double[] ARIxp = {1.0,b1,b2};
		double[] MAIxm = {0.0,a1m,a2m};
		double[] ARIxm = {1.0,b1,b2};
		Matrix Ixp = new Matrix(N,P);
		Matrix Ixm = new Matrix(N,P);
		Matrix Ix = new Matrix(N,P);
		for (int kk=0;kk<P;kk++)
		{
			double[] col = (Ix0.getMatrix(0,N-1,kk,kk)).getColumnPackedCopy();
			Matrix locIxp = new Matrix(filterARMA(MAIxp,ARIxp,col),N);
			Matrix locIxm = new Matrix(flip(filterARMA(MAIxm,ARIxm,flip(col))),N);
			Ixp.setMatrix(0,N-1,kk,kk,locIxp);
			Ixm.setMatrix(0,N-1,kk,kk,locIxm);
		}
		Ix = Ixp.plus(Ixm);
			
		/*----------------------------*/
		/*--     Iy Computation     --*/
		/*----------------------------*/
		double[] MAYp = {0.0, 1.0};
		double[] ARYp = {1.0, b1, b2};
		double[] MAYm = {0.0, 1.0};
		double[] ARYm = {1.0, b1, b2};
		Matrix Yp = new Matrix(N,P);
		Matrix Ym = new Matrix(N,P);
		Matrix Iy0 = new Matrix(N,P);
		for (int kk=0;kk<P;kk++)
		{
			double[] col = (I.getMatrix(0,N-1,kk,kk)).getColumnPackedCopy();
			Matrix locYp = new Matrix(filterARMA(MAYp,ARYp,col),N);
			Matrix locYm = new Matrix(flip(filterARMA(MAYm,ARYm,flip(col))),N);
			Yp.setMatrix(0,N-1,kk,kk,locYp);
			Ym.setMatrix(0,N-1,kk,kk,locYm);
		}
		Iy0 = (Yp.minus(Ym)).times(a1);
		double[] MAIyp = {a0p, a1p};
		double[] ARIyp = {1.0, b1, b2};
		double[] MAIym = {0.0, a1m, a2m};
		double[] ARIym = {1.0, b1, b2};
		Matrix Iyp = new Matrix(N,P);
		Matrix Iym = new Matrix(N,P);
		Matrix Iy = new Matrix(N,P);
		for (int l=0;l<N;l++)
		{
			double[] ligne = (Iy0.getMatrix(l,l,0,P-1)).getRowPackedCopy();
			Matrix locIyp = (new Matrix(filterARMA(MAIyp,ARIyp,ligne),P)).transpose();
			Matrix locIym = (new Matrix(flip(filterARMA(MAIym,ARIym,flip(ligne))),P)).transpose();
			Iyp.setMatrix(l,l,0,P-1,locIyp);
			Iym.setMatrix(l,l,0,P-1,locIym);
		}
		Iy = Iyp.plus(Iym);
		
		/*-- Output computation --*/
		double[][] locIy = Iy.getArray();
		double[][] locIx = Ix.getArray();
		double[][] Ifinal = new double[N][P];
		for (int i=0;i<N;i++)
		{
			for (int j=0;j<P;j++) Ifinal[i][j] = Math.sqrt(locIx[i][j]*locIx[i][j] + locIy[i][j]*locIy[i][j]);
		}

		/*-- Enhance Output --*/
		double Maxi = Ifinal[0][0];
		double Mini = Ifinal[0][0];
		for (int i=0;i<N;i++)
		{
			for (int j=0;j<P;j++)
			{
				if (Ifinal[i][j]>Maxi) Maxi = Ifinal[i][j];
				if (Ifinal[i][j]<Mini) Mini = Ifinal[i][j];
			}
		}

		for (int i=0;i<N;i++)
		{
			for (int j=0;j<P;j++)
			{
				Ifinal[i][j] = (Ifinal[i][j] - Mini)*255.0/(Maxi-Mini);
			}
		}
		return Ifinal;
	}

	/*-------------------------------------------------*/
	/*--   Extended mathematical functions from      --*/
	/*--        com.visualnumerics.javagrande        --*/
	/*-------------------------------------------------*/

	/*
	 *	Returns the value of x with the sign of y.
	 */
	static private double sign(double x, double y)
	{
		double abs_x = ((x < 0) ? -x : x);
		return (y < 0.0) ? -abs_x : abs_x;
	}

	/*
	 *	Evaluate a Chebyschev series
	 */
	static double csevl(double x, double coef[])
	{
		double	b0, b1, b2, twox;
		int		i;
	        b1 = 0.0;
        	b0 = 0.0;
		b2 = 0.0;
	        twox = 2.0*x;
       	 	for (i = coef.length-1;  i >= 0;  i--) {
			b2 = b1;
            		b1 = b0;
            		b0 = twox*b1 - b2 + coef[i];
		}
	        return 0.5*(b0-b2);
	}

	private static final double	SINH_COEF[] = {
		0.1730421940471796,
		0.08759422192276048,
		0.00107947777456713,
		0.00000637484926075,
		0.00000002202366404,
		0.00000000004987940,
		0.00000000000007973,
		0.00000000000000009};
	/**
	 * Returns the inverse (arc) hyperbolic sine of a double.
	 * @param x	A double value.
	 * @return  The arc hyperbolic sine of x.
	 * If x is NaN or less than one, the result is NaN.
	 */
	static public double sinh(double x)
	{
		double	ans;
		double	y = Math.abs(x);
		
		if (Double.isNaN(x)) {
			ans = Double.NaN;
		} else if (Double.isInfinite(y)) {
			return x;
		} else if (y < 2.58096e-08) {
			// 2.58096e-08 = Math.sqrt(6.0*EPSILON_SMALL)
			ans = x;
		} else if (y <= 1.0) {
			ans = x*(1.0+csevl(2.0*x*x-1.0,SINH_COEF));
		} else {
			y = Math.exp(y);
			if (y >= 94906265.62) {
				// 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
				ans = sign(0.5*y,x);
			} else {
				ans = sign(0.5*(y-1.0/y),x);
			}
		}
		return ans;
	}

	/**
	 * Returns the hyperbolic cosine of a double.
	 * @param x	A double value.
	 * @return  The hyperbolic cosine of x.
	 * If x is NaN, the result is NaN.
	 */
	static public double cosh(double x)
	{
		double	ans;
		double	y = Math.exp(Math.abs(x));

		if (Double.isNaN(x)) {
			ans = Double.NaN;
		} else if (Double.isInfinite(x)) {
			ans = x;
		} else if (y < 94906265.62) {
			// 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
			ans = 0.5*(y+1.0/y);
		} else {
			ans = 0.5*y;
		}
		return ans;
	}
}