     
 // Main function
     void tvdenoising(const float* src, float* dst, float lambda,
-            const float* spacing, const long* image_size, int maxIter)
-    {
+            const float* spacing, const long* image_size, int maxIter){
         
         // Prepare for MultiGPU
         int deviceCount = 0;

            // As mem_auxiliary is not expected to be a large value (for a 2000^3 image is around 28Mbytes), lets for now assume we need it all
            size_t mem_free=mem_GPU_global-mem_auxiliary;
            
            splits=(unsigned int)(ceil(((float)(3*mem_size_image)/(float)(deviceCount))/mem_free));
            // Now, there is an overhead here, as each splits should have 2 slices more, to accoutn for overlap of images.
            // lets make sure these 2 slices fit, if they do not, add 1 to splits.
            slices_per_split=(image_size[2]+deviceCount*splits-1)/(deviceCount*splits);
            mem_img_each_GPU=(mem_slice_image*(slices_per_split+buffer_length*2));
            
            // if the new stuff does not fit in the GPU, it measn we are in the edge case where adding that extra slice will overflow memory
            if (mem_GPU_global< 3*mem_img_each_GPU+mem_auxiliary){
                // one more splot shoudl do the job, as its an edge case.
                splits++;
                //recompute for later
                slices_per_split=(image_size[2]+deviceCount*splits-1)/(deviceCount*splits); // amountf of slices that fit on a GPU. Later we add 2 to these, as we need them for overlap
                mem_img_each_GPU=(mem_slice_image*(slices_per_split+buffer_length*2));
            }

            // How many EXTRA buffer slices shoudl be able to fit in here??!?!
            mem_free=mem_GPU_global-(3*mem_img_each_GPU+mem_auxiliary);
            unsigned int extra_buff=(mem_free/mem_slice_image); 
            buffer_length=(extra_buff/2)/3; // we need double whatever this results in, rounded down.
            mem_img_each_GPU=(mem_slice_image*(slices_per_split+buffer_length*2));

            // Assert
            if (mem_GPU_global< 3*mem_img_each_GPU+mem_auxiliary){
                mexErrMsgIdAndTxt("minimizeTV:POCS_TV:GPU","Bad assert. Logic behind spliting flawed! Please tell: ander.biguri@gmail.com\n");
            }
        }
        
        
        
        
        // Init params
        size_t total_pixels = image_size[0] * image_size[1]  * image_size[2] ;
        size_t mem_size = sizeof(float) * total_pixels;
        
        
        
        float *d_src, *d_u, *d_px, *d_py, *d_pz;
        
        // F
        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, d_src, mem_size, cudaMemcpyDeviceToDevice);
        cudaCheckErrors("Memory Malloc and Memset: U");
        // PX
        cudaMalloc(&d_px, mem_size);
        cudaMemset(d_px, 0, mem_size);
        cudaCheckErrors("Memory Malloc and Memset: PX");
        // PY
        cudaMalloc(&d_py, mem_size);
        cudaMemset(d_py, 0, mem_size);
        cudaCheckErrors("Memory Malloc and Memset: PY");
        // PZ
        cudaMalloc(&d_pz, mem_size);
        cudaMemset(d_pz, 0, mem_size);
        cudaCheckErrors("Memory Malloc and Memset: PZ");
        
        // bdim and gdim
        dim3 block(10, 10, 10);
        dim3 grid((image_size[0]+block.x-1)/block.x, (image_size[1]+block.y-1)/block.y, (image_size[2]+block.z-1)/block.z);
        
        int i = 0;
        
        float tau2, tau1;
        
        for ( i = 0; i < maxIter; i++ )
        {
            tau2 = 0.3f + 0.02f * i;
            tau1 = (1.f/tau2) * ((1.f/6.f) - (5.f/(15.f+i)));
            
            update_u<<<grid, block>>>(d_src, d_pz, d_py, d_px, d_u, tau1, lambda,
                    image_size[2], image_size[1],image_size[0],
                    spacing[2], spacing[1], spacing[0]);
            
            update_p<<<grid, block>>>(d_u, d_pz, d_py, d_px, tau2,
                    image_size[2], image_size[1], image_size[0],
                    spacing[2], spacing[1], spacing[0]);
            
        }
        
        cudaCheckErrors("TV minimization");
        

        // Decide how are we handling the distribution of computation
        size_t mem_img_each_GPU;
        
        unsigned int buffer_length=1;
        //Does everything fit in the GPU?
        unsigned int slices_per_split;
        unsigned int splits=1; // if the number does not fit in an uint, you have more serious trouble than this.
        if(mem_GPU_global> 5*mem_size_image+5*mem_slice_image*2){
            // We only need to split if we have extra GPUs
            slices_per_split=(image_size[2]+deviceCount-1)/deviceCount;
            mem_img_each_GPU=mem_slice_image*((image_size[2]+buffer_length*2+deviceCount-1)/deviceCount);
