/* SHARED SKETCH ALGORITHM FOR LEARNING ACTIVE BASIS */
# include <stdio.h>
# include <stdlib.h>
# include "mex.h"        
# include "math.h"
# define PI 3.1415926
# define ABS(x) ((x)>0? (x):(-(x)))
# define MAX(x, y) ((x)>(y)? (x):(y))
# define MIN(x, y) ((x)<(y)? (x):(y))
# 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); 
 }
/* key input and output variables */
int numOrient, numElement; /* number of orientations and number of Gabors */   
double **allSymbol; /* symbols of Gabors */              
int halfFilterSize; /* filter size = 2*halfFilterSize + 1 */ 
int sizex, sizey, subsample; /* size of template and subsample rate */ 
double *selectedOrient, *selectedx, *selectedy, *selectedPolar; /* parameters of learned active basis */  
float *commonTemplate; /* templates of active basis */    

/* plot the bar for Gabor at (mx, my, mo) */
void DrawElement(float *Template, int mo, int mx, int my, double w) 
{
  int x, y, here; 
  float a; 
          
  for (x=mx-halfFilterSize; x<=mx+halfFilterSize; x++)
     for (y=my-halfFilterSize; y<=my+halfFilterSize; y++)
      if ((x>=1)&&(x<=sizex)&&(y>=1)&&(y<=sizey)) 
        {
         a = allSymbol[mo][px(x-mx+halfFilterSize+1, y-my+halfFilterSize+1, 
                              2*halfFilterSize+1, 2*halfFilterSize+1)]*w; 
         here = px(x, y, sizex, sizey); 
         if (Template[here]<a)
             Template[here] = a;
        }
}

/* draw the common template */
void DrawTemplate()
{
   int t, x, y; 
   double w; 

   for (x=1; x<=sizex; x++)
       for (y=1; y<=sizey; y++)
           commonTemplate[px(x, y, sizex, sizey)] = 0.; 
   t = 0; 
   do  
   { 
     if (selectedPolar[t]>0.)
         w = 1.; 
     else 
         w = .1; 
     DrawElement(commonTemplate, ROUND(selectedOrient[t]), ROUND(selectedx[t])*subsample, ROUND(selectedy[t])*subsample, w);     
     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, c; 
 mxArray *f;  
 
 c = 0; 
 numOrient = ROUND(mxGetScalar(prhs[c++]));  
 numElement = ROUND(mxGetScalar(prhs[c++])); 
 sizex = ROUND(mxGetScalar(prhs[c++])); 
 sizey = ROUND(mxGetScalar(prhs[c++]));  
 subsample = ROUND(mxGetScalar(prhs[c++]));  
 halfFilterSize = ROUND(mxGetScalar(prhs[c++]));     
 allSymbol = mxCalloc(numOrient, sizeof(double*));    
 for (orient=0; orient<numOrient; orient++)
     {  
       f = mxGetCell(prhs[c], orient); 
       allSymbol[orient] = mxGetPr(f);       
     }
 c++; 
 selectedOrient = mxGetPr(prhs[c++]);                 
 selectedx = mxGetPr(prhs[c++]);         
 selectedy = mxGetPr(prhs[c++]);         
 selectedPolar = mxGetPr(prhs[c++]);  
 commonTemplate = mxGetPr(prhs[c++]);     

 DrawTemplate(); 
}



