/*--- formatted by Jindent 2.1, (www.c-lab.de/~jindent) ---*/

/*
 * The package jm.cepstrum is part of the Renal Function Project
 * for analysis of dynamic contrast medium evalutions MRT-Images
 * of the Kidneys.
 * 
 * Copyright (C) 1999 / 2000 Jens Martin
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License 2
 * as published by the Free Software Foundation.
 * 
 * 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.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */

package jm.cepstrum;


/**
 * Diese Klasse realisiert die bikubische Interpolation.
 * Sie findet Verwendung zur Skalierung und Verschiebung
 * von Bildern.
 * 
 * @see jm.extension.ExImageCanvas
 * @see jm.cepstrum.CepstrumAnalysis
 * @version  2.0, 14/02/2000
 * @author   Jens Martin
 */
public class BicubicInterpolation {


	/**
	 * Array fuer das Bild
	 */
	private int						image[] = null;


	/**
	 * Breite des Bildes
	 */
	private int						width = 0;


	/**
	 * Hoehe des Bildes
	 */
	private int						height = 0;


	/**
	 * Benoetigte Konstant zur Berechnung der Interpolation.
	 * Der lt. Literatur optimale Koeffizient
	 */
	private static double a = -0.75;


	/**
	 * Hilfsarray fuer Berechnung
	 */
	private int						p[] = null;


	/**
	 * Hilfsvariable fuer Berechnung
	 */
	private double				tt;


	/**
	 * Der Konstruktor. Wird ueberladen, damit die Interpolation
	 * von 8, 16, und 32-Bit-Grauwertbildern moeglich ist.
	 * @param    img    byte-Array der Bildpunkte des Bildes
	 * @param    w      Breite des Bildes
	 * @param    h      Hoehe der Bildes
	 */
	public BicubicInterpolation(byte img[], int w, int h) {
		width = w;
		height = h;
		p = new int[16];
		image = new int[width * height];
		for (int i = 0; i < width * height; i++) {
			image[i] = (int) img[i];
		} 

		// System.arraycopy(img,0,image,0,width*height);
	}


	/**
	 * Der Konstruktor. Wird ueberladen, damit die Interpolation
	 * von 8, 16, und 32-Bit-Grauwertbildern moeglich ist.
	 * @param    img    short-Array der Bildpunkte des Bildes
	 * @param    w      Breite des Bildes
	 * @param    h      Hoehe der Bildes
	 */
	public BicubicInterpolation(short img[], int w, int h) {
		width = w;
		height = h;
		p = new int[16];
		image = new int[width * height];
		for (int i = 0; i < width * height; i++) {
			image[i] = (int) img[i];
		} 

		// System.arraycopy(img,0,image,0,width*height);
	}


	/**
	 * Der Konstruktor. Wird ueberladen, damit die Interpolation
	 * von 8, 16, und 32-Bit-Grauwertbildern moeglich ist.
	 * @param    img    int-Array der Bildpunkte des Bildes
	 * @param    w      Breite des Bildes
	 * @param    h      Hoehe der Bildes
	 */
	public BicubicInterpolation(int img[], int w, int h) {
		width = w;
		height = h;
		p = new int[16];
		image = new int[width * height];

		/*
		 * for (int i = 0; i < width*height; i++)
		 * {
		 * image[i] = (int)img[i];
		 * }
		 */
		System.arraycopy(img, 0, image, 0, width * height);
	}


	/**
	 * Einzige Methode der Klasse. Ermoeglicht das Auslesen von
	 * Grauwerten an beliebigen (Zwischen-)Rasterkoordinaten.
	 * @param x x-Position (reellwertig)
	 * @param y y-Position (reellwertig)
	 * @return Der berechnete Grauwert. Es wird immer ein Integer geliefert.
	 */
	public int getPixel(double x, double y) {
		int result = 0;

		if (x >= (double) width || x < 0.0 || y >= (double) height || y < 0.0) {
			return result;
		} 
		double	ix = x - (double) Math.floor(x);
		double	iy = y - (double) Math.floor(y);
		int			posx = ((int) Math.floor(x)) - 1;
		int			posy = ((int) Math.floor(y)) - 1;
		int			m, n;
		int			i, j;

		for (j = 0; j <= 3; j++) {
			n = posy + j;
			if (n < 0) {
				n = 1;
			} 
			if (n == height) {
				n = height - 2;
			} 
			if (n > height) {
				n = height - 3;
			} 
			for (i = 0; i <= 3; i++) {
				m = posx + i;
				if (m < 0) {
					m = 1;
				} 
				if (m == width) {
					m = width - 2;
				} 
				if (m > width) {
					m = width - 3;

					// p[(j*4)+i] = image[(n*width)+m];
				} 
				p[(j << 2) + i] = image[(n * width) + m];
			} 
		} 
		double	h3, h2, h1, h0;
		double	c3x = c3(ix);
		double	c2x = c2(ix);
		double	c1x = c1(ix);
		double	c0x = c0(ix);

		h3 = (p[0] * c3x) + (p[1] * c2x) + (p[2] * c1x) + (p[3] * c0x);
		h2 = (p[4] * c3x) + (p[5] * c2x) + (p[6] * c1x) + (p[7] * c0x);
		h1 = (p[8] * c3x) + (p[9] * c2x) + (p[10] * c1x) + (p[11] * c0x);
		h0 = (p[12] * c3x) + (p[13] * c2x) + (p[14] * c1x) + (p[15] * c0x);
		result = (int) ((h3 * c3(iy)) + (h2 * c2(iy)) + (h1 * c1(iy)) + (h0 * c0(iy)));
		return result;
	} 


	/**
	 * HIlfsmethode zur schnelleren Berechnung.
	 * @param t double-Wert als Einagbe
	 * @return Das berechnete Ergebnis.
	 */
	private double c0(double t) {

		// return ( (-1.0) * a * Math.pow(t,3.0) ) + ( a * Math.pow(t,2.0) );
		tt = t * t;
		return ((-1.0) * a * t * tt) + (a * tt);
	} 


	/**
	 * HIlfsmethode zur schnelleren Berechnung.
	 * @param t double-Wert als Einagbe
	 * @return Das berechnete Ergebnis.
	 */
	private double c1(double t) {

		// return ( (-1.0) * (a + 2.0) * Math.pow(t,3.0) ) + ( ((2.0 *a) + 3.0) * Math.pow(t,2.0) ) - ( a * t );
		tt = t * t;
		return ((-1.0) * (a + 2.0) * t * tt) + (((2.0 * a) + 3.0) * tt) - (a * t);
	} 


	/**
	 * HIlfsmethode zur schnelleren Berechnung.
	 * @param t double-Wert als Einagbe
	 * @return Das berechnete Ergebnis.
	 */
	private double c2(double t) {

		// return ( (a + 2.0) * Math.pow(t,3.0) ) -  ( (a + 3.0) * Math.pow(t,2.0) ) + 1;
		tt = t * t;
		return ((a + 2.0) * t * tt) - ((a + 3.0) * tt) + 1;
	} 


	/**
	 * HIlfsmethode zur schnelleren Berechnung.
	 * @param t double-Wert als Einagbe
	 * @return Das berechnete Ergebnis.
	 */
	private double c3(double t) {

		// return ( a * Math.pow(t,3.0) ) - ( 2.0 * a * Math.pow(t,2.0) ) + ( a * t );
		tt = t * t;
		return (a * t * tt) - (2.0 * a * tt) + (a * t);
	} 

}


/*--- formatting done in "My Own Convention" style on 08-04-2000 ---*/

