 // SIM is the product of the gradient ("S") and image magnitude ("IM"). The
 // metric considers sharp edgegs which are dark to be favorable. Since "dark" 
 // has a smaller numerical value, the image brightness, or local magnitude, is
-// subtracted from the expected maximum intensity so that (intmax - localmag)
+// subtracted from the expected maximum intensity so that 
+//
+//   (maximum_intensity - image_magnitude)
+//
 // gives a large value for favorable dark pixels. The SIM is then
 //
-//   SIM = sfiltmag * (intmax - localmag).
+//   SIM = gradient_magnitude * (maximum_intensity - image_magnitude).
 //
 // sim: output array of the SIM focus metric
-// sfiltmag: magnitude of the steerable filter (a gradient filter), applied
-//           to the reconstructed image
-// localmag: magnitude of a single slice's reconstructed image
-// intmax: maximum intensity expected over all slices; should be a constant
-//           across slices
-// M: number of rows in sim, sfiltmag, and localmag
-// N: number of cols in sim, sfiltmag, and localmag
-static __global__ void computeSIMkernel(float* sim, float* sfiltmag,
-                                        float* localmag, const float intmax, 
+// gradient_magnitude: magnitude of the a gradient filter which has been applied
+//           to the complex-valued reconstructed image
+// image_magnitude: magnitude of a single slice's reconstructed image
+// maximum_intensity: maximum intensity expected over all slices; should be a
+//           constant across slices
+// M: number of rows in sim, gradient_magnitude, and image_magnitude
+// N: number of cols in sim, gradient_magnitude, and image_magnitude
+static __global__ void computeSIMkernel(float* sim, 
+	                                    float* gradient_magnitude,
+                                        float* image_magnitude,
+                                        const float maximum_intensity, 
 										const int M, const int N) {
 	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
 	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
 	if (xidx < N && yidx < M) {
 		int idx = yidx * N + xidx;
-		sim[idx] = sfiltmag[idx] * (intmax - localmag[idx]);
-	}
-}
-
-
-
-static __global__ void addShortKernel(short* dst, const short value,
+		sim[idx] = gradient_magnitude[idx] *
+		          (maximum_intensity - image_magnitude[idx]);
+	}
+}
+
+
+// In-place addition of shorts. The kernel is used to add value to all elements
+// of dest. The dest array has M rows and N columns.
+static __global__ void AddShortsKernel(short* dest, const short value,
                                       const int M, const int N) {
 	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
 	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
-	if (xidx<N && yidx<M){
-		dst[yidx * N + xidx] += value;
-	}
-}
-
-
-
-
+	if (xidx < N && yidx < M){
+		dest[yidx * N + xidx] += value;
+	}
+}
+
+
+
+// TODO(nloomis): docs
 //NB: the source and destination are the same location!
 //TODO: this might be a good use for atomic operators.
-static __global__ void powerKernel(float2* fsrc, const int xpix, const int ypix,
-								   const float pow, const float offset, const float sigma){
+static __global__ void PowerKernel(float2* fsrc,
+                                   const int N, const int M,
+								   const float power_order,
+								   const float power_filter_offset,
+								   const float sigma) {
 	int xidx = threadIdx.x + blockDim.x * blockIdx.x; //pixel index that this thread deals with
 	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
-	if (xidx < xpix && yidx < ypix){
+	if (xidx < N && yidx < M){
+		int idx = yidx * N + xidx;
 
 		//get the FFTW-shifted u,v coordinates (in numbers of samples)
 		//option 4) use normalized freqs (just the magnitude) - note that uidx->u here in the nomenclature.
-		int halfx = xpix / 2;
-		int halfy = ypix / 2;
-		float u = (xidx >= halfx ? (float) (xpix - xidx) / (float) halfx :
-		                           (float) xidx / (float)halfx);
-		float v = (yidx >= halfy ? (float) (ypix - yidx) / (float)halfy :
-		                           (float)yidx / (float)halfy);
-
-		int idx = yidx * ypix + xidx; //the linear index of the data
+		int halfx = N / 2;
+		int halfy = M / 2;
+		float u = (xidx >= halfx ? (float) (N - xidx) / (float) halfx :
+		                           (float) xidx / (float) halfx);
+		float v = (yidx >= halfy ? (float) (M - yidx) / (float) halfy :
+		                           (float)yidx / (float) halfy);
+
 		float2 value = fsrc[idx]; //the value to work with
 
 		//do some operation on value using u,v
 		//here, i'm using a power filter to keep the lower freqs but kill the high freqs.
-		float fkeep = __expf( - ( __powf(u,pow) + __powf(v,pow) ) / sigma);
+		float fkeep = __expf(-( __powf(u, power_order) +
+		                        __powf(v, power_order)) /
+		                      sigma);
 		value.x = value.x * fkeep;
 		value.y = value.y * fkeep; 
-		//the offset is the fraction to keep no matter what; only (1-offset) of the value
+		//the power_filter_offset is the fraction to keep no matter what; only (1-power_filter_offset) of the value
 		//is allowed to be changed.
+
+		// TODO(nloomis): where is power_filter_offset used?!?
 
 		//store the result back to the destination
 		fsrc[idx] = value;
 	}
 }
 
-static __global__ void fillSfiltxKernel(float2* sfiltx, const float sigma, 
+// Computes the convolution kernel for a steerable filter in the x-direction. A
+// steerable filter has two important characteristics:
+//  1) the filter includes some smoothing, so the effects of random noise are
+//     reduced, and
+//  2) the gradient magnitude is independent of the orientation of an edge in
+//     an image.
+// 
+// TODO(nloomis): finish docs.
+// note: use float2 so that fft2 can be taken immediately
+static __global__ void fillSfiltxKernel(float2* gradient_filter_in_x, const float sigma, 
+                                        const int N, const int M) {
+	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
+	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
+	if (xidx < N && yidx < M){
+        int idx = yidx * N + xidx;
+        
+		int halfx = N / 2;
+		int halfy = M / 2;
+		float u = (xidx >= halfx ? -(float)(N - xidx) : (float) xidx);
+		float v = (yidx >= halfy ? -(float)(M - yidx) : (float) yidx);
+		
+		//set the value
+		gradient_filter_in_x[idx].x = -2.0f * u * __expf(-(u * u + v * v) /
+		                             (2.0f * sigma * sigma));
+        // The imaginary component is zero.
+        gradient_filter_in_x[idx].y = 0;
+
+	}
+}
+
+
+// Computes the convolution kernel for a steerable filter in the y-direction.
+// See documentation for fillSfiltxKernel above.
+//
+static __global__ void fillSfiltyKernel(float2* gradient_filter_in_y, const float sigma, 
                                         const int xpix, const int ypix) {
 	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
 	int yidx = threadIdx.y + blockDim.y * blockIdx.y;

		float oldreal = max_float[idx];
		max_metric[idx] = newreal > oldreal ? newreal : oldreal;
		index_of_max[idx] = newreal > oldreal ?
			                (short) comparison_index : index_of_max[idx];
	}
}


//finds the larger of the two; this is the overloaded version that stores the intensity at the max sim as well.
// Compares a 
static __global__ void compareRealKernel(float* max_metric, short* index_of_max_metric, float* slice_metric, float* intensity_at_max_metric, float2* slice_intensity,
										  const int slice_index, const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x*blockIdx.x;
	int yidx = threadIdx.y + blockDim.y*blockIdx.y;
	if (xidx<N && yidx<M){
		int idx = yidx*N+xidx;
		float newreal = slice_metric[idx];
		float oldreal = max_metric[idx];
		max_metric[idx] = (newreal > oldreal ? newreal : oldreal);
		index_of_max_metric[idx] = (newreal > oldreal ? (short) slice_index : index_of_max_metric[idx]);
		// If the focus metric is higher at this slice, record the intensity
		// (or, as a proxy, the real component in the x field) of this slice's
		// reconstruction; if not, keep the value the same.
		intensity_at_max_metric[idx] = (newreal > oldreal ?
		    slice_intensity[idx].x :
		    intensity_at_max_metric[idx]);
	}
}




//finds the smaller of the two
static __global__ void compareReal2Kernel(float* rintdev, short* ridxdev, float2* recdev,
										  const int i, const int M, const int N)
{
	int xidx = threadIdx.x + blockDim.x*blockIdx.x;
	int yidx = threadIdx.y + blockDim.y*blockIdx.y;
	if (xidx<N && yidx<M){
		int idx = yidx*N+xidx;
		float newreal = recdev[idx].x;
		float oldreal = rintdev[idx];
		rintdev[idx] = (newreal < oldreal ? newreal : oldreal);
		ridxdev[idx] = (newreal < oldreal ? (short)i : ridxdev[idx]);
	}
}



// In-place computation of the magnitude of a complex number. The arrray of
// complex numbers, src, has M rows and N columns. The magnitude of src is
// stored back to the real part of src, and the imaginary part is set to zero.
static __global__ void computeMagKernel(float2* src, const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M) {
		int idx = yidx * N + xidx;
		src[idx].x = complexMagnitude(src[idx]);
		src[idx].y = 0;
	}
}




static __global__ void complexMultKernel(float2* dst,
                                         float2* src1, float2* src2,
                                         const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M) {
		int idx = yidx * N + xidx;
		dst[idx] = complexMult32(src1[idx], src2[idx]);
	}
}



// sfiltmag: magnitude of a steerable filter; output
// sx: x-direction response of image convolved with steerable filter
// sy: y-direction response of image convolved with steerable filter
// M: number of rows in sfiltmag, sx, and sy
// N: number of cols in sfiltmag, sx, and sy
static __global__ void sfiltMagKernel(float* sfiltmag, float2* sx, float2* sy, 
	                                  const int M, const int N){
	int xidx = threadIdx.x + blockDim.x*blockIdx.x;
	int yidx = threadIdx.y + blockDim.y*blockIdx.y;
	if (xidx<N && yidx<M) {
		int idx = yidx*N + xidx;
		float sxmag = complexMagnitude(sx[idx]);
		float symag = complexMagnitude(sy[idx]);
		sfiltmag[idx] = sqrtf( sxmag * sxmag + symag * symag );
	}
}	


//     // occur otherwise.
//     if (xidx < N && yidx < M) {
//       int idx = yidx * N + xidx;
//       // Compute something using data from source[idx] and storing back to
//       // dest[idx].
//     }
//   }
//


// Converts a float to a complex float, storing the result into dest. The real
// part of dest is copied from source, while the imaginary part is set to zero. The
// size of the source array is M rows and N columns.
static __global__ void float_to_complexKernel(float2* dest, float* source,
                                              const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M){
		// Index of the data in the array.
		int idx = yidx * N + xidx;
		// dest.x is the real component.
		dest[idx].x = source[idx];
		// des.y is the imaginary component.
		dest[idx].y = 0;
	}
}


// In-place multiplication. The dest is a pointer to an array with size M rows
// and N columns; each value in dest is multipled by the constant scalar, sc.
static __global__ void scaleFloatKernel(float* dest, const float sc,
                                        const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M){
		int idx = yidx * N + xidx;
		dest[idx] = dest[idx] * sc;
	}
}


// Copies the real part of the value in csrc and stores it to dest. Both csrc
// and dest have M rows and N columns.
static __global__ void copyRealKernel(float* dest, float2* csrc,
                                      const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M){
		int idx = yidx * N + xidx;
		dest[idx] = csrc[idx].x;
	}
}

// Computes the magnitude of the complex-valued csrc pixel and stores it to dest.
// Both csrc and dest have M rows and N columns.
static __global__ void copyMagnitudeKernel(float* dest, float2* csrc,
                                           const int M, const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M){
		int idx = yidx * N + xidx;
		dest[idx] = complexMagnitude(csrc[idx]);
	}
}

