/* *  -*- 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. */voidGetDual ( 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(). */voidFreeCo ( 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. */voidGetReal ( 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. */voidFLWT2D ( 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] );                }            }