 	dev_normalize<<<gridm, 128>>>(d_resp, d_pop, nrX*nrY, nrDirs);
 	cutilCheckMsg("dev_normalize() execution failed\n");
 
+	// FIXME: i have no idea where the numbers 1e6, 5e5, and 1e5 are coming from...
+	// FIXME: also, there's no need to de-grayscale the inputs, this has already been done dev_split()
+	// but the idea is to normalize the responses and linearly map them onto [0,50 Hz]
 	for (int scale=0; scale<nrScales; scale++)
 		dev_scale<<<iDivUp(nrX*nrY*nrDirs,128), 128>>>(&d_resp[scale*nrX*nrY*nrDirs], (scale==0?1000000:(scale==1?500000:100000))/255.0/255*50, nrX*nrY*nrDirs);
 

	printf("used GPU memory %ld\n", total - free);
*/
}


#include <cutil_inline.h>
#define IMUL(a, b) __mul24(a, b)

// 28 space-time orientations of V1 simple cells
#define nrDirs 28

// each of these is a unit vector (e.g. d_v1Dirs[0][1..3] => (-0.6559)^2+(0.7246)^2+(0.2113)^2 == 1
// all vectors lie on the surface of a dome (hemisphere) in 3D Fourier space
__constant__  float d_v1Dirs[3][nrDirs] = {{-0.6559, -0.1019, 0.6240, -0.7797, 0.9692, -0.2312, -0.9151, 0.4207, -0.9533, 0.8175, 0.2398, 0.8810, -0.4430, 0.0588, -0.5384, 0.5644, 0.7931, 0.5142, -0.7680, -0.0669, -0.6670, -0.2747, 0.5034, 0.5042, 0.1580, 0.1332, -0.5159, -0.3549},
                                           { 0.7246, -0.9718, 0.7496, -0.5837, -0.0810, 0.9439, 0.3203, -0.8712, -0.1593, -0.5142, 0.9304, 0.3737, -0.8031, -0.8126, 0.6004, -0.5738, 0.0024, 0.5969, 0.1436, 0.7757, -0.4004, -0.5108, 0.2375, -0.2221, -0.5140, 0.5194, -0.0870, 0.3838},
                                           { 0.2113, 0.2126, 0.2210, 0.2266, 0.2327, 0.2359, 0.2451, 0.2529, 0.2567, 0.2593, 0.2772, 0.2902, 0.3984, 0.5799, 0.5913, 0.5935, 0.6091, 0.6160, 0.6241, 0.6275, 0.6283, 0.8146, 0.8308, 0.8345, 0.8431, 0.8441, 0.8522, 0.8525}};


// this filter is used for the 3 different scales in space/time:
// first scale  == original image resolution
// second scale == first scale blurred with Gaussian (width 1), so resolution should scale down sqrt(2)
// third scale  == second scale blurred with Gaussian (width 1)
// FIXME: not sure why this guy is not normalized
#define scalingFiltSize 5
__constant__ float d_scalingFilt[scalingFiltSize] = {0.0884, 0.3536, 0.5303, 0.3536, 0.0884};
float* scalingFilt;

// d_v1Gaus defines the 1D receptive field of a V1 unit, which is then used for all three dimensions (X,Y and T)
// this guy can be reproduced in matlab with g=normpdf(-4:4,0,1.25);
// it is the same as provided by the S&H matlab code
#define v1GausSize 9
__constant__ float d_v1Gaus[v1GausSize] = {0.0007, 0.0155, 0.0903, 0.2345, 0.3179, 0.2345, 0.0903, 0.0155, 0.0007};
float* v1Gaus;

// d_complexV1Filt is the spacial filter for complex cells; it averages over "simple" V1 cells
// all simple cells must have the same space-time orientation and phase
// this guy can be reproduced in matlab with g=normpdf(-5:5,0,1.6);
// it is the same as provided by the S&H matlab code
#define complexV1FiltSize 11
__constant__ float d_complexV1Filt[complexV1FiltSize] = {0.0019, 0.0110, 0.0430, 0.1142, 0.2052, 0.2495, 0.2052, 0.1142, 0.0430, 0.0110, 0.0019};
float* complexV1Filt;

// d_normV1filt is the spatial filter used complex cell normalization
// this guy can be reproduced in matlab with: g=normpdf(-10:10,0,3.35);
#define normV1filtSize 21
__constant__ float d_normV1filt[normV1filtSize] = {0.0013, 0.0031, 0.0067, 0.0132, 0.0237, 0.0389, 0.0584, 0.0800, 0.1001, 0.1146, 0.1199, 0.1146, 0.1001, 0.0800, 0.0584, 0.0389, 0.0237, 0.0132, 0.0067, 0.0031, 0.0013};
float* normV1filt;

// difference operator for taking the first-order derivative
#define diff1filtSize 3
__constant__ float d_diff1filt[diff1filtSize] = {-1/2.0, 0, 1/2.0};
float* diff1filt;

// difference operator for taking the second-order derivative
#define diff2filtSize 3
__constant__ float d_diff2filt[diff2filtSize] = {1, -2, 1};
float* diff2filt;

// difference operator for taking the third-order derivative
// the grayscale values of our input stimuli will be convolved with d_scalingFilt in 3D
#define diff3filtSize 5
__constant__ float d_diff3filt[diff3filtSize] = {-1/2.0, 1, 0, -1, 1/2.0};
float* diff3filt;

// 3 different spatio-temporal scales (these must correspond to the 3 blobs along the x-, y-, or t-axis in Fig. 3 in
// the Simoncelli & Heeger paper... because they are scales in the x,y,t-space they cannot be temporal frequency or
// "speed" alone)
#define nrScales 3

// number of time steps to be considered in computation
// nrT = 5*(3-1)-1 = 9
#define nrT (scalingFiltSize*(nrScales-1)-1)//(scalingFiltSize*nrScales-2)

int stimBufX, stimBufY;
float* d_resp; // will be the returned ME responses
float* d_stimBuf;
float* d_scalingStimBuf; // the temporary matrix that will be filtered once for each scale...
float* d_v1GausBuf;
