 
 __global__ void cuMedianFilter (im_type* src, im_type* dest, hist_type * hist, int32_t rows, int32_t cols, int32_t r, int32_t medPos)
 {
-    __shared__ int32_t H[MF_HIST_SIZE];
-    histogramClearAllColmuns(hist,cols,MF_HIST_SIZE);
-    ClearDest(dest,rows,cols);
-    syncthreads();
+    __shared__ hist_type H[MF_HIST_SIZE];
+    __shared__ hist_type Hscan[MF_HIST_SIZE];
+    __shared__ im_type retval;
+
 
     for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
     	hist[j*MF_HIST_SIZE+src[j]]=r+2;
     }
     syncthreads();
 
-    for(int i=0; i< r; i++){
+    for(int i=1; i< r; i++){
     	int32_t pos=min(i,rows-1);
         for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
-//        	atomicAdd(hist+j*MF_HIST_SIZE+src[pos*cols+j],1);
         	hist[j*MF_HIST_SIZE+src[pos*cols+j]]++;
         }
         syncthreads();
-
-    }
+    }
+
 
     for(int i=0; i< rows; i++){
-        histogramClear(H, MF_HIST_SIZE);
+        histogramClear(H); syncthreads();
 
         int32_t possub=max(0,i-r-1);
         int32_t posadd=min(rows-1,i+r);
-        syncthreads();
 
         for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
-//        	atomicSub(hist+MF_HIST_SIZE*j+ src[possub*cols+j],1);
-//        	atomicAdd(hist+MF_HIST_SIZE*j+ src[posadd*cols+j],1);
-
         	hist[MF_HIST_SIZE*j+ src[possub*cols+j] ]--;
         	hist[MF_HIST_SIZE*j+ src[posadd*cols+j] ]++;
         }
         syncthreads();
 
         for(int32_t j=0;j<(2*r);j++){
-        	histogramAdd(H,hist+MF_HIST_SIZE*j,MF_HIST_SIZE);
-            syncthreads();
+        	histogramAdd(H,hist+MF_HIST_SIZE*j); syncthreads();
         }
 
         for(int32_t j=r;j<cols-r;j++){
             int32_t possub=max(j-r,0);
             int32_t posadd=min(j+r,cols-1);
-            histogramAdd(H, hist+posadd*MF_HIST_SIZE, MF_HIST_SIZE);
+            histogramAdd(H, hist+posadd*MF_HIST_SIZE); syncthreads();
+//            histogramMedian(H,MF_HIST_SIZE,medPos, &retval);
+//            histogramMedianPar(H,Hscan,MF_HIST_SIZE,medPos, 8,&retval);
+            histogramMedianPar256(H,Hscan,MF_HIST_SIZE,medPos, 8,&retval);
+
+//            return;
             syncthreads();
-            int32_t retval;
-            historgramMedian(H,MF_HIST_SIZE,medPos, &retval);
-            syncthreads();
 
             if(threadIdx.x==0){
-           	dest[i*cols+j]=retval;
-            }
-
-            histogramSub(H, hist+possub*MF_HIST_SIZE, MF_HIST_SIZE);
-            syncthreads();
-
+            	dest[i*cols+j]=retval;
+           }
+                syncthreads();
+
+            histogramSub(H, hist+possub*MF_HIST_SIZE); syncthreads();
         }
         syncthreads();
-
-    }
-
-}
-
-
-//int array
+    }
+}
+
+
+__global__ void cuMedianFilterMultiBlock (im_type* src, im_type* dest, hist_type * histPar, int32_t rows, int32_t cols, int32_t r, int32_t medPos)
+{
+    __shared__ hist_type H[MF_HIST_SIZE];
+    __shared__ hist_type Hscan[32];
+    __shared__ im_type retval;
+
+    int32_t extraRowThread=rows%gridDim.x;
+    int32_t doExtraRow=blockIdx.x<extraRowThread;
+    int32_t startRow=0, stopRow=0;
+    int32_t rowsPerBlock= rows/gridDim.x+doExtraRow;
+    if(doExtraRow){
+        startRow=rowsPerBlock*blockIdx.x;
+        stopRow=min(rows, startRow+rowsPerBlock);
+    }
+    else{
+        startRow=(rowsPerBlock+1)*extraRowThread+(rowsPerBlock)*(blockIdx.x-extraRowThread);    
+        stopRow=min(rows, startRow+rowsPerBlock);        
+    }
+
+    hist_type* hist=histPar+cols*MF_HIST_SIZE*blockIdx.x;
+    if (blockIdx.x==(gridDim.x-1))
+    	stopRow=rows;
+    syncthreads();
+    int32_t initNeeded=0, initVal, initStartRow, initStopRow;
+
+    if(blockIdx.x==0){
+    	initNeeded=1; initVal=r+2; initStartRow=1;	initStopRow=r;
+    }
+    else if (startRow<(r+2)){
+    	//initNeeded=1; initVal=r+2-startRow-1; initStartRow=1+startRow;	initStopRow=r+startRow+1;
+    	initNeeded=1; initVal=r+2-startRow; initStartRow=1;	initStopRow=r+stopRow-startRow;
+    }
+    else{
+    	initNeeded=0; initVal=0; initStartRow=startRow-(r+1);	initStopRow=r+startRow;    	
+    }
+   syncthreads();
+//   int counter=0;
+    if (initNeeded){
+		for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
+			hist[j*MF_HIST_SIZE+src[j]]=initVal;
+		}
+//		counter+=initVal;
+    }
+
+    syncthreads();
+    for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
+		for(int i=initStartRow; i<initStopRow; i++){
+			int32_t pos=min(i,rows-1);
+				hist[j*MF_HIST_SIZE+src[pos*cols+j]]++;
+			}
+	}
+    
+/*    
+	for(int i=initStartRow; i<initStopRow; i++){
+		int32_t pos=min(i,rows-1);
+		    for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
+				hist[j*MF_HIST_SIZE+src[pos*cols+j]]++;
+		      }
+//			counter++;  
+	}
+*/    syncthreads();
+//     if(threadIdx.x==0 && initNeeded)
+//      printf("%d, %d, %d, %d \n",blockIdx.x, startRow,stopRow, counter);
+    int32_t inc=blockDim.x*MF_HIST_SIZE;
+     for(int i=startRow; i< stopRow; i++){
+         histogramClear(H);
+         int32_t possub=max(0,i-r-1);
+         int32_t posadd=min(rows-1,i+r);
+	  int32_t possubMcols=possub*cols;
+	  int32_t posaddMcols=posadd*cols;
+         syncthreads();
+         int32_t histPos=threadIdx.x*MF_HIST_SIZE;
+         for (int32_t j=threadIdx.x; j<cols; j+=blockDim.x){
+//          	hist[histPos+ src[possubMcols+j] ]--;
+//          	hist[histPos+ src[posaddMcols+j] ]++;
+          	histPos+=inc;
+	  syncthreads();
+         }
+
+         
+         histogramMultipleAdd(H,hist, 2*r+1);         
+
+         syncthreads();         	
+         int32_t rowpos=i*cols;
+         int32_t cols_m_1=cols-1;
+         for(int32_t j=r;j<cols-r;j++){
+             int32_t possub=max(j-r,0);
+             int32_t posadd=min(j+1+r,cols_m_1);
+             //histogramMedianPar32Multi(H,Hscan,MF_HIST_SIZE,medPos, 8,&retval);
+             histogramMedianPar32WorkInefficient(H,Hscan,MF_HIST_SIZE,medPos, 8,&retval);
+             syncthreads();
+             if (threadIdx.x==0)
+            	 dest[rowpos+j]=retval;
+
+             //histogramAddAndSub(H, hist+posadd*MF_HIST_SIZE,hist+possub*MF_HIST_SIZE);                 
+             histogramAddAndSub(H, hist+(int)(posadd<<8),hist+(int)(possub<<8));                 
+             syncthreads();
+        }
+         syncthreads();
+     }
+
+}
+
 
 void readImage(char* filename,im_type* imread, int32_t rows, int32_t cols);
 float psnr(im_type* im1,im_type* im2, int32_t rows, int32_t cols);

}


__device__ void historgramMedian(hist_type* H, const int32_t size, const int32_t medPos, int32_t* retval){
	int32_t tx=threadIdx.x;
	if(tx==0){
		int32_t sum=0;
		*retval=size;
		for(int32_t i=0; i<size; i++){
			sum+=H[i];
			if(sum>=medPos){
				*retval=i;
				return;
			}
		}

	}
}

__device__ void ClearDest(im_type* dest,int32_t rows,int32_t cols){
	int32_t tx=threadIdx.x;
	int32_t totalPixels=rows*cols;
	for(; tx<totalPixels;tx+=blockDim.x)
		dest[tx]=0;
}

/////////////////////////////////////////
/////////////////////////////////////////

	cudaThreadSynchronize();
	cudaEventElapsedTime(&time, start, stop);


	// Copying filtered image back from the device to the host.
	CUDA(cudaMemcpy(hostDest,devDest,memBytesImage,cudaMemcpyDeviceToHost));

	// Loading reference image for comparison purposes.
	char refFileName[]="ref.csv";
	readImage(refFileName, hostRef, rows, cols);
	float psnrval = psnr(hostDest, hostRef, rows, cols);
	printf("Time: %f  (secs) \t\t PSNR: %f \n",time/1000,psnrval);


	for(int gridSize=32; gridSize<512; gridSize+=32)
	//for(int gridSize=32; gridSize<64; gridSize+=32)
	{
		// Setting CUDA kernel properties.
		dim3 gridDim; gridDim.x=gridSize;
		dim3 blockDim; blockDim.x=64;
		hist_type* devHistMulti;
		int32_t memBytesHistMulti=sizeof(hist_type)*cols*MF_HIST_SIZE*gridDim.x;
		CUDA(cudaMalloc((void**)(&devHistMulti), memBytesHistMulti));
		CUDA(cudaMemset(devDest,0,memBytesImage));
		CUDA(cudaMemset(devHistMulti,0,memBytesHistMulti));

		cudaDeviceSetCacheConfig(cudaFuncCachePreferL1);
		
		cudaEvent_t start,stop; float multiTime;
		cudaEventCreate(&start); cudaEventCreate(&stop);
		cudaEventRecord(start,0);
		cuMedianFilterMultiBlock<<<gridDim,blockDim>>>(devSrc, devDest, devHistMulti, rows, cols, r, medPos);
		cudaEventRecord(stop,0);
		cudaEventSynchronize(stop);
		cudaThreadSynchronize();
		cudaEventElapsedTime(&multiTime, start, stop);

		cudaFree(devHistMulti);

		// Copying filtered image back from the device to the host.
		CUDA(cudaMemcpy(hostDest,devDest,memBytesImage,cudaMemcpyDeviceToHost));

		float psnrval = psnr(hostDest, hostRef, rows, cols);
		printf("Time: %f  (secs) \t\t PSNR: %f \t\t Speedup %f\n",multiTime/1000,psnrval,time/multiTime);

	}

	// Deallocating host and device memory.
	cudaFree(devHist);
	cudaFree(devSrc);
	cudaFree(devDest);
