improved png compression
[swftools.git] / lib / gocr / otsu.c
1 /*
2 This is a Optical-Character-Recognition program
3 Copyright (C) 2000-2007  Joerg Schulenburg
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18
19  see README for EMAIL-address
20
21  the following code was send by Ryan Dibble <dibbler@umich.edu>
22  
23   The algorithm is very simple but works good hopefully.
24  
25   Compare the grayscale histogram with a mass density diagram:
26   I think the algorithm is a kind of
27   divide a body into two parts in a way that the mass 
28   centers have the largest distance from each other,
29   the function is weighted in a way that same masses have a advantage
30   
31   - otsu algorithm is failing on diskrete multi color images
32   
33   TODO:
34     RGB: do the same with all colors (CMYG?) seperately 
35
36     test: hardest case = two colors
37        bbg: test done, using a two color gray file. Output:
38        # threshold: Value = 43 gmin=43 gmax=188
39
40  my changes:
41    - float -> double
42    - debug option added (vvv & 1..2)
43    - **image => *image,  &image[i][1] => &image[i*cols+1]
44    - do only count pixels near contrast regions
45      this makes otsu much better for shadowed fonts or multi colored text 
46      on white background
47
48  (m) Joerg Schulenburg (see README for email address)
49
50  ToDo:
51    - measure contrast
52    - detect low-contrast regions
53
54  */
55
56 #include <stdio.h>
57 #include <string.h>
58
59 #define Abs(x) ((x<0)?-(x):x)
60
61 /*======================================================================*/
62 /* global thresholding routine                                     */
63 /*   takes a 2D unsigned char array pointer, number of rows, and        */
64 /*   number of cols in the array. returns the value of the threshold    */
65 /*======================================================================*/
66 int
67 otsu (unsigned char *image, int rows, int cols,
68       int x0, int y0, int dx, int dy, int vvv) {
69
70   unsigned char *np;    // pointer to position in the image we are working with
71   unsigned char op1, op2;   // predecessor of pixel *np (start value)
72   int maxc=0;           // maximum contrast (start value) 
73   int thresholdValue=1; // value we will threshold at
74   int ihist[256];       // image histogram
75   int chist[256];       // contrast histogram
76
77   int i, j, k;          // various counters
78   int is, i1, i2, ns, n1, n2, gmin, gmax;
79   double m1, m2, sum, csum, fmax, sb;
80
81   // zero out histogram ...
82   memset(ihist, 0, sizeof(ihist));
83   memset(chist, 0, sizeof(chist));
84   op1=op2=0;
85
86   gmin=255; gmax=0; k=dy/512+1;
87   // v0.43 first get max contrast, dont do it together with next step
88   //  because it failes if we have pattern as background (on top)
89   for (i =  0; i <  dy ; i+=k) {
90     np = &image[(y0+i)*cols+x0];
91     for (j = 0; j < dx ; j++) {
92       ihist[*np]++;
93       if(*np > gmax) gmax=*np;
94       if(*np < gmin) gmin=*np;
95       if (Abs(*np-op1)>maxc) maxc=Abs(*np-op1); /* new maximum contrast */
96       if (Abs(*np-op2)>maxc) maxc=Abs(*np-op2); /* new maximum contrast */
97       /* we hope that maxc will be find its maximum very fast */
98       op2=op1;    /* shift old pixel to next older */
99       op1=*np;     /* store old pixel for contrast check */
100       np++;       /* next pixel */
101     }
102   }
103
104   // generate the histogram
105   // Aug06 images with large white or black homogeneous
106   //   areas give bad results, so we only add pixels on contrast edges
107   for (i =  0; i <  dy ; i+=k) {
108     np = &image[(y0+i)*cols+x0];
109     for (j = 0; j < dx ; j++) {
110       if (Abs(*np-op1)>maxc/4
111        || Abs(*np-op2)>maxc/4)
112          chist[*np]++; // count only relevant pixels
113       op2=op1;    /* shift old pixel to next older */
114       op1=*np;     /* store old pixel for contrast check */
115       np++;       /* next pixel */
116     }
117   }
118
119   // set up everything
120   sum = csum = 0.0;
121   ns = 0;
122   is = 0;
123   
124   for (k = 0; k <= 255; k++) {
125     sum += (double) k * (double) chist[k];  /* x*f(x) cmass moment */
126     ns  += chist[k];                        /*  f(x)    cmass      */
127     is  += ihist[k];                        /*  f(x)    imass      */
128     // Debug: output to out_hist.dat?
129     // fprintf(stderr,"\chistogram %3d %6d (brightness weight)", k, ihist[k]);
130   }
131
132   if (!ns) {
133     // if n has no value we have problems...
134     fprintf (stderr, "NOT NORMAL, thresholdValue = 160\n");
135     return (160);
136   }
137
138   // ToDo: only care about extremas in a 3 pixel environment
139   //       check if there are more than 2 mass centers (more colors)
140   //       return object colors and color radius instead of threshold value
141   //        also the reagion, where colored objects are found
142   //       what if more than one background color? no otsu at all?
143   //       whats background? box with lot of other boxes in it
144   //       threshold each box (examples/invers.png,colors.png) 
145   //  get maximum white and minimum black pixel color (possible range)
146   //    check range between them for low..high contrast ???
147   // typical scenes (which must be covered): 
148   //    - white page with text of different colors (gray values)
149   //    - binear page: background (gray=1) + black text (gray=0)
150   //    - text mixed with big (dark) images
151   //  ToDo: recursive clustering for maximum multipol moments?
152   //  idea: normalize ihist to max=1024 before otsu?
153   
154   // do the otsu global thresholding method
155
156   if ((vvv&1)) // Debug
157      fprintf(stderr,"# threshold: value ihist chist mass_dipol_moment\n");
158   fmax = -1.0;
159   n1 = 0;
160   for (k = 0; k < 255; k++) {
161     n1 += chist[k];          // left  mass (integration)
162     if (!n1) continue;       // we need at least one foreground pixel
163     n2 = ns - n1;            // right mass (num pixels - left mass)
164     if (n2 == 0) break;      // we need at least one background pixel
165     csum += (double) k *chist[k];  // left mass moment
166     m1 =        csum  / n1;        // left  mass center (black chars)
167     m2 = (sum - csum) / n2;        // right mass center (white background)
168     // max. dipol moment?
169     // orig: sb = (double) n1 *(double) n2 * (m1 - m2) * (m1 - m2);
170     sb = (double) n1 *(double) n2 * (m2 - m1); // seems to be better Aug06
171     /* bbg: note: can be optimized. */
172     if (sb > fmax) {
173       fmax = sb;
174       thresholdValue = k + 1;
175       // thresholdValue = (m1 + 3 * m2) / 4;
176     }
177     if ((vvv&1) && ihist[k]) // Debug
178      fprintf(stderr,"# threshold: %3d %6d %6d %8.2f\n",
179        k, ihist[k], chist[k],
180        sb/(dx*dy));  /* normalized dipol moment */
181   }
182   // ToDo: error = left/right point where sb is 90% of maximum?
183   // now we count all pixels for background detection
184   i1 = 0;
185   for (k = 0; k < thresholdValue; k++) {
186     i1 += ihist[k];          // left  mass (integration)
187   }
188   i2 = is - i1;            // right mass (num pixels - left mass)
189
190   // at this point we have our thresholding value
191   // black_char: value<cs,  white_background: value>=cs
192   
193   // can it happen? check for sureness
194   if (thresholdValue >  gmax) {
195     fprintf(stderr,"# threshold: Value >gmax\n");
196     thresholdValue = gmax;
197   }
198   if (thresholdValue <= gmin) {
199     fprintf(stderr,"# threshold: Value<=gmin\n");
200     thresholdValue = gmin+1;
201   }
202
203   // debug code to display thresholding values
204   if ( vvv & 1 )
205   fprintf(stderr,"# threshold: Value = %d gmin=%d gmax=%d cmax=%d"
206                  " i= %d %d\n",
207      thresholdValue, gmin, gmax, maxc, i1, i2);
208
209   if (i1>=4*i2) { // black>=4*white, obviously black is background
210     if ( vvv & 1 )
211       fprintf(stderr,"# threshold: invert the image\n");
212     // we do inversion here (no data lost)
213     for (i =  0; i <  dy ; i++) {
214       np = &image[(y0+i)*cols+x0];
215       for (j = 0; j < dx ; j++) {
216         *np=255-*np;
217         np++;       /* next pixel */
218       }
219     }
220     thresholdValue=255-thresholdValue+1;
221   }
222
223   return(thresholdValue); 
224   /* range: 0 < thresholdValue <= 255, example: 1 on b/w images */
225   /* 0..threshold-1 is foreground */
226   /* threshold..255 is background */
227   /* ToDo:  min=blackmasscenter/2,thresh,max=(whitemasscenter+255)/2 */
228 }
229
230 /*======================================================================*/
231 /* thresholding the image  (set threshold to 128+32=160=0xA0)           */
232 /*   now we have a fixed thresholdValue good to recognize on gray image */
233 /*   - so lower bits can used for other things (bad design?)            */
234 /* ToDo: different foreground colors, gray on black/white background    */
235 /*======================================================================*/
236 int
237 thresholding (unsigned char *image, int rows, int cols,
238   int x0, int y0, int dx, int dy, int thresholdValue) {
239
240   unsigned char *np; // pointer to position in the image we are working with
241
242   int i, j;          // various counters
243   int gmin=255,gmax=0;
244   int nmin=255,nmax=0;
245
246   // calculate min/max (twice?)
247   for (i = y0 + 1; i < y0 + dy - 1; i++) {
248     np = &image[i*cols+x0+1];
249     for (j = x0 + 1; j < x0 + dx - 1; j++) {
250       if(*np > gmax) gmax=*np;
251       if(*np < gmin) gmin=*np;
252       np++; /* next pixel */
253     }
254   }
255   
256   /* allowed_threshold=gmin+1..gmax v0.43 */
257   if (thresholdValue<=gmin || thresholdValue>gmax){
258     thresholdValue=(gmin+gmax+1)/2; /* range=0..1 -> threshold=1 */
259     fprintf(stderr,"# thresholdValue out of range %d..%d, reset to %d\n",
260      gmin, gmax, thresholdValue);
261   }
262
263   /* b/w: min=0,tresh=1,max=1 v0.43 */
264   // actually performs the thresholding of the image...
265   //  later: grayvalues should also be used, only rescaling threshold=160=0xA0
266   for (i = y0; i < y0+dy; i++) {
267     np = &image[i*cols+x0];
268     for (j = x0; j < x0+dx; j++) {
269       *np = (unsigned char) (*np >= thresholdValue ?
270          (255-(gmax - *np)* 80/(gmax - thresholdValue + 1)) :
271          (  0+(*np - gmin)*150/(thresholdValue - gmin    )) );
272       if(*np > nmax) nmax=*np;
273       if(*np < nmin) nmin=*np;
274       np++;
275     }
276   }
277
278   // fprintf(stderr,"# thresholding:  nmin=%d nmax=%d\n", nmin, nmax);
279
280   return(128+32); // return the new normalized threshold value
281   /*   0..159 is foreground */
282   /* 160..255 is background */
283 }
284