1 // This file is part of OpenCV project.
2 // It is subject to the license terms in the LICENSE file found in the top-level
3 // directory of this distribution and at http://opencv.org/license.html.
4
5 #include "precomp.hpp"
6 #include "msd_pyramid.hpp"
7
8 namespace cv
9 {
10 namespace xfeatures2d
11 {
12
13 class TBMR_Impl CV_FINAL : public TBMR
14 {
15 public:
16 struct Params
17 {
Paramscv::xfeatures2d::CV_FINAL::Params18 Params(int _min_area = 60, float _max_area_relative = 0.01,
19 float _scale = 1.5, int _n_scale = -1)
20 {
21 CV_Assert(_min_area >= 0);
22 CV_Assert(_max_area_relative >=
23 std::numeric_limits<float>::epsilon());
24
25 minArea = _min_area;
26 maxAreaRelative = _max_area_relative;
27 scale = _scale;
28 n_scale = _n_scale;
29 }
30
31 uint minArea;
32 float maxAreaRelative;
33 int n_scale;
34 float scale;
35 };
36
TBMR_Impl(const Params & _params)37 explicit TBMR_Impl(const Params &_params) : params(_params) {}
38
~TBMR_Impl()39 virtual ~TBMR_Impl() CV_OVERRIDE {}
40
setMinArea(int minArea)41 virtual void setMinArea(int minArea) CV_OVERRIDE
42 {
43 params.minArea = std::max(minArea, 0);
44 }
getMinArea() const45 int getMinArea() const CV_OVERRIDE { return params.minArea; }
46
setMaxAreaRelative(float maxAreaRelative)47 virtual void setMaxAreaRelative(float maxAreaRelative) CV_OVERRIDE
48 {
49 params.maxAreaRelative =
50 std::max(maxAreaRelative, std::numeric_limits<float>::epsilon());
51 }
getMaxAreaRelative() const52 virtual float getMaxAreaRelative() const CV_OVERRIDE
53 {
54 return params.maxAreaRelative;
55 }
setScaleFactor(float scale_factor)56 virtual void setScaleFactor(float scale_factor) CV_OVERRIDE
57 {
58 params.scale = std::max(scale_factor, 1.f);
59 }
getScaleFactor() const60 virtual float getScaleFactor() const CV_OVERRIDE { return params.scale; }
setNScales(int n_scales)61 virtual void setNScales(int n_scales) CV_OVERRIDE
62 {
63 params.n_scale = n_scales;
64 }
getNScales() const65 virtual int getNScales() const CV_OVERRIDE { return params.n_scale; }
66
67 virtual void detect(InputArray image,
68 CV_OUT std::vector<KeyPoint> &keypoints,
69 InputArray mask = noArray()) CV_OVERRIDE;
70
71 virtual void detect(InputArray image,
72 CV_OUT std::vector<Elliptic_KeyPoint> &keypoints,
73 InputArray mask = noArray()) CV_OVERRIDE;
74
75 virtual void
76 detectAndCompute(InputArray image, InputArray mask,
77 CV_OUT std::vector<Elliptic_KeyPoint> &keypoints,
78 OutputArray descriptors,
79 bool useProvidedKeypoints = false) CV_OVERRIDE;
80
zfindroot(uint * parent,uint p)81 CV_INLINE uint zfindroot(uint *parent, uint p)
82 {
83 if (parent[p] == p)
84 return p;
85 else
86 return parent[p] = zfindroot(parent, parent[p]);
87 }
88
89 // Calculate the Component tree. Based on the order of S, it will be a
90 // min or max tree.
calcMinMaxTree(Mat ima)91 void calcMinMaxTree(Mat ima)
92 {
93 int rs = ima.rows;
94 int cs = ima.cols;
95 uint imSize = (uint)rs * cs;
96
97 std::array<int, 4> offsets = {
98 -ima.cols, -1, 1, ima.cols
99 }; // {-1,0}, {0,-1}, {0,1}, {1,0} yx
100 std::array<Vec2i, 4> offsetsv = { Vec2i(0, -1), Vec2i(-1, 0),
101 Vec2i(1, 0), Vec2i(0, 1) }; // xy
102 AutoBuffer<uint> zparb(imSize);
103 AutoBuffer<uint> rootb(imSize);
104 AutoBuffer<uint> rankb(imSize);
105 memset(rankb.data(), 0, imSize * sizeof(uint));
106 uint* zpar = zparb.data();
107 uint *root = rootb.data();
108 uint *rank = rankb.data();
109 parent = Mat(rs, cs, CV_32S); // unsigned
110 AutoBuffer<bool> dejaVub(imSize);
111 memset(dejaVub.data(), 0, imSize * sizeof(bool));
112 bool* dejaVu = dejaVub.data();
113
114 const uint *S_ptr = S.ptr<const uint>();
115 uint *parent_ptr = parent.ptr<uint>();
116 Vec<uint, 6> *imaAttribute = imaAttributes.ptr<Vec<uint, 6>>();
117
118 for (int i = imSize - 1; i >= 0; --i)
119 {
120 uint p = S_ptr[i];
121
122 Vec2i idx_p(p % cs, p / cs);
123 // make set
124 {
125 parent_ptr[p] = p;
126 zpar[p] = p;
127 root[p] = p;
128 dejaVu[p] = true;
129 imaAttribute[p][0] = 1; // area
130 imaAttribute[p][1] = idx_p[0]; // sum_x
131 imaAttribute[p][2] = idx_p[1]; // sum_y
132 imaAttribute[p][3] = idx_p[0] * idx_p[1]; // sum_xy
133 imaAttribute[p][4] = idx_p[0] * idx_p[0]; // sum_xx
134 imaAttribute[p][5] = idx_p[1] * idx_p[1]; // sum_yy
135 }
136
137 uint x = p; // zpar of p
138 for (unsigned k = 0; k < offsets.size(); ++k)
139 {
140 uint q = p + offsets[k];
141
142 Vec2i q_idx = idx_p + offsetsv[k];
143 bool inBorder = q_idx[0] >= 0 && q_idx[0] < ima.cols &&
144 q_idx[1] >= 0 &&
145 q_idx[1] < ima.rows; // filter out border cases
146
147 if (inBorder && dejaVu[q]) // remove first check
148 // obsolete
149 {
150 uint r = zfindroot(zpar, q);
151 if (r != x) // make union
152 {
153 parent_ptr[root[r]] = p;
154 // accumulate information
155 imaAttribute[p][0] += imaAttribute[root[r]][0]; // area
156 imaAttribute[p][1] += imaAttribute[root[r]][1]; // sum_x
157 imaAttribute[p][2] += imaAttribute[root[r]][2]; // sum_y
158 imaAttribute[p][3] +=
159 imaAttribute[root[r]][3]; // sum_xy
160 imaAttribute[p][4] +=
161 imaAttribute[root[r]][4]; // sum_xx
162 imaAttribute[p][5] +=
163 imaAttribute[root[r]][5]; // sum_yy
164
165 if (rank[x] < rank[r])
166 {
167 // we merge p to r
168 zpar[x] = r;
169 root[r] = p;
170 x = r;
171 }
172 else if (rank[r] < rank[p])
173 {
174 // merge r to p
175 zpar[r] = p;
176 }
177 else
178 {
179 // same height
180 zpar[r] = p;
181 rank[p] += 1;
182 }
183 }
184 }
185 }
186 }
187 }
188
calculateTBMRs(const Mat & image,std::vector<Elliptic_KeyPoint> & tbmrs,const Mat & mask,float scale,int octave)189 void calculateTBMRs(const Mat &image, std::vector<Elliptic_KeyPoint> &tbmrs,
190 const Mat &mask, float scale, int octave)
191 {
192 uint imSize = image.cols * image.rows;
193 uint maxArea =
194 static_cast<uint>(params.maxAreaRelative * imSize * scale);
195 uint minArea = static_cast<uint>(params.minArea * scale);
196
197 if (parent.empty() || parent.size != image.size)
198 parent = Mat(image.rows, image.cols, CV_32S);
199
200 if (imaAttributes.empty() || imaAttributes.size != image.size)
201 imaAttributes = Mat(image.rows, image.cols, CV_32SC(6));
202
203 calcMinMaxTree(image);
204
205 const Vec<uint, 6> *imaAttribute =
206 imaAttributes.ptr<const Vec<uint, 6>>();
207 const uint8_t *ima_ptr = image.ptr<const uint8_t>();
208 const uint *S_ptr = S.ptr<const uint>();
209 uint *parent_ptr = parent.ptr<uint>();
210
211 // canonization
212 for (uint i = 0; i < imSize; ++i)
213 {
214 uint p = S_ptr[i];
215 uint q = parent_ptr[p];
216 if (ima_ptr[parent_ptr[q]] == ima_ptr[q])
217 parent_ptr[p] = parent_ptr[q];
218 }
219
220 // TBMRs extraction
221 //------------------------------------------------------------------------
222 // small variant of the given algorithm in the paper. For each
223 // critical node having more than one child, we check if the
224 // largest region containing this node without any change of
225 // topology is above its parent, if not, discard this critical
226 // node.
227 //
228 // note also that we do not select the critical nodes themselves
229 // as final TBMRs
230 //--------------------------------------------------------------------------
231
232 AutoBuffer<uint> numSonsb(imSize);
233 memset(numSonsb.data(), 0, imSize * sizeof(uint));
234 uint* numSons = numSonsb.data();
235 uint vecNodesSize = imaAttribute[S_ptr[0]][0]; // area
236 AutoBuffer<uint> vecNodesb(vecNodesSize);
237 memset(vecNodesb.data(), 0, vecNodesSize * sizeof(uint));
238 uint *vecNodes = vecNodesb.data(); // area
239 uint numNodes = 0;
240
241 // leaf to root propagation to select the canonized nodes
242 for (int i = imSize - 1; i >= 0; --i)
243 {
244 uint p = S_ptr[i];
245 if (parent_ptr[p] == p || ima_ptr[p] != ima_ptr[parent_ptr[p]])
246 {
247 vecNodes[numNodes++] = p;
248 if (imaAttribute[p][0] >= minArea) // area
249 numSons[parent_ptr[p]]++;
250 }
251 }
252
253 AutoBuffer<bool> isSeenb(imSize);
254 memset(isSeenb.data(), 0, imSize * sizeof(bool));
255 bool *isSeen = isSeenb.data();
256
257 // parent of critical leaf node
258 AutoBuffer<bool> isParentofLeafb(imSize);
259 memset(isParentofLeafb.data(), 0, imSize * sizeof(bool));
260 bool* isParentofLeaf = isParentofLeafb.data();
261
262 for (uint i = 0; i < vecNodesSize; i++)
263 {
264 uint p = vecNodes[i];
265 if (numSons[p] == 0 && numSons[parent_ptr[p]] == 1)
266 isParentofLeaf[parent_ptr[p]] = true;
267 }
268
269 uint numTbmrs = 0;
270 AutoBuffer<uint> vecTbmrsb(numNodes);
271 uint* vecTbmrs = vecTbmrsb.data();
272 for (uint i = 0; i < vecNodesSize; i++)
273 {
274 uint p = vecNodes[i];
275 if (numSons[p] == 1 && !isSeen[p] && imaAttribute[p][0] <= maxArea)
276 {
277 uint num_ancestors = 0;
278 uint pt = p;
279 uint po = pt;
280 while (numSons[pt] == 1 && imaAttribute[pt][0] <= maxArea)
281 {
282 isSeen[pt] = true;
283 num_ancestors++;
284 po = pt;
285 pt = parent_ptr[pt];
286 }
287 if (!isParentofLeaf[p] || num_ancestors > 1)
288 {
289 vecTbmrs[numTbmrs++] = po;
290 }
291 }
292 }
293 // end of TBMRs extraction
294 //------------------------------------------------------------------------
295
296 // compute best fitting ellipses
297 //------------------------------------------------------------------------
298 for (uint i = 0; i < numTbmrs; i++)
299 {
300 uint p = vecTbmrs[i];
301 double area = static_cast<double>(imaAttribute[p][0]);
302 double sum_x = static_cast<double>(imaAttribute[p][1]);
303 double sum_y = static_cast<double>(imaAttribute[p][2]);
304 double sum_xy = static_cast<double>(imaAttribute[p][3]);
305 double sum_xx = static_cast<double>(imaAttribute[p][4]);
306 double sum_yy = static_cast<double>(imaAttribute[p][5]);
307
308 // Barycenter:
309 double x = sum_x / area;
310 double y = sum_y / area;
311
312 double i20 = sum_xx - area * x * x;
313 double i02 = sum_yy - area * y * y;
314 double i11 = sum_xy - area * x * y;
315 double n = i20 * i02 - i11 * i11;
316 if (n != 0)
317 {
318 double a = (i02 / n) * (area - 1) / 4;
319 double b = (-i11 / n) * (area - 1) / 4;
320 double c = (i20 / n) * (area - 1) / 4;
321
322 // filter out some non meaningful ellipses
323 double a1 = a;
324 double b1 = b;
325 double c1 = c;
326 uint ai = 0;
327 uint bi = 0;
328 uint ci = 0;
329 if (a > 0)
330 {
331 if (a < 0.00005)
332 a1 = 0;
333 else if (a < 0.0001)
334 {
335 a1 = 0.0001;
336 }
337 else
338 {
339 ai = (uint)(10000 * a);
340 a1 = (double)ai / 10000;
341 }
342 }
343 else
344 {
345 if (a > -0.00005)
346 a1 = 0;
347 else if (a > -0.0001)
348 a1 = -0.0001;
349 else
350 {
351 ai = (uint)(10000 * (-a));
352 a1 = -(double)ai / 10000;
353 }
354 }
355
356 if (b > 0)
357 {
358 if (b < 0.00005)
359 b1 = 0;
360 else if (b < 0.0001)
361 {
362 b1 = 0.0001;
363 }
364 else
365 {
366 bi = (uint)(10000 * b);
367 b1 = (double)bi / 10000;
368 }
369 }
370 else
371 {
372 if (b > -0.00005)
373 b1 = 0;
374 else if (b > -0.0001)
375 b1 = -0.0001;
376 else
377 {
378 bi = (uint)(10000 * (-b));
379 b1 = -(double)bi / 10000;
380 }
381 }
382
383 if (c > 0)
384 {
385 if (c < 0.00005)
386 c1 = 0;
387 else if (c < 0.0001)
388 {
389 c1 = 0.0001;
390 }
391 else
392 {
393 ci = (uint)(10000 * c);
394 c1 = (double)ci / 10000;
395 }
396 }
397 else
398 {
399 if (c > -0.00005)
400 c1 = 0;
401 else if (c > -0.0001)
402 c1 = -0.0001;
403 else
404 {
405 ci = (uint)(10000 * (-c));
406 c1 = -(double)ci / 10000;
407 }
408 }
409 double v =
410 (a1 + c1 -
411 std::sqrt(a1 * a1 + c1 * c1 + 4 * b1 * b1 - 2 * a1 * c1)) /
412 2;
413
414 double l1 = 1. / std::sqrt((a + c +
415 std::sqrt(a * a + c * c +
416 4 * b * b - 2 * a * c)) /
417 2);
418 double l2 = 1. / std::sqrt((a + c -
419 std::sqrt(a * a + c * c +
420 4 * b * b - 2 * a * c)) /
421 2);
422 double minAxL = std::min(l1, l2);
423 double majAxL = std::max(l1, l2);
424
425 if (minAxL >= 1.5 && v != 0 &&
426 (mask.empty() ||
427 mask.at<uchar>(cvRound(y), cvRound(x)) != 0))
428 {
429 double theta = 0;
430 if (b == 0)
431 if (a < c)
432 theta = 0;
433 else
434 theta = CV_PI / 2.;
435 else
436 theta = CV_PI / 2. + 0.5 * std::atan2(2 * b, (a - c));
437
438 float size = (float)majAxL;
439
440 // not sure if we should scale or not scale x,y,axes,size
441 // (as scale is stored in si)
442 Elliptic_KeyPoint ekp(
443 Point2f((float)x, (float)y) * scale, (float)theta,
444 cv::Size2f((float)majAxL, (float)minAxL) * scale,
445 size * scale, scale);
446 ekp.octave = octave;
447 tbmrs.push_back(ekp);
448 }
449 }
450 }
451 //---------------------------------------------
452 }
453
454 Mat tempsrc;
455
456 // component tree representation (parent,S): see
457 // https://ieeexplore.ieee.org/document/6850018
458 Mat parent;
459 Mat S;
460 // moments: compound type of: (area, x, y, xy, xx, yy)
461 Mat imaAttributes;
462
463 Params params;
464 };
465
detect(InputArray _image,std::vector<KeyPoint> & keypoints,InputArray _mask)466 void TBMR_Impl::detect(InputArray _image, std::vector<KeyPoint> &keypoints,
467 InputArray _mask)
468 {
469 std::vector<Elliptic_KeyPoint> kp;
470 detect(_image, kp, _mask);
471 keypoints.resize(kp.size());
472 for (size_t i = 0; i < kp.size(); ++i)
473 keypoints[i] = kp[i];
474 }
475
detect(InputArray _image,std::vector<Elliptic_KeyPoint> & keypoints,InputArray _mask)476 void TBMR_Impl::detect(InputArray _image,
477 std::vector<Elliptic_KeyPoint> &keypoints,
478 InputArray _mask)
479 {
480 Mat mask = _mask.getMat();
481 Mat src = _image.getMat();
482
483 keypoints.clear();
484
485 if (src.empty())
486 return;
487
488 if (!mask.empty())
489 {
490 CV_Assert(mask.type() == CV_8UC1);
491 CV_Assert(mask.size == src.size);
492 }
493
494 if (!src.isContinuous())
495 {
496 src.copyTo(tempsrc);
497 src = tempsrc;
498 }
499
500 CV_Assert(src.depth() == CV_8U);
501
502 if (src.channels() != 1)
503 cv::cvtColor(src, src, cv::COLOR_BGR2GRAY);
504
505 int m_cur_n_scales =
506 params.n_scale > 0
507 ? params.n_scale
508 : 1 /*todo calculate optimal scale factor from image size*/;
509 float m_scale_factor = params.scale;
510
511 // track and eliminate duplicates introduced with multi scale position ->
512 // (size)
513 Mat dupl(src.rows / 4, src.cols / 4, CV_32F, cv::Scalar::all(0));
514 float *dupl_ptr = dupl.ptr<float>();
515
516 std::vector<Mat> pyr;
517 MSDImagePyramid scaleSpacer(src, m_cur_n_scales, m_scale_factor);
518 pyr = scaleSpacer.getImPyr();
519
520 int oct = 0;
521 for (auto &s : pyr)
522 {
523 float scale = ((float)s.cols) / pyr.begin()->cols;
524 std::vector<Elliptic_KeyPoint> kpts;
525
526 // append max tree tbmrs
527 sortIdx(s.reshape(1, 1), S,
528 SortFlags::SORT_ASCENDING | SortFlags::SORT_EVERY_ROW);
529 calculateTBMRs(s, kpts, mask, scale, oct);
530
531 // reverse instead of sort
532 flip(S, S, -1);
533 calculateTBMRs(s, kpts, mask, scale, oct);
534
535 if (oct == 0)
536 {
537 for (const auto &k : kpts)
538 {
539 dupl_ptr[(int)(k.pt.x / 4) +
540 (int)(k.pt.y / 4) * (src.cols / 4)] = k.size;
541 }
542 keypoints.insert(keypoints.end(), kpts.begin(), kpts.end());
543 }
544 else
545 {
546 for (const auto &k : kpts)
547 {
548 float &sz = dupl_ptr[(int)(k.pt.x / 4) +
549 (int)(k.pt.y / 4) * (src.cols / 4)];
550 // we hereby add only features that are at least 4 pixels away
551 // or have a significantly different size
552 if (std::abs(k.size - sz) / std::max(k.size, sz) >= 0.2f)
553 {
554 sz = k.size;
555 keypoints.push_back(k);
556 }
557 }
558 }
559
560 oct++;
561 }
562 }
563
detectAndCompute(InputArray image,InputArray mask,CV_OUT std::vector<Elliptic_KeyPoint> & keypoints,OutputArray descriptors,bool useProvidedKeypoints)564 void TBMR_Impl::detectAndCompute(
565 InputArray image, InputArray mask,
566 CV_OUT std::vector<Elliptic_KeyPoint> &keypoints, OutputArray descriptors,
567 bool useProvidedKeypoints)
568 {
569 // We can use SIFT to compute descriptors for the extracted keypoints...
570 auto sift = SIFT::create();
571 auto dac = AffineFeature2D::create(this, sift);
572 dac->detectAndCompute(image, mask, keypoints, descriptors,
573 useProvidedKeypoints);
574 }
575
create(int _min_area,float _max_area_relative,float _scale,int _n_scale)576 Ptr<TBMR> TBMR::create(int _min_area, float _max_area_relative, float _scale,
577 int _n_scale)
578 {
579 return cv::makePtr<TBMR_Impl>(
580 TBMR_Impl::Params(_min_area, _max_area_relative, _scale, _n_scale));
581 }
582
583 } // namespace xfeatures2d
584 } // namespace cv