/* *  -*- 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 voidFLWTFilter1D ( Vector, const int, const int, const int, const Flt );static voidFLWTFilter1DC ( 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. */BetaFilterQuery ( 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. */voidFilterAPI ( 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). */voidFLWTFilter2D ( 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 voidFLWTFilter1D ( 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 voidFLWTFilter1DC ( 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;    }}