ShenCastan.cpp

1  /***********************************************************************************  
2      ImAnalyse : software in image processing and image analysis 
3     
4      Copyright (C) 10 juillet 2008  <Vincent MORARD> 
5    Version: 2.1 
6      Contact: vincent<POINT>morard<AROBAS>cpe<POINT>fr 
7    Website: http://ImAnalyse.free.fr 
8    Website: http://pistol.petesampras.free.fr 
9   
10      This program is free software: you can redistribute it and/or modify 
11      it under the terms of the GNU General Public License as published by 
12      the Free Software Foundation, either version 3 of the License, or 
13      (at your option) any later version. 
14   
15      This program is distributed in the hope that it will be useful, 
16      but WITHOUT ANY WARRANTY; without even the implied warranty of 
17      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
18      GNU General Public License for more details. 
19   
20      You should have received a copy of the GNU General Public License 
21      along with this program.  If not, see <http://www.gnu.org/licenses/ 
22  **********************************************************************************/ 
23   
24  //****************************************************************************** 
25  //Shen_Castan.cpp 
26  //Ce fichier regroupe les fonctions permettant de trouver les contours selon les  
27  //critères de Shen-Castan.  
28  //OUTLINE corespond au nombre de pixels que l'on ne prendra pas en compte lors  
29  //de la convolution. 
30  //Aussi, pour pouvoir effectuer la détection de contour à tous les pixels, y  
31  //compris ceux du bord de l'image, on augmentera la taille de l'image de  
32  //2*OUTLINE. C'est le rôle des fonctions embded et debed 
33  //WINDOW_SIZE est un paramètre qui determine la taille de la fenêtre pour  
34  //calculer le gradient adaptatif. 
35  // 
36  //plus b est faible plus il y a de détail (plus de bruit aussi) 
37  //****************************************************************************** 
38  #include "../CImage.h" 
39  #include "../AdvancEdge.h" 
40   
41  #define WINDOW_SIZE    7 
42  #define OUTLINE     25 
43  double b; 
44   
45   
46  //************************************************************************* 
47  //ShenCastan: Detection des contours optimaux: 
48  //Procedure: On augmente d'abord la taille de l'image source en on convertit  
49  //l'image initiale en float. On envoie cette  
50  //nouvelle image à la fonction Shen avec ses nouvelles dimensions 
51  //Une fois la fin du traitement, on remet les dimensions initiales et l'on place  
52  //l'image obtenue dans le buffer de sortie. 
53  //************************************************************************* 
54  bool CImage::ShenCastan(CImage *ImgDest,double B) 
55  { 
56    float **Img=0; 
57    int Largeur,Hauteur; 
58   
59    b=B; 
60    if(hBmp==0){ 
61      MessageBox(NULL,"Shen-Castan : L'image source est vide", 
62        NULL,MB_OK|MB_ICONWARNING); 
63      return 0; 
64    } 
65   
66    if(ImgDest != 0 && ImgDest != this) 
67      ImgDest->Copy(this); 
68   
69   
70   
71    if(ImgDest == 0) 
72      ImgDest=this; 
73   
74    GetBitmapBits(hBmp,Width*Height*4,ucBits); 
75   
76    //Ajout d'une bordure de OUTLINE pixel autours de l'image 
77    Img=embed(OUTLINE,&Largeur,&Hauteur); 
78   
79    Shen(Img,Largeur,Hauteur); 
80   
81    //Retrait de la bordure 
82    ImgDest->debed(Img,OUTLINE,Largeur,Hauteur); 
83   
84    SetBitmapBits(ImgDest->hBmp,Width*Height*4,ImgDest->ucBits); 
85    ImgDest->ImgType=GRAY; 
86   
87    DesAllocT_float(Img,Largeur); 
88    return 1; 
89  } 
90   
91   
92  //****************************************************************************** 
93  //Shen:  
94  //C'est dans cette fonction que tout le traitement est effectué 
95  //On reçoit en paramètre d'entrée l'image directement accessible en pixel  
96  //ainsi que ses dimensions. 
97  //L'image de sortie sera dans la variable Img 
98  //****************************************************************************** 
99  void Shen(float **Img,int Largeur,int Hauteur) 
100  { 
101    float **BufFiltrer=0; 
102    int   **ImgBli=0; 
103   
104    BufFiltrer=AllocT_float(Largeur,Hauteur); 
105   
106    //On filtre le bruit en appliquant l'algo ISEF 
107    ComputeISEF(Img,BufFiltrer,Largeur,Hauteur); 
108   
109    //On trouve les pixels de l'image qui on un Laplacian positif  
110    ImgBli=ComputeBli(BufFiltrer,Img,Largeur,Hauteur); 
111   
112    //Detection des coutours à partir de l'ImgBli et de l'image filtrer 
113    LocateZeroCrossings(Img,BufFiltrer,ImgBli,Largeur,Hauteur); 
114   
115    MaxContraste(Img,Largeur,Hauteur); 
116   
117    //Desallocation de la memoire allouer 
118    DesAllocT_float(BufFiltrer,Largeur); 
119    DesAllocT_int(ImgBli,Largeur); 
120   
121   
122   
123  } 
124   
125  //****************************************************************************** 
126  //ISEF  (Infinite Symetrical Exponential Filter) 
127  //Filtrage de l'image horizontalement et verticalement. L'image filtrer sera  
128  //placée dans la variable y 
129  //****************************************************************************** 
130  void ComputeISEF(float **x,float **y,int Largeur,int Hauteur) 
131  { 
132    float **Causal=0,**AntiCausal=0; 
133   
134    Causal=AllocT_float(Largeur,Hauteur); 
135    AntiCausal=AllocT_float(Largeur,Hauteur); 
136   
137    //On applique d'abord le filtre dans la direction verticale 
138    ApplyISEF_Vertical(x,y,Causal,AntiCausal,Largeur,Hauteur); 
139    ApplyISEF_Horizontal(y,y,Causal,AntiCausal,Largeur,Hauteur); 
140   
141   
142   
143    //Libération de la mémoire 
144    DesAllocT_float(Causal,Largeur); 
145    DesAllocT_float(AntiCausal,Largeur); 
146   
147   
148  } 
149   
150   
151  //******************************************************************************* 
152  //ApplyISEF_Vertical 
153  //Filtrage vertical : Calcul des composantes causales et anticausales 
154  //******************************************************************************* 
155  void ApplyISEF_Vertical(float **x,float **y,float **Causal,float **AntiCausal, 
156              int Largeur,int Hauteur) 
157  { 
158    int i,j; 
159    float b1,b2; 
160    b1 = (float)((1.0-b)/(1.0+b)); 
161    b2 = (float)(b*b1); 
162   
163    //Methode recurssive donc on calcule les pixels des bords de l'image 
164    for(j=0;j<Hauteur;j++) 
165    { 
166      Causal[0][j]=b1*x[0][j]; 
167      AntiCausal[Largeur-1][j]=b2*x[Largeur-1][j]; 
168    } 
169   
170    //Calcul des composantes causales 
171    for(i=1;i<Largeur;i++) 
172      for(j=0;j<Hauteur;j++) 
173        Causal[i][j]=(float)(b1*x[i][j]+b*Causal[i-1][j]); 
174   
175    //Calcul des composantes anti-causales 
176    for(i=Largeur-2;i>=0;i--) 
177      for(j=0;j<Hauteur;j++) 
178        AntiCausal[i][j]=(float)(b2*x[i][j]+b*AntiCausal[i+1][j]); 
179   
180    //on calcule les pixels des bords de l'image de sortie 
181    for(j=0;j<Hauteur-1;j++) 
182      y[Largeur-1][j]=Causal[Largeur-1][j]; 
183   
184   
185    //On calcule l'image de sortie du premier filtre que l'on place dans la variable y 
186    //Correspond a la somme des composantes causal et anti-causales 
187    for(i=0;i<Largeur-2;i++) 
188      for(j=0;j<Hauteur-1;j++) 
189        y[i][j]=Causal[i][j]+AntiCausal[i+1][j]; 
190   
191  } 
192   
193  //******************************************************************************* 
194  //ApplyISEF_Horizontal 
195  //Filtrage vertical : Calcul des composantes causales et anticausales 
196  //******************************************************************************* 
197  void ApplyISEF_Horizontal(float **x,float **y,float **Causal,float **AntiCausal, 
198                int Largeur,int Hauteur) 
199  { 
200    int i,j; 
201    float b1,b2; 
202    b1 = (float)((1.0-b)/(1.0+b)); 
203    b2 = (float)(b*b1); 
204   
205    //on calcule les pixels des bords de l'image 
206    for(i=0;i<Largeur;i++) 
207    { 
208      Causal[i][0]=b1*x[i][0]; 
209      AntiCausal[i][Hauteur-2]=b2*x[i][Hauteur-2]; 
210    } 
211   
212    //Calcul des composantes causales 
213    for(j=1;j<Hauteur;j++) 
214      for(i=0;i<Largeur;i++) 
215        Causal[i][j]=(float)(b1*x[i][j]+b*Causal[i][j-1]); 
216   
217    //Compute anti causal component 
218    for(j=Hauteur-3;j>=0;j--) 
219      for(i=0;i<Largeur;i++) 
220        AntiCausal[i][j]=(float)(b2*x[i][j]+b*AntiCausal[i][j+1]); 
221   
222    //Calcul des composantes anti-causales 
223    for(i=0;i<Largeur-1;i++) 
224      y[i][Hauteur-1]=Causal[i][Hauteur-1]; 
225   
226    //On calcule l'image de sortie du premier filtre que l'on place dans la variable y 
227    //Correspond a la somme des composantes causal et anti-causales 
228    for(i=0;i<Largeur-1;i++) 
229      for(j=0;j<Hauteur-2;j++) 
230        y[i][j]=Causal[i][j]+AntiCausal[i][j+1]; 
231   
232  } 
233   
234  //******************************************************************************* 
235  //ComputeBli: 
236  //On fait la différence des deux images et on compare le resultat à 0 
237  //ImgBli est donc une image composé uniquement de 0 et de 1 
238  //******************************************************************************* 
239  int **ComputeBli(float **ImgFiltrer,float **ImgBuf,int nRows,int nCols) 
240  { 
241    int Row,Col; 
242    int **ImgBli=0; 
243   
244    ImgBli=AllocT_int(nRows,nCols); 
245   
246    //On prend la difference entre l'image lisse et l'image originale. 
247    //On calcule l'ImgBli en mettant a 1 tous les pixels ou le Laplacian est positif. 0 sinon 
248    for(Row=0;Row<nRows;Row++) 
249    { 
250      for(Col=0;Col<nCols;Col++) 
251      { 
252        ImgBli[Row][Col]=0; 
253        if(Row<OUTLINE || Row>=nRows-OUTLINE || Col<OUTLINE || Col>=nCols-OUTLINE) 
254          continue; 
255        ImgBli[Row][Col]=((ImgFiltrer[Row][Col]-ImgBuf[Row][Col])>0.0); 
256   
257      } 
258    } 
259    return ImgBli; 
260   
261  } 
262   
263   
264  //************************************************************************************* 
265  //LocateZeroCrossings 
266  //Cette fonction permettra de déterminer pour tous les pixels de l'image s'il est un pixel 
267  //appartenant à un contour ou non. 
268  //************************************************************************************** 
269  void LocateZeroCrossings(float **Orig,float **BufFiltrer,int **ImgBli,int nRows,int nCols) 
270  { 
271    int Row,Col; 
272   
273    for(Row=0;Row<nRows;Row++) 
274      for(Col=0;Col<nCols;Col++) 
275      { 
276        //On ignore les pixels que l'on a ajoute pour le calcule 
277        if(Row<OUTLINE || Row>=nRows-OUTLINE || Col<OUTLINE || Col>=nCols-OUTLINE) 
278          Orig[Row][Col]=0.0; 
279   
280        //On verifie si ce pixel est un "zero crossing" pour le Laplacian 
281        else if(IsCandidateEdge(ImgBli,BufFiltrer,Row,Col)) 
282        { 
283          //On Calcule le gradian adaptatif 
284          Orig[Row][Col]=ComputeAdaptativeGradient(ImgBli,BufFiltrer,Row,Col); 
285   
286        } 
287        else 
288          Orig[Row][Col]=0.0; 
289      } 
290  } 
291   
292  //************************************************************************************** 
293  //IsCandidateEdge 
294  //On regarde le pixel voisin en on regarde s'il y a un franchissement de zero. 
295  //On effectue donc la multiplication des 2 pixels et on compare à 0 
296  //**************************************************************************************  
297  bool IsCandidateEdge(int **Buff,float **Orig,int Row,int Col) 
298  { 
299    //a positive z-c must have a positive 1st derivative,where positive z-c 
300    //means the second derivative goes from + to - as we cross the edge 
301   
302    if(Buff[Row][Col]==1 && Buff[Row+1][Col]==0) 
303      return (Orig[Row+1][Col]-Orig[Row-1][Col]>0 ? TRUE:FALSE); //positive z-c 
304   
305    else if(Buff[Row][Col]==1 && Buff[Row][Col+1]==0) 
306      return (Orig[Row][Col+1]-Orig[Row][Col-1]>0 ? TRUE:FALSE); //positive z-c 
307   
308    else if(Buff[Row][Col]==1 && Buff[Row-1][Col]==0) 
309      return (Orig[Row+1][Col]-Orig[Row-1][Col]<0 ? TRUE:FALSE); //negative z-c 
310   
311    else if(Buff[Row][Col]==1 && Buff[Row][Col-1]==0) 
312      return (Orig[Row][Col+1]-Orig[Row][Col-1]<0 ? TRUE:FALSE); //negative z-c 
313   
314   
315    return FALSE;    //not a z-c 
316   
317  } 
318   
319  //************************************************************************************* 
320  //ComputeAdaptativeGradient 
321  //On calcule un seuil pour chaque pixel. Le seuil sera determiner grace aux pixels voisins 
322  //présent dans la fenêtre WINDOW_SIZE 
323  //************************************************************************************* 
324  float ComputeAdaptativeGradient(int **Bli,float **Orig,int Row,int Col) 
325  { 
326    int i,j; 
327    float SumOn,SumOff; 
328    float AvgOn,AvgOff; 
329    int NbOn,NbOff; 
330   
331    SumOn=0.0; 
332    SumOff=0.0; 
333    NbOn=0; 
334    NbOff=0; 
335   
336    //On regarde par rapport aux pixels voisins, le nombre de pixel a 1 et a 0 
337    for(i=(-WINDOW_SIZE/2);i<=(WINDOW_SIZE/2);i++) 
338      for(j=(-WINDOW_SIZE/2);j<=(WINDOW_SIZE/2);j++) 
339      { 
340        if(Bli[Row+i][Col+j]) 
341        { 
342          SumOn+=Orig[Row+i][Col+j]; 
343          NbOn++; 
344        } 
345        else 
346        { 
347          SumOff+=Orig[Row+i][Col+j]; 
348          NbOff++; 
349        } 
350      } 
351   
352      if(SumOff) 
353        AvgOff=SumOff/(float)NbOff; 
354      else 
355        AvgOff=0.0; 
356   
357      if(SumOn) 
358        AvgOn=SumOn/(float)NbOn; 
359      else 
360        AvgOn=0.0; 
361      return (AvgOff-AvgOn); 
362  } 
363   
364   
365   
366  //******************************************************************************** 
367  //embed 
368  //Cette fonction retourne un pointeur sur un float qui de l'image initial ou l'on a  
369  //ajouté des bordure de largeur et de hauteur de Width 
370  //******************************************************************************** 
371  float **CImage::embed(int taille,int *L,int *H) 
372  { 
373    int i,j,I,J; 
374   
375    float **Img=0; 
376   
377   
378    taille+=2; 
379    *L=Width+2*taille; 
380    *H=Height+2*taille; 
381    Img=AllocT_float(*L,*H); 
382   
383    for(i=0;i<*L;i++) 
384      for(j=0;j<*H;j++) 
385      { 
386        I=(i-taille+Width)%Width; 
387        J=(j-taille+Height)%Height; 
388        Img[i][j]=(float)GetPixel(I,J,GRAY); 
389      } 
390   
391      return Img; 
392  } 
393   
394  //******************************************************************************** 
395  //debed 
396  //Cette fonction retourne inscrit l'image obtenue dans le buffer Dest en enlevant  
397  //les bordures ajoutées 
398  //******************************************************************************** 
399  void CImage::debed(float **Img,int taille,int L,int H) 
400  { 
401    int i,j; 
402    taille+=2; 
403    int Val; 
404    for(i=taille;i<L-taille;i++) 
405      for(j=taille;j<H-taille;j++) 
406      { 
407        Val=Limit((int)Img[i][j]); 
408   
409        SetPixel(i-taille,j-taille,Val,Val,Val); 
410      } 
411  } 
412   
413   
414  //Etire le contraste au maximum. 
415  //Attention tous les pixels negatifs sont mis a 0 
416  void MaxContraste(float **Img,int nRows,int nCols) 
417  { 
418    int i,j; 
419    float X,Scale,Vmin,Vmax; 
420    Vmin=Img[50][50]; 
421    Vmax=Vmin; 
422    for(i=0;i<nRows;i++) 
423    { 
424      for(j=0;j<nCols;j++) 
425      { 
426        if(i<OUTLINE || i>=nRows-OUTLINE || j<OUTLINE || j>=nCols-OUTLINE) 
427          continue; 
428        X=Img[i][j]; 
429        if(X<0){Img[i][j]=0;X=0;} 
430        if(Vmin>X) Vmin=X; 
431        if(Vmax<X) Vmax=X; 
432      } 
433    } 
434   
435    Scale = (float)(256.0/(Vmax-Vmin+1)); 
436   
437    for(i=0;i<nRows;i++) 
438      for(j=0;j<nCols;j++) 
439      { 
440        if(i<OUTLINE || i>=nRows-OUTLINE || j<OUTLINE || j>=nCols-OUTLINE) 
441          continue; 
442   
443        Img[i][j]=((Img[i][j]-Vmin)*Scale); 
444   
445      } 
446   
447   
448   
449  } 

lien1 lien2 lien3