hist_isect_c.c

来自「论文Beyond Bags of Features源代码」· C语言 代码 · 共 98 行

/******************************************************************************

Compute matrix of histogram intersection kernel values between two sets of vectors.
If one argument is a matrix of size m x o (where m is the number of vectors and
o is the dimensionality of each vector) and the second one is of size n x o,
the output is a matrix of size m x n.

From MATLAB, compile this mex function with the following command: 
mex hist_isect_c.c -lm 

Adapted from the svm_v0.55 toolbox: http://theoval.sys.uea.ac.uk/~gcc/svm/toolbox

******************************************************************************/

#include <math.h>
#include "mex.h"

#define min(a,b) (((a)<(b))?(a):(b))

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
   double *x1, *x2, *y;

   int row, i, j, k, m, n, o, N;

   /* check number of input and output arguments */

   if (nrhs != 2)
   {
      mexErrMsgTxt("Wrong number of input arguments.");
   }
   else if (nlhs > 1)
   {
      mexErrMsgTxt("Too many output arguments.");
   }

   /* get input arguments */

   if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]))
   {
      mexErrMsgTxt("x1 must be a double matrix.");
   }

   m  = mxGetM(prhs[0]);
   x1 = mxGetPr(prhs[0]);

   if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]))
   {
      mexErrMsgTxt("x2 must be a double matrix.");
   }

   n  = mxGetM(prhs[1]);
   o  = mxGetN(prhs[1]);
   x2 = mxGetPr(prhs[1]);

   /* allocate and initialise output matrix */

   plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);

   y = mxGetPr(plhs[0]);

   /* compute kernel matrix */

if (n > m) 
{
   for (i = 0; i < m; i++)
   {
      for (k = 0; k < o; k++)
      {
	 if (x1[i+k*m]==0) continue;

         for (j = 0; j < n; j++)
         {
            y[i+j*m] += min(x1[i+k*m], x2[j+k*n]);
         }
      }
   }
}
else
{
  for (j = 0; j < n; j++)
    {
      for (k = 0; k < o; k++)
	{
	  if (x2[j+k*n]==0) continue;

	  for (i = 0; i < m; i++)
	    {
	      y[i+j*m] += min(x1[i+k*m], x2[j+k*n]);
	    }
	}
    }
}

}