1 /**********************************************************************
2 * File: otsuthr.cpp
3 * Description: Simple Otsu thresholding for binarizing images.
4 * Author: Ray Smith
5 * Created: Fri Mar 07 12:31:01 PST 2008
6 *
7 * (C) Copyright 2008, Google Inc.
8 ** Licensed under the Apache License, Version 2.0 (the "License");
9 ** you may not use this file except in compliance with the License.
10 ** You may obtain a copy of the License at
11 ** http://www.apache.org/licenses/LICENSE-2.0
12 ** Unless required by applicable law or agreed to in writing, software
13 ** distributed under the License is distributed on an "AS IS" BASIS,
14 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 ** See the License for the specific language governing permissions and
16 ** limitations under the License.
17 *
18 **********************************************************************/
19
20 #include <string.h>
21 #include "otsuthr.h"
22
23 namespace tesseract {
24
25 // Compute the Otsu threshold(s) for the given image rectangle, making one
26 // for each channel. Each channel is always one byte per pixel.
27 // Returns an array of threshold values and an array of hi_values, such
28 // that a pixel value >threshold[channel] is considered foreground if
29 // hi_values[channel] is 0 or background if 1. A hi_value of -1 indicates
30 // that there is no apparent foreground. At least one hi_value will not be -1.
31 // Delete thresholds and hi_values with delete [] after use.
OtsuThreshold(const unsigned char * imagedata,int bytes_per_pixel,int bytes_per_line,int left,int top,int width,int height,int ** thresholds,int ** hi_values)32 void OtsuThreshold(const unsigned char* imagedata,
33 int bytes_per_pixel, int bytes_per_line,
34 int left, int top, int width, int height,
35 int** thresholds, int** hi_values) {
36 // Of all channels with no good hi_value, keep the best so we can always
37 // produce at least one answer.
38 int best_hi_value = 1;
39 int best_hi_index = 0;
40 bool any_good_hivalue = false;
41 double best_hi_dist = 0.0;
42 *thresholds = new int[bytes_per_pixel];
43 *hi_values = new int[bytes_per_pixel];
44
45 for (int ch = 0; ch < bytes_per_pixel; ++ch) {
46 (*thresholds)[ch] = -1;
47 (*hi_values)[ch] = -1;
48 // Compute the histogram of the image rectangle.
49 int histogram[kHistogramSize];
50 HistogramRect(imagedata + ch, bytes_per_pixel, bytes_per_line,
51 left, top, width, height, histogram);
52 int H;
53 int best_omega_0;
54 int best_t = OtsuStats(histogram, &H, &best_omega_0);
55 if (best_omega_0 == 0 || best_omega_0 == H) {
56 // This channel is empty.
57 continue;
58 }
59 // To be a convincing foreground we must have a small fraction of H
60 // or to be a convincing background we must have a large fraction of H.
61 // In between we assume this channel contains no thresholding information.
62 int hi_value = best_omega_0 < H * 0.5;
63 (*thresholds)[ch] = best_t;
64 if (best_omega_0 > H * 0.75) {
65 any_good_hivalue = true;
66 (*hi_values)[ch] = 0;
67 } else if (best_omega_0 < H * 0.25) {
68 any_good_hivalue = true;
69 (*hi_values)[ch] = 1;
70 } else {
71 // In case all channels are like this, keep the best of the bad lot.
72 double hi_dist = hi_value ? (H - best_omega_0) : best_omega_0;
73 if (hi_dist > best_hi_dist) {
74 best_hi_dist = hi_dist;
75 best_hi_value = hi_value;
76 best_hi_index = ch;
77 }
78 }
79 }
80 if (!any_good_hivalue) {
81 // Use the best of the ones that were not good enough.
82 (*hi_values)[best_hi_index] = best_hi_value;
83 }
84 }
85
86 // Compute the histogram for the given image rectangle, and the given
87 // channel. (Channel pointed to by imagedata.) Each channel is always
88 // one byte per pixel.
89 // Bytes per pixel is used to skip channels not being
90 // counted with this call in a multi-channel (pixel-major) image.
91 // Histogram is always a kHistogramSize(256) element array to count
92 // occurrences of each pixel value.
HistogramRect(const unsigned char * imagedata,int bytes_per_pixel,int bytes_per_line,int left,int top,int width,int height,int * histogram)93 void HistogramRect(const unsigned char* imagedata,
94 int bytes_per_pixel, int bytes_per_line,
95 int left, int top, int width, int height,
96 int* histogram) {
97 int bottom = top + height;
98 memset(histogram, 0, sizeof(*histogram) * kHistogramSize);
99 const unsigned char* pixels = imagedata +
100 top * bytes_per_line +
101 left * bytes_per_pixel;
102 for (int y = top; y < bottom; ++y) {
103 for (int x = 0; x < width; ++x) {
104 ++histogram[pixels[x * bytes_per_pixel]];
105 }
106 pixels += bytes_per_line;
107 }
108 }
109
110 // Compute the Otsu threshold(s) for the given histogram.
111 // Also returns H = total count in histogram, and
112 // omega0 = count of histogram below threshold.
OtsuStats(const int * histogram,int * H_out,int * omega0_out)113 int OtsuStats(const int* histogram, int* H_out, int* omega0_out) {
114 int H = 0;
115 double mu_T = 0.0;
116 for (int i = 0; i < kHistogramSize; ++i) {
117 H += histogram[i];
118 mu_T += i * histogram[i];
119 }
120
121 // Now maximize sig_sq_B over t.
122 // http://www.ctie.monash.edu.au/hargreave/Cornall_Terry_328.pdf
123 int best_t = -1;
124 int omega_0, omega_1;
125 int best_omega_0 = 0;
126 double best_sig_sq_B = 0.0;
127 double mu_0, mu_1, mu_t;
128 omega_0 = 0;
129 mu_t = 0.0;
130 for (int t = 0; t < kHistogramSize - 1; ++t) {
131 omega_0 += histogram[t];
132 mu_t += t * static_cast<double>(histogram[t]);
133 if (omega_0 == 0)
134 continue;
135 omega_1 = H - omega_0;
136 if (omega_1 == 0)
137 break;
138 mu_0 = mu_t / omega_0;
139 mu_1 = (mu_T - mu_t) / omega_1;
140 double sig_sq_B = mu_1 - mu_0;
141 sig_sq_B *= sig_sq_B * omega_0 * omega_1;
142 if (best_t < 0 || sig_sq_B > best_sig_sq_B) {
143 best_sig_sq_B = sig_sq_B;
144 best_t = t;
145 best_omega_0 = omega_0;
146 }
147 }
148 if (H_out != NULL) *H_out = H;
149 if (omega0_out != NULL) *omega0_out = best_omega_0;
150 return best_t;
151 }
152
153 } // namespace tesseract.
154