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

	Program:		RITK Plugins
	Module:			CUDASeparableNormalizedConvolutionImageFilter.txx
	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

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

#include "CUDASeparableNormalizedConvolutionImageFilter.h"


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::CUDASeparableNormalizedConvolutionImageFilter()
{
	// First, create our second output. We simply increment the output index from itk::ProcessObject
	m_MaskImageOutputIdx = ProcessObject::GetNumberOfOutputs();
	MaskImagePointer OutputMaskImage = static_cast<MaskImageType*>(MakeOutput(m_MaskImageOutputIdx).GetPointer());
	ProcessObject::SetNumberOfRequiredOutputs(ProcessObject::GetNumberOfOutputs()+1); 
	ProcessObject::SetNthOutput(m_MaskImageOutputIdx, OutputMaskImage.GetPointer());

	// Invalid CUDA memory by default
	m_CurrentSize[0] = 0;
	m_CurrentSize[1] = 0;
	m_InputImageDataContainer = NULL;
	m_OutputImageDataContainer = NULL;
	m_InputMaskDataContainer = NULL;
	m_OutputMaskDataContainer = NULL;

	// Default variance and kernel size
	m_Sigma = 50;
	m_KernelSize = 5;
	CudaBlockedMethodCallMacro(ComputeGaussianKernel);

	m_Plan = NULL;

	// Note that we can't set the CUDA pixel container here due to the itk::ImportImageContainer::Reserve() implementation
	//this->GetOutput()->SetPixelContainer(ritk::CudaRegularMemoryImportImageContainerF::New());
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::~CUDASeparableNormalizedConvolutionImageFilter()
{
	// CUDA memory should be cleaned up automatically, however, we need to destroy our plan
	if ( m_Plan )
		delete m_Plan;
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
typename CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>::DataObjectPointer
CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>
::MakeOutput(unsigned int idx)
{
	// Create the appropriate output
	if ( idx == m_MaskImageOutputIdx )
		return static_cast<itk::DataObject*>(MaskImageType::New().GetPointer());
	else
		return Superclass::MakeOutput(idx);
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
typename CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>::MaskImagePointer
CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>
::GetMaskImageOutput()
{
	return static_cast<MaskImageType*>(ProcessObject::GetOutput(m_MaskImageOutputIdx));
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
typename CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>::MaskImageConstPointer
CUDASeparableNormalizedConvolutionImageFilter<TInput, TMaskImage, TOutput>
::GetMaskImageOutput() const
{
	return static_cast<const MaskImageType*>(ProcessObject::GetOutput(m_MaskImageOutputIdx));
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::SetSigma(double Sigma)
{
	m_Sigma = Sigma;
	CudaBlockedMethodCallMacro(ComputeGaussianKernel);
	this->Modified();
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
double
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::GetSigma() const
{
	return m_Sigma;
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::SetKernelSize(int Size)
{
	m_KernelSize = Size;
	CudaBlockedMethodCallMacro(ComputeGaussianKernel);
	this->Modified();
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
int
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::GetKernelSize()
{
	return m_KernelSize;
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::ComputeGaussianKernel()
{
	// Allocate our container
	m_KernelContainer = KernelContainerType::New();
	KernelContainerType::SizeType Size;
	Size[0] = m_KernelSize*2+1;
	m_KernelContainer->SetContainerSize(Size);
	m_KernelContainer->Reserve(m_KernelSize*2+1);

	// Compute the gaussian
	for ( int i = 0; i < m_KernelSize*2+1; i++ )
	{
		double d1 = 1.0 / (sqrt(2.0*3.1415)*m_Sigma);
		double d2 = (i-m_KernelSize)/m_Sigma;
		double d3 = d1 * exp(-0.5*d2*d2);
		m_KernelContainer->GetBufferPointer()[i] = d3;
	}

	// Copy to the device
	SetConvolutionKernel(m_KernelContainer->GetBufferPointer(), m_KernelSize);
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::AllocateCudaMemory()
{
	// Size in bytes of the image
	typename InputImageType::SizeType Size = this->GetInput()->GetBufferedRegion().GetSize();

	// Allocate CUDA device memory
	if ( m_CurrentSize != Size )
	{
		// Create new input image container
		m_InputImageDataContainer = ImageDataContainerType::New();
		m_InputImageDataContainer->SetContainerSize(Size);
		m_InputImageDataContainer->Reserve(Size[0]*Size[1]);

		// Create new mask image container
		m_InputMaskDataContainer = MaskDataContainerType::New();
		m_InputMaskDataContainer->SetContainerSize(Size);
		m_InputMaskDataContainer->Reserve(Size[0]*Size[1]);

		// Create output container
		m_OutputImageDataContainer = ImageDataContainerType::New();
		m_OutputImageDataContainer->SetContainerSize(Size);
		m_OutputImageDataContainer->Reserve(Size[0]*Size[1]);
		this->GetOutput()->SetPixelContainer(m_OutputImageDataContainer);

		// Create mask output container
		m_OutputMaskDataContainer = MaskDataContainerType::New();
		m_OutputMaskDataContainer->SetContainerSize(Size);
		m_OutputMaskDataContainer->Reserve(Size[0]*Size[1]);
		GetMaskImageOutput()->SetPixelContainer(m_OutputMaskDataContainer);

		// Update the size member
		m_CurrentSize = Size;
	}	
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::DoFiltering()
{
	// Create a plan on demand
	if ( !m_Plan || m_Plan->NX != m_CurrentSize[0] || m_Plan->NY != m_CurrentSize[1] )
	{
		if ( m_Plan )
			delete m_Plan;
		m_Plan = new SeparableNormalizedConvolutionPlan(m_CurrentSize[0], m_CurrentSize[1]);
	}

	// Copy the data
	m_InputImageDataContainer->DeepCopy(this->GetInput()->GetPixelContainer());
	m_InputMaskDataContainer->DeepCopy(GetMaskImageInput()->GetPixelContainer());

	// CUDA invocation
	SeparableNormalizedConvolution
	(
		m_KernelSize,
		*m_Plan,
		m_OutputImageDataContainer->GetCudaMemoryPointer(),
		m_InputImageDataContainer->GetCudaMemoryPointer(),
		m_InputImageDataContainer->GetPitch(),
		m_OutputMaskDataContainer->GetCudaMemoryPointer(),
		m_InputMaskDataContainer->GetCudaMemoryPointer(),
		m_InputMaskDataContainer->GetPitch()
	);
}


//----------------------------------------------------------------------------
template<class TInput, class TMaskImage, class TOutput>
void
CUDASeparableNormalizedConvolutionImageFilter<TInput,TMaskImage,TOutput>
::GenerateData()
{
	// Sanity checks
	if ( !this->GetInput() )
		itkExceptionMacro(<<"No input provided!");
	if ( !GetMaskImageInput() )
		itkExceptionMacro(<<"No mask provided!");
	if ( this->GetInput()->GetRequestedRegion() != GetMaskImageInput()->GetRequestedRegion() )
		itkExceptionMacro(<<"Region mismatch detected: Input and mask must have the same region!");

	// Ensure that we have valid outputs
	this->GetOutput()->SetRequestedRegion(this->GetInput()->GetRequestedRegion());
	this->GetOutput()->SetBufferedRegion(this->GetInput()->GetBufferedRegion());
	this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
	this->GetMaskImageOutput()->SetRequestedRegion(this->GetMaskImageInput()->GetRequestedRegion());
	this->GetMaskImageOutput()->SetBufferedRegion(this->GetMaskImageInput()->GetBufferedRegion());
	this->GetMaskImageOutput()->SetLargestPossibleRegion(this->GetMaskImageInput()->GetLargestPossibleRegion());

	// Allocate host memory
	this->GetOutput()->Allocate();
	this->GetMaskImageOutput()->Allocate();

	// Allocate device memory
	AllocateCudaMemory();

	// Do the actual filtering
	CudaBlockedMethodCallMacro(DoFiltering);
}







