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