/* SUM-MAX MAPS FOR TEMPLATE MATCHING INFERENCE */
# include <stdio.h>
# include <stdlib.h>
# include "mex.h"        
# include "math.h"
# define PI 3.1415926
# define ROUND(x) (floor((x)+.5))
# define NEGMAX -1e10
/* 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); 
 }
 /* variables */
int numResolution; /* number of resolutions and number of Gabors */   
double *allSizex, *allSizey;  /* sizes of images at multiple resolutions */
float **SUM2map, *MAX2score;    
int sizex, sizey, subsample, sizexSubsample, sizeySubsample; 
int sizeTemplatex, sizeTemplatey, sizeTemplatexSubsample, sizeTemplateySubsample, sizexEnlarge, sizeyEnlarge; 
double *allFx, *allFy; /* detected location over multiple resolutions */
/* compute the MAX2 score for each resolution */
void ComputeMAX2score(int resolution)
{
   int x, y, Fx, Fy;  
   float F, r, maxResponse; 

   F = NEGMAX; 
   for (x=-sizeTemplatexSubsample+1; x<=sizexSubsample; x++)
         for (y=-sizeTemplateySubsample+1; y<=sizeySubsample; y++) 
          {       
             r = SUM2map[resolution][px(x+sizeTemplatexSubsample, y+sizeTemplateySubsample, sizexEnlarge, sizeyEnlarge)]; 
             if (F<r) 
             {
                 F = r; Fx = x; Fy = y; 
             }
          }  
   MAX2score[resolution] = F; allFx[resolution] = Fx; allFy[resolution] = Fy; 
} 
           
/* read in the input and output variables */
void mexFunction(int nlhs, mxArray *plhs[], 
                 int nrhs, const mxArray *prhs[])                
{
 int resolution, orient, c; 
 mxArray *f;  
 
 c = 0; 
 numResolution = ROUND(mxGetScalar(prhs[c++]));
 allSizex = mxGetPr(prhs[c++]);
 allSizey = mxGetPr(prhs[c++]);
 subsample = ROUND(mxGetScalar(prhs[c++]));
 sizeTemplatex = ROUND(mxGetScalar(prhs[c++]));  
 sizeTemplatey = ROUND(mxGetScalar(prhs[c++])); 
 SUM2map = mxCalloc(numResolution, sizeof(float*)); 
 for (resolution=0; resolution<numResolution; resolution++)
   {
     f = mxGetCell(prhs[c], resolution); 
     SUM2map[resolution] = mxGetPr(f);  
  } 
 c++; 
 MAX2score = mxGetPr(prhs[c++]);
 allFx =  mxGetPr(prhs[c++]);
 allFy =  mxGetPr(prhs[c++]);
 
 sizeTemplatexSubsample = floor((double)sizeTemplatex/subsample); 
 sizeTemplateySubsample = floor((double)sizeTemplatey/subsample); 
 for (resolution=0; resolution<numResolution; resolution++)
 {
 sizex = ROUND(allSizex[resolution]); 
 sizey = ROUND(allSizey[resolution]); 
 sizexSubsample = floor((double)sizex/subsample); 
 sizeySubsample = floor((double)sizey/subsample); 
 sizexEnlarge = floor((double)(sizeTemplatex+sizex)/subsample); 
 sizeyEnlarge = floor((double)(sizeTemplatey+sizey)/subsample);   
 ComputeMAX2score(resolution);  
 }
}


     



 

                    