Fawkes API  Fawkes Development Version
rht_circle.cpp
1 
2 /***************************************************************************
3  * rht_circle.cpp - Implementation of a circle shape finder
4  * with Randomized Hough Transform
5  *
6  * Created: Tue Jun 28 00:00:00 2005
7  * Copyright 2005 Hu Yuxiao <Yuxiao.Hu@rwth-aachen.de>
8  * Tim Niemueller [www.niemueller.de]
9  *
10  ****************************************************************************/
11 
12 /* This program is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version. A runtime exception applies to
16  * this software (see LICENSE.GPL_WRE file mentioned below for details).
17  *
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU Library General Public License for more details.
22  *
23  * Read the full text in the LICENSE.GPL_WRE file in the doc directory.
24  */
25 
26 #include <fvmodels/shape/rht_circle.h>
27 
28 #include <cmath>
29 #include <sys/time.h>
30 
31 using namespace std;
32 using namespace fawkes;
33 
34 #define TEST_IF_IS_A_PIXEL(x) ((x)>230)
35 
36 #define TBY_SQUARED_DIST(x1,y1,x2,y2) \
37  (((x1)-(x2))*((x1)-(x2))+((y1)-(y2))*((y1)-(y2)))
38 #define TBY_RADIUS_DIFF(x1, y1, x2, y2, r) \
39  (sqrt(((x1)-(x2))*((x1)-(x2))+((y1)-(y2))*((y1)-(y2))-(r)*(r)))
40 
41 
42 namespace firevision {
43 #if 0 /* just to make Emacs auto-indent happy */
44 }
45 #endif
46 
47 const float RhtCircleModel::RHT_MIN_RADIUS = 40.0;
48 const float RhtCircleModel::RHT_MAX_RADIUS = 500.0;
49 
50 /** @class RhtCircleModel <fvmodels/shape/rht_circle.h>
51  * Randomized Hough-Transform circle model.
52  */
53 
54 /** Constructor. */
55 RhtCircleModel::RhtCircleModel(void)
56 {
57 }
58 
59 
60 /** Destructor. */
61 RhtCircleModel::~RhtCircleModel(void)
62 {
63  m_Circles.clear();
64 }
65 
66 /**************************************************************
67  * In this function I implement the circle detection algorithm
68  * from the following literature
69  * Randomized Hough Transform
70  **************************************************************/
71 int
72 RhtCircleModel::parseImage( unsigned char* buf,
73  ROI *roi )
74 {
75 
76  unsigned char *buffer = roi->get_roi_buffer_start(buf);
77 
78  unsigned int x, y;
79  vector<upoint_t> pixels;
80  struct timeval start, end;
81 
82  // clear the accumulator
83  accumulator.reset();
84 
85  // clear all the remembered circles
86  m_Circles.clear();
87 
88  // The following constants are used as stopping criteria
89  const int RHT_MAX_TIME = 10000; // = 10ms (is given in microseconds)
90  const int RHT_MAX_ITER = 1000; // Maximal number of iterations.
91 
92  // The following constant is used for eliminating circles with too few votes
93  const float RHT_MIN_VOTE_RATE = 0.0f;
94 
95  // The following constants are used for RHT accumulator precision
96  const int RHT_XY_SCALE = 8;
97  const int RHT_RADIUS_SCALE= 8;
98 
99  // All the pixels with a less distance difference than below
100  // are taken into account for circle fitting.
101  const float RHT_FITTING_DIST_DIFF = 15.0f;
102 
103  // The following constant is used for the size of the hollow window in the ROI.
104  const float ROI_HOLLOW_RATE = 0.70f;
105 
106  const unsigned int roi_hollow_top = (int)(roi->height * ((1.0f - ROI_HOLLOW_RATE) / 2));
107  const unsigned int roi_hollow_bottom = roi->height - roi_hollow_top;
108  const unsigned int roi_hollow_left = (int)(roi->width * ((1.0f - ROI_HOLLOW_RATE) / 2));
109  const unsigned int roi_hollow_right = roi->width - roi_hollow_left;
110 
111  // First, find all the pixels on the edges,
112  // and store them in the 'pixels' vector.
113  // NEW: excluding the hollow window
114  unsigned char *line_start = buffer;
115 
116  gettimeofday(&start, NULL);
117  end.tv_usec = start.tv_usec;
118 
119  // top "1/3"
120  for (y = 0; y < roi_hollow_top; ++y) {
121  for (x = 0; x < roi->width; ++x) {
122  if (TEST_IF_IS_A_PIXEL(*buffer)) {
123  upoint_t pt={x, y};
124  pixels.push_back(pt);
125  }
126  // NOTE: this assumes roi->pixel_step == 1
127  ++buffer;
128  }
129  line_start += roi->line_step;
130  buffer = line_start;
131  }
132  // middle "1/3"
133  for (y = roi_hollow_top; y < roi_hollow_bottom; ++y) {
134  for (x = 0; x < roi_hollow_left; ++x) {
135  if (TEST_IF_IS_A_PIXEL(*buffer)) {
136  upoint_t pt={x, y};
137  pixels.push_back(pt);
138  }
139  // NOTE: this assumes roi->pixel_step == 1
140  ++buffer;
141  }
142  buffer+=(roi_hollow_right - roi_hollow_left);
143  for (x = roi_hollow_right; x < roi->width; ++x) {
144  if (TEST_IF_IS_A_PIXEL(*buffer)) {
145  upoint_t pt={x, y};
146  pixels.push_back(pt);
147  }
148  // NOTE: this assumes roi->pixel_step == 1
149  ++buffer;
150  }
151  line_start += roi->line_step;
152  buffer = line_start;
153  }
154  // bottom "1/3"
155  for (y = roi_hollow_bottom; y < roi->height; ++y) {
156  for (x = 0; x < roi->width; ++x) {
157  if (TEST_IF_IS_A_PIXEL(*buffer)) {
158  upoint_t pt={x, y};
159  pixels.push_back(pt);
160  }
161  // NOTE: this assumes roi->pixel_step == 1
162  ++buffer;
163  }
164  line_start += roi->line_step;
165  buffer = line_start;
166  }
167 
168  // Then perform the RHT algorithm
169  upoint_t p[3];
170  center_in_roi_t center;
171  float radius;
172  vector< upoint_t >::iterator pos;
173  int num_iter = 0;
174  int num_points = (int) pixels.size();
175  if (num_points == 0) {
176  // No pixels found => no edge => no circle
177  return 0;
178  }
179 
180  while( (num_iter++ < RHT_MAX_ITER) &&
181  ( ((end.tv_usec - start.tv_usec) < RHT_MAX_TIME) ||
182  ((end.tv_usec + 1000000 - start.tv_usec) < RHT_MAX_TIME) )
183  // this only works if time constraint smaller than 500ms
184  ) {
185 
186  // Pick three points, and move them to the remove_list.
187  for (int i=0; i < 3; ++i) {
188  int ri = rand() % num_points;
189  pos = pixels.begin() + ri;
190  p[i] = *pos; // use * operator of iterator
191  }
192 
193  // Now calculate the center and radius
194  // based on the three points.
195  calcCircle(p[0], p[1], p[2], center, radius);
196 
197  // Accumulate this circle to the 3-D space...
198  if (radius > RHT_MIN_RADIUS && radius < RHT_MAX_RADIUS) {
199  accumulator.accumulate((int)(center.x / RHT_XY_SCALE),
200  (int)(center.y / RHT_XY_SCALE),
201  (int)(radius / RHT_RADIUS_SCALE));
202  }
203 
204  gettimeofday(&end, NULL);
205  }
206 
207  // Find the most dense region, and decide on the ball.
208  int max, x_max, y_max, r_max;
209  max = accumulator.getMax(x_max, y_max, r_max);
210 
211  cout << "Max vote is with " << max << " of the " << num_iter << " ones." << endl;
212 
213  // Only when votes exceeds a ratio will it be considered a real circle
214  if (max > num_iter * RHT_MIN_VOTE_RATE)
215  {
216  center_in_roi_t center;
217  center.x = (float)(x_max * RHT_XY_SCALE + RHT_XY_SCALE/2);
218  center.y = (float)(y_max * RHT_XY_SCALE + RHT_XY_SCALE/2);
219  float c_r = (float)(r_max * RHT_RADIUS_SCALE + RHT_RADIUS_SCALE/2);
220 
221  // With circle fitting
222  for(vector< upoint_t >::iterator pos = pixels.begin(); pos != pixels.end(); )
223  {
224  if (TBY_RADIUS_DIFF(pos->x, pos->y, center.x, center.y, c_r)
225  > RHT_FITTING_DIST_DIFF)
226  {
227  pixels.erase(pos);
228  }
229  else
230  {
231  pos++;
232  }
233  }
234 
235  Circle c;
236  c.fitCircle(pixels);
237  c.count = max;
238  m_Circles.push_back(c);
239 
240  /*
241  // Without circle fitting
242  m_Circles.push_back(Circle(center, c_r, max));
243  */
244 
245  return 0;
246  }
247 
248  // Return... here in this algorithm, we only find at most one most likely circle,
249  // (if none found, return-ed above)
250  // so the return value here is always 1
251  return 1;
252 }
253 
254 
255 int
256 RhtCircleModel::getShapeCount(void) const
257 {
258  return m_Circles.size();
259 }
260 
261 Circle *
262 RhtCircleModel::getShape(int id) const
263 {
264  if (id < 0 || (unsigned int)id >= m_Circles.size()) {
265  return NULL;
266  } else {
267  return const_cast<Circle*>(&m_Circles[id]); // or use const Shape* def?!...
268  }
269 }
270 
271 
272 Circle *
273 RhtCircleModel::getMostLikelyShape(void) const
274 {
275  if (m_Circles.size() == 0) {
276  return NULL;
277  } else if (m_Circles.size() == 1) {
278  return const_cast<Circle*>(&m_Circles[0]); // or use const Shape* def?!...
279  } else {
280  int cur=0;
281  for (unsigned int i=1; i < m_Circles.size(); ++i) {
282  if (m_Circles[i].count > m_Circles[cur].count) {
283  cur = i;
284  }
285  }
286  return const_cast<Circle*>(&m_Circles[cur]); // or use const Shape* definition?!...
287  }
288 }
289 
290 
291 void
292 RhtCircleModel::calcCircle(
293  const upoint_t& p1,
294  const upoint_t& p2,
295  const upoint_t& p3,
296  center_in_roi_t& center,
297  float& radius)
298  // Given three points p1, p2, p3,
299  // this function calculates the center and radius
300  // of the circle that is determined
301 {
302  const int &x1=p1.x, &y1=p1.y, &x2=p2.x, &y2=p2.y, &x3=p3.x, &y3=p3.y;
303  float dx, dy;
304  int div = 2*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1));
305 
306  if (div == 0)
307  {
308  // p1, p2 and p3 are in a straight line.
309  radius = -1.0;
310  return;
311  }
312  center.x = ((float)((x2*x2+y2*y2-x1*x1-y1*y1)*(y3-y1)
313  -(x3*x3+y3*y3-x1*x1-y1*y1)*(y2-y1))
314  /div);
315  center.y = ((float)((x2-x1)*(x3*x3+y3*y3-x1*x1-y1*y1)
316  -(x3-x1)*(x2*x2+y2*y2-x1*x1-y1*y1))
317  /div);
318  dx = center.x - x1;
319  dy = center.y - y1;
320  radius = (float)sqrt(dx*dx+dy*dy);
321 }
322 
323 } // end namespace firevision
Fawkes library namespace.
unsigned int y
y coordinate
Definition: types.h:36
STL namespace.
unsigned int x
x coordinate
Definition: types.h:35
unsigned int width
ROI width.
Definition: roi.h:121
Region of interest.
Definition: roi.h:58
float x
x in pixels
Definition: types.h:40
Circle shape.
Definition: circle.h:45
unsigned char * get_roi_buffer_start(unsigned char *buffer) const
Get ROI buffer start.
Definition: roi.cpp:556
Point with cartesian coordinates as unsigned integers.
Definition: types.h:34
void fitCircle(std::vector< fawkes::upoint_t > &points)
Fit circle.
Definition: circle.cpp:75
unsigned int height
ROI height.
Definition: roi.h:123
unsigned int line_step
line step
Definition: roi.h:129
int count
Number of pixels.
Definition: circle.h:64
Center in ROI.
Definition: types.h:39
float y
y in pixels
Definition: types.h:41