/*======================================================================================

	Program:		RITK Plugins
	Module:			CUDASeparableNormalizedConvolutionImageFilter.cu
	Version:		$Revision: 898 $
	Date:			$Date: 2011-04-29 18:51:40 +0200 (Fri, 29 Apr 2011) $

	Copyright (c) 2010-2011
	Pattern Recognition Lab, 
    University Erlangen-Nuremberg, Germany
    All rights reserved.

	This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
	See RITKCopyright.txt for more information.

	http://www5.informatik.uni-erlangen.de/research/software/ritk

========================================================================================*/

// Due to MS Windows SDK and debug built
#ifdef _WIN32
	#define WIN32_LEAN_AND_MEAN
#endif

#include "ritkCudaMacros.h"
#include "CUDASeparableNormalizedConvolutionImageFilterInterfaces.h"
#include "CUDASeparableNormalizedConvolutionImageFilterKernel.h"

// ================================== //
// SeparableNormalizedConvolutionPlan //
// ================================== //

#if 0
//----------------------------------------------------------------------------
SeparableNormalizedConvolutionPlan::SeparableNormalizedConvolutionPlan()
{
	printf("Default constructor not supported...\n");
	NX = 0;
	NY = 0;
	Buffer1 = 0;
	Buffer2 = 0;
	Buffer3 = 0;
}
#endif


//----------------------------------------------------------------------------
SeparableNormalizedConvolutionPlan::SeparableNormalizedConvolutionPlan(unsigned int _NX, unsigned int _NY)
{
	// Set the size
	NX = _NX;
	NY = _NY;

	// Allocate CUDA memory
	ritkCudaSafeCall(cudaMallocPitch((void**)&FloatMask, &Pitch, NX*sizeof(float), NY));
	ritkCudaSafeCall(cudaMallocPitch((void**)&Buffer1, &Pitch, NX*sizeof(float), NY));
	ritkCudaSafeCall(cudaMallocPitch((void**)&Buffer2, &Pitch, NX*sizeof(float), NY));
	ritkCudaSafeCall(cudaMallocPitch((void**)&Buffer3, &Pitch, NX*sizeof(float), NY));

}


//----------------------------------------------------------------------------
SeparableNormalizedConvolutionPlan::~SeparableNormalizedConvolutionPlan()
{
	// Clean up
	if ( FloatMask )
		ritkCudaSafeCall(cudaFree((void*)FloatMask));
	if ( Buffer1 )
		ritkCudaSafeCall(cudaFree((void*)Buffer1));
	if ( Buffer2 )
		ritkCudaSafeCall(cudaFree((void*)Buffer2));
	if ( Buffer3 )
		ritkCudaSafeCall(cudaFree((void*)Buffer3));
}




// ======================== //
// CUDA function interfaces //
// ======================== //


// Copy the kernel data to the GPU.
//----------------------------------------------------------------------------
bool SetConvolutionKernel(float *h_Kernel, int KernelRadius)
{
	// Safety check
	if ( KernelRadius > MAX_KERNEL_RADIUS )
		return false;

	// Copy the data
    cudaMemcpyToSymbol(c_Kernel, h_Kernel, (2*KernelRadius+1) * sizeof(float));
	return true;
}


// 2D Normalized convolution that internally performs two separated convolutions
//----------------------------------------------------------------------------
template<unsigned int KERNEL_RADIUS>
void 
SeparableNormalizedConvolution(SeparableNormalizedConvolutionPlan &Plan, float *ImageDataOutput, const float *ImageDataInput, size_t ImageDataPitch, unsigned char *MaskDataOutput, const unsigned char *MaskDataInput, size_t MaskDataPitch)
{
	// ============ //
	// MASK CASTING //
	// ============ //

	// First of all, cast the uchar mask to a float mask. Note that each thread will read one uchar4 and writes 4 floats
	const dim3 DimBlock_CAST(16,16);
	const dim3 DimGrid_CAST(DivUp(Plan.NX, 4*DimBlock_CAST.x), DivUp(Plan.NY, DimBlock_CAST.y));
	UC2F_CastKernel<16,16><<<DimGrid_CAST, DimBlock_CAST>>>(Plan.NX, Plan.NY, Plan.FloatMask, Plan.Pitch/sizeof(float), (uchar4*)MaskDataInput, MaskDataPitch/sizeof(uchar4));


	// ========================= //
	// MASK-INPUT MULTIPLICATION // 
	// ========================= //

	// Set up the CUDA thread configuration for the multiplication
	const unsigned int TILE_STEPS_MULT = 4;
	const unsigned int BLOCKSIZE_X_MULT = 32;
	const unsigned int BLOCKSIZE_Y_MULT = 8;
	const dim3 DimBlock_MULT(BLOCKSIZE_X_MULT, BLOCKSIZE_Y_MULT);
	const dim3 DimGrid_MULT(DivUp(Plan.NX,TILE_STEPS_MULT*DimBlock_MULT.x), DivUp(Plan.NY,DimBlock_MULT.y));

	// Multiplication
	MultiplicationKernel<BLOCKSIZE_X_MULT, BLOCKSIZE_Y_MULT, TILE_STEPS_MULT><<<DimGrid_MULT, DimBlock_MULT>>>
	(
		Plan.NX,
		Plan.NY,
		Plan.Buffer1,
		Plan.FloatMask,
		ImageDataInput,
		Plan.Pitch/sizeof(float)
	);


	// =============== //
	// ROW CONVOLUTION //
	// =============== //

	// Set up the CUDA thread configuration for the row convolution
	const unsigned int TILE_STEPS_RC = 1;
	const unsigned int BLOCKSIZE_X_RC = 64;
	const unsigned int BLOCKSIZE_Y_RC = 2;
	const dim3 DimBlock_RC(BLOCKSIZE_X_RC, BLOCKSIZE_Y_RC);
	const dim3 DimGrid_RC(DivUp(Plan.NX,TILE_STEPS_RC*DimBlock_RC.x), DivUp(Plan.NY,DimBlock_RC.y));

	// Image/mask data convolution
	CombinedRowConvolutionKernel<BLOCKSIZE_X_RC, BLOCKSIZE_Y_RC, TILE_STEPS_RC, KERNEL_RADIUS><<<DimGrid_RC, DimBlock_RC>>>
	(
		Plan.NX,
		Plan.NY,
		Plan.Buffer2,
		Plan.Buffer3,
		Plan.Buffer1,
		Plan.FloatMask,
		Plan.Pitch/sizeof(float)
	);


	// ================== //
	// COLUMN CONVOLUTION //
	// ================== //

	// Set up the CUDA thread configuration for the column convolution
	const unsigned int TILE_STEPS_CC = 2;
	const unsigned int BLOCKSIZE_X_CC = 16;
	const unsigned int BLOCKSIZE_Y_CC = 8;
	const dim3 DimBlock_CC(BLOCKSIZE_X_CC, BLOCKSIZE_Y_CC);
	const dim3 DimGrid_CC(DivUp(Plan.NX,DimBlock_CC.x), DivUp(Plan.NY,TILE_STEPS_CC*DimBlock_CC.y));

	// Image data convolution
	ColumnConvolutionKernel<BLOCKSIZE_X_CC, BLOCKSIZE_Y_CC, TILE_STEPS_CC, KERNEL_RADIUS><<<DimGrid_CC, DimBlock_CC>>>
	(
		Plan.NX,
		Plan.NY,
		Plan.Buffer1,
		Plan.Buffer2,
		Plan.Pitch/sizeof(float)
	);

	// Mask convolution
	ColumnConvolutionKernel<BLOCKSIZE_X_CC, BLOCKSIZE_Y_CC, TILE_STEPS_CC, KERNEL_RADIUS><<<DimGrid_CC, DimBlock_CC>>>
	(
		Plan.NX,
		Plan.NY,
		Plan.Buffer2,
		Plan.Buffer3,
		Plan.Pitch/sizeof(float)
	);


	// =============================== //
	// NORMALIZATION AND MASK CREATION //
	// =============================== //

	// Division. Use the same cuda configuration as for the multiplication
	DivisionKernel<BLOCKSIZE_X_MULT, BLOCKSIZE_Y_MULT, TILE_STEPS_MULT><<<DimGrid_MULT, DimBlock_MULT>>>
	(
		Plan.NX,
		Plan.NY,
		Plan.Buffer3,
		Plan.Buffer1,
		Plan.Buffer2,
		Plan.Pitch/sizeof(float)
	);

	// ============== //
	// FINAL ASSEMBLY //
	// ============== //

	// Assembly. Use the same cuda configuration as for the multiplication
	FinalAssemblyKernel<BLOCKSIZE_X_MULT, BLOCKSIZE_Y_MULT, TILE_STEPS_MULT><<<DimGrid_MULT, DimBlock_MULT>>>
	(
		Plan.NX,
		Plan.NY,
		ImageDataOutput,
		Plan.Buffer1,
		ImageDataInput,
		Plan.Buffer3,
		Plan.FloatMask,
		Plan.Pitch/sizeof(float)
	);

	// Cast back
	F2UC_CastKernel<16,16><<<DimGrid_CAST, DimBlock_CAST>>>(Plan.NX, Plan.NY, (uchar4*)MaskDataOutput, MaskDataPitch/sizeof(uchar4), Plan.Buffer1, Plan.Pitch/sizeof(float));
}


// Separable normalized convolution. This method will resolve the KernelRadius argument to a templated version
//----------------------------------------------------------------------------
bool 
SeparableNormalizedConvolution(unsigned int KernelRadius, SeparableNormalizedConvolutionPlan &Plan, float *ImageDataOutput, const float *ImageDataInput, size_t ImageDataPitch, unsigned char *MaskDataOutput, const unsigned char *MaskDataInput, size_t MaskDataPitch)
{
	// Safety check
	if ( KernelRadius > MAX_KERNEL_RADIUS )
		return false;

	// Depending on the kernel radius
	switch ( KernelRadius )
	{
		case 1:
			SeparableNormalizedConvolution<1>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;	
		case 2:
			SeparableNormalizedConvolution<2>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 3:
			SeparableNormalizedConvolution<3>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 4:
			SeparableNormalizedConvolution<4>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 5:
			SeparableNormalizedConvolution<5>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 6:
			SeparableNormalizedConvolution<6>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 7:
			SeparableNormalizedConvolution<7>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 8:
			SeparableNormalizedConvolution<8>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 9:
			SeparableNormalizedConvolution<9>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 10:
			SeparableNormalizedConvolution<10>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 11:
			SeparableNormalizedConvolution<11>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 12:
			SeparableNormalizedConvolution<12>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 13:
			SeparableNormalizedConvolution<13>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 14:
			SeparableNormalizedConvolution<14>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		case 15:
			SeparableNormalizedConvolution<15>(Plan, ImageDataOutput, ImageDataInput, ImageDataPitch, MaskDataOutput, MaskDataInput, MaskDataPitch);
			return true;
		default:
			return false;
	}
}







