/* *  -*- 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 MatrixGetMatrix ( Hedge *__graph, const int __s, const int __x, const int __y );static voidReconstructPackets ( Hedge *__graph, const int __x, const int __y,                                     const int __N, const int __nTilde,                                     const int __l );static VectorGetVector ( Hedge *__graph, const int __s, const int __x );static voidReconstructPackets1D ( 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. */voidIFLWTPackets ( 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. */voidIFLWTPackets1D ( 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 MatrixGetMatrix ( 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 voidReconstructPackets ( 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 VectorGetVector ( 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 voidReconstructPackets1D ( 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--;    }}