#define MAX_DISTANCE        (8.0f)
#define MIN_DELTA           (0.001f)
#define MAX_REFLECTIONS     (8)

#define MOD(x, y)           ((x) - (y) * floor((x) / (y));
#ifndef M_PI
#define M_PI                (3.1415926535897932384626433832795f)
#endif

float2 rand2n(float2 *seed)
{
    *seed += (float2)(-1.0f, 1.0f);
    // implementation based on: lumina.sourceforge.net/Tutorials/Noise.html
    float s = 0.0f;
    return (float2)(fract(sin(dot((*seed).xy, (float2)(12.9898f, 78.233f))) * 43758.5453f, &s),
        fract(cos(dot((*seed).xy, (float2)(4.898f, 7.23f))) * 23421.631f, &s));
}

float3 ortho(float3 v)
{
    //  See : http://lolengine.net/blog/2013/09/21/picking-orthogonal-vector-combing-coconuts
    return fabs(v.x) > fabs(v.z) ? (float3)(-v.y, v.x, 0.0f) : (float3)(0.0, -v.z, v.y);
}

float3 getSampleBiased(float3 dir, float power, float2 *seed)
{
    dir = normalize(dir);
    float3 o1 = normalize(ortho(dir));
    float3 o2 = normalize(cross(dir, o1));
    float2 r = rand2n(seed);
    r.x = r.x * 2.0f * M_PI;
    r.y = pow(r.y, 1.0f / (power + 1.0f));
    float oneminus = sqrt(1.0f - r.y*r.y);
    return cos(r.x)*oneminus*o1 + sin(r.x)*oneminus*o2 + r.y*dir;
}

float3 getSample(float3 dir, float2 *seed)
{
    return getSampleBiased(dir, 0.0f, seed); // <- unbiased!
}

float3 getCosineWeightedSample(float3 dir, float2 *seed)
{
    return getSampleBiased(dir, 1.0f, seed);
}

float3 getConeSample(float3 dir, float extent, float2 *seed)
{
    // Formula 34 in GI Compendium
    dir = normalize(dir);
    float3 o1 = normalize(ortho(dir));
    float3 o2 = normalize(cross(dir, o1));
    float2 r = rand2n(seed);
    r.x = r.x * 2.0f * M_PI;
    r.y = 1.0 - r.y*extent;
    float oneminus = sqrt(1.0f - r.y*r.y);
    return cos(r.x)*oneminus*o1 + sin(r.x)*oneminus*o2 + r.y*dir;
}

float3 vRotateX(float3 p, float angle)
{
    float c = cos(angle);
    float s = sin(angle);
    return (float3)(p.x, c*p.y + s*p.z, -s*p.y + c*p.z);
}

float3 vRotateY(float3 p, float angle)
{
    float c = cos(angle);
    float s = sin(angle);
    return (float3)(c*p.x - s*p.z, p.y, s*p.x + c*p.z);
}

float3 vRotateZ(float3 p, float angle)
{
    float c = cos(angle);
    float s = sin(angle);
    return (float3)(c*p.x + s*p.y, -s*p.x + c*p.y, p.z);
}

__kernel void mainimage(
    __global float4 *tile,
    int width,
    int height,
    int tile_size,
    int offset_x,
    int offset_y,
    int sample)
{
    int gid_x = get_global_id(0);
    int gid_y = get_global_id(1);

    float fsample = (float)(sample + 1);

    float2 screen = (float2)(gid_x + offset_x, gid_y + offset_y);
    float2 seed = screen * fsample * 1.2256f;
    float2 offset = rand2n(&seed);
    float2 p = (screen + offset) / (float)height;

    float4 final = (float4)(
        p.x,
        p.y,
        atan2(p.x, p.y),
        1.0f);

    float a = -3.0f * p.x + 3.0f - p.y;
    if (a > 0.0 && a < 1.0)
        final = (float4)(1.0, 1.0, 1.0, 1.0);

    tile[gid_x + gid_y * tile_size] = mix(
        tile[gid_x + gid_y * tile_size],
        final,
        1.0f / fsample);
}