/ocr/ocrservice/jni/opticalflow/feature_detector.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: Andrew Harp
 18//
 19// Various feature detecting functions.
 20
 21#include <float.h>
 22
 23#include "utils.h"
 24#include "time_log.h"
 25
 26#include "image.h"
 27#include "feature_detector.h"
 28
 29// Threshold for pixels to be considered different.
 30#define FAST_DIFF_AMOUNT 10
 31
 32// How far from edge of frame to stop looking for FAST features.
 33#define FAST_BORDER_BUFFER 20
 34
 35// Minimum enforced distance between detected features.
 36// Default
 37#define MIN_FEATURE_DIST_NORMAL 24
 38
 39// Regions selected as "interesting" (aka have annotations) can have denser
 40// feature coverage.
 41#define MIN_FEATURE_DIST_INTEREST 12
 42
 43// How many FAST qualifying pixels must be connected to a pixel for it to be
 44// considered a candidate feature for Harris filtering.
 45#define MIN_NUM_CONNECTED 8
 46
 47// Size of the window to integrate over for Harris filtering.
 48// Compare to WINDOW_SIZE in optical_flow.h.
 49#define HARRIS_WINDOW_SIZE 2
 50
 51// Arbitrary parameter for how picky Harris filter is, the higher the more
 52// discriminating.
 53#define SENSITIVITY 0.2f
 54
 55namespace flow {
 56
 57void scoreFeatures(const Image<int32>& I_x, const Image<int32>& I_y,
 58                   const int32 num_candidates,
 59                   Point2D* const candidate_features) {
 60  // Score all the features
 61  for (int32 i = 0; i < num_candidates; ++i) {
 62    Point2D* const feature = candidate_features + i;
 63    feature->score = harrisFilter(I_x, I_y, feature->x, feature->y);
 64  }
 65}
 66
 67// Quicksorts detected features by score and then selects them such that
 68// they are separated by a minimum distance.
 69int32 sortAndSelect(const int32 num_candidates, const int32 max_features,
 70                    const Image<bool>& interest_map,
 71                    Point2D* const candidate_features,
 72                    Point2D* const final_features,
 73                    Image<uint8>* const best_feature_map) {
 74  qsort(candidate_features, num_candidates);
 75
 76#ifdef SANITY_CHECKS
 77  // Verify that the array got sorted.
 78  float32 last_score = -FLT_MAX;
 79  for (int32 i = 0; i < num_candidates; ++i) {
 80    const float32 curr_score = (candidate_features + i)->score;
 81
 82    // Scores should be monotonically increasing.
 83    CHECK(last_score <= curr_score,
 84          "Quicksort failure! %d: %.5f > %d: %.5f",
 85          i - 1, last_score, i, curr_score);
 86
 87    last_score = curr_score;
 88  }
 89#endif
 90
 91  best_feature_map->clear(false);
 92
 93  int32 num_features = 0;
 94
 95  for (int32 i = num_candidates - 1; i >= 0; --i) {
 96    const Point2D& candidate = candidate_features[i];
 97
 98    // Since features are sorted, the first 0 or less value means we can stop
 99    // looking.
100    if (candidate.score <= 0.0f) {
101      break;
102    }
103
104    // Lookup whether this feature is in an interest region.  If so, other
105    // features may appear closer to it than normal.
106    const int32 distance = interest_map.getPixel(candidate.x, candidate.y) ?
107        MIN_FEATURE_DIST_INTEREST : MIN_FEATURE_DIST_NORMAL;
108
109    if (markImage(candidate.x, candidate.y, distance, best_feature_map)) {
110      final_features[num_features] = candidate;
111      num_features++;
112
113      if (num_features >= max_features) {
114        break;
115      }
116    }
117  }
118
119  return num_features;
120}
121
122// Walks along the given circle checking for pixels above or below the center.
123// Returns a score, or 0 if the feature did not pass the criteria.
124//
125// Parameters:
126//  circle_perimeter: the circumference in pixels of the circle.
127//  threshold: the minimum number of contiguous pixels that must be above or
128//             below the center value.
129//  center_ptr: the location of the center pixel in memory
130//  offsets: the relative offsets from the center pixel of the edge pixels.
131inline int32 testCircle(const int32 circle_perimeter, const int32 threshold,
132                        const uint8* const center_ptr,
133                        const int32* offsets) {
134  // Get the actual value of the center pixel for easier reference later on.
135  const int32 center_value = static_cast<int32>(*center_ptr);
136
137  // Number of total pixels to check.  Have to wrap around some in case
138  // the contiguous section is split by the array edges.
139  const int32 num_total = circle_perimeter + threshold - 1;
140
141  int32 num_above = 0;
142  int32 above_diff = 0;
143
144  int32 num_below = 0;
145  int32 below_diff = 0;
146
147  // Used to tell when this is definitely not going to meet the threshold so we
148  // can early abort.
149  int32 minimum_by_now = threshold - num_total + 1;
150
151  // Go through every pixel along the perimeter of the circle, and then around
152  // again a little bit.
153  for (int32 i = 0; i < num_total; ++i) {
154    // This should be faster than mod.
155    const int32 perim_index = i < circle_perimeter ? i : i - circle_perimeter;
156
157    // This gets the value of the current pixel along the perimeter by using
158    // a precomputed offset.
159    const int32 curr_value =
160        static_cast<int32>(center_ptr[offsets[perim_index]]);
161
162    const int32 difference = curr_value - center_value;
163
164    if (difference > FAST_DIFF_AMOUNT) {
165      above_diff += difference;
166      ++num_above;
167
168      num_below = 0;
169      below_diff = 0;
170
171      if (num_above >= threshold) {
172        return above_diff;
173      }
174    } else if (difference < -FAST_DIFF_AMOUNT) {
175      below_diff += difference;
176      ++num_below;
177
178      num_above = 0;
179      above_diff = 0;
180
181      if (num_below >= threshold) {
182        return below_diff;
183      }
184    } else {
185      num_above = 0;
186      num_below = 0;
187      above_diff = 0;
188      below_diff = 0;
189    }
190
191    // See if there's any chance of making the threshold.
192    if (max(num_above, num_below) < minimum_by_now) {
193      // Didn't pass.
194      return 0;
195    }
196    ++minimum_by_now;
197  }
198
199  // Didn't pass.
200  return 0;
201}
202
203// Creates features in a regular grid, regardless of image contents.
204int32 seedFeatures(const Image<uint8>& frame,
205                   const int32 num_x, const int32 num_y,
206                   const float32 left, const float32 top,
207                   const float32 right, const float32 bottom,
208                   const int32 type, Point2D* const features) {
209  int32 num_features = 0;
210
211  const float32 step_x = ((right - left) / (num_x - 1));
212  const float32 step_y = ((bottom - top) / (num_y - 1));
213
214  for (int32 x = 0; x < num_x; ++x) {
215    for (int32 y = 0; y < num_y; ++y) {
216      const int32 x_pos = x * step_x + left;
217      const int32 y_pos = y * step_y + top;
218      if (inRange(x_pos, 0, frame.width_less_one_) &&
219          inRange(y_pos, 0, frame.height_less_one_)) {
220        Point2D* const feature = features + num_features;
221        feature->x = x_pos;
222        feature->y = y_pos;
223        feature->type = type;
224
225        ++num_features;
226      }
227    }
228  }
229
230  return num_features;
231}
232
233
234// Returns how likely a point in the image is to be a corner.
235float32 harrisFilter(const Image<int32>& I_x, const Image<int32>& I_y,
236                     const int32 x, const int32 y) {
237  // Image gradient matrix.
238  float32 G[] = { 0, 0, 0, 0 };
239  calculateG(HARRIS_WINDOW_SIZE, x, y, I_x, I_y, G);
240
241  const float32 g_sum = G[0] + G[1] + G[2] + G[3];
242
243  const float32 a = G[0] / g_sum;
244  const float32 b = G[1] / g_sum;
245  const float32 c = G[2] / g_sum;
246  const float32 d = G[3] / g_sum;
247
248  const float32 det = a * d - b * c;
249  const float32 trace = a + d;
250
251  const float32 inner = square(trace) - 4 * det;
252
253  if (inner >= 0.0f) {
254    const float32 square_root_inner = sqrtf(inner);
255    const float32 eig1 = (trace + square_root_inner) / 2.0f;
256    const float32 eig2 = (trace - square_root_inner) / 2.0f;
257    return eig1 * eig2 - SENSITIVITY * square(eig1 + eig2);
258  }
259
260  // Way negative.
261  return -100.0f;
262}
263
264// FAST feature detector.
265int32 findFastFeatures(const Image<uint8>& frame, const int32 max_num_features,
266                       Point2D* const features,
267                       Image<uint8>* const best_feature_map) {
268  /*
269   // Reference for a circle of diameter 7.
270   const int32 circle[] = {0, 0, 1, 1, 1, 0, 0,
271                           0, 1, 0, 0, 0, 1, 0,
272                           1, 0, 0, 0, 0, 0, 1,
273                           1, 0, 0, 0, 0, 0, 1,
274                           1, 0, 0, 0, 0, 0, 1,
275                           0, 1, 0, 0, 0, 1, 0,
276                           0, 0, 1, 1, 1, 0, 0};
277   const int32 circle_offset[] =
278       {2, 3, 4, 8, 12, 14, 20, 21, 27, 28, 34, 36, 40, 44, 45, 46};
279   */
280
281  // Quick test of compass directions.  Any length 16 circle with a break of up
282  // to 4 pixels will have at least 3 of these 4 pixels active.
283  static const int32 short_circle_perimeter = 4;
284  static const int32 short_threshold = 3;
285  static const int32 short_circle_x[] = { -3,  0, +3,  0 };
286  static const int32 short_circle_y[] = {  0, -3,  0, +3 };
287
288  // Precompute image offsets.
289  int32 short_offsets[short_circle_perimeter];
290  for (int i = 0; i < short_circle_perimeter; ++i) {
291    short_offsets[i] = short_circle_x[i] + short_circle_y[i] * frame.getWidth();
292  }
293
294  // Large circle values.
295  static const int32 full_circle_perimeter = 16;
296  static const int32 full_threshold = 12;
297  static const int32 full_circle_x[] =
298      { -1,  0, +1, +2, +3, +3, +3, +2, +1, +0, -1, -2, -3, -3, -3, -2 };
299  static const int32 full_circle_y[] =
300      { -3, -3, -3, -2, -1,  0, +1, +2, +3, +3, +3, +2, +1, +0, -1, -2 };
301
302  // Precompute image offsets.
303  int32 full_offsets[full_circle_perimeter];
304  for (int i = 0; i < full_circle_perimeter; ++i) {
305    full_offsets[i] = full_circle_x[i] + full_circle_y[i] * frame.getWidth();
306  }
307
308  const int frame_width = frame.getWidth();
309
310  const int end_y = frame.getHeight() - FAST_BORDER_BUFFER;
311  const int end_x = frame.getWidth() - FAST_BORDER_BUFFER;
312
313  best_feature_map->clear(0);
314
315  // Loop through once to find FAST feature clumps.
316  for (int32 img_y = FAST_BORDER_BUFFER; img_y < end_y; ++img_y) {
317    const uint8* curr_pixel_ptr =
318        frame.getPixelPtrConst(FAST_BORDER_BUFFER, img_y);
319
320    for (int32 img_x = FAST_BORDER_BUFFER; img_x < end_x; ++img_x) {
321      // Only insert it if it meets the quick minimum requirements test.
322      if (testCircle(short_circle_perimeter, short_threshold,
323                     curr_pixel_ptr, short_offsets) != 0) {
324
325        // Longer test for actual feature score..
326        const int32 fast_score = testCircle(full_circle_perimeter,
327                                            full_threshold,
328                                            curr_pixel_ptr,
329                                            full_offsets);
330
331        // Non-zero score means the feature was found.
332        if (fast_score != 0) {
333          uint8* const center_ptr = best_feature_map->getPixelPtr(img_x, img_y);
334
335          // Increase the feature count on this pixel and the pixels in all
336          // 4 cardinal directions.
337          *center_ptr += 5;
338          *(center_ptr - 1) += 1;
339          *(center_ptr + 1) += 1;
340          *(center_ptr - frame_width) += 1;
341          *(center_ptr + frame_width) += 1;
342        }
343      }
344
345      ++curr_pixel_ptr;
346    }  // x
347  }  // y
348
349  timeLog("Found FAST features");
350
351  int32 num_features = 0;
352  // Loop through again and Harris filter pixels in the center of clumps.
353  // We can shrink the window by 1 pixel on every side.
354  for (int32 img_y = FAST_BORDER_BUFFER + 1; img_y < end_y - 1; ++img_y) {
355    const int32 start_x = FAST_BORDER_BUFFER + 1;
356
357    const uint8* curr_pixel_ptr =
358        best_feature_map->getPixelPtrConst(start_x, img_y);
359
360    for (int32 img_x = start_x; img_x < end_x - 1; ++img_x) {
361      if (*curr_pixel_ptr >= MIN_NUM_CONNECTED) {
362        Point2D* const feature = features + num_features;
363        feature->x = img_x;
364        feature->y = img_y;
365        feature->score = 0;
366        feature->type = FEATURE_FAST;
367
368        ++num_features;
369        if (num_features >= max_num_features) {
370          return num_features;
371        }
372      }
373
374      ++curr_pixel_ptr;
375    }  // x
376  }  // y
377
378  timeLog("Filtered FAST features");
379  return num_features;
380}
381
382}  // namespace flow