flwt.c

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

        lSum += *(lambdaPtr);
    }
}

/* CODE NEEDED FOR THE FLWT PACKETS */

/*
 * DSFLWT2D function: This function does a disjoint Separable two-
 *                    dimensional Fast Lifted Wavelet Transform on the
 *                    given vector in[] up to one level.
 *                    It reorganizes the coefficients in Mallat format
 *                    and returns the vector out[] with four matrices
 *                    with 4 rows and width*height/4 columns.
 *                    Each row represents a region of the decomposition
 *                    in the following order: upper left for row 0,
 *                    upper right for row 1, bottom left for row 2, and
 *                    bottom right for row 3.
 *
 *                                  Mallat               +---------+ +---------+
 *  +----------------+       +---------+---------+       |         | |         |
 *  |                |       |         |         |       |    0    | |    1    |
 *  |                |       |    0    |    1    |       |         | |         |
 *  |    in-place    |       |         |         |       +---------+ +---------+
 *  |    forward     |  ==>  +---------+---------+  ==>         Separated
 *  |     flwt       |       |         |         |       +---------+ +---------+
 *  |                |       |    2    |    3    |       |         | |         |
 *  |                |       |         |         |       |    2    | |    3    |
 *  +----------------+       +---------+---------+       |         | |         |
 *                                                       +---------+ +---------+
 */
Matrix *
DSFLWT2D ( Matrix in, const int width, const int height,
                      const int N, const int nTilde )
{
    int x, y,                   /* counters */
        x_2 = ( width == 1 ) ? 0 : width>>1,
        y_2 = ( height == 1 ) ? 0 : height>>1;
    Matrix *out;
    Vector outPtr1, outPtr2;

    /* Allocate memory for the ouput regions matrices (4, X/2, Y/2) */
    out = MEM_Init3D( 4, (long)x_2, (long)y_2 );
    if ( !out )
        Error( "DSFLWT2D", MEMORY_MATRIX, ABORT);

    /* Apply forward transform */
    if ( ODD(N) ) {
	FLWT2D_Haar( in, width, height, 1, FALSE );
    } else {
        FLWT2D( in, width, height, N, nTilde, 1, FALSE );
    }

    /* Reorganize coefficients Mallat format */
    for ( y=0 ; y<height ; y++ )
        InplaceChange1D( &in[y][0], width, 1, FALSE );
    for ( x=0 ; x<width ; x++ )
        InplaceChange1D( &in[0][x], height, width, FALSE );
	
    /* Separate regions */
    for ( y=0, outPtr1=&out[0][0][0], outPtr2=&out[1][0][0] ; y<y_2 ; y++ ) {
        for ( x=0 ; x<x_2 ; x++ )
            *(outPtr1++) = in[y][x];   /* fill region 0 */
        for ( x=x_2 ; x<width ; x++ )
            *(outPtr2++) = in[y][x];   /* fill region 1 */
    }
    for ( y=y_2, outPtr1=&out[2][0][0], outPtr2=&out[3][0][0] ; y<height ; y++ ) {
        for ( x=0 ; x<x_2 ; x++ )
            *(outPtr1++) = in[y][x];   /* fill region 2 */
        for ( x=x_2 ; x<width ; x++ )
            *(outPtr2++) = in[y][x];   /* fill region 3 */
    }

    return out;
}


/*
 * ISFLWT2D function: This function does an adjoint Separable two-
 *                    dimensional Fast Lifted Wavelet Transform on the
 *                    given matrix in[][] up to one level.
 *                    It takes the vectors in in[] to create an initial
 *                    matrix out[][] of with Mallat format and size
 *                    width*height. Then, it reorganizes it in Lifiting
 *                    format and applies the inverse transform. The result
 *                    is the in-place version of the inverse transform.
 *
 *                                  Mallat               +---------+ +---------+
 *  +----------------+       +---------+---------+       |         | |         |
 *  |                |       |         |         |       |    0    | |    1    |
 *  |                |       |    0    |    1    |       |         | |         |
 *  |    in-place    |       |         |         |       +---------+ +---------+
 *  |    inverse     |  <==  +---------+---------+  <==         Separated
 *  |      flwt      |       |         |         |       +---------+ +---------+
 *  |                |       |    2    |    3    |       |         | |         |
 *  |                |       |         |         |       |    2    | |    3    |
 *  +----------------+       +---------+---------+       |         | |         |
 *                                                       +---------+ +---------+
 */
Matrix
ISFLWT2D ( Matrix *in, const int width, const int height, 
                       const int N, const int nTilde )
{
    int x, y,                /* counters */
        x_2 = width>>1,
        y_2 = height>>1;
    Matrix out;
    Vector inPtr1, inPtr2;

    /* Allocate memory for output coefficients matrix (width, height) */
    out = MEM_Init2D( (long)width, (long)height );
    if ( !out )
        Error( "SFLWT2D", MEMORY_MATRIX, ABORT);

    /* Adjoint regions */
    for ( y=0, inPtr1=&in[0][0][0], inPtr2=&in[1][0][0] ; y<y_2 ; y++ ) {
        for ( x=0 ; x<x_2 ; x++ )
            out[y][x] = *(inPtr1++);   /* fill region 0 */
        for ( x=x_2 ; x<width ; x++ )
            out[y][x] = *(inPtr2++);   /* fill region 1 */
    }
    for ( y=y_2, inPtr1=&in[2][0][0], inPtr2=&in[3][0][0] ; y<height ; y++ ) {
        for ( x=0 ; x<x_2 ; x++ )
            out[y][x] = *(inPtr1++);   /* fill region 2 */
        for ( x=x_2 ; x<width ; x++ )
            out[y][x] = *(inPtr2++);   /* fill region 3 */
    }

    /* Reorganize coefficients in Lifting format */
    for ( x=0 ; x<width ; x++ )
        InplaceChange1D( &out[0][x], height, width, TRUE );
    for ( y=0 ; y<height ; y++ )
        InplaceChange1D( &out[y][0], width, 1, TRUE );

    /* Apply in-place inverse transform */
    if ( ODD(N) ) {
	FLWT2D_Haar( out, width, height, 1, TRUE );
    } else {
        FLWT2D( out, width, height, N, nTilde, 1, TRUE );
    }

    return out;
}

/*
 * IFLWTP2D function: This function does an inverse two-dimensional
 *                    Fast Lifted Wavelet Transform on the given matrix
 *                    in[][] up to one level.
 *                    The matrix in[][] is in Mallat format with size
 *                    width*height. The routine reorganizes it in Lifiting
 *                    format and applies the inverse transform. The result
 *                    is the in-place version of the inverse transform.
 *
 *                                  Mallat        
 *  +----------------+       +---------+---------+
 *  |                |       |         |         |
 *  |                |       |    0    |    1    |
 *  |    in-place    |       |         |         |
 *  |    inverse     |  <==  +---------+---------+
 *  |      flwt      |       |         |         |
 *  |                |       |    2    |    3    |
 *  |                |       |         |         |
 *  +----------------+       +---------+---------+
 */
void
IFLWTP2D ( Matrix in, const int width, const int height,
                      const int N, const int nTilde )
{
    int x, y;   /* counters */

    /* Reorganize coefficients in Lifting format */
    for ( x=0 ; x<width ; x++ )
        InplaceChange1D( &in[0][x], height, width, TRUE );
    for ( y=0 ; y<height ; y++ )
        InplaceChange1D( &in[y][0], width, 1, TRUE );

    /* Apply in-place inverse transform */
    if ( ODD(N) ) {
	FLWT2D_Haar( in, width, height, 1, TRUE );
    } else {
        FLWT2D( in, width, height, N, nTilde, 1, TRUE );
    }
}

/*
 * DSFLWT1D function: This function does a disjoint Separable one-
 *                    dimensional Fast Lifted Wavelet Transform on the
 *                    given vector in up to one level.
 *                    It reorganizes the coefficients in Mallat format
 *                    and returns the vector out[] with two vectors
 *                    with width/2 columns.
 *
 *       In-Place                 Mallat               
 *  +----------------+     +---------+---------+     +---------+ +---------+
 *  |    (forward)   | ==> |         |         | ==> |         | | |
 *  +----------------+     +---------+---------+     +---------+ +---------+
 */
Vector *
DSFLWT1D ( Vector in, const int width, const int N, const int nTilde )
{
    int x,                   /* counter */
        x_2 = ( width == 1 ) ? 0 : width>>1;
    Vector *out;

    /* Allocate memory for the ouput regions vectors (2, X/2) */
    out = MEM_InitVecs( 2, (long)x_2 );
    if ( !out )
        Error( "DSFLWT1D", MEMORY_MATRIX, ABORT);

    /* Apply forward transform */
    FLWT1D( in, width, N, nTilde, 1, FALSE );

    /* Reorganize coefficients Mallat format */
    InplaceChange1D( in, width, 1, FALSE );
	
    /* Separate regions */
    for ( x=0 ; x<x_2 ; x++ ) {
        out[0][x] = in[x];   /* fill region 0 */
        out[1][x] = in[x+x_2];   /* fill region 1 */
    }

    return out;
}


/*
 * ISFLWT1D function: This function does an adjoint Separable one-
 *                    dimensional Fast Lifted Wavelet Transform on the
 *                    given vectors in[] up to one level.
 *                    It takes the vectors in in to create an initial
 *                    vector out[] of with Mallat format and size
 *                    width. Then, it reorganizes it in Lifting
 *                    format and applies the inverse transform. The
 *                    result is the in-place version of the inverse 
 *                    transform.
 *
 *       In-Place                 Mallat               
 *  +----------------+     +---------+---------+     +---------+ +---------+
 *  |    (inverse)   | <== |         |         | <== |         | | |
 *  +----------------+     +---------+---------+     +---------+ +---------+
 */
Vector
ISFLWT1D ( Vector *in, const int width, const int N, const int nTilde )
{
    int x,                /* counter */
        x_2 = width>>1;
    Vector out;

    /* Allocate memory for output coefficients vector (width) */
    out = vector( 0, (long)(width-1) );
    if ( !out )
        Error( "ISFLWT1D", MEMORY_MATRIX, ABORT);

    /* Adjoint regions */
    for ( x=0 ; x<x_2 ; x++ ) {
        out[x] = in[0][x];   /* fill region 0 */
        out[x+x_2] = in[1][x];   /* fill region 1 */
    }

    /* Reorganize coefficients in Lifting format */
    InplaceChange1D( out, width, 1, TRUE );

    /* Apply in-place inverse transform */
    FLWT1D( out, width, N, nTilde, 1, TRUE );

    return out;
}

/*
 * IFLWTP1D function: This function does an inverse one-dimensional
 *                    Fast Lifted Wavelet Transform on the given vector
 *                    in[] up to one level.
 *                    The vector in[] is in Mallat format with size
 *                    width. The routine reorganizes it in Lifting
 *                    format and applies the inverse transform. The
 *                    result is the in-place version of the inverse 
 *                    transform.
 *
 *       In-Place                   Mallat        
 *  +----------------+       +---------+---------+
 *  |    (inverse)   |  <==  |         |         |
 *  +----------------+       +---------+---------+
 */
void
IFLWTP1D ( Vector in, const int width, const int N, const int nTilde )
{
    /* Reorganize coefficients in Lifting format */
    InplaceChange1D( in, width, 1, TRUE );

    /* Apply in-place inverse transform */
    FLWT1D( in, width, N, nTilde, 1, TRUE );
}