     if(x<w && y<h)
     {
         size_t idx = x + (size_t)y*w;
+
         float a = 0.0f;
         float b = 0.0f;
         float c = 0.0f;


    for(int n1=0; n1<N1; n1++)
    {
        for(int n2=0; n2<N2; n2++)
        {
            DiamondDotProduct<<<grid, block>>>(dP, dD, IMG_WIDTH, IMG_HEIGHT);
            cudaDeviceSynchronize();

            ComputeU<<<grid, block>>>(dU, dV, dD, IMG_WIDTH, IMG_HEIGHT, theta);
            cudaDeviceSynchronize();
            
            DiamondOperator<<<grid, block>>>(dU, dDiaD, IMG_WIDTH, IMG_HEIGHT);

    cudaMemcpy(dU, fInDepth, (size_t)IMG_WIDTH*IMG_HEIGHT*sizeof(float), cudaMemcpyHostToDevice); CUDA_CHECK;
    cudaMemcpy(dV, dU, (size_t)IMG_WIDTH*IMG_HEIGHT*sizeof(float), cudaMemcpyDeviceToDevice); CUDA_CHECK;
    cudaMemset(dD, 0, (size_t)IMG_WIDTH*IMG_HEIGHT*sizeof(float)); CUDA_CHECK;
    //cudaMemset(dP, 0, (size_t)3*IMG_WIDTH*IMG_HEIGHT*sizeof(float)); CUDA_CHECK;
    cudaMemset(dDiaD, 0, (size_t)3*IMG_WIDTH*IMG_HEIGHT*sizeof(float)); CUDA_CHECK;

    // Init block and grid sizes
