tesseract  4.0.0-1-g2a2b
otsuthr.cpp
Go to the documentation of this file.
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 "otsuthr.h"
21 
22 #include <cstring>
23 #include "allheaders.h"
24 #include "helpers.h"
25 #include "openclwrapper.h"
26 
27 
28 namespace tesseract {
29 
30 // Computes the Otsu threshold(s) for the given image rectangle, making one
31 // for each channel. Each channel is always one byte per pixel.
32 // Returns an array of threshold values and an array of hi_values, such
33 // that a pixel value >threshold[channel] is considered foreground if
34 // hi_values[channel] is 0 or background if 1. A hi_value of -1 indicates
35 // that there is no apparent foreground. At least one hi_value will not be -1.
36 // Delete thresholds and hi_values with delete [] after use.
37 // The return value is the number of channels in the input image, being
38 // the size of the output thresholds and hi_values arrays.
39 int OtsuThreshold(Pix* src_pix, int left, int top, int width, int height,
40  int** thresholds, int** hi_values) {
41  int num_channels = pixGetDepth(src_pix) / 8;
42  // Of all channels with no good hi_value, keep the best so we can always
43  // produce at least one answer.
44  PERF_COUNT_START("OtsuThreshold")
45  int best_hi_value = 1;
46  int best_hi_index = 0;
47  bool any_good_hivalue = false;
48  double best_hi_dist = 0.0;
49  *thresholds = new int[num_channels];
50  *hi_values = new int[num_channels];
51 
52  // only use opencl if compiled w/ OpenCL and selected device is opencl
53 #ifdef USE_OPENCL
54  // all of channel 0 then all of channel 1...
55  int* histogramAllChannels = new int[kHistogramSize * num_channels];
56 
57  // Calculate Histogram on GPU
58  OpenclDevice od;
59  if (od.selectedDeviceIsOpenCL() && (num_channels == 1 || num_channels == 4) &&
60  top == 0 && left == 0) {
61  od.HistogramRectOCL(pixGetData(src_pix), num_channels,
62  pixGetWpl(src_pix) * 4, left, top, width, height,
63  kHistogramSize, histogramAllChannels);
64 
65  // Calculate Threshold from Histogram on cpu
66  for (int ch = 0; ch < num_channels; ++ch) {
67  (*thresholds)[ch] = -1;
68  (*hi_values)[ch] = -1;
69  int *histogram = &histogramAllChannels[kHistogramSize * ch];
70  int H;
71  int best_omega_0;
72  int best_t = OtsuStats(histogram, &H, &best_omega_0);
73  if (best_omega_0 == 0 || best_omega_0 == H) {
74  // This channel is empty.
75  continue;
76  }
77  // To be a convincing foreground we must have a small fraction of H
78  // or to be a convincing background we must have a large fraction of H.
79  // In between we assume this channel contains no thresholding information.
80  int hi_value = best_omega_0 < H * 0.5;
81  (*thresholds)[ch] = best_t;
82  if (best_omega_0 > H * 0.75) {
83  any_good_hivalue = true;
84  (*hi_values)[ch] = 0;
85  } else if (best_omega_0 < H * 0.25) {
86  any_good_hivalue = true;
87  (*hi_values)[ch] = 1;
88  } else {
89  // In case all channels are like this, keep the best of the bad lot.
90  double hi_dist = hi_value ? (H - best_omega_0) : best_omega_0;
91  if (hi_dist > best_hi_dist) {
92  best_hi_dist = hi_dist;
93  best_hi_value = hi_value;
94  best_hi_index = ch;
95  }
96  }
97  }
98  } else {
99 #endif
100  for (int ch = 0; ch < num_channels; ++ch) {
101  (*thresholds)[ch] = -1;
102  (*hi_values)[ch] = -1;
103  // Compute the histogram of the image rectangle.
104  int histogram[kHistogramSize];
105  HistogramRect(src_pix, ch, left, top, width, height, histogram);
106  int H;
107  int best_omega_0;
108  int best_t = OtsuStats(histogram, &H, &best_omega_0);
109  if (best_omega_0 == 0 || best_omega_0 == H) {
110  // This channel is empty.
111  continue;
112  }
113  // To be a convincing foreground we must have a small fraction of H
114  // or to be a convincing background we must have a large fraction of H.
115  // In between we assume this channel contains no thresholding information.
116  int hi_value = best_omega_0 < H * 0.5;
117  (*thresholds)[ch] = best_t;
118  if (best_omega_0 > H * 0.75) {
119  any_good_hivalue = true;
120  (*hi_values)[ch] = 0;
121  } else if (best_omega_0 < H * 0.25) {
122  any_good_hivalue = true;
123  (*hi_values)[ch] = 1;
124  } else {
125  // In case all channels are like this, keep the best of the bad lot.
126  double hi_dist = hi_value ? (H - best_omega_0) : best_omega_0;
127  if (hi_dist > best_hi_dist) {
128  best_hi_dist = hi_dist;
129  best_hi_value = hi_value;
130  best_hi_index = ch;
131  }
132  }
133  }
134 #ifdef USE_OPENCL
135  }
136  delete[] histogramAllChannels;
137 #endif // USE_OPENCL
138 
139  if (!any_good_hivalue) {
140  // Use the best of the ones that were not good enough.
141  (*hi_values)[best_hi_index] = best_hi_value;
142  }
144  return num_channels;
145 }
146 
147 // Computes the histogram for the given image rectangle, and the given
148 // single channel. Each channel is always one byte per pixel.
149 // Histogram is always a kHistogramSize(256) element array to count
150 // occurrences of each pixel value.
151 void HistogramRect(Pix* src_pix, int channel,
152  int left, int top, int width, int height,
153  int* histogram) {
154  PERF_COUNT_START("HistogramRect")
155  int num_channels = pixGetDepth(src_pix) / 8;
156  channel = ClipToRange(channel, 0, num_channels - 1);
157  int bottom = top + height;
158  memset(histogram, 0, sizeof(*histogram) * kHistogramSize);
159  int src_wpl = pixGetWpl(src_pix);
160  l_uint32* srcdata = pixGetData(src_pix);
161  for (int y = top; y < bottom; ++y) {
162  const l_uint32* linedata = srcdata + y * src_wpl;
163  for (int x = 0; x < width; ++x) {
164  int pixel = GET_DATA_BYTE(linedata, (x + left) * num_channels + channel);
165  ++histogram[pixel];
166  }
167  }
169 }
170 
171 // Computes the Otsu threshold(s) for the given histogram.
172 // Also returns H = total count in histogram, and
173 // omega0 = count of histogram below threshold.
174 int OtsuStats(const int* histogram, int* H_out, int* omega0_out) {
175  int H = 0;
176  double mu_T = 0.0;
177  for (int i = 0; i < kHistogramSize; ++i) {
178  H += histogram[i];
179  mu_T += static_cast<double>(i) * histogram[i];
180  }
181 
182  // Now maximize sig_sq_B over t.
183  // http://www.ctie.monash.edu.au/hargreave/Cornall_Terry_328.pdf
184  int best_t = -1;
185  int omega_0, omega_1;
186  int best_omega_0 = 0;
187  double best_sig_sq_B = 0.0;
188  double mu_0, mu_1, mu_t;
189  omega_0 = 0;
190  mu_t = 0.0;
191  for (int t = 0; t < kHistogramSize - 1; ++t) {
192  omega_0 += histogram[t];
193  mu_t += t * static_cast<double>(histogram[t]);
194  if (omega_0 == 0)
195  continue;
196  omega_1 = H - omega_0;
197  if (omega_1 == 0)
198  break;
199  mu_0 = mu_t / omega_0;
200  mu_1 = (mu_T - mu_t) / omega_1;
201  double sig_sq_B = mu_1 - mu_0;
202  sig_sq_B *= sig_sq_B * omega_0 * omega_1;
203  if (best_t < 0 || sig_sq_B > best_sig_sq_B) {
204  best_sig_sq_B = sig_sq_B;
205  best_t = t;
206  best_omega_0 = omega_0;
207  }
208  }
209  if (H_out != nullptr) *H_out = H;
210  if (omega0_out != nullptr) *omega0_out = best_omega_0;
211  return best_t;
212 }
213 
214 } // namespace tesseract.
int OtsuStats(const int *histogram, int *H_out, int *omega0_out)
Definition: otsuthr.cpp:174
#define PERF_COUNT_START(FUNCT_NAME)
const int kHistogramSize
Definition: otsuthr.h:27
void HistogramRect(Pix *src_pix, int channel, int left, int top, int width, int height, int *histogram)
Definition: otsuthr.cpp:151
#define PERF_COUNT_END
T ClipToRange(const T &x, const T &lower_bound, const T &upper_bound)
Definition: helpers.h:111
int OtsuThreshold(Pix *src_pix, int left, int top, int width, int height, int **thresholds, int **hi_values)
Definition: otsuthr.cpp:39