cvemd.cpp

opencv库在TI DM6437上的移植,目前包括两个库cv.lib和cxcore.lib的工程

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

/*
    Partially based on Yossi Rubner code:
    =========================================================================
    emd.c

    Last update: 3/14/98

    An implementation of the Earth Movers Distance.
    Based of the solution for the Transportation problem as described in
    "Introduction to Mathematical Programming" by F. S. Hillier and 
    G. J. Lieberman, McGraw-Hill, 1990.

    Copyright (C) 1998 Yossi Rubner
    Computer Science Department, Stanford University
    E-Mail: [e-mail]   URL: http://vision.stanford.edu/~rubner
    ==========================================================================
*/
#include "_cv.h"

#define MAX_ITERATIONS 500
#define CV_EMD_INF   ((float)1e20)
#define CV_EMD_EPS   ((float)1e-5)

/* CvNode1D is used for lists, representing 1D sparse array */
typedef struct CvNode1D
{
    float val;
    struct CvNode1D *next;
}
CvNode1D;

/* CvNode2D is used for lists, representing 2D sparse matrix */
typedef struct CvNode2D
{
    float val;
    struct CvNode2D *next[2];  /* next row & next column */
    int i, j;
}
CvNode2D;


typedef struct CvEMDState
{
    int ssize, dsize;

    float **cost;
    CvNode2D *_x;
    CvNode2D *end_x;
    CvNode2D *enter_x;
    char **is_x;

    CvNode2D **rows_x;
    CvNode2D **cols_x;

    CvNode1D *u;
    CvNode1D *v;

    int* idx1;
    int* idx2;

    /* find_loop buffers */
    CvNode2D **loop;
    char *is_used;

    /* russel buffers */
    float *s;
    float *d;
    float **delta;

    float weight, max_cost;
    char *buffer;
}
CvEMDState;

/* static function declaration */
static CvStatus icvInitEMD( const float *signature1, int size1,
                            const float *signature2, int size2,
                            int dims, CvDistanceFunction dist_func, void *user_param,
                            const float* cost, int cost_step,
                            CvEMDState * state, float *lower_bound,
                            char *local_buffer, int local_buffer_size );

static CvStatus icvFindBasicVariables( float **cost, char **is_x,
                                       CvNode1D * u, CvNode1D * v, int ssize, int dsize );

static float icvIsOptimal( float **cost, char **is_x,
                           CvNode1D * u, CvNode1D * v,
                           int ssize, int dsize, CvNode2D * enter_x );

static void icvRussel( CvEMDState * state );


static CvStatus icvNewSolution( CvEMDState * state );
static int icvFindLoop( CvEMDState * state );

static void icvAddBasicVariable( CvEMDState * state,
                                 int min_i, int min_j,
                                 CvNode1D * prev_u_min_i,
                                 CvNode1D * prev_v_min_j,
                                 CvNode1D * u_head );

static float icvDistL2( const float *x, const float *y, void *user_param );
static float icvDistL1( const float *x, const float *y, void *user_param );
static float icvDistC( const float *x, const float *y, void *user_param );

/* The main function */
CV_IMPL float
cvCalcEMD2( const CvArr* signature_arr1,
            const CvArr* signature_arr2,
            int dist_type,
            CvDistanceFunction dist_func,
            const CvArr* cost_matrix,
            CvArr* flow_matrix,
            float *lower_bound,
            void *user_param )
{
    char local_buffer[16384];
    char *local_buffer_ptr = (char *)cvAlignPtr(local_buffer,16);
    CvEMDState state;
    float emd = 0;

    CV_FUNCNAME( "cvCalcEMD2" );

    memset( &state, 0, sizeof(state));

    __BEGIN__;

    double total_cost = 0;
    CvStatus result = CV_NO_ERR;
    float eps, min_delta;
    CvNode2D *xp = 0;
    CvMat sign_stub1, *signature1 = (CvMat*)signature_arr1;
    CvMat sign_stub2, *signature2 = (CvMat*)signature_arr2;
    CvMat cost_stub, *cost = &cost_stub;
    CvMat flow_stub, *flow = (CvMat*)flow_matrix;
    int dims, size1, size2;

    CV_CALL( signature1 = cvGetMat( signature1, &sign_stub1 ));
    CV_CALL( signature2 = cvGetMat( signature2, &sign_stub2 ));

    if( signature1->cols != signature2->cols )
        CV_ERROR( CV_StsUnmatchedSizes, "The arrays must have equal number of columns (which is number of dimensions but 1)" );

    dims = signature1->cols - 1;
    size1 = signature1->rows;
    size2 = signature2->rows;

    if( !CV_ARE_TYPES_EQ( signature1, signature2 ))
        CV_ERROR( CV_StsUnmatchedFormats, "The array must have equal types" );

    if( CV_MAT_TYPE( signature1->type ) != CV_32FC1 )
        CV_ERROR( CV_StsUnsupportedFormat, "The signatures must be 32fC1" );

    if( flow )
    {
        CV_CALL( flow = cvGetMat( flow, &flow_stub ));

        if( flow->rows != size1 || flow->cols != size2 )
            CV_ERROR( CV_StsUnmatchedSizes,
            "The flow matrix size does not match to the signatures' sizes" );

        if( CV_MAT_TYPE( flow->type ) != CV_32FC1 )
            CV_ERROR( CV_StsUnsupportedFormat, "The flow matrix must be 32fC1" );
    }

    cost->data.fl = 0;
    cost->step = 0;

    if( dist_type < 0 )
    {
        if( cost_matrix )
        {
            if( dist_func )
                CV_ERROR( CV_StsBadArg,
                "Only one of cost matrix or distance function should be non-NULL in case of user-defined distance" );

            if( lower_bound )
                CV_ERROR( CV_StsBadArg,
                "The lower boundary can not be calculated if the cost matrix is used" );

            CV_CALL( cost = cvGetMat( cost_matrix, &cost_stub ));
            if( cost->rows != size1 || cost->cols != size2 )
                CV_ERROR( CV_StsUnmatchedSizes,
                "The cost matrix size does not match to the signatures' sizes" );

            if( CV_MAT_TYPE( cost->type ) != CV_32FC1 )
                CV_ERROR( CV_StsUnsupportedFormat, "The cost matrix must be 32fC1" );
        }
        else if( !dist_func )
            CV_ERROR( CV_StsNullPtr, "In case of user-defined distance Distance function is undefined" );
    }
    else
    {
        if( dims == 0 )
            CV_ERROR( CV_StsBadSize,
            "Number of dimensions can be 0 only if a user-defined metric is used" );
        user_param = (void *) (size_t)dims;
        switch (dist_type)
        {
        case CV_DIST_L1:
            dist_func = (CvDistanceFunction)icvDistL1;
            break;
        case CV_DIST_L2:
            dist_func = (CvDistanceFunction)icvDistL2;
            break;
        case CV_DIST_C:
            dist_func = (CvDistanceFunction)icvDistC;
            break;
        default:
            CV_ERROR( CV_StsBadFlag, "Bad or unsupported metric type" );
        }
    }

    IPPI_CALL( result = icvInitEMD( signature1->data.fl, size1,
                                    signature2->data.fl, size2,
                                    dims, dist_func, user_param,
                                    cost->data.fl, cost->step,
                                    &state, lower_bound, local_buffer_ptr,
                                    sizeof( local_buffer ) - 16 ));

    if( result > 0 && lower_bound )
    {
        emd = *lower_bound;
        EXIT;
    }

    eps = CV_EMD_EPS * state.max_cost;

    /* if ssize = 1 or dsize = 1 then we are done, else ... */
    if( state.ssize > 1 && state.dsize > 1 )
    {
        int itr;

        for( itr = 1; itr < MAX_ITERATIONS; itr++ )
        {
            /* find basic variables */
            result = icvFindBasicVariables( state.cost, state.is_x,
                                            state.u, state.v, state.ssize, state.dsize );
            if( result < 0 )
                break;

            /* check for optimality */
            min_delta = icvIsOptimal( state.cost, state.is_x,
                                      state.u, state.v,
                                      state.ssize, state.dsize, state.enter_x );

            if( min_delta == CV_EMD_INF )
            {
                CV_ERROR( CV_StsNoConv, "" );
            }

            /* if no negative deltamin, we found the optimal solution */
            if( min_delta >= -eps )
                break;

            /* improve solution */
            IPPI_CALL( icvNewSolution( &state ));
        }
    }

    /* compute the total flow */
    for( xp = state._x; xp < state.end_x; xp++ )
    {
        float val = xp->val;
        int i = xp->i;
        int j = xp->j;
        int ci = state.idx1[i];
        int cj = state.idx2[j];

        if( xp != state.enter_x && ci >= 0 && cj >= 0 )
        {
            total_cost += (double)val * state.cost[i][j];
            if( flow )
                ((float*)(flow->data.ptr + flow->step*ci))[cj] = val;
        }
    }

    emd = (float) (total_cost / state.weight);

    __END__;

    if( state.buffer && state.buffer != local_buffer_ptr )
        cvFree( &state.buffer );

    return emd;
}


/************************************************************************************\
*          initialize structure, allocate buffers and generate initial golution      *
\************************************************************************************/
static CvStatus
icvInitEMD( const float* signature1, int size1,
            const float* signature2, int size2,
            int dims, CvDistanceFunction dist_func, void* user_param,
            const float* cost, int cost_step,
            CvEMDState* state, float* lower_bound,
            char* local_buffer, int local_buffer_size )
{
    float s_sum = 0, d_sum = 0, diff;
    int i, j;
    int ssize = 0, dsize = 0;
    int equal_sums = 1;
    int buffer_size;
    float max_cost = 0;
    char *buffer, *buffer_end;

    memset( state, 0, sizeof( *state ));
    assert( cost_step % sizeof(float) == 0 );
    cost_step /= sizeof(float);

    /* calculate buffer size */
    buffer_size = (size1+1) * (size2+1) * (sizeof( float ) +    /* cost */
                                   sizeof( char ) +     /* is_x */
                                   sizeof( float )) +   /* delta matrix */
        (size1 + size2 + 2) * (sizeof( CvNode2D ) + /* _x */
                           sizeof( CvNode2D * ) +  /* cols_x & rows_x */
                           sizeof( CvNode1D ) + /* u & v */
                           sizeof( float ) + /* s & d */
                           sizeof( int ) + sizeof(CvNode2D*)) +  /* idx1 & idx2 */ 
        (size1+1) * (sizeof( float * ) + sizeof( char * ) + /* rows pointers for */
                 sizeof( float * )) + 256;      /*  cost, is_x and delta */

    if( buffer_size < (int) (dims * 2 * sizeof( float )))
    {
        buffer_size = dims * 2 * sizeof( float );
    }

    /* allocate buffers */
    if( local_buffer != 0 && local_buffer_size >= buffer_size )
    {
        buffer = local_buffer;
    }
    else
    {
        buffer = (char*)cvAlloc( buffer_size );
        if( !buffer )
            return CV_OUTOFMEM_ERR;
    }

    state->buffer = buffer;
    buffer_end = buffer + buffer_size;

    state->idx1 = (int*) buffer;
    buffer += (size1 + 1) * sizeof( int );

    state->idx2 = (int*) buffer;
    buffer += (size2 + 1) * sizeof( int );