/*************************************************************************
Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

package com.one;

public class CubicSpline
{
	/*************************************************************************
	Ресэмплирование бикубическим сплайном.

		Процедура получает значения  функции  на  сетке  OldWidth*OldHeight  и
	путем интерполяции бикубическим  сплайном  вычисляет  значения  функции  в
	узлах Декартовой сетки размером NewWidth*NewHeight. Новая сетка может быть
	как более, так и менее плотная, чем старая.

	Входные параметры:
		A           - массив значений функции на старой сетке.
					  Нумерация элементов [0..OldHeight-1,
					  0..OldWidth-1]
		OldHeight   - старый размер сетки, OldHeight>1
		OldWidth    - старый размер сетки, OldWidth>1
		NewHeight   - новый размер сетки, NewHeight>1
		NewWidth    - новый размер сетки, NewWidth>1

	Выходные параметры:
		B           - массив значений функции на новой сетке.
					  Нумерация элементов [0..NewHeight-1,
					  0..NewWidth-1]

	  -- ALGLIB routine --
		 15 May, 2007
		 Copyright by Bochkanov Sergey
	*************************************************************************/
	public static float[][] bicubicResampleCartesian(float[][] a, int oldheight, int oldwidth, int newheight, int newwidth)
	{
		float[][] buf;
		float[][] b;
		
		float[] x;
		float[] y;
		float[] c;
		
		int i = 0;
		int j = 0;
		int mw = 0;
		int mh = 0;

		//
		// Prepare
		//
		
		mw = Math.max(oldwidth, newwidth);
		mh = Math.max(oldheight, newheight);
		
		buf = new float[oldheight][newwidth];
		b = new float[newheight][newwidth];
		
		x = new float[Math.max(mw, mh)];
		y = new float[Math.max(mw, mh)];

		//
		// Horizontal interpolation
		//
		
		for(i = 0; i < oldheight; i++)
		{

			//
			// Fill X, Y
			//
			
			for(j = 0; j < oldwidth; j++)
			{
				x[j] = (float)j / (float)(oldwidth - 1);
				y[j] = a[i][j];
			}

			//
			// Interpolate and place result into temporary matrix
			//
			
			c = CubicSpline.buildCubicSpline(x, y, oldwidth, 0, 0.0f, 0, 0.0f);
			
			for(j = 0; j < newwidth; j++)
			{
				buf[i][j] = CubicSpline.splineInterpolation(c, (float)j / (float)(newwidth - 1));
			}
		}

		//
		// Vertical interpolation
		//
		
		for(j = 0; j < newwidth; j++)
		{

			//
			// Fill X, Y
			//
			
			for(i = 0; i < oldheight; i++)
			{
				x[i] = (float)i / (float)(oldheight - 1);
				y[i] = buf[i][j];
			}

			//
			// Interpolate and place result into B
			//
			
			c = CubicSpline.buildCubicSpline(x, y, oldheight, 0, 0.0f, 0, 0.0f);
			
			for(i = 0; i < newheight; i++)
			{
				b[i][j] = CubicSpline.splineInterpolation(c, (float)i / (float)(newheight - 1));
			}
		}
		
		return b;
	}
	
	/*************************************************************************
	Построение таблицы коэффициентов кубического сплайна

	Входные параметры:
		X           -   абсциссы, массив с нумерацией элементов [0..N-1]
		Y           -   значения функции,
						массив с нумерацией элементов [0..N-1]
		N           -   число точек, N>=2
		BoundLType  -   тип граничного условия (левая граница)
		BoundL      -   значение первой (или второй, в зависимости от
						BoundType) производной сплайна на левой границе
		BoundRType  -   тип граничного условия (правая граница)
		BoundR      -   значение первой (или второй, в зависимости от
						BoundType) производной сплайна на правой границе

	Выходные параметры:
		C           -   таблица коэффициентов сплайна для использования в
						подпрограмме SplineInterpolation

	Параметры BoundLType/BoundRType задают тип граничного условия и  принимают
	следующие значения:
		* 0, что соответствует граничному условию типа "сплайн, завершающийся
		  параболой" (при этом значения BoundL/BoundR игнорируются).
		* 1, что соответствует граничному условию для первой производной
		* 2, что соответствует граничному условию для второй производной

	  -- ALGLIB PROJECT --
		 Copyright 23.06.2007 by Bochkanov Sergey
	*************************************************************************/
	public static float[] buildCubicSpline(float[] x, float[] y, int n, int boundltype, float boundl, int boundrtype, float boundr)
	{
		if(n < 2 || boundltype < 0 || boundltype > 2 || boundrtype < 0 || boundrtype > 2)
		{
			throw new IllegalArgumentException();
		}
		
		float[] d;
		int i = 0;
		int tblsize = 0;
		float delta = 0;
		float delta2 = 0;
		float delta3 = 0;

		x = cloneArray(x);
		y = cloneArray(y);

		float[] a1 = new float[n];
		float[] a2 = new float[n];
		float[] a3 = new float[n];
		float[] b = new float[n];

		//
		// Special case:
		// * N=2
		// * parabolic terminated boundary condition on both ends
		//
		
		if(n == 2 & boundltype == 0 & boundrtype == 0)
		{

			//
			// Change task type
			//
			
			boundltype = 2;
			boundl = 0;
			boundrtype = 2;
			boundr = 0;
		}

		//
		// Sort points
		//
		
		heapSortPoints(x, y, n);

		//
		// Left boundary conditions
		//
		
		if(boundltype == 0)
		{
			a1[0] = 0;
			a2[0] = 1;
			a3[0] = 1;
			b[0] = 2 * (y[1] - y[0]) / (x[1] - x[0]);
		}
		
		if(boundltype == 1)
		{
			a1[0] = 0;
			a2[0] = 1;
			a3[0] = 0;
			b[0] = boundl;
		}
		
		if(boundltype == 2)
		{
			a1[0] = 0;
			a2[0] = 2;
			a3[0] = 1;
			b[0] = 3 * (y[1] - y[0]) / (x[1] - x[0]) - 0.5f * boundl * (x[1] - x[0]);
		}

		//
		// Central conditions
		//
		
		for(i = 1; i <= n - 2; i++)
		{
			a1[i] = x[i + 1] - x[i];
			a2[i] = 2 * (x[i + 1] - x[i - 1]);
			a3[i] = x[i] - x[i - 1];
			b[i] = 3 * (y[i] - y[i - 1]) / (x[i] - x[i - 1]) * (x[i + 1] - x[i]) + 3 * (y[i + 1] - y[i]) / (x[i + 1] - x[i]) * (x[i] - x[i - 1]);
		}

		//
		// Right boundary conditions
		//
		
		if(boundrtype == 0)
		{
			a1[n - 1] = 1;
			a2[n - 1] = 1;
			a3[n - 1] = 0;
			b[n - 1] = 2 * (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]);
		}
		
		if(boundrtype == 1)
		{
			a1[n - 1] = 0;
			a2[n - 1] = 1;
			a3[n - 1] = 0;
			b[n - 1] = boundr;
		}
		
		if(boundrtype == 2)
		{
			a1[n - 1] = 1;
			a2[n - 1] = 2;
			a3[n - 1] = 0;
			b[n - 1] = 3 * (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]) + 0.5f * boundr * (x[n - 1] - x[n - 2]);
		}

		//
		// Solve
		//
		
		d = solveTridiagonal(a1, a2, a3, b, n);

		//
		// Now problem is reduced to the cubic Hermite spline
		//
		
		return buildHermiteSpline(x, y, d, n);
	}

	/*************************************************************************
	Построение таблицы коэффициентов сплайна Эрмита

	Входные параметры:
		X           -   абсциссы, массив с нумерацией элементов [0..N-1]
		Y           -   значения функции,
						массив с нумерацией элементов [0..N-1]
		D           -   значения производной,
						массив с нумерацией элементов [0..N-1]
		N           -   число точек, N>=2

	Выходные параметры:
		C           -   таблица коэффициентов сплайна для использования в
						подпрограмме SplineInterpolation

	  -- ALGLIB PROJECT --
		 Copyright 23.06.2007 by Bochkanov Sergey
	*************************************************************************/
	public static float[] buildHermiteSpline(float[] x, float[] y, float[] d, int n)
	{
		if(n < 2)
		{
			throw new IllegalArgumentException();
		}

		int i = 0;
		int tblsize = 0;
		float delta = 0;
		float delta2 = 0;
		float delta3 = 0;

		x = cloneArray(x);
		y = cloneArray(y);
		d = cloneArray(d);

		//
		// Sort points
		//
		
		heapSortDPoints(x, y, d, n);

		//
		// Fill C:
		//  C[0]            -   length(C)
		//  C[1]            -   type(C):
		//                      3 - general cubic spline
		//  C[2]            -   N
		//  C[3]...C[3+N-1] -   x[i], i = 0...N-1
		//  C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
		//
		
		tblsize = 3 + n + (n - 1) * 4;
		
		float[] c = new float[tblsize];
		
		c[0] = tblsize;
		c[1] = 3;
		c[2] = n;
		
		for(i = 0; i < n; i++)
		{
			c[3 + i] = x[i];
		}
		
		for(i = 0; i < n - 1; i++)
		{
			delta = x[i + 1] - x[i];
			delta2 = delta * delta;
			delta3 = delta * delta2;
			
			c[3 + n + 4 * i + 0] = y[i];
			c[3 + n + 4 * i + 1] = d[i];
			c[3 + n + 4 * i + 2] = (3 * (y[i + 1] - y[i]) - 2 * d[i] * delta - d[i + 1] * delta) / delta2;
			c[3 + n + 4 * i + 3] = (2 * (y[i] - y[i + 1]) + d[i] * delta + d[i + 1] * delta) / delta3;
		}
		
		return c;
	}

	/*************************************************************************
	Вычисление интерполирующего сплайна

	Входные параметры:
		C   -   массив коэффициентов, вычисленный подпрограммой для
				построения сплайна.
		X   -   точка, в которой вычисляется значение сплайна

	Результат:
		значение сплайна в точке X

	  -- ALGLIB PROJECT --
		 Copyright 23.06.2007 by Bochkanov Sergey
	*************************************************************************/
	public static float splineInterpolation(float[] c, float x)
	{
		if((int)c[1] != 3)
		{
			throw new IllegalArgumentException();
		}
		
		int n = 0;
		int l = 0;
		int r = 0;
		int m = 0;

		n = (int)c[2];

		//
		// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
		//
		
		l = 3;
		r = 3 + n - 2 + 1;
		
		while(l != r - 1)
		{
			m = (l + r) / 2;
			
			if(c[m] >= x)
			{
				r = m;
			}
			else
			{
				l = m;
			}
		}

		//
		// Interpolation
		//
		
		x = x - c[l];
		m = 3 + n + 4 * (l - 3);
		
		return c[m] + x * (c[m + 1] + x * (c[m + 2] + x * c[m + 3]));
	}

	/*************************************************************************
	Internal subroutine. Heap sort.
	*************************************************************************/
	protected static void heapSortPoints(float[] x, float[] y, int n)
	{
		int i = 0;
		int j = 0;
		int k = 0;
		int t = 0;
		
		float tmp = 0;
		
		boolean isascending = true;
		boolean isdescending = true;

		//
		// Test for already sorted set
		//
		
		for(i = 1; i < n; i++)
		{
			isascending &= x[i] > x[i - 1];
			isdescending &= x[i] < x[i - 1];
		}
		
		if(isascending)
		{
			return;
		}
		
		if(isdescending)
		{
			for(i = 0; i < n; i++)
			{
				j = n - 1 - i;
				
				if(j <= i)
				{
					break;
				}
				
				tmp = x[i];
				x[i] = x[j];
				x[j] = tmp;
				
				tmp = y[i];
				y[i] = y[j];
				y[j] = tmp;
			}
			
			return;
		}

		//
		// Special case: N=1
		//
		
		if(n == 1)
		{
			return;
		}

		//
		// General case
		//
		
		i = 2;
		
		do
		{
			t = i;
			
			while(t != 1)
			{
				k = t / 2;
				
				if(x[k - 1] >= x[t - 1])
				{
					t = 1;
				}
				else
				{
					tmp = x[k - 1];
					x[k - 1] = x[t - 1];
					x[t - 1] = tmp;
					
					tmp = y[k - 1];
					y[k - 1] = y[t - 1];
					y[t - 1] = tmp;
					
					t = k;
				}
			}
		}
		while(i++ <= n);
		
		i = n - 1;
		
		do
		{
			tmp = x[i];
			x[i] = x[0];
			x[0] = tmp;
			
			tmp = y[i];
			y[i] = y[0];
			y[0] = tmp;
			
			t = 1;
			
			while(t != 0)
			{
				k = 2 * t;
				
				if(k > i)
				{
					t = 0;
				}
				else
				{
					if(k < i)
					{
						if(x[k] > x[k - 1])
						{
							k = k + 1;
						}
					}
					
					if(x[t - 1] >= x[k - 1])
					{
						t = 0;
					}
					else
					{
						tmp = x[k - 1];
						x[k - 1] = x[t - 1];
						x[t - 1] = tmp;
						
						tmp = y[k - 1];
						y[k - 1] = y[t - 1];
						y[t - 1] = tmp;
						
						t = k;
					}
				}
			}
		}
		while(i-- >= 1);
	}


	/*************************************************************************
	Internal subroutine. Heap sort.
	*************************************************************************/
	protected static void heapSortDPoints(float[] x, float[] y, float[] d, int n)
	{
		int i = 0;
		int j = 0;
		int k = 0;
		int t = 0;
		
		float tmp = 0;
		
		boolean isascending = true;
		boolean isdescending = true;

		//
		// Test for already sorted set
		//
		
		for(i = 1; i < n; i++)
		{
			isascending &= x[i] > x[i - 1];
			isdescending &= x[i] < x[i - 1];
		}
		
		if(isascending)
		{
			return;
		}
		
		if(isdescending)
		{
			for(i = 0; i < n; i++)
			{
				j = n - 1 - i;
				
				if(j <= i)
				{
					break;
				}
				
				tmp = x[i];
				x[i] = x[j];
				x[j] = tmp;
				
				tmp = y[i];
				y[i] = y[j];
				y[j] = tmp;
				
				tmp = d[i];
				d[i] = d[j];
				d[j] = tmp;
			}
			
			return;
		}

		//
		// Special case: N=1
		//
		
		if(n == 1)
		{
			return;
		}

		//
		// General case
		//
		
		i = 2;
		
		do
		{
			t = i;
			
			while(t != 1)
			{
				k = t / 2;
				
				if(x[k - 1] >= x[t - 1])
				{
					t = 1;
				}
				else
				{
					tmp = x[k - 1];
					x[k - 1] = x[t - 1];
					x[t - 1] = tmp;
					
					tmp = y[k - 1];
					y[k - 1] = y[t - 1];
					y[t - 1] = tmp;
					
					tmp = d[k - 1];
					d[k - 1] = d[t - 1];
					d[t - 1] = tmp;
					
					t = k;
				}
			}
		}
		while(i++ <= n);
		
		i = n - 1;
		
		do
		{
			tmp = x[i];
			x[i] = x[0];
			x[0] = tmp;
			
			tmp = y[i];
			y[i] = y[0];
			y[0] = tmp;
			
			tmp = d[i];
			d[i] = d[0];
			d[0] = tmp;
			
			t = 1;
			
			while(t != 0)
			{
				k = 2 * t;
				
				if(k > i)
				{
					t = 0;
				}
				else
				{
					if(k < i)
					{
						if(x[k] > x[k - 1])
						{
							k = k + 1;
						}
					}
					
					if(x[t - 1] >= x[k - 1])
					{
						t = 0;
					}
					else
					{
						tmp = x[k - 1];
						x[k - 1] = x[t - 1];
						x[t - 1] = tmp;
						
						tmp = y[k - 1];
						y[k - 1] = y[t - 1];
						y[t - 1] = tmp;
						
						tmp = d[k - 1];
						d[k - 1] = d[t - 1];
						d[t - 1] = tmp;
						
						t = k;
					}
				}
			}
		}
		while(i-- >= 1);
	}

	/*************************************************************************
	Internal subroutine. Tridiagonal solver.
	*************************************************************************/
	protected static float[] solveTridiagonal(float[] a, float[] b, float[] c, float[] d, int n)
	{
		int k = 0;
		float t = 0;

		a = cloneArray(a);
		b = cloneArray(b);
		c = cloneArray(c);
		d = cloneArray(d);

		float[] x = new float[n];
		
		a[0] = 0;
		c[n - 1] = 0;
		
		for(k = 1; k <= n - 1; k++)
		{
			t = a[k] / b[k - 1];
			
			b[k] = b[k] - t * c[k - 1];
			d[k] = d[k] - t * d[k - 1];
		}
		
		x[n - 1] = d[n - 1] / b[n - 1];
		
		for(k = n - 2; k >= 0; k--)
		{
			x[k] = (d[k] - c[k] * x[k + 1]) / b[k];
		}
		
		return x;
	}

	/*************************************************************************
	Internal subroutine. Three-point differentiation
	*************************************************************************/
	protected static float diffThreePoint(float t, float x0, float f0, float x1, float f1, float x2, float f2)
	{
		float a = 0;
		float b = 0;

		t = t - x0;
		x1 = x1 - x0;
		x2 = x2 - x0;
		a = (f2 - f0 - x2 / x1 * (f1 - f0)) / (x2 * x2 - x1 * x2);
		b = (f1 - f0 - a * x1 * x1) / x1;
		
		return 2 * a * t + b;
	}
	
	/*************************************************************************
	Internal subroutine. Array cloning
	*************************************************************************/
	protected static float[] cloneArray(float[] array)
	{
		float[] clone = new float[array.length];
		System.arraycopy(array, 0, clone, 0, array.length);
		
		return clone;
	}
}
