/* 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);
}
- 1
- 2
- 3
前往页