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