#ifdef __DEVICE_EMULATION__
#include <stdio.h>
#endif

#ifndef _TEXTURE3D_KERNEL_H_
#define _TEXTURE3D_KERNEL_H_

#include "cutil_math.h"

#define tex3Dv(t,p) tex3D(t, p.x, p.y, p.z)

// a bit of matrix algebra
typedef struct {
    float4 m[3];
} float3x4;
__device__ float3 operator*(const float3x4 &M, const float3 &v) { return make_float3(dot(v, make_float3(M.m[0])), dot(v, make_float3(M.m[1])), dot(v, make_float3(M.m[2]))); }
__device__ float4 operator*(const float3x4 &M, const float4 &v) { return make_float4(dot(v, M.m[0]), dot(v, M.m[1]), dot(v, M.m[2]), 1.0f); }

struct Ray {
	float3 o;	// origin
	float3 d;	// direction
};

// intersect ray with a box
// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm
__device__ int intersectBox(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    float3 invR = make_float3(1.0f) / r.d;
    float3 tbot = invR * (boxmin - r.o);
    float3 ttop = invR * (boxmax - r.o);

    // re-order intersections to find smallest and largest on each axis
    float3 tmin = fminf(ttop, tbot);
    float3 tmax = fmaxf(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
    float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));

    *tnear = largest_tmin;
    *tfar = smallest_tmax;

    return smallest_tmax > largest_tmin;
}

// color conversion functions
__device__ uint rgbaFloatToInt(float4 rgba) { return (uint(__saturatef(rgba.w)*255)<<24) | (uint(__saturatef(rgba.z)*255)<<16) | (uint(__saturatef(rgba.y)*255)<<8) | uint(__saturatef(rgba.x)*255); }
__device__ uint rFloatToInt(float r)
{
    r = __saturatef(r);
    return (uint(r*255)<<24) | (uint(r*255)<<16) | (uint(r*255)<<8) | uint(r*255);
}

////////////////////////////////////////////////////////////////////////////////
// constants

#define RENDER_AMBIENT   1
#define RENDER_DIFFUSE   2
#define RENDER_SPECULAR  4

// the inverse view matrix (needed for raytracing)
__constant__ float3x4 dInvViewMatrix;

// these are set by host code before the kernels are run
__constant__ int dFieldDimension; // dimension of vector / scalar fields
__constant__ float3 dLightPos; // light position
__constant__ int dRenderFlags; // parameters
__constant__ float dShiny = 7.0f; // for specular lighting
__constant__ float dIso;

// "pressure" texture
texture<float, 3, cudaReadModeElementType> pressureTex;

////////////////////////////////////////////////////////////////////////////////
// spectrum function (from my lab 7, and from some other code of mine before that)

__device__ const float cycler[] = {0.0f,1.0f,1.0f,1.0f,0.0f,0.0f,0.0f};
__device__ const float cycleg[] = {0.0f,0.0f,0.5f,1.0f,1.0f,1.0f,0.0f};
__device__ const float cycleb[] = {0.0f,0.0f,0.0f,0.0f,0.0f,1.0f,1.0f};

__device__ void InterpolateSpectrum(float position, float* rvalue, float* gvalue, float* bvalue)
{
    position *= 6;
    float positionFraction = position - floor(position);
    int index = (int)floor(position);
    *rvalue = cycler[index] * (1.0f - positionFraction) + cycler[index + 1] * positionFraction;
    *gvalue = cycleg[index] * (1.0f - positionFraction) + cycleg[index + 1] * positionFraction;
    *bvalue = cycleb[index] * (1.0f - positionFraction) + cycleb[index + 1] * positionFraction;
}

////////////////////////////////////////////////////////////////////////////////
// display kernel

__device__ float3 PressureTexNormal(float3 pos)
{
    float3 normal;
    float d = 0.5f / dFieldDimension;
    normal.x = (tex3D(pressureTex, pos.x+d, pos.y, pos.z) - tex3D(pressureTex, pos.x-d, pos.y, pos.z));
    normal.y = (tex3D(pressureTex, pos.x, pos.y+d, pos.z) - tex3D(pressureTex, pos.x, pos.y-d, pos.z));
    normal.z = (tex3D(pressureTex, pos.x, pos.y, pos.z+d) - tex3D(pressureTex, pos.x, pos.y, pos.z-d));
    return normalize(normal);
}

__global__ void RenderPressureIsosurfaceK(unsigned long* output, unsigned long imageW, unsigned long imageH)
{
    // amount to step by
    float tstep = 0.003f;
    int maxSteps = 2000;

    unsigned long x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned long y = blockIdx.y * blockDim.y + threadIdx.y;

    float u = (x / (float) imageW) * 2.0f - 1.0f;
    float v = (y / (float) imageH) * 2.0f - 1.0f;

    // the index of the pixel to be written by this call
    long outputColorIndex = y * imageW + x;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(dInvViewMatrix * make_float4(0.0f, 0.0f, 0.0f, 1.0f));
    eyeRay.d = dInvViewMatrix * normalize(make_float3(u, v, -2.0f));

    // find intersection with box
    float tnear, tfar;

    // return if not intersecting
    if (!intersectBox(eyeRay, make_float3(-2.0f,-2.0f,-2.0f), make_float3(2.0f,2.0f,2.0f), &tnear, &tfar))
    {
        output[outputColorIndex] = 0;
        return;
    }

    // clamp to near plane
    if (tnear < 0.0f) tnear = 0.0f;

    // start stepping through space
    float t = tnear;
    float3 pos, texPos;
    for (int i = 0; (i < maxSteps) && (t <= tfar); i++)
    {
        texPos = (eyeRay.o + eyeRay.d * t) * 0.25f + 0.5f;
        if (tex3Dv(pressureTex, texPos) > dIso) break;
        t += tstep;
    }

    // recalculate pos in world space (not texture space)
    pos = (eyeRay.o + eyeRay.d * t);

    // this pixel isn't part of the surface? then it's black
    if (t > tfar)
    {
        output[outputColorIndex] = 0;
        return;
    }

    // get the normal vector
    float3 normal = PressureTexNormal(pos * 0.25f + 0.5f);
    if ((x < imageW) && (y < imageH))
    {
        // calculate the color based on whatever
        float4 color, objectColor;

        switch (dRenderFlags & 0x18)
        {
          case 0x00:
            objectColor = make_float4(cos(2 * pos.x),cos(2 * pos.y),cos(2 * pos.z),0.0f); // position spectrum
            objectColor *= objectColor;
            break;
          case 0x08:
            objectColor = make_float4((normal + 1.0f) / 2, 0.0f); // normals
            break;
          case 0x10:
            objectColor = make_float4(1.0f,0.0f,0.0f,0.0f); // red
            break;
          case 0x18:
            InterpolateSpectrum(length(pos) / 1.733f, &objectColor.x, &objectColor.y, &objectColor.z);
        }

        // if nothing else, it's black
        color = make_float4(0.0f, 0.0f, 0.0f, 1.0f);

        // ambient lighting...
        if (dRenderFlags & RENDER_AMBIENT)
            color = objectColor * 0.05f;

        // diffuse lighting...
        if (dRenderFlags & RENDER_DIFFUSE)
        {
            float NdotL = max(dot(normal, normalize(pos - dLightPos)), 0.0f);
            if (NdotL > 0.0f) color += objectColor * NdotL;
        }

        // specular lighting...
        if (dRenderFlags & RENDER_SPECULAR)
        {
            float RdotEye = max(dot(normalize(reflect(dLightPos - pos, normal)), normalize(eyeRay.d)), 0.0f);
            if (RdotEye > 0.0f)
            {
                float4 specColor = make_float4(1.0f, 1.0f, 1.0f, 0.0f);
                color += specColor * pow(RdotEye, dShiny);
            }
        }

        // finally, write output color
        output[outputColorIndex] = rgbaFloatToInt(color);
    }
}

#endif // #ifndef _TEXTURE3D_KERNEL_H_