summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/mser/src/matlab/mser.mex.c')
-rwxr-xr-xSD-VBS/benchmarks/mser/src/matlab/mser.mex.c815
1 files changed, 0 insertions, 815 deletions
diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c b/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c
deleted file mode 100755
index 8473afe..0000000
--- a/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c
+++ /dev/null
@@ -1,815 +0,0 @@
1/* file: mser.mex.c
2** description: Maximally Stable Extremal Regions
3** author: Andrea Vedaldi
4**/
5
6/* AUTORIGHTS
7Copyright (C) 2006 Regents of the University of California
8All rights reserved
9
10Written by Andrea Vedaldi (UCLA VisionLab).
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met
14
15 * Redistributions of source code must retain the above copyright
16 notice, this list of conditions and the following disclaimer.
17 * Redistributions in binary form must reproduce the above copyright
18 notice, this list of conditions and the following disclaimer in the
19 documentation and/or other materials provided with the distribution.
20 * Neither the name of the University of California, Berkeley nor the
21 names of its contributors may be used to endorse or promote products
22 derived from this software without specific prior written permission.
23
24THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
25EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
26WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
27DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
28DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
31ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
32(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34*/
35
36/** @file
37 ** @brief Maximally Stable Extremal Regions - MEX implementation
38 **/
39
40#include<mexutils.c>
41#include<stdio.h>
42#include<stdlib.h>
43#include<math.h>
44#include<string.h>
45#include<assert.h>
46
47#define MIN(x,y) (((x)<(y))?(x):(y))
48#define MAX(x,y) (((x)>(y))?(x):(y))
49
50#define BUCKETS 256
51
52#define USE_BUCKET_SORT
53/*#define USE_RANK_UNION
54*/
55
56typedef char unsigned val_t ;
57typedef int unsigned idx_t ;
58typedef long long int unsigned acc_t ;
59
60/* pairs are used to sort the pixels */
61typedef struct
62{
63 val_t value ;
64 idx_t index ;
65} pair_t ;
66
67/* forest node */
68typedef struct
69{
70 idx_t parent ; /**< parent pixel */
71 idx_t shortcut ; /**< shortcut to the root */
72 idx_t region ; /**< index of the region */
73 int area ; /**< area of the region */
74#ifdef USE_RANK_UNION
75 int height ; /**< node height */
76#endif
77} node_t ;
78
79/* extremal regions */
80typedef struct
81{
82 idx_t parent ; /**< parent region */
83 idx_t index ; /**< index of root pixel */
84 val_t value ; /**< value of root pixel */
85 int area ; /**< area of the region */
86 int area_top ; /**< area of the region DELTA levels above */
87 int area_bot ; /**< area of the region DELTA levels below */
88 float variation ; /**< variation */
89 int maxstable ; /**< max stable number (=0 if not maxstable) */
90} region_t ;
91
92/* predicate used to sort pixels by increasing intensity */
93int
94cmp_pair(void const* a, void const* b)
95{
96 pair_t* pa = (pair_t*) a;
97 pair_t* pb = (pair_t*) b;
98 return pa->value - pb->value ;
99}
100
101/* advance N-dimensional subscript */
102void
103adv(int const* dims, int ndims, int* subs_pt)
104{
105 int d = 0 ;
106 while(d < ndims) {
107 if( ++subs_pt[d] < dims[d] ) return ;
108 subs_pt[d++] = 0 ;
109 }
110}
111
112/* driver */
113void
114mexFunction(int nout, mxArray *out[],
115 int nin, const mxArray *in[])
116{
117 enum {IN_I=0, IN_DELTA} ;
118 enum {OUT_REGIONS=0, OUT_ELL, OUT_PARENTS, OUT_AREA} ;
119
120 idx_t i ;
121 idx_t rindex = 0 ;
122 int k ;
123
124 /* configuration */
125 int verbose = 0 ; /* be verbose */
126 int small_cleanup= 1 ; /* remove very small regions */
127 int big_cleanup = 1 ; /* remove very big regions */
128 int bad_cleanup = 0 ; /* remove very bad regions */
129 int dup_cleanup = 1 ; /* remove duplicates */
130 val_t delta ; /* stability delta */
131
132 /* node value denoting a void node */
133 idx_t const node_is_void = 0xffffffff ;
134
135 int* subs_pt ; /* N-dimensional subscript */
136 int* nsubs_pt ; /* diff-subscript to point to neigh. */
137 idx_t* strides_pt ; /* strides to move in image array */
138 idx_t* visited_pt ; /* flag */
139
140 int nel ; /* number of image elements (pixels) */
141 int ner = 0 ; /* number of extremal regions */
142 int nmer = 0 ; /* number of maximally stable */
143 int ndims ; /* number of dimensions */
144 int const* dims ; /* dimensions */
145 int njoins = 0 ; /* number of join ops */
146
147 val_t const* I_pt ; /* source image */
148 pair_t* pairs_pt ; /* scratch buffer to sort pixels */
149 node_t* forest_pt ; /* the extremal regions forest */
150 region_t* regions_pt ; /* list of extremal regions found */
151
152 /* ellipses fitting */
153 acc_t* acc_pt ; /* accumulator to integrate region moments */
154 acc_t* ell_pt ; /* ellipses parameters */
155 int gdl ; /* number of parameters of an ellipse */
156 idx_t* joins_pt ; /* sequence of joins */
157
158 /** -----------------------------------------------------------------
159 ** Check the arguments
160 ** -------------------------------------------------------------- */
161 if (nin != 2) {
162 mexErrMsgTxt("Two arguments required.") ;
163 } else if (nout > 4) {
164 mexErrMsgTxt("Too many output arguments.");
165 }
166
167 if(mxGetClassID(in[IN_I]) != mxUINT8_CLASS) {
168 mexErrMsgTxt("I must be of class UINT8") ;
169 }
170
171 if(!uIsScalar(in[IN_DELTA])) {
172 mexErrMsgTxt("DELTA must be scalar") ;
173 }
174
175 delta = 0 ;
176 switch(mxGetClassID(in[IN_DELTA])) {
177 case mxUINT8_CLASS :
178 delta = * (val_t*) mxGetData(in[IN_DELTA]) ;
179 break ;
180
181 case mxDOUBLE_CLASS :
182 {
183 double x = *mxGetPr(in[IN_DELTA]) ;
184 if(x < 0.0) {
185 mexErrMsgTxt("DELTA must be non-negative") ;
186 }
187 delta = (val_t) x ;
188 }
189 break ;
190
191 default :
192 mexErrMsgTxt("DELTA must be of class DOUBLE or UINT8") ;
193 }
194
195 /* get dimensions */
196 nel = mxGetNumberOfElements(in[IN_I]) ;
197 ndims = mxGetNumberOfDimensions(in[IN_I]) ;
198 dims = mxGetDimensions(in[IN_I]) ;
199 I_pt = mxGetData(in[IN_I]) ;
200
201 /* allocate stuff */
202 subs_pt = mxMalloc( sizeof(int) * ndims ) ;
203 nsubs_pt = mxMalloc( sizeof(int) * ndims ) ;
204 strides_pt = mxMalloc( sizeof(idx_t) * ndims ) ;
205 visited_pt = mxMalloc( sizeof(idx_t) * nel ) ;
206 regions_pt = mxMalloc( sizeof(region_t) * nel ) ;
207 pairs_pt = mxMalloc( sizeof(pair_t) * nel ) ;
208 forest_pt = mxMalloc( sizeof(node_t) * nel ) ;
209 joins_pt = mxMalloc( sizeof(idx_t) * nel ) ;
210
211 /* compute strides to move into the N-dimensional image array */
212 strides_pt [0] = 1 ;
213 for(k = 1 ; k < ndims ; ++k) {
214 strides_pt [k] = strides_pt [k-1] * dims [k-1] ;
215 }
216
217 /* sort pixels by increasing intensity*/
218 verbose && mexPrintf("Sorting pixels ... ") ;
219
220#ifndef USE_BUCKET_SORT
221 for(i = 0 ; i < nel ; ++i) {
222 pairs_pt [i].value = I_pt [i] ;
223 pairs_pt [i].index = i ;
224 }
225 qsort(pairs_pt, nel, sizeof(pair_t), cmp_pair) ;
226#else
227 {
228 int unsigned buckets [BUCKETS] ;
229 int unsigned v;
230 memset(buckets, 0, sizeof(int unsigned)*BUCKETS) ;
231 for(i = 0 ; i < nel ; ++i) {
232 v = (unsigned int)I_pt [i] ;
233 ++buckets[v] ;
234 }
235 for(i = 1 ; i < BUCKETS ; ++i) {
236 buckets[i] += buckets[i-1] ;
237 }
238 for(i = nel ; i >= 1 ; ) {
239 v = I_pt [--i] ;
240 idx_t j = --buckets [v] ;
241 pairs_pt [j].value = v ;
242 pairs_pt [j].index = i ;
243 }
244 }
245#endif
246 verbose && mexPrintf("done\n") ;
247
248 /* initialize the forest with all void nodes */
249 for(i = 0 ; i < nel ; ++i) {
250 forest_pt [i].parent = node_is_void ;
251 }
252
253 /* number of ellipse free parameters */
254 gdl = ndims*(ndims+1)/2 + ndims ;
255
256 /* -----------------------------------------------------------------
257 * Compute extremal regions tree
258 * -------------------------------------------------------------- */
259 verbose && mexPrintf("Computing extremal regions ... ") ;
260 for(i = 0 ; i < nel ; ++i) {
261
262 /* pop next node xi */
263 idx_t index = pairs_pt [i].index ;
264 val_t value = pairs_pt [i].value ;
265
266
267 /* this will be needed later */
268 rindex = index ;
269
270 /* push it into the tree */
271 forest_pt [index] .parent = index ;
272 forest_pt [index] .shortcut = index ;
273 forest_pt [index] .area = 1 ;
274#ifdef USE_RANK_UNION
275 forest_pt [index] .height = 1 ;
276#endif
277
278 /* convert index into a subscript sub; also initialize nsubs
279 to (-1,-1,...,-1) */
280 {
281 idx_t temp = index ;
282 for(k = ndims-1 ; k >=0 ; --k) {
283 nsubs_pt [k] = -1 ;
284 subs_pt [k] = temp / strides_pt [k] ;
285 temp = temp % strides_pt [k] ;
286 }
287 }
288
289 /* process neighbors of xi */
290 while( true ) {
291 int good = true ;
292 idx_t nindex = 0 ;
293
294 /* compute NSUBS+SUB, the correspoinding neighbor index NINDEX
295 and check that the pixel is within image boundaries. */
296 for(k = 0 ; k < ndims && good ; ++k) {
297 int temp = nsubs_pt [k] + subs_pt [k] ;
298 good &= 0 <= temp && temp < dims[k] ;
299 nindex += temp * strides_pt [k] ;
300 }
301
302 /* keep going only if
303 1 - the neighbor is within image boundaries;
304 2 - the neighbor is indeed different from the current node
305 (this happens when nsub=(0,0,...,0));
306 3 - the nieghbor is already in the tree, meaning that
307 is a pixel older than xi.
308 */
309 if(good &&
310 nindex != index &&
311 forest_pt[nindex].parent != node_is_void ) {
312
313 idx_t nrindex = 0, nvisited ;
314 val_t nrvalue = 0 ;
315
316#ifdef USE_RANK_UNION
317 int height = forest_pt [ rindex] .height ;
318 int nheight = forest_pt [nrindex] .height ;
319#endif
320
321 /* RINDEX = ROOT(INDEX) might change as we merge trees, so we
322 need to update it after each merge */
323
324 /* find the root of the current node */
325 /* also update the shortcuts */
326 nvisited = 0 ;
327 while( forest_pt[rindex].shortcut != rindex ) {
328 visited_pt[ nvisited++ ] = rindex ;
329 rindex = forest_pt[rindex].shortcut ;
330 }
331 while( nvisited-- ) {
332 forest_pt [ visited_pt[nvisited] ] .shortcut = rindex ;
333 }
334
335 /* find the root of the neighbor */
336 nrindex = nindex ;
337 nvisited = 0 ;
338 while( forest_pt[nrindex].shortcut != nrindex ) {
339 visited_pt[ nvisited++ ] = nrindex ;
340 nrindex = forest_pt[nrindex].shortcut ;
341 }
342 while( nvisited-- ) {
343 forest_pt [ visited_pt[nvisited] ] .shortcut = nrindex ;
344 }
345
346 /*
347 Now we join the two subtrees rooted at
348
349 RINDEX = ROOT(INDEX) and NRINDEX = ROOT(NINDEX).
350
351 Only three things can happen:
352
353 a - ROOT(INDEX) == ROOT(NRINDEX). In this case the two trees
354 have already been joined and we do not do anything.
355
356 b - I(ROOT(INDEX)) == I(ROOT(NRINDEX)). In this case index
357 is extending an extremal region with the same
358 value. Since ROOT(NRINDEX) will NOT be an extremal
359 region of the full image, ROOT(INDEX) can be safely
360 addedd as children of ROOT(NRINDEX) if this reduces
361 the height according to union rank.
362
363 c - I(ROOT(INDEX)) > I(ROOT(NRINDEX)) as index is extending
364 an extremal region, but increasing its level. In this
365 case ROOT(NRINDEX) WILL be an extremal region of the
366 final image and the only possibility is to add
367 ROOT(NRINDEX) as children of ROOT(INDEX).
368 */
369
370 if( rindex != nrindex ) {
371 /* this is a genuine join */
372
373 nrvalue = I_pt [nrindex] ;
374 if( nrvalue == value
375#ifdef USE_RANK_UNION
376 && height < nheight
377#endif
378 ) {
379 /* ROOT(INDEX) becomes the child */
380 forest_pt[rindex] .parent = nrindex ;
381 forest_pt[rindex] .shortcut = nrindex ;
382 forest_pt[nrindex].area += forest_pt[rindex].area ;
383
384#ifdef USE_RANK_UNION
385 forest_pt[nrindex].height = MAX(nheight, height+1) ;
386#endif
387
388 joins_pt[njoins++] = rindex ;
389
390 } else {
391 /* ROOT(index) becomes parent */
392 forest_pt[nrindex] .parent = rindex ;
393 forest_pt[nrindex] .shortcut = rindex ;
394 forest_pt[rindex] .area += forest_pt[nrindex].area ;
395
396#ifdef USE_RANK_UNION
397 forest_pt[rindex].height = MAX(height, nheight+1) ;
398#endif
399 if( nrvalue != value ) {
400 /* nrindex is extremal region: save for later */
401 forest_pt[nrindex].region = ner ;
402 regions_pt [ner] .index = nrindex ;
403 regions_pt [ner] .parent = ner ;
404 regions_pt [ner] .value = nrvalue ;
405 regions_pt [ner] .area = forest_pt [nrindex].area ;
406 regions_pt [ner] .area_top = nel ;
407 regions_pt [ner] .area_bot = 0 ;
408 ++ner ;
409 }
410
411 /* annote join operation for post-processing */
412 joins_pt[njoins++] = nrindex ;
413 }
414 }
415
416 } /* neighbor done */
417
418 /* move to next neighbor */
419 k = 0 ;
420 while(++ nsubs_pt [k] > 1) {
421 nsubs_pt [k++] = -1 ;
422 if(k == ndims) goto done_all_neighbors ;
423 }
424 } /* next neighbor */
425 done_all_neighbors : ;
426 } /* next pixel */
427
428
429 /* the root of the last processed pixel must be a region */
430 forest_pt [rindex].region = ner ;
431 regions_pt [ner] .index = rindex ;
432 regions_pt [ner] .parent = ner ;
433 regions_pt [ner] .value = I_pt [rindex] ;
434 regions_pt [ner] .area = forest_pt [rindex] .area ;
435 regions_pt [ner] .area_top = nel ;
436 regions_pt [ner] .area_bot = 0 ;
437 ++ner ;
438
439 verbose && mexPrintf("done\nExtremal regions: %d\n", ner) ;
440
441 /* -----------------------------------------------------------------
442 * Compute region parents
443 * -------------------------------------------------------------- */
444 for( i = 0 ; i < ner ; ++i) {
445 idx_t index = regions_pt [i].index ;
446 val_t value = regions_pt [i].value ;
447 idx_t j = i ;
448
449 while(j == i) {
450 idx_t pindex = forest_pt [index].parent ;
451 val_t pvalue = I_pt [pindex] ;
452
453 /* top of the tree */
454 if(index == pindex) {
455 j = forest_pt[index].region ;
456 break ;
457 }
458
459 /* if index is the root of a region, either this is still
460 i, or it is the parent region we are looking for. */
461 if(value < pvalue) {
462 j = forest_pt[index].region ;
463 }
464
465 index = pindex ;
466 value = pvalue ;
467 }
468 regions_pt[i]. parent = j ;
469 }
470
471 /* -----------------------------------------------------------------
472 * Compute areas of tops and bottoms
473 * -------------------------------------------------------------- */
474
475 /* We scan the list of regions from the bottom. Let x0 be the current
476 region and be x1 = PARENT(x0), x2 = PARENT(x1) and so on.
477
478 Here we do two things:
479
480 1) Look for regions x for which x0 is the BOTTOM. This requires
481 VAL(x0) <= VAL(x) - DELTA < VAL(x1).
482 We update AREA_BOT(x) for each of such x found.
483
484 2) Look for the region y which is the TOP of x0. This requires
485 VAL(y) <= VAL(x0) + DELTA < VAL(y+1)
486 We update AREA_TOP(x0) as soon as we find such y.
487
488 */
489
490 for( i = 0 ; i < ner ; ++i) {
491 /* fix xi as the region, then xj are the parents */
492 idx_t parent = regions_pt [i].parent ;
493 int val0 = regions_pt [i].value ;
494 int val1 = regions_pt [parent].value ;
495 int val = val0 ;
496 idx_t j = i ;
497
498 while(true) {
499 int valp = regions_pt [parent].value ;
500
501 /* i is the bottom of j */
502 if(val0 <= val - delta && val - delta < val1) {
503 regions_pt [j].area_bot =
504 MAX(regions_pt [j].area_bot, regions_pt [i].area) ;
505 }
506
507 /* j is the top of i */
508 if(val <= val0 + delta && val0 + delta < valp) {
509 regions_pt [i].area_top = regions_pt [j].area ;
510 }
511
512 /* stop if going on is useless */
513 if(val1 <= val - delta && val0 + delta < val)
514 break ;
515
516 /* stop also if j is the root */
517 if(j == parent)
518 break ;
519
520 /* next region upward */
521 j = parent ;
522 parent = regions_pt [j].parent ;
523 val = valp ;
524 }
525 }
526
527 /* -----------------------------------------------------------------
528 * Compute variation
529 * -------------------------------------------------------------- */
530 for(i = 0 ; i < ner ; ++i) {
531 int area = regions_pt [i].area ;
532 int area_top = regions_pt [i].area_top ;
533 int area_bot = regions_pt [i].area_bot ;
534 regions_pt [i].variation =
535 (float)(area_top - area_bot) / (float)area ;
536
537 /* initialize .mastable to 1 for all nodes */
538 regions_pt [i].maxstable = 1 ;
539 }
540
541 /* -----------------------------------------------------------------
542 * Remove regions which are NOT maximally stable
543 * -------------------------------------------------------------- */
544 nmer = ner ;
545 for(i = 0 ; i < ner ; ++i) {
546 idx_t parent = regions_pt [i] .parent ;
547 float var = regions_pt [i] .variation ;
548 float pvar = regions_pt [parent] .variation ;
549 idx_t loser ;
550
551 /* decide which one to keep and put that in loser */
552 if(var < pvar) loser = parent ; else loser = i ;
553
554 /* make loser NON maximally stable */
555 if(regions_pt [loser].maxstable) --nmer ;
556 regions_pt [loser].maxstable = 0 ;
557 }
558
559 verbose && mexPrintf("Maximally stable regions: %d (%.1f%%)\n",
560 nmer, 100.0 * (double) nmer / ner) ;
561
562 /* -----------------------------------------------------------------
563 * Remove more regions
564 * -------------------------------------------------------------- */
565
566 /* it is critical for correct duplicate detection to remove regions
567 from the bottom (smallest one first) */
568
569 if( big_cleanup || small_cleanup || bad_cleanup || dup_cleanup ) {
570 int nbig = 0 ;
571 int nsmall = 0 ;
572 int nbad = 0 ;
573 int ndup = 0 ;
574
575 /* scann all extremal regions */
576 for(i = 0 ; i < ner ; ++i) {
577
578 /* process only maximally stable extremal regions */
579 if(! regions_pt [i].maxstable) continue ;
580
581 if( bad_cleanup && regions_pt[i].variation >= 1.0f ) {
582 ++nbad ;
583 goto remove_this_region ;
584 }
585
586 if( big_cleanup && regions_pt[i].area > nel/2 ) {
587 ++nbig ;
588 goto remove_this_region ;
589 }
590
591 if( small_cleanup && regions_pt[i].area < 25 ) {
592 ++nsmall ;
593 goto remove_this_region ;
594 }
595
596 /*
597 * Remove duplicates
598 */
599 if( dup_cleanup ) {
600 idx_t parent = regions_pt [i].parent ;
601 int area, parea ;
602 float change ;
603
604 /* the search does not apply to root regions */
605 if(parent != i) {
606
607 /* search for the maximally stable parent region */
608 while(! regions_pt[parent].maxstable) {
609 idx_t next = regions_pt[parent].parent ;
610 if(next == parent) break ;
611 parent = next ;
612 }
613
614 /* compare with the parent region; if the current and parent
615 regions are too similar, keep only the parent */
616 area = regions_pt [i].area ;
617 parea = regions_pt [parent].area ;
618 change = (float)(parea - area)/area ;
619
620 if(change < 0.5) {
621 ++ndup ;
622 goto remove_this_region ;
623 }
624
625 } /* drop duplicates */
626 }
627 continue ;
628 remove_this_region :
629 regions_pt[i].maxstable = false ;
630 --nmer ;
631 } /* next region to cleanup */
632
633 if(verbose) {
634 mexPrintf(" Bad regions: %d\n", nbad ) ;
635 mexPrintf(" Small regions: %d\n", nsmall ) ;
636 mexPrintf(" Big regions: %d\n", nbig ) ;
637 mexPrintf(" Duplicated regions: %d\n", ndup ) ;
638 }
639 }
640
641 verbose && mexPrintf("Cleaned-up regions: %d (%.1f%%)\n",
642 nmer, 100.0 * (double) nmer / ner) ;
643
644 /* -----------------------------------------------------------------
645 * Fit ellipses
646 * -------------------------------------------------------------- */
647
648 ell_pt = 0 ;
649 if (nout >= 1) {
650 int midx = 1 ;
651 int d, index, j ;
652
653 verbose && mexPrintf("Fitting ellipses...\n") ;
654
655 /* enumerate maxstable regions */
656 for(i = 0 ; i < ner ; ++i) {
657 if(! regions_pt [i].maxstable) continue ;
658 regions_pt [i].maxstable = midx++ ;
659 }
660
661 /* allocate space */
662 acc_pt = mxMalloc(sizeof(acc_t) * nel) ;
663 ell_pt = mxMalloc(sizeof(acc_t) * gdl * nmer) ;
664
665 /* clear accumulators */
666 memset(ell_pt, 0, sizeof(int) * gdl * nmer) ;
667
668 /* for each gdl */
669 for(d = 0 ; d < gdl ; ++d) {
670 /* initalize parameter */
671 memset(subs_pt, 0, sizeof(int) * ndims) ;
672
673 if(d < ndims) {
674 verbose && mexPrintf(" mean %d\n",d) ;
675 for(index = 0 ; index < nel ; ++ index) {
676 acc_pt[index] = subs_pt[d] ;
677 adv(dims, ndims, subs_pt) ;
678 }
679
680 } else {
681
682 /* decode d-ndims into a (i,j) pair */
683 i = d-ndims ;
684 j = 0 ;
685 while(i > j) {
686 i -= j + 1 ;
687 j ++ ;
688 }
689
690 verbose && mexPrintf(" corr (%d,%d)\n",i,j) ;
691
692 /* add x_i * x_j */
693 for(index = 0 ; index < nel ; ++ index){
694 acc_pt[index] = subs_pt[i]*subs_pt[j] ;
695 adv(dims, ndims, subs_pt) ;
696 }
697 }
698
699 /* integrate parameter */
700 for(i = 0 ; i < njoins ; ++i) {
701 idx_t index = joins_pt[i] ;
702 idx_t parent = forest_pt [ index ].parent ;
703 acc_pt[parent] += acc_pt[index] ;
704 }
705
706 /* save back to ellpises */
707 for(i = 0 ; i < ner ; ++i) {
708 idx_t region = regions_pt [i].maxstable ;
709
710 /* skip if not extremal region */
711 if(region-- == 0) continue ;
712 ell_pt [d + gdl*region] = acc_pt [ regions_pt[i].index ] ;
713 }
714
715 /* next gdl */
716 }
717 mxFree(acc_pt) ;
718 }
719
720
721 /* -----------------------------------------------------------------
722 * Save back and exit
723 * -------------------------------------------------------------- */
724
725 /*
726 * Save extremal regions
727 */
728 {
729 int dims[2] ;
730 int unsigned * pt ;
731 dims[0] = nmer ;
732 out[OUT_REGIONS] = mxCreateNumericArray(1,dims,mxUINT32_CLASS,mxREAL);
733 pt = mxGetData(out[OUT_REGIONS]) ;
734 for (i = 0 ; i < ner ; ++i) {
735 if( regions_pt[i].maxstable ) {
736 /* adjust for MATLAB index compatibility */
737 *pt++ = regions_pt[i].index + 1 ;
738 }
739 }
740 }
741
742 /*
743 * Save fitted ellipses
744 */
745 if(nout >= 2) {
746 int dims[2], d, j, index ;
747 double * pt ;
748 dims[0] = gdl ;
749 dims[1] = nmer ;
750
751 out[OUT_ELL] = mxCreateNumericArray(2,dims,mxDOUBLE_CLASS,mxREAL) ;
752 pt = mxGetData(out[OUT_ELL]) ;
753
754 for(index = 0 ; index < nel ; ++index) {
755
756 idx_t region = regions_pt [index] .maxstable ;
757 int N = regions_pt [index] .area ;
758
759 if(region-- == 0) continue ;
760
761 for(d = 0 ; d < gdl ; ++d) {
762
763 pt[d] = (double) ell_pt[gdl*region + d] / N ;
764
765 if(d < ndims) {
766 /* adjust for MATLAB coordinate frame convention */
767 pt[d] += 1 ;
768 } else {
769 /* remove squared mean from moment to get variance */
770 i = d - ndims ;
771 j = 0 ;
772 while(i > j) {
773 i -= j + 1 ;
774 j ++ ;
775 }
776 pt[d] -= (pt[i]-1)*(pt[j]-1) ;
777 }
778 }
779 pt += gdl ;
780 }
781 mxFree(ell_pt) ;
782 }
783
784 if(nout >= 3) {
785 int unsigned * pt ;
786 out[OUT_PARENTS] = mxCreateNumericArray(ndims,dims,mxUINT32_CLASS,mxREAL) ;
787 pt = mxGetData(out[OUT_PARENTS]) ;
788 for(i = 0 ; i < nel ; ++i) {
789 *pt++ = forest_pt[i].parent ;
790 }
791 }
792
793 if(nout >= 4) {
794 int dims[2] ;
795 int unsigned * pt ;
796 dims[0] = 3 ;
797 dims[1]= ner ;
798 out[OUT_AREA] = mxCreateNumericArray(2,dims,mxUINT32_CLASS,mxREAL);
799 pt = mxGetData(out[OUT_AREA]) ;
800 for( i = 0 ; i < ner ; ++i ) {
801 *pt++ = regions_pt [i]. area_bot ;
802 *pt++ = regions_pt [i]. area ;
803 *pt++ = regions_pt [i]. area_top ;
804 }
805 }
806
807 /* free stuff */
808 mxFree( forest_pt ) ;
809 mxFree( pairs_pt ) ;
810 mxFree( regions_pt ) ;
811 mxFree( visited_pt ) ;
812 mxFree( strides_pt ) ;
813 mxFree( nsubs_pt ) ;
814 mxFree( subs_pt ) ;
815}
generated by cgit v1.2.2 (git 2.25.0) at 2025-04-30 19:38:01 -0400