Commit 9db0bcc9 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Add mumford-shah computation on the tree.

	*  apps/CMakeLists.txt,
	*  apps/attributes/gradient_magnitude.hpp: Add gradient on 0F too.
	*  apps/mumford_shah_on_tree/CMakeLists.txt,
	*  apps/mumford_shah_on_tree/mumford_shah_on_tree.cpp,
	*  apps/mumford_shah_on_tree/mumford_shah_on_tree.hpp,
	*  apps/mumford_shah_on_tree/mumford_shah_on_tree_full.cpp: New.
parent be1cd23f
......@@ -6,6 +6,8 @@ add_subdirectory(maxtree_comparison)
add_subdirectory(attributes)
add_subdirectory(simplification)
add_subdirectory(clattice)
add_subdirectory(g2)
add_subdirectory(mumford_shah_on_tree)
add_subdirectory(tosgui)
add_subdirectory(saliency)
add_subdirectory(supervised-gui)
......
......@@ -18,7 +18,7 @@ namespace mln
mln_entering("mln::compute_gradient_magnitude");
image2d<float> grad;
resize(grad, f);
resize(grad, f).init(1.0);
box2d dom = f.domain();
sbox2d d2f = {dom.pmin - point2d{2,2}, dom.pmax, point2d{2,2}};
......@@ -29,6 +29,7 @@ namespace mln
{
grad.at(p + point2d{0,1}) = 1 - l1norm(f.at(p) - f.at(p + point2d{0,2})) / lmax;
grad.at(p + point2d{1,0}) = 1 - l1norm(f.at(p) - f.at(p + point2d{2,0})) / lmax;
grad.at(p + point2d{1,1}) = 1 - l1norm(f.at(p) - f.at(p + point2d{2,2})) / lmax;
}
auto res = compute_attribute_on_contour(tree, grad, accu::features::mean<> ());
......
link_libraries(-lfreeimage)
add_executable(mumford_shah_on_tree mumford_shah_on_tree.cpp)
add_executable(mumford_shah_on_tree_full mumford_shah_on_tree_full.cpp
${CMAKE_SOURCE_DIR}/apps/g2/satmaxtree.cpp
${CMAKE_SOURCE_DIR}/apps/g2/compute_g2.cpp)
target_link_libraries(mumford_shah_on_tree_full ${TBB_LIBRARIES})
\ No newline at end of file
#include <mln/core/image/image2d.hpp>
#include <mln/morpho/component_tree/io.hpp>
#include <mln/accu/accumulators/variance.hpp>
#include <mln/accu/accumulators/accu_as_it.hpp>
#include <mln/io/imread.hpp>
#include <mln/io/imsave.hpp>
#include <apps/tos/Kinterpolate.hpp>
#include "mumford_shah_on_tree.hpp"
typedef mln::morpho::component_tree<unsigned, mln::image2d<unsigned> > tree_t;
tree_t::vertex_id_t
find_canonical(mln::property_map<tree_t, tree_t::vertex_id_t>& newpar,
tree_t::vertex_id_t x)
{
if (newpar[x] == tree_t::npos()) // x is alive
return x;
else
return newpar[x] = find_canonical(newpar, newpar[x]);
}
int main(int argc, char** argv)
{
if (argc < 6) {
std::cerr << "Usage: " << argv[0] << " input.ppm input.tree lambda grain output.ppm\n"
;
std::exit(1);
}
const char* input_path = argv[1];
const char* tree_path = argv[2];
float lambda = std::atof(argv[3]);
unsigned grain = std::atoi(argv[4]);
const char* output_path = argv[5];
using namespace mln;
image2d<rgb8> f, F;
tree_t tree;
io::imread(input_path, f);
morpho::load(tree_path, tree);
F = Kadjust_to(f, tree._get_data()->m_pmap.domain());
// Filter the grain first
{
accu::accumulators::accu_if< accu::accumulators::count<>,
K1::is_face_2_t,
point2d > counter;
auto areamap = morpho::paccumulate(tree, F, counter);
auto pred = make_functional_property_map<tree_t::vertex_id_t>([&areamap, grain](tree_t::vertex_id_t x) {
return areamap[x] > grain;
});
morpho::filter_direct_inplace(tree, pred);
tree.shrink_to_fit();
}
mumford_shah_on_tree(tree, F, lambda);
F = Kadjust_to(F, f.domain());
io::imsave(F, output_path);
}
#ifndef MUMFORD_SHAH_ON_TREE_HPP
# define MUMFORD_SHAH_ON_TREE_HPP
# include <mln/morpho/component_tree/accumulate.hpp>
# include <mln/morpho/component_tree/reconstruction.hpp>
# include <mln/morpho/component_tree/filtering.hpp>
# include <mln/accu/accumulators/variance.hpp>
# include <mln/accu/accumulators/mean.hpp>
# include <mln/accu/accumulators/accu_as_it.hpp>
# include <apps/attributes/gradient_magnitude.hpp>
/// \brief Perform the MS simplification
///
/// \param tree The input tree
/// \param[inout] F The input image (interpolated to match tree's domain)
/// \param lambda Mumford-shash regularization parameter
template <class T, class I>
void
mumford_shah_on_tree(const T& tree,
I& F,
double lambda);
/***************************************/
/*** Implementation ****/
/***************************************/
namespace internal
{
typedef mln::morpho::component_tree<unsigned, mln::image2d<unsigned> > tree_t;
tree_t::vertex_id_t
find_canonical(mln::property_map<tree_t, tree_t::vertex_id_t>& newpar,
tree_t::vertex_id_t x)
{
if (newpar[x] == tree_t::npos()) // x is alive
return x;
else
return newpar[x] = find_canonical(newpar, newpar[x]);
}
}
template <class T, class I>
void
mumford_shah_on_tree(T& tree,
I& F,
double lambda)
{
using namespace mln;
typedef mln::morpho::component_tree<unsigned, mln::image2d<unsigned> > tree_t;
auto A = morpho::vaccumulate_proper(tree, F,
accu::accumulators::accu_as_it<
accu::accumulators::variance< rgb8, rgb<double>, rgb<double> >
> ());
auto B = compute_attribute_on_contour(tree, F, accu::features::count<>());
// Sort the node by gradient
auto mgrad = compute_gradient_magnitude(tree, F); // note: minima = strong gradient
std::vector<tree_t::vertex_id_t> S;
S.reserve(tree.size());
mln_foreach(tree_t::node_type x, tree.nodes_without_root())
S.push_back(x.id());
std::sort(S.begin(), S.end(), [&mgrad](tree_t::vertex_id_t x, tree_t::vertex_id_t y) {
return mgrad[x] > mgrad[y];
});
// Energy function
auto denergy = [lambda, &A, &B](tree_t::node_type x, tree_t::node_type q) {
auto tmp = A[q];
double before = (tmp.to_result() * accu::extractor::count(tmp));
tmp.take(A[x]);
double after = (tmp.to_result() * accu::extractor::count(tmp));
//std::cout << "before: " << before << " after: " << after << std::endl;
double delta = -lambda * B[x] + (after - before);
//std::cout << "delta: " << delta << std::endl;
return delta;
};
// Minimization glutone
std::cout << "Before: " << tree.size() << std::endl;
unsigned k = tree.size();
property_map<tree_t, bool> alive(tree, true);
property_map<tree_t, tree_t::vertex_id_t> newpar(tree, tree_t::npos());
for (tree_t::vertex_id_t i : S)
{
tree_t::node_type x = tree.get_node(i);
tree_t::node_type q = tree.get_node(::internal::find_canonical(newpar, x.get_parent_id()));
if (denergy(x,q) < 0)
{
alive[x] = false;
A[q].take(A[x]);
newpar[x] = q.id();
--k;
}
}
std::cout << "After: " << k << std::endl;
// Reconstruction
{
morpho::filter_direct_inplace(tree, alive);
tree.shrink_to_fit();
auto vmap = morpho::vaccumulate_proper(tree, F, accu::features::mean<> ());
morpho::reconstruction(tree, vmap, F);
}
}
#endif // ! MUMFORD_SHAH_ON_TREE_HPP
#include <mln/core/image/image2d.hpp>
#include <mln/core/neighb2d.hpp>
#include <mln/morpho/tos/ctos.hpp>
#include <mln/morpho/component_tree/compute_depth.hpp>
#include <mln/morpho/component_tree/reconstruction.hpp>
#include <mln/morpho/component_tree/accumulate.hpp>
#include <tbb/parallel_for.h>
#include <tbb/task_scheduler_init.h>
#include <mln/io/imread.hpp>
#include <mln/io/imsave.hpp>
#include <mln/core/dontcare.hpp>
#include <boost/graph/dag_shortest_paths.hpp>
#include <boost/graph/transpose_graph.hpp>
#include <boost/property_map/function_property_map.hpp>
#include <apps/tos/addborder.hpp>
#include <apps/tos/Kinterpolate.hpp>
#include <apps/tos/croutines.hpp>
#include <apps/g2/compute_g2.hpp>
#include <apps/g2/satmaxtree.hpp>
#include "mumford_shah_on_tree.hpp"
// Compute the depth attribute of each graph node
boost::vector_property_map<unsigned>
compute_graph_depth(const MyGraph& g)
{
mln_entering("Compute graph depth");
boost::vector_property_map<unsigned> depth(boost::num_vertices(g));
auto one = [](mln::dontcare_t) -> int{ return 1; };
auto w = boost::make_function_property_map<MyGraph::edge_descriptor, int, decltype(one)>(one);
MyGraph::vertex_descriptor root = boost::vertex(0, g);
depth[root] = 0;
MyGraph gT;
boost::transpose_graph(g, gT);
boost::dag_shortest_paths(gT, root, boost::weight_map(w)
.distance_map(depth)
.distance_compare(std::greater<int> ())
.distance_inf(-1)
.distance_zero(0)
);
mln_exiting();
return depth;
}
// Compute the per-pixel attribute and reconstruct
template <class ValueMap>
void
write_vmap_to_image(const tree_t* t, const tlink_t* tlink,
const ValueMap& vmap, mln::image2d<mln::uint16>& out)
{
mln_foreach(auto px, out.pixels())
{
unsigned w = 0;
for (int i = 0; i < NTREE; ++i)
{
tree_t::node_type tnode = t[i].get_node_at(px.index());
MyGraph::vertex_descriptor gnode = tlink[i][tnode];
w = std::max(w, vmap[gnode]);
}
px.val() = w;
}
}
/// \brief Remove non-2F from the tree
template <class P>
mln::morpho::component_tree<P, mln::image2d<P> >
tree_keep_2F(const mln::morpho::component_tree<P, mln::image2d<P> >& tree)
{
using namespace mln;
morpho::component_tree<P, image2d<P> > out;
auto newdata = out._get_data();
auto olddata = tree._get_data();
// 1. Copy the point2node map
box2d olddom = olddata->m_pmap.domain();
box2d dom;
dom.pmin = olddom.pmin / 2;
dom.pmax = (olddom.pmax + 1) / 2;
newdata->m_pmap.resize(dom);
copy(olddata->m_pmap | sbox2d{olddom.pmin, olddom.pmax, {2,2}},
newdata->m_pmap);
// 2. Copy the node
newdata->m_nodes = olddata->m_nodes;
// 3. Copy the point set and update node first point index/
newdata->m_S.resize(dom.size());
unsigned j = 0;
for (unsigned i = 0; i < olddata->m_S.size(); ++i)
{
P p = olddata->m_S[i];
point2d pt = olddata->m_pmap.point_at_index(p);
if (K1::is_face_2(pt))
{
newdata->m_S[j] = newdata->m_pmap.index_of_point(pt/2);
auto node = tree.get_node_at(p);
if (node.get_first_point_id() == i)
newdata->m_nodes[node.id()].m_point_index = j;
++j;
}
}
// 4. Do not forget the sentinel
newdata->m_nodes[out.npos()].m_point_index = j;
return out.get_subtree(tree.get_root_id());
}
void usage(char** argv)
{
std::cout << "Usage: " << argv[0] << " input[rgb] α₀ α₁ λ output[rgb]\n"
"α₀\tGrain filter size before merging trees (0 to disable)\n"
"α₁\tGrain filter size on the color ToS (0 to disable)\n"
\tMumford-shah regularisation weight (e.g. 5000)\n";
std::exit(1);
}
int main(int argc, char** argv)
{
if (argc < 5)
usage(argv);
const char* input_path = argv[1];
int a0 = std::atoi(argv[2]);
int a1 = std::atoi(argv[3]);
int lambda = std::atoi(argv[4]);
const char* output_path = argv[5];
tbb::task_scheduler_init init;
// 1. Compute the individual ToS
using namespace mln;
typedef rgb8 V;
image2d<V> ima;
io::imread(input_path, ima);
image2d<V> f = addborder(ima, lexicographicalorder_less<value_t>());
image2d<V> F = interpolate_k1(f);
// COmpute the color tree of shapes T
tree_t T;
{
// 1. Compute the marginal ToS and filter them if necessary.
tree_t t[NTREE];
tbb::parallel_for(0, (int)NTREE, [&t,&f,a0](int i){
t[i] = morpho::cToS(imtransform(f, [i](value_t x) { return x[i]; }), c4);
if (a0 > 0) {
grain_filter_inplace(t[i], a0);
t[i].shrink_to_fit();
}
});
// 2. Compute the Gos.
MyGraph g2;
std::array<property_map<tree_t, typename MyGraph::vertex_descriptor>, NTREE> tlink;
std::tie(g2, tlink) = compute_g2<NTREE>(t);
// 3. Compute the depth image
boost::vector_property_map<unsigned> gdepth = compute_graph_depth(g2);
image2d<uint16> imdepth = imchvalue<uint16>(F);
write_vmap_to_image(t, &tlink[0], gdepth, imdepth);
// debug
// io::imsave(imdepth, "depth.tiff");
// 4. Compute the saturated maxtree
std::tie(T, std::ignore) = satmaxtree(imdepth);
T = tree_keep_2F(T);
}
// Mumford-shah + simplification
{
// Grain filter
if (a1 > 0)
{
grain_filter_inplace(T, a1);
T.shrink_to_fit();
}
std::cout << "Simplification (number of llines): " << std::endl;
mumford_shah_on_tree(T, F, lambda);
F = Kadjust_to(F, ima.domain());
io::imsave(F, output_path);
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment