1 /** @file mser.c
2 ** @brief MSER - Definition
3 ** @author Andrea Vedaldi
4 **/
5
6 /*
7 Copyright (C) 2007-13 Andrea Vedaldi and Brian Fulkerson.
8 All rights reserved.
9
10 This file is part of the VLFeat library and is made available under
11 the terms of the BSD license (see the COPYING file).
12 */
13
14 /**
15 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
16 @page mser Maximally Stable Extremal Regions (MSER)
17 @author Andrea Vedaldi
18 @tableofcontents
19 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
20
21 @ref mser.h implements the *Maximally Stable Extremal Regions* (MSER)
22 local feature detector of @cite{matas03robust}. This detector extracts
23 as features the the connected components of the level sets of the
24 input intensity image. Among all such regions, the ones that are
25 locally maximally stable are selected. MSERs are affine co-variant, as
26 well as largely co-variant to generic diffeomorphic transformations.
27
28 See @ref mser-starting for an introduction on how to use the detector
29 from the C API. For further details refer to:
30
31 - @subpage mser-fundamentals - MSER definition and parameters.
32
33 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
34 @section mser-starting Getting started with the MSER detector
35 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
36
37 Running the MSER filter usually involves the following steps:
38
39 - Initialize the MSER filter by ::vl_mser_new(). The
40 filter can be reused for images of the same size.
41 - Compute the MSERs by ::vl_mser_process().
42 - Optionally fit ellipsoids to the MSERs by ::vl_mser_ell_fit().
43 - Retrieve the results by ::vl_mser_get_regions() (and optionally ::vl_mser_get_ell()).
44 - Optionally retrieve filter statistics by ::vl_mser_get_stats().
45 - Delete the MSER filter by ::vl_mser_delete().
46
47 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
48 @page mser-fundamentals MSER fundamentals
49 @tableofcontents
50 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
51
52 The *extermal regions* of an image are the connected components of the
53 level sets $S_l = \{ x : I(x) \leq l \}, l \in \real$ of the image
54 $I(x)$. Consider a discretization of the intensity levels $l$
55 consisting of $M$ samples $\mathcal{L}=\{0,\dots,M-1\}$. The extremal
56 regions $R_l \subset S_l$ of the level sets $S_l, l \in \mathcal{L}$
57 can be arranged in a tree, where a region $R_l$ is a children of a
58 region $R_{l+1}$ if $R_l \subset R_{l+1}$. The following figures shows
59 a 1D example where the regions are denoted by dark thick lines:
60
61 @image html mser-tree.png "Connected components of the image level sets arranged in a tree."
62
63 Note that, depending on the image, regions at different levels can be
64 identical as sets:
65
66 @image html mser-er-step.png "Connected components when the image contains step changes."
67
68 A *stable extremal region* is an extremal region that does not change
69 much as the index $l$ is varied. Here we use a criterion which is
70 similar but not identical to the original paper. This definition is
71 somewhat simpler both to understand and code.
72
73 Let $B(R_l)=(R_l,R_{l+1},\dots,R_{l+\Delta})$ be the branch of the
74 tree $R_l \subset R_{l+1} \subset \dots \subset R_{l + \Delta}$
75 rooted at $R_l$. We associate to the branch the (in)stability score
76
77 @f[
78 v(R_l) = \frac{|R_{l+\Delta} - R_l|}{|R_l|}.
79 @f]
80
81 This score is a relative measure of how much $R_l$ changes as the
82 index is increased from $l$ to $l+\Delta$, as illustrated in the
83 following figure.
84
85 @image html mser-er.png "Stability is measured by looking at how much a region changes with the intensity level."
86
87 The score is low if the regions along the branch have similar area
88 (and thus similar shape). We aim to select maximally stable
89 branches; then a maximally stable region is just a representative
90 region selected from a maximally stable branch (for simplicity we
91 select $R_l$, but one could choose for example
92 $R_{l+\Delta/2}$).
93
94 Roughly speaking, a branch is maximally stable if it is a local
95 minimum of the (in)stability score. More accurately, we start by
96 assuming that all branches are maximally stable. Then we consider
97 each branch $B(R_{l})$ and its parent branch
98 $B(R_{l+1}):R_{l+1}\supset R_l$ (notice that, due to the
99 discrete nature of the calculations, they might be geometrically
100 identical) and we mark as unstable the less stable one, i.e.:
101
102 - if $v(R_l)<v(R_{l+1})$, mark $R_{l+1}$ as unstable;
103 - if $v(R_l)>v(R_{l+1})$, mark $R_{l}$ as unstable;
104 - otherwise, do nothing.
105
106 This criterion selects among nearby regions the ones that are more
107 stable. We optionally refine the selection by running (starting
108 from the bigger and going to the smaller regions) the following
109 tests:
110
111 - $a_- \leq |R_{l}|/|R_{\infty}| \leq a_+$: exclude MSERs too
112 small or too big ($|R_{\infty}|$ is the area of the image).
113
114 - $v(R_{l}) < v_+$: exclude MSERs too unstable.
115
116 - For any MSER $R_l$, find the parent MSER $R_{l'}$ and check
117 if
118 $|R_{l'} - R_l|/|R_l'| < d_+$: remove duplicated MSERs.
119
120 <table>
121 <tr>
122 <td>parameter</td>
123 <td>alt. name</td>
124 <td>standard value</td>
125 <td>set by</td>
126 </tr>
127 <tr>
128 <td>$\Delta$</td>
129 <td>@c delta</td>
130 <td>5</td>
131 <td>::vl_mser_set_delta()</td>
132 </tr>
133 <tr>
134 <td>$a_+$</td>
135 <td>@c max_area</td>
136 <td>0.75</td>
137 <td>::vl_mser_set_max_area()</td>
138 </tr>
139 <tr>
140 <td>$a_-$</td>
141 <td>@c min_area</td>
142 <td>3.0/$|R_\infty|$</td>
143 <td>::vl_mser_set_min_area()</td>
144 </tr>
145 <tr>
146 <td>$v_+$</td>
147 <td>@c max_var</td>
148 <td>0.25</td>
149 <td>::vl_mser_set_max_variation()</td>
150 </tr>
151 <tr>
152 <td>$d_+$</td>
153 <td>@c min_diversity</td>
154 <td>0.2</td>
155 <td>::vl_mser_set_min_diversity()</td>
156 </tr>
157 </table>
158
159 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
160 @section mser-vol Volumetric images
161 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
162
163 The code supports images of arbitrary dimension. For instance, it
164 is possible to find the MSER regions of volumetric images or time
165 sequences. See ::vl_mser_new() for further details
166
167 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
168 @section mser-ell Ellipsoids
169 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
170
171 Usually extremal regions are returned as a set of ellipsoids
172 fitted to the actual regions (which have arbitrary shape). The fit
173 is done by calculating the mean and variance of the pixels
174 composing the region:
175 @f[
176 \mu_l = \frac{1}{|R_l|}\sum_{x\in R_l}x,
177 \qquad
178 \Sigma_l = \frac{1}{|R_l|}\sum_{x\in R_l} (x-\mu_l)^\top(x-\mu_l)
179 @f]
180 Ellipsoids are fitted by ::vl_mser_ell_fit(). Notice that for a
181 <em>n</em> dimensional image, the mean has <em>n</em> components
182 and the variance has <em>n(n+1)/2</em> independent components. The
183 total number of components is obtained by ::vl_mser_get_ell_dof()
184 and the total number of fitted ellipsoids by
185 ::vl_mser_get_ell_num(). A matrix with an ellipsoid per column is
186 returned by ::vl_mser_get_ell(). The column is the stacking of the
187 mean and of the independent components of the variance, in the
188 order <em>(1,1),(1,2),..,(1,n), (2,2),(2,3)...</em>. In the
189 calculations, the pixel coordinate $x=(x_1,...,x_n)$ use the
190 standard index order and ranges.
191
192 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
193 @section mser-algo Algorithm
194 <!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
195
196 The algorithm is quite efficient. While some details may be
197 tricky, the overall idea is easy to grasp.
198
199 - Pixels are sorted by increasing intensity.
200 - Pixels are added to a forest by increasing intensity. The forest has the
201 following properties:
202 - All the descendent of a certain pixels are subset of an extremal region.
203 - All the extremal regions are the descendants of some pixels.
204 - Extremal regions are extracted from the region tree and the extremal regions tree is
205 calculated.
206 - Stable regions are marked.
207 - Duplicates and other bad regions are removed.
208
209 @remark The extremal region tree which is calculated is a subset
210 of the actual extremal region tree. In particular, it does not
211 contain redundant entries extremal regions that coincide as
212 sets. So, for example, in the calculated extremal region tree, the
213 parent $R_q$ of an extremal region $R_{l}$ may or may
214 <em>not</em> correspond to $R_{l+1}$, depending whether
215 $q\leq l+1$ or not. These subtleties are important when
216 calculating the stability tests.
217
218 **/
219
220 #include "mser.h"
221 #include<stdlib.h>
222 #include<string.h>
223 #include<assert.h>
224
225 /** -------------------------------------------------------------------
226 ** @brief Advance N-dimensional subscript
227 **
228 ** The function increments by one the subscript @a subs indexing an
229 ** array the @a ndims dimensions @a dims.
230 **
231 ** @param ndims number of dimensions.
232 ** @param dims dimensions.
233 ** @param subs subscript to advance.
234 **/
235
236 VL_INLINE void
adv(int ndims,int const * dims,int * subs)237 adv(int ndims, int const *dims, int *subs)
238 {
239 int d = 0 ;
240 while(d < ndims) {
241 if( ++subs[d] < dims[d] ) return ;
242 subs[d++] = 0 ;
243 }
244 }
245
246 /** -------------------------------------------------------------------
247 ** @brief Climb the region forest to reach aa root
248 **
249 ** The function climbs the regions forest @a r starting from the node
250 ** @a idx to the corresponding root.
251 **
252 ** To speed-up the operation, the function uses the
253 ** VlMserReg::shortcut field to quickly jump to the root. After the
254 ** root is reached, all the used shortcut are updated.
255 **
256 ** @param r regions' forest.
257 ** @param idx stating node.
258 ** @return index of the reached root.
259 **/
260
261 VL_INLINE vl_uint
climb(VlMserReg * r,vl_uint idx)262 climb (VlMserReg* r, vl_uint idx)
263 {
264
265 vl_uint prev_idx = idx ;
266 vl_uint next_idx ;
267 vl_uint root_idx ;
268
269 /* move towards root to find it */
270 while (1) {
271
272 /* next jump to the root */
273 next_idx = r [idx] .shortcut ;
274
275 /* recycle shortcut to remember how we came here */
276 r [idx] .shortcut = prev_idx ;
277
278 /* stop if the root is found */
279 if( next_idx == idx ) break ;
280
281 /* next guy */
282 prev_idx = idx ;
283 idx = next_idx ;
284 }
285
286 root_idx = idx ;
287
288 /* move backward to update shortcuts */
289 while (1) {
290
291 /* get previously visited one */
292 prev_idx = r [idx] .shortcut ;
293
294 /* update shortcut to point to the new root */
295 r [idx] .shortcut = root_idx ;
296
297 /* stop if the first visited node is reached */
298 if( prev_idx == idx ) break ;
299
300 /* next guy */
301 idx = prev_idx ;
302 }
303
304 return root_idx ;
305 }
306
307 /** -------------------------------------------------------------------
308 ** @brief Create a new MSER filter
309 **
310 ** Initializes a new MSER filter for images of the specified
311 ** dimensions. Images are @a ndims -dimensional arrays of dimensions
312 ** @a dims.
313 **
314 ** @param ndims number of dimensions.
315 ** @param dims dimensions.
316 **/
317 VL_EXPORT
318 VlMserFilt*
vl_mser_new(int ndims,int const * dims)319 vl_mser_new (int ndims, int const* dims)
320 {
321 VlMserFilt* f ;
322 int *strides, k ;
323
324 f = vl_calloc (sizeof(VlMserFilt), 1) ;
325
326 f-> ndims = ndims ;
327 f-> dims = vl_malloc (sizeof(int) * ndims) ;
328 f-> subs = vl_malloc (sizeof(int) * ndims) ;
329 f-> dsubs = vl_malloc (sizeof(int) * ndims) ;
330 f-> strides = vl_malloc (sizeof(int) * ndims) ;
331
332 /* shortcuts */
333 strides = f-> strides ;
334
335 /* copy dims to f->dims */
336 for(k = 0 ; k < ndims ; ++k) {
337 f-> dims [k] = dims [k] ;
338 }
339
340 /* compute strides to move into the N-dimensional image array */
341 strides [0] = 1 ;
342 for(k = 1 ; k < ndims ; ++k) {
343 strides [k] = strides [k-1] * dims [k-1] ;
344 }
345
346 /* total number of pixels */
347 f-> nel = strides [ndims-1] * dims [ndims-1] ;
348
349 /* dof of ellipsoids */
350 f-> dof = ndims * (ndims + 1) / 2 + ndims ;
351
352 /* more buffers */
353 f-> perm = vl_malloc (sizeof(vl_uint) * f-> nel) ;
354 f-> joins = vl_malloc (sizeof(vl_uint) * f-> nel) ;
355 f-> r = vl_malloc (sizeof(VlMserReg) * f-> nel) ;
356
357 f-> er = 0 ;
358 f-> rer = 0 ;
359 f-> mer = 0 ;
360 f-> rmer = 0 ;
361 f-> ell = 0 ;
362 f-> rell = 0 ;
363
364 /* other parameters */
365 f-> delta = 5 ;
366 f-> max_area = 0.75 ;
367 f-> min_area = 3.0 / f-> nel ;
368 f-> max_variation = 0.25 ;
369 f-> min_diversity = 0.2 ;
370
371 return f ;
372 }
373
374 /** -------------------------------------------------------------------
375 ** @brief Delete MSER filter
376 **
377 ** The function releases the MSER filter @a f and all its resources.
378 **
379 ** @param f MSER filter to be deleted.
380 **/
381 VL_EXPORT
382 void
vl_mser_delete(VlMserFilt * f)383 vl_mser_delete (VlMserFilt* f)
384 {
385 if(f) {
386 if(f-> acc ) vl_free( f-> acc ) ;
387 if(f-> ell ) vl_free( f-> ell ) ;
388
389 if(f-> er ) vl_free( f-> er ) ;
390 if(f-> r ) vl_free( f-> r ) ;
391 if(f-> joins ) vl_free( f-> joins ) ;
392 if(f-> perm ) vl_free( f-> perm ) ;
393
394 if(f-> strides) vl_free( f-> strides) ;
395 if(f-> dsubs ) vl_free( f-> dsubs ) ;
396 if(f-> subs ) vl_free( f-> subs ) ;
397 if(f-> dims ) vl_free( f-> dims ) ;
398
399 if(f-> mer ) vl_free( f-> mer ) ;
400 vl_free (f) ;
401 }
402 }
403
404
405 /** -------------------------------------------------------------------
406 ** @brief Process image
407 **
408 ** The functions calculates the Maximally Stable Extremal Regions
409 ** (MSERs) of image @a im using the MSER filter @a f.
410 **
411 ** The filter @a f must have been initialized to be compatible with
412 ** the dimensions of @a im.
413 **
414 ** @param f MSER filter.
415 ** @param im image data.
416 **/
417 VL_EXPORT
418 void
vl_mser_process(VlMserFilt * f,vl_mser_pix const * im)419 vl_mser_process (VlMserFilt* f, vl_mser_pix const* im)
420 {
421 /* shortcuts */
422 vl_uint nel = f-> nel ;
423 vl_uint *perm = f-> perm ;
424 vl_uint *joins = f-> joins ;
425 int ndims = f-> ndims ;
426 int *dims = f-> dims ;
427 int *subs = f-> subs ;
428 int *dsubs = f-> dsubs ;
429 int *strides = f-> strides ;
430 VlMserReg *r = f-> r ;
431 VlMserExtrReg *er = f-> er ;
432 vl_uint *mer = f-> mer ;
433 int delta = f-> delta ;
434
435 int njoins = 0 ;
436 int ner = 0 ;
437 int nmer = 0 ;
438 int nbig = 0 ;
439 int nsmall = 0 ;
440 int nbad = 0 ;
441 int ndup = 0 ;
442
443 int i, j, k ;
444
445 /* delete any previosuly computed ellipsoid */
446 f-> nell = 0 ;
447
448 /* -----------------------------------------------------------------
449 * Sort pixels by intensity
450 * -------------------------------------------------------------- */
451
452 {
453 vl_uint buckets [ VL_MSER_PIX_MAXVAL ] ;
454
455 /* clear buckets */
456 memset (buckets, 0, sizeof(vl_uint) * VL_MSER_PIX_MAXVAL ) ;
457
458 /* compute bucket size (how many pixels for each intensity
459 value) */
460 for(i = 0 ; i < (int) nel ; ++i) {
461 vl_mser_pix v = im [i] ;
462 ++ buckets [v] ;
463 }
464
465 /* cumulatively add bucket sizes */
466 for(i = 1 ; i < VL_MSER_PIX_MAXVAL ; ++i) {
467 buckets [i] += buckets [i-1] ;
468 }
469
470 /* empty buckets computing pixel ordering */
471 for(i = nel ; i >= 1 ; ) {
472 vl_mser_pix v = im [ --i ] ;
473 vl_uint j = -- buckets [v] ;
474 perm [j] = i ;
475 }
476 }
477
478 /* initialize the forest with all void nodes */
479 for(i = 0 ; i < (int) nel ; ++i) {
480 r [i] .parent = VL_MSER_VOID_NODE ;
481 }
482
483 /* -----------------------------------------------------------------
484 * Compute regions and count extremal regions
485 * -------------------------------------------------------------- */
486 /*
487 In the following:
488
489 idx : index of the current pixel
490 val : intensity of the current pixel
491 r_idx : index of the root of the current pixel
492 n_idx : index of the neighbors of the current pixel
493 nr_idx : index of the root of the neighbor of the current pixel
494
495 */
496
497 /* process each pixel by increasing intensity */
498 for(i = 0 ; i < (int) nel ; ++i) {
499
500 /* pop next node xi */
501 vl_uint idx = perm [i] ;
502 vl_mser_pix val = im [idx] ;
503 vl_uint r_idx ;
504
505 /* add the pixel to the forest as a root for now */
506 r [idx] .parent = idx ;
507 r [idx] .shortcut = idx ;
508 r [idx] .area = 1 ;
509 r [idx] .height = 1 ;
510
511 r_idx = idx ;
512
513 /* convert the index IDX into the subscript SUBS; also initialize
514 DSUBS to (-1,-1,...,-1) */
515 {
516 vl_uint temp = idx ;
517 for(k = ndims - 1 ; k >= 0 ; --k) {
518 dsubs [k] = -1 ;
519 subs [k] = temp / strides [k] ;
520 temp = temp % strides [k] ;
521 }
522 }
523
524 /* examine the neighbors of the current pixel */
525 while (1) {
526 vl_uint n_idx = 0 ;
527 vl_bool good = 1 ;
528
529 /*
530 Compute the neighbor subscript as NSUBS+SUB, the
531 corresponding neighbor index NINDEX and check that the
532 neighbor is within the image domain.
533 */
534 for(k = 0 ; k < ndims && good ; ++k) {
535 int temp = dsubs [k] + subs [k] ;
536 good &= (0 <= temp) && (temp < dims [k]) ;
537 n_idx += temp * strides [k] ;
538 }
539
540 /*
541 The neighbor should be processed if the following conditions
542 are met:
543
544 1. The neighbor is within image boundaries.
545
546 2. The neighbor is indeed different from the current node
547 (the opposite happens when DSUB=(0,0,...,0)).
548
549 3. The neighbor is already in the forest, meaning that it has
550 already been processed.
551 */
552 if (good &&
553 n_idx != idx &&
554 r [n_idx] .parent != VL_MSER_VOID_NODE ) {
555
556 vl_mser_pix nr_val = 0 ;
557 vl_uint nr_idx = 0 ;
558 int hgt = r [ r_idx] .height ;
559 int n_hgt = r [nr_idx] .height ;
560
561 /*
562 Now we join the two subtrees rooted at
563
564 R_IDX = ROOT( IDX)
565 NR_IDX = ROOT(N_IDX).
566
567 Note that R_IDX = ROOT(IDX) might change as we process more
568 neighbors, so we need keep updating it.
569 */
570
571 r_idx = climb(r, idx) ;
572 nr_idx = climb(r, n_idx) ;
573
574 /*
575 At this point we have three possibilities:
576
577 (A) ROOT(IDX) == ROOT(NR_IDX). In this case the two trees
578 have already been joined and we do not do anything.
579
580 (B) I(ROOT(IDX)) == I(ROOT(NR_IDX)). In this case the pixel
581 IDX is extending an extremal region with the same
582 intensity value. Since ROOT(NR_IDX) will NOT be an
583 extremal region of the full image, ROOT(IDX) can be
584 safely added as children of ROOT(NR_IDX) if this
585 reduces the height according to the union rank
586 heuristic.
587
588 (C) I(ROOT(IDX)) > I(ROOT(NR_IDX)). In this case the pixel
589 IDX is starting a new extremal region. Thus ROOT(NR_IDX)
590 WILL be an extremal region of the final image and the
591 only possibility is to add ROOT(NR_IDX) as children of
592 ROOT(IDX), which becomes parent.
593 */
594
595 if( r_idx != nr_idx ) { /* skip if (A) */
596
597 nr_val = im [nr_idx] ;
598
599 if( nr_val == val && hgt < n_hgt ) {
600
601 /* ROOT(IDX) becomes the child */
602 r [r_idx] .parent = nr_idx ;
603 r [r_idx] .shortcut = nr_idx ;
604 r [nr_idx] .area += r [r_idx] .area ;
605 r [nr_idx] .height = VL_MAX(n_hgt, hgt+1) ;
606
607 joins [njoins++] = r_idx ;
608
609 } else {
610
611 /* cases ROOT(IDX) becomes the parent */
612 r [nr_idx] .parent = r_idx ;
613 r [nr_idx] .shortcut = r_idx ;
614 r [r_idx] .area += r [nr_idx] .area ;
615 r [r_idx] .height = VL_MAX(hgt, n_hgt + 1) ;
616
617 joins [njoins++] = nr_idx ;
618
619 /* count if extremal */
620 if (nr_val != val) ++ ner ;
621
622 } /* check b vs c */
623 } /* check a vs b or c */
624 } /* neighbor done */
625
626 /* move to next neighbor */
627 k = 0 ;
628 while(++ dsubs [k] > 1) {
629 dsubs [k++] = -1 ;
630 if(k == ndims) goto done_all_neighbors ;
631 }
632 } /* next neighbor */
633 done_all_neighbors : ;
634 } /* next pixel */
635
636 /* the last root is extremal too */
637 ++ ner ;
638
639 /* save back */
640 f-> njoins = njoins ;
641
642 f-> stats. num_extremal = ner ;
643
644 /* -----------------------------------------------------------------
645 * Extract extremal regions
646 * -------------------------------------------------------------- */
647
648 /*
649 Extremal regions are extracted and stored into the array ER. The
650 structure R is also updated so that .SHORTCUT indexes the
651 corresponding extremal region if any (otherwise it is set to
652 VOID).
653 */
654
655 /* make room */
656 if (f-> rer < ner) {
657 if (er) vl_free (er) ;
658 f->er = er = vl_malloc (sizeof(VlMserExtrReg) * ner) ;
659 f->rer = ner ;
660 } ;
661
662 /* save back */
663 f-> nmer = ner ;
664
665 /* count again */
666 ner = 0 ;
667
668 /* scan all regions Xi */
669 for(i = 0 ; i < (int) nel ; ++i) {
670
671 /* pop next node xi */
672 vl_uint idx = perm [i] ;
673
674 vl_mser_pix val = im [idx] ;
675 vl_uint p_idx = r [idx] .parent ;
676 vl_mser_pix p_val = im [p_idx] ;
677
678 /* is extremal ? */
679 vl_bool is_extr = (p_val > val) || idx == p_idx ;
680
681 if( is_extr ) {
682
683 /* if so, add it */
684 er [ner] .index = idx ;
685 er [ner] .parent = ner ;
686 er [ner] .value = im [idx] ;
687 er [ner] .area = r [idx] .area ;
688
689 /* link this region to this extremal region */
690 r [idx] .shortcut = ner ;
691
692 /* increase count */
693 ++ ner ;
694 } else {
695 /* link this region to void */
696 r [idx] .shortcut = VL_MSER_VOID_NODE ;
697 }
698 }
699
700 /* -----------------------------------------------------------------
701 * Link extremal regions in a tree
702 * -------------------------------------------------------------- */
703
704 for(i = 0 ; i < ner ; ++i) {
705
706 vl_uint idx = er [i] .index ;
707
708 do {
709 idx = r[idx] .parent ;
710 } while (r[idx] .shortcut == VL_MSER_VOID_NODE) ;
711
712 er[i] .parent = r[idx] .shortcut ;
713 er[i] .shortcut = i ;
714 }
715
716 /* -----------------------------------------------------------------
717 * Compute variability of +DELTA branches
718 * -------------------------------------------------------------- */
719 /* For each extremal region Xi of value VAL we look for the biggest
720 * parent that has value not greater than VAL+DELTA. This is dubbed
721 * `top parent'. */
722
723 for(i = 0 ; i < ner ; ++i) {
724
725 /* Xj is the current region the region and Xj are the parents */
726 int top_val = er [i] .value + delta ;
727 int top = er [i] .shortcut ;
728
729 /* examine all parents */
730 while (1) {
731 int next = er [top] .parent ;
732 int next_val = er [next] .value ;
733
734 /* Break if:
735 * - there is no node above the top or
736 * - the next node is above the top value.
737 */
738 if (next == top || next_val > top_val) break ;
739
740 /* so next could be the top */
741 top = next ;
742 }
743
744 /* calculate branch variation */
745 {
746 int area = er [i ] .area ;
747 int area_top = er [top] .area ;
748 er [i] .variation = (float) (area_top - area) / area ;
749 er [i] .max_stable = 1 ;
750 }
751
752 /* Optimization: since extremal regions are processed by
753 * increasing intensity, all next extremal regions being processed
754 * have value at least equal to the one of Xi. If any of them has
755 * parent the parent of Xi (this comprises the parent itself), we
756 * can safely skip most intermediate node along the branch and
757 * skip directly to the top to start our search. */
758 {
759 int parent = er [i] .parent ;
760 int curr = er [parent] .shortcut ;
761 er [parent] .shortcut = VL_MAX (top, curr) ;
762 }
763 }
764
765 /* -----------------------------------------------------------------
766 * Select maximally stable branches
767 * -------------------------------------------------------------- */
768
769 nmer = ner ;
770 for(i = 0 ; i < ner ; ++i) {
771 vl_uint parent = er [i ] .parent ;
772 vl_mser_pix val = er [i ] .value ;
773 float var = er [i ] .variation ;
774 vl_mser_pix p_val = er [parent] .value ;
775 float p_var = er [parent] .variation ;
776 vl_uint loser ;
777
778 /*
779 Notice that R_parent = R_{l+1} only if p_val = val + 1. If not,
780 this and the parent region coincide and there is nothing to do.
781 */
782 if(p_val > val + 1) continue ;
783
784 /* decide which one to keep and put that in loser */
785 if(var < p_var) loser = parent ; else loser = i ;
786
787 /* make loser NON maximally stable */
788 if(er [loser] .max_stable) {
789 -- nmer ;
790 er [loser] .max_stable = 0 ;
791 }
792 }
793
794 f-> stats. num_unstable = ner - nmer ;
795
796 /* -----------------------------------------------------------------
797 * Further filtering
798 * -------------------------------------------------------------- */
799 /* It is critical for correct duplicate detection to remove regions
800 * from the bottom (smallest one first). */
801 {
802 float max_area = (float) f-> max_area * nel ;
803 float min_area = (float) f-> min_area * nel ;
804 float max_var = (float) f-> max_variation ;
805 float min_div = (float) f-> min_diversity ;
806
807 /* scan all extremal regions (intensity value order) */
808 for(i = ner-1 ; i >= 0L ; --i) {
809
810 /* process only maximally stable extremal regions */
811 if (! er [i] .max_stable) continue ;
812
813 if (er [i] .variation >= max_var ) { ++ nbad ; goto remove ; }
814 if (er [i] .area > max_area) { ++ nbig ; goto remove ; }
815 if (er [i] .area < min_area) { ++ nsmall ; goto remove ; }
816
817 /*
818 * Remove duplicates
819 */
820 if (min_div < 1.0) {
821 vl_uint parent = er [i] .parent ;
822 int area, p_area ;
823 float div ;
824
825 /* check all but the root mser */
826 if((int) parent != i) {
827
828 /* search for the maximally stable parent region */
829 while(! er [parent] .max_stable) {
830 vl_uint next = er [parent] .parent ;
831 if(next == parent) break ;
832 parent = next ;
833 }
834
835 /* Compare with the parent region; if the current and parent
836 * regions are too similar, keep only the parent. */
837 area = er [i] .area ;
838 p_area = er [parent] .area ;
839 div = (float) (p_area - area) / (float) p_area ;
840
841 if (div < min_div) { ++ ndup ; goto remove ; }
842 } /* remove dups end */
843
844 }
845 continue ;
846 remove :
847 er [i] .max_stable = 0 ;
848 -- nmer ;
849 } /* check next region */
850
851 f-> stats .num_abs_unstable = nbad ;
852 f-> stats .num_too_big = nbig ;
853 f-> stats .num_too_small = nsmall ;
854 f-> stats .num_duplicates = ndup ;
855 }
856 /* -----------------------------------------------------------------
857 * Save the result
858 * -------------------------------------------------------------- */
859
860 /* make room */
861 if (f-> rmer < nmer) {
862 if (mer) vl_free (mer) ;
863 f->mer = mer = vl_malloc( sizeof(vl_uint) * nmer) ;
864 f->rmer = nmer ;
865 }
866
867 /* save back */
868 f-> nmer = nmer ;
869
870 j = 0 ;
871 for (i = 0 ; i < ner ; ++i) {
872 if (er [i] .max_stable) mer [j++] = er [i] .index ;
873 }
874 }
875
876 /** -------------------------------------------------------------------
877 ** @brief Fit ellipsoids
878 **
879 ** @param f MSER filter.
880 **
881 ** @sa @ref mser-ell
882 **/
883
884 VL_EXPORT
885 void
vl_mser_ell_fit(VlMserFilt * f)886 vl_mser_ell_fit (VlMserFilt* f)
887 {
888 /* shortcuts */
889 int nel = f-> nel ;
890 int dof = f-> dof ;
891 int *dims = f-> dims ;
892 int ndims = f-> ndims ;
893 int *subs = f-> subs ;
894 int njoins = f-> njoins ;
895 vl_uint *joins = f-> joins ;
896 VlMserReg *r = f-> r ;
897 vl_uint *mer = f-> mer ;
898 int nmer = f-> nmer ;
899 vl_mser_acc *acc = f-> acc ;
900 vl_mser_acc *ell = f-> ell ;
901
902 int d, index, i, j ;
903
904 /* already fit ? */
905 if (f->nell == f->nmer) return ;
906
907 /* make room */
908 if (f->rell < f->nmer) {
909 if (f->ell) vl_free (f->ell) ;
910 f->ell = vl_malloc (sizeof(float) * f->nmer * f->dof) ;
911 f->rell = f-> nmer ;
912 }
913
914 if (f->acc == 0) {
915 f->acc = vl_malloc (sizeof(float) * f->nel) ;
916 }
917
918 acc = f-> acc ;
919 ell = f-> ell ;
920
921 /* -----------------------------------------------------------------
922 * Integrate moments
923 * -------------------------------------------------------------- */
924
925 /* for each dof */
926 for(d = 0 ; d < f->dof ; ++d) {
927
928 /* start from the upper-left pixel (0,0,...,0) */
929 memset (subs, 0, sizeof(int) * ndims) ;
930
931 /* step 1: fill acc pretending that each region has only one pixel */
932 if(d < ndims) {
933 /* 1-order ................................................... */
934
935 for(index = 0 ; index < nel ; ++ index) {
936 acc [index] = subs [d] ;
937 adv(ndims, dims, subs) ;
938 }
939 }
940 else {
941 /* 2-order ................................................... */
942
943 /* map the dof d to a second order moment E[x_i x_j] */
944 i = d - ndims ;
945 j = 0 ;
946 while(i > j) {
947 i -= j + 1 ;
948 j ++ ;
949 }
950 /* initialize acc with x_i * x_j */
951 for(index = 0 ; index < nel ; ++ index){
952 acc [index] = subs [i] * subs [j] ;
953 adv(ndims, dims, subs) ;
954 }
955 }
956
957 /* step 2: integrate */
958 for(i = 0 ; i < njoins ; ++i) {
959 vl_uint index = joins [i] ;
960 vl_uint parent = r [index] .parent ;
961 acc [parent] += acc [index] ;
962 }
963
964 /* step 3: save back to ellpises */
965 for(i = 0 ; i < nmer ; ++i) {
966 vl_uint idx = mer [i] ;
967 ell [d + dof*i] = acc [idx] ;
968 }
969
970 } /* next dof */
971
972 /* -----------------------------------------------------------------
973 * Compute central moments
974 * -------------------------------------------------------------- */
975
976 for(index = 0 ; index < nmer ; ++index) {
977 float *pt = ell + index * dof ;
978 vl_uint idx = mer [index] ;
979 float area = r [idx] .area ;
980
981 for(d = 0 ; d < dof ; ++d) {
982
983 pt [d] /= area ;
984
985 if(d >= ndims) {
986 /* remove squared mean from moment to get variance */
987 i = d - ndims ;
988 j = 0 ;
989 while(i > j) {
990 i -= j + 1 ;
991 j ++ ;
992 }
993 pt [d] -= pt [i] * pt [j] ;
994 }
995
996 }
997 }
998
999 /* save back */
1000 f-> nell = nmer ;
1001 }
1002