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

	Program:		Range Imaging Toolkit (RITK)
	Module:			ritkCudaRayCaster.txx
	Version:		$Revision: 2882 $
	Date:			$Date: 2012-04-02 17:00:07 +0200 (Mo, 02 Apr 2012) $

	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 "ritkCuda3DArrayImportImageContainer.h"
#include "ritkCudaPitchedMemoryImportImageContainer.h"
#include "ritkCudaRayCasterInterfaces.h"
#include "ritkCudaRayCaster.h"
#include <vnl/vnl_double_3.h>

//#define PERFORMANCE_CHECK

namespace ritk
{

//----------------------------------------------------------------------------
template<class TInputVolume, class TOutputImage>
CudaRayCaster<TInputVolume, TOutputImage>
::CudaRayCaster()
{
	// Default size
	Superclass::m_OutputSize[0] = 512;
	Superclass::m_OutputSize[1] = 512;
}


//----------------------------------------------------------------------------
template<class TInputVolume, class TOutputImage>
CudaRayCaster<TInputVolume,TOutputImage>
::~CudaRayCaster()
{
}


//----------------------------------------------------------------------------
template<class TInputVolume, class TOutputImage>
void
CudaRayCaster<TInputVolume,TOutputImage>
::GenerateOutputInformation()
{
	// Different dimensions are handles by the ray caster base class
	Superclass::GenerateOutputInformation();

	// Set the correct output container.
	typedef ritk::CudaPitchedMemoryImportImageContainer<unsigned long, typename OutputImageType::PixelType> OutputCudaContainerType;
	typename OutputCudaContainerType::Pointer OutputCudaContainer = dynamic_cast<OutputCudaContainerType*>(this->GetOutput()->GetPixelContainer());
	if ( !OutputCudaContainer )
	{
		OutputCudaContainer = OutputCudaContainerType::New();
		this->GetOutput()->SetPixelContainer(OutputCudaContainer);
	}
	OutputCudaContainer->SetContainerSize(Superclass::m_OutputSize);
}


//----------------------------------------------------------------------------
template<class TInputVolume, class TOutputImage>
void
CudaRayCaster<TInputVolume,TOutputImage>
::GenerateData()
{
	// Allocate outputs as a first step
	this->AllocateOutputs();

#if 0
	// Debug information
	std::cout << "InputVolume" << std::endl;
	std::cout << "===========" << std::endl;
	InputImageConstPointer Input = this->GetInput();
	std::cout << Input << std::endl;
	std::cout << std::endl << std::endl;

	std::cout << "OutputImage" << std::endl;
	std::cout << "===========" << std::endl;
	OutputImagePointer Output = this->GetOutput();
	std::cout << Output << std::endl;
	std::cout << std::endl << std::endl;

	if ( this->GetProjectionMatrixInput() )
	{
		std::cout << "ProjectionMatrix" << std::endl;
		std::cout << "================" << std::endl;
		std::cout << this->GetProjectionMatrixInput()->Get() << std::endl;
		std::cout << std::endl << std::endl;
	}
#endif

	// Sanity check for correct input cuda container
	typedef ritk::Cuda3DArrayImportImageContainer<unsigned long, typename InputVolumeType::PixelType> InputCudaContainerType;
	typename InputCudaContainerType::ConstPointer InputCudaContainer = dynamic_cast<const InputCudaContainerType*>(this->GetInput()->GetPixelContainer());
	if ( !InputCudaContainer )
	{
		itkExceptionMacro(<<"Invalid input image container detected. Must be an ritk::Cuda3DArrayImportImageContainer<unsigned long, float>");
	}

	// Sanity check for correct output cuda container
	typedef ritk::CudaPitchedMemoryImportImageContainer<unsigned long, typename OutputImageType::PixelType> OutputCudaContainerType;
	typename OutputCudaContainerType::Pointer OutputCudaContainer = dynamic_cast<OutputCudaContainerType*>(this->GetOutput()->GetPixelContainer());
	if ( !OutputCudaContainer )
	{
		itkExceptionMacro(<<"Invalid output image container detected. Must be an ritk::CudaPitchedMemoryImportImageContainer<unsigned long, float>");
	}


#ifdef PERFORMANCE_CHECK
	cudaEvent_t start_event, stop_event;
    cudaEventCreateWithFlags(&start_event, cudaEventBlockingSync);
    cudaEventCreateWithFlags(&stop_event, cudaEventBlockingSync);
	cudaEventRecord(start_event, 0);
	int NumRuns = 50;
	for ( int i = 0; i < NumRuns; i++ )
	{
#endif


	// Set up the volume parameters
	InputVolumeConstPointer InputImage = this->GetInput();
	ritkCudaRayCasterVolumeParameters VP;
	VP.Size.x = InputImage->GetRequestedRegion().GetSize()[0];
	VP.Size.y = InputImage->GetRequestedRegion().GetSize()[1];
	VP.Size.z = InputImage->GetRequestedRegion().GetSize()[2];
	VP.Origin.x = InputImage->GetOrigin()[0];
	VP.Origin.y = InputImage->GetOrigin()[1];
	VP.Origin.z = InputImage->GetOrigin()[2];
	VP.Spacing.x = InputImage->GetSpacing()[0];
	VP.Spacing.y = InputImage->GetSpacing()[1];
	VP.Spacing.z = InputImage->GetSpacing()[2];
	VP.Center.x = 0.5*InputImage->GetRequestedRegion().GetSize()[0]*InputImage->GetSpacing()[0];
	VP.Center.y = 0.5*InputImage->GetRequestedRegion().GetSize()[1]*InputImage->GetSpacing()[1];
	VP.Center.z = 0.5*InputImage->GetRequestedRegion().GetSize()[2]*InputImage->GetSpacing()[2];
	VP.ChannelDesc = &InputCudaContainer->GetChannelDesc();
	VP.DataArray = InputCudaContainer->GetCudaMemoryArray();

	// Set up the projection image parameters. Note that the pitch is in elements NOT bytes
	OutputImagePointer OutputImage = this->GetOutput();
	ritkCudaRayCasterProjectionImageParameters<typename OutputImageType::PixelType> IP;
	IP.Size.x = OutputImage->GetRequestedRegion().GetSize()[0];
	IP.Size.y = OutputImage->GetRequestedRegion().GetSize()[1];
	IP.Origin.x = OutputImage->GetOrigin()[0];
	IP.Origin.y = OutputImage->GetOrigin()[1];
	IP.Spacing.x = OutputImage->GetSpacing()[0];
	IP.Spacing.y = OutputImage->GetSpacing()[1];
	IP.Pitch = OutputCudaContainer->GetPitch()/sizeof(typename OutputImageType::PixelType);
	IP.DataPtr = OutputCudaContainer->GetCudaMemoryPointer();

	// Set up the projection parameters
	ritkCudaRayCasterProjectionParameters PP;
	// Inverse projection matrix
	ProjectionMatrixType P = this->GetProjectionMatrixInput()->Get();
	vnl_double_3x3 tmp;
	tmp(0,0) = P(0,0); tmp(0,1) = P(0,1); tmp(0,2) = P(0,2);
	tmp(1,0) = P(1,0); tmp(1,1) = P(1,1); tmp(1,2) = P(1,2);
	tmp(2,0) = P(2,0); tmp(2,1) = P(2,1); tmp(2,2) = P(2,2);
	vnl_svd<double> SVD(tmp);
	ProjectionMatrixInverseType IPM = SVD.pinverse();
	for ( unsigned int row = 0; row < 3; row++ )
		for ( unsigned int col = 0; col < 3; col++ )
			PP.InverseProjectionMatrix[row*3+col] = IPM(row,col);
	// Extract the translation part out of the projection matrix
	vnl_double_3 q;
	q(0) = -1*P(0,3);
	q(1) = -1*P(1,3);
	q(2) = -1*P(2,3);
	q = IPM*q;
	PP.OpticalCenter.x = q[0];
	PP.OpticalCenter.y = q[1];
	PP.OpticalCenter.z = q[2];

	// Call the kernel
	KernelInvocationHelper<typename InputVolumeType::PixelType, typename OutputImageType::PixelType>(VP,IP,PP);

	// Synchronize the host
	OutputCudaContainer->SynchronizeHost();

#ifdef PERFORMANCE_CHECK
	}
	cudaEventRecord(stop_event, 0);
    cudaEventSynchronize(stop_event);   
	float Time = 0;
    cudaEventElapsedTime(&Time, start_event, stop_event);
	std::cout << "Elapsed time for rendering: " << (Time/NumRuns) << std::endl;
#endif

	
}

} // namespace ritk
