[bwlabel] calculs de statistiques dans des surfaces seuillees

[KrokroAus]

Bonjour a tous,

Je travaille avec des images satellites afin de cartographier la vegetation dans le desert Australien. J'ai un programme qui me permet de calculer des statistiques (superficie, moyenne, valeur max...) dans une region precise, delimité par un seuil et une localisation (i.e je garde les valeurs > seuil et les plus pres de la position definie par les coordonnées X,Y).

Ce programme utilise la fonction bwlabel.
Je voudrais effectuer les memes calculs mais en considerant l'ensemble des "paquets/zones" superieur a mon seuil (sans garder uniquement la zone la plus proche du point X,Y)

J'ai essaye de remplacer dans le code :

Code :

1
2
VEGETcarto = Data.*(L==idx)

par un

Code :

VEGETcarto = Data.*(L)

Les surface delimitees semble bonnes mais pas les valeurs calculees (la carte affichee et les statistiques calculees sont fausses...

Quelqu'un a t'il une idee d'ou vient l'erreur??

Merci beaucoup,
Caroline

Code :

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
function [VEGETcarto,MAXVEGET,MINVEGET,MOYENNEVEGET,MEDIANVEGET,SUPERFICIEVEGET,VALEmb]=Sat_Stat(Data,Seuil,dx,dy,LonX,LatY,yemb)

  %========================================================
 % Calcul de la hauteur et largeur des pixels Modis (en mètre)
 % Rt : rayon terre 
 %========================================================
   Rt=6378140;
   dym=Rt*dy*pi/180./1000.;
   dxm=Rt*dx*pi/180.*cos(pi*yemb/180.)/1000.;

%======================================================== 
% Délimitation des surfaces et calculs des indices statistiques
%========================================================

% Délimitation des surface
%--------------------------

   BW = Data > Seuil;                % matrice avec des 1 qd l inegalite est vraie
   [L, num] = bwlabel(BW);           % On separe les differents paquets
   %figure
   %fig=figure(100);clf 
   %pcolor(L); shading interp;      % Pour observer les régions sélectionnées 

  reg = regionprops(L);             % Proprietes des paquets (regions)
       
   for n=1:num  
      [y,x] = find(L == n);          % Positions associees a une region
      d(n) = min(sqrt((LonX-x).^2+(LatY-y).^2)); % Calcul distance / (LonX,LongY)
      %disp(['n : ' int2str(n) '  d : ' num2str(d(n))]) % Pour afficher les différentes distances calculées
   end
   %idx = find(d<=Dist);
   minim=min(d); % permet de selectionner la "surface" la plus proche d'un point fixe
   idx = find(d<=minim); 
   [pasbesoin,idx2] = sort([reg(idx).Area]);
   idx = idx(idx2(end));

   
%calcul des indices satistiques dans la surface selectionnée
%------------------------------------------------------------
  
%    VEGETcarto = Data.*(L==idx) 
%    VEGET= Data.*(L==idx); 
%    VEGET(find(VEGET==0))=NaN;
%    VEGET=VEGET(finite(VEGET)); 
%    VALEmb = Data(LatY,LonX); 
%    MAXVEGET = max(VEGET(:));
%    MINVEGET = min(VEGET(:));
%    MOYENNEVEGET = mean(mean(VEGET));
%    MEDIANVEGET  = median(median(VEGET));
%    SUPERFICIEVEGET = sum(L(:)==idx(:));
%    disp(['     Nbre Pixels VEGET : ' int2str(SUPERFICIEVEGET)])
%    SUPERFICIEVEGET=SUPERFICIEVEGET*dxm*dym;

%calcul des indices satistiques dans l'ensemble des surfaces
%------------------------------------------------------------

   VEGETcarto = Data.*(L) 
   VEGET= Data.*(L) 
   VEGET(find(VEGET==0))=NaN;
   VEGET=VEGET(finite(VEGET)); 
   VALEmb = Data(LatY,LonX); 
   MAXVEGET = max(VEGET(:));
   MINVEGET = min(VEGET(:));
   MOYENNEVEGET = mean(mean(VEGET));
   MEDIANVEGET  = median(median(VEGET));
   SUPERFICIEVEGET = sum(L(:));
   disp(['     Nbre Pixels VEGET : ' int2str(SUPERFICIEVEGET)])
   SUPERFICIEVEGET=SUPERFICIEVEGET*dxm*dym;
    
   
end

[Dut]

Pour ne prendre en compte toutes les régions d'intérêt et effectuer les calculs d'un seul coup, tu fais "simplement" :

Code :

VEGETcarto = Data.*(L>0)

Comme ça, tu appliques le masque binaire dont les valeurs 1 (vraies) sont celles qui correspondent aux valeurs non nulles de L (donc au numéro des régions d'intérêt détectées)

[KrokroAus]

Super! Merci enormement pour cette reponse rapide!

Pour le calcul de ma surface totale c'est pareil? En ecrivant :

Code :

1
2
3
4
   VEGETcarto = Data.*(L>0) 
   VEGET= Data.*(L>0) 
   SUPERFICIEVEGET = sum(L(:)>0);

J'ai bien la superficie totale i.e de l'ensemble de mes surfaces repondant a la condition valeurImage > seuil ??

merci en tout cas!

[Dut]

Envoyé par KrokroAus

Pour le calcul de ma surface totale c'est pareil? En ecrivant :

Code :

J'ai bien la superficie totale i.e de l'ensemble de mes surfaces repondant a la condition valeurImage > seuil ??

Oui... sauf que les deux premières lignes ne servent à rien dans ce calcul

Sinon dans ton code, tu peux faire quelques améliorations...

Envoyé par KrokroAus

Code :

reg = regionprops(L);

Comme tu ne calcules que les surfaces :

Code :

reg = regionprops(L,'Area');

Envoyé par KrokroAus

Code :

revient à faire "simplement"

Code :

[pasbesoin,idx] = min(d);

Envoyé par KrokroAus

Code :

   VEGET(find(VEGET==0))=NaN;

Avec l'indexage logique :

Code :

VEGET(VEGET==0) = NaN;

Envoyé par KrokroAus

Code :

Comme pour la fonction SUM :

Code :

1
2
   MOYENNEVEGET = mean(VEGET(:));
   MEDIANVEGET  = median(VEGET(:));

Voila, j'espère ne pas avoir écrit de bêtises

[KrokroAus]

Super! Merci beaucoup pour ton aide!
Une derniere question : est t'il possible d'extraire les coordonnees geographiques correspondant aux surface selectionnees, ie repondant a la condition > seuil?

[Dut]

Envoyé par KrokroAus

extraire les coordonnees geographiques correspondant aux surface selectionnees

Qu'entends-tu par "coordonnées géographiques" ?

Ceux des pixels à l'intérieur de la zone ?
Les sommets du polygone frontière ?
Dans quel système d'unité (pixel, métrique, longitude/latitude) ?

[KrokroAus]

J'entends par "coordonnees" les coordonnes (en lat/lon) des pixels a l'interieur des surfaces definies par la condition valeur > seuil.

[Dut]

Pour récupérer les coordonnées des pixels à l'intérieur de chaque région, tu peux faire comme ceci (à adapter à ton problème bien entendu) :

Code :

1
2
3
4
5
[L, num] = bwlabel(...);

for n = 1:num
   [row,col] = find(L==n);
end

Ensuite pour la conversion vers le système de coordonnées final, cela dépend des paramètres de ton image initiale (système de coordonnées, projection, origine, ...)

Note déjà que FIND renvoie d'abord les indices des lignes puis ceux des colonnes... dans un système cartésien, les lignes correspondent aux ordonnées (axe y), et les colonnes aux abscisses (axe x)

[KrokroAus]

Merci pour ton aide et desole pour la reponse tardive!
Je vais tester ca!