 	}
 
 	//reconstruct hologram and find mins
-	reconstruct(SIM, index_of_max_focus_metric, min_intensity_cpu, index_of_min_intensity, hologram_cpu, z, nz, M, N, wavelength, pixel_size, 
-		power_filter_order, power_filter_offset, sigma, maximum_intensity, intensity_at_best_focus_cpu, rintmode, record_intensity_at_best_focus);
+	Reconstruct(sim_focus_metric, index_of_max_focus_metric,
+                min_intensity_cpu, index_of_min_intensity,
+                hologram_cpu,
+                reconstruction_z_depths, num_z_depths,
+                M, N,
+                wavelength, pixel_size, 
+		        power_filter_order, power_filter_offset,
+                sigma,
+                maximum_intensity, intensity_at_best_focus_cpu,
+                record_min_intensity,
+                record_intensity_at_best_focus);
 
 	// If the hologram was originally passed in as a double or unsigned chars,
 	// hologram_cpu would have been malloc'd. Free it now, before exiting.

	//C is column-wise ordering (column-major)
	
	//retrieve the depths at which to do reconstructions
	z = mxGetPr(prhs[1]);
	nz = mxGetM(prhs[1]) * mxGetN(prhs[1]);
	
	//mexPrintf("M: %i, N: %i, nz: %i\n", M, N, nz);

	//get optional parameters; use defaults otherwise
	if (nrhs>=3)

	}
}

// Hint: it would have been much better to create structs to hold the parameters
// and options, reducing the number of variables passed in.
// Hint: the mode flags would have been better as enums than as bools.

void Reconstruct(float* focus_metric,
                 short* index_of_max_focus_metric,
                 float* min_intensity_cpu,
                 short* index_of_min_intensity,
                 float* hologram_cpu, 
			     double* reconstruction_z_depths,
                 const int num_z_depths,
                 const int M,
                 const int N, 
				 const double wavelength,
                 const double pixel_size,
				 const float power_filter_order,
                 const float power_filter_offset,
				 const float steerable_filter_sigma,
                 const float maximum_intensity,
				 float* intensity_at_best_focus_cpu,
                 bool record_min_intensity,
                 bool record_intensity_at_best_focus) {
	float* hologram_gpu;
    float* min_intensity_gpu;
    float* max_sim_focus_metric_gpu;
    float* gradient_magnitude_gpu;
    float* image_magnitude;
    float* sim;
    float* intensity_at_best_focus_gpu;
	short* index_of_min_intensity_gpu;
    short* index_of_max_focus_metric_gpu;
	float2* complex_hologram_gpu;
    float2* reconstruction_gpu;
    float2* gradient_filter_in_x;
    float2* gradient_filter_in_y;
    float2* gradient_response_in_x_gpu;
    float2* gradient_response_in_y_gpu;
	float adudu, cutoff2, fftscale;
	float2 offsetPhaseExp;
	float sigma;
    //...for both x and y dims in the block
    const int kBlockSize = 16;

    // A block is the number of threads which are launched at each round. A size
    // of 16x16 is used to step through the image 256 threads at a time. For
    // CUDA (as of 2010), threads are launched internally in multiples of 16,
    // so I used a multiple of 16 here for better processing. The exact size and
    // shape were not optimized.
	dim3 myblock(kBlockSize, kBlockSize);

    // A grid describes to the GPU how many times it should launch thread
    // blocks of a given size. There are enough grids so that one thread is
    // launched for every pixel. If N or M are not evenly divisible by the block
    // size, there will be more threads launched than pixels (due to the limited
    // granularity the GPU gives in its thread blocks), so the kernels all
    // include tests to make sure its thread corresponds to a pixel.
	dim3 mygrid( roundUpDiv(N, kBlockSize), roundUpDiv(M, kBlockSize) );

	// Upload the hologram to the GPU's memory.
    mexPrintf("Uploading hologram\n");
	cudaMalloc((void**) &hologram_gpu, sizeof(float) * M * N);
    cudaMemcpy(hologram_gpu,
               hologram_cpu,
               sizeof(float) * M * N,
               cudaMemcpyHostToDevice);

    // Convert the hologram from real-valued floats to complex-valued floats.
    // This is needed because we need the complex-valued FFT for computing the
    // physical wave propagation.
	cudaMalloc((void**) &complex_hologram_gpu, sizeof(float2) * M * N);
	FloatToComplexKernel<<<mygrid, myblock>>>(complex_hologram_gpu,
                                              hologram_gpu,
                                              M, N);
    // After the FFT is computed, the original source data does not need to be
    // maintained on the GPU.
	cudaFree(hologram_gpu);

	// Allocate memory for the remaining variables in the GPU.
    mexPrintf("Creating memory on GPU\n");
	cudaMalloc((void**) &max_sim_focus_metric_gpu, sizeof(float) * M * N);
	cudaMalloc((void**) &index_of_max_focus_metric_gpu, sizeof(short) * M * N);
	cudaMemset(index_of_max_focus_metric_gpu, 0, sizeof(short) * M * N);
	if (record_min_intensity){
		cudaMalloc((void**) &min_intensity_gpu, sizeof(float) * M * N);
		cudaMalloc((void**) &index_of_min_intensity_gpu, sizeof(short) * M * N);
		cudaMemset(index_of_min_intensity_gpu, 0, sizeof(short) * M * N);
	}

	// Create an FFT plan. The plan tells the GPU how to compute a Fourier
    // transform, its dimensionality, and the complexity. There is a TINY
    // amount of overhead while CUDA "plans" (optimizes) the computation, so
    // that it will work extremely fast later. The plan is meant to be re-used.
    mexPrintf("Creating FFT plan\n");
	cufftHandle plan;
    // CUFFT_C2C indicates that the transform is from complex-valued data and
    // that the results should be stored as complex-valued.
	cufftPlan2d(&plan, N, M, CUFFT_C2C);

    // Compute the 2D FFT of the hologram data, in place, yielding the spatial
    // frequencies of the image.
	mexPrintf("Executing FFT\n");
    // CUFFT_FORWARD indicates that this transform is a forward Fourier
    // transform (not an inverese FFT).
	cufftExecC2C(plan,
                 complex_hologram_gpu,
                 complex_hologram_gpu,
                 CUFFT_FORWARD);

// TODO(nloomis): explain sigma on power filter; finish docs for power filter

    // Holograms may be captured using a monochrome sensor or a sensor with a
    // Bayer filter. In the Bayer'ed holograms, the effect of the color filter
    // manifests itself as replicated spectra at the Nyquist frequencies. A
    // power filter acts as a low-pass filter with a smooth transition at a
    // controllable rate and location in the frequency domain.
	mexPrintf("Applying power filter\n");
    // Cheap trick: if power_filter_order is negative, the filter isn't applied.
	if (power_filter_order > 0) {
		// mexPrintf("Applying power kernel\n");
		sigma = -1.0f / log(power_filter_offset);
		PowerKernel<<<mygrid, myblock>>>(complex_hologram_gpu, N, M,
			power_filter_order, power_filter_offset, sigma);
	}

    // Compute the steerable filter, a gradient filter with smoothing. The
    // gradient filter is applied to each reconstructed slice.
    mexPrintf("Computing steerable filter\n");
	cudaMalloc((void**) &gradient_filter_in_x, sizeof(float2) * M * N);
	cudaMalloc((void**) &gradient_filter_in_y, sizeof(float2) * M * N);
	ComputeSteerableFilterXCoefsKernel<<<mygrid, myblock>>>(
        gradient_filter_in_x,
        steerable_filter_sigma,
        M, N);
	ComputeSteerableFilterYCoefsKernel<<<mygrid, myblock>>>(
        gradient_filter_in_y,
        steerable_filter_sigma,
        M, N);

    // Compute the filter's frequency-domain response. This makes it easy to
    // apply the filter later: instead of a slow spatial convolution, we can do
    // a point-wise multiplication in the Fourier domain and take the inverse
    // 2D FFT.
    // The forward 2D FFTs are performed in-place.
	cufftExecC2C(plan,
                 gradient_filter_in_x,
                 gradient_filter_in_x,
                 CUFFT_FORWARD);
	cufftExecC2C(plan,
                 gradient_filter_in_y,
                 gradient_filter_in_y,
                 CUFFT_FORWARD);
    // Reserve memory for the gradient response in each direction and the 
    // magnitude of the gradient response.
	cudaMalloc((void**) &gradient_response_in_x_gpu, sizeof(float2) * M * N);
	cudaMalloc((void**) &gradient_response_in_y_gpu, sizeof(float2) * M * N);
	cudaMalloc((void**) &gradient_magnitude_gpu, sizeof(float) * M * N);


	// Reserve memory for a single reconstructed slice of the volume.
	cudaMalloc( (void**)&reconstruction_gpu, sizeof(float2) * M * N);


