diff --git a/generic/imgraph.c b/generic/imgraph.c index 36d6aed..342feae 100644 --- a/generic/imgraph.c +++ b/generic/imgraph.c @@ -665,6 +665,9 @@ static int imgraph_(mergetree)(lua_State *L) { t->cs = cs; t->rs = rs; t->altitudes = altitudes; + + + return 1; } @@ -714,6 +717,21 @@ static int imgraph_(filtertree)(lua_State *L) { return 1; } + +static int imgraph_(cuttree)(lua_State *L) { + // get args + MergeTree *t = lua_toMergeTree(L, 1); + int table_weights = 2; // arg 2 + + + //calling the labeling method on the merge tree + list * cut; + cut = MSF_Kruskal(t); + + // done + return 1; +} + static int imgraph_(weighttree)(lua_State *L) { // get args MergeTree *t = lua_toMergeTree(L, 1); @@ -1376,6 +1394,7 @@ static const struct luaL_Reg imgraph_(methods__) [] = { {"tree2components", imgraph_(tree2components)}, {"adjacency", imgraph_(adjacency)}, {"segm2components", imgraph_(segm2components)}, + {"cuttree", imgraph_(cuttree)}, {NULL, NULL} }; diff --git a/init.lua b/init.lua index cc65a11..ee08570 100644 --- a/init.lua +++ b/init.lua @@ -336,6 +336,28 @@ function imgraph.weighttree(...) torch.Tensor().imgraph.weighttree(tree, weights) end +---------------------------------------------------------------------- +-- computes a cut in a tree +-- +function imgraph.cuttree(...) + --get args + local args = {...} + local tree = args[1] + local weights = args[2] + + -- usage + if not tree or not weights then + print(xlua.usage('imgraph.cuttree', + 'computes a cut in a tree', nil, + {type='imgraph.MergeTree', help='merge tree to be weighted', req=true}, + {type='table', help='a list of weights t[k] = w, with k the index of the node to be weighted', req=true})) + xlua.error('incorrect arguments', 'imgraph.cuttree') + end + + -- compute cut + torch.Tensor().imgraph.cuttree(tree, weights) +end + ---------------------------------------------------------------------- -- transform a merge tree back into a graph, for visualization -- @@ -752,14 +774,14 @@ imgraph._example = [[ local mt = imgraph.mergetree(graph) -- (7) display results - image.display{image=inputimg, legend='input image'} - image.display{image=cc, legend='thresholded graph'} - image.display{image=watershed, legend='watershed on the graph'} - image.display{image=watershedcc, legend='components of watershed'} - image.display{image=mstsegmcolor, legend='segmented graph, using min-spanning tree'} - image.display{image=pool, legend='original imaged hist-pooled by segmentation'} - image.display{image=hierarchy, legend='raw edge-weighted graph watershed'} - image.display{image=filteredhierarchy, legend='filtered edge-weighted graph watershed'} + image.display{image=inputimg, legend='input image'} + image.display{image=cc, legend='thresholded graph'} + image.display{image=watershed, legend='watershed on the graph'} + image.display{image=watershedcc, legend='components of watershed'} + image.display{image=mstsegmcolor, legend='segmented graph, using min-spanning tree'} + image.display{image=pool, legend='original imaged hist-pooled by segmentation'} + image.display{image=hierarchy, legend='raw edge-weighted graph watershed'} + image.display{image=filteredhierarchy, legend='filtered edge-weighted graph watershed'} ]] function imgraph.testme(usrimg) local inputimg diff --git a/mergetree.h b/mergetree.h index bcc2e4f..d3e9248 100644 --- a/mergetree.h +++ b/mergetree.h @@ -14,6 +14,13 @@ typedef struct { int rs; } MergeTree; + typedef struct list + { + int index; + struct list *next; + } list ; + + static MergeTree *lua_toMergeTree (lua_State *L, int index) { MergeTree *mt = (MergeTree *)lua_touserdata(L, index); diff --git a/pink/CMakeLists.txt b/pink/CMakeLists.txt index f35f55e..352e122 100644 --- a/pink/CMakeLists.txt +++ b/pink/CMakeLists.txt @@ -23,7 +23,9 @@ add_library(pink mcindic.c mclifo.c mcfifo.c - mcunionfind.c) + mcunionfind.c + mtrand64.c + kruskal.c) set (CMAKE_C_FLAGS "-fpic -g -DUNIXIO") target_link_libraries(pink) diff --git a/pink/kruskal.c b/pink/kruskal.c new file mode 100644 index 0000000..39115a0 --- /dev/null +++ b/pink/kruskal.c @@ -0,0 +1,420 @@ +/* +Kruskal algorithm for maximum spanning forest computation +author: Camille Couprie +21 oct. 2011 +*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define false 0 +#define true 1 + +typedef struct { + mtree *tree; + RAG *rag; + struct xvimage *labels; + int32_t *altitudes; + float *weights; + int cs; + int rs; +} MergeTree; + + typedef struct list + { + int index; + struct list *next; + } list ; + + +void Insert(list **sl, int index) + { + list *tmp = NULL; + list *csl = *sl; + list *elem = malloc(sizeof(list)); + if(!elem) exit(EXIT_FAILURE); + elem->index = index; + while(csl) + { + tmp = csl; + csl = csl->next; + } + elem->next = csl; + if(tmp) tmp->next = elem; + else *sl = elem; + } + +void PrintList(list *sl); +int element_link( int x,int y, uint32_t *Rnk, uint32_t *Fth); +int element_find(int x, uint32_t *Fth ); +void TriRapideStochastique_dec (float * A, uint32_t *I, long p, long r); +long PartitionStochastique_dec (float *A, uint32_t * I, long p, long r); +int nb_neighbors(int x, JCctree *CT, int nb_leafs); +int neighbor(int x, int k, JCctree *CT, int nb_leafs, int * SeededNodes); + +/*=====================================================================================*/ + list * MSF_Kruskal(MergeTree * MT) // max_weight /* maximum weight value */ ) +/*=====================================================================================*/ +/*Segment a tree into two components. + Returns a list of nodes correspunding to the Max Spanning Forest cut, + computed using Kruskal's algorithm */ +{ + int i, j, k, x, y, e1, e2; + int nb_markers; int nb_leafs; + long N, M; + + mtree * T= MT->tree; + float * W = MT->weights; + // mergeTreePrint(T); + JCctree *CT = T->CT; + int root_node = CT->root; + + //nb nodes + M = CT->nbnodes; + + // nb_edges + nb_leafs = 0; + for (i = 0; i < M; i++) + if (CT->tabnodes[i].nbsons == 0) + nb_leafs++; + + nb_markers = nb_leafs+1; + N=M+nb_markers; + M=N-1; + printf("Nb nodes:%d Nb edges: %d Nb leafs :%d \n", N, M, nb_leafs); + + // indexes of edges : son's nodes indexes + //Memory allocation of temporary arrays for Krukal's algorithm + Lifo * LIFO; + LIFO = CreeLifoVide(M); + if (LIFO == NULL) { fprintf(stderr, "kruskal : CreeLifoVide failed\n"); exit(0); } + + int * Mrk = (int*)calloc(N ,sizeof(int)); + if (Mrk == NULL) { fprintf(stderr, "kruskal : malloc failed\n"); exit(0); } + + uint32_t * SeededNodes = (uint32_t*)malloc(nb_markers*sizeof(uint32_t)); + if (SeededNodes == NULL) { fprintf(stderr, "kruskal : malloc failed\n"); exit(0); } + + // markers + + SeededNodes[0]= M; + j=1; + for (i = 0; i < CT->nbnodes; i++) + if (CT->tabnodes[i].nbsons == 0) + { + SeededNodes[j]= i+CT->nbnodes; + Mrk[SeededNodes[j]] = 1; + printf("seeded node %d \n", SeededNodes[j]); + j++; + } + Mrk[M] = 0; + + + uint32_t * Rnk = (uint32_t*)calloc(N, sizeof(uint32_t)); + if (Rnk == NULL) { fprintf(stderr, "kruskal : malloc failed\n"); exit(0); } + uint32_t * Fth = (uint32_t*)malloc(N*sizeof(uint32_t)); + if (Fth == NULL) { fprintf(stderr, "kruskal : malloc failed\n"); exit(0); } + for(k=0;knbnodes+k]=0.5; + fprintf(stderr,"%f\n", sorted_weights[CT->nbnodes+k]); + } + TriRapideStochastique_dec(sorted_weights,Es, 0, M-1); + free(sorted_weights); + + long nb_arete = 0; + int e_max, root; + long cpt_aretes = 0; + + // beginning of main loop + while (nb_arete < N-nb_markers+1) + { + e_max=Es[cpt_aretes]; + cpt_aretes=cpt_aretes+1; + e1= e_max; // e1 = Edges[0][e_max]; + if (e_maxnbnodes) e2= CT->tabnodes[e_max].father; + else e2= e_max-CT->nbnodes; + if (e2==-1)e2=M; //e2 = Edges[1][e_max]; + x = element_find(e1, Fth ); + y = element_find(e2, Fth ); + if ((x != y) && (!(Mrk[x]>=1 && Mrk[y]>=1))) + { + root = element_link( x,y, Rnk, Fth); + nb_arete=nb_arete+1; + if ( Mrk[x]>=1) Mrk[root]= Mrk[x]; + else if ( Mrk[y]>=1) Mrk[root]= Mrk[y]; + } + } + + //building the labeling for each individual markers in map + // (find the root vertex of each tree) + int * Map2 = (int *)malloc(N*sizeof(int)); + int * Map = (int *)malloc(N*sizeof(int)); + for (i=0; inbnodes; i++) + { + // nodes having a different value than their father are in the cut + if( Map[CT->tabnodes[i].father] != Map[i]) + Insert(&cut, i); + // leafs having the same label as the root are in the cut + if ((CT->tabnodes[i].nbsons == 0) && (Map[i]==0)) + Insert(&cut, i); + } + + PrintList(cut); + + + LifoTermine(LIFO); + free(Mrk); + free(Rnk); + free(Es); + free(Map); + free(Map2); + free(Fth); + return cut; +} +/*================================================*/ +void PrintList(list *sl) +/*================================================*/ +{ +while(sl) + { + printf("%d\n",sl->index); + sl = sl->next; + } +} + +/*================================================*/ +int nb_neighbors(int x, JCctree *CT, int nb_leafs) +/*================================================*/ +{ + int tmp; + if (xnbnodes) + { + tmp = CT->tabnodes[x].nbsons; + if (tmp==0) tmp=1; + return tmp+1; + } + else if (xnbnodes+nb_leafs) + return 1; //CT->tabnodes[CT->root].nbsons; + else return 1; +} + + +/*================================================*/ +int neighbor(int x, int k, JCctree *CT, int nb_leafs, int * SeededNodes) +/*================================================*/ +{ + JCsoncell *s;int i, tmp; + if (xnbnodes) + { + tmp = CT->tabnodes[x].nbsons; + if (tmp==0) + { if (k==0) return CT->tabnodes[x].father; + return SeededNodes[x+1]; + } + else if (k<=tmp) + { + if (k==tmp) return CT->tabnodes[x].father; + s = CT->tabnodes[x].sonlist; + for (i=0;inext; // INEFFICACE A REFAIRE + return s->son; + } + } + else if (xnbnodes+nb_leafs) + return x-CT->nbnodes; + else + { + return CT->root; + } +} + +/*================================================*/ +int element_link( int x,int y, uint32_t *Rnk, uint32_t *Fth) +/*================================================*/ +{ + if( Rnk[x] > Rnk[y]) + { + int t; + t=x; + x=y; + y=t; + } + if( Rnk[x] == Rnk[y]) + { + Rnk[y]=Rnk[y]+1; + } + Fth[x] = y; + return y; +} + + + +/*===============================*/ +int element_find(int x, uint32_t *Fth ) +/*===============================*/ +{ + if (Fth[x] != x) + Fth[x] = element_find(Fth[x], Fth); + return Fth[x]; +} + + /* // Printing the merge tree + int32_t i; + JCsoncell *s; + JCctree *CT = MT->CT; + printf("root = %d ; nbnodes: %d ; nbsoncells: %d\n", CT->root, CT->nbnodes, CT->nbsoncells); + for (i = 0; i < CT->nbnodes; i++) + { + printf("node: %d ; level %d ; nbsons: %d ; father: %d ; mergeEdge %d", + i, CT->tabnodes[i].data, CT->tabnodes[i].nbsons, CT->tabnodes[i].father, MT->mergeEdge[i]); + if (CT->tabnodes[i].nbsons > 0) + { + printf("sons: "); + for (s = CT->tabnodes[i].sonlist; s != NULL; s = s->next) + printf("%d ", s->son); + } + printf("\n"); + } + */ + + + +/* =============================================================== */ +long Partitionner_dec(float *A, uint32_t * I, long p, long r) +/* =============================================================== */ +/* + partitionne les elements de A entre l'indice p (compris) et l'indice r (compris) + en deux groupes : ceux <= A[p] et les autres. +*/ +{ + float t; + int t1; + float x = A[p]; + long i = p - 1; + long j = r + 1; + while (1) + { + do j--; while (A[j] < x); + do i++; while (A[i] > x); + if (i < j) + { + t = A[i]; + A[i] = A[j]; + A[j] = t; + t1 = I[i]; + I[i] = I[j]; + I[j] = t1; + } + else return j; + } /* while (1) */ +} /* Partitionner() */ + +/* =============================================================== */ +long PartitionStochastique_dec (float *A, uint32_t * I, long p, long r) +/* =============================================================== */ +/* + partitionne les elements de A entre l'indice p (compris) et l'indice r (compris) + en deux groupes : ceux <= A[q] et les autres, avec q tire au hasard dans [p,r]. +*/ +{ + float t; + int t1; + long q; + + q = p + (genrand64_int64() % (r - p + 1)); /* rand must be 64-bit safe, should be OK now */ + // assert((q >= p) && (q <= r)) ; /* you'd be surprised */ + t = A[p]; /* echange A[p] et A[q] */ + A[p] = A[q]; + A[q] = t; + + t1 = I[p]; /* echange I[p] et I[q] */ + I[p] = I[q]; + I[q] = t1; + + return Partitionner_dec(A, I, p, r); +} /* PartitionStochastique() */ + + +/* =============================================================== */ +void TriRapideStochastique_dec (float * A, uint32_t *I, long p, long r) +/* =============================================================== */ +/* + trie les valeurs du tableau A de l'indice p (compris) a l'indice r (compris) + par ordre decroissant +*/ +{ + long q; + if (p < r) + { + q = PartitionStochastique_dec(A, I, p, r); + TriRapideStochastique_dec (A, I, p, q) ; + TriRapideStochastique_dec (A, I, q+1, r) ; + } +} /* TriRapideStochastique() */ diff --git a/pink/mtrand64.c b/pink/mtrand64.c new file mode 100644 index 0000000..e9bc516 --- /dev/null +++ b/pink/mtrand64.c @@ -0,0 +1,167 @@ +/* + A C-program for MT19937-64 (2004/9/29 version). + Coded by Takuji Nishimura and Makoto Matsumoto. + + This is a 64-bit version of Mersenne Twister pseudorandom number + generator. + + Before using, initialize the state by using init_genrand64(seed) + or init_by_array64(init_key, key_length). + + Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT OWNER OR + 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. + + References: + T. Nishimura, ``Tables of 64-bit Mersenne Twisters'' + ACM Transactions on Modeling and + Computer Simulation 10. (2000) 348--357. + M. Matsumoto and T. Nishimura, + ``Mersenne Twister: a 623-dimensionally equidistributed + uniform pseudorandom number generator'' + ACM Transactions on Modeling and + Computer Simulation 8. (Jan. 1998) 3--30. + + Any feedback is very welcome. + http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces) +*/ + + +#include +#include "mtrand64.h" + +#define NN 312 +#define MM 156 +#define MATRIX_A 0xB5026F5AA96619E9ULL +#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */ +#define LM 0x7FFFFFFFULL /* Least significant 31 bits */ + + +/* The array for the state vector */ +static unsigned long long mt[NN]; +/* mti==NN+1 means mt[NN] is not initialized */ +static int mti=NN+1; + +/* initializes mt[NN] with a seed */ +void init_genrand64(unsigned long long seed) +{ + mt[0] = seed; + for (mti=1; mti> 62)) + mti); +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +void init_by_array64(unsigned long long init_key[], + unsigned long long key_length) +{ + unsigned long long i, j, k; + init_genrand64(19650218ULL); + i=1; j=0; + k = (NN>key_length ? NN : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 3935559000370003845ULL)) + + init_key[j] + j; /* non linear */ + i++; j++; + if (i>=NN) { mt[0] = mt[NN-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=NN-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 2862933555777941757ULL)) + - i; /* non linear */ + i++; + if (i>=NN) { mt[0] = mt[NN-1]; i=1; } + } + + mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0, 2^64-1]-interval */ +unsigned long long genrand64_int64(void) +{ + int i; + unsigned long long x; + static unsigned long long mag01[2]={0ULL, MATRIX_A}; + + if (mti >= NN) { /* generate NN words at one time */ + + /* if init_genrand64() has not been called, */ + /* a default initial seed is used */ + if (mti == NN+1) + init_genrand64(5489ULL); + + for (i=0;i>1) ^ mag01[(int)(x&1ULL)]; + } + for (;i>1) ^ mag01[(int)(x&1ULL)]; + } + x = (mt[NN-1]&UM)|(mt[0]&LM); + mt[NN-1] = mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)]; + + mti = 0; + } + + x = mt[mti++]; + + x ^= (x >> 29) & 0x5555555555555555ULL; + x ^= (x << 17) & 0x71D67FFFEDA60000ULL; + x ^= (x << 37) & 0xFFF7EEE000000000ULL; + x ^= (x >> 43); + + return x; +} + +/* generates a random number on [0, 2^63-1]-interval */ +long long genrand64_int63(void) +{ + return (long long)(genrand64_int64() >> 1); +} + +/* generates a random number on [0,1]-real-interval */ +double genrand64_real1(void) +{ + return (genrand64_int64() >> 11) * (1.0/9007199254740991.0); +} + +/* generates a random number on [0,1)-real-interval */ +double genrand64_real2(void) +{ + return (genrand64_int64() >> 11) * (1.0/9007199254740992.0); +} + +/* generates a random number on (0,1)-real-interval */ +double genrand64_real3(void) +{ + return ((genrand64_int64() >> 12) + 0.5) * (1.0/4503599627370496.0); +} diff --git a/pink/mtrand64.h b/pink/mtrand64.h new file mode 100644 index 0000000..0a1873c --- /dev/null +++ b/pink/mtrand64.h @@ -0,0 +1,80 @@ +/* + A C-program for MT19937-64 (2004/9/29 version). + Coded by Takuji Nishimura and Makoto Matsumoto. + + This is a 64-bit version of Mersenne Twister pseudorandom number + generator. + + Before using, initialize the state by using init_genrand64(seed) + or init_by_array64(init_key, key_length). + + Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT OWNER OR + 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. + + References: + T. Nishimura, ``Tables of 64-bit Mersenne Twisters'' + ACM Transactions on Modeling and + Computer Simulation 10. (2000) 348--357. + M. Matsumoto and T. Nishimura, + ``Mersenne Twister: a 623-dimensionally equidistributed + uniform pseudorandom number generator'' + ACM Transactions on Modeling and + Computer Simulation 8. (Jan. 1998) 3--30. + + Any feedback is very welcome. + http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces) +*/ + + +/* initializes mt[NN] with a seed */ +void init_genrand64(unsigned long long seed); + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +void init_by_array64(unsigned long long init_key[], + unsigned long long key_length); + +/* generates a random number on [0, 2^64-1]-interval */ +unsigned long long genrand64_int64(void); + + +/* generates a random number on [0, 2^63-1]-interval */ +long long genrand64_int63(void); + +/* generates a random number on [0,1]-real-interval */ +double genrand64_real1(void); + +/* generates a random number on [0,1)-real-interval */ +double genrand64_real2(void); + +/* generates a random number on (0,1)-real-interval */ +double genrand64_real3(void);