 			            yIndex * cols +            
 						xIndex;
 
-		output[idx]              = blue;
-		output[idx +   chanDist] = green;
-		output[idx + 2*chanDist] = red;
+		output[idx]                = blue;
+		output[idx +     chanDist] = green;
+		output[idx + 2 * chanDist] = red;
 	}
 }
 
 // Math to add two intermediate steps of mean & stddev 
 // See http://www.johndcook.com/blog/skewness_kurtosis/
-__device__ void combine_running_totals(float &M1_1, float M1_2, float &M2_1, float M2_2, unsigned int &n_1, unsigned int n_2)
+__device__ void combine_running_totals(float &M1_1, const float M1_2, float &M2_1, const float M2_2, unsigned int &n_1, const unsigned int n_2)
 {
 	unsigned int combined_n = n_1 + n_2;
 
-	const float delta = M1_2 - M1_1;
+	const float delta  = M1_2 - M1_1;
 	const float delta2 = delta * delta;
 
 	float combined_M1 = (n_1 * M1_1 + n_2 * M1_2) / combined_n;
 	float combined_M2 = M2_1 + M2_2 + delta2 * n_1 * n_2 / combined_n;
 
-	n_1 = combined_n;
+	n_1  = combined_n;
 	M1_1 = combined_M1;
 	M2_1 = combined_M2;
 }



// For 24x24 reductions, there are 4 entries x 3 channels
// left over after a 16x16 thread is reduced (i.e. 24x24 fills
// 4 thread blocks). Reduce those 4 values down to one
__global__ void mean_stddev_reduction_kernel24(const float *M1In, const float *M2In, const unsigned int *nIn,
		float *mean, float *stddev)
{

	const int zIndex = blockIdx.z * blockDim.z + threadIdx.z;

	//Only valid threads perform memory I/O
	if((xIndex < input[zIndex].cols) && (yIndex < input[zIndex].rows))
	{
		// xIndex * 3 since col has a blue green and red component
		const float blue	= input[zIndex](yIndex, 3*xIndex);
