do_sp2.c

Standard model object recognition matlab code

/**********************************************/
/* do_sp2.c                                   */
/* M. Kouh                                    */
/* August, 2005                               */
/**********************************************/

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

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
  double *x1, *x2, *w1, *w2, *w_idx, *mn_size, *shift;
  double *y1_out, *y2_out;
  double *x1_orig, *x2_orig;
  double y1, y2;
  int wm, wn, ym, shift_val;
  int xm_max, xn_max, wm_max, wn_max;
  int ym_max, yn_max;
  int i, k, m, n, idx;

  /* Check for proper number of arguments. */
  if (nrhs!=7) {
    mexErrMsgTxt("error: Nine inputs required.");
  } else if (nlhs!=2) {
    mexErrMsgTxt("error: Number of output arguments");
  }

  /* Assign pointers to each input. */
  x1      = mxGetPr(prhs[0]); /* input (image to be convolved) */
  x2      = mxGetPr(prhs[1]); /* input (image to be convolved) */
  w1      = mxGetPr(prhs[2]); /* weights (filters for convolution) */
  w2      = mxGetPr(prhs[3]); /* weights (filters for convolution) */
  w_idx   = mxGetPr(prhs[4]); /* index of nonzero weights */
  mn_size = mxGetPr(prhs[5]); /* size of image and filter */
  shift   = mxGetPr(prhs[6]); /* shift */

  /* get size of the input and calculate the output size */
  wm = mxGetM(prhs[2]); /* number of old features times x times y */
  wn = mxGetN(prhs[2]); /* number of new features */
  xm_max = mn_size[0];  /* x-dim size of input image */
  xn_max = mn_size[1];  /* y-dim size of input image */
  wm_max = mn_size[2];  /* x-dim size of filter */
  wn_max = mn_size[3];  /* y-dim size of filter */
  shift_val = *shift;

  /* The size of x has been either cropped or zero-padded to ensure
     the convolution would work.  ie. ym_max and yn_max must be
     integer values greater than (>) 0. */
  ym_max = floor((xm_max-wm_max)/shift_val)+1;
  yn_max = floor((xn_max-wn_max)/shift_val)+1;
  ym = ym_max*yn_max;

  /* if input is smaller than the filter, just return 0 */
  /* Better size check is done with crop_or_zeropad.m */
  if (ym_max<1 | yn_max<1) { 
    plhs[0] = mxCreateDoubleMatrix(1,wn,mxREAL);
    y1_out = mxGetPr(plhs[0]); /* output */
    y2_out = mxGetPr(plhs[1]); /* output */
  }

  else {
    plhs[0] = mxCreateDoubleMatrix(ym,wn,mxREAL);
    plhs[1] = mxCreateDoubleMatrix(ym,wn,mxREAL);
    y1_out = mxGetPr(plhs[0]); /* output */
    y2_out = mxGetPr(plhs[1]); /* output */

    x1_orig = x1;
    x2_orig = x2;

    /* cycle thru the features */
    /* w_idx has one dimension more than w. */
    for (i=0; i<wn; i++, w1+=wm, w2+=wm, w_idx+=wm+1) {

      /* move along the x(m) and y(n) dimension of the image */
      /* and shift by the prescribed amount. */
      for (n=0; n<yn_max; n++) {
	x1 = x1_orig+n*shift_val*xm_max;
	x2 = x2_orig+n*shift_val*xm_max;

	for (m=0; m<ym_max; m++, x1+=shift_val, x2+=shift_val) { 
	  y1=0;
	  y2=0;

	  /* go thru the filter, and stop the loop when idx<0 */
	  for (k=0, idx=*w_idx; idx>=0; k++, idx=*(w_idx+k)) { 
	    y1 += *(w1+k) * *(x1+idx);
	    y2 += *(w2+k) * *(x2+idx);
	  }
	  *(y1_out++) = y1;
	  *(y2_out++) = y2;
	}
      }
    }
  }
}