/* *  -*- 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. */voidFLWTChangeFormat ( 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) );}