/**************************************************
 Mex C Code for computing the sigmoid transformation
***************************************************/
# 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 images */
 int N;                      /* number of orientations */
 double **fI;                /* filtered images */
 int sx, sy;                 /* sizes of image */  
 double **mI;                /* locally maximized images */
 int **sinsh, **cossh;       /* store the shift to avoid repeated sin and cosin computation */
 int L, ore, sub;                 /* allowed ranges of shifting in location and orientation */
 
/* 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 = 0.; 
   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(m); 
}

 
 
void Cmm()
{  
   int x, y, here, ind, i, ii, rm[5]; 
   
   mI = double_matrix(N, sx*sy);  
   
   for (i=0; i<n; i++)
   { 
    for (ind=0; ind<N; ind++)
    {                   
       for (x=L*sub+1; x<sx-L*sub-1; x++)
           for (y=L*sub+1; y<sy-L*sub-1; y++)
            {
               here = px(x, y, sx, sy);   
               mI[ind][here] = shiftmax(i, ind, x, y, rm); 
            }
    }
   for (ind=0; ind<N; ind++)
    {                 
       ii = ind*n+i;     
       for (x=L*sub+1; x<sx-L*sub-1; x++)
           for (y=L*sub+1; y<sy-L*sub-1; y++)
            {
               here = px(x, y, sx, sy);   
               fI[ii][here] = mI[ind][here]; 
            }
    }   
  }
   
  free_matrix(mI, N, sx*sy);
}
   


/* 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, i, c; 
 mxArray *f;  
 
 c = 0; /* counter for input pointers and scalars */
 n = floor(mxGetScalar(prhs[c++])+.5);  /* number of images */
 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++; 
 
 sx = floor(mxGetScalar(prhs[c++])+.5);
 sy = floor(mxGetScalar(prhs[c++])+.5);   /* size of images */
 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 */

 storeshift();
 Cmm();
}

     



 

                    