#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);
}

// ------------------------------------------------------------------------------------------------
float3 pal(float t, float3 a, float3 b, float3 c, float3 d )
{
    return a + b*cos(6.28318f * (c * t + d));
}

// ------------------------------------------------------------------------------------------------
#define DE2_SCALE   (1.5f)
#define DE2_OFFSET  (2.0f)
#define DE2_ITER    (2*48)
float DE2(float3 z, float *orbit)
{
    float r;
    float d;

    float min_dist = 1e9f;

    z = vRotateZ(z, M_PI / 2.0f);

    for (int n = 0; n < DE2_ITER; n++)
    {
        z = vRotateX(z, 0.31);
        z = vRotateY(z, 0.31);
        if(z.x + z.y < 0.0f) { z.xy = -z.yx; } // fold 1
        if(z.x + z.z < 0.0f) { z.xz = -z.zx; } // fold 2
        if(z.y + z.z < 0.0f) { z.zy = -z.yz; } // fold 3  
        //z.xy = fabs(z.xy);
        z = z * DE2_SCALE - DE2_OFFSET * (DE2_SCALE - 1.0f);

        d = length(z);
        min_dist = min(min_dist, d);
    }

    *orbit = min_dist;
    return d * pow(DE2_SCALE, -(float)DE2_ITER);
}

// ------------------------------------------------------------------------------------------------
__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 = 2.0f * (screen + offset) / (float)height;
    p -= (float2)((float)width / (float)height, 1.0f);

    float4 final = (float4)(0.0f, 0.0f, 0.0f, 1.0f);

    float orbit = 0.0f;
    float val = DE2((float3)(p.x, p.y, 0.0f), &orbit);
    final.xyz = pal(
        2.0f * orbit,
        (float3)(0.8,0.5,0.4),
        (float3)(0.2,0.4,0.2),
        (float3)(2.0,1.0,1.0),
        (float3)(0.0,0.25,0.25));

    //if (val <= 0.0001f)
    //  final.xyz = (float3)(1.0f, 1.0f, 1.0f);

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