 	float2* pOutputPixel = &( outputGrid( gx, gy, gz ) );
 	atomicAdd( &( pOutputPixel->x ), lum );
 	atomicAdd( &( pOutputPixel->y ), 1 );
+
+    #else
+	int gx = blockIdx.x * blockDim.x + threadIdx.x;
+	int gy = blockIdx.y * blockDim.y + threadIdx.y;
+    
+    for (int ry = 0; ry < samplingSpatial; ry++) {
+        for (int rx = 0; rx < samplingSpatial; rx++) {
+            int x = gx * samplingSpatial + rx - samplingSpatial/2;
+            int y = gy * samplingSpatial + ry - samplingSpatial/2;
+                
+        	if( x >= inputImage.width ||
+        		y >= inputImage.height )
+        	{
+        		continue;
+        	}
+        	
+            float lum = inputImage( x, y );
+        	lum = ( lum - inputMin ) / inputRange;
+        	int gz = roundToInt( lum / samplingRange ) + paddingZ;
+            
+            float *pOutputPixel, *pOutputNorm;
+            #if 0
+            if (gz % 2 == 0) {
+        	    pOutputPixel = &( outputGrid( gx, gy, gz/2 ).x );
+                pOutputNorm = &( outputGrid( gx, gy, gz ).x );
+            } else {
+        	    pOutputPixel = &( outputGrid( gx, gy, gz/2 ).y );
+                pOutputNorm = &( outputGrid( gx, gy, gz ).y );
+            }
+            #else
+    	    pOutputPixel = &( outputGrid( gx, gy, gz ).x );
+            pOutputNorm = &( outputGrid( gx, gy, gz ).y );
+            #endif
+            
+            #if 0
+        	atomicAdd( pOutputPixel, lum );
+            atomicAdd( pOutputNorm, 1 );
+            #else
+        	*pOutputPixel += lum;
+            *pOutputNorm += 1;
+            #endif
+        }
+    }
+    #endif
 }
 
 __global__

					  int paddingXY, int paddingZ,
					  DeviceArray3D< float2 > outputGrid )
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;


#include <MathUtils.h>

#include <cutil_math.h>

#define JRK_SCHEDULE 1

// TODO: be very careful with this variable
texture< float2, cudaTextureType2D, cudaReadModeElementType > g_gridTexture( 0, cudaFilterModeLinear );
