summaryrefslogtreecommitdiffstats
path: root/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c')
-rwxr-xr-xSD-VBS/common/toolbox/MultiNcut/multiAffinityic.c216
1 files changed, 216 insertions, 0 deletions
diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c
new file mode 100755
index 0000000..6919db9
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c
@@ -0,0 +1,216 @@
1/*================================================================
2* function w = affinityic(emag,ephase,pi,pj,subgrid,nrSubgrid,ncSubgrid,subpi,sigma)
3* Input:
4* emag = edge strength at each pixel
5* ephase = edge phase at each pixel
6* [pi,pj] = index pair representation for MALTAB sparse matrices, pi index in original grid, size(pj)=nrW*ncW+1
7* subgrid= index of the subgrid in the original image (first pixel's index is one!)
8* [nrSubgrid,ncSubgrid]= number of rows et colums of the subgrid
9* subpi = pi but with index in subgrid
10* sigma = sigma for IC energy
11* Output:
12* w = affinity with IC at [pi,pj]
13*
14
15% test sequence
16f = synimg(10);
17[i,j] = cimgnbmap(size(f),2);
18[ex,ey,egx,egy] = quadedgep(f);
19a = affinityic(ex,ey,egx,egy,i,j)
20show_dist_w(f,a);
21
22* Stella X. Yu, Nov 19, 2001.
23*=================================================================*/
24
25# include "mex.h"
26# include "math.h"
27
28void mexFunction(
29 int nargout,
30 mxArray *out[],
31 int nargin,
32 const mxArray *in[]
33)
34{
35 /* declare variables */
36 int nr, nc, np, nW, total;
37 int i, j, k, t, ix, iy, jx, jy, ii, jj, iip1, jjp1, iip2, jjp2, step,nrSubgrid, ncSubgrid;
38 double sigma, di, dj, a, z, maxori, phase1, phase2, slope;
39 int *ir, *jc;
40 /* unsigned long *pi, *pj, *subpi; */
41 unsigned int *pi, *pj, *subpi;
42 double *emag, *ephase, *w,*tmp,*subgrid;
43
44
45 /* check argument */
46 if (nargin<8) {
47 mexErrMsgTxt("Eight input arguments required");
48 }
49 if (nargout>1) {
50 mexErrMsgTxt("Too many output arguments");
51 }
52
53 /* get edgel information */
54 nr = mxGetM(in[0]);
55 nc = mxGetN(in[0]);
56 if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) {
57 mexErrMsgTxt("Edge magnitude and phase shall be of the same image size");
58 }
59 emag = mxGetPr(in[0]);
60 ephase = mxGetPr(in[1]);
61 np = nr * nc;
62
63 /*get subgrid information*/
64
65 tmp = mxGetData(in[5]);
66 nrSubgrid = (int)tmp[0];
67 tmp = mxGetData(in[6]);
68 ncSubgrid = (int)tmp[0];
69
70 /* printf("nrSubgrid=%d\n",nrSubgrid);
71 printf("ncSubgrid=%d\n",ncSubgrid); */
72 if (nrSubgrid* ncSubgrid != mxGetM(in[4])*mxGetN(in[4])) {
73 mexErrMsgTxt("Error in the size of the subgrid");
74 }
75 subgrid = mxGetData(in[4]);
76 nW = nrSubgrid * ncSubgrid;
77
78 /* get new index pair */
79 if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) {
80 mexErrMsgTxt("Index pair shall be of type UINT32");
81 }
82 if (mxGetM(in[3]) * mxGetN(in[3]) != nW + 1) {
83 mexErrMsgTxt("Wrong index representation");
84 }
85 pi = mxGetData(in[2]);
86 pj = mxGetData(in[3]);
87 subpi = mxGetData(in[7]);
88 /*{printf("pi[50] = %d\n",pi[50]);}
89 {printf("subpi[5] = %d\n",subpi[5]);}
90 {printf("subpi[6] = %d\n",subpi[6]);}
91 {printf("subpi[4] = %d\n",subpi[4]);}*/
92
93 /* create output */ /*!!!!!!!!!!!!!!!!!!!!!!!changer taille output!!!!!!!!!!*/
94 out[0] = mxCreateSparse(nW,nW,pj[nW],mxREAL);
95 if (out[0]==NULL) {
96 mexErrMsgTxt("Not enough memory for the output matrix");
97 }
98 w = mxGetPr(out[0]);
99 ir = mxGetIr(out[0]);
100 jc = mxGetJc(out[0]);
101
102 /* find my sigma */
103 if (nargin<9) {
104 sigma = 0;
105 for (k=0; k<np; k++) {
106 if (emag[k]>sigma) { sigma = emag[k]; }
107 }
108 sigma = sigma / 6;
109 printf("sigma = %6.5f",sigma);
110 } else {
111 sigma = mxGetScalar(in[8]);
112 }
113 a = 0.5 / (sigma * sigma);
114
115 /* computation */
116 total = 0;
117
118 for (j=0; j<nW; j++){
119 t= (int)subgrid[j]-1; /*on parcourt tous les pixels de la sous-grille dans la grille d'origine*/
120
121 jc[j] = total; /* total represente le nombre voisins du pixel j*/
122 jx = t / nr; /* col */
123 jy = t % nr; /* row */
124
125 for (k=pj[j]; k<pj[j+1]; k++) { /*k represente les indices correspondant au pixel j dans pi*/
126
127 i = pi[k]-1; /*i est un voisin de j a considerer*/
128
129 if (i==j) {
130 maxori = 1;
131
132 } else {
133
134 ix = i / nr;
135 iy = i % nr;
136
137 /* scan */
138 di = (double) (iy - jy);
139 dj = (double) (ix - jx);
140
141 maxori = 0.;
142 phase1 = ephase[j];
143
144
145 /* sample in i direction */
146 if (abs(di) >= abs(dj)) {
147 slope = dj / di;
148 step = (iy>=jy) ? 1 : -1;
149
150 iip1 = jy;
151 jjp1 = jx;
152
153
154 for (ii=0;ii<abs(di);ii++){
155 iip2 = iip1 + step;
156 jjp2 = (int)(0.5 + slope*(iip2-jy) + jx);
157
158 phase2 = ephase[iip2+jjp2*nr];
159
160 if (phase1 != phase2) {
161 z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
162 if (z > maxori){
163 maxori = z;
164 }
165 }
166
167 iip1 = iip2;
168 jjp1 = jjp2;
169 phase1 = phase2;
170 }
171
172 /* sample in j direction */
173 } else {
174 slope = di / dj;
175 step = (ix>=jx) ? 1: -1;
176
177 jjp1 = jx;
178 iip1 = jy;
179
180
181 for (jj=0;jj<abs(dj);jj++){
182 jjp2 = jjp1 + step;
183 iip2 = (int)(0.5+ slope*(jjp2-jx) + jy);
184
185 phase2 = ephase[iip2+jjp2*nr];
186
187 if (phase1 != phase2){
188 z = (emag[iip1+jjp1*nr] + emag[iip2+jjp2*nr]);
189 if (z > maxori){