/*
 * Copyright (c) 2011, Norma Kuehn
 * All rights reserved.
 *
 */

// Pseudo random number generator 'ran1' used in the stimulation program.

#include <mex.h>
#include <memory.h>
#include <matrix.h>
#define _USE_MATH_DEFINES 
#include <math.h>

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

#define squareWidth 800
#define squareHeight 800

int iset = 0;
double gset = 0;


// function declarations
double gauss(long &idum);

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    if (nrhs < 2) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
    if (nrhs > 3) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
    if (nlhs != 1) mexErrMsgTxt("Wrong number of output arguments.");
    
    // required input arguments
    int bgStixel, filterSTDV;
    
    // defined constants
    //int squareWidth = 800;
    //int squareHeight = 800;
    
    long seed = -10000;
    
    float meanintensity = 0.5;
    float contrast = 1; // 0.3;

    //get input data
	if (mxIsDouble(prhs[0])) bgStixel = (int)mxGetScalar(prhs[0]);
	else mexErrMsgTxt("Invalid class of input argument 'bgStixel'.");
    
	if (mxIsDouble(prhs[1])) filterSTDV = (int)mxGetScalar(prhs[1]);
	else mexErrMsgTxt("Invalid class of input argument 'filterSTDV'.");

    // if third argument is given, it's a new seed
    if (nrhs == 3) {
	  if (mxIsDouble(prhs[2])) {
          seed = (int)mxGetScalar(prhs[2]);
          if (seed > 0) seed = -seed;
      }
      else mexErrMsgTxt("Invalid class of input argument 'seed'.");
    }

    // create pointer array
    int xstixels = squareWidth/bgStixel;
    int ystixels = squareHeight/bgStixel;
    double *stimulus = new double[xstixels*ystixels];
    
	// generate array of random numbers
    int filterWidth = filterSTDV / bgStixel * 3;
    int noiseLim[2] = { ceil(float(squareHeight) / bgStixel), ceil(float(squareWidth) / bgStixel) };
    static double noiseField[squareHeight][squareWidth];
    static double Gfilter[squareHeight][squareWidth];
    float norm = 0.;
    for (int i = 0; i < noiseLim[0]; i++) {
        for (int j = 0; j < noiseLim[1]; j++) {
            noiseField[i][j] = meanintensity + meanintensity*contrast*gauss(seed);
            if ((i < (filterWidth * 2 + 1)) && (j < (filterWidth * 2 + 1))) {
                Gfilter[i][j] = exp2f(-(pow(float(i - filterWidth), 2) + pow(float(j - filterWidth), 2)) / (2. * pow(float(filterSTDV) / bgStixel, 2)));
                norm += Gfilter[i][j];
            }
        }
    }
    for (int i = 0; i < noiseLim[0]; i++) {
        for (int j = 0; j < noiseLim[1]; j++) {
            double help = 0.;
            for (int ki = 0; ki < (filterWidth * 2 + 1); ki++) {
                for (int kj = 0; kj < (filterWidth * 2 + 1); kj++) {
                    help += noiseField[(i - (ki - filterWidth) + noiseLim[0]) % noiseLim[0]][(j - (kj - filterWidth) + noiseLim[1]) % noiseLim[1]] * Gfilter[ki][kj];
                }
            }
            stimulus[j*ystixels + i] = (help / norm - meanintensity) *filterSTDV / bgStixel + meanintensity;
        }
    }
    
    // allocate space for output mxArray
    int dims[2] = {ystixels, xstixels};
	plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
    // plhs[0] = mxCreateDoubleMatrix(xstixels*ystixels*nframes, 1, mxREAL);

    memcpy(mxGetPr(plhs[0]), stimulus, xstixels*ystixels*sizeof(double));
    
    // delete pointer array
    delete [] stimulus;
}


double ran1(long &idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if (idum <= 0 || !iy) {
		if (-(idum) < 1) idum=1;
		else idum = -(idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(idum)/IQ;
			idum=IA*(idum-k*IQ)-IR*k;
			if (idum < 0) idum += IM;
			if (j < NTAB) iv[j] = idum;
		}
		iy=iv[0];
	}
	k=(idum)/IQ;
	idum=IA*(idum-k*IQ)-IR*k;
	if (idum < 0) idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}

double gauss(long &idum)
{
	double fac,rsq,v1,v2;

	if (idum < 0) iset=0;
	if  (iset == 0) {
		do {
			v1=2.0*ran1(idum)-1.0;
			v2=2.0*ran1(idum)-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} 
    else {
		iset=0;
		return gset;
	}
}