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

	Program:		Range Imaging Toolkit (RITK)
	Module:			ritkRGBAVolumeCudaRayCasterKernel.h
	Version:		$Revision: 2871 $
	Date:			$Date: 2012-03-31 17:42:14 +0200 (Sa, 31 Mrz 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

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


#ifndef RITKRGBAVOLUMECUDARAYCASTERKERNEL_H__
#define RITKRGBAVOLUMECUDARAYCASTERKERNEL_H__

#include "ritkCudaRayCasterInterfaces.h"
#include "ritkCudaRayCasterCommon.h"

//#define USE_JITTER

/// Our volume texture
/// ==================
texture<uchar4,3,cudaReadModeElementType> VolumeTexture_uchar;
texture<float4,3,cudaReadModeElementType> VolumeTexture_float;


/// The actual ray casting kernel
/// =============================
template<int NUMSAMPLES>
__global__ void
RGBAVolumeCudaRayCasterKernel_uchar(const ritkCudaRayCasterVolumeParameters VP, const ritkCudaRayCasterProjectionImageParameters<uchar4> IP, const ritkCudaRayCasterProjectionParameters PP, float3 WorldToTex, float FixedScale)
{
	// Get the output image indices for this thread
	int IX = blockIdx.x*blockDim.x + threadIdx.x;
	int IY = blockIdx.y*blockDim.y + threadIdx.y;

	// Check for out-of-bounds
	if ( IX >= IP.Size.x || IY >= IP.Size.y )
		return;

	// The linear output index for this threaed
	int pos = IY*IP.Pitch + IX;

#ifdef USE_JITTER
	// Multi-ray cast
	int Steps = 1;
	float delta = 0.5f;
	float dx = delta * IP.Spacing.x;
	float dy = delta * IP.Spacing.y;
	float4 Accu = make_float4(0,0,0,0);
	for ( int jx = -Steps; jx <= Steps; jx++ )
	{
		for ( int jy = -Steps; jy <= Steps; jy++ )
		{
			// Compute physical (aka continous) image indices from the pixel indices for this thread
			float PX = (IP.Origin.x/IP.Spacing.x) + IX + jx*dx;
			float PY = (IP.Origin.y/IP.Spacing.y) + IY + jy*dy;

			// Set up the ray for this position in the image plan
			float3 RayDirection;
			ComputeRayDirection(PX,PY,PP.InverseProjectionMatrix,&RayDirection);

			// Compute intersection with the volume (results are stored in the ray)
			float2 Alpha;
			IntersectVolume(&PP.OpticalCenter, &RayDirection, &VP, &Alpha); 

			// Sanity check
			if ( Alpha.y <= Alpha.x || fabsf(Alpha.y-Alpha.x) < 1e-3f )
			{
				IP.DataPtr[pos] = make_uchar4(0,0,0,0);
				return;
			}

			// Initialize step length
			float3 Delta = RayDirection * (Alpha.y-Alpha.x)/(NUMSAMPLES-1);

			// Texture position (SamplePosition-VolumeOrigin)*WorldToTex and incremental StepLength*WorldToTex
			float3 TexCoord = (PP.OpticalCenter + Alpha.x * RayDirection - VP.Origin) * WorldToTex;
			float3 TexCoordDelta = Delta*WorldToTex;

			// Evaluate the ray
			float4 IntensityAccu = make_float4(0.f,0.f,0.f,0.f);
			for ( int i = 0; i < NUMSAMPLES; i++ )
			{
				// Update position
				TexCoord += TexCoordDelta;

				// Get value
				uchar4 VolumeValue = tex3D(VolumeTexture,TexCoord.x,TexCoord.y,TexCoord.z);

				// Accumulate
				IntensityAccu.x += VolumeValue.x;
				IntensityAccu.y += VolumeValue.y;
				IntensityAccu.z += VolumeValue.z;
				IntensityAccu.w += VolumeValue.w;

				// Termination
				if ( IntensityAccu.x > 1e-3f )
					break;
				if ( IntensityAccu.y > 1e-3f )
					break;
				if ( IntensityAccu.z > 1e-3f )
					break;
				if ( IntensityAccu.w > 1e-3f )
					break;
			}

			// Accumulate the multi-ray accu
			Accu += IntensityAccu;
		}
	}

	// Normalize multi-ray accu
	Accu /= (Steps*2+1)*(Steps*2+1);

	// Set values in the output image
	uchar4 Result;
	Result.x = Accu.x;
	Result.y = Accu.y;
	Result.z = Accu.z;
	Result.w = Accu.w;
	IP.DataPtr[pos] = Result;
#else
	// Compute physical (aka continous) image indices from the pixel indices for this thread
	float PX = (IP.Origin.x/IP.Spacing.x) + IX;
	float PY = (IP.Origin.y/IP.Spacing.y) + IY;

	// Set up the ray for this position in the image plan
	float3 RayDirection;
	ComputeRayDirection(PX,PY,PP.InverseProjectionMatrix,&RayDirection);

	// Compute intersection with the volume (results are stored in the ray)
	float2 Alpha;
	IntersectVolume(&PP.OpticalCenter, &RayDirection, &VP, &Alpha); 

	// Sanity check
	if ( Alpha.y <= Alpha.x || fabsf(Alpha.y-Alpha.x) < 1e-3f )
	{
		IP.DataPtr[pos] = make_uchar4(0,0,0,0);
		return;
	}

	// Initialize step length
	float3 Delta = RayDirection * (Alpha.y-Alpha.x)/(NUMSAMPLES-1);
			
	// Texture position (SamplePosition-VolumeOrigin)*WorldToTex and incremental StepLength*WorldToTex
	float3 TexCoord = (PP.OpticalCenter + Alpha.x * RayDirection - VP.Origin) * WorldToTex;
	float3 TexCoordDelta = Delta*WorldToTex;

	// Evaluate the ray
	float4 IntensityAccu = make_float4(0.f,0.f,0.f,0.f);
	for ( int i = 0; i < NUMSAMPLES; i++ )
	{
		// Update position
		TexCoord += TexCoordDelta;

		// Get value
		uchar4 VolumeValue = tex3D(VolumeTexture_uchar,TexCoord.x,TexCoord.y,TexCoord.z);

		// Accumulate
		IntensityAccu.x += VolumeValue.x;
		IntensityAccu.y += VolumeValue.y;
		IntensityAccu.z += VolumeValue.z;
		IntensityAccu.w += VolumeValue.w;

		// Termination
		if ( IntensityAccu.x > 1e-3f )
			break;
		if ( IntensityAccu.y > 1e-3f )
			break;
		if ( IntensityAccu.z > 1e-3f )
			break;
		if ( IntensityAccu.w > 1e-3f )
			break;
	}

	// Set values in the output image
	uchar4 Result;
	Result.x = IntensityAccu.x;
	Result.y = IntensityAccu.y;
	Result.z = IntensityAccu.z;
	Result.w = IntensityAccu.w;
	IP.DataPtr[pos] = Result;
#endif
}


/// The actual ray casting kernel
/// =============================
template<int NUMSAMPLES>
__global__ void
RGBAVolumeCudaRayCasterKernel_float(const ritkCudaRayCasterVolumeParameters VP, const ritkCudaRayCasterProjectionImageParameters<float4> IP, const ritkCudaRayCasterProjectionParameters PP, float3 WorldToTex, float FixedScale)
{
	// Get the output image indices for this thread
	int IX = blockIdx.x*blockDim.x + threadIdx.x;
	int IY = blockIdx.y*blockDim.y + threadIdx.y;

	// Check for out-of-bounds
	if ( IX >= IP.Size.x || IY >= IP.Size.y )
		return;

	// The linear output index for this threaed
	int pos = IY*IP.Pitch + IX;

	// Compute physical (aka continous) image indices from the pixel indices for this thread
	float PX = (IP.Origin.x/IP.Spacing.x) + IX;
	float PY = (IP.Origin.y/IP.Spacing.y) + IY;

	// Set up the ray for this position in the image plan
	float3 RayDirection;
	ComputeRayDirection(PX,PY,PP.InverseProjectionMatrix,&RayDirection);

	// Compute intersection with the volume (results are stored in the ray)
	float2 Alpha;
	IntersectVolume(&PP.OpticalCenter, &RayDirection, &VP, &Alpha); 

	// Sanity check
	if ( Alpha.y <= Alpha.x || fabsf(Alpha.y-Alpha.x) < 1e-3f )
	{
		IP.DataPtr[pos] = make_float4(0,0,0,0);
		return;
	}

	// Initialize step length
	float3 Delta = RayDirection * (Alpha.y-Alpha.x)/(NUMSAMPLES-1);
			
	// Texture position (SamplePosition-VolumeOrigin)*WorldToTex and incremental StepLength*WorldToTex
	float3 TexCoord = (PP.OpticalCenter + Alpha.x * RayDirection - VP.Origin) * WorldToTex;
	float3 TexCoordDelta = Delta*WorldToTex;

	// Evaluate the ray
	float4 RGBAccu = make_float4(0.f,0.f,0.f,0.f);
	float RemainingOpacity = 1.5f;
	float Weights = 0.f;

	for ( int i = 0; i < NUMSAMPLES; i++ )
	{
		// Update position
		TexCoord += TexCoordDelta;

		// Get the volume value
		float4 VolumeValue = tex3D(VolumeTexture_float,TexCoord.x,TexCoord.y,TexCoord.z);

		// Compute the effective opacity according to previous attenuations
		float Opacity = (VolumeValue.w > 32.f) ? VolumeValue.w : 0.0f;
		Opacity = fminf(Opacity/255.f, RemainingOpacity);

		// Update the RGB and weights accus and update the remaining opacity
		RGBAccu += Opacity*VolumeValue;
		Weights += Opacity;
		RemainingOpacity -= Opacity;

		// Termination
		if ( fabsf(RemainingOpacity) < 1e-6f )
			break;
	}


	// Set values in the output image
	float4 Result;
	float Scale = (Weights > 1e-6f) ? 1.f/Weights : 0.0f;
	Scale /= 255.f;
	Result.x = RGBAccu.x*Scale;
	Result.y = RGBAccu.y*Scale;
	Result.z = RGBAccu.z*Scale;
	Result.w = RGBAccu.w*Scale;
	IP.DataPtr[pos] = Result;
}


#endif // RITKRGBAVOLUMECUDARAYCASTERKERNEL_H__