From bec858191621e724b9d05acaff35eb489293c5b5 Mon Sep 17 00:00:00 2001 From: Brecht Van Lommel Date: Thu, 19 Dec 2024 01:09:05 +0100 Subject: [PATCH] Refactor: Remove Boost from Quadriflow Extract Boykov-Kolmogorov max flow algorithm from boost::graph. As part of work towards removing Boost as a Blender dependency. Together with #133347 this removes the last direct Blender dependency on Boost. Pull Request: https://projects.blender.org/blender/blender/pulls/132142 --- CMakeLists.txt | 1 - extern/quadriflow/CMakeLists.txt | 3 +- extern/quadriflow/patches/blender.patch | 166 ++++- .../patches/boykov_kolmogorov_max_flow.hpp | 691 ++++++++++++++++++ extern/quadriflow/src/config.hpp | 5 - extern/quadriflow/src/flow.hpp | 101 +-- extern/quadriflow/src/post-solver.cpp | 2 + intern/quadriflow/CMakeLists.txt | 2 - 8 files changed, 879 insertions(+), 92 deletions(-) create mode 100644 extern/quadriflow/patches/boykov_kolmogorov_max_flow.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a9061d06ec..6f22f8b4a22 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1278,7 +1278,6 @@ set_and_warn_dependency(WITH_PYTHON WITH_MOD_FLUID OFF) # otherwise if the user disabled set_and_warn_dependency(WITH_BOOST WITH_OPENVDB OFF) -set_and_warn_dependency(WITH_BOOST WITH_QUADRIFLOW OFF) set_and_warn_dependency(WITH_BOOST WITH_USD OFF) if(WITH_CYCLES) set_and_warn_dependency(WITH_BOOST WITH_CYCLES_OSL OFF) diff --git a/extern/quadriflow/CMakeLists.txt b/extern/quadriflow/CMakeLists.txt index 8c0938ca730..f42a68ab246 100644 --- a/extern/quadriflow/CMakeLists.txt +++ b/extern/quadriflow/CMakeLists.txt @@ -51,11 +51,11 @@ set(INC ) set(INC_SYS - ${BOOST_INCLUDE_DIR} ${EIGEN3_INCLUDE_DIRS} ) set(SRC + patches/boykov_kolmogorov_max_flow.hpp src/adjacent-matrix.cpp src/adjacent-matrix.hpp src/compare-key.hpp @@ -91,7 +91,6 @@ set(SRC ) set(LIB - ${BOOST_LIBRARIES} ) blender_add_lib(extern_quadriflow "${SRC}" "${INC}" "${INC_SYS}" "${LIB}") diff --git a/extern/quadriflow/patches/blender.patch b/extern/quadriflow/patches/blender.patch index ade60d28d02..cc7a21adaa8 100644 --- a/extern/quadriflow/patches/blender.patch +++ b/extern/quadriflow/patches/blender.patch @@ -189,22 +189,6 @@ index 8de74ede8a9..f9861f39169 100644 num = length - shift; while (num--) { -diff --git a/extern/quadriflow/src/config.hpp b/extern/quadriflow/src/config.hpp -index 842b885a209..bf597ad0f39 100644 ---- a/extern/quadriflow/src/config.hpp -+++ b/extern/quadriflow/src/config.hpp -@@ -1,6 +1,11 @@ - #ifndef CONFIG_H_ - #define CONFIG_H_ - -+/* Workaround a bug in boost 1.68, until we upgrade to a newer version. */ -+#if defined(__clang__) && defined(WIN32) -+ #include -+ using namespace boost; -+#endif - // Move settings to cmake to make CMake happy :) - - // #define WITH_SCALE diff --git a/extern/quadriflow/src/hierarchy.cpp b/extern/quadriflow/src/hierarchy.cpp index c333256a139..70a9628320f 100644 --- a/extern/quadriflow/src/hierarchy.cpp @@ -255,3 +239,153 @@ index aa27066e6e4..5b9d717db71 100644 std::size_t operator()(const obj_vertex &v) const { size_t hash = std::hash()(v.p); hash = hash * 37 + std::hash()(v.uv); +diff --git a/extern/quadriflow/src/flow.hpp b/extern/quadriflow/src/flow.hpp +index ab4a01c..a77f7ae 100644 +--- a/extern/quadriflow/src/flow.hpp ++++ b/extern/quadriflow/src/flow.hpp +@@ -7,17 +7,12 @@ + #include + + #include "config.hpp" +- +-#include +-#include +-#include +-#include ++#include "../patches/boykov_kolmogorov_max_flow.hpp" + + #include + #include + #include + +-using namespace boost; + using namespace Eigen; + + namespace qflow { +@@ -34,78 +29,52 @@ class MaxFlowHelper { + + class BoykovMaxFlowHelper : public MaxFlowHelper { + public: +- typedef int EdgeWeightType; +- typedef adjacency_list_traits Traits; +- // clang-format off +- typedef adjacency_list < vecS, vecS, directedS, +- property < vertex_name_t, std::string, +- property < vertex_index_t, long, +- property < vertex_color_t, boost::default_color_type, +- property < vertex_distance_t, long, +- property < vertex_predecessor_t, Traits::edge_descriptor > > > > >, +- +- property < edge_capacity_t, EdgeWeightType, +- property < edge_residual_capacity_t, EdgeWeightType, +- property < edge_reverse_t, Traits::edge_descriptor > > > > Graph; +- // clang-format on +- +- public: +- BoykovMaxFlowHelper() { rev = get(edge_reverse, g); } +- void resize(int n, int m) { +- vertex_descriptors.resize(n); +- for (int i = 0; i < n; ++i) vertex_descriptors[i] = add_vertex(g); +- } +- int compute() { +- EdgeWeightType flow = +- boykov_kolmogorov_max_flow(g, vertex_descriptors.front(), vertex_descriptors.back()); +- return flow; ++ BoykovMaxFlowHelper() = default; ++ void resize(int n, int m) override { ++ num_verts = n; ++ num_edges = 0; ++ flow.resize(num_verts, m * 2); + } +- void addDirectEdge(Traits::vertex_descriptor& v1, Traits::vertex_descriptor& v2, +- property_map::type& rev, const int capacity, +- const int inv_capacity, Graph& g, Traits::edge_descriptor& e1, +- Traits::edge_descriptor& e2) { +- e1 = add_edge(v1, v2, g).first; +- e2 = add_edge(v2, v1, g).first; +- put(edge_capacity, g, e1, capacity); +- put(edge_capacity, g, e2, inv_capacity); +- +- rev[e1] = e2; +- rev[e2] = e1; ++ int compute() override { ++ return flow.max_flow(0, num_verts - 1); + } +- void addEdge(int x, int y, int c, int rc, int v, int cost = 1) { +- Traits::edge_descriptor e1, e2; +- addDirectEdge(vertex_descriptors[x], vertex_descriptors[y], rev, c, rc, g, e1, e2); ++ void addEdge(int x, int y, int c, int rc, int v, int cost = 1) override { ++ const int e1 = num_edges++; ++ const int e2 = num_edges++; ++ flow.set_edge(e1, e2, x, y, c); ++ flow.set_edge(e2, e1, y, x, rc); + if (v != -1) { +- edge_to_variables[e1] = std::make_pair(v, -1); +- edge_to_variables[e2] = std::make_pair(v, 1); ++ edge_to_variables.emplace_back(v, -1); ++ edge_to_variables.emplace_back(v, 1); ++ } ++ else { ++ edge_to_variables.emplace_back(-1, -1); ++ edge_to_variables.emplace_back(-1, -1); + } + } +- void applyTo(std::vector& edge_diff) { +- property_map::type capacity = get(edge_capacity, g); +- property_map::type residual_capacity = +- get(edge_residual_capacity, g); +- +- graph_traits::vertex_iterator u_iter, u_end; +- graph_traits::out_edge_iterator ei, e_end; +- for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) +- for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) +- if (capacity[*ei] > 0) { +- int flow = (capacity[*ei] - residual_capacity[*ei]); ++ void applyTo(std::vector& edge_diff) override { ++ for (int vert = 0; vert < num_verts; vert++) { ++ for (int edge : flow.vertex_out_edges(vert)) { ++ const int capacity = flow.edge_capacity(edge); ++ const int residual_capacity = flow.edge_residual_capacity(edge); ++ if (capacity > 0) { ++ int flow = (capacity - residual_capacity); + if (flow > 0) { +- auto it = edge_to_variables.find(*ei); +- if (it != edge_to_variables.end()) { +- edge_diff[it->second.first / 2][it->second.first % 2] += +- it->second.second * flow; ++ std::pair e2v = edge_to_variables[edge]; ++ if (e2v.first != -1) { ++ edge_diff[e2v.first / 2][e2v.first % 2] += e2v.second * flow; + } + } + } ++ } ++ } + } + + private: +- Graph g; +- property_map::type rev; +- std::vector vertex_descriptors; +- std::map> edge_to_variables; ++ BoykovKolmogorovMaxFlow flow; ++ std::vector> edge_to_variables; ++ int num_verts = 0; ++ int num_edges = 0; + }; + + class NetworkSimplexFlowHelper : public MaxFlowHelper { +diff --git a/extern/quadriflow/src/post-solver.cpp b/extern/quadriflow/src/post-solver.cpp +index 6027ddd..ccefd15 100644 +--- a/extern/quadriflow/src/post-solver.cpp ++++ b/extern/quadriflow/src/post-solver.cpp +@@ -5,7 +5,9 @@ + // Created by Jingwei on 2/5/18. + // + #include ++#ifdef POST_SOLVER + #include ++#endif + #include + #include + #include diff --git a/extern/quadriflow/patches/boykov_kolmogorov_max_flow.hpp b/extern/quadriflow/patches/boykov_kolmogorov_max_flow.hpp new file mode 100644 index 00000000000..96df563f0c4 --- /dev/null +++ b/extern/quadriflow/patches/boykov_kolmogorov_max_flow.hpp @@ -0,0 +1,691 @@ +// SPDX-FileCopyrightText: 2006 Stephan Diederich +// SPDX-FileCopyrightText: 2024 Blender Authors +// SPDX-License-Identifier: MIT +// +// Adapted from boost::graph + +#include +#include +#include +#include +#include +#include + +namespace qflow { + +class BoykovKolmogorovMaxFlow { + // Types + enum class Color { white, black, gray }; + + struct Edge { + int source; + int target; + }; + + const int NULL_VERTEX = -1; + const int NULL_EDGE = -1; + + public: + BoykovKolmogorovMaxFlow() = default; + + void resize(int num_verts, int num_edges) + { + m_graph_edges.resize(num_edges); + m_graph_out_edges.resize(num_verts); + m_cap_map.resize(num_edges); + m_res_cap_map.resize(num_edges); + m_rev_edge_map.resize(num_edges); + m_pre_map.resize(num_verts, 0); + m_tree_map.resize(num_verts, Color::gray); + m_dist_map.resize(num_verts, 0); + m_in_active_list_map.resize(num_verts, false); + m_has_parent_map.resize(num_verts, false); + m_time_map.resize(num_verts, 0); + } + + void set_edge(const int edge, + const int reverse_edge, + const int source, + const int target, + const int capacity) + { + m_graph_edges[edge] = Edge{source, target}; + m_graph_out_edges[source].push_back(edge); + m_rev_edge_map[edge] = reverse_edge; + + // Initialize flow to zero which means initializing + // the residual capacity equal to the capacity + m_cap_map[edge] = capacity; + m_res_cap_map[edge] = capacity; + } + + int edge_capacity(const int edge) + { + return m_cap_map[edge]; + } + + int edge_residual_capacity(const int edge) + { + return m_res_cap_map[edge]; + } + + const std::vector &vertex_out_edges(const int vertex) + { + return m_graph_out_edges[vertex]; + } + + int max_flow(int src, int sink) + { + m_source = src; + m_sink = sink; + + // init the search trees with the two terminals + m_tree_map[m_source] = Color::black; + m_tree_map[m_sink] = Color::white; + m_time_map[m_source] = 1; + m_time_map[m_sink] = 1; + + // augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK + augment_direct_paths(); + // start the main-loop + while (true) { + bool path_found; + int connecting_edge; + std::tie(connecting_edge, path_found) = grow(); // find a path from source to sink + if (!path_found) { + // we're finished, no more paths were found + break; + } + ++m_time; + augment(connecting_edge); // augment that path + adopt(); // rebuild search tree structure + } + return m_flow; + } + + protected: + int lookup_edge(int source, int target) + { + for (const int e : m_graph_out_edges[source]) { + if (m_graph_edges[e].target == target) { + return e; + } + } + return NULL_EDGE; + } + + void augment_direct_paths() + { + // in a first step, we augment all direct paths from + // source->NODE->sink and additionally paths from source->sink. This + // improves especially graphcuts for segmentation, as most of the + // nodes have source/sink connects but shouldn't have an impact on + // other maxflow problems (this is done in grow() anyway) + for (const int ei : m_graph_out_edges[m_source]) { + int from_source = ei; + int current_node = m_graph_edges[from_source].target; + if (current_node == m_sink) { + int cap = m_res_cap_map[from_source]; + m_res_cap_map[from_source] = 0; + m_flow += cap; + continue; + } + const int to_sink = lookup_edge(current_node, m_sink); + if (to_sink != NULL_EDGE) { + int cap_from_source = m_res_cap_map[from_source]; + int cap_to_sink = m_res_cap_map[to_sink]; + if (cap_from_source > cap_to_sink) { + m_tree_map[current_node] = Color::black; + add_active_node(current_node); + set_edge_to_parent(current_node, from_source); + m_dist_map[current_node] = 1; + m_time_map[current_node] = 1; + // add stuff to flow and update residuals. we dont need + // to update reverse_edges, as incoming/outgoing edges + // to/from source/sink don't count for max-flow + m_res_cap_map[from_source] = m_res_cap_map[from_source] - cap_to_sink; + m_res_cap_map[to_sink] = 0; + m_flow += cap_to_sink; + } + else if (cap_to_sink > 0) { + m_tree_map[current_node] = Color::white; + add_active_node(current_node); + set_edge_to_parent(current_node, to_sink); + m_dist_map[current_node] = 1; + m_time_map[current_node] = 1; + // add stuff to flow and update residuals. we dont need + // to update reverse_edges, as incoming/outgoing edges + // to/from source/sink don't count for max-flow + m_res_cap_map[to_sink] = m_res_cap_map[to_sink] - cap_from_source; + m_res_cap_map[from_source] = 0; + m_flow += cap_from_source; + } + } + else if (m_res_cap_map[from_source]) { + // there is no sink connect, so we can't augment this path, + // but to avoid adding m_source to the active nodes, we just + // activate this node and set the approciate things + m_tree_map[current_node] = Color::black; + set_edge_to_parent(current_node, from_source); + m_dist_map[current_node] = 1; + m_time_map[current_node] = 1; + add_active_node(current_node); + } + } + for (const int ei : m_graph_out_edges[m_sink]) { + int to_sink = m_rev_edge_map[ei]; + int current_node = m_graph_edges[to_sink].source; + if (m_res_cap_map[to_sink]) { + m_tree_map[current_node] = Color::white; + set_edge_to_parent(current_node, to_sink); + m_dist_map[current_node] = 1; + m_time_map[current_node] = 1; + add_active_node(current_node); + } + } + } + + /** + * Returns a pair of an edge and a boolean. if the bool is true, the + * edge is a connection of a found path from s->t , read "the link" and + * m_graph_edges[returnVal].source is the end of the path found in the + * source-tree m_graph_edges[returnVal].target is the beginning of the path found + * in the sink-tree + */ + std::pair grow() + { + assert(m_orphans.empty()); + int current_node; + while ((current_node = get_next_active_node()) != NULL_VERTEX) { // if there is one + assert(m_tree_map[current_node] != Color::gray && + (has_parent(current_node) || current_node == m_source || current_node == m_sink)); + + if (m_tree_map[current_node] == Color::black) { + // source tree growing + if (current_node != m_last_grow_vertex) { + m_last_grow_vertex = current_node; + m_last_grow_out_edge = 0; + } + const std::vector &out_edges = m_graph_out_edges[m_last_grow_vertex]; + for (; m_last_grow_out_edge < out_edges.size(); m_last_grow_out_edge++) { + int out_edge = out_edges[m_last_grow_out_edge]; + if (m_res_cap_map[out_edge] > 0) { // check if we have capacity left on this edge + int other_node = m_graph_edges[out_edge].target; + if (m_tree_map[other_node] == Color::gray) { // it's a free node + // aquire other node to our search tree + m_tree_map[other_node] = Color::black; + set_edge_to_parent(other_node, out_edge); // set us as parent + m_dist_map[other_node] = m_dist_map[current_node] + 1; // and update the + // distance-heuristic + m_time_map[other_node] = m_time_map[current_node]; + add_active_node(other_node); + } + else if (m_tree_map[other_node] == Color::black) { + // we do this to get shorter paths. check if we + // are nearer to the source as its parent is + if (is_closer_to_terminal(current_node, other_node)) { + set_edge_to_parent(other_node, out_edge); + m_dist_map[other_node] = m_dist_map[current_node] + 1; + m_time_map[other_node] = m_time_map[current_node]; + } + } + else { + assert(m_tree_map[other_node] == Color::white); + // kewl, found a path from one to the other + // search tree, return + // the connecting edge in src->sink dir + return std::make_pair(out_edge, true); + } + } + } // for all out-edges + } // source-tree-growing + else { + assert(m_tree_map[current_node] == Color::white); + if (current_node != m_last_grow_vertex) { + m_last_grow_vertex = current_node; + m_last_grow_out_edge = 0; + } + const std::vector &out_edges = m_graph_out_edges[m_last_grow_vertex]; + for (; m_last_grow_out_edge < out_edges.size(); m_last_grow_out_edge++) { + int in_edge = m_rev_edge_map[out_edges[m_last_grow_out_edge]]; + if (m_res_cap_map[in_edge] > 0) { // check if there is capacity left + int other_node = m_graph_edges[in_edge].source; + if (m_tree_map[other_node] == Color::gray) { // it's a free node + // aquire other node to our search tree + m_tree_map[other_node] = Color::white; + set_edge_to_parent(other_node, in_edge); // set us as parent + add_active_node(other_node); // activate that node + m_dist_map[other_node] = m_dist_map[current_node] + 1; // set its distance + m_time_map[other_node] = m_time_map[current_node]; // and time + } + else if (m_tree_map[other_node] == Color::white) { + if (is_closer_to_terminal(current_node, other_node)) { + // we are closer to the sink than its parent + // is, so we "adopt" him + set_edge_to_parent(other_node, in_edge); + m_dist_map[other_node] = m_dist_map[current_node] + 1; + m_time_map[other_node] = m_time_map[current_node]; + } + } + else { + assert(m_tree_map[other_node] == Color::black); + // kewl, found a path from one to the other + // search tree, + // return the connecting edge in src->sink dir + return std::make_pair(in_edge, true); + } + } + } // for all out-edges + } // sink-tree growing + + // all edges of that node are processed, and no more paths were + // found. + // remove if from the front of the active queue + finish_node(current_node); + } // while active_nodes not empty + + // no active nodes anymore and no path found, we're done + return std::make_pair(int(), false); + } + + /** + * augments path from s->t and updates residual graph + * m_graph_edges[e].source is the end of the path found in the source-tree + * m_graph_edges[e].target is the beginning of the path found in the sink-tree + * this phase generates orphans on satured edges, if the attached verts + * are from different search-trees orphans are ordered in distance to + * sink/source. first the farest from the source are front_inserted into + * the orphans list, and after that the sink-tree-orphans are + * front_inserted. when going to adoption stage the orphans are + * popped_front, and so we process the nearest verts to the terminals + * first + */ + void augment(int e) + { + assert(m_tree_map[m_graph_edges[e].target] == Color::white); + assert(m_tree_map[m_graph_edges[e].source] == Color::black); + assert(m_orphans.empty()); + + const int bottleneck = find_bottleneck(e); + // now we push the found flow through the path + // for each edge we saturate we have to look for the verts that + // belong to that edge, one of them becomes an orphans now process + // the connecting edge + m_res_cap_map[e] = m_res_cap_map[e] - bottleneck; + assert(m_res_cap_map[e] >= 0); + m_res_cap_map[m_rev_edge_map[e]] = m_res_cap_map[m_rev_edge_map[e]] + bottleneck; + + // now we follow the path back to the source + int current_node = m_graph_edges[e].source; + while (current_node != m_source) { + int pred = get_edge_to_parent(current_node); + m_res_cap_map[pred] = m_res_cap_map[pred] - bottleneck; + assert(m_res_cap_map[pred] >= 0); + m_res_cap_map[m_rev_edge_map[pred]] = m_res_cap_map[m_rev_edge_map[pred]] + bottleneck; + if (m_res_cap_map[pred] == 0) { + set_no_parent(current_node); + m_orphans.push_front(current_node); + } + current_node = m_graph_edges[pred].source; + } + // then go forward in the sink-tree + current_node = m_graph_edges[e].target; + while (current_node != m_sink) { + int pred = get_edge_to_parent(current_node); + m_res_cap_map[pred] = m_res_cap_map[pred] - bottleneck; + assert(m_res_cap_map[pred] >= 0); + m_res_cap_map[m_rev_edge_map[pred]] = m_res_cap_map[m_rev_edge_map[pred]] + bottleneck; + if (m_res_cap_map[pred] == 0) { + set_no_parent(current_node); + m_orphans.push_front(current_node); + } + current_node = m_graph_edges[pred].target; + } + // and add it to the max-flow + m_flow += bottleneck; + } + + /** + * returns the bottleneck of a s->t path (end_of_path is last vertex in + * source-tree, begin_of_path is first vertex in sink-tree) + */ + int find_bottleneck(int e) + { + int minimum_cap = m_res_cap_map[e]; + int current_node = m_graph_edges[e].source; + // first go back in the source tree + while (current_node != m_source) { + int pred = get_edge_to_parent(current_node); + minimum_cap = std::min(minimum_cap, m_res_cap_map[pred]); + current_node = m_graph_edges[pred].source; + } + // then go forward in the sink-tree + current_node = m_graph_edges[e].target; + while (current_node != m_sink) { + int pred = get_edge_to_parent(current_node); + minimum_cap = std::min(minimum_cap, m_res_cap_map[pred]); + current_node = m_graph_edges[pred].target; + } + return minimum_cap; + } + + /** + * rebuild search trees + * empty the queue of orphans, and find new parents for them or just + * drop them from the search trees + */ + void adopt() + { + while (!m_orphans.empty() || !m_child_orphans.empty()) { + int current_node; + if (m_child_orphans.empty()) { + // get the next orphan from the main-queue and remove it + current_node = m_orphans.front(); + m_orphans.pop_front(); + } + else { + current_node = m_child_orphans.front(); + m_child_orphans.pop(); + } + if (m_tree_map[current_node] == Color::black) { + // we're in the source-tree + int min_distance = (std::numeric_limits::max)(); + int new_parent_edge; + for (const int ei : m_graph_out_edges[current_node]) { + const int in_edge = m_rev_edge_map[ei]; + assert(m_graph_edges[in_edge].target == current_node); // we should be the target of + // this edge + if (m_res_cap_map[in_edge] > 0) { + int other_node = m_graph_edges[in_edge].source; + if (m_tree_map[other_node] == Color::black && has_source_connect(other_node)) { + if (m_dist_map[other_node] < min_distance) { + min_distance = m_dist_map[other_node]; + new_parent_edge = in_edge; + } + } + } + } + if (min_distance != (std::numeric_limits::max)()) { + set_edge_to_parent(current_node, new_parent_edge); + m_dist_map[current_node] = min_distance + 1; + m_time_map[current_node] = m_time; + } + else { + m_time_map[current_node] = 0; + for (const int ei : m_graph_out_edges[current_node]) { + int in_edge = m_rev_edge_map[ei]; + int other_node = m_graph_edges[in_edge].source; + if (m_tree_map[other_node] == Color::black && other_node != m_source) { + if (m_res_cap_map[in_edge] > 0) { + add_active_node(other_node); + } + if (has_parent(other_node) && + m_graph_edges[get_edge_to_parent(other_node)].source == current_node) + { + // we are the parent of that node + // it has to find a new parent, too + set_no_parent(other_node); + m_child_orphans.push(other_node); + } + } + } + m_tree_map[current_node] = Color::gray; + } // no parent found + } // source-tree-adoption + else { + // now we should be in the sink-tree, check that... + assert(m_tree_map[current_node] == Color::white); + int new_parent_edge; + int min_distance = (std::numeric_limits::max)(); + for (const int ei : m_graph_out_edges[current_node]) { + const int out_edge = ei; + if (m_res_cap_map[out_edge] > 0) { + const int other_node = m_graph_edges[out_edge].target; + if (m_tree_map[other_node] == Color::white && has_sink_connect(other_node)) { + if (m_dist_map[other_node] < min_distance) { + min_distance = m_dist_map[other_node]; + new_parent_edge = out_edge; + } + } + } + } + if (min_distance != (std::numeric_limits::max)()) { + set_edge_to_parent(current_node, new_parent_edge); + m_dist_map[current_node] = min_distance + 1; + m_time_map[current_node] = m_time; + } + else { + m_time_map[current_node] = 0; + for (const int ei : m_graph_out_edges[current_node]) { + const int out_edge = ei; + const int other_node = m_graph_edges[out_edge].target; + if (m_tree_map[other_node] == Color::white && other_node != m_sink) { + if (m_res_cap_map[out_edge] > 0) { + add_active_node(other_node); + } + if (has_parent(other_node) && + m_graph_edges[get_edge_to_parent(other_node)].target == current_node) + { + // we were it's parent, so it has to find a + // new one, too + set_no_parent(other_node); + m_child_orphans.push(other_node); + } + } + } + m_tree_map[current_node] = Color::gray; + } // no parent found + } // sink-tree adoption + } // while !orphans.empty() + } // adopt + + /** + * return next active vertex if there is one, otherwise a null_vertex + */ + int get_next_active_node() + { + while (true) { + if (m_active_nodes.empty()) { + return NULL_VERTEX; + } + int v = m_active_nodes.front(); + + // if it has no parent, this node can't be active (if its not + // source or sink) + if (!has_parent(v) && v != m_source && v != m_sink) { + m_active_nodes.pop(); + m_in_active_list_map[v] = false; + } + else { + assert(m_tree_map[v] == Color::black || m_tree_map[v] == Color::white); + return v; + } + } + } + + /** + * adds v as an active vertex, but only if its not in the list already + */ + void add_active_node(int v) + { + assert(m_tree_map[v] != Color::gray); + if (m_in_active_list_map[v]) { + if (m_last_grow_vertex == v) { + m_last_grow_vertex = NULL_VERTEX; + } + return; + } + + m_in_active_list_map[v] = true; + m_active_nodes.push(v); + } + + /** + * finish_node removes a node from the front of the active queue (its + * called in grow phase, if no more paths can be found using this node) + */ + void finish_node(int v) + { + assert(m_active_nodes.front() == v); + m_active_nodes.pop(); + m_in_active_list_map[v] = false; + m_last_grow_vertex = NULL_VERTEX; + } + + /** + * returns edge to parent vertex of v; + */ + int get_edge_to_parent(int v) const + { + return m_pre_map[v]; + } + + /** + * returns true if the edge stored in m_pre_map[v] is a valid entry + */ + bool has_parent(int v) const + { + return m_has_parent_map[v]; + } + + /** + * sets edge to parent vertex of v; + */ + void set_edge_to_parent(int v, int f_edge_to_parent) + { + assert(m_res_cap_map[f_edge_to_parent] > 0); + m_pre_map[v] = f_edge_to_parent; + m_has_parent_map[v] = true; + } + + /** + * removes the edge to parent of v (this is done by invalidating the + * entry an additional map) + */ + void set_no_parent(int v) + { + m_has_parent_map[v] = false; + } + + /** + * checks if vertex v has a connect to the sink-vertex (@var m_sink) + * @param v the vertex which is checked + * @return true if a path to the sink was found, false if not + */ + bool has_sink_connect(int v) + { + int current_distance = 0; + int current_vertex = v; + while (true) { + if (m_time_map[current_vertex] == m_time) { + // we found a node which was already checked this round. use + // it for distance calculations + current_distance += m_dist_map[current_vertex]; + break; + } + if (current_vertex == m_sink) { + m_time_map[m_sink] = m_time; + break; + } + if (has_parent(current_vertex)) { + // it has a parent, so get it + current_vertex = m_graph_edges[get_edge_to_parent(current_vertex)].target; + ++current_distance; + } + else { + // no path found + return false; + } + } + current_vertex = v; + while (m_time_map[current_vertex] != m_time) { + m_dist_map[current_vertex] = current_distance; + --current_distance; + m_time_map[current_vertex] = m_time; + current_vertex = m_graph_edges[get_edge_to_parent(current_vertex)].target; + } + return true; + } + + /** + * checks if vertex v has a connect to the source-vertex (@var m_source) + * @param v the vertex which is checked + * @return true if a path to the source was found, false if not + */ + bool has_source_connect(int v) + { + int current_distance = 0; + int current_vertex = v; + while (true) { + if (m_time_map[current_vertex] == m_time) { + // we found a node which was already checked this round. use + // it for distance calculations + current_distance += m_dist_map[current_vertex]; + break; + } + if (current_vertex == m_source) { + m_time_map[m_source] = m_time; + break; + } + if (has_parent(current_vertex)) { + // it has a parent, so get it + current_vertex = m_graph_edges[get_edge_to_parent(current_vertex)].source; + ++current_distance; + } + else { + // no path found + return false; + } + } + current_vertex = v; + while (m_time_map[current_vertex] != m_time) { + m_dist_map[current_vertex] = current_distance; + --current_distance; + m_time_map[current_vertex] = m_time; + current_vertex = m_graph_edges[get_edge_to_parent(current_vertex)].source; + } + return true; + } + + /** + * returns true, if p is closer to a terminal than q + */ + bool is_closer_to_terminal(int p, int q) + { + // checks the timestamps first, to build no cycles, and after that + // the real distance + return (m_time_map[q] <= m_time_map[p] && m_dist_map[q] > m_dist_map[p] + 1); + } + + // Member variables + std::vector m_graph_edges; + std::vector> m_graph_out_edges; + + std::vector m_cap_map; + std::vector m_res_cap_map; + std::vector m_rev_edge_map; + std::vector m_pre_map; // stores paths found in the growth stage + std::vector m_tree_map; // maps each vertex into white, black or gray + std::vector m_dist_map; // stores distance to source/sink nodes + + int m_source; + int m_sink; + + std::queue m_active_nodes; + std::vector m_in_active_list_map; + + std::list m_orphans; + std::queue m_child_orphans; // we use a second queuqe for child orphans, as + // they are FIFO processed + std::vector m_has_parent_map; + + std::vector m_time_map; // timestamp of each node, used for + // sink/source-path calculations + int m_flow = 0; + int m_time = 1; + + int m_last_grow_vertex = NULL_VERTEX; + int m_last_grow_out_edge = 0; +}; + +} // namespace qflow diff --git a/extern/quadriflow/src/config.hpp b/extern/quadriflow/src/config.hpp index bf597ad0f39..842b885a209 100644 --- a/extern/quadriflow/src/config.hpp +++ b/extern/quadriflow/src/config.hpp @@ -1,11 +1,6 @@ #ifndef CONFIG_H_ #define CONFIG_H_ -/* Workaround a bug in boost 1.68, until we upgrade to a newer version. */ -#if defined(__clang__) && defined(WIN32) - #include - using namespace boost; -#endif // Move settings to cmake to make CMake happy :) // #define WITH_SCALE diff --git a/extern/quadriflow/src/flow.hpp b/extern/quadriflow/src/flow.hpp index ab4a01cb4ed..a77f7ae6639 100644 --- a/extern/quadriflow/src/flow.hpp +++ b/extern/quadriflow/src/flow.hpp @@ -7,17 +7,12 @@ #include #include "config.hpp" - -#include -#include -#include -#include +#include "../patches/boykov_kolmogorov_max_flow.hpp" #include #include #include -using namespace boost; using namespace Eigen; namespace qflow { @@ -34,78 +29,52 @@ class MaxFlowHelper { class BoykovMaxFlowHelper : public MaxFlowHelper { public: - typedef int EdgeWeightType; - typedef adjacency_list_traits Traits; - // clang-format off - typedef adjacency_list < vecS, vecS, directedS, - property < vertex_name_t, std::string, - property < vertex_index_t, long, - property < vertex_color_t, boost::default_color_type, - property < vertex_distance_t, long, - property < vertex_predecessor_t, Traits::edge_descriptor > > > > >, - - property < edge_capacity_t, EdgeWeightType, - property < edge_residual_capacity_t, EdgeWeightType, - property < edge_reverse_t, Traits::edge_descriptor > > > > Graph; - // clang-format on - - public: - BoykovMaxFlowHelper() { rev = get(edge_reverse, g); } - void resize(int n, int m) { - vertex_descriptors.resize(n); - for (int i = 0; i < n; ++i) vertex_descriptors[i] = add_vertex(g); + BoykovMaxFlowHelper() = default; + void resize(int n, int m) override { + num_verts = n; + num_edges = 0; + flow.resize(num_verts, m * 2); } - int compute() { - EdgeWeightType flow = - boykov_kolmogorov_max_flow(g, vertex_descriptors.front(), vertex_descriptors.back()); - return flow; + int compute() override { + return flow.max_flow(0, num_verts - 1); } - void addDirectEdge(Traits::vertex_descriptor& v1, Traits::vertex_descriptor& v2, - property_map::type& rev, const int capacity, - const int inv_capacity, Graph& g, Traits::edge_descriptor& e1, - Traits::edge_descriptor& e2) { - e1 = add_edge(v1, v2, g).first; - e2 = add_edge(v2, v1, g).first; - put(edge_capacity, g, e1, capacity); - put(edge_capacity, g, e2, inv_capacity); - - rev[e1] = e2; - rev[e2] = e1; - } - void addEdge(int x, int y, int c, int rc, int v, int cost = 1) { - Traits::edge_descriptor e1, e2; - addDirectEdge(vertex_descriptors[x], vertex_descriptors[y], rev, c, rc, g, e1, e2); + void addEdge(int x, int y, int c, int rc, int v, int cost = 1) override { + const int e1 = num_edges++; + const int e2 = num_edges++; + flow.set_edge(e1, e2, x, y, c); + flow.set_edge(e2, e1, y, x, rc); if (v != -1) { - edge_to_variables[e1] = std::make_pair(v, -1); - edge_to_variables[e2] = std::make_pair(v, 1); + edge_to_variables.emplace_back(v, -1); + edge_to_variables.emplace_back(v, 1); + } + else { + edge_to_variables.emplace_back(-1, -1); + edge_to_variables.emplace_back(-1, -1); } } - void applyTo(std::vector& edge_diff) { - property_map::type capacity = get(edge_capacity, g); - property_map::type residual_capacity = - get(edge_residual_capacity, g); - - graph_traits::vertex_iterator u_iter, u_end; - graph_traits::out_edge_iterator ei, e_end; - for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) - for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) - if (capacity[*ei] > 0) { - int flow = (capacity[*ei] - residual_capacity[*ei]); + void applyTo(std::vector& edge_diff) override { + for (int vert = 0; vert < num_verts; vert++) { + for (int edge : flow.vertex_out_edges(vert)) { + const int capacity = flow.edge_capacity(edge); + const int residual_capacity = flow.edge_residual_capacity(edge); + if (capacity > 0) { + int flow = (capacity - residual_capacity); if (flow > 0) { - auto it = edge_to_variables.find(*ei); - if (it != edge_to_variables.end()) { - edge_diff[it->second.first / 2][it->second.first % 2] += - it->second.second * flow; + std::pair e2v = edge_to_variables[edge]; + if (e2v.first != -1) { + edge_diff[e2v.first / 2][e2v.first % 2] += e2v.second * flow; } } } + } + } } private: - Graph g; - property_map::type rev; - std::vector vertex_descriptors; - std::map> edge_to_variables; + BoykovKolmogorovMaxFlow flow; + std::vector> edge_to_variables; + int num_verts = 0; + int num_edges = 0; }; class NetworkSimplexFlowHelper : public MaxFlowHelper { diff --git a/extern/quadriflow/src/post-solver.cpp b/extern/quadriflow/src/post-solver.cpp index 6027ddd2eb7..ccefd1557d8 100644 --- a/extern/quadriflow/src/post-solver.cpp +++ b/extern/quadriflow/src/post-solver.cpp @@ -5,7 +5,9 @@ // Created by Jingwei on 2/5/18. // #include +#ifdef POST_SOLVER #include +#endif #include #include #include diff --git a/intern/quadriflow/CMakeLists.txt b/intern/quadriflow/CMakeLists.txt index b253d1feb0e..5aef735fc7e 100644 --- a/intern/quadriflow/CMakeLists.txt +++ b/intern/quadriflow/CMakeLists.txt @@ -8,7 +8,6 @@ set(INC set(INC_SYS ../../extern/quadriflow/src - ${BOOST_INCLUDE_DIR} ${EIGEN3_INCLUDE_DIRS} ) @@ -20,7 +19,6 @@ set(SRC set(LIB PRIVATE bf::intern::guardedalloc extern_quadriflow - ${BOOST_LIBRARIES} ) if(WIN32)