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

	Program:		Range Imaging Toolkit (RITK)
	Module:			ritkFloatVolumeCudaRayCasterKernel.h
	Version:		$Revision: 2884 $
	Date:			$Date: 2012-04-02 17:43:53 +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

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


#ifndef RITKFLOATVOLUMECUDARAYCASTERKERNEL_H__
#define RITKFLOATVOLUMECUDARAYCASTERKERNEL_H__

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

/// Our volume texture
/// ==================
texture<float,3,cudaReadModeElementType> VolumeTexture;

/**	@name	The parameters for ray casting */
//@{
__constant__ ritkCudaRayCasterVolumeParameters c_VP;
__constant__ ritkCudaRayCasterProjectionImageParameters<float> c_IP;
__constant__ ritkCudaRayCasterProjectionParameters c_PP;
//@}




/// The actual ray casting kernel
/// =============================
template<int NUMSAMPLES>
__global__ void
FloatVolumeCudaRayCasterKernel(const ritkCudaRayCasterVolumeParameters VP, const ritkCudaRayCasterProjectionImageParameters<float> 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] = 0.f;
		return;
	}

	// Initialize step length
	float3 Delta = RayDirection * (Alpha.y-Alpha.x)/(NUMSAMPLES-1);
	
	// World [0:N] to texture [0:1] conversion
	//float3 WorldToTex = make_float3(1.f)/(VP.Spacing*VP.Size);

	// 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
	float IntensityAccu = 0.f;
#pragma unroll 
	for ( int i = 0; i < NUMSAMPLES; i++ )
	{
		// Update position
		TexCoord += TexCoordDelta;

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

		// Accumulate
		IntensityAccu += VolumeValue;
	}

	// Set values in the output image
	//IntensityAccu += NUMSAMPLES*1024;
	float scale = FixedScale*(Alpha.y-Alpha.x);
	IntensityAccu *= scale;
	//IntensityAccu = fminf(IntensityAccu,0.0004f);
	//float MaxVal = 0.5f;
	//IntensityAccu = fminf(IntensityAccu,MaxVal)/MaxVal;
	IP.DataPtr[pos] = IntensityAccu;
}










#if 0
/// Kernel to compute the rays for each thread.
/// ==========================================
__global__ void
ComputeRays(const ritkCudaRayCasterVolumeParameters VP, const ritkCudaRayCasterProjectionImageParameters IP, const ritkCudaRayCasterProjectionParameters PP, float4 *Output)
{
	// 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);

	// Assemble the ray
	float4 Ray;
	Ray.x = RayDirection.x * Alpha.x;
	Ray.y = RayDirection.y * Alpha.x;
	Ray.z = RayDirection.z * Alpha.x;
	Ray.w = Alpha.y;

	// Write the output
	Output[pos] = Ray;
}


/// The actual ray casting kernel
/// =============================
template<int NUMSAMPLES>
__global__ void
RayCasterKernel(const ritkCudaRayCasterVolumeParameters VP, const ritkCudaRayCasterProjectionImageParameters IP, const ritkCudaRayCasterProjectionParameters PP, float4 *Rays, 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 index for this threaed
	int pos = IY*IP.Pitch + IX;

	// Load the ray for this thread
	float4 Ray = Rays[pos];
	float2 Alpha;
	Alpha.x = sqrtf(Ray.x*Ray.x + Ray.y*Ray.y + Ray.z*Ray.z);
	Alpha.y = Ray.w;
	float3 RayDirection;
	RayDirection.x = Ray.x / Alpha.x;
	RayDirection.y = Ray.y / Alpha.x;
	RayDirection.z = Ray.z / Alpha.x;

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

	// Initialize step length
	float3 Delta = RayDirection * (Alpha.y-Alpha.x)/(NUMSAMPLES-1);
	
	// World [0:N] to texture [0:1] conversion
	//float3 WorldToTex = make_float3(1.f)/(VP.Spacing*VP.Size);

	// 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
	const int NumInnerLoops = 4;
	float IntensityAccu = 0.f;
	for ( int i = 0; i < NUMSAMPLES/NumInnerLoops; i++ )
	{
#pragma unroll 
		for ( int j = 0; j < NumInnerLoops; j++ )
		{
			// Update position
			TexCoord += TexCoordDelta;

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

			/*if ( fabsf(VolumeValue-40) < 40 )
				VolumeValue = 0;
			else
				VolumeValue = -1024;*/

			// Accumulate
			IntensityAccu += VolumeValue;
		}
	}

	// Set values in the output image
	//IntensityAccu += NUMSAMPLES*1024;
	float scale = FixedScale*(Alpha.y-Alpha.x);
	IntensityAccu *= scale;
	IntensityAccu = fminf(IntensityAccu,0.0007f);
	IP.DataPtr[pos] = IntensityAccu;
}
#endif



#endif // RITKFLOATVOLUMECUDARAYCASTERKERNEL_H__