/* TEMPLATE MATCHING INFERENCE */
# include <stdio.h>
# include <stdlib.h>
# include "mex.h"        
# include "math.h"
# define ROUND(x) (floor((x)+.5))
/* Compute pixel index in the vector that stores image */
int px(int x, int y, int lengthx, int lengthy)  /* the image is lengthx*lengthy */
{            
   return (x + (y-1)*lengthx - 1); 
 }
/* key input and output variables */
int numOrient, numImage, numElement; /* number of images and number of Gabors */   
float **MAX1map; 
int sizex, sizey, subsample, sizexSubsample, sizeySubsample; /* MAX1 maps are smaller than SUM1 maps by subsample */ 
double *selectedOrient, *selectedx, *selectedy, *selectedlambda, *selectedLogZ; /* parameters of learned active basis */  
double *SUM2score; /* template matching scores */ 
/* template matching */
void ComputeSUM2score()
{
   int img, t, here, besto, bestx, besty; 
   double bestlambda, bestLogZ; 
   
   for (img=0; img<numImage; img++)
          SUM2score[img] = 0.; 
   t = 0; 
   do  
   {   
     besto = ROUND(selectedOrient[t]); bestx = ROUND(selectedx[t]); besty = ROUND(selectedy[t]);  
     bestlambda = selectedlambda[t]; bestLogZ = selectedLogZ[t];   
     here = px(bestx , besty, sizexSubsample, sizeySubsample); 
     for (img=0; img<numImage; img++)
          SUM2score[img] += (bestlambda*MAX1map[besto*numImage+img][here] - bestLogZ); 
      t++; 
   }
  while (t<numElement);  
}
/* read in input variables and run the algorithm */
void mexFunction(int nlhs, mxArray *plhs[], 
                 int nrhs, const mxArray *prhs[])                
{
 int orient, img, i, c; 
 mxArray *f;  
 
 c = 0; 
 /* about active basis */
 numOrient = ROUND(mxGetScalar(prhs[c++]));  
 subsample = ROUND(mxGetScalar(prhs[c++]));   
 numElement = ROUND(mxGetScalar(prhs[c++])); 
 /* about input images */
 numImage = ROUND(mxGetScalar(prhs[c++])); 
 sizex = ROUND(mxGetScalar(prhs[c++])); 
 sizey = ROUND(mxGetScalar(prhs[c++]));   
 MAX1map = mxCalloc(numImage*numOrient, sizeof(float*));  
 for (img=0; img<numImage; img++)
     for (orient=0; orient<numOrient; orient++)
      {  
       i = orient*numImage+img; 
       f = mxGetCell(prhs[c], i); 
       MAX1map[i] = mxGetPr(f);    
      }
 c++; 
 /* learned parameters of active basis */
 selectedOrient = mxGetPr(prhs[c++]);                 
 selectedx = mxGetPr(prhs[c++]);         
 selectedy = mxGetPr(prhs[c++]);   
 selectedlambda = mxGetPr(prhs[c++]);       
 selectedLogZ = mxGetPr(prhs[c++]);   
 SUM2score = mxGetPr(prhs[c++]); 
 /* MAX1 maps are smaller than SUM1 maps */
 sizexSubsample = floor((double)sizex/subsample+.1); 
 sizeySubsample = floor((double)sizey/subsample+.1);
 /* run the shared sketch algorithm */
 ComputeSUM2score(); 
}



