GaussianFrames.cpp 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137
  1. /*
  2. * Copyright (c) 2011, Norma Kuehn
  3. * All rights reserved.
  4. *
  5. */
  6. // Pseudo random number generator 'ran1' used in the stimulation program.
  7. #include <mex.h>
  8. #include <memory.h>
  9. #include <matrix.h>
  10. #include <math.h>
  11. #define IA 16807
  12. #define IM 2147483647
  13. #define AM (1.0/IM)
  14. #define IQ 127773
  15. #define IR 2836
  16. #define NTAB 32
  17. #define NDIV (1+(IM-1)/NTAB)
  18. #define EPS 1.2e-7
  19. #define RNMX (1.0-EPS)
  20. int iset = 0;
  21. double gset = 0;
  22. // function declarations
  23. double gauss(long &idum);
  24. void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  25. {
  26. if (nrhs < 3) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
  27. if (nrhs > 4) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
  28. if (nlhs != 1) mexErrMsgTxt("Wrong number of output arguments.");
  29. // required input arguments
  30. int xstixels, ystixels, nframes;
  31. long seed = -1000;
  32. float meanintensity = 0.5;
  33. float contrast = 1; // 0.3;
  34. //get input data
  35. if (mxIsDouble(prhs[0])) xstixels = (int)mxGetScalar(prhs[0]);
  36. else mexErrMsgTxt("Invalid class of input argument 'xpixels'.");
  37. if (mxIsDouble(prhs[1])) ystixels = (int)mxGetScalar(prhs[1]);
  38. else mexErrMsgTxt("Invalid class of input argument 'ypixels'.");
  39. if (mxIsDouble(prhs[2])) nframes = (int)mxGetScalar(prhs[2]);
  40. else mexErrMsgTxt("Invalid class of input argument 'nframes'.");
  41. if (nrhs == 4) {
  42. if (mxIsDouble(prhs[3])) {
  43. seed = (int)mxGetScalar(prhs[3]);
  44. if (seed > 0) seed = -seed;
  45. }
  46. else mexErrMsgTxt("Invalid class of input argument 'nframes'.");
  47. }
  48. // create pointer array
  49. double *stimulus = new double[xstixels*ystixels*nframes];
  50. // generate array of random numbers
  51. for (int n=0; n<nframes; n++) {
  52. for (int x=0; x<xstixels; x++) {
  53. for (int y=0; y<ystixels; y++) {
  54. stimulus[n*xstixels*ystixels + x*ystixels + y] = meanintensity*(1 + contrast*gauss( seed ));
  55. }
  56. }
  57. }
  58. // allocate space for output mxArray
  59. int dims[3] = {ystixels, xstixels, nframes};
  60. plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
  61. // plhs[0] = mxCreateDoubleMatrix(xstixels*ystixels*nframes, 1, mxREAL);
  62. memcpy(mxGetPr(plhs[0]), stimulus, xstixels*ystixels*nframes*sizeof(double));
  63. // delete pointer array
  64. delete [] stimulus;
  65. }
  66. double ran1(long &idum)
  67. {
  68. int j;
  69. long k;
  70. static long iy=0;
  71. static long iv[NTAB];
  72. double temp;
  73. if (idum <= 0 || !iy) {
  74. if (-(idum) < 1) idum=1;
  75. else idum = -(idum);
  76. for (j=NTAB+7;j>=0;j--) {
  77. k=(idum)/IQ;
  78. idum=IA*(idum-k*IQ)-IR*k;
  79. if (idum < 0) idum += IM;
  80. if (j < NTAB) iv[j] = idum;
  81. }
  82. iy=iv[0];
  83. }
  84. k=(idum)/IQ;
  85. idum=IA*(idum-k*IQ)-IR*k;
  86. if (idum < 0) idum += IM;
  87. j=iy/NDIV;
  88. iy=iv[j];
  89. iv[j] = idum;
  90. if ((temp=AM*iy) > RNMX) return RNMX;
  91. else return temp;
  92. }
  93. double gauss(long &idum)
  94. {
  95. double fac,rsq,v1,v2;
  96. if (idum < 0) iset=0;
  97. if (iset == 0) {
  98. do {
  99. v1=2.0*ran1(idum)-1.0;
  100. v2=2.0*ran1(idum)-1.0;
  101. rsq=v1*v1+v2*v2;
  102. } while (rsq >= 1.0 || rsq == 0.0);
  103. fac=sqrt(-2.0*log(rsq)/rsq);
  104. gset=v1*fac;
  105. iset=1;
  106. return v2*fac;
  107. }
  108. else {
  109. iset=0;
  110. return gset;
  111. }
  112. }