/ocr/ocrservice/jni/imageutils/blur.cpp

  1/*
  2 * Copyright 2011, Google Inc.
  3 *
  4 * Licensed under the Apache License, Version 2.0 (the "License");
  5 * you may not use this file except in compliance with the License.
  6 * You may obtain a copy of the License at
  7 *
  8 *     http://www.apache.org/licenses/LICENSE-2.0
  9 *
 10 * Unless required by applicable law or agreed to in writing, software
 11 * distributed under the License is distributed on an "AS IS" BASIS,
 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 13 * See the License for the specific language governing permissions and
 14 * limitations under the License.
 15 */
 16
 17// Author: Xiaotao Duan
 18//
 19// This library contains image processing method to detect
 20// image blurriness.
 21//
 22// This library is *not* thread safe because static memory is
 23// used for performance.
 24//
 25// A method to detect whether a given image is blurred or not.
 26// The algorithm is based on H. Tong, M. Li, H. Zhang, J. He,
 27// and C. Zhang. "Blur detection for digital images using wavelet
 28// transform".
 29//
 30// To achieve better performance on client side, the method
 31// is running on four 128x128 portions which compose the 256x256
 32// central area of the given image. On Nexus One, average time
 33// to process a single image is ~5 milliseconds.
 34
 35#include <math.h>
 36
 37#include "blur.h"
 38#include "utils.h"
 39
 40static const int kDecomposition = 3;
 41static const int kThreshold = 35;
 42static const float kMinZero = 0.05;
 43
 44static const int kMaximumWidth = 256;
 45static const int kMaximumHeight = 256;
 46
 47static int32 _smatrix[kMaximumWidth * kMaximumHeight];
 48static int32 _arow[kMaximumWidth > kMaximumHeight ?
 49    kMaximumWidth : kMaximumHeight];
 50
 51// Does Haar Wavelet Transformation in place on a given row of a matrix.
 52// The matrix is in size of matrix_height * matrix_width and represented
 53// in a linear array. Parameter offset_row indicates transformation is
 54// performed on which row. offset_column and num_columns indicate column
 55// range of the given row.
 56inline void Haar1DX(int* matrix, int matrix_height, int matrix_width,
 57    int offset_row, int offset_column, int num_columns) {
 58  int32* ptr_a = _arow;
 59  int32* ptr_matrix = matrix + offset_row * matrix_width + offset_column;
 60  int half_num_columns = num_columns / 2;
 61
 62  int32* a_tmp = ptr_a;
 63  int32* matrix_tmp = ptr_matrix;
 64  for (int j = 0; j < half_num_columns; ++j) {
 65    *a_tmp++ = (matrix_tmp[0] + matrix_tmp[1]) / 2;
 66    matrix_tmp += 2;
 67  }
 68
 69  int32* average = ptr_a;
 70  a_tmp = ptr_a + half_num_columns;
 71  matrix_tmp = ptr_matrix;
 72  for (int j = 0; j < half_num_columns; ++j) {
 73    *a_tmp++ = *matrix_tmp - *average++;
 74    matrix_tmp += 2;
 75  }
 76
 77  memcpy(ptr_matrix, ptr_a, sizeof(int32) * num_columns);
 78}
 79
 80// Does Haar Wavelet Transformation in place on a given column of a matrix.
 81inline void Haar1DY(int* matrix, int matrix_height, int matrix_width,
 82    int offset_column, int offset_row, int num_rows) {
 83  int32* ptr_a = _arow;
 84  int32* ptr_matrix = matrix + offset_row * matrix_width + offset_column;
 85  int half_num_rows = num_rows / 2;
 86  int two_line_width = matrix_width * 2;
 87
 88  int32* a_tmp = ptr_a;
 89  int32* matrix_tmp = ptr_matrix;
 90  for (int j = 0; j < half_num_rows; ++j) {
 91    *a_tmp++ = (matrix_tmp[matrix_width] + matrix_tmp[0]) / 2;
 92    matrix_tmp += two_line_width;
 93  }
 94
 95  int32* average = ptr_a;
 96  a_tmp = ptr_a + half_num_rows;
 97  matrix_tmp = ptr_matrix;
 98  for (int j = 0; j < num_rows; j += 2) {
 99    *a_tmp++ = *matrix_tmp - *average++;
100    matrix_tmp += two_line_width;
101  }
102
103  for (int j = 0; j < num_rows; ++j) {
104    *ptr_matrix = *ptr_a++;
105    ptr_matrix += matrix_width;
106  }
107}
108
109// Does Haar Wavelet Transformation in place for a specified area of
110// a matrix. The matrix size is specified by matrix_width and matrix_height.
111// The area on which the transformation is performed is specified by
112// offset_column, num_columns, offset_row and num_rows.
113void Haar2D(int* matrix, int matrix_height, int matrix_width,
114    int offset_column, int num_columns, int offset_row, int num_rows) {
115  for (int i = offset_row; i < offset_row + num_rows; ++i) {
116    Haar1DX(matrix, matrix_height, matrix_width, i, offset_column, num_columns);
117  }
118
119  for (int i = offset_column; i < offset_column + num_columns; ++i){
120    Haar1DY(matrix, matrix_height, matrix_width, i, offset_row, num_rows);
121  }
122}
123
124// Reads in a given matrix, does first round HWT and outputs result
125// matrix into target array. This function is used for optimization by
126// avoiding a memory copy. The input matrix has height rows and width
127// columns. The transformation is performed on the given area specified
128// by offset_column, num_columns, offset_row, num_rows. After
129// transformation, the output matrix has num_columns columns and
130// num_rows rows.
131void HwtFirstRound(const uint8* const data, int height, int width,
132    int offset_column, int num_columns,
133    int offset_row, int num_rows, int32* matrix) {
134  int32* ptr_a = _arow;
135  const uint8* ptr_data = data + offset_row * width + offset_column;
136  int half_num_columns = num_columns / 2;
137  for (int i = 0; i < num_rows; ++i) {
138    int32* a_tmp = ptr_a;
139    const uint8* data_tmp = ptr_data;
140    for (int j = 0; j < half_num_columns; ++j) {
141      *a_tmp++ = (int32) ((data_tmp[0] + data_tmp[1]) / 2);
142      data_tmp += 2;
143    }
144
145    int32* average = ptr_a;
146    a_tmp = ptr_a + half_num_columns;
147    data_tmp = ptr_data;
148    for (int j = 0; j < half_num_columns; ++j) {
149      *a_tmp++ = *data_tmp - *average++;
150      data_tmp += 2;
151    }
152
153    int32* ptr_matrix = matrix + i * num_columns;
154    a_tmp = ptr_a;
155    for (int j = 0; j < num_columns; ++j) {
156      *ptr_matrix++ = *a_tmp++;
157    }
158
159    ptr_data += width;
160  }
161
162  // Column transformation does not involve input data.
163  for (int i = 0; i < num_columns; ++i) {
164    Haar1DY(matrix, num_rows, num_columns, i, 0, num_rows);
165  }
166}
167
168// Returns the weight of a given point in a certain scale of a matrix
169// after wavelet transformation.
170// The point is specified by k and l which are y and x coordinate
171// respectively. Parameter scale tells in which scale the weight is
172// computed, must be 1, 2 or 3 which stands respectively for 1/2, 1/4
173// and 1/8 of original size.
174int ComputeEdgePointWeight(int* matrix, int width, int height,
175    int k, int l, int scale) {
176  int r = k >> scale;
177  int c = l >> scale;
178  int window_row = height >> scale;
179  int window_column = width >> scale;
180
181  int v_top_right = square(matrix[r * width + c + window_column]);
182  int v_bot_left = square(matrix[(r + window_row) * width + c]);
183  int v_bot_right =
184      square(matrix[(r + window_row) * width + c + window_column]);
185
186  int v = sqrt(v_top_right + v_bot_left + v_bot_right);
187  return v;
188}
189
190// Computes point with maximum weight for a given local window for a
191// given scale.
192// Parameter scaled_width and scaled_height define scaled image size
193// of a certain decomposition level. The window size is defined by
194// window_size. Output value k and l store row (y coordinate) and
195// column (x coordinate) respectively of the point with maximum weight.
196// The maximum weight is returned.
197int ComputeLocalMaximum(int* matrix, int width, int height,
198    int scaled_width, int scaled_height,
199    int top, int left, int window_size, int* k, int* l) {
200  int max = -1;
201  *k = top;
202  *l = left;
203
204  for (int i = 0; i < window_size; ++i) {
205    for (int j = 0; j < window_size; ++j) {
206      int r = top + i;
207      int c = left + j;
208
209      int v_top_right = abs(matrix[r * width + c + scaled_width]);
210      int v_bot_left = abs(matrix[(r + scaled_height) * width + c]);
211      int v_bot_right =
212          abs(matrix[(r + scaled_height) * width + c + scaled_width]);
213      int v = v_top_right + v_bot_left + v_bot_right;
214
215      if (v > max) {
216        max = v;
217        *k = r;
218        *l = c;
219      }
220    }
221  }
222
223  int r = *k;
224  int c = *l;
225  int v_top_right = square(matrix[r * width + c + scaled_width]);
226  int v_bot_left = square(matrix[(r + scaled_height) * width + c]);
227  int v_bot_right =
228      square(matrix[(r + scaled_height) * width + c + scaled_width]);
229  int v = sqrt(v_top_right + v_bot_left + v_bot_right);
230
231  return v;
232}
233
234// Detects blurriness of a transformed matrix.
235// Blur confidence and extent will be returned through blur_conf
236// and blur_extent. 1 is returned while input matrix is blurred.
237int DetectBlur(int* matrix, int width, int height,
238    float* blur_conf, float* blur_extent) {
239  int nedge = 0;
240  int nda = 0;
241  int nrg = 0;
242  int nbrg = 0;
243
244  // For each scale
245  for (int current_scale = kDecomposition; current_scale > 0; --current_scale) {
246    int scaled_width = width >> current_scale;
247    int scaled_height = height >> current_scale;
248    int window_size = 16 >> current_scale;  // 2, 4, 8
249    // For each window
250    for (int r = 0; r + window_size < scaled_height; r += window_size) {
251      for (int c = 0; c + window_size < scaled_width; c += window_size) {
252        int k, l;
253        int emax = ComputeLocalMaximum(matrix, width, height,
254            scaled_width, scaled_height, r, c, window_size, &k, &l);
255        if (emax > kThreshold) {
256          int emax1, emax2, emax3;
257          switch (current_scale) {
258            case 1:
259              emax1 = emax;
260              emax2 = ComputeEdgePointWeight(matrix, width, height,
261                  k << current_scale, l << current_scale, 2);
262              emax3 = ComputeEdgePointWeight(matrix, width, height,
263                  k << current_scale, l << current_scale, 3);
264              break;
265            case 2:
266              emax1 = ComputeEdgePointWeight(matrix, width, height,
267                  k << current_scale, l << current_scale, 1);
268              emax2 = emax;
269              emax3 = ComputeEdgePointWeight(matrix, width, height,
270                  k << current_scale, l << current_scale, 3);
271              break;
272            case 3:
273              emax1 = ComputeEdgePointWeight(matrix, width, height,
274                  k << current_scale, l << current_scale, 1);
275              emax2 = ComputeEdgePointWeight(matrix, width, height,
276                  k << current_scale, l << current_scale, 2);
277              emax3 = emax;
278              break;
279          }
280
281          nedge++;
282          if (emax1 > emax2 && emax2 > emax3) {
283            nda++;
284          }
285          if (emax1 < emax2 && emax2 < emax3) {
286            nrg++;
287            if (emax1 < kThreshold) {
288              nbrg++;
289            }
290          }
291          if (emax2 > emax1 && emax2 > emax3) {
292            nrg++;
293            if (emax1 < kThreshold) {
294              nbrg++;
295            }
296          }
297        }
298      }
299    }
300  }
301
302  // TODO(xiaotao): No edge point at all, blurred or not?
303  float per = nedge == 0 ? 0 : (float)nda / nedge;
304
305  *blur_conf = per;
306  *blur_extent = (float)nbrg / nrg;
307
308  return per < kMinZero;
309}
310
311// Detects blurriness of a given portion of a luminance matrix.
312int IsBlurredInner(const uint8* const luminance,
313    const int width, const int height,
314    const int left, const int top,
315    const int width_wanted, const int height_wanted,
316    float* const blur, float* const extent) {
317  int32* matrix = _smatrix;
318
319  HwtFirstRound(luminance, height, width,
320                left, width_wanted, top, height_wanted, matrix);
321  Haar2D(matrix, height_wanted, width_wanted,
322         0, width_wanted >> 1, 0, height_wanted >> 1);
323  Haar2D(matrix, height_wanted, width_wanted,
324         0, width_wanted >> 2, 0, height_wanted >> 2);
325
326  int blurred = DetectBlur(matrix, width_wanted, height_wanted, blur, extent);
327
328  return blurred;
329}
330
331int IsBlurred(const uint8* const luminance,
332    const int width, const int height, float* const blur, float* const extent) {
333
334  int desired_width = min(kMaximumWidth, width);
335  int desired_height = min(kMaximumHeight, height);
336  int left = (width - desired_width) >> 1;
337  int top = (height - desired_height) >> 1;
338
339  float conf1, extent1;
340  int blur1 = IsBlurredInner(luminance, width, height,
341      left, top, desired_width >> 1, desired_height >> 1, &conf1, &extent1);
342  float conf2, extent2;
343  int blur2 = IsBlurredInner(luminance, width, height,
344      left + (desired_width >> 1), top, desired_width >> 1, desired_height >> 1,
345      &conf2, &extent2);
346  float conf3, extent3;
347  int blur3 = IsBlurredInner(luminance, width, height,
348      left, top + (desired_height >> 1), desired_width >> 1,
349      desired_height >> 1, &conf3, &extent3);
350  float conf4, extent4;
351  int blur4 = IsBlurredInner(luminance, width, height,
352      left + (desired_width >> 1), top + (desired_height >> 1),
353      desired_width >> 1, desired_height >> 1, &conf4, &extent4);
354
355  *blur = (conf1 + conf2 + conf3 + conf4) / 4;
356  *extent = (extent1 + extent2 + extent3 + extent4) / 4;
357  return *blur < kMinZero;
358}