From f618466c25d43f3bae9e40920273bf77de1e1149 Mon Sep 17 00:00:00 2001 From: leochanj105 Date: Mon, 19 Oct 2020 23:09:30 -0400 Subject: initial sd-vbs initial sd-vbs add sd-vbs sd-vbs --- SD-VBS/benchmarks/mser/src/c/mser.c | 714 ++++++++++++++++++ SD-VBS/benchmarks/mser/src/c/mser.h | 83 +++ SD-VBS/benchmarks/mser/src/c/script_mser.c | 120 +++ SD-VBS/benchmarks/mser/src/matlab/Makefile | 123 ++++ SD-VBS/benchmarks/mser/src/matlab/TIMESTAMP | 2 + SD-VBS/benchmarks/mser/src/matlab/erfill.m | 13 + SD-VBS/benchmarks/mser/src/matlab/erfill.mex.c | 223 ++++++ SD-VBS/benchmarks/mser/src/matlab/erfill.mexa64 | Bin 0 -> 10930 bytes SD-VBS/benchmarks/mser/src/matlab/erfill.mexglx | Bin 0 -> 8744 bytes SD-VBS/benchmarks/mser/src/matlab/mexutils.c | 111 +++ SD-VBS/benchmarks/mser/src/matlab/mser.mex.c | 815 +++++++++++++++++++++ SD-VBS/benchmarks/mser/src/matlab/mser.mexa64 | Bin 0 -> 14107 bytes SD-VBS/benchmarks/mser/src/matlab/mser.mexglx | Bin 0 -> 12459 bytes SD-VBS/benchmarks/mser/src/matlab/mser_compile.m | 7 + SD-VBS/benchmarks/mser/src/matlab/mser_demo2.m | 62 ++ SD-VBS/benchmarks/mser/src/matlab/mser_demo3.m | 117 +++ SD-VBS/benchmarks/mser/src/matlab/old/Makefile | 123 ++++ SD-VBS/benchmarks/mser/src/matlab/old/TIMESTAMP | 2 + SD-VBS/benchmarks/mser/src/matlab/old/erfill.m | 13 + SD-VBS/benchmarks/mser/src/matlab/old/erfill.mex.c | 223 ++++++ .../benchmarks/mser/src/matlab/old/erfill.mexa64 | Bin 0 -> 18714 bytes SD-VBS/benchmarks/mser/src/matlab/old/mexutils.c | 111 +++ SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c | 809 ++++++++++++++++++++ SD-VBS/benchmarks/mser/src/matlab/old/mser.mexa64 | Bin 0 -> 28582 bytes .../benchmarks/mser/src/matlab/old/mser_compile.m | 7 + SD-VBS/benchmarks/mser/src/matlab/old/mser_demo2.m | 62 ++ SD-VBS/benchmarks/mser/src/matlab/old/mser_demo3.m | 117 +++ .../benchmarks/mser/src/matlab/old/overview_mser.m | 8 + .../mser/src/matlab/old/script_run_profile.m | 137 ++++ SD-VBS/benchmarks/mser/src/matlab/overview_mser.m | 8 + .../mser/src/matlab/script_run_profile.m | 126 ++++ 31 files changed, 4136 insertions(+) create mode 100644 SD-VBS/benchmarks/mser/src/c/mser.c create mode 100644 SD-VBS/benchmarks/mser/src/c/mser.h create mode 100644 SD-VBS/benchmarks/mser/src/c/script_mser.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/Makefile create mode 100755 SD-VBS/benchmarks/mser/src/matlab/TIMESTAMP create mode 100755 SD-VBS/benchmarks/mser/src/matlab/erfill.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/erfill.mex.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/erfill.mexa64 create mode 100755 SD-VBS/benchmarks/mser/src/matlab/erfill.mexglx create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mexutils.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser.mex.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser.mexa64 create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser.mexglx create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser_compile.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser_demo2.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/mser_demo3.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/Makefile create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/TIMESTAMP create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/erfill.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/erfill.mex.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/erfill.mexa64 create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mexutils.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser.mexa64 create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser_compile.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser_demo2.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/mser_demo3.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/overview_mser.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/old/script_run_profile.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/overview_mser.m create mode 100755 SD-VBS/benchmarks/mser/src/matlab/script_run_profile.m (limited to 'SD-VBS/benchmarks/mser/src') diff --git a/SD-VBS/benchmarks/mser/src/c/mser.c b/SD-VBS/benchmarks/mser/src/c/mser.c new file mode 100644 index 0000000..d886d7a --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/c/mser.c @@ -0,0 +1,714 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +/*** +% MSER Maximally Stable Extremal Regions +% R=MSER(I,DELTA) computes the Maximally Stable Extremal Regions +% (MSER) of image I with stability threshold DELTA. I is any +% array of class UINT8, while DELTA is a scalar of the same class. +% R is an index set (of class UINT32) which enumerates the +% representative pixels of the detected regions. +% +% A region R can be recovered from a representative pixel X as the +% connected component of the level set {Y:I(Y) <= I(X)} which +% contains X. +***/ + + +#include "mser.h" +#include + +/* advance N-dimensional subscript */ +void +adv(iArray *dims, int ndims, iArray *subs_pt) +{ + int d = 0 ; + while(d < ndims) + { + sref(subs_pt,d) = sref(subs_pt,d) + 1; + if( sref(subs_pt,d) < sref(dims,d) ) + return ; + sref(subs_pt,d++) = 0 ; + } +} + +/** driver **/ +I2D* mser(I2D* I, int in_delta, + iArray* subs_pt, iArray* nsubs_pt, iArray* strides_pt, iArray* visited_pt, iArray* dims, + uiArray* joins_pt, + region_t* regions_pt, + pair_t* pairs_pt, + node_t* forest_pt, + ulliArray* acc_pt, ulliArray* ell_pt, + I2D* out) +{ + idx_t i, rindex=0; + int k; + int nout = 1; + + int OUT_REGIONS=0; + int OUT_ELL = 1; + int OUT_PARENTS = 2; + int OUT_AREA = 3; + int BUCKETS = 256; + + //I2D* out; + + int IN_I = 0; + int IN_DELTA = 1; + + /* 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 ; + + //iArray* subs_pt ; /* N-dimensional subscript + //iArray* nsubs_pt ; /* diff-subscript to point to neigh. +// uiArray* strides_pt ; /* strides to move in image array +// uiArray* 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 */ +// iArray* dims ; /* dimensions + int njoins = 0 ; /* number of join ops */ + + I2D* 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 + int regions_pt_size; + int pairs_pt_size; + int forest_pt_size; + + /* ellipses fitting */ + //ulliArray* acc_pt ; /* accumulator to integrate region moments + //ulliArray* ell_pt ; /* ellipses parameters + int gdl ; /* number of parameters of an ellipse */ + //uiArray* joins_pt ; /* sequence of joins + + delta = 0; + delta = in_delta; + + /* get dimensions */ + + nel = I->height*I->width; /* number of elements of src image */ + ndims = 2; + //dims = malloc(sizeof(iArray) + sizeof(int)*ndims); + I_pt = I; + + sref(dims,0) = I->height; + sref(dims,1) = I->width; + + /* allocate stuff */ + //subs_pt = malloc(sizeof(iArray) + sizeof(int)*ndims); + //nsubs_pt = malloc(sizeof(iArray) + sizeof(int)*ndims); + + //strides_pt = malloc(sizeof(uiArray)+sizeof(unsigned int)*ndims); + //visited_pt = malloc(sizeof(uiArray) + sizeof(unsigned int)*nel); + //joins_pt = malloc(sizeof(uiArray) + sizeof(unsigned int)*nel); + + //regions_pt = (region_t*)malloc(sizeof(region_t)*nel); + regions_pt_size = nel; + + //pairs_pt = (pair_t*)malloc(sizeof(pair_t)*nel); + pairs_pt_size = nel; + + //forest_pt = (node_t*)malloc(sizeof(node_t)*nel); + forest_pt_size = nel; + + /* compute strides to move into the N-dimensional image array */ + sref(strides_pt,0) = 1; + for(k = 1 ; k < ndims ; ++k) + { + sref(strides_pt,k) = sref(strides_pt,k-1) * sref(dims,k-1) ; + } + + /* sort pixels in increasing order of intensity: using Bucket Sort */ + { + int unsigned buckets [BUCKETS] ; + memset(buckets, 0, sizeof(int unsigned)*BUCKETS) ; + + for(i = 0 ; i < nel ; ++i) + { + val_t v = asubsref(I_pt,i) ; + ++buckets[v] ; + } + + for(i = 1 ; i < BUCKETS ; ++i) + { + arrayref(buckets,i) += arrayref(buckets,i-1) ; + } + + for(i = nel ; i >= 1 ; ) + { + val_t v = asubsref(I_pt,--i) ; + idx_t j = --buckets[v] ; + pairs_pt[j].value = v ; + pairs_pt[j].index = i ; + } + } + + /* 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 + * -------------------------------------------------------------- */ + + 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) + { + sref(nsubs_pt,k) = -1 ; + sref(subs_pt,k) = temp / sref(strides_pt,k) ; + temp = temp % sref(strides_pt,k) ; + } + } + + /* process neighbors of xi */ + while(1) + { + int good = 1 ; + 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 = sref(nsubs_pt,k) + sref(subs_pt,k) ; + good &= 0 <= temp && temp < sref(dims,k) ; + nindex += temp * sref(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 ) + { + sref(visited_pt,nvisited++) = rindex ; + rindex = forest_pt[rindex].shortcut ; + } + while( nvisited-- ) + { + forest_pt [ sref(visited_pt,nvisited) ] .shortcut = rindex ; + } + + /* find the root of the neighbor */ + nrindex = nindex ; + nvisited = 0 ; + while( forest_pt[nrindex].shortcut != nrindex ) + { + sref(visited_pt, nvisited++) = nrindex ; + nrindex = forest_pt[nrindex].shortcut ; + } + while( nvisited-- ) + { + forest_pt [ sref(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 = asubsref(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 + + sref(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 ; + } + + /* annote join operation for post-processing */ + sref(joins_pt,njoins++) = nrindex ; + } + } + + } /* neighbor done */ + + /* move to next neighbor */ + k = 0 ; + sref(nsubs_pt,k) = sref(nsubs_pt,k) + 1; + while( sref(nsubs_pt, k) > 1) + { + sref(nsubs_pt,k++) = -1 ; + if(k == ndims) goto done_all_neighbors ; + sref(nsubs_pt,k) = sref(nsubs_pt,k) + 1; + } + } /* 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 = asubsref(I_pt,rindex) ; + regions_pt [ner] .area = forest_pt [rindex] .area ; + regions_pt [ner] .area_top = nel ; + regions_pt [ner] .area_bot = 0 ; + ++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 = asubsref(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(1) + { + 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 = (area_top - area_bot) / (area*1.0) ; + + /* 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 ; + } + + + /* ----------------------------------------------------------------- + * 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 = (parea - area)/(area*1.0) ; + + if(change < 0.5) + { + ++ndup ; + goto remove_this_region ; + } + + } /* drop duplicates */ + } + + continue ; + remove_this_region : + regions_pt[i].maxstable = 0 ; + --nmer ; + + } /* next region to cleanup */ + + if(0) + { + printf(" Bad regions: %d\n", nbad ) ; + printf(" Small regions: %d\n", nsmall ) ; + printf(" Big regions: %d\n", nbig ) ; + printf(" Duplicated regions: %d\n", ndup ) ; + } + } + +/* printf("Cleaned-up regions: %d (%.1f%%)\n", + nmer, 100.0 * (double) nmer / ner) ; +*/ + /* ----------------------------------------------------------------- + * Fit ellipses + * -------------------------------------------------------------- */ + //ell_pt = 0 ; + //memset(ell_pt, sizeof(ulliArray) + sizeof(acc_t)*gdl*nmer, 0) ; + if (nout >= 1) + { + int midx = 1 ; + int d, index, j ; + + /* enumerate maxstable regions */ + for(i = 0 ; i < ner ; ++i) + { + if(! regions_pt [i].maxstable) continue ; + regions_pt [i].maxstable = midx++ ; + } + + /* allocate space */ + //acc_pt = malloc(sizeof(ulliArray) + sizeof(acc_t)*nel) ; + //printf("nmer = %d\n", nmer); + //ell_pt = malloc(sizeof(ulliArray) + sizeof(acc_t)*gdl*nmer) ; + + /* clear accumulators */ + for(d=0; d<(gdl*nmer); d++) + sref(ell_pt,d) = 0; + + /* for each gdl */ + for(d = 0 ; d < gdl ; ++d) + { + /* initalize parameter */ + int counter_i; + for(counter_i=0; counter_i j) + { + i -= j + 1 ; + j ++ ; + } + + /* add x_i * x_j */ + for(index = 0 ; index < nel ; ++ index) + { + sref(acc_pt,index) = sref(subs_pt,i) * sref(subs_pt,j) ; + adv(dims, ndims, subs_pt) ; + } + } + + /* integrate parameter */ + for(i = 0 ; i < njoins ; ++i) + { + idx_t index = sref(joins_pt,i); + idx_t parent = forest_pt [ index ].parent ; + sref(acc_pt,parent) += sref(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 ; + sref(ell_pt,d + gdl*region) = sref(acc_pt, regions_pt[i].index) ; + } + + /* next gdl */ + } + //free(acc_pt) ; + //free(ell_pt) ; + } + + /* ----------------------------------------------------------------- + * Save back and exit + * -------------------------------------------------------------- */ + + /* + * Save extremal regions + */ + { + int dims[2], j=0; + I2D* pt ; + dims[0] = nmer ; + //out = iMallocHandle(1, nmer); + out->height = 1; + out->width = nmer; + pt = out; + for (i = 0 ; i < ner ; ++i) + { + if( regions_pt[i].maxstable ) + { + /* adjust for MATLAB index compatibility */ +// *pt++ = regions_pt[i].index + 1 ; + asubsref(pt,j++) = regions_pt[i].index + 1 ; + } + } + } + + /* free stuff */ + //free(dims); + //free( forest_pt ) ; + //free( pairs_pt ) ; + //free( regions_pt ) ; + //free( visited_pt ) ; + //free( strides_pt ) ; + //free( nsubs_pt ) ; + //free( subs_pt ) ; + //free( joins_pt ) ; + + return out; +} + + + diff --git a/SD-VBS/benchmarks/mser/src/c/mser.h b/SD-VBS/benchmarks/mser/src/c/mser.h new file mode 100644 index 0000000..8876311 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/c/mser.h @@ -0,0 +1,83 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#ifndef _MSER_ +#define _MSER_ + +#define sref(a,i) a->data[i] + +#include "sdvbs_common.h" +#define NMER_MAX 756 + +typedef int val_t; + +typedef struct +{ + int width; + int data[]; +}iArray; + +typedef struct +{ + int width; + unsigned int data[]; +}uiArray; + +typedef struct +{ + int width; + long long int unsigned data[]; +}ulliArray; + +#define MIN(a,b) (ab)?a:b + +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 ; + +int script_mser(); +I2D* mser(I2D* I, int in_delta, + iArray* subs_pt, iArray* nsubs_pt, iArray* strides_pt, iArray* visited_pt, iArray* dims, + uiArray* joins_pt, + region_t* regions_pt, + pair_t* pairs_pt, + node_t* forest_pt, + ulliArray* acc_pt, ulliArray* ell_pt, + I2D* out); + +#endif + + diff --git a/SD-VBS/benchmarks/mser/src/c/script_mser.c b/SD-VBS/benchmarks/mser/src/c/script_mser.c new file mode 100644 index 0000000..d4a98cd --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/c/script_mser.c @@ -0,0 +1,120 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#include "mser.h" +#include +#include "extra.h" +#define min(a,b) (ab)?a:b + +int main(int argc, char* argv[]) +{ + SET_UP + int which_image; + int i, j, k; + I2D *idx; + I2D *I; + I2D *It; + I2D *out; + int rows=196, cols=98; + int minVal = 1000; + int maxVal = -1000; + int lev = 10; + + char im1[100], im2[100]; + + iArray *subs_pt, *nsubs_pt, *strides_pt, *visited_pt, *dims; + uiArray* joins_pt; + ulliArray *acc_pt, *ell_pt; + region_t* regions_pt; + pair_t* pairs_pt; + node_t* forest_pt; + + int ndims, nel, gdl, nmer; + + printf("Input Image: "); + scanf("%s", im1); + + I = readImage(im1); + + rows = I->height; + cols = I->width; + + It = readImage(im1); + + k = 0; + for(i=0; iheight * It->width; + gdl = ndims * (ndims+1)/2 + ndims; + nmer = NMER_MAX; + + dims = malloc(sizeof(iArray) + sizeof(int)*ndims); + /* allocate stuff */ + subs_pt = malloc(sizeof(iArray) + sizeof(int)*ndims); + nsubs_pt = malloc(sizeof(iArray) + sizeof(int)*ndims); + strides_pt = malloc(sizeof(uiArray)+sizeof(unsigned int)*ndims); + visited_pt = malloc(sizeof(uiArray) + sizeof(unsigned int)*nel); + joins_pt = malloc(sizeof(uiArray) + sizeof(unsigned int)*nel); + + regions_pt = (region_t*)malloc(sizeof(region_t)*nel); + pairs_pt = (pair_t*)malloc(sizeof(pair_t)*nel); + forest_pt = (node_t*)malloc(sizeof(node_t)*nel); + + acc_pt = malloc(sizeof(ulliArray) + sizeof(acc_t)*nel) ; + ell_pt = malloc(sizeof(ulliArray) + sizeof(acc_t)*gdl*nmer) ; + + + out = iMallocHandle(1, nmer); + printf("start\n"); + for_each_job{ + idx = mser(It, 2, subs_pt, nsubs_pt, strides_pt, visited_pt, dims, + joins_pt, + regions_pt, + pairs_pt, + forest_pt, + acc_pt, ell_pt, + out); + } + printf("end..\n"); + +#ifdef CHECK + /** Self checking - use expected.txt from data directory **/ + { + int tol, ret=0; + tol = 1; +#ifdef GENERATE_OUTPUT + writeMatrix(idx, argv[1]); +#endif + ret = selfCheck(idx, "expected_C.txt", tol); + if (ret == -1) + printf("Error in MSER\n"); + } + /** Self checking done **/ +#endif + free(dims); + free( forest_pt ) ; + free( pairs_pt ) ; + free( regions_pt ) ; + free( visited_pt ) ; + free( strides_pt ) ; + free( nsubs_pt ) ; + free( subs_pt ) ; + free( joins_pt ) ; + free( acc_pt ) ; + free( ell_pt ) ; + iFreeHandle(idx); + iFreeHandle(I); + iFreeHandle(It); + WRITE_TO_FILE + return 0; +} + diff --git a/SD-VBS/benchmarks/mser/src/matlab/Makefile b/SD-VBS/benchmarks/mser/src/matlab/Makefile new file mode 100755 index 0000000..29c0982 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/Makefile @@ -0,0 +1,123 @@ +# file: Makefile +# author: Andrea Vedaldi +# description: Build mex files + +# -------------------------------------------------------------------- +# +# -------------------------------------------------------------------- + +# Determine on the flight the system we are running on +Darwin_ARCH := mac +Linux_ARCH := glx +ARCH := $($(shell uname)_ARCH) + +mac_CFLAGS := -I. -pedantic -Wall -Wno-long-long +mac_MEX_CFLAGS := -g -O CFLAGS='$$CFLAGS $(mac_CFLAGS)' +mac_MEX_SUFFIX := mexmac + +glx_CFLAGS := -I. -pedantic -Wall -Wno-long-long +glx_MEX_CFLAGS := -g -O CFLAGS='$$CFLAGS $(glx_CFLAGS)' +glx_MEX_SUFFIX := mexglx + +MEX_SUFFIX := $($(ARCH)_MEX_SUFFIX) +MEX_CFLAGS := $($(ARCH)_MEX_CFLAGS) + +VER := 0.4 +DIST := mser-$(VER) +BINDIST := $(DIST)-$(ARCH) + +# -------------------------------------------------------------------- +# +# -------------------------------------------------------------------- + +vpath %.mex.c . + +src := $(wildcard *.mex.c) +msrc := $(wildcard *.m) +stem := $(notdir $(basename $(basename $(src)))) +tgt := $(addprefix ./, $(addsuffix .$(MEX_SUFFIX),$(stem))) + +%.$(MEX_SUFFIX) : %.mex.c + mex -I. $(MEX_CFLAGS) $< -output $* + +.PHONY: all +all: $(tgt) + +.PHONY: info +info : + @echo src = $(src) + @echo stem = $(stem) + @echo tgt = $(tgt) + +# PDF documentation +.PHONY: doc +doc: mser.html doc/mser.pdf + +mser.html : $(msrc) + mdoc --output=mser.html . \ + --exclude='.*(_demo|_compile).*.m' + +.PHONY: clean +clean: + rm -f $(tgt) + find . -name '.DS_Store' -exec rm -f \{\} \; + find . -name '.gdb_history' -exec rm -f \{\} \; + find . -name '*~' -exec rm -f \{\} \; + find . -name '*.bak' -exec rm -f \{\} \; + make -C doc/figures clean + +.PHONY: distclean +distclean: clean + rm -f *.mexmac *.mexglx + rm -f mser.html + rm -f mser-*.tar.gz + rm -f doc/*.log + rm -f doc/*.aux + rm -f doc/*.toc + rm -f doc/*.bbl + rm -f doc/*.blg + rm -f doc/*.out + rm -f $(DIST).tar.gz + rm -f $(BINDIST).tar.gz + rm -rf $(BINDIST) + +.PHONY: dist +dist: distclean + echo Version $(VER) >TIMESTAMP + echo Archive created on `date` >>TIMESTAMP + d=$(notdir $(CURDIR)) ; \ + tar chzvf $(DIST).tar.gz \ + --exclude mser_demo4.m \ + --exclude data/seq.avi \ + --exclude results \ + ../$${d} + +.PHONY: bindist +bindist: all + test -e $(BINDIST) || mkdir $(BINDIST) + cp *.$(MEX_SUFFIX) $(BINDIST) + cd $(BINDIST) ; strip -S *.$(MEX_SUFFIX) + tar chzvf $(BINDIST).tar.gz $(BINDIST) + +.PHONY: autorights +autorights: + autorights . \ + --verbose \ + --recursive \ + --template cal \ + --years 2006 \ + --authors "Andrea Vedaldi (UCLA VisionLab)" \ + --program "Video Extremal Regions" + +doc/mser.pdf : doc/*.tex doc/*.bib doc/figures/*.fig + make -C doc/figures all + cd doc ; \ + for k in 1 2 3 ; \ + do \ + pdflatex -file-line-error-style -interaction batchmode \ + mser.tex ; \ + if test "$$k" = '1' ; \ + then \ + bibtex mser.aux ; \ + fi ; \ + done diff --git a/SD-VBS/benchmarks/mser/src/matlab/TIMESTAMP b/SD-VBS/benchmarks/mser/src/matlab/TIMESTAMP new file mode 100755 index 0000000..1de1720 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/TIMESTAMP @@ -0,0 +1,2 @@ +Version 0.4 +Archive created on Wed Feb 7 11:08:47 PST 2007 diff --git a/SD-VBS/benchmarks/mser/src/matlab/erfill.m b/SD-VBS/benchmarks/mser/src/matlab/erfill.m new file mode 100755 index 0000000..6e11fc6 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/erfill.m @@ -0,0 +1,13 @@ +% ERFILL Fill extremal region +% MEMBERS=ERFILL(I,ER) returns the list MEMBERS of the pixels which +% belongs to the extremal region represented by the pixel ER. +% +% The selected region is the one that contains pixel ER and of +% inensity I(ER). +% +% I must be of class UINT8 and ER must be a (scalar) index of the +% region representative point. +% +% See also MSER(). + + diff --git a/SD-VBS/benchmarks/mser/src/matlab/erfill.mex.c b/SD-VBS/benchmarks/mser/src/matlab/erfill.mex.c new file mode 100755 index 0000000..893d346 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/erfill.mex.c @@ -0,0 +1,223 @@ +/* file: erfill.mex.c +** description: Extremal Regions filling +** 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)) + +typedef char unsigned val_t ; +typedef int unsigned idx_t ; +typedef long long int unsigned acc_t ; + +/* 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_ER} ; + enum {OUT_MEMBERS} ; + + idx_t i ; + int k, nel, ndims ; + int const * dims ; + val_t const * I_pt ; + int last = 0 ; + int last_expanded = 0 ; + val_t value = 0 ; + + double const * er_pt ; + + 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 */ + val_t* visited_pt ; /* flag */ + idx_t* members_pt ; /* region members */ + + /** ----------------------------------------------------------------- + ** 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(!uIsRealScalar(in[IN_ER])) { + mexErrMsgTxt("ER must be a DOUBLE scalar.") ; + } + + /* 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(val_t) * nel ) ; + members_pt = mxMalloc( sizeof(idx_t) * nel ) ; + + er_pt = mxGetPr(in[IN_ER]) ; + + /* 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] ; + } + + /* load first pixel */ + memset(visited_pt, 0, sizeof(val_t) * nel) ; + { + idx_t idx = (idx_t) *er_pt ; + if( idx < 1 || idx > nel ) { + char buff[80] ; + snprintf(buff,80,"ER=%d out of range [1,%d]",idx,nel) ; + mexErrMsgTxt(buff) ; + } + members_pt [last++] = idx - 1 ; + } + value = I_pt[ members_pt[0] ] ; + + /* ----------------------------------------------------------------- + * Fill region + * -------------------------------------------------------------- */ + while(last_expanded < last) { + + /* pop next node xi */ + idx_t index = members_pt[last_expanded++] ; + + /* 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] ; + } + + /* process neighbor + 1 - the pixel is within image boundaries; + 2 - the pixel is indeed different from the current node + (this happens when nsub=(0,0,...,0)); + 3 - the pixel has value not greather than val + is a pixel older than xi + 4 - the pixel has not been visited yet + */ + if(good + && nindex != index + && I_pt [nindex] <= value + && ! visited_pt [nindex] ) { + + /* mark as visited */ + visited_pt [nindex] = 1 ; + + /* add to list */ + members_pt [last++] = nindex ; + } + + /* 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 : ; + } /* goto pop next member */ + + /* + * Save results + */ + { + int dims[2] ; + int unsigned * pt ; + dims[0] = last ; + out[OUT_MEMBERS] = mxCreateNumericArray(1,dims,mxUINT32_CLASS,mxREAL); + pt = mxGetData(out[OUT_MEMBERS]) ; + for (i = 0 ; i < last ; ++i) { + *pt++ = members_pt[i] + 1 ; + } + } + + /* free stuff */ + mxFree( members_pt ) ; + mxFree( visited_pt ) ; + mxFree( strides_pt ) ; + mxFree( nsubs_pt ) ; + mxFree( subs_pt ) ; +} diff --git a/SD-VBS/benchmarks/mser/src/matlab/erfill.mexa64 b/SD-VBS/benchmarks/mser/src/matlab/erfill.mexa64 new file mode 100755 index 0000000..bc54d65 Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/erfill.mexa64 differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/erfill.mexglx b/SD-VBS/benchmarks/mser/src/matlab/erfill.mexglx new file mode 100755 index 0000000..8eec110 Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/erfill.mexglx differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/mexutils.c b/SD-VBS/benchmarks/mser/src/matlab/mexutils.c new file mode 100755 index 0000000..0fc664b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/mexutils.c @@ -0,0 +1,111 @@ +/* file: mexutils.c +** author: Andrea Vedaldi +** description: Utility functions to write MEX files. +**/ + +#include"mex.h" + +#undef M_PI +#define M_PI 3.14159265358979 + +/** @brief Is scalar? + ** + ** @return @c true if the array @a A is a scalar. + **/ +int +uIsScalar(const mxArray* A) +{ + return + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + mxGetM(A) == 1 && + mxGetN(A) == 1 ; +} + +/** @brief Is real scalar? + ** + ** @return @c true if the array @a A is a real scalar. + **/ +int +uIsRealScalar(const mxArray* A) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + mxGetM(A) == 1 && + mxGetN(A) == 1 ; +} + +/** @brief Is real matrix? + ** + ** The function checks wether the argument @a A is a real matrix. In + ** addition, if @a M >= 0, it checks wether the number of rows is + ** equal to @a M and, if @a N >= 0, if the number of columns is equal + ** to @a N. + ** + ** @param M number of rows. + ** @param N number of columns. + ** @return @c true if the array is a real matrix with the specified format. + **/ +int +uIsRealMatrix(const mxArray* A, int M, int N) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + ((M>=0)?(mxGetM(A) == M):1) && + ((N>=0)?(mxGetN(A) == N):1) ; +} + +/** @brief Is real vector? + ** + ** The function checks wether the argument @a V is a real vector. By + ** definiton, a matrix is a vector if one of its dimension is one. + ** In addition, if @a D >= 0, it checks wether the dimension of the + ** vecotr is equal to @a D. + ** + ** @param D lenght of the vector. + ** @return @c true if the array is a real vector of the specified dimension. + **/ +int +uIsRealVector(const mxArray* V, int D) +{ + int M = mxGetM(V) ; + int N = mxGetN(V) ; + int is_vector = (N == 1) || (M == 1) ; + + return + mxIsDouble(V) && + !mxIsComplex(V) && + mxGetNumberOfDimensions(V) == 2 && + is_vector && + ( D < 0 || N == D || M == D) ; +} + + +/** @brief Is a string? + ** + ** The function checks wether the array @a S is a string. If + ** @a L is non-negative, it also check wether the strign has + ** length @a L. + ** + ** @return @a c true if S is a string of the specified length. + **/ +int +uIsString(const mxArray* S, int L) +{ + int M = mxGetM(S) ; + int N = mxGetN(S) ; + + return + mxIsChar(S) && + M == 1 && + (L < 0 || N == L) ; +} + +/** + ** + **/ + diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c b/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c new file mode 100755 index 0000000..8473afe --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/mser.mex.c @@ -0,0 +1,815 @@ +/* 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 BUCKETS 256 + +#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 = 0 ; /* 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_BUCKET_SORT + 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 [BUCKETS] ; + int unsigned v; + memset(buckets, 0, sizeof(int unsigned)*BUCKETS) ; + for(i = 0 ; i < nel ; ++i) { + v = (unsigned int)I_pt [i] ; + ++buckets[v] ; + } + for(i = 1 ; i < BUCKETS ; ++i) { + buckets[i] += buckets[i-1] ; + } + for(i = nel ; i >= 1 ; ) { + 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 ; + } + + /* 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 ) ; +} diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser.mexa64 b/SD-VBS/benchmarks/mser/src/matlab/mser.mexa64 new file mode 100755 index 0000000..e3ca56b Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/mser.mexa64 differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser.mexglx b/SD-VBS/benchmarks/mser/src/matlab/mser.mexglx new file mode 100755 index 0000000..b69e1d2 Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/mser.mexglx differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser_compile.m b/SD-VBS/benchmarks/mser/src/matlab/mser_compile.m new file mode 100755 index 0000000..5e3562b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/mser_compile.m @@ -0,0 +1,7 @@ +function mser_compile(type) +% MSER_COMPILE Compile MEX files + +opts = { '-O', '-I.' } ; + +mex('mser.mex.c','-output', 'mser',opts{:}) ; +mex('erfill.mex.c','-output', 'erfill',opts{:}) ; diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser_demo2.m b/SD-VBS/benchmarks/mser/src/matlab/mser_demo2.m new file mode 100755 index 0000000..37a4ed1 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/mser_demo2.m @@ -0,0 +1,62 @@ +% MSER_DEMO2 Demonstrate MSER code + +% 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. + +I = load('clown') ; I = uint8(I.X) ; +figure(1) ; imagesc(I) ; colormap gray; hold on ; + +[M,N] = size(I) ; +i = double(i) ; +j = double(j) ; + +[r,ell] = mser(I,5) ; + +r=double(r) ; + +[i,j]=ind2sub(size(I),r) ; +plot(j,i,'r*') ; + +ell = ell([2 1 5 4 3],:) ; +plotframe(ell); + +figure(2) ; + +clear MOV ; +K = size(ell,2) ; +for k=1:K + clf ; + sel = erfill(I,r(k)) ; + mask = zeros(M,N) ; mask(sel) =1 ; + imagesc(cat(3,I,255*uint8(mask),I)) ; colormap gray ; hold on ; + set(gca,'position',[0 0 1 1]) ; axis off ; axis equal ; + plot(j(k),i(k),'r*') ; + plotframe(ell(:,k),'color','r') ; + MOV(k) = getframe(gca) ; +end diff --git a/SD-VBS/benchmarks/mser/src/matlab/mser_demo3.m b/SD-VBS/benchmarks/mser/src/matlab/mser_demo3.m new file mode 100755 index 0000000..4669437 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/mser_demo3.m @@ -0,0 +1,117 @@ +% MSER_DEMO3 Demonstrates MSER on a volumetric image + +% 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. + +% -------------------------------------------------------------------- +% Create data +% -------------------------------------------------------------------- + +% volumetric coordinate (x,y,z) +x = linspace(-1,1,50) ; +[x,y,z] = meshgrid(x,x,x) ; + +% create funny volumetric image +I = sin(4*x).*cos(4*y).*sin(z) ; +I = I-min(I(:)) ; +I = I/max(I(:)) ; + +% quantize the image in 10 levels +lev = 10 ; +I = lev*I ; +Ir = round(I) ; + +% -------------------------------------------------------------------- +% Compute regions +% -------------------------------------------------------------------- +[idx,ell,p] = mser(uint8(Ir),1); + +% -------------------------------------------------------------------- +% Plots +% -------------------------------------------------------------------- + +% The image is quantized; store in LEV its range. +lev = unique(Ir(idx)) ; + +figure(100); clf; +K=min(length(lev),4) ; + +r=.99 ; + +% one level per time +for k=1:K + tightsubplot(K,k) ; + [i,j,m] = ind2sub(size(I), idx(Ir(idx)==lev(k)) ) ; + + % compute level set of level LEV(k) + Is = double(Ir<=lev(k)) ; + + p1 = patch(isosurface(Is,r), ... + 'FaceColor','blue','EdgeColor','none') ; + p2 = patch(isocaps(Is,r),... + 'FaceColor','interp','EdgeColor','none') ; + isonormals(I,p1) + hold on ; + + view(3); axis vis3d tight + camlight; lighting phong ; + + % find regions that have this level + sel = find( Ir(idx) == lev(k) ) ; + + % plot fitted ellipsoid + for r=sel' + E = ell(:,r) ; + c = E(1:3) ; + A = zeros(3) ; + A(1,1) = E(4) ; + A(1,2) = E(5) ; + A(2,2) = E(6) ; + A(1,3) = E(7) ; + A(2,3) = E(8) ; + A(3,3) = E(9) ; + + A = A + A' - diag(diag(A)) ; + + % correct var. order + perm = [0 1 0 ; 1 0 0 ; 0 0 1] ; + A = perm*A*perm ; + + [V,D] = eig(A) ; + A = 2.5*V*sqrt(D) ; + + [x,y,z]=sphere ; + [P,Q]=size(x) ; + X=A*[x(:)';y(:)';z(:)'] ; + x=reshape(X(1,:),P,Q)+c(2) ; + y=reshape(X(2,:),P,Q)+c(1) ; + z=reshape(X(3,:),P,Q)+c(3) ; + surf(x,y,z,'FaceAlpha',.5) ; + end +end diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/Makefile b/SD-VBS/benchmarks/mser/src/matlab/old/Makefile new file mode 100755 index 0000000..29c0982 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/Makefile @@ -0,0 +1,123 @@ +# file: Makefile +# author: Andrea Vedaldi +# description: Build mex files + +# -------------------------------------------------------------------- +# +# -------------------------------------------------------------------- + +# Determine on the flight the system we are running on +Darwin_ARCH := mac +Linux_ARCH := glx +ARCH := $($(shell uname)_ARCH) + +mac_CFLAGS := -I. -pedantic -Wall -Wno-long-long +mac_MEX_CFLAGS := -g -O CFLAGS='$$CFLAGS $(mac_CFLAGS)' +mac_MEX_SUFFIX := mexmac + +glx_CFLAGS := -I. -pedantic -Wall -Wno-long-long +glx_MEX_CFLAGS := -g -O CFLAGS='$$CFLAGS $(glx_CFLAGS)' +glx_MEX_SUFFIX := mexglx + +MEX_SUFFIX := $($(ARCH)_MEX_SUFFIX) +MEX_CFLAGS := $($(ARCH)_MEX_CFLAGS) + +VER := 0.4 +DIST := mser-$(VER) +BINDIST := $(DIST)-$(ARCH) + +# -------------------------------------------------------------------- +# +# -------------------------------------------------------------------- + +vpath %.mex.c . + +src := $(wildcard *.mex.c) +msrc := $(wildcard *.m) +stem := $(notdir $(basename $(basename $(src)))) +tgt := $(addprefix ./, $(addsuffix .$(MEX_SUFFIX),$(stem))) + +%.$(MEX_SUFFIX) : %.mex.c + mex -I. $(MEX_CFLAGS) $< -output $* + +.PHONY: all +all: $(tgt) + +.PHONY: info +info : + @echo src = $(src) + @echo stem = $(stem) + @echo tgt = $(tgt) + +# PDF documentation +.PHONY: doc +doc: mser.html doc/mser.pdf + +mser.html : $(msrc) + mdoc --output=mser.html . \ + --exclude='.*(_demo|_compile).*.m' + +.PHONY: clean +clean: + rm -f $(tgt) + find . -name '.DS_Store' -exec rm -f \{\} \; + find . -name '.gdb_history' -exec rm -f \{\} \; + find . -name '*~' -exec rm -f \{\} \; + find . -name '*.bak' -exec rm -f \{\} \; + make -C doc/figures clean + +.PHONY: distclean +distclean: clean + rm -f *.mexmac *.mexglx + rm -f mser.html + rm -f mser-*.tar.gz + rm -f doc/*.log + rm -f doc/*.aux + rm -f doc/*.toc + rm -f doc/*.bbl + rm -f doc/*.blg + rm -f doc/*.out + rm -f $(DIST).tar.gz + rm -f $(BINDIST).tar.gz + rm -rf $(BINDIST) + +.PHONY: dist +dist: distclean + echo Version $(VER) >TIMESTAMP + echo Archive created on `date` >>TIMESTAMP + d=$(notdir $(CURDIR)) ; \ + tar chzvf $(DIST).tar.gz \ + --exclude mser_demo4.m \ + --exclude data/seq.avi \ + --exclude results \ + ../$${d} + +.PHONY: bindist +bindist: all + test -e $(BINDIST) || mkdir $(BINDIST) + cp *.$(MEX_SUFFIX) $(BINDIST) + cd $(BINDIST) ; strip -S *.$(MEX_SUFFIX) + tar chzvf $(BINDIST).tar.gz $(BINDIST) + +.PHONY: autorights +autorights: + autorights . \ + --verbose \ + --recursive \ + --template cal \ + --years 2006 \ + --authors "Andrea Vedaldi (UCLA VisionLab)" \ + --program "Video Extremal Regions" + +doc/mser.pdf : doc/*.tex doc/*.bib doc/figures/*.fig + make -C doc/figures all + cd doc ; \ + for k in 1 2 3 ; \ + do \ + pdflatex -file-line-error-style -interaction batchmode \ + mser.tex ; \ + if test "$$k" = '1' ; \ + then \ + bibtex mser.aux ; \ + fi ; \ + done diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/TIMESTAMP b/SD-VBS/benchmarks/mser/src/matlab/old/TIMESTAMP new file mode 100755 index 0000000..1de1720 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/TIMESTAMP @@ -0,0 +1,2 @@ +Version 0.4 +Archive created on Wed Feb 7 11:08:47 PST 2007 diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/erfill.m b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.m new file mode 100755 index 0000000..6e11fc6 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.m @@ -0,0 +1,13 @@ +% ERFILL Fill extremal region +% MEMBERS=ERFILL(I,ER) returns the list MEMBERS of the pixels which +% belongs to the extremal region represented by the pixel ER. +% +% The selected region is the one that contains pixel ER and of +% inensity I(ER). +% +% I must be of class UINT8 and ER must be a (scalar) index of the +% region representative point. +% +% See also MSER(). + + diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mex.c b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mex.c new file mode 100755 index 0000000..893d346 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mex.c @@ -0,0 +1,223 @@ +/* file: erfill.mex.c +** description: Extremal Regions filling +** 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)) + +typedef char unsigned val_t ; +typedef int unsigned idx_t ; +typedef long long int unsigned acc_t ; + +/* 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_ER} ; + enum {OUT_MEMBERS} ; + + idx_t i ; + int k, nel, ndims ; + int const * dims ; + val_t const * I_pt ; + int last = 0 ; + int last_expanded = 0 ; + val_t value = 0 ; + + double const * er_pt ; + + 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 */ + val_t* visited_pt ; /* flag */ + idx_t* members_pt ; /* region members */ + + /** ----------------------------------------------------------------- + ** 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(!uIsRealScalar(in[IN_ER])) { + mexErrMsgTxt("ER must be a DOUBLE scalar.") ; + } + + /* 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(val_t) * nel ) ; + members_pt = mxMalloc( sizeof(idx_t) * nel ) ; + + er_pt = mxGetPr(in[IN_ER]) ; + + /* 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] ; + } + + /* load first pixel */ + memset(visited_pt, 0, sizeof(val_t) * nel) ; + { + idx_t idx = (idx_t) *er_pt ; + if( idx < 1 || idx > nel ) { + char buff[80] ; + snprintf(buff,80,"ER=%d out of range [1,%d]",idx,nel) ; + mexErrMsgTxt(buff) ; + } + members_pt [last++] = idx - 1 ; + } + value = I_pt[ members_pt[0] ] ; + + /* ----------------------------------------------------------------- + * Fill region + * -------------------------------------------------------------- */ + while(last_expanded < last) { + + /* pop next node xi */ + idx_t index = members_pt[last_expanded++] ; + + /* 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] ; + } + + /* process neighbor + 1 - the pixel is within image boundaries; + 2 - the pixel is indeed different from the current node + (this happens when nsub=(0,0,...,0)); + 3 - the pixel has value not greather than val + is a pixel older than xi + 4 - the pixel has not been visited yet + */ + if(good + && nindex != index + && I_pt [nindex] <= value + && ! visited_pt [nindex] ) { + + /* mark as visited */ + visited_pt [nindex] = 1 ; + + /* add to list */ + members_pt [last++] = nindex ; + } + + /* 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 : ; + } /* goto pop next member */ + + /* + * Save results + */ + { + int dims[2] ; + int unsigned * pt ; + dims[0] = last ; + out[OUT_MEMBERS] = mxCreateNumericArray(1,dims,mxUINT32_CLASS,mxREAL); + pt = mxGetData(out[OUT_MEMBERS]) ; + for (i = 0 ; i < last ; ++i) { + *pt++ = members_pt[i] + 1 ; + } + } + + /* free stuff */ + mxFree( members_pt ) ; + mxFree( visited_pt ) ; + mxFree( strides_pt ) ; + mxFree( nsubs_pt ) ; + mxFree( subs_pt ) ; +} diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mexa64 b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mexa64 new file mode 100755 index 0000000..679e972 Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/old/erfill.mexa64 differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mexutils.c b/SD-VBS/benchmarks/mser/src/matlab/old/mexutils.c new file mode 100755 index 0000000..0fc664b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/mexutils.c @@ -0,0 +1,111 @@ +/* file: mexutils.c +** author: Andrea Vedaldi +** description: Utility functions to write MEX files. +**/ + +#include"mex.h" + +#undef M_PI +#define M_PI 3.14159265358979 + +/** @brief Is scalar? + ** + ** @return @c true if the array @a A is a scalar. + **/ +int +uIsScalar(const mxArray* A) +{ + return + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + mxGetM(A) == 1 && + mxGetN(A) == 1 ; +} + +/** @brief Is real scalar? + ** + ** @return @c true if the array @a A is a real scalar. + **/ +int +uIsRealScalar(const mxArray* A) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + mxGetM(A) == 1 && + mxGetN(A) == 1 ; +} + +/** @brief Is real matrix? + ** + ** The function checks wether the argument @a A is a real matrix. In + ** addition, if @a M >= 0, it checks wether the number of rows is + ** equal to @a M and, if @a N >= 0, if the number of columns is equal + ** to @a N. + ** + ** @param M number of rows. + ** @param N number of columns. + ** @return @c true if the array is a real matrix with the specified format. + **/ +int +uIsRealMatrix(const mxArray* A, int M, int N) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + ((M>=0)?(mxGetM(A) == M):1) && + ((N>=0)?(mxGetN(A) == N):1) ; +} + +/** @brief Is real vector? + ** + ** The function checks wether the argument @a V is a real vector. By + ** definiton, a matrix is a vector if one of its dimension is one. + ** In addition, if @a D >= 0, it checks wether the dimension of the + ** vecotr is equal to @a D. + ** + ** @param D lenght of the vector. + ** @return @c true if the array is a real vector of the specified dimension. + **/ +int +uIsRealVector(const mxArray* V, int D) +{ + int M = mxGetM(V) ; + int N = mxGetN(V) ; + int is_vector = (N == 1) || (M == 1) ; + + return + mxIsDouble(V) && + !mxIsComplex(V) && + mxGetNumberOfDimensions(V) == 2 && + is_vector && + ( D < 0 || N == D || M == D) ; +} + + +/** @brief Is a string? + ** + ** The function checks wether the array @a S is a string. If + ** @a L is non-negative, it also check wether the strign has + ** length @a L. + ** + ** @return @a c true if S is a string of the specified length. + **/ +int +uIsString(const mxArray* S, int L) +{ + int M = mxGetM(S) ; + int N = mxGetN(S) ; + + return + mxIsChar(S) && + M == 1 && + (L < 0 || N == L) ; +} + +/** + ** + **/ + diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c b/SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c new file mode 100755 index 0000000..48c788e --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/mser.mex.c @@ -0,0 +1,809 @@ +/* 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 ) ; +} diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mser.mexa64 b/SD-VBS/benchmarks/mser/src/matlab/old/mser.mexa64 new file mode 100755 index 0000000..fb9fb0b Binary files /dev/null and b/SD-VBS/benchmarks/mser/src/matlab/old/mser.mexa64 differ diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mser_compile.m b/SD-VBS/benchmarks/mser/src/matlab/old/mser_compile.m new file mode 100755 index 0000000..5e3562b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/mser_compile.m @@ -0,0 +1,7 @@ +function mser_compile(type) +% MSER_COMPILE Compile MEX files + +opts = { '-O', '-I.' } ; + +mex('mser.mex.c','-output', 'mser',opts{:}) ; +mex('erfill.mex.c','-output', 'erfill',opts{:}) ; diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo2.m b/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo2.m new file mode 100755 index 0000000..37a4ed1 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo2.m @@ -0,0 +1,62 @@ +% MSER_DEMO2 Demonstrate MSER code + +% 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. + +I = load('clown') ; I = uint8(I.X) ; +figure(1) ; imagesc(I) ; colormap gray; hold on ; + +[M,N] = size(I) ; +i = double(i) ; +j = double(j) ; + +[r,ell] = mser(I,5) ; + +r=double(r) ; + +[i,j]=ind2sub(size(I),r) ; +plot(j,i,'r*') ; + +ell = ell([2 1 5 4 3],:) ; +plotframe(ell); + +figure(2) ; + +clear MOV ; +K = size(ell,2) ; +for k=1:K + clf ; + sel = erfill(I,r(k)) ; + mask = zeros(M,N) ; mask(sel) =1 ; + imagesc(cat(3,I,255*uint8(mask),I)) ; colormap gray ; hold on ; + set(gca,'position',[0 0 1 1]) ; axis off ; axis equal ; + plot(j(k),i(k),'r*') ; + plotframe(ell(:,k),'color','r') ; + MOV(k) = getframe(gca) ; +end diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo3.m b/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo3.m new file mode 100755 index 0000000..4669437 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/mser_demo3.m @@ -0,0 +1,117 @@ +% MSER_DEMO3 Demonstrates MSER on a volumetric image + +% 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. + +% -------------------------------------------------------------------- +% Create data +% -------------------------------------------------------------------- + +% volumetric coordinate (x,y,z) +x = linspace(-1,1,50) ; +[x,y,z] = meshgrid(x,x,x) ; + +% create funny volumetric image +I = sin(4*x).*cos(4*y).*sin(z) ; +I = I-min(I(:)) ; +I = I/max(I(:)) ; + +% quantize the image in 10 levels +lev = 10 ; +I = lev*I ; +Ir = round(I) ; + +% -------------------------------------------------------------------- +% Compute regions +% -------------------------------------------------------------------- +[idx,ell,p] = mser(uint8(Ir),1); + +% -------------------------------------------------------------------- +% Plots +% -------------------------------------------------------------------- + +% The image is quantized; store in LEV its range. +lev = unique(Ir(idx)) ; + +figure(100); clf; +K=min(length(lev),4) ; + +r=.99 ; + +% one level per time +for k=1:K + tightsubplot(K,k) ; + [i,j,m] = ind2sub(size(I), idx(Ir(idx)==lev(k)) ) ; + + % compute level set of level LEV(k) + Is = double(Ir<=lev(k)) ; + + p1 = patch(isosurface(Is,r), ... + 'FaceColor','blue','EdgeColor','none') ; + p2 = patch(isocaps(Is,r),... + 'FaceColor','interp','EdgeColor','none') ; + isonormals(I,p1) + hold on ; + + view(3); axis vis3d tight + camlight; lighting phong ; + + % find regions that have this level + sel = find( Ir(idx) == lev(k) ) ; + + % plot fitted ellipsoid + for r=sel' + E = ell(:,r) ; + c = E(1:3) ; + A = zeros(3) ; + A(1,1) = E(4) ; + A(1,2) = E(5) ; + A(2,2) = E(6) ; + A(1,3) = E(7) ; + A(2,3) = E(8) ; + A(3,3) = E(9) ; + + A = A + A' - diag(diag(A)) ; + + % correct var. order + perm = [0 1 0 ; 1 0 0 ; 0 0 1] ; + A = perm*A*perm ; + + [V,D] = eig(A) ; + A = 2.5*V*sqrt(D) ; + + [x,y,z]=sphere ; + [P,Q]=size(x) ; + X=A*[x(:)';y(:)';z(:)'] ; + x=reshape(X(1,:),P,Q)+c(2) ; + y=reshape(X(2,:),P,Q)+c(1) ; + z=reshape(X(3,:),P,Q)+c(3) ; + surf(x,y,z,'FaceAlpha',.5) ; + end +end diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/overview_mser.m b/SD-VBS/benchmarks/mser/src/matlab/old/overview_mser.m new file mode 100755 index 0000000..b7fcf2b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/overview_mser.m @@ -0,0 +1,8 @@ +% OVERVIEW_MSER Maximally Stable Extremal Regions +% This is a MATLAB/MEX implementation of Maximally Stable Extremal +% Regions (MSER). You can: +% +% * Use MSER() to extract the maximally stable extremal regions from +% a given image. +% +% For practical coding example, see MSER_DEMO() and MSER_DEMO3(). diff --git a/SD-VBS/benchmarks/mser/src/matlab/old/script_run_profile.m b/SD-VBS/benchmarks/mser/src/matlab/old/script_run_profile.m new file mode 100755 index 0000000..bdb2c04 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/old/script_run_profile.m @@ -0,0 +1,137 @@ +function script_run_profile(dataDir, resultDir, type, common,toolDir) + +path(path, common); + +% MSER_DEMO Demonstrates MSER + +% 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. + +which_image = 3 ; + +% -------------------------------------------------------------------- +% Create data +% -------------------------------------------------------------------- +switch which_image + case 1 + I = rand(200,200) ; + I = imsmooth(I,10) ; + I = I-min(I(:)) ; + I = I/max(I(:)) ; + lev = 10 ; + I = uint8(round(I*lev)) ; + + case 2 + I = zeros(200,200) ; + I(50:150,50:150)=5 ; + I = imsmooth(I,10) ; + I = uint8(round(I)) ; + + case 3 + Files = dir([dataDir,'/1.bmp']); + I=readImage(fullfile(dataDir,Files(1).name)); + +% if(ndims(Image) > 2) +% +% for i=1:size(Image,1) +% for j=1:size(Image,2) +% I(i,j) = (Image(i,j,1) + Image(i,j,2)*6 + Image(i,j,3)*3)/10; +% end +% end +% +% else +% I = Image; +% end + disp(size(I)); +end + +% -------------------------------------------------------------------- +% Compute MSERs +% -------------------------------------------------------------------- + +%% Self check params +tol = 0.1; +elapsed = zeros(1,2); + +%% Timing +start = photonStartTiming; + +[idx,ell,p,a] = mser(uint8(I), 2) ; + +%% Timing +stop = photonEndTiming; + +temp = photonReportTiming(start, stop); +elapsed(1) = elapsed(1) + temp(1); +elapsed(2) = elapsed(2) + temp(2); + + %% Self checking + writeMatrix(idx, dataDir); + ret = selfCheck(idx, dataDir, tol); + if(ret == -1) + disp('Error in MSER'); + end + +%% Timing +photonPrintTiming(elapsed); + +%% -------------------------------------------------------------------- +%% Plots +%% -------------------------------------------------------------------- +%[i,j] = ind2sub(size(I),idx) ; +% +%figure(100) ; clf ; imagesc(I) ; hold on ; +%set(gca,'Position',[0 0 1 1]) ; +%plot(j,i,'g*') ; colormap gray ; +% +%% swap x with y +%ell = ell([2 1 5 4 3],:) ; +% +%for k=1:size(ell,2) +% E = ell(:,k) ; +% c = E(1:2) ; +% A = zeros(2) ; +% A(1,1) = E(3) ; +% A(1,2) = E(4) ; +% A(2,2) = E(5) ; +% A = A + A' - diag(diag(A)) ; +% +% [V,D] = eig(A) ; +% A = 2.5*V*sqrt(D) ; +% +% X = A*[cos(linspace(0,2*pi,30)) ; sin(linspace(0,2*pi,30)) ;] ; +% X(1,:) = X(1,:) + c(1) ; +% X(2,:) = X(2,:) + c(2) ; +% +% plot(X(1,:),X(2,:),'r-','LineWidth',2) ; +% plot(c(1),c(2),'r.') ; +% plot(j(k),i(k),'g*') ; +%end +% +%line([j'; ell(1,:)],[i'; ell(2,:)],'color','b') ; diff --git a/SD-VBS/benchmarks/mser/src/matlab/overview_mser.m b/SD-VBS/benchmarks/mser/src/matlab/overview_mser.m new file mode 100755 index 0000000..b7fcf2b --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/overview_mser.m @@ -0,0 +1,8 @@ +% OVERVIEW_MSER Maximally Stable Extremal Regions +% This is a MATLAB/MEX implementation of Maximally Stable Extremal +% Regions (MSER). You can: +% +% * Use MSER() to extract the maximally stable extremal regions from +% a given image. +% +% For practical coding example, see MSER_DEMO() and MSER_DEMO3(). diff --git a/SD-VBS/benchmarks/mser/src/matlab/script_run_profile.m b/SD-VBS/benchmarks/mser/src/matlab/script_run_profile.m new file mode 100755 index 0000000..ec06539 --- /dev/null +++ b/SD-VBS/benchmarks/mser/src/matlab/script_run_profile.m @@ -0,0 +1,126 @@ +function script_run_profile(dataDir, resultDir, type, common,toolDir) + +mser_compile; +path(path, common); + +% MSER_DEMO Demonstrates MSER + +% 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. + +which_image = 3 ; + +% -------------------------------------------------------------------- +% Create data +% -------------------------------------------------------------------- +switch which_image + case 1 + I = rand(200,200) ; + I = imsmooth(I,10) ; + I = I-min(I(:)) ; + I = I/max(I(:)) ; + lev = 10 ; + I = uint8(round(I*lev)) ; + + case 2 + I = zeros(200,200) ; + I(50:150,50:150)=5 ; + I = imsmooth(I,10) ; + I = uint8(round(I)) ; + + case 3 + Files = dir([dataDir,'/1.bmp']); + I=readImage(fullfile(dataDir,Files(1).name)); +end + +% -------------------------------------------------------------------- +% Compute MSERs +% -------------------------------------------------------------------- + +%% Self check params +tol = 0.1; +elapsed = zeros(1,2); + +rows = size(I,1); +cols = size(I,2); + +fprintf(1,'Input size\t\t- (%dx%d)\n', rows, cols); + +%% Timing +start = photonStartTiming; + +[idx] = mser(uint8(I), 2) ; + +%% Timing +stop = photonEndTiming; + +temp = photonReportTiming(start, stop); +elapsed(1) = elapsed(1) + temp(1); +elapsed(2) = elapsed(2) + temp(2); + + %% Self checking + writeMatrix(idx, dataDir); + +%% Timing +photonPrintTiming(elapsed); + +%% -------------------------------------------------------------------- +%% Plots +%% -------------------------------------------------------------------- +%[i,j] = ind2sub(size(I),idx) ; +% +%figure(100) ; clf ; imagesc(I) ; hold on ; +%set(gca,'Position',[0 0 1 1]) ; +%plot(j,i,'g*') ; colormap gray ; +% +%% swap x with y +%ell = ell([2 1 5 4 3],:) ; +% +%for k=1:size(ell,2) +% E = ell(:,k) ; +% c = E(1:2) ; +% A = zeros(2) ; +% A(1,1) = E(3) ; +% A(1,2) = E(4) ; +% A(2,2) = E(5) ; +% A = A + A' - diag(diag(A)) ; +% +% [V,D] = eig(A) ; +% A = 2.5*V*sqrt(D) ; +% +% X = A*[cos(linspace(0,2*pi,30)) ; sin(linspace(0,2*pi,30)) ;] ; +% X(1,:) = X(1,:) + c(1) ; +% X(2,:) = X(2,:) + c(2) ; +% +% plot(X(1,:),X(2,:),'r-','LineWidth',2) ; +% plot(c(1),c(2),'r.') ; +% plot(j(k),i(k),'g*') ; +%end +% +%line([j'; ell(1,:)],[i'; ell(2,:)],'color','b') ; -- cgit v1.2.2