|
@@ -0,0 +1,162 @@
|
|
|
+/*
|
|
|
+ * 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;
|
|
|
+ }
|
|
|
+}
|