flwtpack.c

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

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: flwtpack.c,v 1.7 1996/10/14 21:03:35 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/flwtpack.c,v $
 *  Author: Gabriel Fernandez
 *
 *  Using the "Lifting Scheme" and the "Square 2D" method, the coefficients
 *  of "Biorthogonal Wavelet Packets" are found for one and two-dimensional
 *  data sets. 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> */
/* FLWT header files */
#include <cost.h>
#include <data.h>
#include <flwterr.h>
#include <flwt.h>
#include <flwtpack.h>
#include <imgmem.h>
#include <mallat.h>
#include <mem.h>
#include <memchk.h>
#include <util.h>

/* Prototypes to private functions */
static Matrix
GetMatrix ( Hedge *__graph, const int __s, const int __x, const int __y );
static void
ReconstructPackets ( Hedge *__graph, const int __x, const int __y,
                                     const int __N, const int __nTilde,
                                     const int __l );
static Vector
GetVector ( Hedge *__graph, const int __s, const int __x );
static void
ReconstructPackets1D ( Hedge *__graph, const int __x,
                                       const int __N, const int __nTilde,
                                       const int __l );

/* code */

/*
 * DFLWTPackets function: for the forward case, it sets some initial 
 *                        variables needed for the best basis function,
 *                        finds the hedge corresponding to the best
 *                        basis description of the given data, and
 *                        cleans up any left overs. For the inverse case,
 *                        it takes the given array and finds the
 *                        corresponding hedge. Then, it builds the
 *                        original data set using the inverse calculations.
 */
int *
DFLWTPackets ( Matrix Data, const int sizeX, const int sizeY,
                            const int N, const int nTilde, const int maxLevel,
                            int *blocks )
{
    Hedge *graph;        /* hedge with the best basis of Data */
    register int i;      /* counter */
    int *levels,         /* levels in which the regions are located */
        *outlevels,      /* levels for output */
        nX, nY,          /* number of levels in X and Y */
        L,               /* maximum level in decomposition */
        XY;              /* size of contents vector in hedge */
    Flt maxN,            /* maximum number of vanishing moments */
        cost;            /* cost of the best basis tree */
    Vector dataPtr;      /* pointer to the output matrix */
    Matrix bestCoeffs;   /* matrix that will store the best coefficients */

    /* Calculate maximum level */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;
    nX = (sizeX == 1) ? 0 : (int)logBaseN( (Flt)(sizeX-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN(maxLevel, nX);
    nY = (sizeY == 1) ? 0 : (int)logBaseN( (Flt)(sizeY-1)/maxN, (Flt)2 );
    nY = (nY < 0) ? 0 : MIN(maxLevel, nY);
    L = MAX(nX, nY);

    /* Allocate initial hedge */
    *blocks = (1<<2*L) + 1;
    bestCoeffs = (Matrix)MEM_Init2D( (long)ceil((Flt)(sizeX*sizeY)/(Flt)*blocks),
                                     (long)*blocks );
    levels = ivector( 0, (long)(*blocks-1) );
    graph = (Hedge *)makehedge( *blocks, bestCoeffs, levels, NULL );
    graph->blocks = 0;
    
    /* Find best basis */
    cost = bbwp2( graph, Data, sizeX, sizeY, N, nTilde, 0, L );
    
    /* Put results in final image buffer */
    /*
    fprintf( stdout, "\nCost best basis = %g\t\t", cost );
    fflush (stdout);
    */
    XY = sizeX*sizeY;
    outlevels = ivector( 0, (long)(graph->blocks-1) );
    *blocks = graph->blocks;
    dataPtr = &Data[0][0];
    for ( i=0; i<graph->blocks ; i++ ) {
        register int j;
        int length = XY/(1<<(2*graph->levels[i]));

        outlevels[i] = graph->levels[i];
        for ( j=0 ; j<length ; j++, dataPtr++ ) {
            *dataPtr = graph->contents[i][j];
        }
    }

    /* Clean up */
    graph = freehedge( graph );

    return outlevels;
}


/*
 * IFLWTPackets function: performs the inverse operation to the packet
 *                        transform in order to find the original signal.
 */
void
IFLWTPackets ( Matrix Data, const int sizeX, const int sizeY,
               const int blocks, int *levels, const int N, const int nTilde )
{
    Hedge *graph;
    Matrix content;
    register int i;
    int column, XY, length, maxl;

    /* Allocate a hedge for the manipulation of the basis */
    content = (Matrix)malloc( (size_t)(blocks * sizeof( Vector )) );   /* row of new pointers */
    graph   = (Hedge *)makehedge( blocks, content, levels, NULL );

    /* Extract hedge of coefficients */
    column = 0;
    maxl = -999;
    XY = sizeX * sizeY;
    for ( i=0 ; i<graph->blocks ; i++ ) {
        maxl = MAX(maxl, graph->levels[i]);      /* find maximum level */
        length = XY/(1<<(2*graph->levels[i]));   /* length of vector block */
        graph->contents[i] = Data[0] + column;   /* pointer to beginning of vector block */
        column += length;                        /* place for next column pointer */
    }

    /* Reconstruct original image */
    if ( maxl > 0 ) {
        Vector dataPtr;

        ReconstructPackets( graph, sizeX, sizeY, N, nTilde, maxl );
        dataPtr = &Data[0][0];
        for ( i=0; i<graph->blocks ; i++ ) {
            register int j;
            int length = XY/(1<<(2*graph->levels[i]));

            for ( j=0 ; j<length ; j++, dataPtr++ ) {
                *dataPtr = graph->contents[i][j];
            }
        }
    }

    /* Clean up */
    free( (Matrix)content );
    free( (Hedge *)graph );
}

/*
 * DFLWTPackets1D function: for the forward case, it sets some initial 
 *                          variables needed for the best basis function, 
 *                          finds the hedge corresponding to the best 
 *                          basis description of the given data, and 
 *                          cleans up any left overs. For the inverse 
 *                          case, it takes the given array and finds the 
 *                          corresponding hedge. Then, it builds the
 *                          original data set using the inverse
 *                          calculations.
 */
int *
DFLWTPackets1D ( Vector Data, const int sizeX,
                              const int N, const int nTilde, 
                              const int maxLevel, int *blocks )
{
    Hedge *graph;        /* hedge with the best basis of Data */
    register int i;      /* counter */
    int *levels,         /* levels in which the regions are located */
        *outlevels,      /* levels for output */
        nX;              /* number of levels in X */
    Flt maxN,            /* maximum number of vanishing moments */
        cost;            /* cost of the best basis tree */
    Vector dataPtr;      /* pointer to the output matrix */
    Matrix bestCoeffs;   /* matrix that will store the best coefficients */

    /* Calculate maximum level */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;
    nX = (sizeX == 1) ? 0 : (int)logBaseN( (Flt)(sizeX-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN(maxLevel, nX);

    /* Allocate initial hedge */
    *blocks = (1<<nX) + 1;
    bestCoeffs = (Matrix)MEM_Init2D( (long)ceil(sizeX/(Flt)*blocks),
                                     (long)*blocks );
    levels = ivector( 0, (long)(*blocks-1) );
    graph = (Hedge *)makehedge( *blocks, bestCoeffs, levels, NULL );
    graph->blocks = 0;
    
    /* Find best basis */
    cost = bbwp1( graph, Data, sizeX, N, nTilde, 0, nX );
    
    /* Put results in final image buffer */
    outlevels = ivector( 0, (long)(graph->blocks-1) );
    *blocks = graph->blocks;
    dataPtr = Data;
    for ( i=0; i<graph->blocks ; i++ ) {
        register int j;
        int length = sizeX/(1<<(graph->levels[i]));

        outlevels[i] = graph->levels[i];
        for ( j=0 ; j<length ; j++, dataPtr++ )
            *dataPtr = graph->contents[i][j];
    }

    /* Clean up */
    graph = freehedge( graph );

    return outlevels;
}

/*
 * IFLWTPackets1D function: performs the inverse operation to the packet
 *                        transform in order to find the original signal.
 */
void
IFLWTPackets1D ( Vector Data, const int sizeX,
                 const int blocks, int *levels, const int N, 
                 const int nTilde )
{
    Hedge *graph;
    Matrix content;
    register int i;
    int column, length, maxl;

    /* Allocate a hedge for the manipulation of the basis */
    content = (Matrix)malloc( (size_t)(blocks * sizeof( Vector )) );
    /* row of new pointers */
    graph   = (Hedge *)makehedge( blocks, content, levels, NULL );

    /* Extract hedge of coefficients */
    column = 0;
    maxl = -999;
    for ( i=0 ; i<graph->blocks ; i++ ) {
        maxl = MAX(maxl, graph->levels[i]);      /* find maximum level */
        length = sizeX/(1<<(graph->levels[i]));  /* length of vector block */
        graph->contents[i] = Data + column;      /* pointer to beginning of
                                                    vector block */
        column += length;                        /* place for next column 
                                                    pointer */
    }

    /* Reconstruct original image */
    if ( maxl > 0 ) {
        Vector dataPtr;

        ReconstructPackets1D( graph, sizeX, N, nTilde, maxl );
        dataPtr = Data;
        for ( i=0; i<graph->blocks ; i++ ) {
            register int j;
            int length = sizeX/(1<<(graph->levels[i]));

            for ( j=0 ; j<length ; j++, dataPtr++ ) {
                *dataPtr = graph->contents[i][j];
            }
        }
    }

    /* Clean up */
    free( (Matrix)content );
    free( (Hedge *)graph );
}


/*************************/
/**  Private functions  **/
/*************************/

/*
 * GetMatrix function: build a matrix out from the vector in a hedge describing
 *                     a region in the decomposition at scale s.
 */
static Matrix
GetMatrix ( Hedge *graph, const int s, const int rows, const int cols )
{
    int x, y;
    Matrix out;
    Vector hedgePtr;

    /* Allocate memory */
    out = MEM_Init2D( (long)rows, (long)cols );

    /* Copy data into matrix */
    for ( y=0, hedgePtr=graph->contents[s] ; y<rows ; y++ )
        for ( x=0 ; x<cols ; x++ )
            out[y][x] = *(hedgePtr++);

    return out;
}


/*
 * ReconstructPackets fucntion: given a hedge with a packet decomposition, the
 *                              inverse operation is applied to find the 
 *                              original signal. 
 */
static void
ReconstructPackets ( Hedge *graph, const int x, const int y,
                                   const int N, const int nTilde,
                                   const int l )
{
    int i, j,
        maxLevel,
        x1, y1,
        xy,
        length;
    Matrix K[4];
    Matrix out;
    Vector outPtr, outPtr1;

    maxLevel = l;
    xy = x*y;
    while ( maxLevel ) {

        /* Find the nodes in this level and merge them */
        i = 0;
        while ( i < graph->blocks ) {
            int block;

            /* Find a node in this level */
            while ( graph->levels[i] < maxLevel && i < graph->blocks ) { i++; }
            if ( i >= graph->blocks )
                break;
            block = i;
            length = xy/(1<<(2*graph->levels[block]));
            x1 = y1 = (int)sqrt( (Flt)length );   /* assuming all regions are square */
            /* Take all the four siblings and merge them */
            K[0] = GetMatrix( graph, block  , x1, y1 );
            K[1] = GetMatrix( graph, block+1, x1, y1 );
            K[2] = GetMatrix( graph, block+2, x1, y1 );
            K[3] = GetMatrix( graph, block+3, x1, y1 );
            out = ISFLWT2D ( K, x1<<1, y1<<1, N, nTilde );
            MEM_Free2D( K[0] );
            MEM_Free2D( K[1] );
            MEM_Free2D( K[2] );
            MEM_Free2D( K[3] );
            /* Update pointers and levels */
            graph->blocks -= 3;
            graph->levels[block] -= 1;
            i = block + 1;
            for ( j=i ; j<graph->blocks ; j++ ) {
                graph->contents[j] = graph->contents[j+3];
                graph->levels[j] = graph->levels[j+3];
            }
            /* Copy new node in graph->contents */
            length *= 4;
            for ( j=0, outPtr=graph->contents[block], outPtr1=&out[0][0] ; j<length ; j++ )
                *(outPtr++) = *(outPtr1++);
            MEM_Free2D( out );
        }
        maxLevel--;
        
    }
}

/*
 * GetVector function: build a vector out from the vector in a hedge
 *                     describing a region in the decomposition at scale s.
 */
static Vector
GetVector ( Hedge *graph, const int s, const int cols )
{
    int x;
    Vector out;
    Vector hedgePtr;

    /* Allocate memory */
    out = vector( 0, (long)(cols-1) );

    /* Copy data into vector */
    hedgePtr=graph->contents[s];
    for ( x=0 ; x<cols ; x++ )
        out[x] = *(hedgePtr++);

    return out;
}

/*
 * ReconstructPackets1D function: given a hedge with a packet
 *                                decomposition, the inverse 
 *                                operation is applied to find the
 *                                original signal.
 */
static void
ReconstructPackets1D ( Hedge *graph, const int x,
                                     const int N, const int nTilde,
                                     const int l )
{
    int i, j,
        maxLevel,
        length;
    Vector K[2];
    Vector out;
    Vector outPtr, outPtr1;

    maxLevel = l;
    while ( maxLevel ) {

        /* Find the nodes in this level and merge them */
        i = 0;
        while ( i < graph->blocks ) {
            int block;

            /* Find a node in this level */
            while ( graph->levels[i] < maxLevel && i < graph->blocks ) i++;
            if ( i >= graph->blocks )
                break;
            block = i;
            length = x/(1<<(graph->levels[block]));
            /* Take the two siblings and merge them */
            K[0] = GetVector( graph, block  , length );
            K[1] = GetVector( graph, block+1, length );
            out = ISFLWT1D ( K, length<<1, N, nTilde );
            free( K[0] );
            free( K[1] );
            /* Update pointers and levels */
            graph->blocks--;
            graph->levels[block]--;
            i = block + 1;
            for ( j=i ; j<graph->blocks ; j++ ) {
                graph->contents[j] = graph->contents[j+1];
                graph->levels[j] = graph->levels[j+1];
            }
            /* Copy new node in graph->contents */
            length *= 2;
            for ( j=0, outPtr=graph->contents[block], outPtr1=out ; 
                  j<length ; j++ )
                *(outPtr++) = *(outPtr1++);
            free( out );
        }
        maxLevel--;
    }
}