/****************************************************************************
* $Id: wave_linear.c,v 1.2 2006-11-22 12:58:13+01 rolf Exp rolf $
*****************************************************************************/
/*
  2D linear wavelet
  In place transformation

  1st implementation:

  U H U H U H U H	H = Horizontal details
  V D V D V D V D	V = Vertical details
  U H U H U H U H	D = Diagonal details
  V D V D V D V D	U = Update (low-pass)

  1. Step, horizontal prediction
  	H' = H - (U.left + U.right) / 2
  	D' = D - (V.left + V.right) / 2
  2. Step, horizontal update
  	U' = U + (H'.left + H'.right) / 4
  	V' = V + (D'.left + D'.right) / 4
  3. Step, vertical prediction
  	V'' = V' - (U'.up + U'.down) / 2
  	D'' = D' - (H'.up + H'.down) / 2
  4. Step, vertical prediction
  	U'' = U' + (V''.up + V''.down) / 4
  	H'' = H' + (D''.up + D''.down) / 4

  Limitations:
    image size = A*2^N x B*2^N pixels
    --> perform max N wavelet steps
    * less border conditions
    * simple indexing/stepping

  Changes:
	12.October 2000:
		- Changed /2,/4 to >>1,>>2, because this is equivalent to
		  LMI-Forth opeartor /2
	    - As compiler option: INT-DIV
*/

#include <stdio.h>
#include "defines.h"
#include "wave_linear.h"

#define WAV_INT_DIV		1	// TRUE: -3/4=0, FALSE: -3>>2=-1

#if WAV_INT_DIV
  #if (-1 >> 1) != -1
    #error Machine dependency in Wavelet-Transformation, see comment in source file
  #endif
#endif

typedef void wave_proc( sWord *x, size_t step, size_t count );

/************************************************************************************
  Basic procedures for prediction/update steps
   - Go one line in horizontal, vertical or diagonal (down-right/up-right) direction.
   - The starting point is always on the left/top edge
  Parameters:
    x	  - starting address left/top edge of the image
    step  - size of a single step
    count - number of elements in line ( must be >= 2 !!! )
************************************************************************************/
static void prediction( sWord *x, size_t step, size_t count )
{
  sWord lx=0, rx=0;	/* left/right value, I hope this will be optimized in registers */

  rx =*x;		/* init right value as previous x */
  x += step;		/* point to our data value */
  if( count == 2 ) {	/* special case: take left value (Haar) */
    *x -= rx;
    return;
  }

  while( count > 2 ) {	/* loop (count/2-1) times */
    lx = rx;
    rx = x[step];
#if WAV_INT_DIV
    *x -= (lx + rx) / 2;
#else
    *x -= (lx + rx) >> 1;
#endif
    x += 2*step;
    count -= 2;
  }

  /* right border */
#if WAV_INT_DIV
  *x -= ( 3*rx - lx ) / 2;
#else
  *x -= ( 3*rx - lx ) >> 1;
#endif
}

static void undo_prediction( sWord *x, size_t step, size_t count )
{
  sWord lx=0, rx=0;	/* left/right value, I hope this will be optimized in registers */

  rx =*x;		/* init right value as previous x */
  x += step;		/* point to our data value */
  if( count == 2 ) {	/* special case: take left value (Haar) */
    *x += rx;
    return;
  }

  while( count > 2 ) {	/* loop (count/2-1) times */
    lx = rx;
    rx = x[step];
#if WAV_INT_DIV
    *x += (lx + rx) / 2;
#else
    *x += (lx + rx) >> 1;
#endif
    x += 2*step;
    count -= 2;
  }

  /* right border */
#if WAV_INT_DIV
  *x += ( 3*rx - lx ) / 2;
#else
  *x += ( 3*rx - lx ) >> 1;
#endif
}

static void update( sWord *x, size_t step, size_t count )
{
  sWord lx, rx;		/* left/right value, I hope this will be optimized in registers */

  lx =x[step];		/* init right value as previous x */
  if( count == 2 ) {	/* special case: take left value (Haar) */
#if WAV_INT_DIV
    *x += lx / 2;
#else
    *x += lx >> 1;
#endif
    return;
  }

  /* left border */
  rx = x[3*step];
#if WAV_INT_DIV
  *x += ( 3*lx - rx ) / 4;
#else
  *x += ( 3*lx - rx ) >> 2;
#endif
  x += 2*step;

  while( count > 2 ) {	/* loop (count/2-1) times */
    rx = x[step];
#if WAV_INT_DIV
    *x += (lx + rx) / 4;
#else
    *x += (lx + rx) >> 2;
#endif
    lx = rx;
    x += 2*step;
    count -= 2;
  }
}

static void undo_update( sWord *x, size_t step, size_t count )
{
  sWord lx, rx;		/* left/right value, I hope this will be optimized in registers */

  lx =x[step];		/* init right value as previous x */
  if( count == 2 ) {	/* special case: take left value (Haar) */
#if WAV_INT_DIV
    *x -= lx / 2;
#else
    *x -= lx >> 1;
#endif
    return;
  }

  /* left border */
  rx = x[3*step];
#if WAV_INT_DIV
  *x -= ( 3*lx - rx ) / 4;
#else
  *x -= ( 3*lx - rx ) >> 2;
#endif
  x += 2*step;

  while( count > 2 ) {	/* loop (count/2-1) times */
    rx = x[step];
#if WAV_INT_DIV
    *x -= (lx + rx) / 4;
#else
    *x -= (lx + rx) >> 2;
#endif
    lx = rx;
    x += 2*step;
    count -= 2;
  }
}


/* Basic horizontal/vertical loops */
static void horizontal( wave_proc proc, sWord *x, size_t xStep, size_t xCount, size_t yStep, size_t yCount )
{
	for( ; yCount>0; yCount--, x+=yStep ) {	/* for every row */
		proc( x, xStep, xCount );			/* do prediction/update */
	}
}
static void vertical( wave_proc proc, sWord *x, size_t xStep, size_t xCount, size_t yStep, size_t yCount )
{
	for( ; xCount>0; xCount--, x+=xStep ) {	/* for every column */
		proc( x, yStep, yCount );			/* do prediction/update */
	}
}

/* final (exported) wavelet procedures */
void lin_2d_forward( sWord *data, size_t xSize, size_t ySize, int n )
{
  size_t step;
  size_t xCount, yCount;

  step = 1;
  xCount = xSize, yCount = ySize;

  while( n-- > 0 ) {	/* perform max n wavelet levels */

    if( step < xSize ) {
      horizontal( prediction, data, step, xCount, step*xSize, yCount );
      horizontal( update, data, step, xCount, step*xSize, yCount );
    }

    if( step < ySize ) {
      vertical( prediction, data, step, xCount, step*xSize, yCount );
      vertical( update, data, step, xCount, step*xSize, yCount );
    }
    step *= 2;
    if( xCount > 1 ) xCount /= 2;
    if( yCount > 1 ) yCount /= 2;
  }
}

void lin_2d_reverse( sWord *data, size_t xSize, size_t ySize, int n )
{
  size_t step;
  size_t xCount, yCount;

  step = ( 1 << n );
  xCount = xSize / step;
  yCount = ySize / step;

  while( n-- > 0 ) {	/* perform max n wavelet levels */
    step /= 2;
    if( step < ySize ) yCount *= 2;
    if( step < xSize ) xCount *= 2;

    if( step < ySize ) {
      vertical( undo_update, data, step, xCount, step*xSize, yCount );
      vertical( undo_prediction, data, step, xCount, step*xSize, yCount );
    }
    if( step < xSize ) {
      horizontal( undo_update, data, step, xCount, step*xSize, yCount );
      horizontal( undo_prediction, data, step, xCount, step*xSize, yCount );
    }
  }
}
