lift.c

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

/*
 *  -*- Mode ANSI C -*-
 *  $Id: lift.c,v 1.7 1996/11/20 02:38:51 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/lift.c,v $
 *  Author: Gabriel Fernandez
 *
 *  This module finds the lifting coefficients used in the wavelet
 *  transform to update the lambda coefficients using the contribution
 *  of the neighboring gamma values.
 *  The function GetLifting() updates the values of the (integral,moment) tuples.
 *  Then, calculates the corresponding lifting coefficients.
 */
/* do not edit anyhting above this line */

/* System header files */
#include <assert.h>
#include <math.h>
#include <stdio.h>
/* FLWT header files */
#include <flwterr.h>   /* Error() function */
#include <flwtdef.h>   /* Matrix type and error codes definitions */
#include <lift.h>      /* prototypes */
#include <lu.h>        /* LU decomposition functions */
#include <mem.h>       /* utilities used by the LU decomposition functions */
#include <util.h>      /* ODD() macro */

/* Global static variables */
static Matrix momentX = NULL;   /* Integral and Moments in X direction */
static Matrix momentY = NULL;   /* Integral and Moments in Y direction */

/* Prototype definitions */
static void UpdateMoment ( Matrix, const Vector, const int, const int,
                                   const int, const int, const int );
static void UpdateLifting ( Vector, const Matrix, const int, const int,
                                    const int, const int );

/* code */

/*
 * InitMoment function: Initializes the values of the Integral-Moments
 *                      matrices. The moments are equal to k^i, where
 *                      k is the coefficient index and i is the moment
 *                      number (0,1,2,...).
 *                      There is one matrix for each direction in case
 *                      the idmensions of X and Y are different.
 */
boolean
InitMoment ( const int R, const int sizeX, const int sizeY,
             const boolean print )
{
    int row, col;

    /* Check if R > 0 and even */
    if ( ODD(R) || R<0 )
        return FALSE;

    /* Allocate memory for the Integral-Moment matrices */
    momentX = matrix( 0, (long)(sizeX-1), 0, (long)(R-1) );
    momentY = matrix( 0, (long)(sizeY-1), 0, (long)(R-1) );

    /* Initialize Integral and Moments */
    /* Integral is equal to one at the beginning since all the filter */
    /* coefficients must add to one for each gamma coefficient. */
    /* 1st moment = k, 2nd moment = k^3, 3rd moment = k^3, ... */
    for ( row=0 ; row<sizeY ; row++ )    /* for the cols */
        for ( col=0 ; col<R ; col++ ) {
            if ( row==0 && col==0 )
                momentY[row][col] = (Flt)1.0;   /* pow() domain error */
            else
                momentY[row][col] = (Flt)pow( (Flt)row, (Flt)col );
        }
    for ( row=0 ; row<sizeX ; row++ )    /* for the rows */
        for ( col=0 ; col<R ; col++ ) {
            if ( row==0 && col==0 )
                momentX[row][col] = (Flt)1.0;   /* pow() domain error */
            else
                momentX[row][col] = (Flt)pow( (Flt)row, (Flt)col );
        }

    /* Print Moment Coefficients */
    if (print) {
        fprintf (stdout, "Integral-Moments Coefficients for X (R=%2d)\n", R);
        for ( row=0 ; row<sizeX ; row++ ) {
            fprintf (stdout, "Coeff[%d] (%.20g", row, momentX[row][0]);
            for ( col=1 ; col<R ; col++ )
                fprintf (stdout, ",%.20g", momentX[row][col]);
            fprintf (stdout, ")\n");
        }
        fprintf (stdout, "Integral-Moments Coefficients for Y (R=%2d)\n", R);
        for ( row=0 ; row<sizeY ; row++ ) {
            fprintf (stdout, "Coeff[%d] (%.20g", row, momentY[row][0]);
            for ( col=1 ; col<R ; col++ )
                fprintf (stdout, ",%.20g", momentY[row][col]);
            fprintf (stdout, ")\n");
        }
    }

    return TRUE;
}

/*
 * FreeMoment function: Deallocate memory allocated for the moment matrices.
 */
void
FreeMoment ( const int sizeX, const int sizeY, const int R )
{
    /* Free memory allocated for the moment matrices */
    free_matrix( momentX, 0, (long)(sizeX-1), 0, (long)(R-1) );
    free_matrix( momentY, 0, (long)(sizeY-1), 0, (long)(R-1) );
}

/*
 * GetLifting function: Updates corresponding moment matrix (X if dir is FALSE
 *                      and Y if dir is TRUE) and calculates the lifting
 *                      coefficients using the new moments.
 *                      This function is used for the forward transform and
 *                      uses the original filter coefficients.
 */
void
GetLifting ( Vector lifting, const Vector filter, const int step,
                             const int N, const int R, const int L,
                             const boolean dir, const boolean type )
{
    Matrix moment;   /* pointer to the Integral-Moments matrix */
    int len,         /* number of coefficients in this level */
        step2,       /* step size between same coefficients */
        noGammas;    /* number of Gamma coeffs */

    /*******************************/
    /* Check for boolean variables */
    /*******************************/

    /* Check for the lifting matrix */
    if (!lifting)
        Error ("GetLifting", NO_LIFTING_MATRIX, ABORT);

    /* Select the appropriate moment matrix */
    if (!dir)    /* X direction */
        moment = momentX;
    else         /* Y direction */
        moment = momentY;

    /**********************************/
    /* Initialize important constants */
    /**********************************/

    /* Calculate number of coefficients and step size */
    len = CEIL (L, step);
    step2 = step << 1;

    /* Check for number of Gamma coefficients */
    noGammas  = len >> 1;   /* # of Gammas  = floor(len/2) */

    /******************************/
    /* Decide order of operations */
    /******************************/
    if ( type == FALSE ) {   /* Direct transform */
        UpdateMoment ( moment, filter, step2, noGammas, len, N, R );
        UpdateLifting ( lifting, moment, len, step2, noGammas, R );
    } else {   /* Inverse transform */
        UpdateLifting ( lifting, moment, len, step2, noGammas, R );
        UpdateMoment ( moment, filter, step2, noGammas, len, N, R );
    }
}

/*
 * UpdateMoment function: calculates the integral-moment tuple for the current
 *                        level of calculations.
 */
static void
UpdateMoment ( Matrix moment, const Vector filter, const int step2,
                              const int noGammas, const int len,
                              const int N, const int R )
{
    int i, j,              /* counters */
        row, col,          /* indices of the matrices */
        idxL, idxG,        /* pointers to Lambda & Gamma coeffs, resp. */
        top1, top2, top3;  /* number of iterations for L<R, L=R, L>R */
    Vector filterPtr;      /* pointer to filter coefficients */

    /***************************************/
    /* Update Integral-Moments information */
    /***************************************/

    /* Calculate number of iterations for each case */
    top1 = (N>>1) - 1;                 /* L < R */
    top3 = (N>>1) - ODD(len);          /* L > R */
    top2 = noGammas - (top1 + top3);   /* L = R */

    /* Cases where nbr. left Lambdas < nbr. right Lambdas */
    filterPtr = filter;   /* initialize pointer to first row */
    idxG = step2>>1;      /* second coefficient is Gamma */
    for ( row = 1 ; row <= top1 ; row++ ) {
        idxL = 0;   /* the first Lambda is always the first coefficient */
        filterPtr += N;   /* go to next filter row */
        for ( col = 0 ; col < N ; col++ ) {
            /* Update (int,mom_1,mom_2,...) */
            for ( j = 0 ; j < R ; j++ ) {
                moment[idxL][j] += filterPtr[col]*moment[idxG][j];
            }
            /* Jump to next Lambda coefficient */
            idxL += step2;
        }
        idxG += step2;   /* go to next Gamma coefficient */
    }

    /* Cases where nbr. left Lambdas = nbr. right Lambdas */
    filterPtr += N;   /* go to last filter row */
    for ( i = 0 ; i < top2 ; i++ ) {
        idxL = i*step2;
        for ( col = 0 ; col < N ; col++ ) {
            /* Update (int,mom_1,mom_2,...) */
            for ( j = 0 ; j < R ; j++ ) {
                moment[idxL][j] += filterPtr[col]*moment[idxG][j];
            }
            /* Jump to next Lambda coefficient */
            idxL += step2;
        }
        idxG += step2;   /* go to next Gamma coefficient and stay */
                         /* in the same row of filter coefficients */
    }

    /* Cases where nbr. left Lambdas > nbr. right Lambdas */
    for ( row = top3 ; row >= 1 ; row-- ) {
        idxL = (top2-1)*step2;   /* first Lambda is always in this place */
        filterPtr -= N;   /* go to previous filter row */
        for ( col = N-1 ; col >= 0 ; col-- ) {
            /* Update (int,mom_1,mom_2,...) */
            for ( j = 0 ; j < R ; j++ ) {
                moment[idxL][j] += filterPtr[col]*moment[idxG][j];
            }
            /* Jump to next Lambda coefficient and next filter row */
            idxL += step2;
        }
        idxG += step2;   /* go to next Gamma coefficient */
    }
}

/*
 * UpdateLifting function: calculates the lifting coefficients using the given
 *                         integral-moment tuple.
 */
static void
UpdateLifting ( Vector lifting, const Matrix moment, const int len,
                                const int step2, const int noGammas, const int R )
{
    int lcIdx,             /* index of lifting vector */
        i, j,              /* counters */
        row, col,          /* indices of the matrices */
        idxL, idxG,        /* pointers to Lambda & Gamma coeffs, resp. */
        top1, top2, top3;  /* number of iterations for L<R, L=R, L>R */
    Matrix lift;           /* used to find the lifting coefficients */
    Vector b;              /* used to find the lifting coefficients */
    int *indx;             /* used by the LU routines */
    Flt d;                 /* used by the LU routines */

    /**********************************/
    /* Calculate Lifting Coefficients */
    /**********************************/

    lcIdx = 0;         /* first element of lifting vector */

    /* Allocate space for the indx[] vector */
    indx = ivector( 1, (long)R );

    /* Allocate memory for the temporary lifting matrix and b */
    lift = matrix( 1, (long)R, 1, (long)R );   /* temporary matrix to be solved */
    b = vector( 1, (long)R );                  /* temporary solution vector */

    /* Calculate number of iterations for each case */
    top1 = (R>>1) - 1;                 /* L < R */
    top3 = (R>>1) - ODD(len);          /* L > R */
    top2 = noGammas - (top1 + top3);   /* L = R */

    /* Cases where nbr. left Lambdas < nbr. right Lambdas */
    idxG = step2>>1;   /* second coefficient is Gamma */
    for ( i=0 ; i<top1 ; i++ ) {
        idxL = 0;   /* the first Lambda is always the first coefficient */
        for ( col=1 ; col<=R ; col++ ) {
            /* Load temporary matrix to be solved */
            for ( row=1 ; row<=R ; row++ )
                lift[row][col] = moment[idxL][row-1];      /* matrix */
            /* Jump to next Lambda coefficient */
            idxL += step2;
        }
        /* Apply LU decomposition to lift[][] */
        LUdcmp (lift, R, indx, &d);
        /* Load independent vector */
        for ( j=1 ; j<=R ; j++)
            b[j] = moment[idxG][j-1];   /* independent vector */
        /* Apply back substitution to find lifting coefficients */
        LUbksb (lift, R, indx, b);
        for ( col=1 ; col<=R ; col++ )    /* save them in lifting vector */
            lifting[lcIdx++] = b[col];

        idxG += step2;   /* go to next Gamma coefficient */
    }

    /* Cases where nbr. left Lambdas = nbr. right Lambdas */
    for ( i=0 ; i<top2 ; i++ ) {
        idxL = i*step2;
        for ( col=1 ; col<=R ; col++ ) {
            /* Load temporary matrix to be solved */
            for ( row=1 ; row<=R ; row++ )
                lift[row][col] = moment[idxL][row-1];      /* matrix */
            /* Jump to next Lambda coefficient */
            idxL += step2;
        }
        /* Apply LU decomposition to lift[][] */
        LUdcmp (lift, R, indx, &d);
        /* Load independent vector */
        for ( j=1 ; j<=R ; j++)
            b[j] = moment[idxG][j-1];   /* independent vector */
        /* Apply back substitution to find lifting coefficients */
        LUbksb (lift, R, indx, b);
        for ( col=1 ; col<=R ; col++ )    /* save them in lifting vector */
            lifting[lcIdx++] = b[col];

        idxG += step2;   /* go to next Gamma coefficient */
    }

    /* Cases where nbr. left Lambdas > nbr. right Lambdas */
    for ( i=0 ; i<top3 ; i++ ) {
        idxL = (top2-1)*step2;   /* first Lambda is always in this place */
        for ( col=1 ; col<=R ; col++ ) {
            /* Load temporary matrix to be solved */
            for ( row=1 ; row<=R ; row++ )
                lift[row][col] = moment[idxL][row-1];      /* matrix */
            /* Jump to next Lambda coefficient */
            idxL += step2;
        }
        /* Apply LU decomposition to lift[][] */
        LUdcmp (lift, R, indx, &d);
        /* Load independent vector */
        for ( j=1 ; j<=R ; j++)
            b[j] = moment[idxG][j-1];   /* independent vector */
        /* Apply back substitution to find lifting coefficients */
        LUbksb (lift, R, indx, b);
        for ( col=1 ; col<=R ; col++ )    /* save them in lifting vector */
            lifting[lcIdx++] = b[col];

        idxG += step2;   /* go to next Gamma coefficient */
    }

    /* Free memory */
    free_matrix( lift, 1, (long)R, 1, (long)R );
    free_vector( b, 1, (long)R );
    free_ivector( indx, 1, (long)R );
}