maxflow.cpp

markov random field in matlab code

                    if (orphan_last) orphan_last -> next = np;
                    else             orphan_first        = np;
                    orphan_last = np;
                    np -> next = NULL;
                }
            }
        }
        for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
        {
            a0_for = a0_rev -> sister;
            j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
            if (!j->is_sink && (a=j->parent))
            {
                if (a0_for->r_cap) set_active(j);
                if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
                {
                    /* add j to the adoption list */
                    j -> parent = ORPHAN;
                    np = nodeptr_block -> New();
                    np -> ptr = j;
                    if (orphan_last) orphan_last -> next = np;
                    else             orphan_first        = np;
                    orphan_last = np;
                    np -> next = NULL;
                }
            }
        }
    }
}

void Graph::process_sink_orphan(node *i)
{
    node *j;
    arc_forward *a0_for, *a0_for_first, *a0_for_last;
    arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
    arc_forward *a0_min = NULL, *a;
    nodeptr *np;
    int d, d_min = INFINITE_D;

    /* trying to find a new parent */
    a0_for_first = i -> first_out;
    if (IS_ODD(a0_for_first))
    {
        a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
        a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
    }
    else a0_for_last = (i + 1) -> first_out;
    a0_rev_first = i -> first_in;
    if (IS_ODD(a0_rev_first))
    {
        a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
        a0_rev_last  = (arc_reverse *) ((a0_rev_first ++) -> sister);
    }
    else a0_rev_last = (i + 1) -> first_in;


    for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
    if (a0_for->r_cap)
    {
        j = NEIGHBOR_NODE(i, a0_for -> shift);
        if (j->is_sink && (a=j->parent))
        {
            /* checking the origin of j */
            d = 0;
            while ( 1 )
            {
                if (j->TS == TIME)
                {
                    d += j -> DIST;
                    break;
                }
                a = j -> parent;
                d ++;
                if (a==TERMINAL)
                {
                    j -> TS = TIME;
                    j -> DIST = 1;
                    break;
                }
                if (a==ORPHAN) { d = INFINITE_D; break; }
                if (IS_ODD(a))
                    j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                else
                    j = NEIGHBOR_NODE(j, a -> shift);
            }
            if (d<INFINITE_D) /* j originates from the sink - done */
            {
                if (d<d_min)
                {
                    a0_min = a0_for;
                    d_min = d;
                }
                /* set marks along the path */
                for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; )
                {
                    j -> TS = TIME;
                    j -> DIST = d --;
                    a = j->parent;
                    if (IS_ODD(a))
                        j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                    else
                        j = NEIGHBOR_NODE(j, a -> shift);
                }
            }
        }
    }
    for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
    {
        a0_for = a0_rev -> sister;
        if (a0_for->r_rev_cap)
        {
            j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
            if (j->is_sink && (a=j->parent))
            {
                /* checking the origin of j */
                d = 0;
                while ( 1 )
                {
                    if (j->TS == TIME)
                    {
                        d += j -> DIST;
                        break;
                    }
                    a = j -> parent;
                    d ++;
                    if (a==TERMINAL)
                    {
                        j -> TS = TIME;
                        j -> DIST = 1;
                        break;
                    }
                    if (a==ORPHAN) { d = INFINITE_D; break; }
                    if (IS_ODD(a))
                        j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                    else
                        j = NEIGHBOR_NODE(j, a -> shift);
                }
                if (d<INFINITE_D) /* j originates from the sink - done */
                {
                    if (d<d_min)
                    {
                        a0_min = MAKE_ODD(a0_for);
                        d_min = d;
                    }
                    /* set marks along the path */
                    for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; )
                    {
                        j -> TS = TIME;
                        j -> DIST = d --;
                        a = j->parent;
                        if (IS_ODD(a))
                            j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                        else
                            j = NEIGHBOR_NODE(j, a -> shift);
                    }
                }
            }
        }
    }

    if ((i->parent = a0_min))
    {
        i -> TS = TIME;
        i -> DIST = d_min + 1;
    }
    else
    {
        /* no parent is found */
        i -> TS = 0;

        /* process neighbors */
        for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
        {
            j = NEIGHBOR_NODE(i, a0_for -> shift);
            if (j->is_sink && (a=j->parent))
            {
                if (a0_for->r_cap) set_active(j);
                if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)
                {
                    /* add j to the adoption list */
                    j -> parent = ORPHAN;
                    np = nodeptr_block -> New();
                    np -> ptr = j;
                    if (orphan_last) orphan_last -> next = np;
                    else             orphan_first        = np;
                    orphan_last = np;
                    np -> next = NULL;
                }
            }
        }
        for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
        {
            a0_for = a0_rev -> sister;
            j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
            if (j->is_sink && (a=j->parent))
            {
                if (a0_for->r_rev_cap) set_active(j);
                if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
                {
                    /* add j to the adoption list */
                    j -> parent = ORPHAN;
                    np = nodeptr_block -> New();
                    np -> ptr = j;
                    if (orphan_last) orphan_last -> next = np;
                    else             orphan_first        = np;
                    orphan_last = np;
                    np -> next = NULL;
                }
            }
        }
    }
}

/***********************************************************************/

Graph::flowtype Graph::maxflow()
{
    node *i, *j, *current_node = NULL, *s_start, *t_start = NULL;
    captype *cap_middle = NULL, *rev_cap_middle = NULL;
    arc_forward *a_for, *a_for_first, *a_for_last;
    arc_reverse *a_rev, *a_rev_first, *a_rev_last;
    nodeptr *np, *np_next;

    prepare_graph();
    maxflow_init();
    nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);

    while ( 1 )
    {
        if ((i=current_node))
        {
            i -> next = NULL; /* remove active flag */
            if (!i->parent) i = NULL;
        }
        if (!i)
        {
            if (!(i = next_active())) break;
        }

        /* growth */
        s_start = NULL;

        a_for_first = i -> first_out;
        if (IS_ODD(a_for_first))
        {
            a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
            a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
        }
        else a_for_last = (i + 1) -> first_out;
        a_rev_first = i -> first_in;
        if (IS_ODD(a_rev_first))
        {
            a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);
            a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);
        }
        else a_rev_last = (i + 1) -> first_in;

        if (!i->is_sink)
        {
            /* grow source tree */
            for (a_for=a_for_first; a_for<a_for_last; a_for++)
            if (a_for->r_cap)
            {
                j = NEIGHBOR_NODE(i, a_for -> shift);
                if (!j->parent)
                {
                    j -> is_sink = 0;
                    j -> parent = MAKE_ODD(a_for);
                    j -> TS = i -> TS;
                    j -> DIST = i -> DIST + 1;
                    set_active(j);
                }
                else if (j->is_sink)
                {
                    s_start = i;
                    t_start = j;
                    cap_middle     = & ( a_for -> r_cap );
                    rev_cap_middle = & ( a_for -> r_rev_cap );
                    break;
                }
                else if (j->TS <= i->TS &&
                         j->DIST > i->DIST)
                {
                    /* heuristic - trying to make the distance from j to the source shorter */
                    j -> parent = MAKE_ODD(a_for);
                    j -> TS = i -> TS;
                    j -> DIST = i -> DIST + 1;
                }
            }
            if (!s_start)
            for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)
            {
                a_for = a_rev -> sister;
                if (a_for->r_rev_cap)
                {
                    j = NEIGHBOR_NODE_REV(i, a_for -> shift);
                    if (!j->parent)
                    {
                        j -> is_sink = 0;
                        j -> parent = a_for;
                        j -> TS = i -> TS;
                        j -> DIST = i -> DIST + 1;
                        set_active(j);
                    }
                    else if (j->is_sink)
                    {
                        s_start = i;
                        t_start = j;
                        cap_middle     = & ( a_for -> r_rev_cap );
                        rev_cap_middle = & ( a_for -> r_cap );
                        break;
                    }
                    else if (j->TS <= i->TS &&
                             j->DIST > i->DIST)
                    {
                        /* heuristic - trying to make the distance from j to the source shorter */
                        j -> parent = a_for;
                        j -> TS = i -> TS;
                        j -> DIST = i -> DIST + 1;
                    }
                }
            }
        }
        else
        {
            /* grow sink tree */
            for (a_for=a_for_first; a_for<a_for_last; a_for++)
            if (a_for->r_rev_cap)
            {
                j = NEIGHBOR_NODE(i, a_for -> shift);
                if (!j->parent)
                {
                    j -> is_sink = 1;
                    j -> parent = MAKE_ODD(a_for);
                    j -> TS = i -> TS;
                    j -> DIST = i -> DIST + 1;
                    set_active(j);
                }
                else if (!j->is_sink)
                {
                    s_start = j;
                    t_start = i;
                    cap_middle     = & ( a_for -> r_rev_cap );
                    rev_cap_middle = & ( a_for -> r_cap );
                    break;
                }
                else if (j->TS <= i->TS &&
                         j->DIST > i->DIST)
                {
                    /* heuristic - trying to make the distance from j to the sink shorter */
                    j -> parent = MAKE_ODD(a_for);
                    j -> TS = i -> TS;
                    j -> DIST = i -> DIST + 1;
                }
            }
            for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)
            {
                a_for = a_rev -> sister;
                if (a_for->r_cap)
                {
                    j = NEIGHBOR_NODE_REV(i, a_for -> shift);
                    if (!j->parent)
                    {
                        j -> is_sink = 1;
                        j -> parent = a_for;
                        j -> TS = i -> TS;
                        j -> DIST = i -> DIST + 1;
                        set_active(j);
                    }
                    else if (!j->is_sink)
                    {
                        s_start = j;
                        t_start = i;
                        cap_middle     = & ( a_for -> r_cap );
                        rev_cap_middle = & ( a_for -> r_rev_cap );
                        break;
                    }
                    else if (j->TS <= i->TS &&
                             j->DIST > i->DIST)
                    {
                        /* heuristic - trying to make the distance from j to the sink shorter */
                        j -> parent = a_for;
                        j -> TS = i -> TS;
                        j -> DIST = i -> DIST + 1;
                    }
                }
            }
        }

        TIME ++;

        if (s_start)
        {
            i -> next = i; /* set active flag */
            current_node = i;

            /* augmentation */
            augment(s_start, t_start, cap_middle, rev_cap_middle);
            /* augmentation end */

            /* adoption */
            while ((np=orphan_first))
            {
                np_next = np -> next;
                np -> next = NULL;

                while ((np=orphan_first))
                {
                    orphan_first = np -> next;
                    i = np -> ptr;
                    nodeptr_block -> Delete(np);
                    if (!orphan_first) orphan_last = NULL;
                    if (i->is_sink) process_sink_orphan(i);
                    else            process_source_orphan(i);
                }

                orphan_first = np_next;
            }
            /* adoption end */
        }
        else current_node = NULL;
    }

    delete nodeptr_block;

    return flow;
}

/***********************************************************************/

Graph::termtype Graph::what_segment(node_id i)
{
    if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;
    return SINK;
}