Threashold.cpp

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  #include "CImage.h" 
23   
24   
25   
26  bool CImage::IsolerRed(int SeuilR,int SeuilG,int SeuilB,CImage *ImgDest) 
27  { 
28   
29    if(hBmp==0){ 
30      MessageBox(NULL,"Treashold : L'image source est vide", 
31        NULL,MB_OK|MB_ICONWARNING); 
32      return 0; 
33    } 
34   
35    if(ImgDest!=0 && ImgDest!=this) 
36      ImgDest->Copy(this); 
37   
38    if(ImgDest==0)ImgDest=this; 
39   
40    int i; 
41    //recupération des pixels 
42    GetBitmapBits(hBmp,Width*Height*4,ucBits); 
43   
44    for(i=0;i<Width*Height*4;i+=4) 
45      if(ucBits[i]>SeuilB && ucBits[i+1]> SeuilG && ucBits[i+2]> SeuilR){ 
46        ImgDest->ucBits[i]=255; 
47        ImgDest->ucBits[i+1]=255; 
48        ImgDest->ucBits[i+2]=255; 
49      } 
50      else{ 
51        ImgDest->ucBits[i]=0; 
52        ImgDest->ucBits[i+1]=0; 
53        ImgDest->ucBits[i+2]=0; 
54      } 
55      ImgDest->ImgType=BIN; 
56      SetBitmapBits (ImgDest->hBmp,(Width*Height)*4,ImgDest->ucBits); 
57   
58      return 1; 
59  } 
60   
61   
62   
63  bool CImage::Threashold(int Seuil,CImage *ImgDest) 
64  { 
65   
66    if(hBmp==0){ 
67      MessageBox(NULL,"Treashold : L'image source est vide", 
68        NULL,MB_OK|MB_ICONWARNING); 
69      return 0; 
70    } 
71   
72    if(ImgDest!=0 && ImgDest!=this) 
73      ImgDest->Copy(this); 
74   
75    if(ImgDest==0)ImgDest=this; 
76   
77    int i; 
78    //recupération des pixels 
79    GetBitmapBits(hBmp,Width*Height*4,ucBits); 
80   
81    for(i=0;i<Width*Height*4;i+=4) 
82      if((ucBits[i]+ucBits[i+1]+ucBits[i+2])/3 > Seuil){ 
83        ImgDest->ucBits[i]=255; 
84        ImgDest->ucBits[i+1]=255; 
85        ImgDest->ucBits[i+2]=255; 
86      } 
87      else{ 
88        ImgDest->ucBits[i]=0; 
89        ImgDest->ucBits[i+1]=0; 
90        ImgDest->ucBits[i+2]=0; 
91      } 
92      ImgDest->ImgType=BIN; 
93      SetBitmapBits (ImgDest->hBmp,(Width*Height)*4,ImgDest->ucBits); 
94   
95      return 1; 
96  } 
97   
98   
99   
100  typedef struct _Classe 
101  { 
102    int BorneInf; 
103    int BorneSup; 
104   
105  }Classe; 
106   
107   
108  int Max(int x,int y){ return ( x > y ? x : y );} 
109  int Min(int x,int y){ return ( x < y ? x : y );} 
110   
111  //********************************************************************************* 
112  //Permet de calculer l'histogramme cumulé de la classe C 
113  //********************************************************************************* 
114  int HistCumul(Classe C,int *Hist) 
115  { 
116    int i,Somme; 
117    Somme=0; 
118    for(i=C.BorneInf;i<=C.BorneSup;i++) 
119      Somme+=Hist[i]; 
120    return Somme; 
121   
122  } 
123  //********************************************************************************* 
124  //Permet de calculer la moyenne de la classe C  
125  //La paramètre Puissance donne la possibilité d'élever le resultat au carré avant de 
126  //normaliser 
127  //********************************************************************************* 
128  double HistMoy(Classe C,int *Hist,int Puissance) 
129  { 
130    int tot,i; 
131    double Moyenne; 
132   
133    tot=0; 
134    Moyenne=0; 
135   
136    for(i=C.BorneInf;i<=C.BorneSup;i++) 
137    { 
138      Moyenne+=Hist[i]*i; 
139      tot+=Hist[i]; 
140    } 
141    if(tot!=0) 
142      return pow(Moyenne,Puissance)/tot; 
143    else 
144      return 0; 
145  } 
146   
147  //********************************************************************************* 
148  //Permet de calculer l'écart type de la classe C 
149  //********************************************************************************* 
150  double HistSigma(Classe C,int *Hist) 
151  { 
152    int N,i; 
153    double Sigma,M; 
154    Sigma=0.0; 
155   
156    M=HistMoy(C,Hist,1); 
157    if((N=HistCumul(C,Hist))==0) 
158      return 0.0; 
159   
160    for(i=C.BorneInf;i<=C.BorneSup;i++) 
161      Sigma+=Hist[i]*(i-M)*(i-M); 
162   
163    Sigma/=N; 
164    return sqrt(Sigma); 
165   
166   
167  } 
168   
169   
170   
171  //************************************************************************************* 
172  //thrMaxVariance: 
173  //Cette fonction vise à selectionner le seuil qui sépare au mieux les deux classes par  
174  //maximisation de la variance interclasse 
175  //le seuil optimal s* = 0.5*Moy(C1)+ 0.5*Moy(C2) 
176  //La recherche du seuil optimal se fera par dichotomie.  
177  //Le paramètre boolean Mahalanobis permet d'activer ou non la pondération de Mahalanobis 
178  //s* = Sigma2*Moy(C1)+ Sigma1*Moy(C2)/ (Sigma1+Sigma2) 
179  // 
180  //Cette méthode est inadaptée lorsqu'il faut détecter: 
181  //  -des petites classes 
182  //  -des classes d'écart type très différents (la méthode surestime la classe de faible  
183  //   écart type). Ceci est corrigé par la pondération de Mahalanobis 
184  // 
185  //************************************************************************************** 
186  int CImage::thrMaxVariance(BOOL Mahalanobis) 
187  { 
188    int *Hist=0; 
189    int SEUIL,seuil; 
190    double Sigma0,Sigma1; 
191    Classe C[2]; 
192   
193    C[0].BorneInf=0; 
194    C[0].BorneSup=127; 
195    C[1].BorneInf=128; 
196    C[1].BorneSup=255; 
197   
198    if((Hist=Histogramme(GRAY))==0) 
199      return 0; 
200   
201    seuil=128; 
202    if(!Mahalanobis)    //sans la correction de Mahalanobis 
203      SEUIL=(int)((HistMoy(C[0],Hist,1)+HistMoy(C[1],Hist,1))/2); 
204    else 
205    { 
206      Sigma0=HistSigma(C[0],Hist); 
207      Sigma1=HistSigma(C[1],Hist); 
208      if(Sigma1+Sigma0!=0) 
209        SEUIL=(int)((Sigma1*HistMoy(C[0],Hist,1)+Sigma0*HistMoy(C[1],Hist,1))/(Sigma1+Sigma0)); 
210    } 
211   
212    while(SEUIL!=seuil)    //tant que l'on a pas convergé 
213    { 
214   
215      C[0].BorneSup=SEUIL; 
216      C[1].BorneInf=SEUIL+1; 
217      seuil=SEUIL; 
218      if(!Mahalanobis)    //sans la correction de Mahalanobis 
219        SEUIL=(int)((HistMoy(C[0],Hist,1)+HistMoy(C[1],Hist,1))/2); 
220      else 
221      { 
222        Sigma0=HistSigma(C[0],Hist); 
223        Sigma1=HistSigma(C[1],Hist); 
224        if(Sigma1+Sigma0!=0) 
225          SEUIL=(int)((Sigma1*HistMoy(C[0],Hist,1)+Sigma0*HistMoy(C[1],Hist,1))/(Sigma1+Sigma0)); 
226      } 
227    } 
228    delete []Hist; 
229    return seuil; 
230   
231  } 
232   
233   
234  //************************************************************************** 
235  //thrMaxEntropie 
236  //Cette fonction permet de déterminer le seuil qui sépare au mieux les deux  
237  //classes par maximisation de l'entropie. La pondération par le logarithme 
238  //népérien, permet de favoriser les événements d'autant plus qu'ils sont  
239  //rares. Cette méthode est donc particulièrement bien adaptée à la  
240  //segmentation de petites classes. 
241  //On calculera alors pour chaque niveau et chaque classe l'entropie de la  
242  //répartition des niveaux de gris et on maximisera le résultat. 
243  //************************************************************************** 
244  float Entropie(int Seuil,int *Hist) 
245  { 
246    float E1,E2; 
247    int i,Nb_Pixel; 
248   
249    E1=0.0; 
250    E2=0.0; 
251   
252    Nb_Pixel=0; 
253    for(i=0;i<=Seuil;i++) 
254      Nb_Pixel+=Hist[i]; 
255   
256   
257    for(i=0;i<=Seuil;i++) 
258      if(Hist[i]!=0) 
259        E1+=(float)(Hist[i]/Nb_Pixel*log((double)Hist[i]/Nb_Pixel)); 
260   
261   
262    Nb_Pixel=0; 
263    for(i=Seuil+1;i<=255;i++) 
264      Nb_Pixel+=Hist[i]; 
265   
266    for(i=Seuil+1;i<256;i++) 
267      if(Hist[i]!=0) 
268        E2+=(float)(Hist[i]/Nb_Pixel*log((double)Hist[i]/Nb_Pixel)); 
269   
270    return (-E1-E2); 
271  } 
272   
273  int CImage::thrMaxEntropie() 
274  { 
275    int *Hist=0; 
276    int i,Seuil; 
277    float E,Max; 
278   
279    if((Hist=Histogramme(GRAY))==0) 
280      return 0; 
281   
282    Max=-9999999; 
283    for(i=1;i<255;i++)      //On cherche le seuil qui maximise la fonction entropie 
284      if((E=Entropie(i,Hist))>Max) 
285      { 
286        Max=E; 
287        Seuil=i; 
288      } 
289   
290      delete []Hist; 
291      return Seuil; 
292   
293  } 
294   
295   
296  //******************************************************************************************** 
297  //thrLargeurMiHauteur 
298  //On considère que la répartition de l'histogramme est une gausienne de  
299  //valeur moyenne m et d'écart type sigma  
300  //Si la répartition de l'histogramme ne peut pas etre assimiler à un répartition gaussienne 
301  //cette méthode donnera des résultats abérants 
302  //Un coefficient a été ajouté Coeff: il permet de selectionner un seuil d'un coté ou de l'autre 
303  //de la gaussienne (suivant le signe de Coeff). 
304  //******************************************************************************************* 
305  int CImage::thrLargeurMiHauteur(int Coeff) 
306  { 
307    int *Hist=0,Seuil; 
308    double M,Sigma; 
309    Classe C; 
310   
311    C.BorneInf=0; 
312    C.BorneSup=255; 
313    if((Hist=Histogramme(GRAY))==0) 
314      return 0; 
315   
316    M=HistMoy(C,Hist,1); 
317    Sigma=HistSigma(C,Hist); 
318   
319    Seuil=(int)(M+Coeff*Sigma); 
320    Seuil=(Seuil<0 ? 0 : Min(255,Seuil)); 
321   
322    delete []Hist; 
323    return Seuil; 
324   
325  } 
326   
327  //************************************************************************** 
328  //thrOtsu: 
329  //Il s'agit d'une sorte de maximisation de la variance interclasse.  
330  //Il utilise une fonction critère particulière comme mesure de séparation  
331  //statistique. On calcule donc à partir de l'histogramme normaliser,  
332  //l'histogramme moyenne et l'histogramme cumulé pour tous les niveaux de  
333  //gris de l'image. Histn, HistMoy et HistCumul sont les variables qui  
334  //reçoivent ses données. Il reste plus qu’à calculer pour chaque niveau la  
335  //valeur de la fonction caractéristique et de retenir le niveau qui la  
336  //maximise. (Voir le rapport d’AutoSeuil pour connaître la fonction)  
337  //************************************************************************** 
338  int CImage::thrOtsu() 
339  { 
340    int *Hist=0,i,Seuil; 
341    float Histn[256],HistCumul[256],HistMoy[256],s,Max; 
342   
343    if((Hist=Histogramme(GRAY))==0) 
344      return 0; 
345   
346    //normalisation de l'histogramme  
347    for(i=0;i<256;i++) 
348      Histn[i]=(float)Hist[i]/(Width*Height); 
349   
350    HistCumul[0]=Histn[0]; 
351    HistMoy[0]=0; 
352    Max=-999999; 
353    for(i=1;i<256;i++) 
354    { 
355      HistCumul[i]=HistCumul[i-1]+Histn[i]; 
356      HistMoy[i]=HistMoy[i-1]+i*Histn[i]; 
357    } 
358    for(i=1;i<255;i++) 
359    { 
360      s=(int)HistCumul[i]*(1-HistCumul[i])*pow((HistMoy[255]*HistCumul[i]-HistMoy[i]),2); 
361      if(s>Max) 
362      { 
363        Seuil=i; 
364        Max=s; 
365      } 
366    } 
367    delete []Hist; 
368    return Seuil; 
369   
370  } 
371   
372  //************************************************************************** 
373  //thrPourcentage :  
374  //Cette fonction permet de calculer le seuil tel que le nombre de pixel de  
375  //la classe 1 soit égale au coefficient pourcent multiplié par le nombre de  
376  //pixel total. 
377  //Aussi, si pourcent est égale à 50% alors il y aura le même nombre de pixel  
378  //dans les deux classes. 
379  //Pour calculer facilement le seuil, nus calculons l’histogramme cumulé de  
380  //la répartition des niveaux de gris. Il ne reste plus qu’à parcourir le  
381  //nouveau tableau obtenu pour trouver le seuil optimal. 
382  //************************************************************************** 
383  int CImage::thrPourcentage(double Pourcent) 
384  { 
385    int i,NbPixel; 
386    int *Hist=0,HistCumul[256]; 
387   
388    if(Pourcent<0 || Pourcent>100){ 
389      MessageBox(NULL,"thrPourcentage: Le pourcentage doit etre compris entre 0 et 100", 
390        NULL,MB_OK|MB_ICONWARNING); 
391      return -1; 
392    } 
393   
394   
395    if((Hist=Histogramme(GRAY))==0) 
396      return 0; 
397   
398    NbPixel=Height*Width; 
399   
400    HistCumul[0]=Hist[0]; 
401    for(i=1;i<256;i++) 
402      HistCumul[i]=HistCumul[i-1]+Hist[i]; 
403   
404    for(i=0;i<256;i++) 
405      if((double)HistCumul[i]/HistCumul[255]*100 >= Pourcent) 
406        return i; 
407   
408    return 0; 
409  } 
410   
411   
412  //************************************************************************** 
413  //thrJoh: 
414  //Cette méthode est basée sur la minimisation de l'interdépendance entre  
415  //deux groupes de pixels. On minimise la fonction qui se sert de  
416  //l'entropie de l'image. 
417  //************************************************************************** 
418  float entropy(float h) 
419  { 
420    return (float)(h>0.0 ? -h * (float)log((double)h) : 0.0); 
421   
422  } 
423  int CImage::thrJoh() 
424  { 
425   
426    //pt et pq est l'histogramme cumulé 
427    int i,Seuil,Start,End; 
428    int *Hist=0; 
429    float Histn[256],Pt[256],Pq[256],Sb,Sw,Min; 
430   
431    if((Hist=Histogramme(GRAY))==0) 
432      return 0; 
433   
434    //Normalisation 
435    for(i=0;i<256;i++) 
436      Histn[i]=(float)Hist[i]/(Width*Height); 
437   
438    Min=9999999; 
439    //calcul les facteurs 
440    Pt[0]=Histn[0]; 
441    Pq[0]=1-Histn[0]; 
442    for(i=1;i<256;i++) 
443    { 
444      Pt[i]=Pt[i-1]+Histn[i]; 
445      Pq[i]=1-Pt[i-1]; 
446   
447    } 
448    Start=0; 
449    End=255; 
450    while(Histn[Start++]<=0); 
451    while(Histn[End--]<=0); 
452   
453   
454    //On calcule la fonction à minimiser à tous les niveaux. 
455    for(i=Start;i<End;i++) 
456    { 
457      if(Hist[i]<=0)continue; 
458      Sb=(float)(log((double)Pt[i])+(1.0/Pt[i])*(entropy(Histn[i])+entropy(Pt[i-1]))); 
459      Sw=(float)(log((double)Pq[i])+(1.0/Pq[i])*(entropy(Histn[i])+entropy(Pq[i+1]))); 
460   
461   
462      if(Sb+Sw < Min)    //on sauvegarde le minimum de la fonction 
463      { 
464        Seuil=i; 
465        Min=Sb+Sw; 
466      } 
467    } 
468   
469    return Seuil; 
470  } 
471   
472   
473  //************************************************************************** 
474  //thrMinError 
475  //Cette méthode consiste à considérer l'histogramme comme une répartition  
476  //gaussienne. Il est donc préférable de que cette hypothèse soit vérifiée. 
477  //************************************************************************** 
478  int CImage::thrMinError() 
479  { 
480    double Min=999999; 
481    int *Hist=0; 
482    int P1,P2,i,Seuil,Start,End; 
483    float Sig1,Sig2,J; 
484    Classe C1,C2; 
485   
486    if((Hist=Histogramme(GRAY))==0) 
487      return 0; 
488   
489    Start=0; 
490    End=255; 
491    while(Hist[Start++]<=0); 
492    while(Hist[End--]<=0); 
493    C1.BorneInf=Start; 
494    C2.BorneSup=End; 
495    for(i=Start;i<End;i++) 
496    { 
497      C1.BorneSup=i; 
498      C2.BorneInf=i+1; 
499      P1=HistCumul(C1,Hist); 
500      P2=HistCumul(C2,Hist); 
501      Sig1=(float)HistSigma(C1,Hist); 
502      Sig2=(float)HistSigma(C2,Hist); 
503      if(Sig1==0||Sig2==0)continue; 
504      //on minimise la fonction J 
505      J=(float)(1+2*(P1*log((double)Sig1) + P2*log((double)Sig2))-2*(P1*log((double)P1)+P2*log((double)P2))); 
506   
507   
508      if(J<Min) 
509      { 
510        Seuil=i; 
511        Min=J; 
512      } 
513    } 
514   
515    return Seuil; 
516  } 
517   
518   
519  //************************************************************************** 
520  //thrFuzzy :  
521  //Cette méthode est vu comme un problème ensembliste. Il utilise la fonction  
522  //de Shannon  -x * log(x) - (1-x) * log(1-x) .Cette fonction est définie  
523  //entre 0 et 1. Cela signifie que ux doit aussi être compris entre 0 et 1. 
524  //Comme Ux correspond à des probabilités il n’y a pas de problème sauf  
525  //lorsque ux=0 ou ux=1 
526  //************************************************************************** 
527  double Shannon(double x) 
528  { 
529    return -x*log(x)-(1-x)*log(1-x); 
530  } 
531  int CImage::thrFuzzy() 
532  { 
533   
534    int *Hist=0; 
535    int NbPixel,i,k,Contraste,Start,End,Seuil; 
536    double E,Min,Ux,Mu1,Mu2; 
537    Classe C1,C2; 
538   
539    if((Hist=Histogramme(GRAY))==0) 
540      return 0; 
541    NbPixel=Width*Height; 
542    Start=0; 
543    End=255; 
544   
545   
546    while(Hist[Start++]<=0); 
547    while(Hist[End--]<=0); 
548   
549    C1.BorneInf=Start; 
550    C2.BorneSup=End; 
551   
552    Contraste=End-Start; 
553   
554    Min=9999999; 
555    for(i=Start+1;i<End-1;i++) 
556    { 
557      C1.BorneSup=i; 
558      C2.BorneInf=i+1; 
559      Mu1=HistMoy(C1,Hist,1); 
560      Mu2=HistMoy(C2,Hist,1); 
561   
562      E=0; 
563      for(k=Start;k<=End;k++) 
564      { 
565        if(k<=i) 
566          Ux=(double)Contraste/(Contraste+fabs((double)(k-Mu1))); 
567        else 
568          Ux=(double)Contraste/(Contraste+fabs((double)(k-Mu2))); 
569   
570        if(Ux!=0 && Ux!=1) 
571          E+=Shannon(Ux)*Hist[k]; 
572   
573      } 
574      E=E/NbPixel; 
575   
576      if(E<Min) 
577      { 
578        Min=E; 
579        Seuil=i; 
580      } 
581    } 
582   
583    return Seuil; 
584  } 
585   
586   
587  int CImage::thrFisher() 
588  { 
589    int i,Seuil; 
590    int *Hist=0; 
591    Classe C1,C2; 
592    double Moy1,Moy2,J,Max; 
593   
594    if((Hist=Histogramme(GRAY))==0) 
595      return 0; 
596   
597    C1.BorneInf=0; 
598    C2.BorneSup=255; 
599    Max=-9999999; 
600    for(i=1;i<254;i++) 
601    { 
602      C1.BorneSup=i; 
603      C2.BorneInf=i+1; 
604      Moy1=HistMoy(C1,Hist,2); 
605      Moy2=HistMoy(C2,Hist,2); 
606      J=Moy1+Moy2; 
607   
608      if(J>Max) 
609      { 
610        Seuil=i; 
611        Max=J; 
612      } 
613    } 
614   
615    return Seuil; 
616   
617  } 
618   
619   
620  int CImage::thrTwoPeak() 
621  { 
622    int *Hist=0; 
623    int i,IndexMax1,IndexMax2,Max,Min,Seuil; 
624   
625   
626    if((Hist=Histogramme(GRAY))==0) 
627      return 0; 
628   
629    Max=0; 
630    //Détection de la première zone maximale :first peak 
631    for(i=0;i<256;i++) 
632      if(Hist[i]>=Max) 
633      { 
634        Max=Hist[i]; 
635        IndexMax1=i; 
636      } 
637      Max=0; 
638      //Détection de la seconde zone max :second peak; 
639      //On soustrait et on élève au carré pour ne pas etre géné par le premier peak 
640      for(i=0;i<256;i++) 
641        if((i-IndexMax1)*(i-IndexMax1)*Hist[i] > Max) 
642        { 
643          Max=(i-IndexMax1)*(i-IndexMax1)*Hist[i]; 
644          IndexMax2=i; 
645        } 
646   
647   
648        if(IndexMax1>IndexMax2) 
649        { 
650          Min=IndexMax1; 
651          IndexMax1=IndexMax2; 
652          IndexMax2=Min; 
653        } 
654   
655        Min=Hist[IndexMax2]; 
656   
657        for(i=IndexMax1;i<IndexMax2;i++) 
658          if(Hist[i]<Min) 
659          { 
660            Seuil=i; 
661            Min=Hist[i]; 
662          } 
663   
664          return Seuil; 
665  } 

lien1 lien2 lien3