Commit 45c25ccf authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Implement energy on top the new Component Tree Framework.

        *  apps/CMakeLists.txt,
        *  apps/attributes/CMakeLists.txt,
        *  apps/attributes/cMSER.hpp,
        *  apps/attributes/cMeaningFullNess.hpp,
        *  apps/attributes/curvature.cpp,
        *  apps/attributes/curvature.hpp,
        *  apps/attributes/MSER-patternspec-cli.cpp,
        *  apps/attributes/meaningfullness-cli.cpp,
        *  apps/attributes/meaningfullness-patternspec-cli.cpp: New.

Old implementation maintained for bwd compatibility.
        *  apps/attributes/MSER.hpp,
        *  apps/attributes/meaningfullness.hpp,
        *  apps/attributes/meaningfullnessArgparser.cpp,
        *  apps/tests/attributes/meaningfullness.cpp,
parent 7c15ecec
......@@ -3,6 +3,7 @@ add_subdirectory(tests)
add_subdirectory(imview)
add_subdirectory(maxtree_comparison)
add_subdirectory(attributes)
add_subdirectory(simplification)
add_subdirectory(clattice)
add_subdirectory(tosgui)
......
link_libraries(-lfreeimage)
add_executable(MSER-patternspec-cli MSER-patternspec-cli.cpp)
add_executable(meaningfullness-patternspec-cli meaningfullness-patternspec-cli.cpp curvature.cpp)
add_executable(meaningfullness-cli meaningfullness-cli.cpp curvature.cpp)
\ No newline at end of file
#include <mln/core/image/image2d.hpp>
#include <mln/core/algorithm/accumulate.hpp>
#include <mln/accu/accumulators/max.hpp>
#include <mln/morpho/component_tree/io.hpp>
#include <mln/morpho/component_tree/accumulate.hpp>
#include <mln/morpho/component_tree/reconstruction.hpp>
#include <mln/morpho/component_tree/filtering.hpp>
#include <mln/morpho/component_tree/pattern_spectra.hpp>
#include <mln/morpho/algebraic_filter.hpp>
#include <mln/accu/accumulators/accu_if.hpp>
#include <mln/accu/accumulators/count.hpp>
#include <apps/tos/topology.hpp>
#include "cMSER.hpp"
#include <mln/morpho/extinction.hpp>
#include <mln/morpho/datastruct/checktree.hpp>
int main(int argc, char** argv)
{
if (argc < 9)
{
std::cerr << "Usage: " << argv[0] << "tree image.tiff λ ε γ t₁ t₂ spectra.csv\n"
<< " λ:\t Grain filter before anything else (number of 2F) (consider: 20-50, *20*) \n"
<< " ε:\t The delta level considered when fetching the neighborhood of a node in MSER (consider 5-20, *10*)\n"
<< " γ:\t Not used anymore.\n"
<< " t₁:\t Threshold above which node having a MSER value greater that t₁ are removed (consider 0.7-2.0, *1.0*)\n"
<< " t₂:\t Threshold below which node having an extinction value lesser than t₂ are removed (consider 0-t₁, *0.2*).\n";
std::exit(1);
}
using namespace mln;
const char* tree_path = argv[1];
const char* img_path = argv[2];
unsigned grain = std::atoi(argv[3]);
int eps = std::atoi(argv[4]);
int areaAS = std::atoi(argv[5]);
float threshold1 = std::atof(argv[6]);
float threshold2 = std::atof(argv[7]);
const char* out_path = argv[8];
typedef morpho::component_tree<unsigned, image2d<unsigned> > tree_t;
tree_t tree;
{
std::ifstream f(tree_path, std::ios::binary);
morpho::load(f, tree);
}
image2d<uint16> ima;
io::imread(img_path, ima);
if (ima.domain() != tree._get_data()->m_pmap.domain())
{
std::cerr << "Domain between image differ.\n"
<< ima.domain() << " vs "
<< tree._get_data()->m_pmap.domain() << std::endl;
}
auto vmap = morpho::make_attribute_map_from_image(tree, ima);
accu::accumulators::accu_if< accu::accumulators::count<>,
K1::is_face_2_t,
point2d > counter;
auto amap = morpho::paccumulate(tree, ima, counter);
auto pred = make_functional_property_map<tree_t::vertex_id_t>([&amap, grain](tree_t::vertex_id_t x) {
return amap[x] > grain;
});
morpho::filter_direct_inplace(tree, pred);
morpho::internal::checktree(tree);
auto amser = compute_MSER(tree, vmap, amap, eps, MSER_NORM);
// convert to image and filter
auto imser = morpho::make_image(tree, amser);
{
//float maxval = accumulate(imser, accu::features::max<> ());
mln_foreach(float& v, imser.values())
if (v > threshold1)
v = threshold1;
}
auto nmser = morpho::extinction(imser, morpho::tree_neighb_t());
//auto nmser = eval(nmser_ / imser);
//auto nmser = morpho::area_closing(imser, morpho::tree_neighb_t(), areaAS);
{
std::ofstream fout(out_path);
fout << "area,mser,extinction\n";
mln_foreach(auto x, tree.nodes())
fout << amap[x] << "," << amser[x] << "," << nmser(x) << "\n";
fout.close();
}
{
image2d<uint16> tmp;
resize(tmp, ima);
auto& attr = nmser.get_vmap();
auto pred = [&attr, threshold2](const tree_t::node_type& x) { return threshold2 < attr[x]; };
morpho::filter_direct_and_reconstruct(tree, make_functional_property_map<tree_t::node_type>(pred),
vmap, tmp);
io::imsave(tmp, "selection.tiff");
}
// image2d<float> out(200,200);
// morpho::pattern_spectra(tree, amser, amap,
// make_functional_property_map([&](const tree_t::node_type& x) -> float {
// return amap[x];
// }), out, true, true);
// io::imsave(out, argv[4]);
}
......@@ -78,15 +78,27 @@ compute_MSER_attribute(const mln::image2d<V>& f,
// compute area
image2d<unsigned> area;
image2d<unsigned> maxarea; // the max area of the children in Δh
{
resize(area, K).init(0);
resize(maxarea, K).init((unsigned)0);
area[S[0]] = 1;
for (int i = S.size()-1; i > 0; --i)
{
unsigned x = S[i];
++area[x];
area[parent[x]] += area[x];
if (K[x] != K[parent[x]]) { // canonical propagate to parent
unsigned y = parent[x];
unsigned a = area[x];
V v = f[x];
while (parent[y] != y and dist(v, f[y]) < eps) {
y = parent[y];
}
maxarea[y] = std::max(a, maxarea[y]);
}
}
}
......@@ -103,13 +115,19 @@ compute_MSER_attribute(const mln::image2d<V>& f,
V v = f[x];
unsigned y = parent[x];
while (dist(v, f[y]) < eps and y != parent[y])
unsigned smally = parent[x];
while (dist(v, f[y]) < eps and y != parent[y]) {
if (dist(v,f[y]) < (eps / 3))
smally = y;
y = parent[y];
}
switch (amser) {
case MSER_DIFF: mser[x] = area[y] - area[x]; break;
case MSER_RATIO: mser[x] = (float) area[x] / area[y]; break;
case MSER_NORM: mser[x] = (float) (area[y] - area[x]) / area[x]; break;
case MSER_DIFF: mser[x] = area[y] - maxarea[x]; break;
case MSER_RATIO: mser[x] = //(float) area[x] / area[y]; break;
//mser[x] = std::max(10.0f, (float) (area[y] - maxarea[x]) / area[x]); break;
mser[x] = (float) (area[smally] - area[x]) / (area[y] - area[x]); break;
case MSER_NORM: mser[x] = (float) (area[y] - maxarea[x]) / maxarea[x]; break;
}
}
......
#ifndef APPS_ATTRIBUTE_CMSER_HPP
# define APPS_ATTRIBUTE_CMSER_HPP
# include <mln/core/trace.hpp>
# include <mln/morpho/component_tree/component_tree.hpp>
enum eMSER_attribute {
MSER_DIFF, // Val > 0, δΓ = A(Γ⁺) - A(Γ⁻)
MSER_RATIO, // 0 < Val < 1 δΓ = A(Γ⁻) / A(Γ⁺)
MSER_NORM // 0 < Val < 1 δΓ = (A(Γ⁺) - A(Γ⁻)) / A(Γ)
};
enum eMSER_accum_type {
MSER_ABSOLUTE, // Γ⁺ = ⋁ { Γ ⊂ Γ', |V(Γ') - V(Γ)| ≥ ε }
MSER_SUM // Γ⁺ = ⋁ { Γ ⊂ Γ', ∑ |V(par(u)) - V(u)| > ε, u ∈ [Γ → Γ'] }
};
namespace mln
{
template <class P, class Amap,
class ValuePropertyMap,
class AreaPropertyMap,
class Distance>
property_map<morpho::component_tree<P, Amap>, float>
compute_MSER(const morpho::component_tree<P, Amap>& tree,
const ValuePropertyMap& vmap,
const AreaPropertyMap& amap,
float eps,
eMSER_attribute amser = MSER_DIFF,
eMSER_accum_type fsum = MSER_ABSOLUTE,
Distance dist = Distance());
template <class P, class Amap,
class ValuePropertyMap,
class AreaPropertyMap>
property_map<morpho::component_tree<P, Amap>, float>
compute_MSER(const morpho::component_tree<P, Amap>& tree,
const ValuePropertyMap& vmap,
const AreaPropertyMap& amap,
float eps,
eMSER_attribute amser = MSER_DIFF,
eMSER_accum_type fsum = MSER_ABSOLUTE);
/**********************/
/** Implementation **/
/**********************/
template <class P, class Amap,
class ValuePropertyMap,
class AreaPropertyMap,
class Distance>
property_map<morpho::component_tree<P, Amap>, float>
compute_MSER(const morpho::component_tree<P, Amap>& tree,
const ValuePropertyMap& vmap,
const AreaPropertyMap& amap,
float eps,
eMSER_attribute amser,
eMSER_accum_type fsum,
Distance dist)
{
mln_entering("mln::compute_MSER");
typedef morpho::component_tree<P, Amap> tree_t;
typedef typename tree_t::node_type node_t;
property_map<tree_t, unsigned> aplus(tree);
property_map<tree_t, unsigned> aminus(tree, 0);
if (fsum == MSER_ABSOLUTE)
{
mln_reverse_foreach(auto n, tree.nodes())
{
auto x = n;
while (x.get_parent_id() != tree.npos() and dist(vmap[n], vmap[x]) < eps)
x = x.parent();
aminus[x] = std::max<unsigned>(aminus[x], amap[n]);
aplus[n] = amap[x];
}
}
else
{
mln_reverse_foreach(auto n, tree.nodes())
{
auto x = n;
float d = 0;
while (x.id() != x.get_parent_id() and d < eps)
{
d += dist(vmap[x], vmap[x.parent()]);
x = x.parent();
}
aminus[x] = std::max<unsigned>(aminus[x], amap[n]);
aplus[n] = amap[x];
}
}
aplus[tree.get_root_id()] = amap[tree.get_root()];
std::function<float(const node_t&)> f;
switch (amser) {
case MSER_DIFF:
f = [&](const node_t& x) -> float { return aplus[x] - aminus[x]; };
break;
case MSER_RATIO:
f = [&](const node_t& x) -> float { return (float)amap[x] / aplus[x]; };
break;
case MSER_NORM:
f = [&](const node_t& x) -> float { return (float)(aplus[x] - aminus[x]) / amap[x]; };
break;
};
property_map<tree_t, float> out(tree);
mln_foreach(auto x, tree.nodes())
out[x] = f(x);
mln_exiting();
return out;
}
template <class P, class Amap,
class ValuePropertyMap,
class AreaPropertyMap>
property_map<morpho::component_tree<P, Amap>, float>
compute_MSER(const morpho::component_tree<P, Amap>& tree,
const ValuePropertyMap& vmap,
const AreaPropertyMap& amap,
float eps,
eMSER_attribute amser,
eMSER_accum_type fsum)
{
typedef typename ValuePropertyMap::value_type V;
auto dist = [](V x, V y) -> float { return l2norm(x - y); };
return compute_MSER(tree, vmap, amap, eps, amser, fsum, dist);
}
}
#endif // ! APPS_ATTRIBUTE_CMSER_HPP
#ifndef APPS_ATTRIBUTE_CMEANINGFULLNESS_HPP
# define APPS_ATTRIBUTE_CMEANINGFULLNESS_HPP
# include <mln/core/image/image2d.hpp>
# include <mln/core/neighborhood/dyn_neighborhood.hpp>
# include <mln/core/trace.hpp>
# include <mln/accu/accumulator.hpp>
# include <mln/accu/accumulators/variance.hpp>
# include <mln/accu/accumulators/mean.hpp>
# include <mln/morpho/component_tree/component_tree.hpp>
# include <mln/morpho/component_tree/compute_depth.hpp>
# include <apps/tos/topology.hpp>
# include <apps/tos/croutines.hpp>
# include <apps/g2/accu/lca.hpp>
# include "curvature.hpp"
namespace mln
{
/// \brief compute the external regional energy
///
template <class P, class AMap, class I>
property_map<morpho::component_tree<P, AMap>, float>
compute_meaningfullness_external_energy(const morpho::component_tree<P, AMap>& tree,
const Image<I>& vmap,
int eps);
/// \brief generic function to compute a contextual energy
/// This function supposes a tree computed with 0 and 1-faces. The accumulation
/// is performed on 2F only.
/// Considering a ball Β of radius ε, the accumulator is computed for every shape Γ on:
/// * Ext(Γ) = δ(Γ) \ Γ
/// * Int(Γ) = Γ \ ε(Γ)
/// * Reg(Γ) = Ext(Γ) ⋃ Int(Γ)
///
/// \param tree The component tree.
/// \param valuemap Image with values to accumulate (only 2F are used).
/// \param accu The accumulator (-like) to compute.
/// \param eps The radius of the ball used as the SE for shape dilation and erosion.
///
/// \return An attribute map for each node with the 3 values of the
/// accumulation on thes internal/external/whole regions.
template <class P, class I, class AccuLike>
property_map<morpho::component_tree<P, image2d<P> >,
vec<typename accu::result_of<AccuLike, mln_value(I)>::type, 3> >
compute_regional_energy(const morpho::component_tree<P, image2d<P> >& tree,
const Image<I>& valuemap,
const AccumulatorLike<AccuLike>& accu,
int eps);
/*********************/
/** Implementation ***/
/*********************/
template <class P, class I, class AccuLike>
property_map<morpho::component_tree<P, image2d<P> >,
vec<typename accu::result_of<AccuLike, mln_value(I)>::type, 3> >
compute_regional_energy(const morpho::component_tree<P, image2d<P> >& tree,
const Image<I>& valuemap,
const AccumulatorLike<AccuLike>& accu_,
int eps)
{
mln_entering("mln::compute_regional_energy");
typedef morpho::component_tree<P, image2d<P> > tree_t;
typedef typename tree_t::node_type node_type;
typedef typename tree_t::vertex_id_t vertex_id_t;
typedef typename accu::result_of<AccuLike, mln_value(I)>::type R;
const I& vmap = exact(valuemap);
auto acc = accu::make_accumulator(exact(accu_), mln_value(I) ());
typedef decltype(acc) Accu;
// Check that domain matches.
mln_precondition(tree._get_data()->m_pmap.domain() == vmap.domain());
acc.init();
property_map<tree_t, Accu> external(tree, acc);
property_map<tree_t, Accu> internal(tree, acc);
// Create the se (ball only on 2F)
std::vector<point2d> dpoints;
for (int i = -eps; i < eps+1; ++i)
for (int j = -eps; j < eps+1; ++j)
if ( (i != 0 or j != 0) and (i*i + j*j) <= (eps*eps) )
dpoints.emplace_back(2*i,2*j);
dyn_neighborhood<std::vector<point2d>,
dynamic_neighborhood_tag> ball(dpoints);
// Compute the depth image.
property_map<tree_t, unsigned> depth = morpho::compute_depth(tree);
// Compute internal region attribute
{
mln_pixter(px, vmap);
mln_iter(nx, ball(px));
accu::least_common_ancestor<P, image2d<P>> lca(tree, depth);
mln_forall(px)
{
point2d p = px->point();
if (K1::is_face_2(p))
{
lca.init();
lca.take(px->index());
mln_forall(nx)
if (vmap.domain().has(nx->point()))
lca.take(nx->index());
else
break; // border pixel LCA is above the root
node_type x = tree.get_node_at(px->index());
node_type y = nx.finished() ? lca.to_result() : tree.nend();
for (;x != y; x = x.parent())
internal[x].take(px->val());
}
}
}
// Compute external region attribute
{
property_map<tree_t, bool> dejavu(tree, false);
std::vector< std::pair<vertex_id_t, vertex_id_t> > branches;
mln_pixter(px, vmap);
mln_iter(nx, ball(px));
mln_forall(px)
{
point2d p = px->point();
if (K1::is_face_2(p))
{
node_type k = tree.get_node_at(px->index());
branches.clear();
mln_forall(nx)
if (vmap.domain().has(nx->point()))
{
node_type x = tree.get_node_at(nx->index());
node_type y = lca(tree, depth, k, x);
branches.emplace_back(x.id(), y.id());
for (node_type n = x; n != y; n = n.parent())
if (not dejavu[n]) {
external[n].take(px->val());
dejavu[n] = true;
}
}
// undo dejavu
{
vertex_id_t x,y;
mln_foreach(std::tie(x,y), branches)
for (node_type n = tree.get_node(x); n.id() != y; n = n.parent())
dejavu[n] = false;
}
}
}
}
// Return the results.
property_map<tree_t, vec<R, 3> > res(tree);
{
mln_foreach(auto x, tree.nodes())
{
res[x][0] = internal[x].to_result();
res[x][1] = external[x].to_result();
internal[x].take(external[x]);
res[x][2] = internal[x].to_result();
}
}
mln_exiting();
return res;
}
template <class P, class AMap, class I>
property_map<morpho::component_tree<P, AMap>, float>
compute_meaningfullness_external_energy(const morpho::component_tree<P, AMap>& tree,
const Image<I>& vmap,
int eps)
{
auto attrmap = compute_regional_energy(tree, vmap,
accu::features::variance<double> (),
eps);
property_map<morpho::component_tree<P, AMap>, float> energy(tree);
mln_foreach(auto node, tree.nodes())
energy[node] = (attrmap[node][0] + attrmap[node][1]) / attrmap[node][2];
return energy;
}
}
#endif // ! APPS_ATTRIBUTE_CMEANINGFULLNESS_HPP
#include "curvature.hpp"
namespace mln
{
template
image2d<float> compute_curvature<uint8>(const image2d<uint8>&);
template
image2d<float> compute_curvature<uint16>(const image2d<uint16>&);
}
#ifndef APPS_ATTRIBUTE_CURVATURE_HPP
# define APPS_ATTRIBUTE_CURVATURE_HPP
# include <mln/core/image/image2d.hpp>
# include <mln/core/grays.hpp>
# include <mln/core/trace.hpp>
namespace mln
{
/// \brief Add the curvature on a the 1F.
/// The domain of the output image is twice as big as the original one.
template <class V>
image2d<float>
compute_curvature(const image2d<V>& ima);
extern template image2d<float> compute_curvature<uint8>(const image2d<uint8>&);
extern template image2d<float> compute_curvature<uint16>(const image2d<uint16>&);
/***************************/
/**** Implementation ***/
/***************************/
template <class V>
image2d<float>
compute_curvature(const image2d<V>& ima)
{
static_assert(std::is_integral<V>::value,
"V must be an integral type.");
trace::entering("mln::compute_curvature");
auto domain = ima.domain();
domain.pmin *= 2;
domain.pmax = domain.pmax * 2 - 1;
image2d<float> curv(domain, ima.border(), 0);
mln_foreach(const point2d& p, ima.domain())
{
// Right edge
{
float ux = ima.at(p + point2d{0,1}) - ima.at(p);
float uy =
(ima.at(p + point2d{1,0}) - ima.at(p + point2d{-1,0}) +
ima.at(p + point2d{1,1}) - ima.at(p + point2d{-1,1})) / 4.0;
float uxx =
(ima.at(p + point2d{0,-1}) - ima.at(p) -
ima.at(p + point2d{0,1}) + ima.at(p + point2d{0,2})) / 2.0;
float uxy =
(-ima.at(p + point2d{1,0}) + ima.at(p + point2d{-1,0}) +
ima.at(p + point2d{1,1}) - ima.at(p + point2d{-1,1})) / 2.0;
float uyy =
(ima.at(p + point2d{-1,0}) + ima.at(p + point2d{-1,1}) +