CheckerFrames.cpp 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106
  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. // function declarations
  21. double ran1(long &idum);
  22. void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  23. {
  24. if (nrhs != 3) mexErrMsgTxt("Input arguments are stixel width and stixel height."); // number of frames and seed.");
  25. if (nlhs != 1) mexErrMsgTxt("Wrong number of output arguments.");
  26. // required input arguments
  27. int xstixels, ystixels, nframes;
  28. long seed = -10000;
  29. float meanintensity = 0.5;
  30. float contrast = 1;
  31. float lowInt = -1; // (1 - contrast)*meanintensity;
  32. float highInt = 1; // (1 + contrast)*meanintensity;
  33. //get input data
  34. if (mxIsDouble(prhs[0])) xstixels = (int)mxGetScalar(prhs[0]);
  35. else mexErrMsgTxt("Invalid class of input argument 'xpixels'.");
  36. if (mxIsDouble(prhs[1])) ystixels = (int)mxGetScalar(prhs[1]);
  37. else mexErrMsgTxt("Invalid class of input argument 'ypixels'.");
  38. if (mxIsDouble(prhs[2])) nframes = (int)mxGetScalar(prhs[2]);
  39. else mexErrMsgTxt("Invalid class of input argument 'nframes'.");
  40. // create pointer array
  41. double *stimulus = new double[xstixels*ystixels*nframes];
  42. // generate array of random numbers
  43. for (int n=0; n<nframes; n++) {
  44. for (int x=0; x<xstixels; x++) {
  45. for (int y=0; y<ystixels; y++) {
  46. stimulus[n*xstixels*ystixels + x*ystixels + y] = ran1( seed ) < 0.5 ? lowInt : highInt;
  47. // testing of x and y coordinates
  48. // if (x == 2) stimulus[n*xstixels*ystixels + x*ystixels + y] = 1;
  49. }
  50. }
  51. }
  52. // allocate space for output mxArray
  53. int dims[3] = {ystixels, xstixels, nframes};
  54. plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
  55. // plhs[0] = mxCreateDoubleMatrix(xstixels*ystixels*nframes, 1, mxREAL);
  56. memcpy(mxGetPr(plhs[0]), stimulus, xstixels*ystixels*nframes*sizeof(double));
  57. // delete pointer array
  58. delete [] stimulus;
  59. }
  60. double ran1(long &idum)
  61. {
  62. int j;
  63. long k;
  64. static long iy=0;
  65. static long iv[NTAB];
  66. double temp;
  67. if (idum <= 0 || !iy) {
  68. if (-(idum) < 1) idum=1;
  69. else idum = -(idum);
  70. for (j=NTAB+7;j>=0;j--) {
  71. k=(idum)/IQ;
  72. idum=IA*(idum-k*IQ)-IR*k;
  73. if (idum < 0) idum += IM;
  74. if (j < NTAB) iv[j] = idum;
  75. }
  76. iy=iv[0];
  77. }
  78. k=(idum)/IQ;
  79. idum=IA*(idum-k*IQ)-IR*k;
  80. if (idum < 0) idum += IM;
  81. j=iy/NDIV;
  82. iy=iv[j];
  83. iv[j] = idum;
  84. if ((temp=AM*iy) > RNMX) return RNMX;
  85. else return temp;
  86. }