flwt.c

大师写的二代小波经典之作

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: flwt.c,v 1.18 1997/04/03 15:08:04 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/flwt.c,v $
 *  Author: Gabriel Fernandez, Senthil Periaswamy
 *
 *     >>> Fast Lifted Wavelet Transform on the Interval <<<
 *
 *                             .-. 
 *                            / | \
 *                     .-.   / -+- \   .-.
 *                        `-'   |   `-'
 *
 *  Using the "Lifting Scheme" and the "Square 2D" method, the coefficients
 *  of a "Biorthogonal Wavelet" or "Second Generation Wavelet" are found for
 *  a two-dimensional data set. The same scheme is used to apply the inverse
 *  operation.
 */
/* do not edit anything above this line */

/* System header files */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>   /* exit() */
/* FLWT header files */
#include <constant.h>
#include <flwterr.h>
#include <filter.h>
#include <flwt.h>
#include <imgmem.h>
#include <lift.h>
#include <mallat.h>
#include <mem.h>
#include <util.h>

/* Static Global Variables */
static Vector filter = NULL;      /* Vector with filter coefficients */
static Matrix liftingX = NULL;    /* Matrix with lifting coefficients in X */
static Matrix liftingY = NULL;    /* Matrix with lifting coefficients in Y */

/* code */

/* CODE FOR THE FILTER AND LIFTING COEFFICIENTS */

/*
 * GetDual function: Finds the N filter coefficients and allocates a global
 *                   vector with the values to be used in the transform.
 *                   Allocates the global lifting matrices and initialize the
 *                   integral-moments information to calculate the lifting
 *                   coefficients.
 *                   These filter coefficients are used to find the Dual
 *                   Wavelet function.
 */
void
GetDual ( const int width, const int height, const int N, const int nTilde,
          const boolean print, const int levels )
{
    int row, col, noGammas, nX, nY;
    Flt maxN;
    Vector filterPtr;

    /* Calculate Filter Coefficients */
    if ( (filter = GetFilter(N)) == NULL ) {
        Error ("GetDual", NOT_EVEN, RETURN, "N");
        Error ("GetDual", INVALID_RELATIONAL_INT, ABORT, "N", ">", 0);
    } else {
        if (print) {
            int Nrows = (N>>1) + 1;
            filterPtr = filter;
            fprintf (stdout, "\nFilter Coefficients (N=%d)\n", N);
            for (row = 0 ; row < Nrows ; row++) {
                fprintf (stdout, "For %d in the left and %d in the right\n[ ",
                         row, N-row);
                for (col = 0 ; col < N ; col++)
                    fprintf (stdout, "%.18g ", *(filterPtr++));
                fprintf (stdout, "]\n");
            }
        }
    }

    /* Allocate lifting coefficient matrices */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;   /* max vanishing moments */
    noGammas = width >> 1;    /* number of Gammas in X direction */
    if (noGammas > 0) {
        nX = (int)logBaseN ( (Flt)(width-1)/maxN, (Flt)2 );   /* iterations in X */
        nX = (nX < 0) ? 0 : MIN(levels, nX);   /* find lowest level */
        liftingX = matrix( 0, (long)(nX-1), 0, (long)(noGammas*nTilde - 1) );
        for ( row=0 ; row<nX ; row++ )
            for ( col=0 ; col<noGammas*nTilde ; col++ )
                liftingX[row][col] = (Flt)0;
    }
    noGammas = height >> 1;    /* number of Gammas in Y direction */
    if (noGammas > 0) {
        nY = (int)logBaseN ( (Flt)(height-1)/maxN, (Flt)2 );   /* iterations in Y */
        nY = (nY < 0) ? 0 : MIN(levels, nY);   /* find lowest level */
        liftingY = matrix( 0, (long)(nY-1), 0, (long)(noGammas*nTilde - 1) );
        for ( row=0 ; row<nY ; row++ )
            for ( col=0 ; col<noGammas*nTilde ; col++ )
                liftingY[row][col] = (Flt)0;
    }

    /* Initialize Integral-Moments matrices */
    if (!InitMoment (nTilde, width, height, print)) {
        Error ("GetDual", NOT_EVEN, RETURN, "N");
        Error ("GetDual", INVALID_RELATIONAL_INT, ABORT, "N", ">", 0);
    }
}

/*
 * FreeCo function: Frees the memory allocated by the function GetDual().
 */
void
FreeCo ( const int width, const int height, const int N, const int nTilde,
         const int levels )
{
    int noGammas, nRows, nX, nY;
    Flt maxN;

    /* Free Filter Coefficients Vector */
    nRows = (N>>1) + 1;
    free_vector( filter, 0, (long)(nRows*N - 1) );

    /* Free Lifting Coefficient Matrices */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;   /* max vanishing moments */
    if (liftingX) {
        noGammas = width >> 1;
        nX = (int)logBaseN( (Flt)(width-1)/maxN, (Flt)2 );   /* iterations in X */
        nX = (nX < 0) ? 0 : MIN(levels, nX);   /* find lowest level */
        free_matrix( liftingX, 0, (long)(nX-1), 0, (long)(noGammas*nTilde - 1) );
    }
    if (liftingY) {
        noGammas = height >> 1;
        nY = (int)logBaseN ( (Flt)(height-1)/maxN, (Flt)2 );   /* Iterations in Y */
        nY = (nY < 0) ? 0 : MIN(levels, nY);   /* find lowest level */
        free_matrix( liftingY, 0, (long)(nY-1), 0, (long)(noGammas*nTilde - 1) );
    }

    /* Free Integral-Moments Matrices */
    FreeMoment( width, height, nTilde );
}

/*
 * GetReal function: Calculates the moments and lifting coefficients
 *                   needed for the wavelet transform.
 *                   These coefficients preserve the moments of the Real
 *                   Wavelet function.
 */
void
GetReal ( const int width, const int height, const int N, const int nTilde,
          const boolean print, const int levels )
{
    int x, y, k, n, nX, nY, step, max_step, level;
    Flt  maxN;
    Vector liftPtr;
    char buf[BUFSIZ];

    /* Check for filter vector */
    if (!filter) {
        Error ("GetReal", NO_FILTER_VECTOR, ABORT);
    }

    /* Calculate number of iterations n. It is the maximum  */
    /* value such that the following relation is satisfied. */
    /*        (2^n)*(N-1) <= L < (2^(n+1))*(N-1) + 1        */
    /* Where L = max (signal's length in X-Y direction).    */
    /* and N = max (# dual vanish mom & # real vanish mom)  */
    /* Hence, solving for n, we have the following equation */
    /* for all the cases:  n = floor (log_2((L-1)/(N-1))    */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;   /* max vanishing moments */
    /* Iterations in X */
    nX = (width == 1) ? 0 : (int)logBaseN ( (Flt)(width-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN (levels, nX);   /* find lowest level */
    /* Iterations in Y */
    nY = (height == 1) ? 0 : (int)logBaseN ( (Flt)(height-1)/maxN, (Flt)2 );
    nY = (nY < 0) ? 0 : MIN (levels, nY);   /* find lowest level */
    /* Total number of iterations */
    n = MAX(nX, nY);   /* find minimum number of iterations */
    max_step = (n==0) ? 0 : 1<<n;   /* maximum step */

    /* Find lifting coefficients for all levels */
    level = 0;
    for ( step = 1 ; step < max_step ; step <<= 1 ) {
        if (nX-- > 0) {
            /* Get lifting coefficients for rows */
            GetLifting (liftingX[nX], filter, step, N, nTilde, width, FALSE, FALSE);
            if (print) {
                fprintf (stdout, "\nLifting Coefficients for Rows (Level=%d)\n",
                         level);
                liftPtr = liftingX[nX];
                for ( x=0 ; x<(width>>1) ; x++ ) {
                    sprintf (buf, "%.20g", *(liftPtr++));
                    for ( k=1 ; k<nTilde ; k++ )
                        sprintf (buf, "%s,%.20g", buf, *(liftPtr++));
                    fprintf (stdout, "Gamma[%d] (%s)\n", x, buf);
                }
            }
        }

        if (nY-- > 0) {
            /* Get lifting coefficients for columns */
            GetLifting (liftingY[nY], filter, step, N, nTilde, height, TRUE, FALSE);
            if (print) {
                fprintf (stdout,
                         "\nLifting Coefficients for Columns (Level=%d)\n",
                         level);
                liftPtr = liftingY[nY];
                for ( y=0 ; y<(height>>1) ; y++ ) {
                    sprintf (buf, "%.20g", *(liftPtr++));
                    for ( k=1 ; k<nTilde ; k++ )
                        sprintf (buf, "%s,%.20g", buf, *(liftPtr++));
                    fprintf (stdout, "Gamma[%d] (%s)\n", y, buf);
                }
            }
        }
        level++;   /* go to next level */
    }
}

/* CODE FOR THE IN-PLACE FLWT (INTERPOLATION CASE, N EVEN) */

/*
 * FLWT2D function: This is a 2D Direct/Inverse Fast Lifted Wavelet 
 *                  Trasform that uses the "Square 2D" method to find 
 *                  the Gamma coefficients of the wavelet function Phi
 *                  in the Direct case and the Lambda coefficients of the
 *                  original signal in the Inverse case. The "Square 2D"
 *                  method applies the transform first to all the rows
 *                  and then to all the columns. This creates squared
 *                  window transforms for the process. Since the Lifting
 *                  Scheme is used, the final coefficients are organized
 *                  in-place in the original matrix Data[][]. In file
 *                  mallat.c there are functions provided to change this
 *                  organization to the Isotropic format originally
 *                  established by Mallat.
 */
void
FLWT2D ( Matrix Data, const int width, const int height,
                      const int N, const int nTilde,
                      const int levels, const boolean inverse )
{
    int i, k, x, y, /* counters */
        n,          /* maximum number of levels */
        nX, nY,     /* number of iterations in X and Y directions */
        step,       /* step size of current level */
        max_step;   /* maximum and initial step for the direct and */
                    /* inverse transforms, respectively */
    Flt maxN;       /* maximum number of vanishing moments */


    /* Verify for existence of filter and lifting coefficients */
    if (!filter) {
        Error ("FLWT2D", NO_FILTER_VECTOR, ABORT);
    }

    if ( !liftingX && !liftingY ) {
	fprintf( stderr,
	         "\nFLWT2D(): lifting coefficients are not initialized."
		 "\nUse GetDual() and GetReal() before calling this function.\n");
        exit( INPUT_ERROR );
    }

    /* Calculate number of iterations n. It is the maximum  */
    /* value such that the following relation is satisfied. */
    /*        (2^n)*(N-1) <= L < (2^(n+1))*(N-1) + 1        */
    /* Where L = max (signal's length in X-Y direction).    */
    /* and N = max (# dual vanish mom & # real vanish mom)  */
    /* Hence, solving for n, we have the following equation */
    /* for all the cases:  n = floor (log_2((L-1)/(N-1))    */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;   /* max vanishing moments */
    /* Iterations in X */
    nX = (width == 1) ? 0 : (int)logBaseN ( (Flt)(width-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN(levels, nX);   /* find lowest level */
    /* Iterations in Y */
    nY = (height == 1) ? 0 : (int)logBaseN ( (Flt)(height-1)/maxN, (Flt)2 );
    nY = (nY < 0) ? 0 : MIN(levels, nY);   /* find lowest level */
    /* Total number of iterations */
    n = MAX(nX, nY);   /* find minimum number of iterations */


    if ( !inverse ) {    /* Forward transform */
        /* Calculate wavelet coefficients */
        max_step = (n==0) ? 0 : 1<<n;   /* maximum step */
        for ( step=1 ; step<max_step ; step<<=1 ) {
            if (nX-- > 0) {
                /* Apply forward transform to the rows */
                for ( y=0 ; y<height ; y+=step ) {
                    FLWT1D_Predict ( &Data[y][0], (long)width, 1,
                                     (long)step, (long)N );
                    FLWT1D_Update ( &Data[y][0], (long)width, 1,
                                    (long)step, (long)nTilde, liftingX[nX] );
                }
            }