Deriche.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  //Deriche 
26  //Ce fichier regroupe certaines fonctions permettant d'effectuer la détection  
27  //de coutour de Deriche. On ne calculera dans ce fichier que l'image lissée  
28  //par le filtre récursif de Deriche. La localisation des zéros et des contours  
29  //se feront dans le fichier Shen_Castan.sfm 
30  //OUTLINE corespond au nombre de pixel que l'on ne prendra pas en compte lors de  
31  //la convolution. 
32  //Aussi, pour pouvoir effectuer la détection de contour à tous les pixels y  
33  //compris ceux du bord de l'image, on augmentera la taille de l'image de  
34  //2*OUTLINE. C'est le rôle des fonctions embded et debed 
35  //******************************************************************************  
36   
37  #include "../CImage.h" 
38  #include "../AdvancEdge.h" 
39   
40  #define OUTLINE   25 
41  float Alpha=0.5; 
42   
43  //****************************************************************************** 
44  //Deriche 
45  //C'est le point d'entrée de cette méthode, (fonction export). Il prend en  
46  //argument l'image source et l'image destination ainsi que la valeur du  
47  //paramètre alpha du traitement. 
48  //****************************************************************************** 
49  bool CImage::Deriche(CImage *ImgDest,float alpha) 
50  { 
51    float **Img=0; 
52    int Largeur,Hauteur; 
53   
54    Alpha=alpha; 
55    if(hBmp==0){ 
56      MessageBox(NULL,"Deriche : L'image source est vide", 
57        NULL,MB_OK|MB_ICONWARNING); 
58      return 0; 
59    } 
60   
61    if(ImgDest != 0 && ImgDest != this) 
62      ImgDest->Copy(this); 
63   
64   
65    if(ImgDest == 0) 
66      ImgDest=this; 
67   
68    GetBitmapBits(hBmp,Width*Height*4,ucBits); 
69   
70    //On ajoute des bords de largeur et de hauteur OUTLINE pour éviter les effets  
71    //de bord. 
72    Img=embed(OUTLINE,&Largeur,&Hauteur); 
73   
74    //On lance le traitement 
75    Deri(Img,Largeur,Hauteur); 
76   
77    //On supprime les bordures créées. 
78    ImgDest->debed(Img,OUTLINE,Largeur,Hauteur); 
79   
80    SetBitmapBits(ImgDest->hBmp,Width*Height*4,ImgDest->ucBits); 
81    ImgDest->ImgType=GRAY; 
82   
83    DesAllocT_float(Img,Largeur); 
84    return 1; 
85  } 
86   
87   
88  //****************************************************************************** 
89  //Deri 
90  //Tout le traitement est effectué dans cette fonction. On reçoit en entrée  
91  //l'image étendue avec de nouvelles bordures ainsi que la nouvelle Largeur et la  
92  //nouvelle hauteur de l'image. 
93  //On calcule donc l'image filtrée, l'image BLI et on effectue la détection des  
94  //zéros. 
95  //****************************************************************************** 
96  void Deri(float **Img,int Largeur,int Hauteur) 
97  { 
98   
99    float **BufFiltrer=0; 
100    int **ImgBli=0; 
101   
102    BufFiltrer=AllocT_float(Largeur,Hauteur); 
103   
104    ComputeSmoothing(Img,BufFiltrer,Largeur,Hauteur); 
105   
106    //On calule une image qui sera egale a 1 lorsque le Laplacian est positif. 0 sinon 
107    ImgBli=ComputeBli(BufFiltrer,Img,Largeur,Hauteur); 
108   
109    //Detection des contours en localisant les passages a 0 
110    LocateZeroCrossings(Img,BufFiltrer,ImgBli,Largeur,Hauteur); 
111   
112    DesAllocT_float(BufFiltrer,Largeur); 
113    DesAllocT_int(ImgBli,Largeur); 
114   
115   
116  } 
117   
118   
119  //****************************************************************************** 
120  //ComputeSmoothing 
121  //Cette fonction permet de calculer récursivement (au sens traitement du signal)  
122  //l'image Lisser suivant les critères définis par Deriche. 
123  //On calcule donc les composantes causales et anticausales de l'image  
124  //d'origine et on fait la somme de ces deux composantes. 
125  //****************************************************************************** 
126  void ComputeSmoothing(float **x,float **y,int Largeur,int Hauteur) 
127  { 
128    float **Causal,**AntiCausal; 
129   
130    Causal=AllocT_float(Largeur,Hauteur); 
131    AntiCausal=AllocT_float(Largeur,Hauteur); 
132   
133    ApplySmoothing_Horizontal(x,y,Causal,AntiCausal,Largeur,Hauteur); 
134    ApplySmoothing_Vertical(y,y,Causal,AntiCausal,Largeur,Hauteur); 
135   
136    DesAllocT_float(Causal,Largeur); 
137    DesAllocT_float(AntiCausal,Largeur); 
138  } 
139   
140   
141  //****************************************************************************** 
142  //ApplySmoothing_Horizontal 
143  //On calcule ici le flou suivant la direction horizontale. On calcule tout  
144  //d'abord les valeurs des coefficients de la fonction définie par Deriche. 
145  //On doit faire attention à la gestion des pixels des bordures: il faut faire  
146  //leur initialisation avant de lancer la procédure récursive. 
147  //La sortie y correspond à la somme des composantes causales et des  
148  //composantes anticausales. 
149  //****************************************************************************** 
150  void ApplySmoothing_Horizontal(float **x,float **y,float **Causal,float **AntiCausal,int Largeur,int Hauteur) 
151  { 
152    float c1,c2,a1,a2,g; 
153    double s1,s2,s; 
154    int i,j; 
155   
156    g=((1-exp(-Alpha))*(1-exp(-Alpha)))/(1+2*Alpha*exp(-Alpha)-exp(-2*Alpha)); 
157    a1=exp(-Alpha); 
158    a2=exp(-2*Alpha); 
159    c1=Alpha+1; 
160    c2=Alpha-1; 
161   
162    for(j=0;j<Hauteur;j++) 
163    { 
164      s2=x[Largeur-1][j]; 
165      s1=s2; 
166   
167      for(i=0;i<2;i++)  // les 2 derniers pixels (-> on duplique le bord) 
168      { 
169        s=s1; 
170        s1=g*a1*c1*x[Largeur-1][j]-g*a2*x[Largeur-1][j]+2*a1*s1-a2*s2; 
171        s2=s; 
172      } 
173      for(i=Largeur-3;i>=0;i--) 
174      { 
175        s=s1; 
176        s1=g*a1*c1*x[i+1][j]-g*a2*x[i+2][j]+2*a1*s1-a2*s2; //initialisation (lissage anticausal) 
177        s2=s; 
178      } 
179      Causal[0][j]=(float)s1; // recopie du point dans les 2 premieres colonnes. 
180      Causal[1][j]=(float)s2; 
181   
182   
183      for(i=2;i<Largeur;i++) 
184      { 
185        Causal[i][j]=(float)(g*x[i][j]+g*a1*c2*x[i-1][j]+2*a1*Causal[i-1][j]-a2*Causal[i-2][j]); //lissage causal 
186        s1=s; 
187        s1=g*x[i][j]+g*a1*c2*x[i-1][j]+2*a1*s1-a2*s2; //deuxieme lissage causal pour initialisation 
188        s2=s; //de l'anticausal 
189      } 
190   
191      AntiCausal[Largeur-1][j]=(float)s1; // recopie du point d'initialisation 
192      AntiCausal[Largeur-2][j]=(float)s2; // dans les 2 dernieres colonnes. 
193   
194      for(i=Largeur-3;i>=0;i--) 
195        AntiCausal[i][j]=(float)(g*a1*c1*x[i+1][j]-g*a2*x[i+2][j]+2*a1*AntiCausal[i+1][j]-a2*AntiCausal[i+2][j]); //lissage anticausal 
196    } 
197   
198    for(i=0;i<Largeur;i++) 
199      for(j=0;j<Hauteur;j++) 
200        y[i][j]=Causal[i][j]+AntiCausal[i][j]; 
201   
202  } 
203   
204   
205  //************************************************************************************** 
206  //ApplySmoothing_Vertical 
207  //Meme fonction que ApplySmoothing_Horizontal sauf que les composante sont créées suivant 
208  //un axe vertical. 
209  //************************************************************************************** 
210  void ApplySmoothing_Vertical(float **x,float **y,float **Causal,float **AntiCausal,int Largeur,int Hauteur) 
211  { 
212    float c1,c2,a1,a2,g; 
213    double s1,s2,s; 
214    int i,j; 
215   
216    g=((1-exp(-Alpha))*(1-exp(-Alpha)))/(1+2*Alpha*exp(-Alpha)-exp(-2*Alpha)); 
217    a1=exp(-Alpha); 
218    a2=exp(-2*Alpha); 
219    c1=Alpha+1; 
220    c2=Alpha-1; 
221   
222    for(i=0;i<Largeur;i++) 
223    { 
224      s2=x[i][Hauteur-1]; 
225      s1=s2; 
226   
227      for(j=0;j<2;j++)  // les 2 derniers pixels (-> on duplique le bord) 
228      { 
229        s=s1; 
230        s1=g*a1*c1*x[i][Hauteur-1]-g*a2*x[i][Hauteur-1]+2*a1*s1-a2*s2; 
231        s2=s; 
232      } 
233      for(j=Hauteur-3;j>=0;j--) 
234      { 
235        s=s1; 
236        s1=g*a1*c1*x[i][j+1]-g*a2*x[i][j+2]+2*a1*s1-a2*s2; //initialisation (lissage anticausal) 
237        s2=s; 
238      } 
239      Causal[i][0]=(float)s1; // recopie du point dans les 2 premieres colonnes. 
240      Causal[i][1]=(float)s2; 
241   
242   
243      for(j=2;j<Hauteur;j++) 
244      { 
245        Causal[i][j]=(float)(g*x[i][j]+g*a1*c2*x[i][j-1]+2*a1*Causal[i][j-1]-a2*Causal[i][j-2]); //lissage causal 
246        s1=s; 
247        s1=g*x[i][j]+g*a1*c2*x[i][j-1]+2*a1*s1-a2*s2; //deuxieme lissage causal pour initialisation 
248        s2=s; //de l'anticausal 
249      } 
250   
251      AntiCausal[i][Hauteur-1]=(float)s1; // recopie du point d'initialisation 
252      AntiCausal[i][Hauteur-2]=(float)s2; // dans les 2 dernieres colonnes. 
253   
254      for(j=Hauteur-3;j>=0;j--) 
255        AntiCausal[i][j]=(float)(g*a1*c1*x[i][j+1]-g*a2*x[i][j+2]+2*a1*AntiCausal[i][j+1]-a2*AntiCausal[i][j+2]); //lissage anticausal 
256    } 
257   
258    for(i=0;i<Largeur;i++) 
259      for(j=0;j<Hauteur;j++) 
260        y[i][j]=Causal[i][j]+AntiCausal[i][j]; 
261   
262  } 

lien1 lien2 lien3