flwt.c

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

            if (nY-- > 0) {
                /* Apply forward transform to the columns */
                for ( x=0 ; x<width ; x+=step ) {
                    FLWT1D_Predict ( &Data[0][x], (long)height, (long)width,
                                     (long)step, (long)N );
                    FLWT1D_Update ( &Data[0][x], (long)height, (long)width,
                                    (long)step, (long)nTilde, liftingY[nY] );
                }
            }
        }
    } else {   /* Inverse transform */
        /* Switch sign of filter & lifting coefficients */
        for ( i=0 ; i<N*(1+(N>>1)) ; i++ )
            filter[i] = -filter[i];                 /* negate filter */
        for ( x=0 ; x<nX ; x++ )
            for ( k=0 ; k<(width>>1)*nTilde ; k++ )
                liftingX[x][k] = -liftingX[x][k];   /* negate lifting in X */
        for ( y=0 ; y<nY ; y++ )
            for ( k=0 ; k<(height>>1)*nTilde ; k++ )
                liftingY[y][k] = -liftingY[y][k];   /* negate lifting in Y */

        /* 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 ) {
                    FLWT1D_Update ( &Data[0][x], (long)height, (long)width,
                                    (long)step, (long)nTilde, liftingY[nY-n] );
                    FLWT1D_Predict ( &Data[0][x], (long)height, (long)width,
                                     (long)step, (long)N );
                }
            }
    
            if (n <= nX) {
                /* Apply inverse transform to the rows */
                for ( y=0 ; y<height ; y+=step ) {
                    FLWT1D_Update ( &Data[y][0], (long)width, 1,
                                    (long)step, (long)nTilde, liftingX[nX-n] );
                    FLWT1D_Predict ( &Data[y][0], (long)width, 1,
                                     (long)step, (long)N );
                }
            }
    
            n--;   /* one less level to go */
        }

        /* Switch sign of filter & lifting coefficients again */
        for ( i=0 ; i<N*(1+(N>>1)) ; i++ )
            filter[i] = -filter[i];                 /* negate filter */
        for ( x=0 ; x<nX ; x++ )
            for ( k=0 ; k<(width>>1)*nTilde ; k++ )
                liftingX[x][k] = -liftingX[x][k];   /* negate lifting in X */
        for ( y=0 ; y<nY ; y++ )
            for ( k=0 ; k<(height>>1)*nTilde ; k++ )
                liftingY[y][k] = -liftingY[y][k];   /* negate lifting in Y */
    }
}

/*
 * FLWT1D function: This is a 1D Direct/Inverse Fast Lifted Wavelet
 *                  Trasform. Since the Lifting Scheme is used, the final
 *                  coefficients are organized in-place in the original 
 *                  vector Data[]. In file mallat.c there are
 *                  functions provided to change this organization to
 *                  the Isotropic format originally established by Mallat.
 */
void
FLWT1D ( Vector Data, const int width,
                      const int N, const int nTilde,
                      const int levels, const boolean inverse )
{
    int i, k, x, s,       /* counters */
        step, max_step,   /* step size and max./init. step */
        nX;               /* number of iterations in X and Y directions */
    Flt maxN;             /* maximum number of vanishing moments */


    /* Verify for existence of filter and lifting coefficients */
    if (!filter) {
        Error ("FLWT1D", NO_FILTER_VECTOR, ABORT);
    }

    if ( !liftingX && !liftingY ) {
	fprintf( stderr,
	         "\nFLWT1D(): lifting coefficients are not initialized."
		 "\nUse GetDual() and GetReal() before calling this "
                 "function.\n");
        exit( INPUT_ERROR );
    }

    /* 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))    */
    maxN = (Flt)MAX(N, nTilde) - (Flt)1;   /* max vanishing moments */
    /* Iterations in X */
    nX = (width == 1) ? 0 : (int)logBaseN( (Flt)(width-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN(levels, nX);   /* find lowest level */


    if ( !inverse ) {    /* Forward transform */
        /* Calculate wavelet coefficients */
        max_step = (nX==0) ? 0 : 1<<nX;   /* maximum step */
        for ( step=1, s=nX-1 ; step<max_step ; step<<=1, s-- ) {
            /* Apply forward transform to the vector */
            FLWT1D_Predict ( Data, (long)width, 1,
                             (long)step, (long)N );
            FLWT1D_Update ( Data, (long)width, 1,
                            (long)step, (long)nTilde, liftingX[s] );
        }
    } else {   /* Inverse transform */
        /* Switch sign of filter & lifting coefficients */
        for ( i=0 ; i<N*(1+(N>>1)) ; i++ )
            filter[i] = -filter[i];                 /* negate filter */
        for ( x=0 ; x<nX ; x++ )
            for ( k=0 ; k<(width>>1)*nTilde ; k++ )
                liftingX[x][k] = -liftingX[x][k];   /* negate lifting in X */

        /* Find original dataset */
        max_step = (nX==0) ? 0 : 1<<(nX-1);   /* initial step */
        for ( step=max_step, s=0 ; step>=1 ; step>>=1, s++ ) {
            /* Apply inverse transform to the vector */
            FLWT1D_Update ( Data, (long)width, 1,
                            (long)step, (long)nTilde, liftingX[s] );
            FLWT1D_Predict ( Data, (long)width, 1,
                             (long)step, (long)N );
        }

        /* Switch sign of filter & lifting coefficients again */
        for ( i=0 ; i<N*(1+(N>>1)) ; i++ )
            filter[i] = -filter[i];                 /* negate filter */
        for ( x=0 ; x<nX ; x++ )
            for ( k=0 ; k<(width>>1)*nTilde ; k++ )
                liftingX[x][k] = -liftingX[x][k];   /* negate lifting in X */
    }
}

/*
 * FLWT1D_Predict function: The Gamma coefficients are found as an average
 *                          interpolation of their neighbors in order to find
 *                          the "failure to be linear" or "failure to be
 *                          cubic" or the failure to be the order given by
 *                          N-1. This process uses the filter coefficients
 *                          stored in the filter vector and predicts
 *                          the odd samples based on the even ones
 *                          storing the difference in the gammas.
 *                          By doing so, a Dual Wavelet is created.
 */
void
FLWT1D_Predict ( Vector vect, const long size, const long incr,
                              const long step, const long N )
{
    register Vector lambdaPtr,       /* pointer to Lambda coefficients */
                    gammaPtr,        /* pointer to Gamma coefficients */
                    fP, filterPtr;   /* pointers to filter coefficients */
    register long len,               /* number of coefficients at current level */
                  j,                 /* counter for the filter cases */
                  stepIncr;          /* step size between coefficients of the same type */

    long stop1,                      /* number of cases when L < R */
         stop2,                      /* number of cases when L = R */
         stop3,                      /* number of cases when L > R */
         soi;                        /* increment for the middle cases */

    /************************************************/
    /* 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 for each case */
    /************************************************/
    j     = ODD(len);
    stop1 = N >> 1;
    stop3 = stop1 - j;                /* L > R */
    stop2 = (len >> 1) - N + 1 + j;   /* L = R */
    stop1--;                          /* L < R */

    /***************************************************/
    /* Predict Gamma (wavelet) coefficients (odd guys) */
    /***************************************************/

    filterPtr = filter + N;   /* position filter pointer */

    /* Cases where nbr. left Lambdas < nbr. right Lambdas */
    gammaPtr = vect + (stepIncr >> 1);   /* second coefficient is Gamma */
    while(stop1--) {
        lambdaPtr = vect;   /* first coefficient is always first Lambda */
        j = N;
        do {   /* Gamma update (Gamma - predicted value) */
            *(gammaPtr) -= (*(lambdaPtr)*(*(filterPtr++)));
            lambdaPtr   += stepIncr;   /* jump to next Lambda coefficient */
        } while(--j);   /* use all N filter coefficients */
        /* Go to next Gamma coefficient */
        gammaPtr += stepIncr;
    }

    /* Cases where nbr. left Lambdas = nbr. right Lambdas */
    soi = 0;
    while(stop2--) {
        lambdaPtr = vect + soi;   /* first Lambda to be read */
        fP = filterPtr;   /* filter stays in same values for this cases */
        j = N;
        do {   /* Gamma update (Gamma - predicted value) */
            *(gammaPtr) -= (*(lambdaPtr)*(*(fP++)));
            lambdaPtr   += stepIncr;   /* jump to next Lambda coefficient */
        } while(--j);   /* use all N filter coefficients */
        /* Move start point for the Lambdas */
        soi += stepIncr;
        /* Go to next Gamma coefficient */
        gammaPtr += stepIncr;
    }

    /* Cases where nbr. left Lambdas > nbr. right Lambdas */
    fP = filterPtr;   /* start going backwards with the filter coefficients */
    vect += (soi-stepIncr);   /* first Lambda is always in this place */
    while (stop3--) {
        lambdaPtr = vect;   /* position Lambda pointer */
        j = N;
        do {   /* Gamma update (Gamma - predicted value) */
            *(gammaPtr) -= (*(lambdaPtr)*(*(--fP)));
            lambdaPtr   += stepIncr;   /* jump to next Lambda coefficient */
        } while(--j);   /* use all N filter coefficients */
        /* Go to next Gamma coefficient */
        gammaPtr += stepIncr;
    }
}

/*
 * FLWT1D_Update: the Lambda coefficients have to be "lifted" in order
 *                to find the real wavelet coefficients. The new Lambdas
 *                are obtained  by applying the lifting coeffients stored
 *                in the lifting vector together with the gammas found in
 *                the prediction stage. This process assures that the
 *                moments of the wavelet function are preserved.
 */
void
FLWT1D_Update ( Vector vect, const long size, const long incr,
                             const long step, const long nTilde,
                             const Vector lc )
{
    register Vector lcPtr,   /* pointer to lifting coefficient values */
                    vL,      /* pointer to Lambda values */
                    vG;      /* pointer to Gamma values */
    register long j;         /* counter for the lifting cases */

    long len,                /* number of coefficietns at current level */
         stop1,              /* number of cases when L < R */
         stop2,              /* number of cases when L = R */
         stop3,              /* number of cases when L > R */
         noGammas,           /* number of Gamma coefficients at this level */
         stepIncr,           /* step size between coefficients of the same type */
         soi;                /* increment for the middle cases */


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

    /************************************************/
    /* Calculate number of iterations for each case */
    /************************************************/
    j	  = ODD(len);
    stop1 = nTilde >> 1;
    stop3 = stop1 - j;                   /* L > R */
    stop2 = noGammas - nTilde + 1 + j;   /* L = R */
    stop1--;                             /* L < R */

    /**********************************/
    /* Lift Lambda values (even guys) */
    /**********************************/

    lcPtr = lc;   /* position lifting pointer */

    /* Cases where nbr. left Lambdas < nbr. right Lambdas */
    vG = vect + (stepIncr >> 1);   /* second coefficient is Gamma */
    while(stop1--) {
        vL = vect;   /* lambda starts always in first coefficient */
        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;