/**************************************************
Mex-C code for reproducing experiment 2: find and sketch over multiple scales
***************************************************/

# include <stdio.h>
# include <stdlib.h>
# include "mex.h"        /* the algorithm is connect to matlab */
# 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))


double *double_vector(int n)
{
    double *v; 
    
    v = (double*) mxCalloc (n, sizeof(double));
    
    return v; 
}

int *int_vector(int n)
{
    int *v; 
    
    v = (int*) mxCalloc (n, sizeof(int));
    
    return v; 
}

double **double_matrix(int m, int n)
{
    double **mat; 
    int i; 
    
    mat = (double**) mxCalloc(m, sizeof(double*)); 
    for (i=0; i<m; i++)
        mat[i] = double_vector(n); 
    
    return mat; 
}

int **int_matrix(int m, int n)
{
    int **mat; 
    int i; 
    
    mat = (int**) mxCalloc(m, sizeof(int*)); 
    for (i=0; i<m; i++)
        mat[i] = int_vector(n); 
    
    return mat; 
}

void free_matrix(void **a, int nrow, int ncol)
{
	int i;
	for(i = 0; i < nrow; i++)
		mxFree(a[i]);
	mxFree(a);
}


/* getting the index of matlab image, (x,y) location, (sx, sy) sizes */
int px(int x, int y, int bx, int by)    
{            
   return (x + (y-1)*bx - 1); 
 }

 /* variables */
 int N;                      /* number of orientations */
 double **mI;                /* local maximum pooling for filtered images */
 int sx, sy;                 /* current size of image */
 double *Mi, *Mx, *My, *Mm, *Mm1;  /* storing selected bases */
 int sub;                    /* subsampling */
 double *tss;                /* map of log-likelihood score */
 int Totalsketch;            /* total number of bases */

 
 
/* the shared sketch algorithm */
void Cmp()
{
   int x0, y0, x, y, x1, y1, mi, mx, my, t; 
   double resp; 
   
   t = 0; 
   do  
   {
     mi = floor(Mi[t]+.5); mx = floor(Mx[t]+.5); my = floor(My[t]+.5);  
     /* the following is a parallel search */
    for (x=1; x<=floor(sx/sub); x++)
         for (y=1; y<=floor(sy/sub); y++)
          {            
           x0 = x*sub; y0 = y*sub;   
           x1 = x0+mx; y1 = y0+my; 
           if ((x1>=1)&&(x1<=sx)&&(y1>=1)&&(y1<=sy))
               resp = mI[mi][px(x1, y1, sx, sy)];
           else
               resp = 0.; 
              
           tss[px(x0, y0, sx, sy)] +=  (Mm[t]*resp - Mm1[t]);  
          }   
     t++; 
     }
  while (t<Totalsketch);   /* stopping criterion */  
}


/* mex function is used to pass on the pointers and scalars from matlab, 
   so that heavy computation can be done by C, which puts the results into 
   some of the pointers. After that, matlab can then use these results. 
   
   So matlab is very much like a managing platform for organizing the 
   experiments, and mex C is like a work enginee for fast computation. */

void mexFunction(int nlhs, mxArray *plhs[], 
                 int nrhs, const mxArray *prhs[])                
{
 int ind, c; 
 mxArray *f;  
 
 c = 0; /* counter for input pointers and scalars */
 N = floor(mxGetScalar(prhs[c++])+.5);  /* number of orientations */

 mI = mxCalloc(N, sizeof(double*));   /* fitered images */
 for (ind=0; ind<N; ind++)
      {  
       f = mxGetCell(prhs[c], ind); 
       mI[ind] = mxGetPr(f);    /* get pointers to filtered images */       
      }
 c++;
 
 Mi = mxGetPr(prhs[c++]);   /* orientation of selected Gabor */               
 Mx = mxGetPr(prhs[c++]);         
 My = mxGetPr(prhs[c++]);   /* (x, y) position of selected Gabor */
 Mm = mxGetPr(prhs[c++]);   /* lambda */ 
 Mm1 = mxGetPr(prhs[c++]);  /* logZ */
 sub = floor(mxGetScalar(prhs[c++])+.5);   /* sub-sampling of Gabor filters to improve speed */
 Totalsketch = floor(mxGetScalar(prhs[c++])+.5);  /* total number of bases */
 tss = mxGetPr(prhs[c++]); 
 sx = floor(mxGetScalar(prhs[c++])+.5); 
 sy = floor(mxGetScalar(prhs[c++])+.5); 
 
 Cmp();     
}

     



 

                    