 	}
 }
 
+// parallel averaging
 __global__ void dev_ave(float *i1data, float *i2data, float* odata, int len) {
 	const int     tid = IMUL(blockDim.x, blockIdx.x) + threadIdx.x;
 	const int threadN = IMUL(blockDim.x, gridDim.x);

		cutilSafeCall(cudaFree(d_v1GausBuf));
		d_v1GausBuf = tmp;

		for (int orderT=0; orderT<=3; orderT++) {
			cutilSafeCall(cudaMemcpy(diffV1GausBufT, d_v1GausBuf, sizeof(float)*nrX*nrY*v1GausSize, cudaMemcpyDeviceToDevice));
			if (orderT > 0) diffV1GausBufT = diff(diffV1GausBufT, sizes, orderT,2);

			for (int orderY=0; orderY<=3-orderT; orderY++) {
				int orderX = 3-orderY-orderT;

	}
}

// parallel subtraction
__global__ void dev_sub(float *i1data, float *i2data, float* odata, int len) {
	const int     tid = IMUL(blockDim.x, blockIdx.x) + threadIdx.x;
	const int threadN = IMUL(blockDim.x, gridDim.x);
