agondel | Bonjour,
Je travaille dans le cadre d'un stage de recherche sur la localisation des ondes sur une plaque de métal.
J'utilise le logiciel Freefem++ afin de résoudre mes équations aux DPP.
A partir de la fonction trouvée et de sa cartographie (voir la photo que je poste), je dois trouver les maxima (crêtes) de la fonction.
Comme il n'y a pas de méthode "watershed" sur Freefem++, je dois compiler un code C++.
Mon tuteur de stage m'a donné le code C++ que je devais utiliser pour réaliser la fonction "watershed".
Néanmoins, ce code comporte des erreurs, et étant débutante en C++, je n'arrive pas à les résoudre.
Pouvez-vous m'aider s'il-vous-plaît ?
Voici le code :
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
- 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);
|
Et voici le genre de fonction dont je veux trouver les crêtes :
 |
![[:nazztazz:2]](https://forum-images.hardware.fr/images/perso/2/nazztazz.gif "[:nazztazz:2]")