package be.tarsos.dsp.wavelet.lift;

/**
*
 

* @author Ian Kaplan
*/
class PolynomialInterpolation {
	
/** number of polynomial interpolation ponts */
	
private final static int numPts = 4;

	
/** Table for 4-point interpolation coefficients */
	
private float fourPointTable[][];

	
/** Table for 2-point interpolation coefficients */
	
private float twoPointTable[][];

	
/**
	 
* <p>
	 
* The polynomial interpolation algorithm assumes that the known points are
	 
* located at x-coordinates 0, 1,.. N-1. An interpolated point is calculated
	 
* at <b><i>x</i></b>, using N coefficients. The polynomial coefficients for
	 
* the point <b><i>x</i></b> can be calculated staticly, using the Lagrange
	 
* method.
	 
* </p>
	 
*
 

	 
* @param x
	 
*
            
the x-coordinate of the interpolated point
	 
* @param N
	 
*
            
the number of polynomial points.
	 
* @param c
	 
*
            
an array for returning the coefficients
	 
*/

	
private void lagrange(float x, int N, float c[]) {
		
float num, denom;

		
for (int i = 0; i < N; i++) {
			
num = 1;
			
denom = 1;
			
for (int k = 0; k < N; k++) {
				
if (i != k) {
					
num = num * (x - k);
					
denom = denom * (i - k);
				
}
			
} // for k
			
c[i] = num / denom;
		
} // for i
	
} // lagrange

	
/**
	 
* <p>
	 
* For a given N-point polynomial interpolation, fill the coefficient table,
	 
* for points 0.5 ... (N-0.5).
	 
* </p>
	 
*/

	
private void fillTable(int N, float table[][]) {
		
float x;
		
float n = N;
		
int i = 0;

		
for (x = 0.5f; x < n; x = x + 1.0f) {
			
lagrange(x, N, table[i]);
			
i++;
		
}
	
} // fillTable

	
/**
	 
* <p>
	 
* PolynomialWavelets constructor
	 
* </p>
	 
* <p>
	 
* Build the 4-point and 2-point polynomial coefficient tables.
	 
* </p>
	 
*/

	
public PolynomialInterpolation() {

		
// Fill in the 4-point polynomial interplation table
		
// for the points 0.5, 1.5, 2.5, 3.5
		
fourPointTable = new float[numPts][numPts];

		
fillTable(numPts, fourPointTable);

		
// Fill in the 2-point polynomial interpolation table
		
// for 0.5 and 1.5
		
twoPointTable = new float[2][2];

		
fillTable(2, twoPointTable);
	
} // PolynomialWavelets constructor

	
/**
	 
* Print an N x N table polynomial coefficient table
	 
*/

	
private void printTable(float table[][], int N) {
		
System.out.println(N + "-point interpolation table:");
		
double x = 0.5;
		
for (int i = 0; i < N; i++) {
			
System.out.print(x + ": ");
			
for (int j = 0; j < N; j++) {
				
System.out.print(table[i][j]);
				
if (j < N - 1)
					
System.out.print(", ");
			
}
			
System.out.println();
			
x = x + 1.0;
		
}
	
}

	
/**
	 
* Print the 4-point and 2-point polynomial coefficient tables.
	 
*/

	
public void printTables() {
		
printTable(fourPointTable, numPts);
		
printTable(twoPointTable, 2);
	
} // printTables

	
/**
	 
* <p>
	 
* For the polynomial interpolation point x-coordinate <b><i>x</i></b>,
	 
* return the associated polynomial interpolation coefficients.
	 
* </p>
	 
*
 

	 
* @param x
	 
*
            
the x-coordinate for the interpolated pont
	 
* @param n
	 
*
            
the number of polynomial interpolation points
	 
* @param c
	 
*
            
an array to return the polynomial coefficients
	 
*/

	
private void getCoef(float x, int n, float c[]) {
		
float table[][] = null;

		
int j = (int) x;
		
if (j < 0 || j >= n) {
			
System.out.println("PolynomialWavelets::getCoef: n = " + n
					
+ ", bad x value");
		
}

		
if (n == numPts) {
			
table = fourPointTable;
		
} else if (n == 2) {
			
table = twoPointTable;
			
c[2] = 0.0f;
			
c[3] = 0.0f;
		
} else {
			
System.out.println("PolynomialWavelets::getCoef: bad value for N");
		
}

		
if (table != null) {
			
for (int i = 0; i < n; i++) {
				
c[i] = table[j][i];
			
}
		
}
	
} // getCoef

	
/**
	 
* <p>
	 
* Given four points at the x,y coordinates {0,d<sub>0</sub>},
	 
* {1,d<sub>1</sub>}, {2,d<sub>2</sub>}, {3,d<sub>3</sub>} return the
	 
* y-coordinate value for the polynomial interpolated point at
	 
* <b><i>x</i></b>.
	 
* </p>
	 
*
 

	 
* @param x
	 
*
            
the x-coordinate for the point to be interpolated
	 
* @param N
	 
*
            
the number of interpolation points
	 
* @param d
	 
*
            
an array containing the y-coordinate values for the known
	 
*
            
points (which are located at x-coordinates 0..N-1).
	 
* @return the y-coordinate value for the polynomial interpolated point at
	 
*
         
<b><i>x</i></b>.
	 
*/

	
public float interpPoint(float x, int N, float d[]) {
		
float c[] = new float[numPts];
		
float point = 0;

		
int n = numPts;
		
if (N < numPts)
			
n = N;

		
getCoef(x, n, c);

		
if (n == numPts) {
			
point = c[0] * d[0] + c[1] * d[1] + c[2] * d[2] + c[3] * d[3];
		
} else if (n == 2) {
			
point = c[0] * d[0] + c[1] * d[1];
		
}

		
return point;
	
} // interpPoint

} // PolynomialInterpolation