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