cvemd.cpp

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

        if( min_val > temp )
        {
            leave_x = loop[i];
            min_val = temp;
        }
    }

    /* update the loop */
    for( i = 0; i < steps; i += 2 )
    {
        float temp0 = loop[i]->val + min_val;
        float temp1 = loop[i + 1]->val - min_val;

        loop[i]->val = temp0;
        loop[i + 1]->val = temp1;
    }

    /* remove the leaving basic variable */
    i = leave_x->i;
    j = leave_x->j;
    state->is_x[i][j] = 0;

    head.next[0] = state->rows_x[i];
    cur_x = &head;
    while( (next_x = cur_x->next[0]) != leave_x )
    {
        cur_x = next_x;
        assert( cur_x );
    }
    cur_x->next[0] = next_x->next[0];
    state->rows_x[i] = head.next[0];

    head.next[1] = state->cols_x[j];
    cur_x = &head;
    while( (next_x = cur_x->next[1]) != leave_x )
    {
        cur_x = next_x;
        assert( cur_x );
    }
    cur_x->next[1] = next_x->next[1];
    state->cols_x[j] = head.next[1];

    /* set enter_x to be the new empty slot */
    state->enter_x = leave_x;

    return CV_NO_ERR;
}



/****************************************************************************************\
*                                    icvFindLoop                                       *
\****************************************************************************************/
static int
icvFindLoop( CvEMDState * state )
{
    int i, steps = 1;
    CvNode2D *new_x;
    CvNode2D **loop = state->loop;
    CvNode2D *enter_x = state->enter_x, *_x = state->_x;
    char *is_used = state->is_used;

    memset( is_used, 0, state->ssize + state->dsize );

    new_x = loop[0] = enter_x;
    is_used[enter_x - _x] = 1;
    steps = 1;

    do
    {
        if( (steps & 1) == 1 )
        {
            /* find an unused x in the row */
            new_x = state->rows_x[new_x->i];
            while( new_x != 0 && is_used[new_x - _x] )
                new_x = new_x->next[0];
        }
        else
        {
            /* find an unused x in the column, or the entering x */
            new_x = state->cols_x[new_x->j];
            while( new_x != 0 && is_used[new_x - _x] && new_x != enter_x )
                new_x = new_x->next[1];
            if( new_x == enter_x )
                break;
        }

        if( new_x != 0 )        /* found the next x */
        {
            /* add x to the loop */
            loop[steps++] = new_x;
            is_used[new_x - _x] = 1;
        }
        else                    /* didn't find the next x */
        {
            /* backtrack */
            do
            {
                i = steps & 1;
                new_x = loop[steps - 1];
                do
                {
                    new_x = new_x->next[i];
                }
                while( new_x != 0 && is_used[new_x - _x] );

                if( new_x == 0 )
                {
                    is_used[loop[--steps] - _x] = 0;
                }
            }
            while( new_x == 0 && steps > 0 );

            is_used[loop[steps - 1] - _x] = 0;
            loop[steps - 1] = new_x;
            is_used[new_x - _x] = 1;
        }
    }
    while( steps > 0 );

    return steps;
}



/****************************************************************************************\
*                                        icvRussel                                     *
\****************************************************************************************/
static void
icvRussel( CvEMDState * state )
{
    int i, j, min_i = -1, min_j = -1;
    float min_delta, diff;
    CvNode1D u_head, *cur_u, *prev_u;
    CvNode1D v_head, *cur_v, *prev_v;
    CvNode1D *prev_u_min_i = 0, *prev_v_min_j = 0, *remember;
    CvNode1D *u = state->u, *v = state->v;
    int ssize = state->ssize, dsize = state->dsize;
    float eps = CV_EMD_EPS * state->max_cost;
    float **cost = state->cost;
    float **delta = state->delta;

    /* initialize the rows list (ur), and the columns list (vr) */
    u_head.next = u;
    for( i = 0; i < ssize; i++ )
    {
        u[i].next = u + i + 1;
    }
    u[ssize - 1].next = 0;

    v_head.next = v;
    for( i = 0; i < dsize; i++ )
    {
        v[i].val = -CV_EMD_INF;
        v[i].next = v + i + 1;
    }
    v[dsize - 1].next = 0;

    /* find the maximum row and column values (ur[i] and vr[j]) */
    for( i = 0; i < ssize; i++ )
    {
        float u_val = -CV_EMD_INF;
        float *cost_row = cost[i];

        for( j = 0; j < dsize; j++ )
        {
            float temp = cost_row[j];

            if( u_val < temp )
                u_val = temp;
            if( v[j].val < temp )
                v[j].val = temp;
        }
        u[i].val = u_val;
    }

    /* compute the delta matrix */
    for( i = 0; i < ssize; i++ )
    {
        float u_val = u[i].val;
        float *delta_row = delta[i];
        float *cost_row = cost[i];

        for( j = 0; j < dsize; j++ )
        {
            delta_row[j] = cost_row[j] - u_val - v[j].val;
        }
    }

    /* find the basic variables */
    do
    {
        /* find the smallest delta[i][j] */
        min_i = -1;
        min_delta = CV_EMD_INF;
        prev_u = &u_head;
        for( cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next )
        {
            i = (int)(cur_u - u);
            float *delta_row = delta[i];

            prev_v = &v_head;
            for( cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next )
            {
                j = (int)(cur_v - v);
                if( min_delta > delta_row[j] )
                {
                    min_delta = delta_row[j];
                    min_i = i;
                    min_j = j;
                    prev_u_min_i = prev_u;
                    prev_v_min_j = prev_v;
                }
                prev_v = cur_v;
            }
            prev_u = cur_u;
        }

        if( min_i < 0 )
            break;

        /* add x[min_i][min_j] to the basis, and adjust supplies and cost */
        remember = prev_u_min_i->next;
        icvAddBasicVariable( state, min_i, min_j, prev_u_min_i, prev_v_min_j, &u_head );

        /* update the necessary delta[][] */
        if( remember == prev_u_min_i->next )    /* line min_i was deleted */
        {
            for( cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next )
            {
                j = (int)(cur_v - v);
                if( cur_v->val == cost[min_i][j] )      /* column j needs updating */
                {
                    float max_val = -CV_EMD_INF;

                    /* find the new maximum value in the column */
                    for( cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next )
                    {
                        float temp = cost[cur_u - u][j];

                        if( max_val < temp )
                            max_val = temp;
                    }

                    /* if needed, adjust the relevant delta[*][j] */
                    diff = max_val - cur_v->val;
                    cur_v->val = max_val;
                    if( fabs( diff ) < eps )
                    {
                        for( cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next )
                            delta[cur_u - u][j] += diff;
                    }
                }
            }
        }
        else                    /* column min_j was deleted */
        {
            for( cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next )
            {
                i = (int)(cur_u - u);
                if( cur_u->val == cost[i][min_j] )      /* row i needs updating */
                {
                    float max_val = -CV_EMD_INF;

                    /* find the new maximum value in the row */
                    for( cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next )
                    {
                        float temp = cost[i][cur_v - v];

                        if( max_val < temp )
                            max_val = temp;
                    }

                    /* if needed, adjust the relevant delta[i][*] */
                    diff = max_val - cur_u->val;
                    cur_u->val = max_val;

                    if( fabs( diff ) < eps )
                    {
                        for( cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next )
                            delta[i][cur_v - v] += diff;
                    }
                }
            }
        }
    }
    while( u_head.next != 0 || v_head.next != 0 );
}



/****************************************************************************************\
*                                   icvAddBasicVariable                                *
\****************************************************************************************/
static void
icvAddBasicVariable( CvEMDState * state,
                     int min_i, int min_j,
                     CvNode1D * prev_u_min_i, CvNode1D * prev_v_min_j, CvNode1D * u_head )
{
    float temp;
    CvNode2D *end_x = state->end_x;

    if( state->s[min_i] < state->d[min_j] + state->weight * CV_EMD_EPS )
    {                           /* supply exhausted */
        temp = state->s[min_i];
        state->s[min_i] = 0;
        state->d[min_j] -= temp;
    }
    else                        /* demand exhausted */
    {
        temp = state->d[min_j];
        state->d[min_j] = 0;
        state->s[min_i] -= temp;
    }

    /* x(min_i,min_j) is a basic variable */
    state->is_x[min_i][min_j] = 1;

    end_x->val = temp;
    end_x->i = min_i;
    end_x->j = min_j;
    end_x->next[0] = state->rows_x[min_i];
    end_x->next[1] = state->cols_x[min_j];
    state->rows_x[min_i] = end_x;
    state->cols_x[min_j] = end_x;
    state->end_x = end_x + 1;

    /* delete supply row only if the empty, and if not last row */
    if( state->s[min_i] == 0 && u_head->next->next != 0 )
        prev_u_min_i->next = prev_u_min_i->next->next;  /* remove row from list */
    else
        prev_v_min_j->next = prev_v_min_j->next->next;  /* remove column from list */
}


/****************************************************************************************\
*                                  standard  metrics                                     *
\****************************************************************************************/
static float
icvDistL1( const float *x, const float *y, void *user_param )
{
    int i, dims = (int)(size_t)user_param;
    double s = 0;

    for( i = 0; i < dims; i++ )
    {
        double t = x[i] - y[i];

        s += fabs( t );
    }
    return (float)s;
}

static float
icvDistL2( const float *x, const float *y, void *user_param )
{
    int i, dims = (int)(size_t)user_param;
    double s = 0;

    for( i = 0; i < dims; i++ )
    {
        double t = x[i] - y[i];

        s += t * t;
    }
    return cvSqrt( (float)s );
}

static float
icvDistC( const float *x, const float *y, void *user_param )
{
    int i, dims = (int)(size_t)user_param;
    double s = 0;

    for( i = 0; i < dims; i++ )
    {
        double t = fabs( x[i] - y[i] );

        if( s < t )
            s = t;
    }
    return (float)s;
}

/* End of file. */