thinning.cpp
Go to the documentation of this file.
1 #include <tuw_voronoi_map/thinning.h>
2 #include <ros/ros.h>
3 
4 namespace voronoi_map
5 {
6 
7  void sceletonize_iteration(cv::Mat& img, int iter)
8  {
9  CV_Assert(img.channels() == 1);
10  CV_Assert(img.depth() != sizeof(uchar));
11  CV_Assert(img.rows > 3 && img.cols > 3);
12 
13  cv::Mat marker = cv::Mat::zeros(img.size(), CV_8UC1);
14 
15  int nRows = img.rows;
16  int nCols = img.cols;
17 
18  if(img.isContinuous())
19  {
20  nCols *= nRows;
21  nRows = 1;
22  }
23 
24  int x, y;
25  uchar *pAbove;
26  uchar *pCurr;
27  uchar *pBelow;
28  uchar *nw, *no, *ne; // north (pAbove)
29  uchar *we, *me, *ea;
30  uchar *sw, *so, *se; // south (pBelow)
31 
32  uchar *pDst;
33 
34  // initialize row pointers
35  pAbove = NULL;
36  pCurr = img.ptr<uchar>(0);
37  pBelow = img.ptr<uchar>(1);
38 
39  for(y = 1; y < img.rows - 1; ++y)
40  {
41  // shift the rows up by one
42  pAbove = pCurr;
43  pCurr = pBelow;
44  pBelow = img.ptr<uchar>(y + 1);
45 
46  pDst = marker.ptr<uchar>(y);
47 
48  // initialize col pointers
49  no = &(pAbove[0]);
50  ne = &(pAbove[1]);
51  me = &(pCurr[0]);
52  ea = &(pCurr[1]);
53  so = &(pBelow[0]);
54  se = &(pBelow[1]);
55 
56  for(x = 1; x < img.cols - 1; ++x)
57  {
58  // shift col pointers left by one (scan left to right)
59  nw = no;
60  no = ne;
61  ne = &(pAbove[x + 1]);
62  we = me;
63  me = ea;
64  ea = &(pCurr[x + 1]);
65  sw = so;
66  so = se;
67  se = &(pBelow[x + 1]);
68 
69  int A = (*no == 0 && *ne == 1) + (*ne == 0 && *ea == 1) +
70  (*ea == 0 && *se == 1) + (*se == 0 && *so == 1) +
71  (*so == 0 && *sw == 1) + (*sw == 0 && *we == 1) +
72  (*we == 0 && *nw == 1) + (*nw == 0 && *no == 1);
73  int B = *no + *ne + *ea + *se + *so + *sw + *we + *nw;
74  int m1 = iter == 0 ? (*no * *ea * *so) : (*no * *ea * *we);
75  int m2 = iter == 0 ? (*ea * *so * *we) : (*no * *so * *we);
76 
77  if(A == 1 && (B >= 2 && B <= 6) && m1 == 0 && m2 == 0)
78  pDst[x] = 1;
79  }
80  }
81 
82  img &= ~marker;
83  }
84 
85  void sceletonize(const cv::Mat& src, cv::Mat& dst)
86  {
87  dst = src.clone();
88  dst /= 255; // convert to binary image
89 
90  cv::Mat prev = cv::Mat::zeros(dst.size(), CV_8UC1);
91  cv::Mat diff;
92 
93  do
94  {
95  sceletonize_iteration(dst, 0);
96  sceletonize_iteration(dst, 1);
97  cv::absdiff(dst, prev, diff);
98  dst.copyTo(prev);
99  }
100  while(cv::countNonZero(diff) > 0);
101 
102  dst *= 255;
103  }
104 
105 
106  std::queue<Index> q;
107 
108 
109  void greyscale_thinning(const cv::Mat& src, cv::Mat& dst)
110  {
111  int initCosts = 2;
112  int pathCosts = 1;
113 
114  for(int i = 1; i < dst.rows - 1; i++)
115  {
116  for(int j = 1; j < dst.cols - 1; j++)
117  {
118  float pAct = src.at<float_t>(i, j);
119 
120  //skip unknown area
121  if(pAct > 0 && i > 0 && i < src.rows - 1 && j > 0 && j < src.cols - 1)
122  {
123  //Find neighbour size
124  float p[8] =
125  {
126  src.at<float_t>(i, j - 1),
127  src.at<float_t>(i - 1, j - 1),
128  src.at<float_t>(i - 1, j),
129  src.at<float_t>(i - 1, j + 1),
130  src.at<float_t>(i, j + 1),
131  src.at<float_t>(i + 1, j + 1),
132  src.at<float_t>(i + 1, j),
133  src.at<float_t>(i + 1, j - 1)
134  };
135 
136  bool bigger[8] = {p[0] > pAct, p[1] > pAct, p[2] > pAct, p[3] > pAct, p[4] > pAct, p[5] > pAct, p[6] > pAct, p[7] > pAct};
137  bool smaller[8] = {p[0] < pAct, p[1] < pAct, p[2] < pAct, p[3] < pAct, p[4] < pAct, p[5] < pAct, p[6] < pAct, p[7] < pAct};
138 
139 
140 
141 
142  //Find line Saddle points
143  if(bigger[4] && bigger[0] && !bigger[2] && !bigger[6])
144  {
145  //straight saddlepoint
146  // -
147  // + o +
148  // -
149  q.push(Index(i, j, pAct));
150  dst.at<int8_t>(i, j) = initCosts;
151  }
152  else if(!bigger[4] && !bigger[0] && bigger[2] && bigger[6])
153  {
154  //straight saddlepoint
155  // +
156  // - o -
157  // +
158  q.push(Index(i, j, pAct));
159  dst.at<int8_t>(i, j) = initCosts;
160  }
161  else if(bigger[1] && bigger[5] && !(bigger[2] && bigger[3] && bigger[4]) && !(bigger[0] && bigger[7] && bigger[6]))
162  {
163  //diagonal saddlepoint
164  // + - -
165  // - o -
166  // - - +
167  q.push(Index(i, j, pAct));
168  dst.at<int8_t>(i, j) = initCosts;
169  }
170  else if(bigger[3] && bigger[7] && !(bigger[2] && bigger[1] && bigger[0]) && !(bigger[4] && bigger[5] && bigger[6]))
171  {
172  //diagonal saddlepoint
173  // - - +
174  // - o -
175  // + - -
176  q.push(Index(i, j, pAct));
177  dst.at<int8_t>(i, j) = initCosts;
178  }
179 
180  //Find additional saddle points
181 
182  // - - + + - -
183  // - o - - o -
184  // + + ....
185 
186  else if(!bigger[0] && !bigger[1] && !bigger[2] && bigger[3] && !bigger[4] && bigger[6])
187  {
188  q.push(Index(i, j, pAct));
189  dst.at<int8_t>(i, j) = initCosts;
190  }
191  else if(!bigger[0] && bigger[1] && !bigger[2] && !bigger[3] && !bigger[4] && bigger[6])
192  {
193  q.push(Index(i, j, pAct));
194  dst.at<int8_t>(i, j) = initCosts;
195  }
196 
197  else if(!bigger[0] && bigger[2] && !bigger[4] && bigger[5] && !bigger[6] && !bigger[7])
198  {
199  q.push(Index(i, j, pAct));
200  dst.at<int8_t>(i, j) = initCosts;
201  }
202  else if(!bigger[0] && bigger[2] && !bigger[4] && !bigger[5] && !bigger[6] && bigger[7])
203  {
204  q.push(Index(i, j, pAct));
205  dst.at<int8_t>(i, j) = initCosts;
206  }
207 
208  else if(!bigger[0] && bigger[1] && !bigger[2] && bigger[4] && !bigger[6] && !bigger[7])
209  {
210  q.push(Index(i, j, pAct));
211  dst.at<int8_t>(i, j) = initCosts;
212  }
213  else if(!bigger[0] && !bigger[1] && !bigger[2] && bigger[4] && !bigger[6] && bigger[7])
214  {
215  q.push(Index(i, j, pAct));
216  dst.at<int8_t>(i, j) = initCosts;
217  }
218 
219  else if(bigger[0] && !bigger[2] && bigger[3] && !bigger[4] && !bigger[5] && !bigger[6])
220  {
221  q.push(Index(i, j, pAct));
222  dst.at<int8_t>(i, j) = initCosts;
223  }
224  else if(bigger[0] && !bigger[2] && !bigger[3] && !bigger[4] && bigger[5] && !bigger[6])
225  {
226  q.push(Index(i, j, pAct));
227  dst.at<int8_t>(i, j) = initCosts;
228  }
229  //Maxima
230  else if(!bigger[0] && !bigger[1] && !bigger[2] && !bigger[3] && !bigger[4] && !bigger[5] && !bigger[6] && !bigger[7])
231  {
232  dst.at<int8_t>(i, j) = initCosts;
233  q.push(Index(i, j, pAct));
234 
235 
236  }
237 
238  //Find Edges
239 
240 
241  if(smaller[0] && smaller[1] && bigger[3] && smaller[5] && smaller[6] && smaller[7])
242  {
243  //left top edge
244  // - +
245  // - o
246  // - - -
247 
248  q.push(Index(i, j, pAct));
249  dst.at<int8_t>(i, j) = initCosts;
250  }
251  else if(smaller[1] && smaller[2] && smaller[3] && smaller[4] && smaller[5] && bigger[7])
252  {
253  //left bottom edge
254  // - - -
255  // o -
256  // + -
257 
258  q.push(Index(i, j, pAct));
259  dst.at<int8_t>(i, j) = initCosts;
260  }
261  else if(bigger[1] && smaller[3] && smaller[4] && smaller[5] && smaller[6] && smaller[7])
262  {
263  //right top edge
264  // + -
265  // o -
266  // - - -
267 
268  q.push(Index(i, j, pAct));
269  dst.at<int8_t>(i, j) = initCosts;
270  }
271  else if(smaller[0] && smaller[1] && smaller[2] && smaller[3] && bigger[5] && smaller[7])
272  {
273  //right bottom edge
274  // - - -
275  // - o
276  // - +
277 
278  q.push(Index(i, j, pAct));
279  dst.at<int8_t>(i, j) = initCosts;
280  }
281 
282 
283 
284 
285  if(pAct > 1 && i > 1 && j > 1 && i < src.rows - 2 && j < src.cols - 2) //Dont look for duplex saddle points at the edges of the map
286  {
287 
288  //Find duplex saddle points
289  if((!smaller[5] || !smaller[6] || !smaller[4]) && !bigger[7] && !bigger[3])
290  {
291  //Diagonal duplex saddle Point
292  // + + -
293  // + o -
294  // - o +
295  // - + +
296 
297  int m = i - 1;
298  int n = j - 1;
299 
300  float pot = src.at<float_t>(m, n);
301 
302  float pp2[8] =
303  {
304  src.at<float_t>(m, n - 1),
305  src.at<float_t>(m - 1, n - 1),
306  src.at<float_t>(m - 1, n),
307  src.at<float_t>(m - 1, n + 1),
308  src.at<float_t>(m, n + 1),
309  src.at<float_t>(m + 1, n + 1),
310  src.at<float_t>(m + 1, n),
311  src.at<float_t>(m + 1, n - 1)
312  };
313  bool bigger2[8] = {pp2[0] > pot, pp2[1] > pot, pp2[2] > pot, pp2[3] > pot, pp2[4] > pot, pp2[5] > pot, pp2[6] > pot, pp2[7] > pot};
314  bool smaller2[8] = {pp2[0] < pot, pp2[1] < pot, pp2[2] < pot, pp2[3] < pot, pp2[4] < pot, pp2[5] < pot, pp2[6] < pot, pp2[7] < pot};
315 
316  if((!smaller2[0] || !smaller2[1] || !smaller2[2]) && !bigger2[3] && !bigger2[7] && pot > 1)
317  {
318  q.push(Index(i, j, pAct));
319  dst.at<int8_t>(i, j) = initCosts;
320 
321  q.push(Index(m, n, pAct));
322  dst.at<int8_t>(m, n) = initCosts;
323  }
324  }
325 
326  if((!smaller[2] || !smaller[3] || !smaller[4]) && !bigger[1] && !bigger[5])
327  {
328  //Diagonal duplex saddle Point
329  // - + +
330  // - o +
331  // + o -
332  // + + -
333 
334  int m = i + 1;
335  int n = j - 1;
336  float pot = src.at<float_t>(m, n);
337 
338  float pp2[8] =
339  {
340  src.at<float_t>(m, n - 1),
341  src.at<float_t>(m - 1, n - 1),
342  src.at<float_t>(m - 1, n),
343  src.at<float_t>(m - 1, n + 1),
344  src.at<float_t>(m, n + 1),
345  src.at<float_t>(m + 1, n + 1),
346  src.at<float_t>(m + 1, n),
347  src.at<float_t>(m + 1, n - 1)
348  };
349  bool bigger2[8] = {pp2[0] > pot, pp2[1] > pot, pp2[2] > pot, pp2[3] > pot, pp2[4] > pot, pp2[5] > pot, pp2[6] > pot, pp2[7] > pot};
350  bool smaller2[8] = {pp2[0] < pot, pp2[1] < pot, pp2[2] < pot, pp2[3] < pot, pp2[4] < pot, pp2[5] < pot, pp2[6] < pot, pp2[7] < pot};
351 
352  if((!smaller2[0] || !smaller2[7] || !smaller2[6]) && !bigger2[1] && !bigger2[5] && pot > 1)
353  {
354  q.push(Index(i, j, pAct));
355  dst.at<int8_t>(i, j) = initCosts;
356 
357  q.push(Index(m, n, pAct));
358  dst.at<int8_t>(m, n) = initCosts;
359  }
360  }
361 
362  if((!smaller[2] || !smaller[3] || !smaller[1]) && !bigger[0] && !bigger[4])
363  {
364  //Straight duplex saddle Point
365 
366  // + + +
367  // - o -
368  // - o -
369  // + + +
370 
371  int m = i + 1;
372  int n = j;
373 
374  float pot = src.at<float_t>(m, n);
375 
376 
377  float pp2[8] =
378  {
379  src.at<float_t>(m, n - 1),
380  src.at<float_t>(m - 1, n - 1),
381  src.at<float_t>(m - 1, n),
382  src.at<float_t>(m - 1, n + 1),
383  src.at<float_t>(m, n + 1),
384  src.at<float_t>(m + 1, n + 1),
385  src.at<float_t>(m + 1, n),
386  src.at<float_t>(m + 1, n - 1)
387  };
388 
389 
390  bool bigger2[8] = {pp2[0] > pot, pp2[1] > pot, pp2[2] > pot, pp2[3] > pot, pp2[4] > pot, pp2[5] > pot, pp2[6] > pot, pp2[7] > pot};
391  bool smaller2[8] = {pp2[0] < pot, pp2[1] < pot, pp2[2] < pot, pp2[3] < pot, pp2[4] < pot, pp2[5] < pot, pp2[6] < pot, pp2[7] < pot};
392 
393  if((!smaller2[5] || !smaller2[6] || !smaller2[7]) && !bigger2[0] && !bigger2[4] && pot > 1)
394  {
395  q.push(Index(i, j, pAct));
396  dst.at<int8_t>(i, j) = initCosts;
397 
398  q.push(Index(m, n, pAct));
399  dst.at<int8_t>(m, n) = initCosts;
400  }
401  }
402 
403  if((!smaller[7] || !smaller[0] || !smaller[1]) && !bigger[2] && !bigger[6])
404  {
405  //Straight duplex saddle Point
406 
407  // + - - +
408  // + o o +
409  // + - - +
410 
411  int m = i;
412  int n = j + 1;
413  float pot = src.at<float_t>(m, n);
414 
415  float pp2[8] =
416  {
417  src.at<float_t>(m, n - 1),
418  src.at<float_t>(m - 1, n - 1),
419  src.at<float_t>(m - 1, n),
420  src.at<float_t>(m - 1, n + 1),
421  src.at<float_t>(m, n + 1),
422  src.at<float_t>(m + 1, n + 1),
423  src.at<float_t>(m + 1, n),
424  src.at<float_t>(m + 1, n - 1)
425  };
426  bool bigger2[8] = {pp2[0] > pot, pp2[1] > pot, pp2[2] > pot, pp2[3] > pot, pp2[4] > pot, pp2[5] > pot, pp2[6] > pot, pp2[7] > pot};
427  bool smaller2[8] = {pp2[0] < pot, pp2[1] < pot, pp2[2] < pot, pp2[3] < pot, pp2[4] < pot, pp2[5] < pot, pp2[6] < pot, pp2[7] < pot};
428 
429  if((!smaller2[3] || !smaller2[4] || !smaller2[5]) && !bigger2[2] && !bigger2[6] && pot > 1)
430  {
431  q.push(Index(i, j, pAct));
432  dst.at<int8_t>(i, j) = initCosts;
433 
434  q.push(Index(m, n, pAct));
435  dst.at<int8_t>(m, n) = initCosts;
436  }
437  }
438  }
439 
440  }
441  }
442  }
443 
444  while(!q.empty())
445  {
446  Index pixel = q.front();
447  q.pop();
448 
449  int costs = dst.at<uint8_t>(pixel.i, pixel.j);
450  Index next = getMaximumNeighbour(pixel.i, pixel.j, src, dst);
451 
452  if(next.i >= 0 && next.j >= 0 && dst.at<uint8_t>(next.i, next.j) == 0)
453  {
454  if(src.at<float_t>(next.i, next.j) > 0)
455  {
456  dst.at<uint8_t>(next.i, next.j) = costs + pathCosts;
457  q.push(next);
458  }
459  }
460  }
461  }
462 
463  Index getMaximumNeighbour(int i, int j, const cv::Mat& src, cv::Mat& dst)
464  {
465  //Return on boarder
466  if(i == 0 || i == src.rows - 1 || j == 0 || j == src.cols - 1)
467  return Index(-1, -1, -1);
468 
469  int costs = dst.at<uint8_t>(i, j);
470  Index p[8] =
471  {
472  Index(i, j - 1, src.at<float_t>(i, j - 1)),
473  Index(i - 1, j - 1, src.at<float_t>(i - 1, j - 1)),
474  Index(i - 1, j, src.at<float_t>(i - 1, j)),
475  Index(i - 1, j + 1, src.at<float_t>(i - 1, j + 1)),
476  Index(i, j + 1, src.at<float_t>(i, j + 1)),
477  Index(i + 1, j + 1, src.at<float_t>(i + 1, j + 1)),
478  Index(i + 1, j, src.at<float_t>(i + 1, j)),
479  Index(i + 1, j - 1, src.at<float_t>(i + 1, j - 1))
480  };
481 
482  Index act = Index(-1, -1, -1);
483 
484  for(int i = 0; i < 8; i++)
485  {
486  if(p[i].potential > act.potential)
487  {
488  int actCosts = dst.at<uint8_t>(p[i].i, p[i].j);
489 
490  if(actCosts != costs)
491  {
492  act = p[i];
493  }
494 
495  }
496  }
497 
498 
499  return act;
500  }
501 
502 
503 
504 }
505 
voronoi_map::Index
Definition: thinning.h:31
voronoi_map::greyscale_thinning
void greyscale_thinning(const cv::Mat &src, cv::Mat &dst)
Function for finding the ridge of a distance graph.
Definition: thinning.cpp:109
voronoi_map::getMaximumNeighbour
Index getMaximumNeighbour(int i, int j, const cv::Mat &src, cv::Mat &dst)
function for finding the maximum Neighbour ignoring the last detected pixels
Definition: thinning.cpp:463
voronoi_map::Index::i
int i
Definition: thinning.h:44
voronoi_map::sceletonize_iteration
void sceletonize_iteration(cv::Mat &img, int iter)
Perform one thinning iteration.(Normally you wouldn&#39;t call this function directly from your code) ...
Definition: thinning.cpp:7
y
TFSIMD_FORCE_INLINE const tfScalar & y() const
voronoi_map::Index::j
int j
Definition: thinning.h:45
voronoi_map::Index::potential
float potential
Definition: thinning.h:46
voronoi_map::sceletonize
void sceletonize(const cv::Mat &src, cv::Mat &dst)
Function for thinning the given binary image (Paper Zhang-Suen Thinning)
Definition: thinning.cpp:85
voronoi_map
Definition: thinning.h:7
x
TFSIMD_FORCE_INLINE const tfScalar & x() const
voronoi_map::q
std::queue< Index > q
Definition: thinning.cpp:106

tuw_voronoi_graph
Author(s): Benjamin Binder
autogenerated on Mon Jun 10 2019 15:42:44