123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162 |
- /*
- * 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;
- }
- }
|