flwt.c

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

    }

    /* Cases where nbr. left Lambdas = nbr. right Lambdas */
    soi = 0;
    while(stop2--) {
        vL = vect + soi;   /* first Lambda to be read */
        j = nTilde;
        do {
            *(vL) += (*(vG)*(*(lcPtr++)));   /* lift Lambda (Lambda + lifting value) */
            vL    += stepIncr;               /* jump to next Lambda coefficient */
        } while(--j);   /* use all nTilde lifting coefficients */
        /* Go to next Gamma coefficient */
        vG  += stepIncr;
        /* Move start point for the Lambdas */
        soi += stepIncr;
    }

    /* Cases where nbr. left Lambdas = nbr. right Lambdas */
    vect += (soi - stepIncr);   /* first Lambda is always in this place */
    while(stop3--) {
        vL = vect;   /* position Lambda pointer */
        j = nTilde;
        do {
            *(vL) += (*(vG)*(*(lcPtr++)));   /* lift Lambda (Lambda + lifting value) */
            vL    += stepIncr;               /* jump to next Lambda coefficient */
        } while(--j);   /* use all nTilde lifting coefficients */
        /* Go to next Gamma coefficient */
        vG += stepIncr;
    }
}

/* CODE FOR THE IN-PLACE FLWT (AVERAGE INTERPOLATION CASE, N ODD) */

/*
 * FLWT2D_Haar function: calculates the Haar wavelet coefficients for
 *                       a two dimensional signal. This case corresponds
 *                       to average interpolation with N=1 and wavelets
 *                       with nTilde = 1.
 */
void
FLWT2D_Haar( Matrix Data, const int width, const int height,
                          const int levels, const boolean inverse )
{
    int i, k, x, y, /* counters */
        n,          /* maximum number of levels */
        nX, nY,     /* number of iterations in X and Y directions */
        step,       /* step size of current level */
        max_step;   /* maximum and initial step for the direct and */
                    /* inverse transforms, respectively */
    Flt maxN;       /* maximum number of vanishing moments */

    /* Number of iterations */
    nX = (width==1) ? 0 :
        (int)ceil((Flt)logBaseN((Flt)(width),2.0));
    nX = (nX<0) ? 0 : MIN(nX,levels);
    nY = (height==1) ? 0 :
        (int)ceil((Flt)logBaseN((Flt)(height),2.0));
    nY = (nX<0) ? 0 : MIN(nY,levels);
    n = MAX(nX,nY);


    if ( !inverse ) {    /* Forward transform */
        /* Calculate wavelet coefficients */
        max_step = (n==0) ? 0 : 1<<n;   /* maximum step */
        for ( step = 1 ; step < max_step ; step <<= 1 ) {
            if (nX-- > 0) {
                /* Apply forward transform to the rows */
                for ( y=0 ; y<height ; y+=step ) {
                    FLWT_Haar ( &Data[y][0], (long)width, 1, (long)step );
                }
            }
    
            if (nY-- > 0) {
                /* Apply forward transform to the columns */
                for ( x=0 ; x<width ; x+=step ) {
                    FLWT_Haar ( &Data[0][x], (long)height, (long)width,
                                  (long)step );
                }
            }
        }
    } else {   /* Inverse transform */
        /* Find original dataset */
        max_step = (n==0) ? 0 : 1<<(n-1);   /* initial step */
        for ( step = max_step ; step >= 1 ; step >>= 1 ) {
            if (n <= nY) {
                /* Apply inverse transform to the columns */
                for ( x=0 ; x<width ; x+=step ) {
                    IFLWT_Haar ( &Data[0][x], (long)height, (long)width,
                                   (long)step );
                }
            }
    
            if (n <= nX) {
                /* Apply inverse transform to the rows */
                for ( y=0 ; y<height ; y+=step ) {
                    IFLWT_Haar ( &Data[y][0], (long)width, 1, (long)step );
                }
            }
    
            n--;   /* one less level to go */
        }
    }
}

/*
 * FLWT1D_Haar function: calculates the Haar wavelet coefficients for
 *                       a one dimensional signal. This case corresponds
 *                       to average interpolation with N=1 and wavelets
 *                       with nTilde = 1.
 */
void
FLWT1D_Haar ( Vector Data, const int width,
                           const int levels, const boolean inverse )
{
    int i, k, x, s, /* counters */
        nX;         /* number of iterations in X and Y directions */
    Flt maxN;       /* maximum number of vanishing moments */


    /* Number of iterations */
    nX = (width==1) ? 0 :
        (int)ceil((Flt)logBaseN((Flt)(width),2.0));
    nX = (nX<0) ? 0 : MIN(nX,levels);


    if ( !inverse ) {    /* Forward transform */
        /* Calculate wavelet coefficients */
        for ( s = 0 ; s < nX ; s++ ) {
            /* Apply forward transform to the vector */
            FLWT_Haar ( Data, (long)width, 1, (long)(1<<s) );
        }
    } else {   /* Inverse transform */
        /* Find original dataset */
        for ( s = nX-1 ; s >=0 ; s-- ) {
            /* Apply inverse transform to the vector */
            IFLWT_Haar ( Data, (long)width, 1, (long)(1<<s) ); 
        }
    }
}

/*
 * FLWT_Haar function: finds the Haar and Lambda coefficients for
 *                     the given 1D signal. Boundary cases are not
 *                     considered yet. Thus, all signals must be
 *                     of even lengths.
 */
void
FLWT_Haar ( Vector vect, const long size, const long incr,
                         const long step )
{
    register Vector lambdaPtr,   /* pointer to Lambda coefficients */
                    gammaPtr;    /* pointer to Gamma coefficients */
    register long nL,            /* number of Lambda coefficients */
                  nG,            /* number of Gamma coefficients */
                  stepIncr;      /* step size between coefficients of the same type */
    long len;                    /* number of coefficients at current level */
    Flt  lAve,                   /* average of Lambdas at current level */
         lSum = 0;               /* cumulative of Lambdas before last one */


    /************************************************/
    /* Calculate values of some important variables */
    /************************************************/
    len      = CEIL(size, step);   /* number of coefficients at current level */
    stepIncr = (step*incr) << 1;   /* step size betweeen coefficients */

    /**********************************/
    /* Calculate number of iterations */
    /**********************************/
    nG = (len >> 1);
    nL = len - nG;

    /****************************/
    /* Calculate Lambda average */
    /****************************/
    if ( ODD(len) ) {
        long n = len, s = stepIncr >> 1;
        Flt  t = 0;
        Flt *l = vect;

        while (n--) {
            t += *l;
            l += s;
        }
        lAve = t/(Flt)len;   /* average we've to maintain */
    }

    /**********************************************************/
    /* Calculate Gamma (Haar wavelet) coefficients (odd guys) */
    /* and coarser Lambda coefficients (even guys)            */
    /**********************************************************/
    gammaPtr = vect + (stepIncr >> 1);   /* second coefficient is Gamma */
    lambdaPtr = vect;                    /* first Lambda to be read */
    while (nG--) {
        /* New Gamma (Haar) coefficient */
        *(gammaPtr) -= *(lambdaPtr);
	/* New Lambda coefficient */
        *(lambdaPtr) += (*(gammaPtr)/(Flt)2);
        lSum += *(lambdaPtr);

        /* Go to next Gamma coefficient */
        gammaPtr += stepIncr;
        /* Go to next Lambda coefficient */
        lambdaPtr += stepIncr;
    }


    /*****************************************************/
    /* If ODD(len), then we have one additional case for */
    /* the Lambda calculations. This is a boundary case  */
    /* and, therefore, has different filter values.      */
    /*****************************************************/
    if ( ODD(len) ) {
        *(lambdaPtr) = lAve*(Flt)nL - lSum;
        lSum += *(lambdaPtr);
    }
}

/*
 * IFLWT_Haar function: applies the inverse transform to find Lambda 
 *                      coefficients for the given 1D signal. Boundary 
 *                      cases are not considered yet. Thus, all signals 
 *                      must be of even lengths.
 */
void
IFLWT_Haar( Vector vect, const long size, const long incr,
                         const long step )
{
    register Vector lambdaPtr,   /* pointer to Lambda coefficients */
                    gammaPtr;    /* pointer to Gamma coefficients */
    register long nL,            /* number of Lambda coefficients */
                  nG,            /* number of Gamma coefficients */
                  stepIncr;      /* step size between coefficients of the same type */
    long len;                    /* number of coefficients at current level */
    Flt  lAve,                   /* average of Lambdas at current level */
         lSum = 0;               /* cumulative of Lambdas before last one */


    /************************************************/
    /* Calculate values of some important variables */
    /************************************************/
    len      = CEIL(size, step);   /* number of coefficients at current level */
    stepIncr = (step*incr) << 1;   /* step size betweeen coefficients */

    /**********************************/
    /* Calculate number of iterations */
    /**********************************/
    nG = (len >> 1);
    nL = len - nG;

    /****************************/
    /* Calculate Lambda average */
    /****************************/
    if ( ODD(len) ) {
        long n = nL, s = stepIncr;
        Flt  t = 0;
        Flt *l = vect;

        while (n--) {
            t += *l;
            l += s;
        }
        lAve = t/(Flt)nL;   /* average we've to maintain */
    }


    /**********************************************************/
    /* Calculate Gamma (Haar wavelet) coefficients (odd guys) */
    /* and coarser Lambda coefficients (even guys)            */
    /**********************************************************/
    gammaPtr = vect + (stepIncr >> 1);   /* second coefficient is Gamma */
    lambdaPtr = vect;                    /* first Lambda to be read */
    while (nG--) {
	/* New Lambda coefficient */
        *(lambdaPtr) -= (*(gammaPtr)/(Flt)2);
        lSum += *(lambdaPtr);
	/* New Gamma coefficient */
	*(gammaPtr) += *(lambdaPtr);
        lSum += *(gammaPtr);

        /* Go to next Gamma coefficient */
        gammaPtr += stepIncr;
        /* Go to next Lambda coefficient */
        lambdaPtr += stepIncr;
    }
    
    /*****************************************************/
    /* If ODD(len), then we have one additional case for */
    /* the Lambda calculations. This is a boundary case  */
    /* and, therefore, has different filter values.      */
    /*****************************************************/
    if ( ODD(len) ) {
        *(lambdaPtr) = lAve*(Flt)len - lSum;