Added graph cut tools. This is the new min_cut algorithm and find_max_factor_graph_potts()

routine.
This commit is contained in:
Davis King 2012-04-28 16:10:29 -04:00
parent 37bc7b23bd
commit 80e501d8bc
6 changed files with 1929 additions and 0 deletions

13
dlib/graph_cuts.h Normal file
View File

@ -0,0 +1,13 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_GRAPH_CUTs_HEADER_
#define DLIB_GRAPH_CUTs_HEADER_
#include "graph_cuts/min_cut.h"
#include "graph_cuts/general_flow_graph.h"
#include "graph_cuts/find_max_factor_graph_potts.h"
#endif // DLIB_GRAPH_CUTs_HEADER_

View File

@ -0,0 +1,477 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__
#define DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__
#include "find_max_factor_graph_potts_abstract.h"
#include "../matrix.h"
#include "min_cut.h"
namespace dlib
{
// ----------------------------------------------------------------------------------------
namespace impl
{
template <
typename potts_problem,
typename T = void
>
class flows_container
{
/*
This object notionally represents a matrix of flow values. It's
overloaded to represent this matrix efficiently though. In this case
it represents the matrix using a sparse representation.
*/
typedef typename potts_problem::value_type edge_type;
std::vector<std::vector<edge_type> > flows;
public:
void setup(
const potts_problem& p
)
{
flows.resize(p.number_of_nodes());
for (unsigned long i = 0; i < flows.size(); ++i)
{
flows[i].resize(p.number_of_neighbors(i));
}
}
edge_type& operator() (
const long r,
const long c
) { return flows[r][c]; }
const edge_type& operator() (
const long r,
const long c
) const { return flows[r][c]; }
};
// ----------------------------------------------------------------------------------------
template <
typename potts_problem
>
class flows_container<potts_problem,
typename enable_if_c<potts_problem::max_number_of_neighbors!=0>::type>
{
/*
This object notionally represents a matrix of flow values. It's
overloaded to represent this matrix efficiently though. In this case
it represents the matrix using a dense representation.
*/
typedef typename potts_problem::value_type edge_type;
const static unsigned long max_number_of_neighbors = potts_problem::max_number_of_neighbors;
matrix<edge_type,0,max_number_of_neighbors> flows;
public:
void setup(
const potts_problem& p
)
{
flows.set_size(p.number_of_nodes(), max_number_of_neighbors);
}
edge_type& operator() (
const long r,
const long c
) { return flows(r,c); }
const edge_type& operator() (
const long r,
const long c
) const { return flows(r,c); }
};
// ----------------------------------------------------------------------------------------
template <
typename potts_problem
>
class potts_flow_graph
{
public:
typedef typename potts_problem::value_type edge_type;
private:
/*!
This is a utility class used by dlib::min_cut to convert a potts_problem
into the kind of flow graph expected by the min_cut object's main block
of code.
Within this object, we will use the convention that one past
potts_problem::number_of_nodes() is the source node and two past is
the sink node.
!*/
potts_problem& g;
// flows(i,j) == the flow from node id i to it's jth neighbor
flows_container<potts_problem> flows;
// source_flows(i,0) == flow from source to node i,
// source_flows(i,1) == flow from node i to source
matrix<edge_type,0,2> source_flows;
// sink_flows(i,0) == flow from sink to node i,
// sink_flows(i,1) == flow from node i to sink
matrix<edge_type,0,2> sink_flows;
node_label source_label, sink_label;
public:
potts_flow_graph(
potts_problem& g_
) : g(g_)
{
flows.setup(g);
source_flows.set_size(g.number_of_nodes(), 2);
sink_flows.set_size(g.number_of_nodes(), 2);
source_flows = 0;
sink_flows = 0;
source_label = FREE_NODE;
sink_label = FREE_NODE;
// setup flows based on factor potentials
edge_type min_value = std::numeric_limits<edge_type>::max();
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
{
source_flows(i,0) = -g.factor_value(i, true);
sink_flows(i,1) = -g.factor_value(i, false);
source_flows(i,1) = 0;
sink_flows(i,0) = 0;
if (source_flows(i,0) < min_value)
min_value = source_flows(i,0);
if (sink_flows(i,1) < min_value)
min_value = sink_flows(i,1);
for (unsigned long j = 0; j < g.number_of_neighbors(i); ++j)
{
flows(i,j) = g.factor_value_disagreement(i, g.get_neighbor(i,j));
}
}
// The flows can't be negative, so if they are then just adjust them
// so they are positive. Note that doing this doesn't change the optimal
// labeling assignment.
if (min_value < 0)
{
source_flows -= min_value;
sink_flows -= min_value;
}
}
class out_edge_iterator
{
friend class potts_flow_graph;
unsigned long idx; // base node idx
unsigned long cnt; // count over the neighbors of idx
public:
out_edge_iterator(
):idx(0),cnt(0){}
out_edge_iterator(
unsigned long idx_,
unsigned long cnt_
):idx(idx_),cnt(cnt_)
{}
bool operator!= (
const out_edge_iterator& item
) const { return cnt != item.cnt; }
out_edge_iterator& operator++(
)
{
++cnt;
return *this;
}
};
class in_edge_iterator
{
friend class potts_flow_graph;
unsigned long idx; // base node idx
unsigned long cnt; // count over the neighbors of idx
public:
in_edge_iterator(
):idx(0),cnt(0)
{}
in_edge_iterator(
unsigned long idx_,
unsigned long cnt_
):idx(idx_),cnt(cnt_)
{}
bool operator!= (
const in_edge_iterator& item
) const { return cnt != item.cnt; }
in_edge_iterator& operator++(
)
{
++cnt;
return *this;
}
};
unsigned long number_of_nodes (
) const { return g.number_of_nodes() + 2; }
out_edge_iterator out_begin(
const unsigned long& it
) const { return out_edge_iterator(it, 0); }
in_edge_iterator in_begin(
const unsigned long& it
) const { return in_edge_iterator(it, 0); }
out_edge_iterator out_end(
const unsigned long& it
) const
{
if (it >= g.number_of_nodes())
return out_edge_iterator(it, g.number_of_nodes());
else
return out_edge_iterator(it, g.number_of_neighbors(it)+2);
}
in_edge_iterator in_end(
const unsigned long& it
) const
{
if (it >= g.number_of_nodes())
return in_edge_iterator(it, g.number_of_nodes());
else
return in_edge_iterator(it, g.number_of_neighbors(it)+2);
}
template <typename iterator_type>
unsigned long node_id (
const iterator_type& it
) const
{
// if this isn't an iterator over the source or sink nodes
if (it.idx < g.number_of_nodes())
{
const unsigned long num = g.number_of_neighbors(it.idx);
if (it.cnt < num)
return g.get_neighbor(it.idx, it.cnt);
else if (it.cnt == num)
return g.number_of_nodes();
else
return g.number_of_nodes()+1;
}
else
{
return it.cnt;
}
}
edge_type get_flow (
const unsigned long& it1,
const unsigned long& it2
) const
{
if (it1 >= g.number_of_nodes())
{
// if it1 is the source
if (it1 == g.number_of_nodes())
{
return source_flows(it2,0);
}
else // if it1 is the sink
{
return sink_flows(it2,0);
}
}
else if (it2 >= g.number_of_nodes())
{
// if it2 is the source
if (it2 == g.number_of_nodes())
{
return source_flows(it1,1);
}
else // if it2 is the sink
{
return sink_flows(it1,1);
}
}
else
{
return flows(it1, g.get_neighbor_idx(it1, it2));
}
}
edge_type get_flow (
const out_edge_iterator& it
) const
{
if (it.idx < g.number_of_nodes())
{
const unsigned long num = g.number_of_neighbors(it.idx);
if (it.cnt < num)
return flows(it.idx, it.cnt);
else if (it.cnt == num)
return source_flows(it.idx,1);
else
return sink_flows(it.idx,1);
}
else
{
// if it.idx is the source
if (it.idx == g.number_of_nodes())
{
return source_flows(it.cnt,0);
}
else // if it.idx is the sink
{
return sink_flows(it.cnt,0);
}
}
}
edge_type get_flow (
const in_edge_iterator& it
) const
{
return get_flow(node_id(it), it.idx);
}
void adjust_flow (
const unsigned long& it1,
const unsigned long& it2,
const edge_type& value
)
{
if (it1 >= g.number_of_nodes())
{
// if it1 is the source
if (it1 == g.number_of_nodes())
{
source_flows(it2,0) += value;
source_flows(it2,1) -= value;
}
else // if it1 is the sink
{
sink_flows(it2,0) += value;
sink_flows(it2,1) -= value;
}
}
else if (it2 >= g.number_of_nodes())
{
// if it2 is the source
if (it2 == g.number_of_nodes())
{
source_flows(it1,1) += value;
source_flows(it1,0) -= value;
}
else // if it2 is the sink
{
sink_flows(it1,1) += value;
sink_flows(it1,0) -= value;
}
}
else
{
flows(it1, g.get_neighbor_idx(it1, it2)) += value;
flows(it2, g.get_neighbor_idx(it2, it1)) -= value;
}
}
void set_label (
const unsigned long& it,
node_label value
)
{
if (it < g.number_of_nodes())
g.set_label(it, value);
else if (it == g.number_of_nodes())
source_label = value;
else
sink_label = value;
}
node_label get_label (
const unsigned long& it
) const
{
if (it < g.number_of_nodes())
return g.get_label(it);
if (it == g.number_of_nodes())
return source_label;
else
return sink_label;
}
};
}
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
template <
typename potts_model
>
typename potts_model::value_type potts_model_score (
const potts_model& g
)
{
typename potts_model::value_type score = 0;
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
{
const bool label = (g.get_label(i)!=0);
score += g.factor_value(i, label);
}
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
{
for (unsigned long n = 0; n < g.number_of_neighbors(i); ++n)
{
const unsigned long idx2 = g.get_neighbor(i,n);
if (g.get_label(i) != g.get_label(idx2) && i < idx2)
score -= g.factor_value_disagreement(i, idx2);
}
}
return score;
}
// ----------------------------------------------------------------------------------------
template <
typename potts_model
>
void find_max_factor_graph_potts (
potts_model& m
)
{
min_cut mc;
dlib::impl::potts_flow_graph<potts_model> pfg(m);
mc(pfg, m.number_of_nodes(), m.number_of_nodes()+1);
}
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__

View File

@ -0,0 +1,223 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#undef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__
#ifdef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__
#include "../matrix.h"
#include "min_cut_abstract.h"
namespace dlib
{
// ----------------------------------------------------------------------------------------
class potts_problem
{
/*!
WHAT THIS OBJECT REPRESENTS
This object represents a boolean valued factor graph or graphical model
that can be efficiently operated on using graph cuts. In particular, this
object defines the interface a MAP problem on a factor graph must
implement if it is to be solved using the find_max_factor_graph_potts()
routine defined at the bottom of this file.
Note that there is no dlib::potts_problem object. What you are looking
at here is simply the interface definition for a Potts problem. You must
implement your own version of this object for the problem you wish to
solve and then pass it to the find_max_factor_graph_potts() routine.
Note also that a factor graph should not have any nodes which are
neighbors with themselves. Additionally, the graph is undirected. This
mean that if A is a neighbor of B then B must be a neighbor of A for
the MAP problem to be valid.
!*/
public:
unsigned long number_of_nodes (
) const;
/*!
ensures
- returns the number of nodes in the factor graph. Or in other words,
returns the number of variables in the MAP problem/Potts model.
!*/
unsigned long number_of_neighbors (
unsigned long idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns the number of neighbors of node idx.
!*/
// This is an optional variable which specifies a number that is always
// greater than or equal to number_of_neighbors(idx). If you don't know
// the value at compile time then either don't include max_number_of_neighbors
// in your potts_problem object or set it to 0.
const static unsigned long max_number_of_neighbors = 0;
unsigned long get_neighbor (
unsigned long idx,
unsigned long n
) const;
/*!
requires
- idx < number_of_nodes()
- n < number_of_neighbors(idx)
ensures
- returns the node index value of the n-th neighbor of
the node with index value idx.
- The neighbor relationship is reciprocal. That is, if
get_neighbor(A,i)==B then there is a value of j such
that get_neighbor(B,j)==A.
- A node is never its own neighbor. That is, there is
no i such that get_neighbor(idx,i)==idx.
!*/
unsigned long get_neighbor_idx (
unsigned long idx1,
unsigned long idx2
) const;
/*!
requires
- idx1 < number_of_nodes()
- idx2 < number_of_nodes()
ensures
- This function is basically the inverse of get_neighbor().
- returns a number IDX such that:
- get_neighbor(idx1,IDX) == idx2
- IDX < number_of_neighbors(node_idx1)
!*/
void set_label (
const unsigned long& idx,
node_label value
);
/*!
requires
- idx < number_of_nodes()
ensures
- #get_label(idx) == value
!*/
node_label get_label (
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns the current label for the idx-th node. This is a value which is
0 if the node's label is false and is any other value if it is true.
Note that this value is not used by factor_value() or factor_value_disagreement().
It is simply here to provide a mechanism for find_max_factor_graph_potts()
to return its labeled result. Additionally, the reason it returns a
node_label rather than a bool is because doing it this way facilitates
use of a graph cut algorithm for the solution of the MAP problem. For
more of an explanation you should read the paper referenced by the min_cut
object.
!*/
// This typedef should be for a type like int or double. It
// must also be capable of representing signed values.
typedef an_integer_or_real_type value_type;
value_type factor_value (
unsigned long idx,
bool value
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns a value which indicates how "good" it is to assign the idx-node
a label equal to value. The larger the value, the more desirable
the label contained by value.
!*/
value_type factor_value_disagreement (
unsigned long idx1,
unsigned long idx2
) const;
/*!
requires
- idx1 < number_of_nodes()
- idx2 < number_of_nodes()
- idx1 != idx2
- the idx1-th node and idx2-th node are neighbors in the graph. That is,
get_neighbor(idx1,i)==idx2 for some value of i.
ensures
- returns a number >= 0. This is the penalty for giving node idx1 and idx2
different labels. Larger values indicate a larger penalty.
- this function is symmetric. That is, it is true that:
factor_value_disagreement(i,j) == factor_value_disagreement(j,i)
!*/
};
// ----------------------------------------------------------------------------------------
template <
typename potts_problem
>
typename potts_problem::value_type potts_model_score (
const potts_problem& prob
);
/*!
requires
- potts_problem == an object with an interface compatible with the potts_problem
object defined at the top of this file.
- for all valid i and j:
- prob.factor_value_disagreement(i,j) == prob.factor_value_disagreement(j,i)
ensures
- computes the model score for the given potts_problem. To define this
precisely:
- let L(i) == the boolean label of the ith variable in prob. Or in other
words, L(i) == (prob.get_label(i) != 0).
- let F == the sum over valid i of prob.factor_value(i, L(i)).
- Let D == the sum of all values of prob.factor_value_disagreement(i,j)
whenever the following conditions are true about i and j:
- i and j are neighbors in the graph defined by prob, that is,
it is valid to call prob.factor_value_disagreement(i,j).
- L(i) != L(j)
- i < j
(i.e. We want to make sure to only count the edge between i and j once)
- Then this function returns F - D
!*/
// ----------------------------------------------------------------------------------------
template <
typename potts_problem
>
void find_max_factor_graph_potts (
potts_problem& prob
)
/*!
requires
- potts_problem == an object with an interface compatible with the potts_problem
object defined at the top of this file.
- for all valid i and j:
- prob.factor_value_disagreement(i,j) >= 0
- prob.factor_value_disagreement(i,j) == prob.factor_value_disagreement(j,i)
ensures
- This function is a tool for exactly solving the MAP problem in a Potts
model. In particular, this means that this function finds the assignments
to all the labels in prob which maximizes potts_model_score(#prob).
- Note that this routine is a little bit faster if all the values
returned by prob.factor_value() are negative. So if you can arrange for that
to be true without spending any extra CPU cycles then it might be a good idea
to do so.
!*/
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__

View File

@ -0,0 +1,172 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_GENERAL_FLOW_GRaPH_H__
#define DLIB_GENERAL_FLOW_GRaPH_H__
#include "../graph_utils.h"
namespace dlib
{
// ----------------------------------------------------------------------------------------
namespace impl
{
template <
typename directed_graph_type
>
class general_flow_graph
{
/*!
this is a utility class used by dlib::min_cut to convert a directed_graph
into the kind of flow graph expected by the min_cut object's main block
of code.
!*/
directed_graph_type& g;
typedef typename directed_graph_type::node_type node_type;
typedef typename directed_graph_type::type node_label;
public:
general_flow_graph(
directed_graph_type& g_
) : g(g_)
{
}
class out_edge_iterator
{
friend class general_flow_graph;
unsigned long idx; // base node idx
unsigned long cnt; // count over the neighbors of idx
public:
out_edge_iterator(
):idx(0),cnt(0) {}
out_edge_iterator(
unsigned long idx_,
unsigned long cnt_
):idx(idx_),cnt(cnt_)
{}
bool operator!= (
const out_edge_iterator& item
) const { return cnt != item.cnt; }
out_edge_iterator& operator++(
)
{
++cnt;
return *this;
}
};
class in_edge_iterator
{
friend class general_flow_graph;
unsigned long idx; // base node idx
unsigned long cnt; // count over the neighbors of idx
public:
in_edge_iterator(
):idx(0),cnt(0) {}
in_edge_iterator(
unsigned long idx_,
unsigned long cnt_
):idx(idx_),cnt(cnt_)
{}
bool operator!= (
const in_edge_iterator& item
) const { return cnt != item.cnt; }
in_edge_iterator& operator++(
)
{
++cnt;
return *this;
}
};
unsigned long number_of_nodes (
) const { return g.number_of_nodes(); }
out_edge_iterator out_begin(
const unsigned long& it
) const { return out_edge_iterator(it, 0); }
in_edge_iterator in_begin(
const unsigned long& it
) const { return in_edge_iterator(it, 0); }
out_edge_iterator out_end(
const unsigned long& it
) const { return out_edge_iterator(it, g.node(it).number_of_children()); }
in_edge_iterator in_end(
const unsigned long& it
) const { return in_edge_iterator(it, g.node(it).number_of_parents()); }
unsigned long node_id (
const out_edge_iterator& it
) const { return g.node(it.idx).child(it.cnt).index(); }
unsigned long node_id (
const in_edge_iterator& it
) const { return g.node(it.idx).parent(it.cnt).index(); }
typedef typename directed_graph_type::edge_type edge_type;
edge_type get_flow (const unsigned long& it1, const unsigned long& it2) const
{
return edge(g, it1, it2);
}
edge_type get_flow (const out_edge_iterator& it) const
{
return g.node(it.idx).child_edge(it.cnt);
}
edge_type get_flow (const in_edge_iterator& it) const
{
return g.node(it.idx).parent_edge(it.cnt);
}
void adjust_flow (
const unsigned long& it1,
const unsigned long& it2,
const edge_type& value
)
{
edge(g, it1, it2) += value;
edge(g, it2, it1) -= value;
}
void set_label (
const unsigned long& it,
node_label value
)
{
g.node(it).data = value;
}
node_label get_label (
const unsigned long& it
) const
{
return g.node(it).data;
}
};
}
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_GENERAL_FLOW_GRaPH_H__

572
dlib/graph_cuts/min_cut.h Normal file
View File

@ -0,0 +1,572 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_MIN_CuT_H__
#define DLIB_MIN_CuT_H__
#include "min_cut_abstract.h"
#include "../matrix.h"
#include "general_flow_graph.h"
#include "../is_kind.h"
#include <iostream>
#include <fstream>
#include <deque>
// ----------------------------------------------------------------------------------------
namespace dlib
{
typedef unsigned char node_label;
// ----------------------------------------------------------------------------------------
const node_label SOURCE_CUT = 0;
const node_label SINK_CUT = 1;
const node_label FREE_NODE = 2;
// ----------------------------------------------------------------------------------------
template <typename flow_graph>
typename disable_if<is_directed_graph<flow_graph>,typename flow_graph::edge_type>::type
graph_cut_score (
const flow_graph& g
)
{
typedef typename flow_graph::edge_type edge_weight_type;
edge_weight_type score = 0;
typedef typename flow_graph::out_edge_iterator out_edge_iterator;
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
{
if (g.get_label(i) != SOURCE_CUT)
continue;
for (out_edge_iterator n = g.out_begin(i); n != g.out_end(i); ++n)
{
if (g.get_label(g.node_id(n)) != SOURCE_CUT)
{
score += g.get_flow(n);
}
}
}
return score;
}
template <typename directed_graph>
typename enable_if<is_directed_graph<directed_graph>,typename directed_graph::edge_type>::type
graph_cut_score (
const directed_graph& g
)
{
return graph_cut_score(dlib::impl::general_flow_graph<const directed_graph>(g));
}
// ----------------------------------------------------------------------------------------
class min_cut
{
public:
min_cut()
{
}
min_cut( const min_cut& )
{
// Intentionally left empty since all the member variables
// don't logically contribute to the state of this object.
// This copy constructor is here to explicitly void the overhead
// of copying these transient variables.
}
template <
typename directed_graph
>
typename enable_if<is_directed_graph<directed_graph> >::type operator() (
directed_graph& g,
const unsigned long source_node,
const unsigned long sink_node
) const
{
DLIB_ASSERT(graph_contains_length_one_cycle(g) == false,
"\t void min_cut::operator()"
<< "\n\t Invalid arguments were given to this function."
);
DLIB_ASSERT(graph_has_symmetric_edges(g) == true,
"\t void min_cut::operator()"
<< "\n\t Invalid arguments were given to this function."
);
dlib::impl::general_flow_graph<directed_graph> temp(g);
(*this)(temp, source_node, sink_node);
}
template <
typename flow_graph
>
typename disable_if<is_directed_graph<flow_graph> >::type operator() (
flow_graph& g,
const unsigned long source_node,
const unsigned long sink_node
) const
{
#ifdef ENABLE_ASSERTS
DLIB_ASSERT(source_node != sink_node &&
source_node < g.number_of_nodes() &&
sink_node < g.number_of_nodes(),
"\t void min_cut::operator()"
<< "\n\t Invalid arguments were given to this function."
<< "\n\t g.number_of_nodes(): " << g.number_of_nodes()
<< "\n\t source_node: " << source_node
<< "\n\t sink_node: " << sink_node
<< "\n\t this: " << this
);
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
{
typename flow_graph::out_edge_iterator j, end = g.out_end(i);
for (j = g.out_begin(i); j != end; ++j)
{
const unsigned long jj = g.node_id(j);
DLIB_ASSERT(g.get_flow(i,jj) >= 0,
"\t void min_cut::operator()"
<< "\n\t i: "<< i
<< "\n\t jj: "<< jj
<< "\n\t g.get_flow(i,jj): "<< g.get_flow(i,jj)
<< "\n\t this: "<< this
);
}
}
#endif
parent.clear();
active.clear();
orphans.clear();
typedef typename flow_graph::edge_type edge_type;
typedef typename flow_graph::out_edge_iterator out_edge_iterator;
typedef typename flow_graph::in_edge_iterator in_edge_iterator;
// initialize labels
for (unsigned long i = 0; i < g.number_of_nodes(); ++i)
g.set_label(i, FREE_NODE);
g.set_label(source_node, SOURCE_CUT);
g.set_label(sink_node, SINK_CUT);
// used to indicate "no parent"
const unsigned long nil = g.number_of_nodes();
parent.assign(g.number_of_nodes(), nil);
time = 1;
dist.assign(g.number_of_nodes(), 0);
ts.assign(g.number_of_nodes(), time);
active.push_back(source_node);
active.push_back(sink_node);
in_edge_iterator in_begin = g.in_begin(active[0]);
out_edge_iterator out_begin = g.out_begin(active[0]);
unsigned long source_side, sink_side;
while (grow(g,source_side,sink_side, in_begin, out_begin))
{
++time;
ts[source_node] = time;
ts[sink_node] = time;
augment(g, source_node, sink_node, source_side, sink_side);
adopt(g, source_node, sink_node);
}
}
private:
unsigned long distance_to_origin (
const unsigned long nil,
unsigned long p,
unsigned long
) const
{
unsigned long start = p;
unsigned long count = 0;
while (p != nil)
{
if (ts[p] == time)
{
count += dist[p];
unsigned long count_down = count;
// adjust the dist and ts for the nodes on this path.
while (start != p)
{
ts[start] = time;
dist[start] = count_down;
--count_down;
start = parent[start];
}
return count;
}
p = parent[p];
++count;
}
return std::numeric_limits<unsigned long>::max();
}
template <typename flow_graph>
void adopt (
flow_graph& g,
const unsigned long source,
const unsigned long sink
) const
{
typedef typename flow_graph::edge_type edge_type;
typedef typename flow_graph::out_edge_iterator out_edge_iterator;
typedef typename flow_graph::in_edge_iterator in_edge_iterator;
// used to indicate "no parent"
const unsigned long nil = g.number_of_nodes();
while (orphans.size() > 0)
{
const unsigned long p = orphans.front();
orphans.pop_front();
const unsigned char label_p = g.get_label(p);
// Try to find a valid parent for p.
if (label_p == SOURCE_CUT)
{
const in_edge_iterator begin(g.in_begin(p));
const in_edge_iterator end(g.in_end(p));
unsigned long best_dist = std::numeric_limits<unsigned long>::max();
unsigned long best_node = 0;
for(in_edge_iterator q = begin; q != end; ++q)
{
const unsigned long id = g.node_id(q);
if (g.get_label(id) != label_p || g.get_flow(q) <= 0 )
continue;
unsigned long temp = distance_to_origin(nil, id,source);
if (temp < best_dist)
{
best_dist = temp;
best_node = id;
}
}
if (best_dist != std::numeric_limits<unsigned long>::max())
{
parent[p] = best_node;
dist[p] = dist[best_node] + 1;
ts[p] = time;
}
// if we didn't find a parent for p
if (parent[p] == nil)
{
for(in_edge_iterator q = begin; q != end; ++q)
{
const unsigned long id = g.node_id(q);
if (g.get_label(id) != SOURCE_CUT)
continue;
if (g.get_flow(q) > 0)
active.push_back(id);
if (parent[id] == p)
{
parent[id] = nil;
orphans.push_back(id);
}
}
g.set_label(p, FREE_NODE);
}
}
else
{
unsigned long best_node = 0;
unsigned long best_dist = std::numeric_limits<unsigned long>::max();
const out_edge_iterator begin(g.out_begin(p));
const out_edge_iterator end(g.out_end(p));
for(out_edge_iterator q = begin; q != end; ++q)
{
const unsigned long id = g.node_id(q);
if (g.get_label(id) != label_p || g.get_flow(q) <= 0)
continue;
unsigned long temp = distance_to_origin(nil, id,sink);
if (temp < best_dist)
{
best_dist = temp;
best_node = id;
}
}
if (best_dist != std::numeric_limits<unsigned long>::max())
{
parent[p] = best_node;
dist[p] = dist[best_node] + 1;
ts[p] = time;
}
// if we didn't find a parent for p
if (parent[p] == nil)
{
for(out_edge_iterator q = begin; q != end; ++q)
{
const unsigned long id = g.node_id(q);
if (g.get_label(id) != SINK_CUT)
continue;
if (g.get_flow(q) > 0)
active.push_back(id);
if (parent[id] == p)
{
parent[id] = nil;
orphans.push_back(id);
}
}
g.set_label(p, FREE_NODE);
}
}
}
}
template <typename flow_graph>
void augment (
flow_graph& g,
const unsigned long& source,
const unsigned long& sink,
const unsigned long& source_side,
const unsigned long& sink_side
) const
{
typedef typename flow_graph::edge_type edge_type;
typedef typename flow_graph::out_edge_iterator out_edge_iterator;
typedef typename flow_graph::in_edge_iterator in_edge_iterator;
// used to indicate "no parent"
const unsigned long nil = g.number_of_nodes();
unsigned long s = source_side;
unsigned long t = sink_side;
edge_type min_cap = g.get_flow(s,t);
// find the bottleneck capacity on the current path.
// check from source_side back to the source for the min capacity link.
t = s;
while (t != source)
{
s = parent[t];
const edge_type temp = g.get_flow(s, t);
if (temp < min_cap)
{
min_cap = temp;
}
t = s;
}
// check from sink_side back to the sink for the min capacity link
s = sink_side;
while (s != sink)
{
t = parent[s];
const edge_type temp = g.get_flow(s, t);
if (temp < min_cap)
{
min_cap = temp;
}
s = t;
}
// now push the max possible amount of flow though the path
s = source_side;
t = sink_side;
g.adjust_flow(t,s, min_cap);
// trace back towards the source
t = s;
while (t != source)
{
s = parent[t];
g.adjust_flow(t,s, min_cap);
if (g.get_flow(s,t) <= 0)
{
parent[t] = nil;
orphans.push_back(t);
}
t = s;
}
// trace back towards the sink
s = sink_side;
while (s != sink)
{
t = parent[s];
g.adjust_flow(t,s, min_cap);
if (g.get_flow(s,t) <= 0)
{
parent[s] = nil;
orphans.push_back(s);
}
s = t;
}
}
template <typename flow_graph>
bool grow (
flow_graph& g,
unsigned long& source_side,
unsigned long& sink_side,
typename flow_graph::in_edge_iterator& in_begin,
typename flow_graph::out_edge_iterator& out_begin
) const
/*!
ensures
- if (an augmenting path was found) then
- returns true
- (#source_side, #sink_side) == the point where the two trees meet.
#source_side is part of the source tree and #sink_side is part of
the sink tree.
- else
- returns false
!*/
{
typedef typename flow_graph::edge_type edge_type;
typedef typename flow_graph::out_edge_iterator out_edge_iterator;
typedef typename flow_graph::in_edge_iterator in_edge_iterator;
while (active.size() != 0)
{
// pick an active node
const unsigned long A = active[0];
const unsigned char label_A = g.get_label(A);
// process its neighbors
if (label_A == SOURCE_CUT)
{
const out_edge_iterator out_end = g.out_end(A);
for(out_edge_iterator& i = out_begin; i != out_end; ++i)
{
if (g.get_flow(i) > 0)
{
const unsigned long id = g.node_id(i);
const unsigned char label_i = g.get_label(id);
if (label_i == FREE_NODE)
{
active.push_back(id);
g.set_label(id,SOURCE_CUT);
parent[id] = A;
ts[id] = ts[A];
dist[id] = dist[A] + 1;
}
else if (label_A != label_i)
{
source_side = A;
sink_side = id;
return true;
}
else if (is_closer(A, id))
{
parent[id] = A;
ts[id] = ts[A];
dist[id] = dist[A] + 1;
}
}
}
}
else if (label_A == SINK_CUT)
{
const in_edge_iterator in_end = g.in_end(A);
for(in_edge_iterator& i = in_begin; i != in_end; ++i)
{
if (g.get_flow(i) > 0)
{
const unsigned long id = g.node_id(i);
const unsigned char label_i = g.get_label(id);
if (label_i == FREE_NODE)
{
active.push_back(id);
g.set_label(id,SINK_CUT);
parent[id] = A;
ts[id] = ts[A];
dist[id] = dist[A] + 1;
}
else if (label_A != label_i)
{
sink_side = A;
source_side = id;
return true;
}
else if (is_closer(A, id))
{
parent[id] = A;
ts[id] = ts[A];
dist[id] = dist[A] + 1;
}
}
}
}
active.pop_front();
if (active.size() != 0)
{
in_begin = g.in_begin(active[0]);
out_begin = g.out_begin(active[0]);
}
}
return false;
}
inline bool is_closer (
unsigned long p,
unsigned long q
) const
{
// return true if p is closer to a terminal than q
return ts[q] <= ts[p] && dist[q] > dist[p];
}
mutable std::vector<uint32> dist;
mutable std::vector<uint32> ts;
mutable uint32 time;
mutable std::vector<unsigned long> parent;
mutable std::deque<unsigned long> active;
mutable std::deque<unsigned long> orphans;
};
// ----------------------------------------------------------------------------------------
}
// ----------------------------------------------------------------------------------------
#endif // DLIB_MIN_CuT_H__

View File

@ -0,0 +1,472 @@
// Copyright (C) 2012 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#undef DLIB_MIN_CuT_ABSTRACT_H__
#ifdef DLIB_MIN_CuT_ABSTRACT_H__
#include "../graph_utils.h"
// ----------------------------------------------------------------------------------------
namespace dlib
{
typedef unsigned char node_label;
// ----------------------------------------------------------------------------------------
template <
typename directed_graph_type
>
class flow_graph
{
/*!
WHAT THIS OBJECT REPRESENTS
This object represents a flow capacity graph for use with the
min_cut algorithm defined below. In particular, this object
is a kind of directed graph where the edge weights specify the
flow capacities.
Note that there is no dlib::flow_graph object. What you are
looking at here is simply the interface definition for a graph
which can be used with the min_cut algorithm. You must implement
your own version of this object for the graph you wish to work with
and then pass it to the min_cut::operator() routine.
It's also worth pointing out that this graph has symmetric edge
connections. That is, if there is an edge from node A to node B
then there must also be an edge from node B to node A.
!*/
public:
class out_edge_iterator
{
/*!
WHAT THIS OBJECT REPRESENTS
This is a simple forward iterator for iterating over the neighbors
of a node in the graph. It also represents the fact that the neighbors
are on the end of an outgoing edge. That is, the edge represents
the amount of flow which can flow towards the neighbor.
!*/
public:
out_edge_iterator(
);
/*!
ensures
- constructs an iterator in an undefined state. It can't
be used until assigned with a valid iterator.
!*/
out_edge_iterator(
const out_edge_iterator& item
);
/*!
ensures
- #*this is a copy of item
!*/
out_edge_iterator& operator=(
const out_edge_iterator& item
);
/*!
ensures
- #*this is a copy of item
- returns #*this
!*/
bool operator!= (
const out_edge_iterator& item
) const;
/*!
requires
- *this and item are iterators over the neighbors for the
same node.
ensures
- returns false if *this and item both reference the same
node in the graph and true otherwise.
!*/
out_edge_iterator& operator++(
);
/*!
ensures
- advances *this to the next neighbor node.
- returns a reference to the updated *this
(i.e. this is the ++object form of the increment operator)
!*/
};
class in_edge_iterator
{
/*!
WHAT THIS OBJECT REPRESENTS
This is a simple forward iterator for iterating over the neighbors
of a node in the graph. It also represents the fact that the neighbors
are on the end of an incoming edge. That is, the edge represents
the amount of flow which can flow out of the neighbor node.
!*/
public:
in_edge_iterator(
);
/*!
ensures
- constructs an iterator in an undefined state. It can't
be used until assigned with a valid iterator.
!*/
in_edge_iterator(
const in_edge_iterator& item
);
/*!
ensures
- #*this is a copy of item
!*/
in_edge_iterator& operator=(
const in_edge_iterator& item
);
/*!
ensures
- #*this is a copy of item
- returns #*this
!*/
bool operator!= (
const in_edge_iterator& item
) const;
/*!
requires
- *this and item are iterators over the neighbors for the
same node.
ensures
- returns false if *this and item both reference the same
node in the graph and true otherwise.
!*/
in_edge_iterator& operator++(
);
/*!
ensures
- advances *this to the next neighbor node.
- returns a reference to the updated *this
(i.e. this is the ++object form of the increment operator)
!*/
};
unsigned long number_of_nodes (
) const;
/*!
ensures
- returns the number of nodes in the graph.
!*/
out_edge_iterator out_begin(
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns an iterator pointing to the first neighboring node of
the idx-th node. If no such node exists then returns out_end(idx).
- The returned iterator also represents the directed edge going from
node idx to the neighbor.
!*/
in_edge_iterator in_begin(
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns an iterator pointing to the first neighboring node of
the idx-th node. If no such node exists then returns in_end(idx).
- The returned iterator also represents the directed edge going from
the neighbor to node idx.
!*/
out_edge_iterator out_end(
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns an iterator to one past the last neighboring node of
the idx-th node.
!*/
in_edge_iterator in_end(
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns an iterator to one past the last neighboring node of
the idx-th node.
!*/
unsigned long node_id (
const out_edge_iterator& it
) const;
/*!
requires
- it == a valid iterator (i.e. it must be in the range [out_begin(idx), out_end(idx))
for some valid idx)
ensures
- returns a number IDX such that:
- 0 <= IDX < number_of_nodes()
- IDX == The index which uniquely identifies the node pointed to by the
iterator it. This number can be used with any member function in this
object which expect a node index. (e.g. get_label(IDX) == the label for the
node pointed to by it)
!*/
unsigned long node_id (
const in_edge_iterator& it
) const;
/*!
requires
- it == a valid iterator (i.e. it must be in the range [in_begin(idx), in_end(idx))
for some valid idx)
ensures
- returns a number IDX such that:
- 0 <= IDX < number_of_nodes()
- IDX == The index which uniquely identifies the node pointed to by the
iterator it. This number can be used with any member function in this
object which expect a node index. (e.g. get_label(IDX) == the label for the
node pointed to by it)
!*/
// This typedef should be for a type like int or double. It
// must also be capable of representing signed values.
typedef an_integer_or_real_type edge_type;
edge_type get_flow (
const unsigned long& idx1,
const unsigned long& idx2
) const;
/*!
requires
- idx1 < number_of_nodes()
- idx2 < number_of_nodes()
- idx1 and idx2 are neighbors in the graph
ensures
- returns the residual flow capacity from the idx1-th
node to the idx2-th node.
!*/
edge_type get_flow (
const out_edge_iterator& it
) const;
/*!
requires
- it == a valid iterator (i.e. it must be in the range [out_begin(idx), out_end(idx))
for some valid idx)
ensures
- let IDX = node_id(it)
- it represents the directed edge from a node, call it H, to the node IDX. Therefore,
this function returns get_flow(H,IDX)
!*/
edge_type get_flow (
const in_edge_iterator& it
) const;
/*!
requires
- it == a valid iterator (i.e. it must be in the range [in_begin(idx), in_end(idx))
for some valid idx)
ensures
- let IDX = node_id(it)
- it represents the directed edge from node IDX to another node, call it H. Therefore,
this function returns get_flow(IDX,H)
!*/
void adjust_flow (
const unsigned long& idx1,
const unsigned long& idx2,
const edge_type& value
);
/*!
requires
- idx1 < number_of_nodes()
- idx2 < number_of_nodes()
- idx1 and idx2 are neighbors in the graph
ensures
- #get_flow(idx1,idx2) == get_flow(idx1,idx2) + value
- #get_flow(idx2,idx1) == get_flow(idx2,idx1) - value
!*/
void set_label (
const unsigned long& idx,
node_label value
);
/*!
requires
- idx < number_of_nodes()
ensures
- #get_label(idx) == value
!*/
node_label get_label (
const unsigned long& idx
) const;
/*!
requires
- idx < number_of_nodes()
ensures
- returns the label for the idx-th node in the graph.
!*/
};
// ----------------------------------------------------------------------------------------
const node_label SOURCE_CUT = 0;
const node_label SINK_CUT = 1;
const node_label FREE_NODE = 2;
// ----------------------------------------------------------------------------------------
template <
typename flow_graph
>
typename flow_graph::edge_type graph_cut_score (
const flow_graph& g
);
/*!
requires
- flow_graph == an object with an interface compatible with the flow_graph
object defined at the top of this file, or, an implementation of
dlib/directed_graph/directed_graph_kernel_abstract.h.
ensures
- returns the sum of the outgoing edges from nodes with a label of SOURCE_CUT
to nodes with a label != SOURCE_CUT. Note that for a directed_graph object,
the labels are stored in the node's data field.
!*/
// ----------------------------------------------------------------------------------------
class min_cut
{
/*!
WHAT THIS OBJECT REPRESENTS
This is a function object which can be used to find the min cut
on a graph.
The implementation is based on the method described in the following
paper:
An Experimental Comparison of Min-Cut/Max-Flow Algorithms for
Energy Minimization in Vision, by Yuri Boykov and Vladimir Kolmogorov,
in PAMI 2004.
!*/
public:
min_cut(
);
/*!
ensures
- this object is properly initialized
!*/
template <
typename flow_graph
>
void operator() (
flow_graph& g,
const unsigned long source_node,
const unsigned long sink_node
) const;
/*!
requires
- flow_graph == an object with an interface compatible with the flow_graph
object defined at the top of this file.
- source_node != sink_node
- source_node < g.number_of_nodes()
- sink_node < g.number_of_nodes()
- for all valid i and j:
- g.get_flow(i,j) >= 0
(i.e. all the flow capacities/edge weights are positive)
- g does not contain any self loops. That is, no nodes are neighbors with
themselves.
ensures
- Finds the minimum cut on the given graph. That is, this function finds
a labeling of nodes in g such that graph_cut_score(g) would be minimized. Note
that the flow values in #g are modified by this algorithm so if you want
to obtain the min cut score you must call min_cut::operator(), then copy
the flow values back into #g, and then call graph_cut_score(#g). But in most
cases you don't care about the value of the min cut score, rather, you
just want the labels in #g.
- #g.get_label(source_node) == SOURCE_CUT
- #g.get_label(sink_node) == SINK_CUT
- for all valid i:
- #g.get_label(i) == SOURCE_CUT, SINK_CUT, or FREE_NODE
- if (#g.get_label(i) == SOURCE_CUT) then
- The minimum cut of g places node i into the source side of the cut.
- if (#g.get_label(i) == SINK_CUT) then
- The minimum cut of g places node i into the sink side of the cut.
- if (#g.get_label(i) == FREE_NODE) then
- Node i can be labeled SOURCE_CUT or SINK_CUT. Both labelings
result in the same cut score.
- When interpreting g as a graph of flow capacities from the source_node
to the sink_node we can say that the min cut problem is equivalent to
the max flow problem. This equivalent problem is to find out how to push
as much "flow" from the source node to the sink node as possible.
Upon termination, #g will contain the final flow residuals in addition to
the graph cut labels. That is, for all valid i and j:
- #g.get_flow(i,j) == the residual flow capacity left after the max
possible amount of flow is passing from the source node to the sink
node. For example, this means that #g.get_flow(i,j) == 0 whenever
node i and j are in different cuts.
- #g.get_flow(i,j) >= 0
!*/
template <
typename directed_graph
>
void operator() (
directed_graph& g,
const unsigned long source_node,
const unsigned long sink_node
) const;
/*!
requires
- directed_graph == an implementation of dlib/directed_graph/directed_graph_kernel_abstract.h
- directed_graph::type == node_label
- directed_graph::edge_type == and integer or double type
- source_node != sink_node
- source_node < g.number_of_nodes()
- sink_node < g.number_of_nodes()
- for all valid i and j:
- edge(g,i,j) >= 0
(i.e. all the flow capacities/edge weights are positive)
- graph_contains_length_one_cycle(g) == false
- graph_has_symmetric_edges(g) == true
ensures
- This routine simply converts g into a flow graph and calls the version
of operator() defined above. Note that the conversion is done in O(1)
time, it's just an interface adaptor.
- edge weights in g correspond to network flows while the .data field of
each node in g corresponds to the graph node labels.
- upon termination, the flows and labels in g will have been modified
as described in the above operator() routine.
!*/
};
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_MIN_CuT_ABSTRACT_H__