rct.zip_decomposition_modulate2( )_rct_zip

/* MEX file to rank relevant images using WEIGHTED probabilistic distance * * Created by: Minh Do, 7 July 1999 * * To compile: * mex -O -v rankwpd.c * * Note: This is different with rankpd.c in that it takes into the account * on the entropy weights from the query image. */ #include "mex.h" #include <stdlib.h> #include <math.h> /* From Numerical Recipes */ double gammaln(double xx) { double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j<=5;j++) ser += cof[j]/++y; return -tmp+log(2.5066282746310005*ser/x); } /* Symmetric Kullback-Leibler distance between * two generalized Gaussian pdf's */ double kldistS(double a1, double b1, double a2, double b2) { double res; res = pow(a1/a2, b2)*exp(gammaln((b2+1.0)/b1)-gammaln(1.0/b1)) - 1.0/b1 + pow(a2/a1, b1)*exp(gammaln((b1+1.0)/b2)-gammaln(1.0/b2)) - 1.0/b2; return ((res > 0.0) ? res : 0.0); } /* Asymmetric Kullback-Leibler distance between * two generalized Gaussian pdf's * (a1, b1) is the query parameters * (a2, b2) is the candidate image parameters * Result is: * KLD((a1,b1) || (a2,b2)) or KLD(query || image) */ double kldistA(double a1, double b1, double a2, double b2) { double res; res = log((b1 * a2) / (b2 * a1)) + gammaln(1.0/b2) - gammaln(1.0/b1) + pow(a1/a2, b2)*exp(gammaln((b2+1.0)/b1)-gammaln(1.0/b1)) - 1.0/b1; return ((res > 0.0) ? res : 0.0); } /* Kullback-Leibler distance between two feature vectors * fv1 is supposed to be from the query image * fv2 is supposed to be from the candidate image * ws is the weigths (computed from the query image) */ double probdist(double *fv1, double *fv2, double *w, int nb) { /* Assumption on features: a, b, a, b, (next level), ... for total nb pairs corresponding to nb subbands */ int i; double dist = 0.0; for (i = 0; i < nb; i++) dist += w[i] * kldistA(fv1[2*i], fv1[2*i+1], fv2[2*i], fv2[2*i+1]); return dist; } /* function rr = rankwpd(fs, ws, ns) % RANKWPD Rank relevant images using WEIGHTED probabilistic distance % % Input: % fs: feature sets, one column per image % ws: entropy weights, one column per image % ns: (optional), number of subimages per texture class % (default 16) % % Output: % rr: rank of relevant images, one column for each query % % See also: RANKED * **/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *fs; double *ws; double *rr; int nimages, ns, nb, nf; int q, i, c, r; double *dist, dr; double *pr; fs = mxGetPr(prhs[0]); nf = mxGetM(prhs[0]); nb = nf / 2; /* 2 features per band */ nimages = mxGetN(prhs[0]); ws = mxGetPr(prhs[1]); if ((mxGetM(prhs[1]) != nb) || (mxGetN(prhs[1]) != nimages)) mexErrMsgTxt("Matrix dimensions mismatched"); if (nrhs < 3) ns = 16; else ns = (int) *(mxGetPr(prhs[2])); /* Create/allocate output */ plhs[0] = mxCreateDoubleMatrix(ns, nimages, mxREAL); rr = mxGetPr(plhs[0]); /* Array to keep distances */ dist = (double *) malloc(nimages * sizeof(double)); /* Simulate Query by Example process */ for (q = 0; q < nimages; q++) /* each query */ { for (i = 0; i < nimages; i++) /* each image */ dist[i] = probdist(fs + q*nf, fs + i*nf, ws + q*nb, nb); c = q / ns; /* class index of the query */ /* Find the rank of the relevant images */ for (r = 0; r < ns; r++) { pr = rr + q * ns + r; /* pointer to the corresponding rank */ dr = dist[c * ns + r]; /* distance of the relevant image */ *pr = 1; /* restart */ for (i = 0; i < nimages; i++) if (dist[i] < dr) (*pr)++; } } free(dist); }