/**************************************************
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))
# define NEGMAX -1e10

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 images */
 int N;                      /* number of orientations */
 double **fI;                /* filtered images */
 double **mI;                /* local maximum pooling for filtered images */
 double **C;                 /* inhibition coefficients */
 int h;                      /* halfsizes of filters */
 int sx, sy;                 /* current size of image */
 int ax, ay;                 /* half size of bounding box of learned model */ 
 double T;                   /* threshold of selected bases */
 double *Mi, *Mx, *My, *Mm;  /* storing selected bases */
 double **allsymbol;         /* symbol of filters */
 int L, ore;                 /* allowed ranges of shifting in location and orientation */
 int sub;                    /* subsampling */
 double **tss;                /* map of log-likelihood score */
 double Upperbound;          /* satuation level for response */
 int Totalsketch;            /* total number of bases */
 double **Asym;               /* symbolic sketch */
 double *Sx1, *Sy1;           /* sizes of images */
 double *mscale;              /* log-likelihood scores over scales */
 int LIK;   /* If LIK = 1, we use log-likelihood score. If LIK = 0, we use active correlation */
 int **sinsh, **cossh;       /* store the shift to avoid repeated sin and cosin computation */
 
 
 /* store the shift values so that we do not need to repeat the sin and cos computation */ 
void storeshift()
{
    int ind, l;
    double theta; 
    
    sinsh = int_matrix(N, L+L+1); 
    cossh = int_matrix(N, L+L+1); 
    for (ind=0; ind<N; ind++)        
    {
        theta = PI*ind/N; 
        for (l=-L; l<=L; l++)
         {
            sinsh[ind][l+L] = floor(l*sub*sin(theta)+.5); 
            cossh[ind][l+L] = floor(l*sub*cos(theta)+.5); 
        }   
    }
}

/* for Gabor(x, y, orientation = ind), find local maximum in image i */
double shiftmax(int i, int ind, int x, int y, int *rm) 
{
   double m;
   int x1, y1, l, ml, mx, my, de, d, d1, here, mo, md; 

   m = NEGMAX; 
   for (l=-L; l<=L; l++)   
     {
        x1 = x + cossh[ind][l+L]; 
        y1 = y + sinsh[ind][l+L];   /* shifting in normal direction */
        if ((x1>=1)&&(x1<=sx)&&(y1>=1)&&(y1<=sy))
        {
        here = px(x1, y1, sx, sy);
        for (de=-ore; de<=ore; de++)
            {
              d = de+ind;
              d1 = d;
              if (d<0)
                 d1 = d+N;
              if (d>=N)
                  d1 = d-N;   /* shifting in orientation */
              if (fI[d1*n+i][here]>m)
                {
                   m = fI[d1*n+i][here];   /* local maximization */
                   ml = l; mx = x1; my = y1; mo = d1; md = de; 
                }
             }
        }
     }
   rm[0] = ml; rm[1] = mx; rm[2] = my; rm[3] = mo; rm[4] = md; 
   return(MIN(Upperbound, m)); 
}

/* symbolc sketch */
void draw(double *sym, int x0, int y0, int ind, double w)
{
  int x, y; 
  double a; 
          
  for (x=x0-h; x<=x0+h; x++)
     for (y=y0-h; y<=y0+h; y++)
       {
         a = allsymbol[ind][px(x-x0+h+1, y-y0+h+1, 2*h+1, 2*h+1)]*w; 
         if (sym[px(x, y, sx, sy)]<a)
             sym[px(x, y, sx, sy)] = a;
       }
}


/* local maximization */
void Cmm(int i)
{
   int  ind, x, y, rm[5];  

   for (ind=0; ind<N; ind++)
     { 
        for (x=floor((h+1)/sub)+1+L; x<=floor((sx-h)/sub)-1-L; x++)
           for (y=floor((h+1)/sub)+1+L; y<=floor((sy-h)/sub)-1-L; y++)
             {
                   mI[ind*n+i][px(x*sub, y*sub, sx, sy)] = shiftmax(i, ind, x*sub, y*sub, rm);   /* local maximum of active Gabor */  
               }
      }
}


/* the shared sketch algorithm */
void Cmp(int i)
{
   int x0, y0, x, y, mi, mx, my, t, Fx, Fy, rm[5]; 
   double F, mm, ave; 

   t = 0; 
   do  
   {
     mi = Mi[t]; mx = Mx[t]; my = My[t];  
     /* the following is a parallel search */
    for (x=1; x<floor((sx-ax*2)/sub); x++)
         for (y=1; y<floor((sy-ay*2)/sub); y++)
          {            
           x0 = x*sub; y0 = y*sub; 
           if (LIK)
           {
           tss[i][px(x0+ax, y0+ay, sx, sy)] +=  (1. - 1./Mm[t])*mI[mi*n+i][px(x0+mx, y0+my, sx, sy)] - log(Mm[t]); 
           }
           else 
           {
           tss[i][px(x0+ax, y0+ay, sx, sy)] += sqrt(mI[mi*n+i][px(x0+mx, y0+my, sx, sy)])*Mm[t]; 
           }
         }   
     t++; 
     }
  while (t<Totalsketch);   /* stopping criterion */
   
  /* find the center pixel of the best score */ 
  F = NEGMAX; Fx = 0; Fy = 0;
  for (x=1; x<floor((sx-ax*2)/sub); x++)
         for (y=1; y<floor((sy-ay*2)/sub); y++)
          {            
             x0 = x*sub; y0 = y*sub; 
             if (F<tss[i][px(x0+ax, y0+ay, sx, sy)]) 
             {
                 F = tss[i][px(x0+ax, y0+ay, sx, sy)];
                 Fx = x0; 
                 Fy = y0; 
             }
         }  
   mscale[i] = F; 
  
   t = 0; 
   do  
   {
     mi = Mi[t]; mx = Mx[t]; my = My[t];   
     mm = shiftmax(i, mi, Fx+mx, Fy+my, rm);  
     draw(Asym[i], rm[1], rm[2], rm[3], sqrt(mm)); 
     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 i, ind, j, c; 
 mxArray *f;  
 
 c = 0; /* counter for input pointers and scalars */
 n = floor(mxGetScalar(prhs[c++])+.5);
 N = floor(mxGetScalar(prhs[c++])+.5);  /* number of orientations */
  
 fI = mxCalloc(n*N, sizeof(double*));   /* fitered images */
 for (i=0; i<n; i++)
   {
     for (ind=0; ind<N; ind++)
      {  
       f = mxGetCell(prhs[c], ind*n+i); 
       fI[ind*n+i] = mxGetPr(f);    /* get pointers to filtered images */       
      }
    }
 c++;

 mI = mxCalloc(n*N, sizeof(double*));   /* fitered images */
 for (i=0; i<n; i++)
   {
     for (ind=0; ind<N; ind++)
      {  
       f = mxGetCell(prhs[c], ind*n+i); 
       mI[ind*n+i] = mxGetPr(f);    /* get pointers to filtered images */       
      }
    }
 c++;
 
 
 h = floor(mxGetScalar(prhs[c++])+.5);    /* half size of filters */
 
 ax = floor(mxGetScalar(prhs[c++])+.5);
 ay = floor(mxGetScalar(prhs[c++])+.5);   /* size of window */
 
     
 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++]);   /* average response of selected Gabor */ 
 
 allsymbol = mxCalloc(N, sizeof(double*));    
 for (ind=0; ind<N; ind++)
     {  
       f = mxGetCell(prhs[c], ind); 
       allsymbol[ind] = mxGetPr(f);  /* symbols of filters */         
     }
 c++; 
 
 L = floor(mxGetScalar(prhs[c++])+.5);     /* range of shifting along normal direction of Gabor */
 ore = floor(mxGetScalar(prhs[c++])+.5);   /* range of shifting in orientation */
 sub = floor(mxGetScalar(prhs[c++])+.5);   /* sub-sampling of Gabor filters to improve speed */
 Totalsketch = floor(mxGetScalar(prhs[c++])+.5);  /* total number of bases */
 Upperbound = mxGetScalar(prhs[c++]);             /* satuation level */
 
 tss = mxCalloc(n, sizeof(double*)); 
  for (i=0; i<n; i++)
   {
     f = mxGetCell(prhs[c], i); 
     tss[i] = mxGetPr(f);  
  } 
 c++; 
/* map of log-likelihood score */
 Asym = mxCalloc(n, sizeof(double*)); 
  for (i=0; i<n; i++)
   {
     f = mxGetCell(prhs[c], i); 
     Asym[i] = mxGetPr(f);  
  } 
 c++; 
 
 Sx1 = mxGetPr(prhs[c++]);
 Sy1 = mxGetPr(prhs[c++]);
 mscale =  mxGetPr(prhs[c++]);
 
 LIK = floor(mxGetScalar(prhs[c++])+.5);    /* log-likelihood or active correlation? */
 

 storeshift(); 
 for (i=0; i<n; i++)
 {
 sx = floor(Sx1[i]+.5); 
 sy = floor(Sy1[i]+.5); 
 Cmm(i);      /* parallel local maximum pooling */  
 Cmp(i);      /* parallel scoring, like a deformable filter */  
 }
}

     



 

                    