cost.c

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

/*
 *  -*- Mode: ANSI C -*-
 *  $Id: cost.c,v 1.10 1996/09/17 16:10:36 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/cost.c,v $
 *  Author: Gabriel Fernandez
 *
 *  Declaration of cost functions for calculation of best basis.
 */
/* do not edit anything above this line */

/* System header files */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* FLWT header files */
#include <cost.h>
#include <flwterr.h>
#include <flwt.h>
#include <flwtpack.h>
#include <imgmem.h>
#include <mallat.h>
#include <util.h>

/* Pointer to the cost function */
Flt (*infocost)(Vector, const int, const int);
/* Epsilon for the threshold cost function */
static Flt epsilon = (Flt)50;
/* Power for the lpnormp cost function */
static Flt p = (Flt)2.0;

/* code */

/*
 * SetInfoCost function: sets the pointer to the information cost function
 *                       to be used in the calculation of the best packet
 *                       basis. If the value of opt is NULL, then the menu
 *                       of options is displayed. For threshold and lpnorm,
 *                       the values of epsilon and p can be given. A value
 *                       of 999.0 indicates that they have to asked.
 */
void
SetInfoCost ( const int option, const Flt value, const boolean print )
{
    boolean exit = FALSE;
    char input[BUFSIZ];
    int opt;

    /* Is the menu required? */
    if ( option == 0 ) {
        fprintf (stdout, "Select a cost function:\n");
        fprintf (stdout, "a) Aditive.-\n");
        fprintf (stdout, "    [1] Threshold\n");
        fprintf (stdout, "    [2] l1norm\n");
        fprintf (stdout, "    [3] lpnormp\n");
        fprintf (stdout, "    [4] Entropy (ml2logl2)\n");
        fprintf (stdout, "    [5] Distortion (m.s.e.)\n");
        fprintf (stdout, "b) Not Aditive.-\n");
        fprintf (stdout, "    [6] Gauss-Markov (logl2)\n");
        fflush (stdout);
        do {
            fprintf (stdout, "\nOption: ");
            gets (input);
            fflush (stdin);
            opt = (int)atoi(input);
            if ( opt<1 || opt>6 ) {
                Error ("SetInfoCost", INVALID_VALUE_INT, RETURN, "Option", 1, 6);
                exit = FALSE;
            } else {
                exit = TRUE;
            }
        } while (!exit);
    } else {
	opt = option;
    }

    /* Set info cost function pointer */
    switch (opt) {
        case 1:   /* ask for threshold epsilon */
            if ( value != (Flt)999 ) {
		epsilon = value;
	    } else {
                fprintf (stdout, "\nIntroduce the epsilon value: ");
                fflush (stdout);
                gets (input);
	        fflush (stdin);
                epsilon = (Flt)atof(input);
            }
            infocost = &thresh;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [epsilon(%g)]\n",
                         epsilon);
            break;
        case 2:
            infocost = &l1norm;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [l1norm]\n");
            break;
        case 3:   /* ask for pth power for lpnormp */
	    if ( p != (Flt)999 ) {
		p = value;
	    } else {
                do {
                    fprintf (stdout,
		             "\nIntroduce the value of p "
			     "(suggested value is %g): ", p);
                    fflush (stdout);
                    gets (input);
                    fflush (stdin);
                    p = (Flt)atof(input);
                    if ( p <= (Flt)0 ) {
                        Error ("SetInfoCost", INVALID_RELATIONAL_FLT, RETURN,
		               "P", ">", 0);
                        exit = FALSE;
                    } else {
                        exit = TRUE;
                    }
                } while (!exit);
            }
            infocost = &lpnormp;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [lpnormp(p=%g)]\n", p);
            break;
        case 4:
            infocost = &ml2logl2;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [entropy(ml2logl2)]\n");
            break;
        case 5:
            infocost = &mse;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [distortion(m.s.e.)]\n");
            break;
        case 6:
            infocost = &logl2;
            if ( print )
                fprintf (stdout, "\nInformation Cost function [Gauss-Markov(logl2)]\n");
    }
    if ( print )
        fprintf (stdout, "\n");
    fflush (stdout);
}

/**************************/
/* Aditive cost functions */
/**************************/

/*
 * Thresh function: measures the cost of a sequence of coefficients by the
 *                  "number exceeding a threshold".
 */
Flt
thresh ( Vector data, const int least, const int final )
{
    Flt cost;
    int k;

    /* Find threshold 
    Te = epsilon * exp( -ml2logl2( data, least, final ) / lpnormp( data, least, final ) );
    fprintf (stdout, "\nTe(%g)\n", Te);
     */

    /* Threshold information cost functional */
    cost = (Flt)0;
    for ( k=least; k<final ; k++ )
        if ( ABS(data[k]) > epsilon )
            cost += (Flt)1;

    return cost;
}

/*
 * lpnormp functions: calculate the concentration in l^p by summing the pth
 *                    powers of their coefficients, for 0<p<2.
 */
Flt
l1norm ( Vector data, const int least, const int final )
{
    Flt cost;
    int k;

    /* l^1 information cost functional */
    cost = (Flt)0;
    for ( k=least ; k<final ; k++ )
        cost += ABS(data[k]);

    return cost;
}

Flt
lpnormp ( Vector data, const int least, const int final )
{
    Flt cost;
    int k;

    /* l^p information cost functional */
    cost = (Flt)0;
    for ( k=least ; k<final ; k++ ) {
        Flt absdata = ABS(data[k]);
        if ( absdata > (Flt)0 )
            cost += (Flt)exp( p*(Flt)log((Flt)absdata) );
    }

    return cost;
}

/*
 * ml2logl2 function: entropy itself is not an information cost functional, since
 *                    it requires that the sequence be normalized and is thus not
 *                    additive. Instead, we use the -(l^2)log(l^2) "norm" which does
 *                    satisfy the additivity condition required of an informational
 *                    cost functional and is monotonic with the entropy of a sequence.
 */
Flt
ml2logl2 ( Vector data, const int least, const int final )
{
    Flt cost;
    int k;

    /* -(l^2)log(l^2) or "entropy" information cost functional */
    cost = (Flt)0;
    for ( k=least ; k<final ; k++ ) {
        Flt dataSQR = data[k] * data[k];
        if ( dataSQR > (Flt)0 )
            cost -= dataSQR * (Flt)log((Flt)dataSQR);
    }

    return cost;
}


/*
 * distortion function: finds the amount of distorion of the given area using
 *                      the Mean Squared Error of the signal.
 */
Flt
mse ( Vector data, const int least, const int final )
{
    Flt cost;
    int i;

    cost = (Flt)0;
    for ( i=least ; i<final ; i++ )
        cost += data[i] * data[i];
    cost /= (Flt)(final-least+1);

    return cost;
}


/*******************************/
/* Not Additive Cost functions */
/*******************************/

/*
 * logl2 function: measures the entropy of a Gauss-Markov process as the sum of
 *                 the "logarithms of the variances" of its components. If the
 *                 components have zero mean, then these variances may be
 *                 approximated by the energies of the components.
 */
Flt
logl2 ( Flt *data, const int least, const int final )
{
    Flt cost;
    int k;

    /* log(l^2) or Gauss-Markov information cost functional */
    cost = (Flt)0;
    for ( k=least ; k<final ; k++ )
        cost += (Flt)log( (Flt)(data[k] * data[k]) );

    return cost;
}


/*************************/
/* Theoretical dimension */
/*************************/

/*
 * tdim function: using the function ml2logl2() defined above, the "theoretical
 *                dimension" of a sequence is found. Both the energy and the
 *                (l^2)log(l^2) "norm" of the sequence are needed, and they are
 *                computed at once to avoid having to square each array element
 *                twice.
 */ 
Flt
tdim ( Flt *data, const int least, const int final )
{
    Flt md2logd2, energy, theodim;
    int k;

    /* Theoretical dimension of a sequence */
    md2logd2 = (Flt)0;
    energy  = (Flt)0;
    for ( k=least ; k<final ; k++ ) {
        Flt dSQR = data[k] * data[k];
        if ( dSQR > (Flt)0 ) {
            md2logd2 -= dSQR * (Flt)log((Flt)dSQR);
            energy += dSQR;
        }
    }
    if ( energy > (Flt)0 )
        theodim = energy * (Flt)exp( (Flt)(md2logd2 / energy) );
    else
        theodim = (Flt)0;

    return theodim;
}


/**************/
/* Best basis */
/**************/

/*
 * costs2bbasis function: finds the best basis hedge from a BTN tree with cost tags.
 */
Flt
costs2bbasis ( Hedge *graph, BTN *root, const int level )
{
    Flt *bestcost;

    if ( root->left == NULL && root->right ) {
        graph->levels[graph->blocks] = level;
        graph->contents[graph->blocks] = root->content;
        graph->blocks++;
        bestcost = (Flt *)root->tag;
    } else {
        int blocks = graph->blocks;
        Flt cost = (Flt)0, *tag;

        if ( root->left != NULL )
            cost += costs2bbasis( graph, root->left, level+1 );
        if ( root->right != NULL )
            cost += costs2bbasis( graph, root->right, level+1 );
        tag = (Flt *)root->tag;
        if ( *tag > cost )
            bestcost = &cost;
        else {
            bestcost = (Flt *)root->tag;
            graph->blocks = blocks;
            graph->levels[graph->blocks] = level;
            graph->contents[graph->blocks] = root->content;
            graph->blocks++;
         }
    }

    return *bestcost;
}

/*
 * btnt2costs function: adds costs tags to a BTN tree of intervals.
 */
void
btnt2costs ( BTN *root )
{
    if ( root != NULL ) {
        Interval *I = root->content;
        Flt cost = (Flt)(*infocost)( I->origin, I->least, I->final );
        root->tag = &cost;
        btnt2costs( root->left );
        btnt2costs( root->right );
    }

    return;
}

/*
 * bbwp2 function: best basis of separable two-dimensional periodic wavelet packets.
 *                 The input consists of a preallocated Hedge data structure to hold
 *                 the graph, with its first contents member poiting to a preallocated
 *                 output array of length x*y; the number of columns x; the number of
 *                 rows y; the number of vanishing moments N and nTilde; the current
 *                 level s; and the maximum level l.
 *                 The contents array of the hedge should be allocated with (2^l)+1
 *                 memory locations, one more than the maximum number of blocks in
 *                 any l-level graph basis, since it will get one extra pointer to the
 *                 first location after the last block of coefficients. graph.contents[0]
 *                 should be an array with at least x*y locations, since it will be
 *                 filled with the best basis coefficients.
 */
Flt
bbwp2 ( Hedge *graph, Matrix in, const int x, const int y, const N, const nTilde,
                                 const int s, const int l )
{
    register int i, j;
    int xy;
    Flt cost;
    Vector out, inPtr;

    xy = x*y;
    cost = (Flt)(*infocost)( &in[0][0], 0, xy-1 );
    if ( s<l ) {
        int block = graph->blocks;
        Matrix *k = DSFLWT2D( in, x, y, N, nTilde );   /* separate them */
        Flt kcost = (Flt)0;
        for ( j=0 ; j<4 ; j++ )
            kcost += bbwp2( graph, k[j], x>>1, y>>1, N, nTilde, s+1, l );
        if ( kcost < cost ) {
            cost = kcost;
        } else {
            Matrix parent = in;
            /* Get parent back */
            IFLWTP2D( parent, x, y, N, nTilde );
            graph->levels[block] = s;
            out = (Vector)graph->contents[block];
            for ( i=0, inPtr=&parent[0][0] ; i<xy ; i++, inPtr++ ) {
                out[i] = *inPtr;
            }
            graph->blocks = block + 1;
            graph->contents[graph->blocks] = out + xy;
        }
        /* Free temporary memory */
        MEM_Free3D( k, 4 );
    } else {
        graph->levels[graph->blocks] = l;
        out = (Vector)graph->contents[graph->blocks];
        for ( i=0, inPtr=&in[0][0] ; i<xy ; i++, inPtr++ ) {
            out[i] = *inPtr;
        }
        graph->blocks++;
        graph->contents[graph->blocks] = out + xy;
    }

    return cost;
} 

/*
 * bbwp1 function: best basis of separable one-dimensional periodic
 *                 wavelet packets. The input consists of a preallocated 
 *                 Hedge data structure to hold the graph, with its 
 *                 first contents member pointing to a preallocated
 *                 output array of length x; the number of columns x;
 *                 the number of vanishing moments N and nTilde; the
 *                 current level s; and the maximum level l.
 *                 The contents array of the hedge should be allocated
 *                 with (2^l)+1 memory locations, one more than the
 *                 maximum number of blocks in any l-level graph basis,
 *                 since it will get one extra pointer to the  first
 *                 location after the last block of coefficients.
 *                 graph.contents[0] should be an array with at least
 *                 x locations, since it will be filled with the best
 *                 basis coefficients.
 */
Flt
bbwp1 ( Hedge *graph, Vector in, const int x, const N, const nTilde,
                                 const int s, const int l )
{
    register int i, j;
    Flt cost;
    Vector out;

    cost = (Flt)(*infocost)( in, 0, x-1 );
    if ( s<l ) {
        int block = graph->blocks;
        Vector *k = DSFLWT1D( in, x, N, nTilde );

        Flt kcost = (Flt)0;
        for ( j=0 ; j<2 ; j++ )
            kcost += bbwp1( graph, k[j], x>>1, N, nTilde, s+1, l );

        if ( kcost < cost ) {
            cost = kcost;
        } else {
            Vector parent = in;
            IFLWTP1D( parent, x, N, nTilde );
            graph->levels[block] = s;  /* set the lvl for this block */
            out = (Vector)graph->contents[block];
            for ( i=0 ; i<x ; i++ )
                out[i] = parent[i];
            graph->blocks = block + 1; /* there is now one more block */
            /* set contents[] so it points to next block */
            graph->contents[graph->blocks] = out + x;
        }
        MEM_FreeVecs( k, 2 );
    } else {
        graph->levels[graph->blocks] = l;
        out = (Vector)graph->contents[graph->blocks];
        for ( i=0 ; i<x ; i++ )
            out[i] = in[i];
        graph->blocks++;
        graph->contents[graph->blocks] = out + x;
    }

    return cost;
}