LPE (Ligne de Partage des Eaux)

1  /***********************************************************************************  
2      ImAnalyse : software in image processing and image analysis 
3     
4      Copyright (C) 27 avril 2008  <Vincent MORARD> 
5    Version: 2.0 
6      Contact: vincent<POINT>morard<AROBAS>cpe<POINT>fr 
7    Website: http://pistol.petesampras.free.fr 
8   
9      This program is free software: you can redistribute it and/or modify 
10      it under the terms of the GNU General Public License as published by 
11      the Free Software Foundation, either version 3 of the License, or 
12      (at your option) any later version. 
13   
14      This program is distributed in the hope that it will be useful, 
15      but WITHOUT ANY WARRANTY; without even the implied warranty of 
16      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
17      GNU General Public License for more details. 
18   
19      You should have received a copy of the GNU General Public License 
20      along with this program.  If not, see <http://www.gnu.org/licenses/ 
21  **********************************************************************************/ 
22   
23  //*************************************************************************************** 
24  //Classe  LPE : (ligne de partage des eaux)  
25  //Auteur:  MORARD Vincent  
26  //Date :  14 juin 2007 
27  //Principde de S.Beucher 
28  //************************************************************************************** 
29  #include "CImage.h" 
30  #include <queue> 
31  #include "AdvancEdge.h"    //Pour AllocT et desalloc 
32  using namespace std; 
33   
34  #define N 256 
35  #define COLLISION -1 
36  #define F_COLLISION    257 
37   
38  int **ImgLabel,**ImgF; 
39  int Priority_Queue,m_Neighborhood; 
40  queue<POINT>m_FAH[N+2]; 
41   
42   
43   
44   
45  void CImage::LPE_Init(int Label,int Neighborhood) 
46  { 
47    if((ImgF=Convert_2D(GRAY))==NULL) 
48      return; 
49    if((ImgLabel=Img[Label].Convert_2D(GRAY))==NULL) 
50      return; 
51   
52    for(int i=0;i<Width;i++) 
53    { 
54      for(int j=0;j<Height;j++) 
55        if(ImgLabel[i][j]==255) 
56          ImgLabel[i][j]=60000; 
57    } 
58    Labeliser(); 
59    Init_Queue(); 
60    Priority_Queue=0; 
61    if(Neighborhood == 4 || Neighborhood==8) 
62      m_Neighborhood=Neighborhood; 
63    else 
64      m_Neighborhood=4; 
65   
66  } 
67   
68  //****************************************************************************** 
69  //LabelProcessing 
70  //Cette fonction permet d’initialiser l’image des étiquettes. A partir de  
71  //l’image des marqueurs, on labellise tous les pixels connexes par la même  
72  //valeur (Num). Pour atteindre toutes les piles connexes, on utilise une pile.  
73  //On empilera alors tous les pixels voisins qui seront au niveau logique 1. La  
74  //pile simule la récursivité. Cette méthode est bien plus efficace car elle  
75  //évite tous les changements de contexte. On gagne donc largement en temps  
76  //d’exécution 
77  //****************************************************************************** 
78  void CImage::LabelProcessing(int i,int j,int Num) 
79  { 
80    int k,l,m; 
81    POINT Pt; 
82    Pt.x=i; 
83    Pt.y=j; 
84    m_FAH[0].push(Pt); 
85    do 
86    { 
87      Pt=m_FAH[0].front(); 
88      m_FAH[0].pop(); 
89      ImgLabel[Pt.x][Pt.y]=Num; 
90      i=Pt.x; 
91      j=Pt.y; 
92      for(k=-1;k<=1;k++) 
93        for(l=-1;l<=1;l++) 
94        { 
95          Pt.x=i+k; 
96          Pt.y=j+l; 
97          if(Pt.x>=0 && Pt.x<Width && Pt.y>=0 && Pt.y<Height) 
98          { 
99            m=ImgLabel[Pt.x][Pt.y]; 
100            if(m==60000) 
101            { 
102              m_FAH[0].push(Pt); 
103              ImgLabel[Pt.x][Pt.y]=59000; 
104            } 
105          } 
106        } 
107    } 
108    while(!m_FAH[0].empty()); 
109   
110  } 
111   
112   
113  //****************************************************************************** 
114  //Labeliser 
115  //Cette fonction va de paire avec LabelProcessing. En effet dans cette fonction  
116  //on parcourt l’image entière et pour chaque pixel à 255 on appelle la fonction  
117  //précédente. 
118  //On incrémente Num pour signifier une nouvelle zone. 
119  //****************************************************************************** 
120  void CImage::Labeliser() 
121  { 
122    int Num=1; 
123    for(int i=0;i<Width;i++) 
124      for(int j=0;j<Height;j++) 
125        if(ImgLabel[i][j]==60000) 
126        { 
127          LabelProcessing(i,j,Num); 
128          Num+=2; 
129        } 
130   
131  } 
132   
133  //****************************************************************************** 
134  //InitQueue 
135  //Cette fonction permet l’initialisation de la FAH. On va donc empiler tous les  
136  //pixels qui ont été marqués dans les files de niveau correspondant à leur  
137  //niveau de gris. Cette étape marque la fin de la séquence d’initialisation. 
138  //****************************************************************************** 
139  void CImage::Init_Queue() 
140  { 
141    POINT Pt; 
142    for(Pt.x=0 ; Pt.x < Width ;  Pt.x++) 
143      for(Pt.y=0 ; Pt.y < Height ; Pt.y++) 
144        if((ImgLabel[Pt.x][Pt.y])!=0) 
145          m_FAH[ImgF[Pt.x][Pt.y]].push(Pt); 
146   
147  } 
148  bool CImage::IsEmptyFAH() 
149  { 
150    for(int i=0;i<N;i++) 
151    { 
152      if(!m_FAH[i].empty()) 
153        return false; 
154    } 
155      return true; 
156  } 
157  bool CImage::IsEmptyFAH(int num) 
158  { 
159    return m_FAH[num].empty(); 
160  } 
161   
162   
163  //****************************************************************************** 
164  //NeighborProcess 
165  //Cette fonction exécute le traitement du pixel voisin PtNeighbor par rapport au  
166  //pixel (PT) Suivant la valeur des étiquettes et des niveaux de gris  
167  //des pixels, on effectue des actions différentes. 
168  //****************************************************************************** 
169  void CImage::NeighborProcess(POINT PtNeighbor,POINT Pt,int n) 
170  { 
171    int g=ImgLabel[Pt.x][Pt.y]; 
172    int f=ImgF[Pt.x][Pt.y]; 
173    int gn=ImgLabel[PtNeighbor.x][PtNeighbor.y]; 
174    int fn=ImgF[PtNeighbor.x][PtNeighbor.y]; 
175   
176    if(gn==0) 
177    { 
178      //on assigne au voisin y la valeur de l'étiquette de x-1: label temporaire 
179      ImgLabel[PtNeighbor.x][PtNeighbor.y]=g-1; 
180   
181      if(fn > Priority_Queue) 
182      { 
183        //On place dans la file FAH sa priorité correspondant au niveau 
184        //de gris du pixel voisin dans F   
185        m_FAH[fn].push(PtNeighbor); 
186      } 
187      else 
188      { 
189        m_FAH[N].push(PtNeighbor); 
190        ImgF[PtNeighbor.x][PtNeighbor.y]=N+n; 
191   
192      } 
193    } 
194    else 
195    { 
196      if(gn%2==1)          //si le label est impaire 
197        if(f != fn) 
198          if(gn != g-1)  //Si l'inondation ne provient pas du meme bassin versant 
199          { 
200            ImgLabel[PtNeighbor.x][PtNeighbor.y]=COLLISION;//il y a collision 
201            m_FAH[F_COLLISION].push(PtNeighbor); 
202          } 
203    } 
204  } 
205  //****************************************************************************** 
206  //FindNeighbor 
207  //Cette fonction permet de trouver les voisins d’un pixel (i,j) et de lancer la  
208  //fonction permettant le traitement de ce pixel voisin. 
209  //****************************************************************************** 
210  void CImage::FindNeigbor(POINT Pt,int n) 
211  { 
212    POINT PtNeighbor; 
213   
214   
215    for(int i=-1;i<=1;i++) 
216      for(int j=-1;j<=1;j++) 
217      { 
218        PtNeighbor.x=Pt.x+i; 
219        PtNeighbor.y=Pt.y+j; 
220        if(!(PtNeighbor.x < 0 || PtNeighbor.x>=Width 
221          || PtNeighbor.y < 0 || PtNeighbor.y>=Height)) 
222        { 
223          if(m_Neighborhood==8 || (m_Neighborhood==4 && (i!=j && i!=-j))) 
224            NeighborProcess(PtNeighbor,Pt,n); 
225        } 
226      } 
227  } 
228   
229  //****************************************************************************** 
230  //SelectNextQueue 
231  //Cette fonction parcourt la FAH depuis le niveau marqué par PriorityQueue  
232  //jusqu’à ce que l’on ait une file d’attente qui ne soit pas vide. Si toutes les  
233  //files d’attente sont vides, on retourne -1. 
234  //****************************************************************************** 
235  int CImage::SelectNextQueue() 
236  { 
237    for(int i=Priority_Queue;i<N;i++) 
238    { 
239      if(!IsEmptyFAH(i)) 
240        return i; 
241    } 
242    return -1; 
243  } 
244   
245  //****************************************************************************** 
246  //PixelProcess 
247  //Cette fonction permet de traiter le pixel courant. On lui donne alors un label  
248  //définitif et on regarde ses voisins. 
249  //****************************************************************************** 
250  void CImage::PixelProcess(POINT Pt,int n) 
251  { 
252    if(ImgLabel[Pt.x][Pt.y] == COLLISION)      //si ce pixel est en collision 
253    {                    //On ne propage pas  
254      ImgF[Pt.x][Pt.y]=Priority_Queue; 
255      return;                //le pixel en collision sort de la file sans rien faire 
256    } 
257    //on lui donne un label définitif 
258    ImgLabel[Pt.x][Pt.y]++; 
259   
260    //on recherche ses voisins 
261    FindNeigbor(Pt,n); 
262   
263    //S'il ne font pas partie de la meme génération: 
264    if(n>1) 
265      ImgF[Pt.x][Pt.y]=Priority_Queue; 
266   
267   
268  } 
269   
270   
271  //Label = image des marqueurs 
272  bool CImage::LPE(int Label,int ImgDest,int Neighborhood) 
273  { 
274    if(hBmp==0 || Img[Label-1].hBmp==0){ 
275      MessageBox(NULL,"LPE : L'image source est vide", 
276            NULL,MB_OK|MB_ICONWARNING); 
277      return 0; 
278    } 
279   
280    POINT Pt; 
281    int n=1; 
282      Label--;ImgDest--; 
283    if(hBmp==0) 
284      return 0; 
285    if(Img[Label].hBmp==0) 
286      return 0; 
287    if(Width!=Img[Label].Width || Height != Img[Label].Height) 
288      return 0; 
289   
290    LPE_Init(Label,Neighborhood); 
291   
292   
293    int CurrentQueue=SelectNextQueue(); 
294    Priority_Queue=CurrentQueue; 
295    if(CurrentQueue==-1) 
296    { 
297      printf("fonction LPE : pas de pixel empiler lors de l'initialisation\n"); 
298      return 0; 
299    } 
300    while(Priority_Queue < N || ! m_FAH[N].empty()) 
301    { 
302   
303      //Afficher2(m_Label); 
304      Pt=m_FAH[CurrentQueue].front();    //on reccupère le pixel suivant 
305      m_FAH[CurrentQueue].pop();      //on dequeue 
306   
307      //On traite ce pixel 
308      PixelProcess(Pt,n); 
309   
310      //Si la pile courrante est vide: 
311      if(IsEmptyFAH(CurrentQueue)) 
312      { 
313        if(CurrentQueue==N)      //si on a vider la file alpha 
314        { 
315          Priority_Queue=SelectNextQueue(); 
316          CurrentQueue=Priority_Queue; 
317          n=1; 
318        } 
319        else if(IsEmptyFAH(N))          //si la file alpha est vide  
320          { 
321            Priority_Queue=SelectNextQueue(); 
322            CurrentQueue=Priority_Queue; 
323            n=1; 
324          } 
325          else 
326          { 
327            Priority_Queue++; 
328            CurrentQueue=N;        //on va écouler le file alpha 
329            n=n+1; 
330          } 
331      } 
332      if(CurrentQueue==-1 ) 
333      { 
334        if(!IsEmptyFAH(N)) 
335        { 
336          Priority_Queue++; 
337          CurrentQueue=N;        //on va écouler le file alpha 
338          n=n+1; 
339        } 
340        else 
341          Priority_Queue=N; 
342      } 
343   
344    } 
345   
346    Img[Label].PixelCollision(); 
347   
348    Img[ImgDest].Copy(&Img[Label]); 
349    Img[ImgDest].Convert_2D_to_1D(ImgLabel); 
350    Img[ImgDest].ImgType=GRAY; 
351   
352    DesAllocT_int(ImgLabel,Width); 
353    DesAllocT_int(ImgF,Width); 
354    return 1; 
355  } 
356   
357   
358  //****************************************************************************** 
359  //PixelCollision 
360  //A la fin du traitement on va labelliser les pixels qui ont un label COLLISION.  
361  //On va donc dépiler la file d’attente numéro 257 (F_COLLISION) et attribuer un  
362  //label identique à un de ses voisins. 
363  //****************************************************************************** 
364  void CImage::PixelCollision() 
365  { 
366    int i,j,k,l,p,M; 
367    POINT Pt; 
368    bool OK; 
369   
370    while(!m_FAH[F_COLLISION].empty()) 
371    { 
372    Pt=m_FAH[F_COLLISION].front();    //on reccupère le pixel suivant 
373      m_FAH[F_COLLISION].pop();      //on dequeue 
374      i=Pt.x; 
375    j=Pt.y; 
376    OK=false; 
377      for(p=0;p<9;p++) 
378      { 
379    //on utilise une unique boucle for pour que lorsque l’on break,  
380    //on passe directement au pixel suivant.  
381        k=p%3-1; 
382        l=p/3-1; 
383        if(i+k>=0 && i+k<Width && j+l>=0 && j+l<Height) 
384          if((M=ImgLabel[i+k][j+l])!=COLLISION) 
385          { 
386        OK=true; 
387            ImgLabel[i][j]=M; 
388            break; 
389          } 
390      } 
391    if(!OK) 
392      m_FAH[F_COLLISION].push(Pt); 
393    } 
394  } 
395   
396   
397   
398   
399  bool CImage::lpeFrontiere(CImage  *ImgContour,CImage  *ImgDest) 
400  { 
401    if(hBmp==0){ 
402      MessageBox(NULL,"lpeFrontiere : L'image source est vide", 
403            NULL,MB_OK|MB_ICONWARNING); 
404      return 0; 
405    } 
406   
407    if(ImgDest!=0 && ImgDest!=this) 
408      ImgDest->CopyProp(this); 
409    if(ImgDest==0) 
410      ImgDest=this; 
411   
412    int i,j,P1; 
413    for(i=0;i<Width-1;i++) 
414      for(j=0;j<Height-1;j++) 
415      { 
416        if((P1=(int)GetPixel(i,j,GRAY)) != GetPixel(i+1,j,GRAY) ) 
417        { 
418          if(ImgContour->GetPixel(i,j,RED)==255) 
419            ImgDest->SetPixel(i,j,220,0,0); 
420          else 
421            ImgDest->SetPixel(i,j,255,128,0); 
422        } 
423        else if(P1 != GetPixel(i,j+1,GRAY) ) 
424        { 
425          if(ImgContour->GetPixel(i,j,RED)==255) 
426            ImgDest->SetPixel(i,j,220,0,0); 
427          else 
428            ImgDest->SetPixel(i,j,255,128,0); 
429        } 
430        else 
431          ImgDest->SetPixel(i,j,0,0,0); 
432      } 
433   
434    ImgDest->ImgType=RGBi; 
435    SetBitmapBits (ImgDest->hBmp,(Width*Height)*4,ImgDest->ucBits); 
436    return 1; 
437   
438  } 

lien1 lien2 lien3