enh: enhance labelisation_seeded_watershed performance with shortcut (#260)#296
Open
lapertor wants to merge 1 commit into
Open
enh: enhance labelisation_seeded_watershed performance with shortcut (#260)#296lapertor wants to merge 1 commit into
lapertor wants to merge 1 commit into
Conversation
Member
|
Thank you for addressing this. As this is a performance-related issue, can you also provide some benchmark to assess if the optimization really provides an improvement? |
Contributor
Author
|
I ran instrumented benchmarks to isolate the performance of the main loop (excluding Setup:
Results (loop only, in ms):
It shows clear gains when the seed ratio is high enough. For very low seed ratios, the overhead of the counter is negligible (within noise). Benchmark scriptimport higra as hg
import numpy as np
import time
def make_block_seeds(shape, block_shape):
h, w = shape
bh, bw = block_shape
seeds = np.zeros(h * w, dtype=np.int32)
label = 1
for bi in range(0, h, bh):
for bj in range(0, w, bw):
ci = min(bi + bh // 2, h - 1)
cj = min(bj + bw // 2, w - 1)
seeds[ci * w + cj] = label
label += 1
return seeds
def make_block_edge_weights(graph, shape, block_shape):
h, w = shape
bh, bw = block_shape
sources, targets = graph.edge_list()
weights = np.zeros(graph.num_edges(), dtype=np.float32)
for i in range(graph.num_edges()):
u = sources[i]
v = targets[i]
ui, uj = divmod(u, w)
vi, vj = divmod(v, w)
same_block = (ui // bh == vi // bh) and (uj // bw == vj // bw)
weights[i] = 0.0 if same_block else 1.0
return weights
def _benchmark(shape=(2000, 2000), block_shape=(100, 100), n_runs=5):
graph = hg.get_4_adjacency_graph(shape)
seeds = make_block_seeds(shape, block_shape)
edge_weights = make_block_edge_weights(graph, shape, block_shape)
print(f"Shape: {shape}")
print(f"Block shape: {block_shape}")
print(f"Num vertices: {graph.num_vertices()}")
print(f"Num edges: {graph.num_edges()}")
print(f"Num seeds: {np.count_nonzero(seeds)}")
times = []
for _ in range(n_runs):
t0 = time.perf_counter()
labels = hg.labelisation_seeded_watershed(graph, edge_weights, seeds)
t1 = time.perf_counter()
times.append(t1 - t0)
print(f"Times: {[round(t, 4) for t in times]}")
print(f"Mean: {np.mean(times):.4f} s")
print(f"Std: {np.std(times):.4f} s")
print(f"Labels range: {labels.min()} to {labels.max()}")
def benchmark():
shape = (2000, 2000)
for block_shape in [(2000, 2000), (100, 100), (10, 10), (4, 4), (2, 2), (1, 1)]:
print("=" * 60)
_benchmark(shape=shape, block_shape=block_shape)
print()
if __name__ == "__main__":
benchmark()Old implementation/***************************************************************************
* Copyright ESIEE Paris (2018) *
* *
* Contributor(s) : Benjamin Perret *
* *
* Distributed under the terms of the CECILL-B License. *
* *
* The full license is in the file LICENSE, distributed with this software. *
****************************************************************************/
#pragma once
#include "../graph.hpp"
#include "../structure/array.hpp"
#include "higra/structure/unionfind.hpp"
#include "higra/sorting.hpp"
#include <vector>
#include <stack>
#include <chrono>
namespace hg {
/**
* Linear time watershed cut algorithm.
*
* Jean Cousty, Gilles Bertrand, Laurent Najman, Michel Couprie. Watershed Cuts: Minimum Spanning
* Forests and the Drop of Water Principle. IEEE Transactions on Pattern Analysis and Machine
* Intelligence, Institute of Electrical and Electronics Engineers, 2009, 31 (8), pp.1362-1374.
* @tparam graph_t
* @tparam T
* @param graph
* @param xedge_weights
* @return array of labels on graph vertices, numbered from 1 to n with n the number of minima
*/
template<typename graph_t, typename T>
auto
labelisation_watershed(const graph_t &graph, const xt::xexpression<T> &xedge_weights) {
HG_TRACE();
auto &edge_weights = xedge_weights.derived_cast();
hg_assert_edge_weights(graph, edge_weights);
hg_assert_1d_array(edge_weights);
using value_type = typename T::value_type;
using vertex_t = typename graph_traits<graph_t>::vertex_descriptor;
auto fminus = array_1d<value_type>::from_shape({graph.num_vertices()});
for (auto v: vertex_iterator(graph)) {
auto minValue = (std::numeric_limits<value_type>::max)();
for (auto e: out_edge_iterator(v, graph)) {
minValue = (std::min)(minValue, edge_weights(e));
}
fminus[v] = minValue;
}
auto no_label = (std::numeric_limits<index_t>::max)();
auto labels = array_1d<index_t>::from_shape({graph.num_vertices()});
std::fill(labels.begin(), labels.end(), no_label);
auto notInL = array_1d<bool>::from_shape({graph.num_vertices()});
std::fill(notInL.begin(), notInL.end(), true);
std::vector<vertex_t> L;
std::vector<vertex_t> LL;
auto stream = [&L, &LL, &graph, &edge_weights, &fminus, ¬InL, &labels, no_label](vertex_t x) {
L.clear();
LL.clear();
L.push_back(x);
LL.push_back(x);
notInL[x] = false;
while (!LL.empty()) {
auto y = LL[LL.size() - 1];
LL.pop_back();
for (auto e: out_edge_iterator(y, graph)) {
auto adjacent_vertex = target(e, graph);
if (notInL[adjacent_vertex] && edge_weights(e) == fminus[y]) {
if (labels[adjacent_vertex] != no_label) {
return labels[adjacent_vertex];
} else if (fminus[adjacent_vertex] < fminus[y]) {
L.push_back(adjacent_vertex);
notInL[adjacent_vertex] = false;
LL.clear();
LL.push_back(adjacent_vertex);
break; // stop breadth_first
} else {
L.push_back(adjacent_vertex);
notInL[adjacent_vertex] = false;
LL.push_back(adjacent_vertex);
}
}
}
}
return no_label;
};
index_t num_labs = 0;
for (auto v: vertex_iterator(graph)) {
if (labels[v] == no_label) {
auto res = stream(v);
if (res == no_label) {
num_labs++;
for (auto x: L) {
labels[x] = num_labs;
notInL[x] = true;
}
} else {
for (auto x: L) {
labels[x] = res;
notInL[x] = true;
}
}
}
}
return labels;
};
template<typename graph_t, typename T1, typename T2>
auto labelisation_seeded_watershed(
const graph_t &graph,
const xt::xexpression<T1> &xedge_weights,
const xt::xexpression<T2> &xvertex_seeds,
const typename T2::value_type background_label = 0) {
HG_TRACE();
auto &edge_weights = xedge_weights.derived_cast();
auto &vertex_seeds = xvertex_seeds.derived_cast();
hg_assert_edge_weights(graph, edge_weights);
hg_assert_node_weights(graph, vertex_seeds);
hg_assert_1d_array(edge_weights);
hg_assert_1d_array(vertex_seeds);
using label_type = typename T2::value_type;
array_1d<index_t> sorted_edges_indices = stable_arg_sort(edge_weights);
auto t0 = std::chrono::high_resolution_clock::now();
index_t num_nodes = num_vertices(graph);
index_t num_edges = sorted_edges_indices.size();
union_find uf(num_nodes);
array_1d<label_type> labels = vertex_seeds;
for (index_t i = 0; i < num_edges; i++) {
auto ei = sorted_edges_indices[i];
auto e = edge_from_index(ei, graph);
auto c1 = uf.find(source(e, graph));
auto c2 = uf.find(target(e, graph));
if (c1 != c2 && (labels(c1) == background_label || labels(c2) == background_label)) {
if (labels(c1) == background_label) {
labels(c1) = labels(c2);
} else {
labels(c2) = labels(c1);
}
uf.link(c1, c2);
}
}
for (index_t i = 0; i < num_nodes; i++) {
if (labels(i) == background_label) {
labels(i) = labels(uf.find(i));
}
}
auto t1 = std::chrono::high_resolution_clock::now();
std::cerr << "watershed loop: "
<< std::chrono::duration<double, std::milli>(t1 - t0).count()
<< " ms" << std::endl;
return labels;
};
}Output of tests on the old implementationNew implementation/***************************************************************************
* Copyright ESIEE Paris (2018) *
* *
* Contributor(s) : Benjamin Perret *
* *
* Distributed under the terms of the CECILL-B License. *
* *
* The full license is in the file LICENSE, distributed with this software. *
****************************************************************************/
#pragma once
#include "../graph.hpp"
#include "../structure/array.hpp"
#include "higra/structure/unionfind.hpp"
#include "higra/sorting.hpp"
#include <vector>
#include <stack>
#include <chrono>
namespace hg {
/**
* Linear time watershed cut algorithm.
*
* Jean Cousty, Gilles Bertrand, Laurent Najman, Michel Couprie. Watershed Cuts: Minimum Spanning
* Forests and the Drop of Water Principle. IEEE Transactions on Pattern Analysis and Machine
* Intelligence, Institute of Electrical and Electronics Engineers, 2009, 31 (8), pp.1362-1374.
* @tparam graph_t
* @tparam T
* @param graph
* @param xedge_weights
* @return array of labels on graph vertices, numbered from 1 to n with n the number of minima
*/
template<typename graph_t, typename T>
auto
labelisation_watershed(const graph_t &graph, const xt::xexpression<T> &xedge_weights) {
HG_TRACE();
auto &edge_weights = xedge_weights.derived_cast();
hg_assert_edge_weights(graph, edge_weights);
hg_assert_1d_array(edge_weights);
using value_type = typename T::value_type;
using vertex_t = typename graph_traits<graph_t>::vertex_descriptor;
auto fminus = array_1d<value_type>::from_shape({graph.num_vertices()});
for (auto v: vertex_iterator(graph)) {
auto minValue = (std::numeric_limits<value_type>::max)();
for (auto e: out_edge_iterator(v, graph)) {
minValue = (std::min)(minValue, edge_weights(e));
}
fminus[v] = minValue;
}
auto no_label = (std::numeric_limits<index_t>::max)();
auto labels = array_1d<index_t>::from_shape({graph.num_vertices()});
std::fill(labels.begin(), labels.end(), no_label);
auto notInL = array_1d<bool>::from_shape({graph.num_vertices()});
std::fill(notInL.begin(), notInL.end(), true);
std::vector<vertex_t> L;
std::vector<vertex_t> LL;
auto stream = [&L, &LL, &graph, &edge_weights, &fminus, ¬InL, &labels, no_label](vertex_t x) {
L.clear();
LL.clear();
L.push_back(x);
LL.push_back(x);
notInL[x] = false;
while (!LL.empty()) {
auto y = LL[LL.size() - 1];
LL.pop_back();
for (auto e: out_edge_iterator(y, graph)) {
auto adjacent_vertex = target(e, graph);
if (notInL[adjacent_vertex] && edge_weights(e) == fminus[y]) {
if (labels[adjacent_vertex] != no_label) {
return labels[adjacent_vertex];
} else if (fminus[adjacent_vertex] < fminus[y]) {
L.push_back(adjacent_vertex);
notInL[adjacent_vertex] = false;
LL.clear();
LL.push_back(adjacent_vertex);
break; // stop breadth_first
} else {
L.push_back(adjacent_vertex);
notInL[adjacent_vertex] = false;
LL.push_back(adjacent_vertex);
}
}
}
}
return no_label;
};
index_t num_labs = 0;
for (auto v: vertex_iterator(graph)) {
if (labels[v] == no_label) {
auto res = stream(v);
if (res == no_label) {
num_labs++;
for (auto x: L) {
labels[x] = num_labs;
notInL[x] = true;
}
} else {
for (auto x: L) {
labels[x] = res;
notInL[x] = true;
}
}
}
}
return labels;
};
template<typename graph_t, typename T1, typename T2>
auto labelisation_seeded_watershed(
const graph_t &graph,
const xt::xexpression<T1> &xedge_weights,
const xt::xexpression<T2> &xvertex_seeds,
const typename T2::value_type background_label = 0) {
HG_TRACE();
auto &edge_weights = xedge_weights.derived_cast();
auto &vertex_seeds = xvertex_seeds.derived_cast();
hg_assert_edge_weights(graph, edge_weights);
hg_assert_node_weights(graph, vertex_seeds);
hg_assert_1d_array(edge_weights);
hg_assert_1d_array(vertex_seeds);
using label_type = typename T2::value_type;
array_1d<index_t> sorted_edges_indices = stable_arg_sort(edge_weights);
auto t0 = std::chrono::high_resolution_clock::now();
index_t num_nodes = num_vertices(graph);
index_t num_edges = sorted_edges_indices.size();
union_find uf(num_nodes);
array_1d<label_type> labels = vertex_seeds;
index_t num_background_components = 0;
for (index_t i = 0; i < num_nodes; i++) {
if (labels(i) == background_label) {
num_background_components++;
}
}
if (num_background_components > 0) {
for (index_t i = 0; i < num_edges; i++) {
auto ei = sorted_edges_indices[i];
auto e = edge_from_index(ei, graph);
auto c1 = uf.find(source(e, graph));
auto c2 = uf.find(target(e, graph));
if (c1 != c2 && (labels(c1) == background_label || labels(c2) == background_label)) {
if (labels(c1) == background_label) {
labels(c1) = labels(c2);
} else {
labels(c2) = labels(c1);
}
uf.link(c1, c2);
num_background_components--;
if (num_background_components == 0) {
break;
}
}
}
}
for (index_t i = 0; i < num_nodes; i++) {
if (labels(i) == background_label) {
labels(i) = labels(uf.find(i));
}
}
auto t1 = std::chrono::high_resolution_clock::now();
std::cerr << "watershed loop: "
<< std::chrono::duration<double, std::milli>(t1 - t0).count()
<< " ms" << std::endl;
return labels;
};
}Output of tests on the new implementation |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
This PR addresses Issue #260 by implementing the first proposed optimization (early stopping) for
labelisation_seeded_watershed.Changes made
num_background_componentsto track the number of vertices initially set to thebackground_label.breakin the main edge processing loop as soon asnum_background_components == 0.if (num_background_components > 0)to handle the edge case where all vertices are already seeded or no background exists.Partially fixes #260, let's explore idea 2 now (replacing the relabeling).