watershed problem

[agondel]

Hello

I'm currently working on the localization of elastic waves for a research semester at a university.
I use the software Freefem++ in order to solve equations.
Then, I have to find the crests of the solution found thanks to his graph (for example this kind of graph: http://www.casimages.com/i/141022093627498097.png).

Nevertheless, as there's not "watershed" method on Freefem++, I have to compile a C++ code.

My boss gave me the code below, but the latter doesn't work. As I don't know how a C++ code works, it's very hard for me to understand what I have to change.

Do you have an idea ?

Code :

#include "ff++.hpp"
// #ifndef WITH_NO_INIT
// #include "ff++.hpp"
// #include "AFunction_ext.hpp"
// #endif
// using namespace std;
#include <set>
#include <vector>
#include <map>
#include <algorithm>
//#include "msh3.hpp"
// #include <iostream>
using namespace  Fem2D;
// FreeFem glue
class WATERSHED_P1_Op : public E_F0mps
{
public:
    Expression eTh,eff,eret;
   
    static const int n_name_param = 1;
    static basicAC_F0::name_and_type name_param[n_name_param];
    Expression nargs[n_name_param];
public:
    WATERSHED_P1_Op(const basicAC_F0 &  args,Expression tth, Expression fff,Expression rrr)
    : eTh(tth),eff(fff),eret(rrr)
    {
        args.SetNameParam(n_name_param,name_param,nargs);
    }
    AnyType operator()(Stack stack) const;
   
private:
    template<typename T>
    T arg(int i, Stack stack, T a) const {
        return nargs[i]
        ? GetAny< T >( (*nargs[i])(stack) )
        : a;
    }
};
basicAC_F0::name_and_type WATERSHED_P1_Op::name_param[]= {
    {  "eps",  &typeid(double)}
};
// algorithm
typedef int triangle_t;
typedef int vertex_t;
typedef int color_t;
struct fat_vertex_t {
    vertex_t vertex;
    triangle_t triangle;
    int edge;
   
    fat_vertex_t(vertex_t v, triangle_t t, int e)
    : vertex(v), triangle(t), edge(e) {}
    friend bool operator<(fat_vertex_t const& a, fat_vertex_t const& b)
    { return a.vertex < b.vertex; }
   
    friend bool operator==(fat_vertex_t const& a, fat_vertex_t const& b)
    { return a.vertex == b.vertex; }
};
typedef std::vector<fat_vertex_t> vertices_t;
typedef std::pair<fat_vertex_t, double> ver_val_t;
struct cmp_t {
    bool operator()(ver_val_t const& t1, ver_val_t const& t2) const {
        return t1.second < t2.second;
    }
};
typedef std::priority_queue<ver_val_t, std::vector<ver_val_t>, cmp_t> queue_t;
typedef KNM<long> ret_type;
template<typename Func>
void for_each_triangle(Mesh const& Th, triangle_t const triangle0, int const edge0, Func func) {
    int const vertex = Th(triangle0, edge0);
   
    if( !func( triangle0 ) )
      return;
    int edge = edge0;
    int triangle = triangle0;
    for(;;) {
        edge = (edge + 1) % 3;
        if( Th(triangle, edge) == vertex )
            edge = (edge + 1) % 3;
        triangle = Th.ElementAdj( triangle, edge );
        if( triangle == triangle0 )
            return;
       
        if( triangle < 0 )
            break;
       
        if( !func( triangle ) )
            return;
    }
    triangle = triangle0;
    edge = edge0;
    for(;;) {
        edge = (edge - 1) % 3;
        if( Th(triangle, edge) == vertex )
            edge = (edge - 1) % 3;
        triangle = Th.ElementAdj( triangle, edge );
        if( triangle == triangle0 )
            return;
       
        if( triangle < 0 )
            break;
       
        if( !func( triangle ) )
            return;
    }
}
template<typename Func>
struct for_each_neighbor_helper {
    Func func;
    Mesh const& Th;
    bool operator()(triangle_t triangle) {
        for(int e = 0; e < 3; ++e)
            if(! func( Th(triangle, e), triangle, e ) )
                return false;
        return true;
    }
};
template<typename Func>
void for_each_neighbor(Mesh const& Th, triangle_t const triangle0, int const edge0, Func func) {
    for_each_neighbor_helper<Func> help = { func, Th };
    // check adjacent triangles
    for_each_triangle(Th, triangle0, edge0, help);
}
template<typename Cont>
void erase_unique(Cont& cont) {
    std::sort(cont.begin(), cont.end());
    cont.erase(
        std::unique(cont.begin(), cont.end()),
        cont.end()
    );
}
struct maxima_helper {
    KN<double> const& tff;
    double& maxval;
    bool& is_max;
    bool operator()(vertex_t vertex, triangle_t triangle, int edge) const {
        double val = tff[ vertex ];
        if(val > maxval) {
            is_max = false;
            return false;
        }
        return true;
    }
};
static void maxima(Mesh const& Th, KN<double> const& tff, vertices_t& vertices, double epsr)
{
    const int nbt=Th.nt; // nombre de triangles
    // loop over vertices
    for(int it = 0; it < nbt; ++it) {
        int maxiv = 0;
        double maxval = tff[ Th(it,0) ];
        int iv;
        for(iv=1; iv < 3; ++iv) {
            int i = Th(it,iv);
            double val = tff[i];
            if(val > maxval) {
                maxiv = iv;
                maxval = val;
            }
        }
        iv = maxiv;
       
        if(std::abs(maxval) < epsr)
            continue;
        bool is_max = true;
       
        maxima_helper helper = { tff, maxval, is_max };
       
        for_each_neighbor(Th, it, iv, helper);
        if(!is_max)
            continue;
//         std::cout << "FOUND " << it << ' ' << maxiv << ' ' << Th(it, maxiv) << ' ' << maxval << std::endl;
        vertices.push_back(fat_vertex_t( Th(it,maxiv), it, maxiv ));
    }
    erase_unique(vertices);
}
#if 0
static void maxima(Mesh const& Th, KN<double> const& tff, queue_t& roots, std::vector<color_t>& colors, double epsr)
{
    const int nbt=Th.nt; // nombre de triangles
    const int nbv=Th.nv; // nombre de vertices
    enum pixel_type {
        MAXIMUM,
        PLATEAU,
        NON_MAXIMUM
    };
    // the one that increments current_color
    // shall push to roots
    color_t current_color = 1;
    std::vector<bool> visited ( nbv, false );
    auto analyse_neighbors = [&](vertex_t const vertex0, triangle_t const triangle0, int edge0) {
        pixel_type pxl = MAXIMUM;
        for_each_neighbor(Th, triangle0, edge0,
          [&](vertex_t vertex, triangle_t triangle, int edge) {
            if( vertex == vertex0 )
                return true;
            if( tff[vertex] >  tff[vertex0] ) {
                pxl = NON_MAXIMUM;
                return false;
            }
            if( tff[vertex] == tff[vertex0] )
                pxl = PLATEAU;
            return true;
        });
        return pxl;
    };
    auto analyse_plateau = [&](vertex_t const vertex0, triangle_t const triangle0, int edge0) {
        colors[vertex0] = current_color;
        // early exit
        color_t new_label = current_color;
        // do not forget marked nodes
        std::deque<fat_vertex_t> queue;
        queue.push_back({ vertex0, triangle0, edge0 });
        auto it = queue.begin();
        auto const end = queue.end();
        for(; it != end; ++it ) {
            fat_vertex_t const& vv = *it;
            for_each_neighbor(Th, vv.triangle, vv.edge,
              [&](vertex_t vertex, triangle_t triangle, int edge) {
                if( colors[vertex] == -1 && tff[vertex] == tff[vertex0] ) {
                    colors[vertex] = current_color;
                    queue.push_back({ vertex, triangle, edge });
                    visited[vertex] = true;
                }
                else if( tff[vertex] > tff[vertex0] )
                    new_label = -1;
                return true;
            });
        }
        if( new_label == -1 )
            for(fat_vertex_t const& vv : queue)
                colors[vv.vertex] = -1;
        else {
            ++current_color;
            roots.push({ { vertex0, triangle0, edge0 }, tff[vertex0] });
        }
    };
    // loop over vertices
    for(triangle_t triangle = 0; triangle < nbt; ++triangle)
    for(int edge = 0; edge < 3; ++edge) {
        vertex_t vertex = Th( triangle, edge );
        if( visited[vertex] )
            continue;
        pixel_type pxl = analyse_neighbors(vertex, triangle, edge);
        if( pxl == MAXIMUM ) {
            for_each_neighbor(Th, triangle, edge,
              [&](vertex_t vertex2, int,int) {
                ffassert( tff[vertex2] <= tff[vertex] );
                return true;
            });
            colors[vertex] = current_color++;
            roots.push({{ vertex, triangle, edge }, tff[vertex] });
        }
//         else if( pxl == PLATEAU )
//             analyse_plateau(vertex, triangle, edge);
        visited[vertex] = true;
    }
    ffassert( roots.size() == current_color-1 );
}
#endif 

[agondel]

 struct color_one_neighbor {
    KN<double> const& tff;
    fat_vertex_t const& current;
    std::vector<color_t>& colors;
    color_t const current_color;
    queue_t& queue;
    bool operator()(vertex_t vertex, triangle_t triangle, int edge) {
       
        fat_vertex_t vv ( vertex, triangle, edge );
        if(vertex == current.vertex)
            return true;
        color_t& color = colors[vertex];
        if( color == -1 ) {
            color = current_color;
            queue.push(ver_val_t( vv, tff[vertex] ));
        }
        else if( color != current_color ) {
            color = 0;
            frontier.push_back( vv );
//                 ffassert( tff[vertex] <= tff[current.vertex] ); // TODO ça explose ici
//                 std::cout << "FOUND " << vertex << " -> " << color << std::endl;
        }
        return true;
    }
};
AnyType WATERSHED_P1_Op::operator()(Stack stack) const
{
    MeshPoint *mp(MeshPointStack(stack));
    ret_type& ret = *GetAny<ret_type* >( (*eret)(stack) );
    Mesh* pTh = GetAny<Mesh *>( (*eTh)(stack) );
   
    ffassert(pTh);
    double  epsr = arg(0,stack,1e-5);
    Mesh const& Th = *pTh;
    const int nbv=Th.nv; // nombre de sommet
    const int nbt=Th.nt; // nombre de triangles
    const int nbe=Th.neb; // nombre d'aretes fontiere
    const double unset = -1e-100;
    KN<double> tff(nbv, unset);
    // loop over triangle
    for(int it=0; it < nbt; ++it) {
        for(int iv=0; iv<3; ++iv) {
            int i = Th(it,iv);
            if(tff[i]==unset) {
                mp->setP(pTh,it,iv);
                tff[i]=GetAny<double>((*eff)(stack));
            }
        }
    }
    queue_t queue;
    std::vector<color_t> colors ( nbv, -1 );
    vertices_t frontier;
   
    // prefill
    {
        vertices_t roots;
        maxima(Th, tff, roots, epsr);
        color_t color = 1;
       
        vertices_t::iterator it = roots.begin(), en = roots.end();
        for(; it != en; ++it) {
            fat_vertex_t const& current = *it;
            colors[current.vertex] = color++;
            queue.push(ver_val_t( current, tff[current.vertex] ));
        }
    }
    // loop
    while( !queue.empty() ) {
        fat_vertex_t const current = queue.top().first; queue.pop();
        color_t const current_color = colors[current.vertex];
        ffassert( current_color != -1 );
        if( current_color == 0 )
            continue;
        // check adjacent triangles
        for_each_neighbor(
            Th, current.triangle, current.edge,
            color_one_neighbor(tff, current, colors, current_color, queue)
        );
    }
    erase_unique(frontier);
    std::cout << "OUT " << frontier.size() << std::endl;
    ret.resize(2, frontier.size());
    for(int k = 0; k < frontier.size(); ++k) {
        fat_vertex_t const& vv = frontier[k];
        ret(0, k) = vv.triangle;
        ret(1, k) = vv.edge;
    }
    return 0l;
}
class  WATERSHED_P1: public OneOperator { public: 
    typedef Mesh *pmesh;
    typedef std::pair<FEbase<double, v_fes>*, int> fem_t;
    WATERSHED_P1() : OneOperator(atype<long>(),atype<pmesh>(),atype<double>(), atype<ret_type*>() ) {}
   
    E_F0 * code(const basicAC_F0 & args) const
    {
        return  new WATERSHED_P1_Op( args,
                                  t[0]->CastTo(args[0]),
                                  t[1]->CastTo(args[1]),
                                  t[2]->CastTo(args[2]) );
    }
};
void finit()
{
    Global.Add("watershed","(",new WATERSHED_P1);
}
LOADFUNC(finit);

[agondel]

And the mistakes are :

watershed.cxx03.cpp: In member function ‘bool color_one_neighbor::operator()(vertex_t, triangle_t, int)’:
watershed.cxx03.cpp:353:13: error: ‘frontier’ was not declared in this scope
frontier.push_back( vv );
^
watershed.cxx03.cpp: In member function ‘virtual AnyType WATERSHED_P1_Op::operator()(Stack) const’:
watershed.cxx03.cpp:427:74: error: no matching function for call to ‘color_one_neighbor::color_one_neighbor(KN<double>&, const fat_vertex_t&, std::vector<int>&, const color_t&, queue_t& )’
color_one_neighbor(tff, current, colors, current_color, queue)
^
watershed.cxx03.cpp:427:74: note: candidates are:
watershed.cxx03.cpp:330:8: note: color_one_neighbor::color_one_neighbor()
struct color_one_neighbor {
^
watershed.cxx03.cpp:330:8: note: candidate expects 0 arguments, 5 provided
watershed.cxx03.cpp:330:8: note: color_one_neighbor::color_one_neighbor(const color_one_neighbor& )
watershed.cxx03.cpp:330:8: note: candidate expects 1 argument, 5 provided

[norm b]

error: ‘frontier’ was not declared in this scope

You need to pass "frontier" into the struct color_one_neighbor.

As for the other errors, take a look at how maxima_helper is coded and implemented.

I have never used Freefem++ but try the following:

In color_one_neighbor:

struct color_one_neighbor {
    KN<double> const& tff;
    fat_vertex_t const& current;
    std::vector<color_t>& colors;
    color_t const current_color;
    queue_t& queue;
    vertices_t& frontier;  // added line
    bool operator()(vertex_t vertex, triangle_t triangle, int edge) {
    . . .

In WATERSHED_P1_Op::operator():

    . . .
    while( !queue.empty() ) {   // line 68 above
        fat_vertex_t const current = queue.top().first; queue.pop();
        color_t const current_color = colors[current.vertex];
        ffassert( current_color != -1 );
        if( current_color == 0 )
            continue;
        // check adjacent triangles
        color_one_neighbor neighbor = { ttf, current, colors, current_color, queue, frontier }; // added
        for_each_neighbor(Th, current.triangle, current.edge, neighbor); // changed
    }
    . . .

Caveat utilitor!

HTH

[agondel]

hello,

thank you for you help !

normally everything works now !