bgTexture.cpp 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162
  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. #define _USE_MATH_DEFINES
  11. #include <math.h>
  12. #define IA 16807
  13. #define IM 2147483647
  14. #define AM (1.0/IM)
  15. #define IQ 127773
  16. #define IR 2836
  17. #define NTAB 32
  18. #define NDIV (1+(IM-1)/NTAB)
  19. #define EPS 1.2e-7
  20. #define RNMX (1.0-EPS)
  21. #define squareWidth 800
  22. #define squareHeight 800
  23. int iset = 0;
  24. double gset = 0;
  25. // function declarations
  26. double gauss(long &idum);
  27. void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  28. {
  29. if (nrhs < 2) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
  30. if (nrhs > 3) mexErrMsgTxt("Input arguments are stixel width, stixel height, number of frames [and seed].");
  31. if (nlhs != 1) mexErrMsgTxt("Wrong number of output arguments.");
  32. // required input arguments
  33. int bgStixel, filterSTDV;
  34. // defined constants
  35. //int squareWidth = 800;
  36. //int squareHeight = 800;
  37. long seed = -10000;
  38. float meanintensity = 0.5;
  39. float contrast = 1; // 0.3;
  40. //get input data
  41. if (mxIsDouble(prhs[0])) bgStixel = (int)mxGetScalar(prhs[0]);
  42. else mexErrMsgTxt("Invalid class of input argument 'bgStixel'.");
  43. if (mxIsDouble(prhs[1])) filterSTDV = (int)mxGetScalar(prhs[1]);
  44. else mexErrMsgTxt("Invalid class of input argument 'filterSTDV'.");
  45. // if third argument is given, it's a new seed
  46. if (nrhs == 3) {
  47. if (mxIsDouble(prhs[2])) {
  48. seed = (int)mxGetScalar(prhs[2]);
  49. if (seed > 0) seed = -seed;
  50. }
  51. else mexErrMsgTxt("Invalid class of input argument 'seed'.");
  52. }
  53. // create pointer array
  54. int xstixels = squareWidth/bgStixel;
  55. int ystixels = squareHeight/bgStixel;
  56. double *stimulus = new double[xstixels*ystixels];
  57. // generate array of random numbers
  58. int filterWidth = filterSTDV / bgStixel * 3;
  59. int noiseLim[2] = { ceil(float(squareHeight) / bgStixel), ceil(float(squareWidth) / bgStixel) };
  60. static double noiseField[squareHeight][squareWidth];
  61. static double Gfilter[squareHeight][squareWidth];
  62. float norm = 0.;
  63. for (int i = 0; i < noiseLim[0]; i++) {
  64. for (int j = 0; j < noiseLim[1]; j++) {
  65. noiseField[i][j] = meanintensity + meanintensity*contrast*gauss(seed);
  66. if ((i < (filterWidth * 2 + 1)) && (j < (filterWidth * 2 + 1))) {
  67. Gfilter[i][j] = exp2f(-(pow(float(i - filterWidth), 2) + pow(float(j - filterWidth), 2)) / (2. * pow(float(filterSTDV) / bgStixel, 2)));
  68. norm += Gfilter[i][j];
  69. }
  70. }
  71. }
  72. for (int i = 0; i < noiseLim[0]; i++) {
  73. for (int j = 0; j < noiseLim[1]; j++) {
  74. double help = 0.;
  75. for (int ki = 0; ki < (filterWidth * 2 + 1); ki++) {
  76. for (int kj = 0; kj < (filterWidth * 2 + 1); kj++) {
  77. help += noiseField[(i - (ki - filterWidth) + noiseLim[0]) % noiseLim[0]][(j - (kj - filterWidth) + noiseLim[1]) % noiseLim[1]] * Gfilter[ki][kj];
  78. }
  79. }
  80. stimulus[j*ystixels + i] = (help / norm - meanintensity) *filterSTDV / bgStixel + meanintensity;
  81. }
  82. }
  83. // allocate space for output mxArray
  84. int dims[2] = {ystixels, xstixels};
  85. plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
  86. // plhs[0] = mxCreateDoubleMatrix(xstixels*ystixels*nframes, 1, mxREAL);
  87. memcpy(mxGetPr(plhs[0]), stimulus, xstixels*ystixels*sizeof(double));
  88. // delete pointer array
  89. delete [] stimulus;
  90. }
  91. double ran1(long &idum)
  92. {
  93. int j;
  94. long k;
  95. static long iy=0;
  96. static long iv[NTAB];
  97. double temp;
  98. if (idum <= 0 || !iy) {
  99. if (-(idum) < 1) idum=1;
  100. else idum = -(idum);
  101. for (j=NTAB+7;j>=0;j--) {
  102. k=(idum)/IQ;
  103. idum=IA*(idum-k*IQ)-IR*k;
  104. if (idum < 0) idum += IM;
  105. if (j < NTAB) iv[j] = idum;
  106. }
  107. iy=iv[0];
  108. }
  109. k=(idum)/IQ;
  110. idum=IA*(idum-k*IQ)-IR*k;
  111. if (idum < 0) idum += IM;
  112. j=iy/NDIV;
  113. iy=iv[j];
  114. iv[j] = idum;
  115. if ((temp=AM*iy) > RNMX) return RNMX;
  116. else return temp;
  117. }
  118. double gauss(long &idum)
  119. {
  120. double fac,rsq,v1,v2;
  121. if (idum < 0) iset=0;
  122. if (iset == 0) {
  123. do {
  124. v1=2.0*ran1(idum)-1.0;
  125. v2=2.0*ran1(idum)-1.0;
  126. rsq=v1*v1+v2*v2;
  127. } while (rsq >= 1.0 || rsq == 0.0);
  128. fac=sqrt(-2.0*log(rsq)/rsq);
  129. gset=v1*fac;
  130. iset=1;
  131. return v2*fac;
  132. }
  133. else {
  134. iset=0;
  135. return gset;
  136. }
  137. }