flwtfilt.c

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

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: flwtfilt.c,v 1.5 1996/09/17 16:10:09 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/flwtfilt.c,v $
 *  Author: Gabriel Fernandez, Senthil Periaswamy
 *
 *  A filter operator is included to process images. The idea is
 *  to modify the predicted values of the Gammas by some factor in
 *  order to smooth or enhance the original dataset.
 *
 *  The factor is called beta. If it is less than 1, the original
 *  image is smoothed. If it is greater than 1, the image is enhanced.
 */
/* do not edit anything above this line */

/* System header files */
#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* FLWT header files */
#include <flwterr.h>
#include <flwtfilt.h>
#include <util.h>

/* Static declarations */
/* static void (*pFilterFunction)(); */
static void
FLWTFilter1D ( Vector, const int, const int, const int, const Flt );
static void
FLWTFilter1DC ( Vector, const int, const int, const int,
                        const Flt, const Flt,
			const int, const int );

/* code */

/*
 * FilterQuery function: asks the user about the kind of filter that
 *                       has to be applied to the input dataset.
 */
Beta
FilterQuery ( void )
{
    Beta beta;
    int opt;
    char input[BUFSIZ];
    boolean exit;
    
    fprintf (stdout, "\n>>> Filter Operator Query <<<\n");
    fprintf (stdout, "The following questions are useful to determine the\n"
                     "parameters needed for the filter operator.\n");
    /* Type of filter operator */
    do {
        fprintf (stdout, "\nType of filter: [L]inear, [C]ircular, "
	                 "[R]ectangular, [F]unction ? ");
        opt = tolower ( getc (stdin) );
	fflush (stdin);
	exit = TRUE;
	if ( opt == 'l' ) {
	    beta.mode = LINEAR;
	    /* pFilterFunction = (&FLWTFilter1D); */
	} else if ( opt == 'c' ) {
	    beta.mode = CIRCULAR;
	    /* pFilterFunction = NULL; */
        } else if ( opt == 'r' ) {
	    beta.mode = RECTANGULAR;
	    /* pFilterFunction = NULL; */
        } else if ( opt == 'f' ) {
	    beta.mode = FUNCTION;
	    /* pFilterFunction = NULL; */
	} else {
            Error ("FilterQuery", INVALID_ANSWER, RETURN);
            exit = FALSE;
        }
    } while (!exit);
    if ( opt == 'l' ) {
	do {
	    fprintf (stdout, "\nDifferent beta factors for X and Y [Y, N]? ");
	    opt = tolower ( getc(stdin) );
	    fflush (stdin);
            exit = TRUE;
            if ( opt == 'n' ) {
	        do {
	            fprintf (stdout, "\nGive value of beta: ");
	            gets (input);
		    fflush (stdin);
	            beta.val1 = beta.val2 = (Flt)atof(input);
	            if ( beta.val1 <= (Flt)0 ) {
		        Error ("FilterQuery", INVALID_RELATIONAL_FLT, RETURN,
                               "Beta", ">", 0.0);
		        exit = FALSE;
	            } else {
		        exit = TRUE;
		    }
	        } while (!exit);
            } else if ( opt == 'y' ) {
	        do {
	            fprintf (stdout, "\nGive value of beta for X: ");
	            gets (input);
		    fflush (stdin);
	            beta.val1 = (Flt)atof(input);
	            if ( beta.val1 <= (Flt)0 ) {
		        Error ("FilterQuery", INVALID_RELATIONAL_FLT, RETURN,
                               "Beta", ">", 0.0);
		        exit = FALSE;
	            } else {
		        exit = TRUE;
		    }
	        } while (!exit);
		do {
	            fprintf (stdout, "\nGive value of beta for Y: ");
	            gets (input);
		    fflush (stdin);
	            beta.val2 = (Flt)atof(input);
	            if ( beta.val2 <= (Flt)0 ) {
		        Error ("FilterQuery", INVALID_RELATIONAL_FLT, RETURN,
                               "Beta", ">", 0.0);
		        exit = FALSE;
	            } else {
		        exit = TRUE;
		    }
	        } while (!exit);
            } else {
	        Error ("FilterQuery", INVALID_ANSWER, RETURN);
                exit = FALSE;
	    }
	} while (!exit);
    } else if ( opt == 'c' || opt == 'r' ) {
	do {
	    fprintf (stdout, "\nGive value of beta in the middle: ");
	    gets (input);
	    fflush (stdin);
	    beta.val1 = (Flt)atof(input);
	    if ( beta.val1 <= (Flt)0 ) {
                Error ("FilterQuery", INVALID_RELATIONAL_FLT, RETURN,
                       "Beta", ">", 0.0);
		exit = FALSE;
	    } else {
		exit = TRUE;
	    }
	} while (!exit);
	do {
	    fprintf (stdout, "\nGive value of beta in the border: ");
	    gets (input);
	    fflush (stdin);
	    beta.val2 = (Flt)atof(input);
	    if ( beta.val2 <= (Flt)0 ) {
                Error ("FilterQuery", INVALID_RELATIONAL_FLT, RETURN,
                       "Beta", ">", 0.0);
		exit = FALSE;
	    } else {
		exit = TRUE;
            }
	} while (!exit);
    } else if ( opt == 'f' ) {
        Error ("FilterQuery", CUSTOM, ABORT,
                "Function filter is not supported yet. Exiting...\n");
    } else {
        Error ("FilterQuery", FATAL_ERROR, ABORT);
        /* NOTREACHED */
    }

    /* Return the beta structure filled up */
    return beta;
}

/*
 * FilterAPI function: manages the beta structure to see what kind of
 *                     operations have to be done.
 */
void
FilterAPI ( Matrix Data, const int sizeX, const int sizeY,
                         const int N, const int nTilde,
                         const Beta beta, const int levels )
{
    if ( beta.mode == LINEAR ) {
	FLWTFilter2D ( Data, sizeX, sizeY, N, nTilde, beta, levels );
    } else if ( beta.mode == CIRCULAR ) {
	FLWTFilter2D ( Data, sizeX, sizeY, N, nTilde, beta, levels );
    } else if ( beta.mode == RECTANGULAR ) {
	fprintf (stderr, "Not implemented rectangular option yet."
	                 "Nothing is done.\n");
    } else {
	Error ("FilterAPI", FATAL_ERROR, RETURN);
        /* NOTREACHED */
    }
}

/*
 * FLWTFilter2D function: returns an enhanced image by multiplying
 *                        each wavelet detail coefficient at a particular
 *                        scale s by some weight w_s. The weights are
 *                        different for each scale. Data is the image you
 *                        wish to enhance, and Beta is a base weight factor
 *                        you wish to apply to the detail coefficients in
 *                        each scale. The original beta value is applied
 *                        the coarsest level, (beta)^2 is applied to the
 *                        previous level, and so on, until (beta)^max_level
 *                        is applied to the finest leveli (level 1).
 */
void
FLWTFilter2D ( Matrix Data, const int sizeX, const int sizeY,
                            const int N, const int nTilde,
                            const Beta beta, const int levels )
{
    int x, y,
        n,          /* maximum number levels */
        nX, nY,     /* number of iterations in X and Y directions */
        step,       /* step size of current level */
        max_step;   /* maximum step for coefficients in last level */
    Flt maxN,       /* maximum number of vanishing moments */
        fBetaX,     /* beta factor for the X direction */
        fBetaY;     /* beta factor for the Y direction */

    /* 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 = (sizeX == 1) ? 0 : (int)logBaseN ( (Flt)(sizeX-1)/maxN, (Flt)2 );
    nX = (nX < 0) ? 0 : MIN (levels, nX);   /* find lowest level */
    /* Iterations in Y */
    nY = (sizeY == 1) ? 0 : (int)logBaseN ( (Flt)(sizeY-1)/maxN, (Flt)2 );
    nY = (nY < 0) ? 0 : MIN (levels, nY);   /* find lowest level */
    /* Total number of iterations */
    n  = MAX(nX, nY);   /* find minimum number of iterations */

    /* Calculate the maximum beta factor for each direction */
    fBetaX = (Flt)1;
    for ( x=0 ; x<nX ; x++ )
        fBetaX *= beta.val1;
    fBetaY = (Flt)1;
    for ( y=0 ; y<nY ; y++ )
        fBetaY *= beta.val2;

    /* Apply beta operator to rows and columns */
    max_step = (n==0) ? 0 : 1<<n;   /* maximum step */
    for ( step=1 ; step<max_step ; step<<=1 ) {
        /* For the rows */
        if (nX-- > 0) {
            for ( y=0 ; y<sizeY ; y+=step ) {
	        if ( beta.mode == LINEAR ) {
                    FLWTFilter1D ( &Data[y][0], sizeX, 1, step, fBetaX );
		} else if ( beta.mode == CIRCULAR ) {
		    FLWTFilter1DC ( &Data[y][0], sizeX, 1, step,
		                    beta.val1, beta.val2, y, nX );
		}
                /* (*pFilterFunction) ( &Data[y][0], sizeY, 1, step, fBetaX ); */
            }
	    fBetaX /= beta.val1;
        }
        /* For the columns */
        if (nY-- > 0) {
            for ( x=0 ; x<sizeX ; x+=step ) {
	        if ( beta.mode == LINEAR ) {
                    FLWTFilter1D ( &Data[0][x], sizeY, sizeX, step, fBetaY );
		} else if ( beta.mode == CIRCULAR ) {
		    FLWTFilter1DC ( &Data[0][x], sizeY, sizeX, step,
		                    beta.val1, beta.val2, x, nY );
		}
                /* (*pFilterFunction) ( &Data[0][x], sizeY, 1, step, fBetaY ); */
            }
	    fBetaY /= beta.val2;
        }
    }
}

/*
 * FLWTFilter1D function: applies the given factor beta to the given set of
 *                        coefficients.
 */
static void
FLWTFilter1D ( Vector Data, const int len, const int incr, const int step,
                            const Flt beta )
{
    register int gammaIdx;
    int stepIncr, size;

    /* Initialize some variable */
    stepIncr = (step*incr) << 1;
    gammaIdx = stepIncr >> 1;
    size = len*incr;

    /* Step through Gammas and apply Beta factor */
    while (gammaIdx < size) {
        Data[gammaIdx] *= beta;
        gammaIdx += stepIncr;
    }
}

/*
 * FLWTFilter1DC function: applies a circular filter to the given set of
 *                         wavelet coefficients.
 */
static void
FLWTFilter1DC ( Vector Data, const int len, const int incr, const int step,
                             const Flt betaStart, const Flt betaStop,
			     const int xy, const int levels )
{
    register Vector gammaPtr;
    int stepIncr, x, step2, k;
    Flt i, j, s, beta, interval = betaStop - betaStart;

    if ( xy > len ) {
        Error ( "FLWTFilter1DC", INVALID_RELATIONAL_INT, ABORT, "xy", "<=", len );
    }

    /* Initialize some variables */
    x = 0;
    step2 = step << 1;
    stepIncr = step2*incr;
    gammaPtr = Data + (stepIncr >> 1);

    /* Step through Gammas and apply Beta factor */
    for ( x=0 ; x<len ; x+=step2 ) {
        Flt dist;

	i = (Flt)x / (Flt)len - (Flt)0.5;
	j = (Flt)xy / (Flt)len - (Flt)0.5;
        dist = j*j + i*i;
	s = dist*(Flt)len/(Flt)255;   /* s lies b/w 0 and 1 */
	if ( betaStart <= betaStop ) {
	    s = s*interval + betaStart; 
        } else {
            s = (Flt)1 - s;
            s = betaStop - s*interval;
        }
	beta = s;
	for ( k=0 ; k<levels ; k++ )
	    beta *= s;
        *(gammaPtr) *= beta;
        gammaPtr += stepIncr;
    }
}