Discussion :

[vince67]

Bonjour à tous,
je suis étudiant et dans le cadre d'un projet en MNT (modèle numérique de terrain), je dois sélectionner des points sur un échantillon disposé en grille.
En fait, je pars de points répartis en grille, donc les x et y sont connus. Ils ont aussi une valeur d'altitude. Seulement, cet échantillon est trop important. Je dois donc rédiger un algorithme permettant d'éliminer les points les moins significatifs.
Cet algorithme fonctionne sous le principe suivant:
- on supprime un point,
- on génère une triangulation de Delaunay sur tout l'échantillon, sans ce point.
- A partir de la triangulation, on interpole (linéairement) l'altitude de ce point.
- Puis on calcule la différence d'altitude entre ce point et le point interpolé.
- On le remet alors dans l'échantillon et on fait la même chose à tous les autres points de l'échantillon.
-Une fois effectué, on regarde quel point a la plus petite différence et on le supprime de l'échantillon.

Et le truc c'est qu'alors on prend cet échantillon modifié, et on recommence le processus (génération de triangulation de delaunay etc).

J'ai plusieurs questions:
- je pensais partir avec une matrice, mais comme je vais devoir supprimer des données, comment faire une matrice à trous?
- j'ai essayé la fonction de triangulation de delaunay dans Matlab mais je ne comprends pas la signification de la matrice "Tri".
En gros, comment choisir le triangle qu'il faut pour interpoler l'altitude du point.

Merci d'avance de vos réponses

Vincent

[Jean Dumoncel]

Bonjour,

Ta méthode devrait fonctionner mais tu n'as pas de contraintes sur les temps de calcul? Car ça risque d'être un peu long toutes ces triangulations...

Ne pourrais-tu passer par un fit moyen de la surface et étudier les résidus pour détecter les points les plus éloigné de ce fit?

Si tu cherches une méthode plus optimisé, tu peux peut-être poser la question dans le forum Algorithmes.

- je pensais partir avec une matrice, mais comme je vais devoir supprimer des données, comment faire une matrice à trous?

Tout d'abord, pour l'interpolation regarde la fonctions interp2.
Je mettrais des valeurs NaN pour les points supprimés. Tu peux essayer ces contributions du fex, elles sont rapport avec ce que tu souhaites faire :
inpaintnans
surface fitting using gridfit

- j'ai essayé la fonction de triangulation de delaunay dans Matlab mais je ne comprends pas la signification de la matrice "Tri".
En gros, comment choisir le triangle qu'il faut pour interpoler l'altitude du point.

Avec la fonction Delaunay, Tri contient les triplets de points composant chaque triangle : chaque valeur correspond aux indices dans les grilles x et y.

[vince67]

Hello,
merci pour tes réponses!
A priori, je n'ai pas vraiment de contrainte de temps, il faut juste que ça marche!
j'ai effectivement réussi à trouver les sommets dans "tri". Mais la question que je me pose, c'est comment savoir dans quel triangle se trouve le point que j'ai supprimé? une fois que la nouvelle triangulation de delaunay a été faite?je fais "simplement" une comparaison avec des boucles if pour savoir si les coordonnées sont à l'intérieur?
Et pour générer cette nouvelle triangulation, il "suffit" que je mette NaN aux coordonnées du sommet supprimé?

[Jean Dumoncel]

Perso, voici ce que je ferais :
je conserverais :

• d'un côté les coordonnées de tes points sous forme de grille dans laquelle je mettrai à NaN les points qui seront exclus. Cette grille sera utilisé pour l'interpolation. (cela permet de garder sous forme de grid ,ce qui est utile pour pouvoir interpoler les données)
• et une liste des coordonnées des points dans laquelle seront supprimés les points exclus. Cette liste sera utilisé pour le calcul de la triangulation.

Essaie déjà de faire un code sur un échantillon de tes données.

[Jerome Briot]

A la place des NaN, tu peux aussi remplacer les coordonnées du point à supprimer par un point quelconque situé en dehors de ton maillage ou par un point du maillage dont tu sais a priori que tu ne le supprimeras pas. Tu "empiles" donc en quelque sorte successivement ces points vers une localisation "poubelle".

Ensuite pour savoir dans quel triangle se trouve le point, commence par regarder ton maillage depuis la direction z... tu peux donc simplifier ton problème en projettant ton maillage sur le plan (xoy) en ne prenant que les deux premières composantes des coordonnées (x,y) pour le calcul en 2D

Sinon, quelle est ta version de MATLAB ?
Les fonctions de triangulation ont pas mal variée dans les dernières version de MATLAB (avec l'introduction de la bibliothèque CGAL entre autre)

[ol9245]

Ca devrait être une base pour ton algo (vérifier si il faut pas gérer les cas frontière, etc...).

Pour pouvoir itérer facilement, je déplace le point à supprimer au fond du tableau des XY (et je modifie la triangulation en conséquence). De cette manière, l'algo revient à itérer la suppression du dernier point.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
close all
clc
clear
N = 100 ;
XY = rand(N, 2) ;
T1 = delaunay(XY(:,1), XY(:,2)) ;
 
% on enlève systématiquement l'élément n° 10
% ici plugger ton code qui sélectionne le meilleur point à supprimer
d = 10 ;
 
while N>50
    hold off
    triplot(T1, XY(1:N,1), XY(1:N,2), 'k')
    hold on
    % mettre cet élément à la fin de la triangulation
    r = [1:d-1, d+1:N, d] ;
    XY = XY(r,:) ;
    D = T1==d ;
    P = T1>d ;
    T1(P) = T1(P) - 1 ;
    T1(D) = N ;
    triplot(T1, XY(:,1), XY(:,2), 'b')
    plot(XY(end,1), XY(end,2), 'ob')
 
    % On est ramené à supprimer le dernier élément
 
    % refaire la triangulation
    % Note : il existe un théorème qui dit que les modifications de la triangulation ne peuvent être que locales
    % donc en pratique, il est inutile de refaire toute la triangulation. il sufit de virer les triangles locaux, puis de "réparer" en reconstruisant
    % une triangulation à l'intérieure de la zone nettoyée.
    T2 = delaunay(XY(1:end-1,1), XY(1:end-1,2)) ;
    triplot(T2, XY(1:end-1,1), XY(1:end-1,2), 'r') % OK. Les deux triangulations se superposent
 
    % chercher le point P dans la nouvelle triangulation
    K = dsearch(XY(1:end-1,1), XY(1:end-1,2), T2, XY(end,1), XY(end,2)) ;
    Tcible = find(T2(:,1)==K | T2(:,2)==K | T2(:,3)==K) ;
    for t = 1:numel(Tcible)
        tt = Tcible(t) ;
        s = T2(tt,:).' ; % les sommets du triangle potentiellement vainqueur
        p = patch(XY(s,1), XY(s,2), [0,1,0]) ;
        set(p, 'FaceAlpha', 0.2, 'Edgecolor', 'none')
        if inpolygon(XY(end,1), XY(end,2), XY(s,1), XY(s,2))
            fprintf(1, 'triangle n°%d défini par les points n°%d, %d, %d\n', tt, s(1), s(2), s(3))
            set(p, 'FaceAlpha', 1) ;
        end
    end
 
    % comme ton point à jeter est le dernier, tu peux boucler l'algo sur les N-1 points restants
    N = N-1 ;
    T1 = T2 ;
end

Voila voila.

Sinon, je note quand même une contradiction dans ton message d'origine. D'un côté tu dis que tu as tellement de puissance de calcul que tu te moques d'optimiser ton code, d'un autre coté ton problème (réduire le nombre de points sans réduire la quantité d'information) est typiquement motivé par un problème de puissance de calcul. hum... pas bien logique tout ça....

deuxième point : Je suppose que si tu travailles sur cet algo, c'est parceque tu vas avoir à traiter des MNT très lourds, genre LIDAR. Je partage donc à &00% la perplexité des camarades sur le gaspillage de temps de calcul de ton algo. Construire une triangulation est en N.LogN, Pour supprimer un point, tu va construire N triangulations et, finalement, le nombre de points à supprimer est proportionnel à N. Ton algo est donc en N^3.Log(N). Les algos en N^3 sont, par nature et sans même regarder le code, inutilisables sur de grands échantillons. Avant de commencer à coder, je te suggère donc de trouver un algo qui soit en N.Log(N) : les seuls qui fonctionnent sur des très grands échantillons.

Encore une chose : Matlab a une légère tendance à se scratcher quand on lui demande une triangulation sur un grand nombre (millions) de points.

Et finalement : tu as fait quoi comme biblio sur le traitement des MNT LIDAR avant de te mettre à coder ? je sais, je me mêle de ce qui me regarde pas....

en tous cas, le code que je t'ai fait devrait t'aider à répondre à la question posée. Pour le reste .... reviens nous voir quand tu veux



EDIT : une suggestion cosmétique qui transforme ton algo en NLogN :
1) prendre 1000 points au hasard, les supprimer, calculer l'erreur introduite. Ces 1000 erreurs te permettent de déterminer la loi statistique.
2) A partir de cette loi, tu définis un seuil de tolérance : toute erreur inférieure à ce seuil est tolérée et le point sera supprimé.
3) Pour chaque point un à un : le supprimer -> la triangulation T1 se transforme en T2, calculer l'erreur introduite, si acceptable : T1 = T2, sinon, continuer. -> tu gagnes un facteur N , tu n'es plus que en N²LogN.

En ne travaillant les triangulations que sur des voisinages locaux (voir commentaires dans le code sur les théorèmes existants), tu regagnes un autre facteur N et là, tu as enfin un algo semi linéaire (NLogN).

[vince67]

Bonjour Bonjour!
je reviens de week end où je n'étais pas chez moi, et n'avais pas internet.
Merci pour les réponses, ca devrait effectivement m'aider! j'ai la version matlab mac 7.6 (release 2008)
En fait, c'est un excercice de TP, donc dans le cas de l'exo on a juste à garder 15 points sur une matrice 7x7. Donc c'est evident que l'optimisation n'a aucune importance pour l'instant. Mais je vais tout de même regarder et essayer de faire un truc propre.
Je reviens vers vous si j'ai de nouveaux problèmes (je me fais pas d'illusions, ca va surement arriver..)

EDIT:
Bon je me suis mis dedans la et j'ai réussi à générer mes triangulations de Delaunay en supprimant alternativement un point.
J'ai aussi programmé la partie interpolation. c'est assez simple parce que je ne m'occupe pas des bords!
Maintenant, l'étape ultime pour joindre les deux: comment savoir à partir de quel triangle je dois interpoler le point que j'ai supprimé? En gros, comment à partir du résultat "TRI" trouver le le triangle dans lequel est inclus le point que je viens de supprimer?

[vince67]

Bonjour à tous de nouveau!
Pour localiser les triangles, j'ai utilisé la fonction "tsearch". Cependant elle n'existe plus dans les nuvelles versions de matlab (à partir de la version 8 je crois), à cause de bugs parfois.
Quelqu'un pourrait-il me donner son avis sur le code? ce serait ultra sympa!

Merci d'avance

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
 
function[Grille_Lab3mod]=delaunay(Grille_Lab2,dataset)
% intitialisation des données
  x2 = dataset(:,1); 
  y2 = dataset(:,2);
  x3=zeros(48,1);
  y3=zeros(48,1);
  x4=zeros(48,1);
  y4=zeros(48,1);
  datasetz=[Grille_Lab2(:,1);Grille_Lab2(:,2);Grille_Lab2(:,3);Grille_Lab2(:,4);Grille_Lab2(:,5);Grille_Lab2(:,6);Grille_Lab2(:,7)];
 
  for k=1:48
 
        for i=1:48 % initialisation de la matrice x3 et y3
                 x3(i,1)=x2(i,1);
                 y3(i,1)=y2(i,1);
                 x4(i,1)=x2(i,1);
                 y4(i,1)=y2(i,1);
                 z2(i,1)=datasetz(i,1);
        end
% coordonnées de l'élément supprimé
  xsupp=x2(k,1)
  ysupp=y2(k,1)
  zsupp=datasetz(k,1)
 
 
        for i =k:48 %48 on enlève à chaque fois un élément de la matrice x3, y3 depuis la matrice x2,y2
                 x3(i,1)=x2(i+1,1);
                 y3(i,1)=y2(i+1,1);
                 z2(i,1)=datasetz(i+1,1);
        end
 
tri = delaunay(x3,y3); % on exécute la triangulation de Delaunay avec le point supprimé
 
        % figure;
        %triplot(tri,x3,y3);
        %trimesh(tri,x3,y3,z2);
        %hold all;
 
 
    % On s'apprete ici à  interpoler l'altitude du point surpprimé
 
        if  xsupp>1 & xsupp<7 & ysupp<7
 
                T = tsearch(x3,y3,tri,xsupp,ysupp);
 
                h(1,1)= x3(tri(T,1),1);
                h(2,1)= y3(tri(T,1),1);
                h(3,1)=Grille_Lab2(h(1,1),h(2,1));
 
                h(1,2)= x3(tri(T,2),1);
                h(2,2)= y3(tri(T,2),1);
                h(3,2)=Grille_Lab2(h(1,1),h(2,1));
 
                h(1,3)= x3(tri(T,3),1);
                h(2,3)= y3(tri(T,3),1);
                h(3,3)=Grille_Lab2(h(1,1),h(2,1));
 
 
 
 
                for f=1:3 % calcul des paramètres de l'équation du plan
                    matA(f,1)=1;
                    matA(f,2)=h(1,f);
                    matA(f,3)=h(2,f);
                    matF(f,1)=h(3,f);
                end
 
                matx=inv(matA)*matF;  % les paramètres sont sotckés dans la matrice 'matx'
 
                Z=matx(1,1)+matx(2,1)*xsupp+matx(3,1)*ysupp; % Calcul de la valeur interpolée.
 
                diff(xsupp,ysupp)= zsupp-Z;% difference entre la valeur originale et la valeur interpolée.
 
                matZ(xsupp,ysupp)=Z; % on stocke les valeurs interpolées dans une matrice
 
        end
 
 
 
 
  end
 
  % Sélection des points à supprimer
 
for i=2:6
    absdiff(i-1,:)=abs(diff(i,:));
    Grille_Lab2mod(i-1,:)=Grille_Lab2(i,:);
 
end 
 
 
Grille_Lab3mod=[Grille_Lab2mod(:,2),Grille_Lab2mod(:,3),Grille_Lab2mod(:,4),Grille_Lab2mod(:,5),Grille_Lab2mod(:,6)];
absdiff2=[absdiff(:,2),absdiff(:,3),absdiff(:,4),absdiff(:,5),absdiff(:,6)]
 
for s=1:10 % 10= nombre d'éléments à supprimer
 
    [C,I]= min(absdiff2);
    [D,V]= min(C);
 
    for j=1:5
        if (C(j)~=D)
            E=I(V);
            F=V;
            Grille_Lab3mod(E,F)=NaN;
            absdiff2(E,F)=NaN;
        end
    end
 
end