mallat.c

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

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: mallat.c,v 1.5 1996/10/06 09:24:54 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/mallat.c,v $
 *  Author: Gabriel Fernandez, Senthil Periaswamy
 *
 *  This file contains routines to manipulate the results of an in-place
 *  transform in order to re-arrange the coefficients in Mallat format.
 */
/* do not edit anything above this line */

/* System header files */
#include <stdio.h>
#include <math.h>
/* FLWT header files */
#include <flwtdef.h>
#include <mallat.h>
#include <mem.h>
#include <util.h>

/* code */

/*
 * FLWTChangeFormat function: changes a given image from the Lifting format
 *                            to the Mallat format, and viceversa.
 */
void
FLWTChangeFormat ( Matrix Data, const int width, const int height,
                                const int N, const int nTilde,
                                const int levels, const boolean inverse )
{
    int i, x, y,
        lenX, lenX2,
	lenY, lenY2, /* widths and heights at each level of calculations */
        n,           /* maximum number levels */
        nX, nY;      /* number of iterations in X and Y directions */
    Flt maxN;        /* maximum number of vanishing moments */

    /* 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))    */

    if ( N==1 || nTilde==1 ) {   /* Haar wavelets */
        /* Iterations in X */
        nX = (width==1) ? 0 :
             (int)ceil((Flt)logBaseN((Flt)width,2.0));
        /* Iterations is Y */
        nY = (height==1) ? 0 :
             (int)ceil((Flt)logBaseN((Flt)height,2.0));
    } else {   /* Biorthogonal wavelets */
        /* Max vanishing moments */
        maxN = (Flt)MAX(N, nTilde) - (Flt)1;
        /* Iterations in X */
        nX = (width==1) ? 0 :
             (int)logBaseN((Flt)(width-1)/maxN,2.0);
        /* Iterations in Y */
        nY = (height==1) ? 0 :
             (int)logBaseN((Flt)(height-1)/maxN,2.0);
    }
    /* Find lowest levels */
    nX = (nX < 0) ? 0 : MIN (levels, nX);
    nY = (nY < 0) ? 0 : MIN (levels, nY);
    /* Total number of iterations */
    n  = MAX(nX, nY);

    if ( !inverse ) {

        /* Initial lengths */
        lenX = lenX2 = width;
        lenY = lenY2 = height;

        /* Apply Mallat ordering to rows and columns */
        for ( i=0 ; i<n ; i++ ) {
            /* For the rows */
            if (nX-- > 0) {
                for ( y=0 ; y<lenY2 ; y++ ) {
                    InplaceChange1D( &Data[y][0], lenX, 1, inverse );
                }
	        lenX  -= (lenX  >> 1);
	        lenY2 -= (lenY2 >> 1);
            }
            /* For the columns */
            if (nY-- > 0) {
                for ( x=0 ; x<lenX2 ; x++ ) {
                    InplaceChange1D( &Data[0][x], lenY, width, inverse );
                }
	        lenY  -= (lenY  >> 1);
	        lenX2 -= (lenX2 >> 1);
            }
        }

    } else {   /* inverse */

        /* Apply Lifting ordering to rows and columns */
        for ( i=n-1 ; i>=0 ; i-- ) {
            int j;

            /* Initial lengths */
            lenX = lenX2 = width;
            lenY = lenY2 = height;

            /* For the columns */
            if (nY-- > 0) {
                /* Calculate new length */
                for ( j=0 ; j<i ; j++ ) {
                    lenY  -= (lenY  >> 1);
                    lenX2 -= (lenX2 >> 1);
                }
                /* Reorganize coeffs */
                for ( x=0 ; x<lenX2 ; x++ ) {
                    InplaceChange1D( &Data[0][x], lenY, width, inverse );
                }
            }
            /* For the rows */
            if (nX-- > 0) {
                /* Calculate new length */
                for ( j=0 ; j<i ; j++ ) {
                    lenX  -= (lenX  >> 1);
                    lenY2 -= (lenY2 >> 1);
                }
                /* Reorganize coeffs */
                for ( y=0 ; y<lenY2 ; y++ ) {
                    InplaceChange1D( &Data[y][0], lenX, 1, inverse );
                }
            }
        }

    }
 }

/*
 * InplaceChange1D function: arranges the wavelet coefficients in the given
 *                           vector to put them in Mallat or Lifting format,
 *                           depending on the value of the inverse variable.
 */
void InplaceChange1D ( Vector Data, const int len, const int offset,
                       const boolean inverse )
{
    register int i;
    int noLambdas, noGammas;
    Vector lPtr, gPtr, bufPtr;
    Vector buffer;

    /* Allocate memory for buffer vector */
    buffer = vector( 0, (long)(len-1) );

    /* Number of Gamma coefficients */
    noGammas  = len>>1;
    noLambdas = len - noGammas;
    
    /* Check whether to do Lifting format or Mallat format */
    if ( !inverse ) {
        /* Set lambda and gamma pointers */
        lPtr = buffer;            /* point at the beginning of vector */
        gPtr = buffer + (noLambdas);   /* point at half of the vector */

        /* Move all lambdas to the left and all gammas to the right */
        for ( i=0 ; i<len-2 ; i+=2 ) {
	    *(lPtr) = Data[i*offset];
            *(gPtr) = Data[(i+1)*offset];
	    lPtr ++;
	    gPtr ++;
        }
        if ( ODD(len) ) {
            i = len-1;
            *(lPtr) = Data[i*offset];
        } else {
            i = len-2;
	    *(lPtr) = Data[i*offset];
            *(gPtr) = Data[(i+1)*offset];
        }
    } else {
        /* Set lambda and gamma pointers */
        lPtr = Data;                   /* point at the beginning of vector */
        gPtr = Data + (noLambdas)*offset;   /* point at half of the vector */

        /* Alternate lambdas and gammas */
        for ( i=0 ; i<len-2 ; i+=2 ) {
	    buffer[i]   = *(lPtr);
            buffer[i+1] = *(gPtr);
	    lPtr += offset;
	    gPtr += offset;
        }
        if ( ODD(len) ) {
            i = len-1;
	    buffer[i] = *(lPtr);
        } else {
            i = len-2;
	    buffer[i]   = *(lPtr);
            buffer[i+1] = *(gPtr);
        }
    }
    
    /* Write results in original vector */
    bufPtr = buffer;
    for ( i=0 ; i<len*offset ; i+=offset ) {
	Data[i] = *(bufPtr++);
    }
    
    /* Free allocated memory */
    free_vector( buffer, 0, (long)(len-1) );
}