# include <stdio.h>
# include <stdlib.h>
# include "mex.h"        
# include "math.h"
# define ROUND(x) (floor((x)+.5))
/* 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 */
float *I, *J, average; 
double saturation;               
int sx, sy, sx1, sy1, Fx, Fy; 
/* sigmoid transformation */
float Sigmoid(float r)
{
   return(saturation*(2./(1.+exp(-2.*r/saturation))-1.)); 
}
/* copy images */
void CopyImage()
{
   
   int x, y, here, there, x1, y1; 
   
   for (x=1; x<=sx; x++)
       for (y=1; y<=sy; y++)
       {
        here = px(x, y, sx, sy); 
        x1 = Fx+x; y1 = Fy+y; 
        if ((x1>=1)&&(x1<=sx1)&&(y1>=1)&&(y1<=sy1))
        {
           there = px(x1, y1, sx1, sy1); 
           I[here] = Sigmoid(J[there]/average);
        }
        else
            I[here] = 0.; 
        }
}
/* read in images */  
void mexFunction(int nlhs, mxArray *plhs[], 
                 int nrhs, const mxArray *prhs[])                
{
 int c; 
 
 c = 0;
 I = mxGetPr(prhs[c++]);
 J = mxGetPr(prhs[c++]);
 Fx = ROUND(mxGetScalar(prhs[c++]));
 Fy = ROUND(mxGetScalar(prhs[c++])); 
 sx = ROUND(mxGetScalar(prhs[c++]));
 sy = ROUND(mxGetScalar(prhs[c++]));   
 sx1 = ROUND(mxGetScalar(prhs[c++]));
 sy1 = ROUND(mxGetScalar(prhs[c++]));   
 average = mxGetScalar(prhs[c++]);
 saturation = mxGetScalar(prhs[c++]);
 CopyImage();
}

     



 

                    