     py[idx] = q[1] / norm;
     px[idx] = q[2] / norm;
 }
-
+// https://stackoverflow.com/questions/7548370/cuda-maximum-reduction-algorithm-not-working
+__global__ void reduceMax(int n, float *g_idata)
+{
+    extern __shared__ volatile float sdata[];
+    unsigned int tid = threadIdx.x;
+    unsigned int i = blockIdx.x*(blockDim.x) + tid;
+    unsigned int gridSize = blockDim.x*gridDim.x;
+
+    float val = g_idata[i];
+    i += gridSize;
+    while (i < n) {
+        val = fmax(g_idata[i],val);
+        i += gridSize;
+    }
+    sdata[tid] = val;
+    __syncthreads();
+
+    // This versions uses a single warp for the shared memory 
+    // reduction
+# pragma unroll
+    for(int i=(tid+WARPSIZE); ((tid<WARPSIZE)&&(i<blockDim.x)); i+=32)
+        sdata[tid] = fmax(sdata[tid], sdata[i]);
+
+    if (tid < 16) sdata[tid] = fmax(sdata[tid], sdata[tid + 16]);
+    if (tid < 8)  sdata[tid] = fmax(sdata[tid], sdata[tid + 8]);
+    if (tid < 4)  sdata[tid] = fmax(sdata[tid], sdata[tid + 4]);
+    if (tid < 2)  sdata[tid] = fmax(sdata[tid], sdata[tid + 2]);
+    if (tid == 0) d_max = fmax(sdata[tid], sdata[tid + 1]);
+}
+// https://stackoverflow.com/questions/7548370/cuda-maximum-reduction-algorithm-not-working
+__global__ void reduceMin(int n, float *g_idata)
+{
+    extern __shared__ volatile float sdata[];
+    unsigned int tid = threadIdx.x;
+    unsigned int i = blockIdx.x*(blockDim.x) + tid;
+    unsigned int gridSize = blockDim.x*gridDim.x;
+
+    float val = g_idata[i];
+    i += gridSize;
+    while (i < n) {
+        val = fmin(g_idata[i],val);
+        i += gridSize;
+    }
+    sdata[tid] = val;
+    __syncthreads();
+
+    // This versions uses a single warp for the shared memory 
+    // reduction
+# pragma unroll
+    for(int i=(tid+WARPSIZE); ((tid<WARPSIZE)&&(i<blockDim.x)); i+=WARPSIZE)
+        sdata[tid] = fmin(sdata[tid], sdata[i]);
+
+    if (tid < 16) sdata[tid] = fmin(sdata[tid], sdata[tid + 16]);
+    if (tid < 8)  sdata[tid] = fmin(sdata[tid], sdata[tid + 8]);
+    if (tid < 4)  sdata[tid] = fmin(sdata[tid], sdata[tid + 4]);
+    if (tid < 2)  sdata[tid] = fmin(sdata[tid], sdata[tid + 2]);
+    if (tid == 0) d_min = fmin(sdata[tid], sdata[tid + 1]);
+}
 
 // Main function
 void tvdenoising(const float* src, float* dst, float lambda,
                  const float* spacing, const long* image_size, int maxIter)
 {
     // Init params
-    size_t total = image_size[0] * image_size[1]  * image_size[2] ;
-    size_t mem_size = sizeof(float) * total;
+    size_t total_pixels = image_size[0] * image_size[1]  * image_size[2] ;
+    size_t mem_size = sizeof(float) * total_pixels;
 
     // Init cuda memory
     // BEFORE DOING ANYTHING: Use the proper CUDA enabled GPU: Tesla K40c or Gforce GT 740M

    cudaMalloc(&d_src, mem_size);
    cudaMemcpy(d_src, src, mem_size, cudaMemcpyHostToDevice);
    cudaCheckErrors("Memory Malloc and Memset: SRC");
    // U
    cudaMalloc(&d_u, mem_size);
    cudaMemcpy(d_u, src, mem_size, cudaMemcpyHostToDevice);
    cudaCheckErrors("Memory Malloc and Memset: U");
    // PX
    cudaMalloc(&d_px, mem_size);

        } \
} while (0)

__device__ float d_max;
__device__ float d_min;

__device__ __inline__
float divergence(const float* pz, const float* py, const float* px,
                 long z, long y, long x, long depth, long rows, long cols,
