From 0efc775370d2ff91927d1b383a99eab78dc5538f Mon Sep 17 00:00:00 2001 From: leochanj Date: Wed, 21 Oct 2020 01:52:54 -0400 Subject: debug libextra and remove matlab FLUSH_CACHES --- SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c | 809 ----------------------- 1 file changed, 809 deletions(-) delete mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c (limited to 'SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c') 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 @@ -/* file: mser.mex.c -** description: Maximally Stable Extremal Regions -** author: Andrea Vedaldi -**/ - -/* AUTORIGHTS -Copyright (C) 2006 Regents of the University of California -All rights reserved - -Written by Andrea Vedaldi (UCLA VisionLab). - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met - - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. - * Neither the name of the University of California, Berkeley nor the - names of its contributors may be used to endorse or promote products - derived from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY -EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY -DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -/** @file - ** @brief Maximally Stable Extremal Regions - MEX implementation - **/ - -#include -#include -#include -#include -#include -#include - -#define MIN(x,y) (((x)<(y))?(x):(y)) -#define MAX(x,y) (((x)>(y))?(x):(y)) - -#define USE_BUCKET_SORT -/*#define USE_RANK_UNION -*/ - -typedef char unsigned val_t ; -typedef int unsigned idx_t ; -typedef long long int unsigned acc_t ; - -/* pairs are used to sort the pixels */ -typedef struct -{ - val_t value ; - idx_t index ; -} pair_t ; - -/* forest node */ -typedef struct -{ - idx_t parent ; /**< parent pixel */ - idx_t shortcut ; /**< shortcut to the root */ - idx_t region ; /**< index of the region */ - int area ; /**< area of the region */ -#ifdef USE_RANK_UNION - int height ; /**< node height */ -#endif -} node_t ; - -/* extremal regions */ -typedef struct -{ - idx_t parent ; /**< parent region */ - idx_t index ; /**< index of root pixel */ - val_t value ; /**< value of root pixel */ - int area ; /**< area of the region */ - int area_top ; /**< area of the region DELTA levels above */ - int area_bot ; /**< area of the region DELTA levels below */ - float variation ; /**< variation */ - int maxstable ; /**< max stable number (=0 if not maxstable) */ -} region_t ; - -/* predicate used to sort pixels by increasing intensity */ -int -cmp_pair(void const* a, void const* b) -{ - pair_t* pa = (pair_t*) a; - pair_t* pb = (pair_t*) b; - return pa->value - pb->value ; -} - -/* advance N-dimensional subscript */ -void -adv(int const* dims, int ndims, int* subs_pt) -{ - int d = 0 ; - while(d < ndims) { - if( ++subs_pt[d] < dims[d] ) return ; - subs_pt[d++] = 0 ; - } -} - -/* driver */ -void -mexFunction(int nout, mxArray *out[], - int nin, const mxArray *in[]) -{ - enum {IN_I=0, IN_DELTA} ; - enum {OUT_REGIONS=0, OUT_ELL, OUT_PARENTS, OUT_AREA} ; - - idx_t i ; - idx_t rindex = 0 ; - int k ; - - /* configuration */ - int verbose = 1 ; /* be verbose */ - int small_cleanup= 1 ; /* remove very small regions */ - int big_cleanup = 1 ; /* remove very big regions */ - int bad_cleanup = 0 ; /* remove very bad regions */ - int dup_cleanup = 1 ; /* remove duplicates */ - val_t delta ; /* stability delta */ - - /* node value denoting a void node */ - idx_t const node_is_void = 0xffffffff ; - - int* subs_pt ; /* N-dimensional subscript */ - int* nsubs_pt ; /* diff-subscript to point to neigh. */ - idx_t* strides_pt ; /* strides to move in image array */ - idx_t* visited_pt ; /* flag */ - - int nel ; /* number of image elements (pixels) */ - int ner = 0 ; /* number of extremal regions */ - int nmer = 0 ; /* number of maximally stable */ - int ndims ; /* number of dimensions */ - int const* dims ; /* dimensions */ - int njoins = 0 ; /* number of join ops */ - - val_t const* I_pt ; /* source image */ - pair_t* pairs_pt ; /* scratch buffer to sort pixels */ - node_t* forest_pt ; /* the extremal regions forest */ - region_t* regions_pt ; /* list of extremal regions found */ - - /* ellipses fitting */ - acc_t* acc_pt ; /* accumulator to integrate region moments */ - acc_t* ell_pt ; /* ellipses parameters */ - int gdl ; /* number of parameters of an ellipse */ - idx_t* joins_pt ; /* sequence of joins */ - - /** ----------------------------------------------------------------- - ** Check the arguments - ** -------------------------------------------------------------- */ - if (nin != 2) { - mexErrMsgTxt("Two arguments required.") ; - } else if (nout > 4) { - mexErrMsgTxt("Too many output arguments."); - } - - if(mxGetClassID(in[IN_I]) != mxUINT8_CLASS) { - mexErrMsgTxt("I must be of class UINT8") ; - } - - if(!uIsScalar(in[IN_DELTA])) { - mexErrMsgTxt("DELTA must be scalar") ; - } - - delta = 0 ; - switch(mxGetClassID(in[IN_DELTA])) { - case mxUINT8_CLASS : - delta = * (val_t*) mxGetData(in[IN_DELTA]) ; - break ; - - case mxDOUBLE_CLASS : - { - double x = *mxGetPr(in[IN_DELTA]) ; - if(x < 0.0) { - mexErrMsgTxt("DELTA must be non-negative") ; - } - delta = (val_t) x ; - } - break ; - - default : - mexErrMsgTxt("DELTA must be of class DOUBLE or UINT8") ; - } - - /* get dimensions */ - nel = mxGetNumberOfElements(in[IN_I]) ; - ndims = mxGetNumberOfDimensions(in[IN_I]) ; - dims = mxGetDimensions(in[IN_I]) ; - I_pt = mxGetData(in[IN_I]) ; - - /* allocate stuff */ - subs_pt = mxMalloc( sizeof(int) * ndims ) ; - nsubs_pt = mxMalloc( sizeof(int) * ndims ) ; - strides_pt = mxMalloc( sizeof(idx_t) * ndims ) ; - visited_pt = mxMalloc( sizeof(idx_t) * nel ) ; - regions_pt = mxMalloc( sizeof(region_t) * nel ) ; - pairs_pt = mxMalloc( sizeof(pair_t) * nel ) ; - forest_pt = mxMalloc( sizeof(node_t) * nel ) ; - joins_pt = mxMalloc( sizeof(idx_t) * nel ) ; - - /* compute strides to move into the N-dimensional image array */ - strides_pt [0] = 1 ; - for(k = 1 ; k < ndims ; ++k) { - strides_pt [k] = strides_pt [k-1] * dims [k-1] ; - } - - /* sort pixels by increasing intensity*/ - verbose && mexPrintf("Sorting pixels ... ") ; - -#ifndef USE_BUCKETSORT - for(i = 0 ; i < nel ; ++i) { - pairs_pt [i].value = I_pt [i] ; - pairs_pt [i].index = i ; - } - qsort(pairs_pt, nel, sizeof(pair_t), cmp_pair) ; -#else - { - int unsigned buckets [256] ; - memset(buckets, 0, sizeof(int unsigned)*256) ; - for(i = 0 ; i < nel ; ++i) { - val_t v = I_pt [i] ; - ++ buckets[v] ; - } - for(i = 1 ; i < 256 ; ++i) { - buckets[i] += buckets[i-1] ; - } - for(i = nel ; i >= 1 ; ) { - val_t v = I_pt [--i] ; - idx_t j = -- buckets [v] ; - pairs_pt [j].value = v ; - pairs_pt [j].index = i ; - } - } -#endif - verbose && mexPrintf("done\n") ; - - /* initialize the forest with all void nodes */ - for(i = 0 ; i < nel ; ++i) { - forest_pt [i].parent = node_is_void ; - } - - /* number of ellipse free parameters */ - gdl = ndims*(ndims+1)/2 + ndims ; - - /* ----------------------------------------------------------------- - * Compute extremal regions tree - * -------------------------------------------------------------- */ - verbose && mexPrintf("Computing extremal regions ... ") ; - for(i = 0 ; i < nel ; ++i) { - - /* pop next node xi */ - idx_t index = pairs_pt [i].index ; - val_t value = pairs_pt [i].value ; - - /* this will be needed later */ - rindex = index ; - - /* push it into the tree */ - forest_pt [index] .parent = index ; - forest_pt [index] .shortcut = index ; - forest_pt [index] .area = 1 ; -#ifdef USE_RANK_UNION - forest_pt [index] .height = 1 ; -#endif - - /* convert index into a subscript sub; also initialize nsubs - to (-1,-1,...,-1) */ - { - idx_t temp = index ; - for(k = ndims-1 ; k >=0 ; --k) { - nsubs_pt [k] = -1 ; - subs_pt [k] = temp / strides_pt [k] ; - temp = temp % strides_pt [k] ; - } - } - - /* process neighbors of xi */ - while( true ) { - int good = true ; - idx_t nindex = 0 ; - - /* compute NSUBS+SUB, the correspoinding neighbor index NINDEX - and check that the pixel is within image boundaries. */ - for(k = 0 ; k < ndims && good ; ++k) { - int temp = nsubs_pt [k] + subs_pt [k] ; - good &= 0 <= temp && temp < dims[k] ; - nindex += temp * strides_pt [k] ; - } - - /* keep going only if - 1 - the neighbor is within image boundaries; - 2 - the neighbor is indeed different from the current node - (this happens when nsub=(0,0,...,0)); - 3 - the nieghbor is already in the tree, meaning that - is a pixel older than xi. - */ - if(good && - nindex != index && - forest_pt[nindex].parent != node_is_void ) { - - idx_t nrindex = 0, nvisited ; - val_t nrvalue = 0 ; - -#ifdef USE_RANK_UNION - int height = forest_pt [ rindex] .height ; - int nheight = forest_pt [nrindex] .height ; -#endif - - /* RINDEX = ROOT(INDEX) might change as we merge trees, so we - need to update it after each merge */ - - /* find the root of the current node */ - /* also update the shortcuts */ - nvisited = 0 ; - while( forest_pt[rindex].shortcut != rindex ) { - visited_pt[ nvisited++ ] = rindex ; - rindex = forest_pt[rindex].shortcut ; - } - while( nvisited-- ) { - forest_pt [ visited_pt[nvisited] ] .shortcut = rindex ; - } - - /* find the root of the neighbor */ - nrindex = nindex ; - nvisited = 0 ; - while( forest_pt[nrindex].shortcut != nrindex ) { - visited_pt[ nvisited++ ] = nrindex ; - nrindex = forest_pt[nrindex].shortcut ; - } - while( nvisited-- ) { - forest_pt [ visited_pt[nvisited] ] .shortcut = nrindex ; - } - - /* - Now we join the two subtrees rooted at - - RINDEX = ROOT(INDEX) and NRINDEX = ROOT(NINDEX). - - Only three things can happen: - - a - ROOT(INDEX) == ROOT(NRINDEX). In this case the two trees - have already been joined and we do not do anything. - - b - I(ROOT(INDEX)) == I(ROOT(NRINDEX)). In this case index - is extending an extremal region with the same - value. Since ROOT(NRINDEX) will NOT be an extremal - region of the full image, ROOT(INDEX) can be safely - addedd as children of ROOT(NRINDEX) if this reduces - the height according to union rank. - - c - I(ROOT(INDEX)) > I(ROOT(NRINDEX)) as index is extending - an extremal region, but increasing its level. In this - case ROOT(NRINDEX) WILL be an extremal region of the - final image and the only possibility is to add - ROOT(NRINDEX) as children of ROOT(INDEX). - */ - - if( rindex != nrindex ) { - /* this is a genuine join */ - - nrvalue = I_pt [nrindex] ; - if( nrvalue == value -#ifdef USE_RANK_UNION - && height < nheight -#endif - ) { - /* ROOT(INDEX) becomes the child */ - forest_pt[rindex] .parent = nrindex ; - forest_pt[rindex] .shortcut = nrindex ; - forest_pt[nrindex].area += forest_pt[rindex].area ; - -#ifdef USE_RANK_UNION - forest_pt[nrindex].height = MAX(nheight, height+1) ; -#endif - - joins_pt[njoins++] = rindex ; - - } else { - /* ROOT(index) becomes parent */ - forest_pt[nrindex] .parent = rindex ; - forest_pt[nrindex] .shortcut = rindex ; - forest_pt[rindex] .area += forest_pt[nrindex].area ; - -#ifdef USE_RANK_UNION - forest_pt[rindex].height = MAX(height, nheight+1) ; -#endif - if( nrvalue != value ) { - /* nrindex is extremal region: save for later */ - forest_pt[nrindex].region = ner ; - regions_pt [ner] .index = nrindex ; - regions_pt [ner] .parent = ner ; - regions_pt [ner] .value = nrvalue ; - regions_pt [ner] .area = forest_pt [nrindex].area ; - regions_pt [ner] .area_top = nel ; - regions_pt [ner] .area_bot = 0 ; - ++ner ; -/* printf("ner = %d\n", ner);*/ - } - - /* annote join operation for post-processing */ - joins_pt[njoins++] = nrindex ; - } - } - - } /* neighbor done */ - - /* move to next neighbor */ - k = 0 ; - while(++ nsubs_pt [k] > 1) { - nsubs_pt [k++] = -1 ; - if(k == ndims) goto done_all_neighbors ; - } - } /* next neighbor */ - done_all_neighbors : ; - } /* next pixel */ - - /* the root of the last processed pixel must be a region */ - forest_pt [rindex].region = ner ; - regions_pt [ner] .index = rindex ; - regions_pt [ner] .parent = ner ; - regions_pt [ner] .value = I_pt [rindex] ; - regions_pt [ner] .area = forest_pt [rindex] .area ; - regions_pt [ner] .area_top = nel ; - regions_pt [ner] .area_bot = 0 ; - ++ner ; - - verbose && mexPrintf("done\nExtremal regions: %d\n", ner) ; - - /* ----------------------------------------------------------------- - * Compute region parents - * -------------------------------------------------------------- */ - for( i = 0 ; i < ner ; ++i) { - idx_t index = regions_pt [i].index ; - val_t value = regions_pt [i].value ; - idx_t j = i ; - - while(j == i) { - idx_t pindex = forest_pt [index].parent ; - val_t pvalue = I_pt [pindex] ; - - /* top of the tree */ - if(index == pindex) { - j = forest_pt[index].region ; - break ; - } - - /* if index is the root of a region, either this is still - i, or it is the parent region we are looking for. */ - if(value < pvalue) { - j = forest_pt[index].region ; - } - - index = pindex ; - value = pvalue ; - } - regions_pt[i]. parent = j ; - } - - /* ----------------------------------------------------------------- - * Compute areas of tops and bottoms - * -------------------------------------------------------------- */ - - /* We scan the list of regions from the bottom. Let x0 be the current - region and be x1 = PARENT(x0), x2 = PARENT(x1) and so on. - - Here we do two things: - - 1) Look for regions x for which x0 is the BOTTOM. This requires - VAL(x0) <= VAL(x) - DELTA < VAL(x1). - We update AREA_BOT(x) for each of such x found. - - 2) Look for the region y which is the TOP of x0. This requires - VAL(y) <= VAL(x0) + DELTA < VAL(y+1) - We update AREA_TOP(x0) as soon as we find such y. - - */ - - for( i = 0 ; i < ner ; ++i) { - /* fix xi as the region, then xj are the parents */ - idx_t parent = regions_pt [i].parent ; - int val0 = regions_pt [i].value ; - int val1 = regions_pt [parent].value ; - int val = val0 ; - idx_t j = i ; - - while(true) { - int valp = regions_pt [parent].value ; - - /* i is the bottom of j */ - if(val0 <= val - delta && val - delta < val1) { - regions_pt [j].area_bot = - MAX(regions_pt [j].area_bot, regions_pt [i].area) ; - } - - /* j is the top of i */ - if(val <= val0 + delta && val0 + delta < valp) { - regions_pt [i].area_top = regions_pt [j].area ; - } - - /* stop if going on is useless */ - if(val1 <= val - delta && val0 + delta < val) - break ; - - /* stop also if j is the root */ - if(j == parent) - break ; - - /* next region upward */ - j = parent ; - parent = regions_pt [j].parent ; - val = valp ; - } - } - - /* ----------------------------------------------------------------- - * Compute variation - * -------------------------------------------------------------- */ - for(i = 0 ; i < ner ; ++i) { - int area = regions_pt [i].area ; - int area_top = regions_pt [i].area_top ; - int area_bot = regions_pt [i].area_bot ; - regions_pt [i].variation = - (float)(area_top - area_bot) / (float)area ; - - /* initialize .mastable to 1 for all nodes */ - regions_pt [i].maxstable = 1 ; - } - - /* ----------------------------------------------------------------- - * Remove regions which are NOT maximally stable - * -------------------------------------------------------------- */ - nmer = ner ; - for(i = 0 ; i < ner ; ++i) { - idx_t parent = regions_pt [i] .parent ; - float var = regions_pt [i] .variation ; - float pvar = regions_pt [parent] .variation ; - idx_t loser ; - - /* decide which one to keep and put that in loser */ - if(var < pvar) loser = parent ; else loser = i ; - - /* make loser NON maximally stable */ - if(regions_pt [loser].maxstable) --nmer ; - regions_pt [loser].maxstable = 0 ; - } - - verbose && mexPrintf("Maximally stable regions: %d (%.1f%%)\n", - nmer, 100.0 * (double) nmer / ner) ; - - /* ----------------------------------------------------------------- - * Remove more regions - * -------------------------------------------------------------- */ - - /* it is critical for correct duplicate detection to remove regions - from the bottom (smallest one first) */ - - if( big_cleanup || small_cleanup || bad_cleanup || dup_cleanup ) { - int nbig = 0 ; - int nsmall = 0 ; - int nbad = 0 ; - int ndup = 0 ; - - /* scann all extremal regions */ - for(i = 0 ; i < ner ; ++i) { - - /* process only maximally stable extremal regions */ - if(! regions_pt [i].maxstable) continue ; - - if( bad_cleanup && regions_pt[i].variation >= 1.0f ) { - ++nbad ; - goto remove_this_region ; - } - - if( big_cleanup && regions_pt[i].area > nel/2 ) { - ++nbig ; - goto remove_this_region ; - } - - if( small_cleanup && regions_pt[i].area < 25 ) { - ++nsmall ; - goto remove_this_region ; - } - - /* - * Remove duplicates - */ - if( dup_cleanup ) { - idx_t parent = regions_pt [i].parent ; - int area, parea ; - float change ; - - /* the search does not apply to root regions */ - if(parent != i) { - - /* search for the maximally stable parent region */ - while(! regions_pt[parent].maxstable) { - idx_t next = regions_pt[parent].parent ; - if(next == parent) break ; - parent = next ; - } - - /* compare with the parent region; if the current and parent - regions are too similar, keep only the parent */ - area = regions_pt [i].area ; - parea = regions_pt [parent].area ; - change = (float)(parea - area)/area ; - - if(change < 0.5) { - ++ndup ; - goto remove_this_region ; - } - - } /* drop duplicates */ - } - continue ; - remove_this_region : - regions_pt[i].maxstable = false ; - --nmer ; - } /* next region to cleanup */ - - if(verbose) { - mexPrintf(" Bad regions: %d\n", nbad ) ; - mexPrintf(" Small regions: %d\n", nsmall ) ; - mexPrintf(" Big regions: %d\n", nbig ) ; - mexPrintf(" Duplicated regions: %d\n", ndup ) ; - } - } - - verbose && mexPrintf("Cleaned-up regions: %d (%.1f%%)\n", - nmer, 100.0 * (double) nmer / ner) ; - - /* ----------------------------------------------------------------- - * Fit ellipses - * -------------------------------------------------------------- */ - ell_pt = 0 ; - if (nout >= 1) { - int midx = 1 ; - int d, index, j ; - - verbose && mexPrintf("Fitting ellipses...\n") ; - - /* enumerate maxstable regions */ - for(i = 0 ; i < ner ; ++i) { - if(! regions_pt [i].maxstable) continue ; - regions_pt [i].maxstable = midx++ ; - } - - /* allocate space */ - acc_pt = mxMalloc(sizeof(acc_t) * nel) ; - ell_pt = mxMalloc(sizeof(acc_t) * gdl * nmer) ; - - /* clear accumulators */ - memset(ell_pt, 0, sizeof(int) * gdl * nmer) ; - - /* for each gdl */ - for(d = 0 ; d < gdl ; ++d) { - /* initalize parameter */ - memset(subs_pt, 0, sizeof(int) * ndims) ; - - if(d < ndims) { - verbose && mexPrintf(" mean %d\n",d) ; - for(index = 0 ; index < nel ; ++ index) { - acc_pt[index] = subs_pt[d] ; - adv(dims, ndims, subs_pt) ; - } - - } else { - - /* decode d-ndims into a (i,j) pair */ - i = d-ndims ; - j = 0 ; - while(i > j) { - i -= j + 1 ; - j ++ ; - } - - verbose && mexPrintf(" corr (%d,%d)\n",i,j) ; - - /* add x_i * x_j */ - for(index = 0 ; index < nel ; ++ index){ - acc_pt[index] = subs_pt[i]*subs_pt[j] ; - adv(dims, ndims, subs_pt) ; - } - } - - /* integrate parameter */ - for(i = 0 ; i < njoins ; ++i) { - idx_t index = joins_pt[i] ; - idx_t parent = forest_pt [ index ].parent ; - acc_pt[parent] += acc_pt[index] ; - } - - /* save back to ellpises */ - for(i = 0 ; i < ner ; ++i) { - idx_t region = regions_pt [i].maxstable ; - - /* skip if not extremal region */ - if(region-- == 0) continue ; - ell_pt [d + gdl*region] = acc_pt [ regions_pt[i].index ] ; - } - - /* next gdl */ - } - mxFree(acc_pt) ; - } - - /* ----------------------------------------------------------------- - * Save back and exit - * -------------------------------------------------------------- */ - - /* - * Save extremal regions - */ - { - int dims[2] ; - int unsigned * pt ; - dims[0] = nmer ; - out[OUT_REGIONS] = mxCreateNumericArray(1,dims,mxUINT32_CLASS,mxREAL); - pt = mxGetData(out[OUT_REGIONS]) ; - for (i = 0 ; i < ner ; ++i) { - if( regions_pt[i].maxstable ) { - /* adjust for MATLAB index compatibility */ - *pt++ = regions_pt[i].index + 1 ; - } - } - } - - /* - * Save fitted ellipses - */ - if(nout >= 2) { - int dims[2], d, j, index ; - double * pt ; - dims[0] = gdl ; - dims[1] = nmer ; - - out[OUT_ELL] = mxCreateNumericArray(2,dims,mxDOUBLE_CLASS,mxREAL) ; - pt = mxGetData(out[OUT_ELL]) ; - - for(index = 0 ; index < nel ; ++index) { - - idx_t region = regions_pt [index] .maxstable ; - int N = regions_pt [index] .area ; - - if(region-- == 0) continue ; - - for(d = 0 ; d < gdl ; ++d) { - - pt[d] = (double) ell_pt[gdl*region + d] / N ; - - if(d < ndims) { - /* adjust for MATLAB coordinate frame convention */ - pt[d] += 1 ; - } else { - /* remove squared mean from moment to get variance */ - i = d - ndims ; - j = 0 ; - while(i > j) { - i -= j + 1 ; - j ++ ; - } - pt[d] -= (pt[i]-1)*(pt[j]-1) ; - } - } - pt += gdl ; - } - mxFree(ell_pt) ; - } - - if(nout >= 3) { - int unsigned * pt ; - out[OUT_PARENTS] = mxCreateNumericArray(ndims,dims,mxUINT32_CLASS,mxREAL) ; - pt = mxGetData(out[OUT_PARENTS]) ; - for(i = 0 ; i < nel ; ++i) { - *pt++ = forest_pt[i].parent ; - } - } - - if(nout >= 4) { - int dims[2] ; - int unsigned * pt ; - dims[0] = 3 ; - dims[1]= ner ; - out[OUT_AREA] = mxCreateNumericArray(2,dims,mxUINT32_CLASS,mxREAL); - pt = mxGetData(out[OUT_AREA]) ; - for( i = 0 ; i < ner ; ++i ) { - *pt++ = regions_pt [i]. area_bot ; - *pt++ = regions_pt [i]. area ; - *pt++ = regions_pt [i]. area_top ; - } - } - - /* free stuff */ - mxFree( forest_pt ) ; - mxFree( pairs_pt ) ; - mxFree( regions_pt ) ; - mxFree( visited_pt ) ; - mxFree( strides_pt ) ; - mxFree( nsubs_pt ) ; - mxFree( subs_pt ) ; -} -- cgit v1.2.2