package be.tarsos.dsp.wavelet.lift;

/**
 * <p>
 * Line (with slope) wavelet
 * </p>
 * 
 * <p>
 * The wavelet Lifting Scheme "LineWavelet" wavelet approximates the data set
 * using a LineWavelet with with slope (in contrast to the HaarWavelet wavelet
 * where a LineWavelet has zero slope is used to approximate the data).
 * </p>
 * 
 * <p>
 * The predict stage of the LineWavelet wavelet "predicts" that an odd point
 * will lie midway between its two neighboring even points. That is, that the
 * odd point will lie on a LineWavelet between the two adjacent even points. The
 * difference between this "prediction" and the actual odd value replaces the
 * odd element.
 * </p>
 * 
 * <p>
 * The update stage calculates the average of the odd and even element pairs,
 * although the method is indirect, since the predict phase has over written the
 * odd value.
 * </p>
 * 
 * <h4>
 * Copyright and Use</h4>
 * 
 * <p>
 * You may use this source code without limitation and without fee as long as
 * you include:
 * </p>
 * <blockquote> This software was written and is copyrighted by Ian Kaplan, Bear
 * Products International, www.bearcave.com, 2001. </blockquote>
 * <p>
 * This software is provided "as is", without any warrenty or claim as to its
 * usefulness. Anyone who uses this source code uses it at their own risk. Nor
 * is any support provided by Ian Kaplan and Bear Products International.
 * <p>
 * Please send any bug fixes or suggested source changes to:
 * 
 * <pre>
 *      iank@bearcave.com
 * </pre>
 * 
 * @author Ian Kaplan
 */
public class LineWavelet extends LiftingSchemeBaseWavelet {

	/**
	 * <p>
	 * Calculate an extra "even" value for the LineWavelet wavelet algorithm at
	 * the end of the data series. Here we pretend that the last two values in
	 * the data series are at the x-axis coordinates 0 and 1, respectively. We
	 * then need to calculate the y-axis value at the x-axis coordinate 2. This
	 * point lies on a LineWavelet running through the points at 0 and 1.
	 * </p>
	 * <p>
	 * Given two points, x<sub>1</sub>, y<sub>1</sub> and x<sub>2</sub>,
	 * y<sub>2</sub>, where
	 * </p>
	 * 
	 * <pre>
	 *         x<sub>1</sub> = 0
	 *         x<sub>2</sub> = 1
	 * </pre>
	 * <p>
	 * calculate the point on the LineWavelet at x<sub>3</sub>, y<sub>3</sub>,
	 * where
	 * </p>
	 * 
	 * <pre>
	 *         x<sub>3</sub> = 2
	 * </pre>
	 * <p>
	 * The "two-point equation" for a LineWavelet given x<sub>1</sub>,
	 * y<sub>1</sub> and x<sub>2</sub>, y<sub>2</sub> is
	 * </p>
	 * 
	 * <pre>
	 *      .          y<sub>2</sub> - y<sub>1</sub>
	 *      (y - y<sub>1</sub>) = -------- (x - x<sub>1</sub>)
	 *      .          x<sub>2</sub> - x<sub>1</sub>
	 * </pre>
	 * <p>
	 * Solving for y
	 * </p>
	 * 
	 * <pre>
	 *      .    y<sub>2</sub> - y<sub>1</sub>
	 *      y = -------- (x - x<sub>1</sub>) + y<sub>1</sub>
	 *      .    x<sub>2</sub> - x<sub>1</sub>
	 * </pre>
	 * <p>
	 * Since x<sub>1</sub> = 0 and x<sub>2</sub> = 1
	 * </p>
	 * 
	 * <pre>
	 *      .    y<sub>2</sub> - y<sub>1</sub>
	 *      y = -------- (x - 0) + y<sub>1</sub>
	 *      .    1 - 0
	 * </pre>
	 * <p>
	 * or
	 * </p>
	 * 
	 * <pre>
	 *      y = (y<sub>2</sub> - y<sub>1</sub>)*x + y<sub>1</sub>
	 * </pre>
	 * <p>
	 * We're calculating the value at x<sub>3</sub> = 2, so
	 * </p>
	 * 
	 * <pre>
	 *      y = 2*y<sub>2</sub> - 2*y<sub>1</sub> + y<sub>1</sub>
	 * </pre>
	 * <p>
	 * or
	 * </p>
	 * 
	 * <pre>
	 *      y = 2*y<sub>2</sub> - y<sub>1</sub>
	 * </pre>
	 */
	private float new_y(float y1, float y2) {
		float y = 2 * y2 - y1;
		return y;
	}

	/**
	 * <p>
	 * Predict phase of LineWavelet Lifting Scheme wavelet
	 * </p>
	 * 
	 * <p>
	 * The predict step attempts to "predict" the value of an odd element from
	 * the even elements. The difference between the prediction and the actual
	 * element is stored as a wavelet coefficient.
	 * </p>
	 * <p>
	 * The "predict" step takes place after the split step. The split step will
	 * move the odd elements (b<sub>j</sub>) to the second half of the array,
	 * leaving the even elements (a<sub>i</sub>) in the first half
	 * </p>
	 * 
	 * <pre>
	 *     a<sub>0</sub>, a<sub>1</sub>, a<sub>1</sub>, a<sub>3</sub>, b<sub>0</sub>, b<sub>1</sub>, b<sub>2</sub>, b<sub>2</sub>,
	 * </pre>
	 * <p>
	 * The predict step of the LineWavelet wavelet "predicts" that the odd
	 * element will be on a LineWavelet between two even elements.
	 * </p>
	 * 
	 * <pre>
	 *     b<sub>j+1,i</sub> = b<sub>j,i</sub> - (a<sub>j,i</sub> + a<sub>j,i+1</sub>)/2
	 * </pre>
	 * <p>
	 * Note that when we get to the end of the data series the odd element is
	 * the last element in the data series (remember, wavelet algorithms work on
	 * data series with 2<sup>n</sup> elements). Here we "predict" that the odd
	 * element will be on a LineWavelet that runs through the last two even
	 * elements. This can be calculated by assuming that the last two even
	 * elements are located at x-axis coordinates 0 and 1, respectively. The odd
	 * element will be at 2. The <i>new_y()</i> function is called to do this
	 * simple calculation.
	 * </p>
	 */
	protected void predict(float[] vec, int N, int direction) {
		int half = N >> 1;
		float predictVal;

		for (int i = 0; i < half; i++) {
			int j = i + half;
			if (i < half - 1) {
				predictVal = (vec[i] + vec[i + 1]) / 2;
			} else if (N == 2) {
				predictVal = vec[0];
			} else {
				// calculate the last "odd" prediction
				predictVal = new_y(vec[i - 1], vec[i]);
			}

			if (direction == forward) {
				vec[j] = vec[j] - predictVal;
			} else if (direction == inverse) {
				vec[j] = vec[j] + predictVal;
			} else {
				System.out.println("predictline::predict: bad direction value");
			}
		}
	} // predict

	/**
	 * <p>
	 * The predict phase works on the odd elements in the second half of the
	 * array. The update phase works on the even elements in the first half of
	 * the array. The update phase attempts to preserve the average. After the
	 * update phase is completed the average of the even elements should be
	 * approximately the same as the average of the input data set from the
	 * previous iteration. The result of the update phase becomes the input for
	 * the next iteration.
	 * </p>
	 * <p>
	 * In a HaarWavelet wavelet the average that replaces the even element is
	 * calculated as the average of the even element and its associated odd
	 * element (e.g., its odd neighbor before the split). This is not possible
	 * in the LineWavelet wavelet since the odd element has been replaced by the
	 * difference between the odd element and the mid-point of its two even
	 * neighbors. As a result, the odd element cannot be recovered.
	 * </p>
	 * <p>
	 * The value that is added to the even element to preserve the average is
	 * calculated by the equation shown below. This equation is given in Wim
	 * Sweldens' journal articles and his tutorial (<i>Building Your Own
	 * Wavelets at Home</i>) and in <i>Ripples in Mathematics</i>. A somewhat
	 * more complete derivation of this equation is provided in <i>Ripples in
	 * Mathematics</i> by A. Jensen and A. la Cour-Harbo, Springer, 2001.
	 * </p>
	 * <p>
	 * The equation used to calculate the average is shown below for a given
	 * iteratin <i>i</i>. Note that the predict phase has already completed, so
	 * the odd values belong to iteration <i>i+1</i>.
	 * </p>
	 * 
	 * <pre>
	 *   even<sub>i+1,j</sub> = even<sub>i,j</sub> op (odd<sub>i+1,k-1</sub> + odd<sub>i+1,k</sub>)/4
	 * </pre>
	 * <p>
	 * There is an edge problem here, when i = 0 and k = N/2 (e.g., there is no
	 * k-1 element). We assume that the odd<sub>i+1,k-1</sub> is the same as
	 * odd<sub>k</sub>. So for the first element this becomes
	 * 
	 * <pre>
	 *       (2 * odd<sub>k</sub>)/4
	 * </pre>
	 * <p>
	 * or
	 * </p>
	 * 
	 * <pre>
	 *       odd<sub>k</sub>/2
	 * </pre>
	 */
	protected void update(float[] vec, int N, int direction) {
		int half = N >> 1;

		for (int i = 0; i < half; i++) {
			int j = i + half;
			float val;

			if (i == 0) {
				val = vec[j] / 2.0f;
			} else {
				val = (vec[j - 1] + vec[j]) / 4.0f;
			}
			if (direction == forward) {
				vec[i] = vec[i] + val;
			} else if (direction == inverse) {
				vec[i] = vec[i] - val;
			} else {
				System.out.println("update: bad direction value");
			}
		} // for
	}

} // LineWavelet
