flwtapi.c

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

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: flwtapi.c,v 1.33 1996/11/05 18:17:58 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/flwtapi.c,v $
 *  Author: Gabriel Fernandez
 *
 *  This file has the functions needed to manage all the
 *  options included in the param variable.
 */
/* do not edit anything above this line */

/* System header files */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/* FLWT header files */
#include <constant.h>    /* constants declarations */
#include <flwtdef.h>     /* types and general definitions */
#include <file.h>        /* file utility functions */
#include <flwt.h>        /* wavelet transform functions */
#include <flwtapi.h>     /* parameters structure */
#include <flwtpack.h>    /* wavelet packet transform functions */
#include <image.h>       /* image structure and definitions */
#include <imgmem.h>      /* image manipulation routines */
#include <mallat.h>      /* Mallat format functions */
#include <util.h>        /* useful functions */
#include <flwterr.h>     /* error handling function */


/* Prototypes for Private Static Functions */
static void PrintMatrix ( const Matrix, const int, const int );
static void Point ( Image *, const int, const int, const int, const int );
static void Scale ( Matrix, const int, const int, const int,
                            const int, const int );
static void MinMax ( const Matrix, const int, const int, Flt *, Flt * );

/* code */

/*
 * FLWTAPI function: the name stands for Fast Lifted Wavelet Transform
 *                   Application Program Interface. This function takes
 *                   the param structure and determines what kind of
 *                   operations have to be done depending on the given
 *                   options.
 */
void
FLWTAPI ( ParamStruct param )
{
    int        i;
    Image      image;          /* structure with image and useful info */
    FileFormat filetype;       /* type of input/output file */
    clock_t    startclock,     /* used by the clock() function to measure */
               stopclock;      /* execution time */
    Flt        diff, total=(Flt)0;

    /*******************************************************/
    /* Read image from infile or create dataset if WAVELET */
    /* option is chosen.                                   */
    /*******************************************************/
    if (param.mode != WAVELET) {

        /* Read file type */

        fprintf (stdout, "Read file type ");
        fflush (stdout);
        startclock = clock ();
            filetype = ReadFileType (param.infile);
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        if ( filetype == FMT_ERROR ) {
            Error ("FLWTAPI", NO_FILE, ABORT, param.infile);
            /* NOTREACHED */
        } else if ( filetype == FMT_UNKNOWN ) {
            Error ("FLWTAPI", INVALID_FILE_FORMAT, ABORT, param.infile);
            /* NOTREACHED */
        } else if ( filetype < FMT_ERROR ) {
            Error ("FLWTAPI", NO_READABLE_FORMAT, ABORT, param.infile);
            /* NOTREACHED */
        }
        fprintf (stdout, "\t\t\t\t[%3.2f secs]\n", diff);
        total += diff;

        /* At this point, filetype is a known, readable format */
        /* Read input file with picture to be processed */

        fprintf (stdout, "Read input file");
        fflush (stdout);
        startclock = clock ();
            i = ReadFile (param.infile, filetype, &image);
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        if ( !i ) {
            Error ("FLWTAPI", NO_FILE, ABORT, param.infile);
            /* NOTREACHED */
        }
        param.sizeX = image.width;
        param.sizeY = image.height;
        if ( image.width == 0 || image.height == 0 ) {   /* shouldn't happen */
            Error ("FLWTAPI", INVALID_RELATIONAL_INT, ABORT, 
                   "Dimensions", ">", 0);
            /* NOTREACHED */
        }
        fprintf (stdout, "\t[%4dx%-4d]\t\t[%3.2f secs]\n",
                 image.width, image.height, diff);
        total += diff;
        fprintf (stdout, "[%s]\n", image.fullInfo);

    } else {   /* create initial set to generate Biorthogonal Wavelet */

        fprintf (stdout, "Create Point Image\t\t\t");
        fflush (stdout);
        filetype = FMT_PBM;
        startclock = clock ();
            Point (&image, param.sizeX, param.sizeY, param.x, param.y);
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;

    }

    /* Now that we've read the file or created the point image,
       let's set the output type */
    SetOutputType( param.frmType, filetype, &image, param.packets );

    /***************************/
    /* Print out original data */
    /***************************/
    
    if (param.print) {
        fprintf (stdout, "Original Data Set\n");
        for ( i=0 ; i<image.colorPlanes ; i++ ) {
            fprintf (stdout, "Band #%d\n", i);
            PrintMatrix (image.band[i], image.width, image.height);
        }
    }
    

    /****************************/
    /* Find filter coefficients */
    /****************************/
    if ( EVEN(param.N) ) {
	fprintf (stdout, "\nFilter Coefficients\t\t\t");
	fflush (stdout);
	startclock = clock ();
	    GetDual (image.width, image.height, param.N, param.nTilde,
		     param.print, param.levels);
	stopclock = clock ();
	diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
	fprintf (stdout, "[%3.2f secs]\n", diff);
	total += diff;
    }

    /*****************************/
    /* Find lifting coefficients */
    /*****************************/
    if ( EVEN(param.N) ) {
        fprintf (stdout, "Lifting Coefficients\t\t\t");
        fflush (stdout);
        startclock = clock ();
            GetReal (image.width, image.height, param.N, param.nTilde,
                     param.print, param.levels);
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;
    }

    /***************************/
    /* Apply forward transform */
    /***************************/
    if ((param.mode == FORWARD || param.mode == FILTER) && !param.packets) {
        fprintf (stdout, "\nForward Wavelet Transform\t\t");
        fflush (stdout);
        startclock = clock ();
            for ( i=0 ; i<image.colorPlanes ; i++ ) {
	        if ( ODD(param.N) ) {
		    FLWT2D_Haar( image.band[i], image.width, image.height,
                                 param.levels, FALSE );
		} else {
                    FLWT2D( image.band[i], image.width, image.height,
                            param.N, param.nTilde, param.levels, FALSE );
		}
            }
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;
    }
    if (param.packets && param.mode != INVERSE) {
/*         int rows, cols; */
	int sum, blocks, j;
	int buff[BUFSIZ];

        fprintf (stdout, "\nForward Wavelet Packets Transform\t");
        fflush (stdout);
/*
        if ( image.width == 1 ) { rows = 1; cols = image.height; }
        else if ( image.height == 1 ) { rows = 1; cols = image.width; }
        else { rows = cols = 0; }
*/
        startclock = clock ();
            for ( i=0 ; i<image.colorPlanes ; i++ ) {
/*
                if ( rows == 1 ) {
                    Matrix m = image.band[i];
                    image.levels = DFLWTPackets1D( &m[0][0], cols,
                                                   param.N, param.nTilde,
                                                   param.levels,
                                                   &image.blocks );
                } else {
*/
                    image.levels = DFLWTPackets( image.band[i],
                                                 image.width, image.height,
                                                 param.N, param.nTilde,
                                                 param.levels,
                                                 &image.blocks );
/*
                }
*/
            }
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;
	
	/* Show the distribution of the levels */
        sum = 1; blocks = 1;
        for ( i=1, j=0 ; i<image.blocks ; i++ ) {
	    if ( image.levels[i] == image.levels[i-1] )
		sum++;
	    else {
		buff[j++] = sum;
		buff[j++] = image.levels[i-1];
		sum = 1;
		blocks++;
	    }
	}
	buff[j++] = sum;
	buff[j++] = image.levels[i-1];
	for ( i=0 ; i<2*blocks ; i+=2 )
	    fprintf (stdout, "%d -> %d\n", buff[i], buff[i+1]);
	fflush(stdout);
    }

    /**********************************/
    /* Print out wavelet coefficients */
    /**********************************/
    if (param.print && param.mode != INVERSE) {
        fprintf (stdout, "Wavelet Coefficients\n");
        for ( i=0 ; i<image.colorPlanes ; i++ ) {
            fprintf (stdout, "Band #%d\n", i);
            PrintMatrix (image.band[i], image.width, image.height);
        }
    }

    /**************************************/
    /* Apply filter operator, if required */
    /**************************************/
    if (param.mode == FILTER && !param.packets) {
        fprintf (stdout, "Apply Filter Operator\t\t\t");
        fflush (stdout);
        startclock = clock ();
            for ( i=0 ; i<image.colorPlanes ; i++ ) {
                FilterAPI (image.band[i], image.width, image.height,
                           param.N, param.nTilde, param.beta, param.levels);
            }
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;
    } else if (param.mode == FILTER && param.packets) {
        Error( "FLWTAPI", CUSTOM, RETURN, "Filter operation not supported for "
	       "packets yet.\n" );
    }

    /***********************************************/
    /* Organize data in Mallat format, if required */
    /***********************************************/
    if ((param.mode == FORWARD || param.mode == FILTER) && param.mallat) {
	fprintf (stdout, "Organize in Mallat format\t\t");
	fflush (stdout);
        startclock = clock ();
            for ( i=0 ; i<image.colorPlanes ; i++ ) {
                FLWTChangeFormat (image.band[i], image.width, image.height,
                                  param.N, param.nTilde, param.levels,
				  FALSE);