 				//set output to zeros.
 				checkCudaErrors(cudaMemset(d_out, 0, output.totalSize()));
 				timer.Start();
-				//perform non maximum suppression. 
-				nonMaximumSuppression << <grid_size, block_size >> >(d_gradient, d_theta, d_nMax_out, numRows, numCols);
+				//perform non maximum suppression using second derivative. Need to 
+				//find zero crossings in the resulting dConvOut
+				convolve << <grid_size, block_size >> >(d_in, dConvOut, d_kernel, numRows, numCols, blurSize);
 				timer.Stop();
 				float nonMaxMs = timer.Elapsed();
 #if PRINT_INFO
-				printf("Non-maximum kernel took %f ms\n", nonMaxMs);
+				printf("Convolution kernel took %f ms\n", nonMaxMs);
 #endif
 				cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
+
+				//copy data back to convOut
+				checkCudaErrors(cudaMemcpy(convOut, dConvOut, sizeof(float)*numRows*numCols, cudaMemcpyDeviceToHost));
+				//TODO: Find zero crossings.
+				
 				uint16_t max = std::numeric_limits<uint16_t>::max();
 				timer.Start();
 				hysteresisThresholding << <grid_size, block_size >> >(d_nMax_out, d_out, d_theta, upperThresh, lowerThresh, numRows, numCols, max);

				output.gpuAlloc();
				d_out = static_cast<uint16_t*>(output.deviceData());

				//create non maximum suppression output.
				ArrayUint16 nonMaxOut(output.rows, output.cols);
				//allocate on GPU.

	}
	//set the output value. 
	d_out[offset] = static_cast<T>(result);
}

/**
* General convolution kernel.
* @param d_in input data array.
* @param d_out output data array.
* @param kernel the kernel to convolve with the data.
* @param inRows rows in the input data and output
* @param inCols cols in the input data and output
* @param kernelSize size of the kernel, note that this must be odd. 
*/
template<typename T>
__global__ void convolve(T* d_in, float* d_out, float* kernel, int inRows, int inCols, int kernelSize) {
	assert(kernelSize % 2 == 1);
	//get row and column in blcok
	int r = threadIdx.y + blockIdx.y*blockDim.y;
	int c = threadIdx.x + blockIdx.x*blockDim.x;
	//get unique point in image by finding position in grid.
	int index = c + r*blockDim.x*gridDim.x;
	if (index < inRows*inCols) {
		float result = 0.0f;
		for (int filter_r = -kernelSize / 2; filter_r <= kernelSize / 2; ++filter_r) {
			for (int filter_c = -kernelSize / 2; filter_c <= kernelSize / 2; ++filter_c) {
				int rowCompare = r + filter_r >= 0 ? r + filter_r : 0;
				int colCompare = c + filter_c >= 0 ? c + filter_c : 0;

				int image_r = rowCompare <= static_cast<int>(inRows - 1) ? rowCompare : static_cast<int>(inRows - 1);
				int image_c = colCompare <= static_cast<int>(inCols - 1) ? colCompare : static_cast<int>(inCols - 1);

				float image_value = static_cast<float>(d_in[image_r*inCols + image_c]);
				float filter_value = kernel[(filter_r + kernelSize / 2)*kernelSize + filter_c + kernelSize / 2];
				result += image_value*filter_value;
			}
		}
		d_out[index] = result;
	}
}


