flwtapi.c

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

            }
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        total += diff;
    }

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

    /***************************/
    /* Apply inverse transform */
    /***************************/
    if (param.mode != FORWARD) {
        /* Change to Lifting format, if necessary */
        if (param.mallat) {
	    fprintf (stdout, "Organize in Lifting 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,
				      TRUE);
                }
            stopclock = clock ();
            diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
            fprintf (stdout, "[%3.2f secs]\n", diff);
            total += diff;
        }
        /* Apply inverse now */
        if (!param.packets) {
            fprintf (stdout, "\nInverse 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, TRUE );
		    } else {
                        FLWT2D (image.band[i], image.width, image.height,
                                param.N, param.nTilde, param.levels, TRUE);
		    }
                }
            stopclock = clock ();
            diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
            fprintf (stdout, "[%3.2f secs]\n", diff);
            total += diff;
        } else {   /* packets */
/*             int rows, cols; */
	    int sum, blocks, j;
	    int buff[BUFSIZ];

	    /* Show the distribution of the levels */
	    sum = 1; blocks = 1;
	    fprintf (stdout, "\n");
	    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);

	    /* Calculate inverse packets */
            fprintf (stdout, "Inverse Wavelet Packet 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 ) {
                        IFLWTPackets1D( image.band[i][0], cols,
                                        image.blocks, image.levels,
                                        param.N, param.nTilde );
                    } else {
*/
                        IFLWTPackets( image.band[i],
                                      image.width, image.height,
                                      image.blocks, image.levels,
                                      param.N, param.nTilde );
/*
                    }
*/
                }
            stopclock = clock ();
            diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
            fprintf (stdout, "[%3.2f secs]\n", diff);
            total += diff;
        }
    }

    /***********************************************/
    /* Scale values between 0 and 255, if required */
    /***********************************************/
    /* if (param.mode == WAVELET || param.scale) { */
    if (param.scale) {
        fprintf (stdout, "Scale image ");
        fflush (stdout);
        startclock = clock ();
            for ( i=0 ; i<image.colorPlanes ; i++ ) {
                Scale (image.band[i], image.width, image.height,
                       0, 255, param.gamma);
            }
        stopclock = clock ();
        diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
        fprintf (stdout, "[%3.2f secs]\n", diff);
        fflush (stdout);
        total += diff;
    }

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

    /************************/
    /* Write processed data */
    /************************/
    fprintf (stdout, "\nWrite output file");
    fflush (stdout);
    startclock = clock ();
        i = WriteFile (param.outfile, image);
    stopclock = clock ();
    if ( !i ) {
        Error ("FLWTAPI", WRITE_ERROR, ABORT);
        /* NOTREACHED */
    }
    diff = ((Flt)stopclock-(Flt)startclock)/(Flt)CLOCKS_PER_SEC;
    fprintf (stdout, "\t\t[%3.2f secs]\n", diff);
    total += diff;

    fprintf (stdout, "\t\t\t\t       -------------\n");
    fprintf (stdout, "Total execution time\t\t\t[%3.2f secs]\n", total);
    fflush (stdout);

    /*************************/
    /* Free allocated memory */
    /*************************/
    if ( EVEN(param.N) ) {
        FreeCo (param.sizeX, param.sizeY, param.N, param.nTilde, param.levels);
    }
    IMG_FreeMem (image);
}

/** Private functions **/

/* Print out content of matrix A */
static void
PrintMatrix ( const Matrix A, const int sizeX, const int sizeY )
{
    int i, j;

    for ( i=0 ; i<sizeY ; i++ ) {
        fprintf (stdout, "[");
        for ( j=0 ; j<sizeX-1 ; j++ ) {
            fprintf (stdout, "%3.2f, ", A[i][j]);
        }
        fprintf (stdout, "%3.2f]\n", A[i][j]);
    }
}

/* Generate a data set with a point in the middle */
static void
Point ( Image *image, const int width, const int height,
                      const int x, const int y )
{
    int i, j, rv;
    long filesize;
    Matrix A;

    image->width       = width;
    image->height      = height;
    image->band[RED]   = (Matrix)NULL;
    image->band[GREEN] = (Matrix)NULL;
    image->band[BLUE]  = (Matrix)NULL;
    image->comment     = (char *)NULL;
    image->levels      = (int *)NULL;

    /* Allocate memory for the image */
    rv = (int)IMG_InitMem( image, (long)width, (long)height, 1, 8 );
    if ( !rv ) {
        Error ("Point", MEMORY_IMAGE, ABORT);
        /* NOTREACHED */
    }
    image->colType = 1;   /* Greyscale */
    image->frmType = F_PBMRAW;
    filesize = (long)(width * height);
    sprintf(image->fullInfo, "WAV, raw format.  (%ld bytes)", filesize);
    sprintf(image->shrtInfo, "%dx%d WAV.", image->width, image->height);

    /* Fill with zeros and put a one in given position */
    A = (Matrix)(image->band[0]);
    for ( i=0 ; i<height ; i++ ) {
        for ( j=0 ; j<width ; j++ ) {
            if ( i == y && j == x )
                A[i][j] = (Flt)255.0;
            else
                A[i][j] = (Flt)0.0;
        }
    }
}

/* Scale value of image between Min and Max */
static void
Scale ( Matrix A, const int width, const int height,
                  const int Min  , const int Max  , const int gamma )
{
    int i, j;
    Flt  min, max, minT, maxT;

    MinMax (A, width, height, &min, &max);
    fprintf (stdout, "Min(%4.2f) Max(%4.2f)\t", min, max);
    if (min < (Flt)Min) {
        for ( i=0 ; i<height ; i++ )
            for ( j=0 ; j<width ; j++ )
                A[i][j] += (-min);        /* shift numbers to a positive domain */
        MinMax (A, width, height, &minT, &maxT);
        fprintf (stdout, "\n\t    Min (%4.2f) Max (%4.2f)\t", minT, maxT);
    }
    if ( max>(Flt)Max || gamma>1 ) {
        if ( gamma>1 )
            fprintf (stdout, "\n\t    Gamma (%d)\t", gamma);
        else
            fprintf (stdout, "\n\t    (No gamma correction)\t");
        fflush (stdout);
        for ( i=0 ; i<height ; i++ ) {
            for ( j=0 ; j<width ; j++ ) {
                A[i][j] = A[i][j]/(max-min); /* scale numbers between 0 and 1 */
                if ( gamma>1 && A[i][j]!=(Flt)0 )   /* apply gamma correction */
                    A[i][j] = (Flt)pow((Flt)A[i][j], (Flt)1/(Flt)gamma);
                A[i][j] *= (Flt)255;           /* scale numbers between 0 and 255 */
            }
        }
    }
    if ( gamma>1 ) {
        MinMax (A, width, height, &min, &max);
        fprintf (stdout, "\n\t    Min (%4.2f) Max (%4.2f)\t", min, max);
    }
}

/* Find min and max values of a given vector */
static void
MinMax ( const Matrix A, const int width, const int height, Flt *min, Flt *max )
{
    int i, j;

    *min = MAXFLOAT;
    *max = MINFLOAT;

    for ( i=0 ; i<height ; i++ ) {
        for ( j=0 ; j<width ; j++ ) {
            *min = MIN(*min, A[i][j]);
            *max = MAX(*max, A[i][j]);
            /*
            fprintf (stdout, "\nA[%d][%d](%g) minimum(%g) maximum(%g)\n",
                     i, j, A[i][j], *min, *max);
            */
        }
    }
}