#include <windows.h>

#include "string.h"
#include "glmath.h"

// ----------------------------------------------------------------------------------------------------------------------------

class CCamera
{
public:
    vec3 X, Y, Z, Reference, Position;
    mat4x4 Vin, Pin, VPin, RayMatrix;

public:
    CCamera();
    ~CCamera();

    void CalculateRayMatrix();
    void Look(const vec3 &Position, const vec3 &Reference, bool RotateAroundReference = false);
    bool OnKeyDown(UINT nChar);
    void OnMouseMove(int dx, int dy);
    void OnMouseWheel(short zDelta);
};

// ----------------------------------------------------------------------------------------------------------------------------

class CTriangle
{
public:
    float D, D1, D2, D3, lab, lbc, lca;
    vec3 a, b, c, ab, bc, ca, Color, N, N1, N2, N3;

public:
    CTriangle();
    CTriangle(const vec3 &a, const vec3 &b, const vec3 &c, const vec3 &Color);

    bool Inside(float x, float y, float z);
    bool Inside(const vec3 &Point);
    bool Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance, float &Distance, vec3 &Point);
    bool Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance, float &Distance);
    bool Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance);
};

// ----------------------------------------------------------------------------------------------------------------------------

class RTData
{
public:
    float Distance, TestDistance;
    vec3 Color, Point, TestPoint;
    CTriangle *Triangle;

public:
    RTData();
};

// ----------------------------------------------------------------------------------------------------------------------------

class CVoxel
{
public:
    CTriangle **Triangles;
    int TrianglesCount;

protected:
    int MaxTrianglesCount;
    float Size;
    vec3 Min, Max, MinE, MaxE;

public:
    CVoxel();
    ~CVoxel();

    void Add(CTriangle *Triangle);
    void Delete();
    bool Inside(const vec3 &Point);
    bool Intersect(CTriangle *Triangle);
    bool IntersectEdgesX(CTriangle *Triangle, float x, float y1, float y2, float z1, float z2);
    bool IntersectEdgesY(CTriangle *Triangle, float y, float x1, float x2, float z1, float z2);
    bool IntersectEdgesZ(CTriangle *Triangle, float z, float x1, float x2, float y1, float y2);
    bool IntersectFacesX(CTriangle *Triangle, float D1, float D2);
    bool IntersectFacesY(CTriangle *Triangle, float D1, float D2);
    bool IntersectFacesZ(CTriangle *Triangle, float D1, float D2);
    void Set(const vec3 &Min, float Size);
};

// ----------------------------------------------------------------------------------------------------------------------------

class CUniformGrid
{
protected:
    vec3 Min, Max;

protected:
    int X, Y, Z, Xm1, Ym1, Zm1, XY, XYZ;
    float VoxelSize;
    CVoxel *Voxels;

public:
    CUniformGrid();
    ~CUniformGrid();

public:
    void Delete();
    void Generate(CTriangle *Triangles, int TrianglesCount, float VoxelSize = 1.0f);
    vec3 Traverse(const vec3 &Voxel, const vec3 &Origin, const vec3 &Ray);

    float VoxelToWorldX(float x);
    float VoxelToWorldY(float y);
    float VoxelToWorldZ(float z);
    vec3 VoxelToWorld(const vec3 &Voxel);
    int WorldToVoxelX(float x);
    int WorldToVoxelY(float y);
    int WorldToVoxelZ(float z);
    vec3 WorldToVoxel(const vec3 &World);
};

// ----------------------------------------------------------------------------------------------------------------------------

class CRayTracer
{
private:
    BYTE *ColorBuffer;
    BITMAPINFO ColorBufferInfo;
    int Width, LineWidth, Height;

protected:
    CTriangle *Triangles;
    int TrianglesCount;

private:
    CUniformGrid UniformGrid;

public:
    bool SuperSampling;

public:
    CRayTracer();
    ~CRayTracer();

    bool Init();
    void RayTrace(int x, int y);
    void Resize(int Width, int Height);
    void Destroy();

    void ClearColorBuffer();
    void SwapBuffers(HDC hDC);

protected:
    virtual bool InitScene() = 0;
};

// ----------------------------------------------------------------------------------------------------------------------------

class CMyRayTracer : public CRayTracer
{
protected:
    bool InitScene();
};

// ----------------------------------------------------------------------------------------------------------------------------

class CWnd
{
protected:
    char *WindowName;
    HWND hWnd;
    int Width, Height, x, y, LastX, LastY;

public:
    CWnd();
    ~CWnd();

    bool Create(HINSTANCE hInstance, char *WindowName, int Width, int Height);
    void RePaint();
    void Show(bool Maximized = false);
    void MsgLoop();
    void Destroy();

    void OnKeyDown(UINT Key);
    void OnMouseMove(int X, int Y);
    void OnMouseWheel(short zDelta);
    void OnPaint();
    void OnRButtonDown(int X, int Y);
    void OnSize(int Width, int Height);
};

// ----------------------------------------------------------------------------------------------------------------------------

LRESULT CALLBACK WndProc(HWND hWnd, UINT uiMsg, WPARAM wParam, LPARAM lParam);
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR sCmdLine, int iShow);

#include "cpu_ray_tracer.h"

// ----------------------------------------------------------------------------------------------------------------------------

CString ErrorLog;

// ----------------------------------------------------------------------------------------------------------------------------

CCamera::CCamera()
{
    X = vec3(1.0, 0.0, 0.0);
    Y = vec3(0.0, 1.0, 0.0);
    Z = vec3(0.0, 0.0, 1.0);

    Reference = vec3(0.0, 0.0, 0.0);
    Position = vec3(0.0, 0.0, 5.0);
}

CCamera::~CCamera()
{
}

void CCamera::CalculateRayMatrix()
{
    Vin[0] = X.x; Vin[4] = Y.x; Vin[8] = Z.x;
    Vin[1] = X.y; Vin[5] = Y.y; Vin[9] = Z.y;
    Vin[2] = X.z; Vin[6] = Y.z; Vin[10] = Z.z;

    RayMatrix = Vin * Pin * BiasMatrixInverse * VPin;
}

void CCamera::Look(const vec3 &Position, const vec3 &Reference, bool RotateAroundReference)
{
    this->Reference = Reference;
    this->Position = Position;

    Z = normalize(Position - Reference);
    X = normalize(cross(vec3(0.0f, 1.0f, 0.0f), Z));
    Y = cross(Z, X);

    if(!RotateAroundReference)
    {
        this->Reference = this->Position;
        this->Position += Z * 0.05f;
    }

    CalculateRayMatrix();
}

bool CCamera::OnKeyDown(UINT nChar)
{
    float Distance = 0.125f;

    if(GetKeyState(VK_CONTROL) & 0x80) Distance *= 0.5f;
    if(GetKeyState(VK_SHIFT) & 0x80) Distance *= 2.0f;

    vec3 Up(0.0f, 1.0f, 0.0f);
    vec3 Right = X;
    vec3 Forward = cross(Up, Right);

    Up *= Distance;
    Right *= Distance;
    Forward *= Distance;

    vec3 Movement;

    if(nChar == 'W') Movement += Forward;
    if(nChar == 'S') Movement -= Forward;
    if(nChar == 'A') Movement -= Right;
    if(nChar == 'D') Movement += Right;
    if(nChar == 'R') Movement += Up;
    if(nChar == 'F') Movement -= Up;

    Reference += Movement;
    Position += Movement;

    return Movement.x != 0.0f || Movement.y != 0.0f || Movement.z != 0.0f;
}

void CCamera::OnMouseMove(int dx, int dy)
{
    float sensitivity = 0.25f;

    float hangle = (float)dx * sensitivity;
    float vangle = (float)dy * sensitivity;

    Position -= Reference;

    Y = rotate(Y, vangle, X);
    Z = rotate(Z, vangle, X);

    if(Y.y < 0.0f)
    {
        Z = vec3(0.0f, Z.y > 0.0f ? 1.0f : -1.0f, 0.0f);
        Y = cross(Z, X);
    }

    X = rotate(X, hangle, vec3(0.0f, 1.0f, 0.0f));
    Y = rotate(Y, hangle, vec3(0.0f, 1.0f, 0.0f));
    Z = rotate(Z, hangle, vec3(0.0f, 1.0f, 0.0f));

    Position = Reference + Z * length(Position);

    CalculateRayMatrix();
}

void CCamera::OnMouseWheel(short zDelta)
{
    Position -= Reference;

    if(zDelta < 0 && length(Position) < 500.0f)
    {
        Position += Position * 0.1f;
    }

    if(zDelta > 0 && length(Position) > 0.05f)
    {
        Position -= Position * 0.1f;
    }

    Position += Reference;
}

// ----------------------------------------------------------------------------------------------------------------------------

CCamera Camera;

// ----------------------------------------------------------------------------------------------------------------------------

CTriangle::CTriangle()
{
}

CTriangle::CTriangle(const vec3 &a, const vec3 &b, const vec3 &c, const vec3 &Color) : a(a), b(b), c(c), Color(Color)
{
    ab = b - a; bc = c - b; ca = a - c;

    N = normalize(cross(ab, -ca));
    D = -dot(N, a);

    N1 = normalize(cross(N, ab));
    D1 = -dot(N1, a);

    N2 = normalize(cross(N, bc));
    D2 = -dot(N2, b);

    N3 = normalize(cross(N, ca));
    D3 = -dot(N3, c);

    lab = length(ab); ab /= lab;
    lbc = length(bc); bc /= lbc;
    lca = length(ca); ca /= lca;
}

bool CTriangle::Inside(float x, float y, float z)
{
    if(N1.x * x + N1.y * y + N1.z * z + D1 < 0.0f) return false;
    if(N2.x * x + N2.y * y + N2.z * z + D2 < 0.0f) return false;
    if(N3.x * x + N3.y * y + N3.z * z + D3 < 0.0f) return false;

    return true;
}

bool CTriangle::Inside(const vec3 &Point)
{
    if(dot(N1, Point) + D1 < 0.0f) return false;
    if(dot(N2, Point) + D2 < 0.0f) return false;
    if(dot(N3, Point) + D3 < 0.0f) return false;

    return true;
}

bool CTriangle::Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance, float &Distance, vec3 &Point)
{
    float NdotR = -dot(N, Ray);

    if(NdotR > 0.0f)
    {
        Distance = (dot(N, Origin) + D) / NdotR;

        if(Distance >= 0.0f && Distance < MaxDistance)
        {
            Point = Ray * Distance + Origin;

            return Inside(Point);
        }
    }

    return false;
}

bool CTriangle::Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance, float &Distance)
{
    float NdotR = -dot(N, Ray);

    if(NdotR > 0.0f)
    {
        Distance = (dot(N, Origin) + D) / NdotR;

        if(Distance >= 0.0f && Distance < MaxDistance)
        {
            return Inside(Ray * Distance + Origin);
        }
    }

    return false;
}

bool CTriangle::Intersect(const vec3 &Origin, const vec3 &Ray, float MaxDistance)
{
    float NdotR = -dot(N, Ray);

    if(NdotR > 0.0f)
    {
        float Distance = (dot(N, Origin) + D) / NdotR;

        if(Distance >= 0.0f && Distance < MaxDistance)
        {
            return Inside(Ray * Distance + Origin);
        }
    }

    return false;
}

// ----------------------------------------------------------------------------------------------------------------------------

RTData::RTData()
{
    Distance = 1048576.0f;
    Triangle = NULL;
}

// ----------------------------------------------------------------------------------------------------------------------------

CVoxel::CVoxel()
{
    Triangles = NULL;
    TrianglesCount = 0;
    MaxTrianglesCount = 1;
    Size = 0.0f;
}

CVoxel::~CVoxel()
{
}

void CVoxel::Add(CTriangle *Triangle)
{
    if(TrianglesCount % MaxTrianglesCount == 0)
    {
        CTriangle **OldTriangles = Triangles;

        MaxTrianglesCount *= 2;

        Triangles = new CTriangle*[MaxTrianglesCount];

        for(int i = 0; i < TrianglesCount; i++)
        {
            Triangles[i] = OldTriangles[i];
        }

        if(OldTriangles != NULL)
        {
            delete [] OldTriangles;
        }
    }

    Triangles[TrianglesCount] = Triangle;

    TrianglesCount++;
}

void CVoxel::Delete()
{
    if(Triangles != NULL)
    {
        delete [] Triangles;
        Triangles = NULL;
        TrianglesCount = 0;
        MaxTrianglesCount = 1;
        Size = 0.0f;
        Min = Max = MinE = MaxE = vec3(0.0f);
    }
}

bool CVoxel::Inside(const vec3 &Point)
{
    if(MinE.x < Point.x && Point.x < MaxE.x)
    {
        if(MinE.y < Point.y && Point.y < MaxE.y)
        {
            if(MinE.z < Point.z && Point.z < MaxE.z)
            {
                return true;
            }
        }
    }

    return false;
}

bool CVoxel::IntersectEdgesX(CTriangle *Triangle, float x, float y1, float y2, float z1, float z2)
{
    float NdotR = -Triangle->N.x;

    if(NdotR != 0.0f)
    {
        vec3 Origin = vec3(x, y1, z1);

        float Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x + Distance, Origin.y, Origin.z);
        }

        Origin.z = z2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x + Distance, Origin.y, Origin.z);
        }

        Origin.y = y2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x + Distance, Origin.y, Origin.z);
        }

        Origin.z = z1;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x + Distance, Origin.y, Origin.z);
        }
    }

    return false;
}

bool CVoxel::IntersectEdgesY(CTriangle *Triangle, float y, float x1, float x2, float z1, float z2)
{
    float NdotR = -Triangle->N.y;

    if(NdotR != 0.0f)
    {
        vec3 Origin = vec3(x1, y, z1);

        float Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y + Distance, Origin.z);
        }

        Origin.x = x2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y + Distance, Origin.z);
        }

        Origin.z = z2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y + Distance, Origin.z);
        }

        Origin.x = x1;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y + Distance, Origin.z);
        }
    }

    return false;
}

bool CVoxel::IntersectEdgesZ(CTriangle *Triangle, float z, float x1, float x2, float y1, float y2)
{
    float NdotR = -Triangle->N.z;

    if(NdotR != 0.0f)
    {
        vec3 Origin = vec3(x1, y1, z);

        float Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y, Origin.z + Distance);
        }

        Origin.x = x2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y, Origin.z + Distance);
        }

        Origin.y = y2;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y, Origin.z + Distance);
        }

        Origin.x = x1;

        Distance = (dot(Triangle->N, Origin) + Triangle->D) / NdotR;

        if(Distance >= 0.0f && Distance <= Size)
        {
            return Triangle->Inside(Origin.x, Origin.y, Origin.z + Distance);
        }
    }

    return false;
}

bool CVoxel::IntersectFacesX(CTriangle *Triangle, float D1, float D2)
{
    vec3 *Origin = (vec3*)&Triangle->a;
    vec3 *Ray = (vec3*)&Triangle->ab;
    float *Length = &Triangle->lab;

    float NdotR, d, y, z;

    for(int i = 0; i < 3; i++)
    {
        NdotR = -Ray->x;

        if(NdotR != 0.0f)
        {
            d = (Origin->x - D1) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                y = Ray->y * d + Origin->y;

                if(Min.y <= y && y <= Max.y)
                {
                    z = Ray->z * d + Origin->z;

                    if(Min.z <= z && z <= Max.z)
                    {
                        return true;
                    }
                }
            }

            d = (Origin->x - D2) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                y = Ray->y * d + Origin->y;

                if(Min.y <= y && y <= Max.y)
                {
                    z = Ray->z * d + Origin->z;

                    if(Min.z <= z && z <= Max.z)
                    {
                        return true;
                    }
                }
            }
        }

        Ray++;
        Origin++;
        Length++;
    }

    return false;
}

bool CVoxel::IntersectFacesY(CTriangle *Triangle, float D1, float D2)
{
    vec3 *Origin = (vec3*)&Triangle->a;
    vec3 *Ray = (vec3*)&Triangle->ab;
    float *Length = &Triangle->lab;

    float NdotR, d, x, z;

    for(int i = 0; i < 3; i++)
    {
        NdotR = -Ray->y;

        if(NdotR != 0.0f)
        {
            d = (Origin->y - D1) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                x = Ray->x * d + Origin->x;

                if(Min.x <= x && x <= Max.x)
                {
                    z = Ray->z * d + Origin->z;

                    if(Min.z <= z && z <= Max.z)
                    {
                        return true;
                    }
                }
            }

            d = (Origin->y - D2) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                x = Ray->x * d + Origin->x;

                if(Min.x <= x && x <= Max.x)
                {
                    z = Ray->z * d + Origin->z;

                    if(Min.z <= z && z <= Max.z)
                    {
                        return true;
                    }
                }
            }
        }

        Ray++;
        Origin++;
        Length++;
    }

    return false;
}

bool CVoxel::IntersectFacesZ(CTriangle *Triangle, float D1, float D2)
{
    vec3 *Origin = (vec3*)&Triangle->a;
    vec3 *Ray = (vec3*)&Triangle->ab;
    float *Length = &Triangle->lab;

    float NdotR, d, x, y;

    for(int i = 0; i < 3; i++)
    {
        NdotR = -Ray->z;

        if(NdotR != 0.0f)
        {
            d = (Origin->z - D1) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                x = Ray->x * d + Origin->x;

                if(Min.x <= x && x <= Max.x)
                {
                    y = Ray->y * d + Origin->y;

                    if(Min.y <= y && y <= Max.y)
                    {
                        return true;
                    }
                }
            }

            d = (Origin->z - D2) / NdotR;

            if(0.0f <= d && d <= *Length)
            {
                x = Ray->x * d + Origin->x;

                if(Min.x <= x && x <= Max.x)
                {
                    y = Ray->y * d + Origin->y;

                    if(Min.y <= y && y <= Max.y)
                    {
                        return true;
                    }
                }
            }
        }

        Ray++;
        Origin++;
        Length++;
    }

    return false;
}

bool CVoxel::Intersect(CTriangle *Triangle)
{
    if(Inside(Triangle->a)) return true;
    if(Inside(Triangle->b)) return true;
    if(Inside(Triangle->c)) return true;

    if(IntersectFacesX(Triangle, Min.x, Max.x)) return true;
    if(IntersectFacesY(Triangle, Min.y, Max.y)) return true;
    if(IntersectFacesZ(Triangle, Min.z, Max.z)) return true;

    if(IntersectEdgesX(Triangle, Min.x, Min.y, Max.y, Min.z, Max.z)) return true;
    if(IntersectEdgesY(Triangle, Min.y, Min.x, Max.x, Min.z, Max.z)) return true;
    if(IntersectEdgesZ(Triangle, Min.z, Min.x, Max.x, Min.y, Max.y)) return true;

    return false;
}

void CVoxel::Set(const vec3 &Min, float Size)
{
    this->Size = Size;
    this->Min = Min;
    this->Max = this->Min + Size;
    this->MinE = this->Min - 0.001f;
    this->MaxE = this->Max + 0.001f;
}

// ----------------------------------------------------------------------------------------------------------------------------

CUniformGrid::CUniformGrid()
{
    X = Y = Z = Xm1 = Ym1 = Zm1 = XY = XYZ = 0;
    VoxelSize = 0.0f;
    Voxels = NULL;
}

CUniformGrid::~CUniformGrid()
{
}

void CUniformGrid::Delete()
{
    if(Voxels != NULL)
    {
        for(int i = 0; i < XYZ; i++)
        {
            Voxels[i].Delete();
        }

        Min = Max = vec3(0.0f);
        X = Y = Z = Xm1 = Ym1 = Zm1 = XY = XYZ = 0;
        VoxelSize = 0.0f;
        delete [] Voxels;
        Voxels = NULL;
    }
}

void CUniformGrid::Generate(CTriangle *Triangles, int TrianglesCount, float VoxelSize)
{
    Delete();

    if(Triangles != NULL && TrianglesCount > 0)
    {
        this->VoxelSize = VoxelSize;

        CTriangle *LastTriangle = Triangles + TrianglesCount;

        Min = Max = Triangles->a;

        for(CTriangle *Triangle = Triangles; Triangle < LastTriangle; Triangle++)
        {
            if(Triangle->a.x < Min.x) Min.x = Triangle->a.x;
            if(Triangle->a.y < Min.y) Min.y = Triangle->a.y;
            if(Triangle->a.z < Min.z) Min.z = Triangle->a.z;

            if(Triangle->b.x < Min.x) Min.x = Triangle->b.x;
            if(Triangle->b.y < Min.y) Min.y = Triangle->b.y;
            if(Triangle->b.z < Min.z) Min.z = Triangle->b.z;

            if(Triangle->c.x < Min.x) Min.x = Triangle->c.x;
            if(Triangle->c.y < Min.y) Min.y = Triangle->c.y;
            if(Triangle->c.z < Min.z) Min.z = Triangle->c.z;

            if(Triangle->a.x > Max.x) Max.x = Triangle->a.x;
            if(Triangle->a.y > Max.y) Max.y = Triangle->a.y;
            if(Triangle->a.z > Max.z) Max.z = Triangle->a.z;

            if(Triangle->b.x > Max.x) Max.x = Triangle->b.x;
            if(Triangle->b.y > Max.y) Max.y = Triangle->b.y;
            if(Triangle->b.z > Max.z) Max.z = Triangle->b.z;

            if(Triangle->c.x > Max.x) Max.x = Triangle->c.x;
            if(Triangle->c.y > Max.y) Max.y = Triangle->c.y;
            if(Triangle->c.z > Max.z) Max.z = Triangle->c.z;
        }

        Min /= VoxelSize; Max /= VoxelSize;

        Min.x = floor(Min.x); Max.x = ceil(Max.x);
        Min.y = floor(Min.y); Max.y = ceil(Max.y);
        Min.z = floor(Min.z); Max.z = ceil(Max.z);

        if(Min.x == Max.x) Max.x += 1.0f;
        if(Min.y == Max.y) Max.y += 1.0f;
        if(Min.z == Max.z) Max.z += 1.0f;

        X = (int)(Max.x - Min.x); Xm1 = X - 1;
        Y = (int)(Max.y - Min.y); Ym1 = Y - 1;
        Z = (int)(Max.z - Min.z); Zm1 = Z - 1;

        XY = X * Y;
        XYZ = XY * Z;

        Min *= VoxelSize; Max *= VoxelSize;

        Voxels = new CVoxel[XYZ];

        for(int z = 0; z < Z; z++)
        {
            for(int y = 0; y < Y; y++)
            {
                for(int x = 0; x < X; x++)
                {
                    Voxels[XY * z + X * y + x].Set(VoxelToWorld(vec3((float)x, (float)y, (float)z)), VoxelSize);
                }
            }
        }

        for(CTriangle *Triangle = Triangles; Triangle < LastTriangle; Triangle++)
        {
            vec3 tmin = Triangle->a, tmax = tmin;

            if(Triangle->b.x < tmin.x) tmin.x = Triangle->b.x;
            if(Triangle->b.y < tmin.y) tmin.y = Triangle->b.y;
            if(Triangle->b.z < tmin.z) tmin.z = Triangle->b.z;

            if(Triangle->b.x > tmax.x) tmax.x = Triangle->b.x;
            if(Triangle->b.y > tmax.y) tmax.y = Triangle->b.y;
            if(Triangle->b.z > tmax.z) tmax.z = Triangle->b.z;

            if(Triangle->c.x < tmin.x) tmin.x = Triangle->c.x;
            if(Triangle->c.y < tmin.y) tmin.y = Triangle->c.y;
            if(Triangle->c.z < tmin.z) tmin.z = Triangle->c.z;

            if(Triangle->c.x > tmax.x) tmax.x = Triangle->c.x;
            if(Triangle->c.y > tmax.y) tmax.y = Triangle->c.y;
            if(Triangle->c.z > tmax.z) tmax.z = Triangle->c.z;

            int vminx = WorldToVoxelX(tmin.x), vmaxx = WorldToVoxelX(tmax.x);
            int vminy = WorldToVoxelY(tmin.y), vmaxy = WorldToVoxelY(tmax.y);
            int vminz = WorldToVoxelZ(tmin.z), vmaxz = WorldToVoxelZ(tmax.z);

            if(vminx >= X) vminx = Xm1; if(vmaxx >= X) vmaxx = Xm1;
            if(vminy >= Y) vminy = Ym1; if(vmaxy >= Y) vmaxy = Ym1;
            if(vminz >= Z) vminz = Zm1; if(vmaxz >= Z) vmaxz = Zm1;

            for(int z = vminz; z <= vmaxz; z++)
            {
                for(int y = vminy; y <= vmaxy; y++)
                {
                    for(int x = vminx; x <= vmaxx; x++)
                    {
                        CVoxel *Voxel = Voxels + (XY * z + X * y + x);

                        if(Voxel->Intersect(Triangle))
                        {
                            Voxel->Add(Triangle);
                        }
                    }
                }
            }
        }
    }
}

vec3 CUniformGrid::Traverse(const vec3 &Voxel, const vec3 &Origin, const vec3 &Ray)
{
    vec3 voxel = Voxel, step, out, t, delta = VoxelSize / Ray;

    if(Ray.x < 0.0f)
    {
        step.x = -1.0f;
        out.x = voxel.x <= 0.0f ? voxel.x - 1.0f : -1.0f;
        t.x = (VoxelToWorldX(voxel.x) - Origin.x) / Ray.x;
    }
    else
    {
        step.x = 1.0f;
        out.x = voxel.x >= X ? voxel.x + 1 : X;
        t.x = (VoxelToWorldX(voxel.x + 1.0f) - Origin.x) / Ray.x;
    }

    if(Ray.y < 0.0f)
    {
        step.y = -1.0f;
        out.y = voxel.y <= 0.0f ? voxel.y - 1.0f : -1.0f;
        t.y = (VoxelToWorldY(voxel.y) - Origin.y) / Ray.y;
    }
    else
    {
        step.y = 1.0f;
        out.y = voxel.y >= Y ? voxel.y + 1 : Y;
        t.y = (VoxelToWorldY(voxel.y + 1.0f) - Origin.y) / Ray.y;
    }

    if(Ray.z < 0.0f)
    {
        step.z = -1.0f;
        out.z = voxel.z <= 0.0f ? voxel.z - 1.0f : -1.0f;
        t.z = (VoxelToWorldZ(voxel.z) - Origin.z) / Ray.z;
    }
    else
    {
        step.z = 1.0f;
        out.z = voxel.z >= Z ? voxel.z + 1 : Z;
        t.z = (VoxelToWorldZ(voxel.z + 1.0f) - Origin.z) / Ray.z;
    }

    delta *= step;

    while(1)
    {
        int x = (int)voxel.x, y = (int)voxel.y, z = (int)voxel.z;

        if(x >= 0 && x < X && y >= 0 && y < Y && z >= 0 && z < Z)
        {
            CVoxel *Voxel = Voxels + (XY * z + X * y + x);
            CTriangle **Triangles = Voxel->Triangles, *Triangle;
            int TrianglesCount = Voxel->TrianglesCount;

            RTData rtdata;

            for(int i = 0; i < TrianglesCount; i++)
            {
                Triangle = Triangles[i];

                if(Triangle->Intersect(Origin, Ray, rtdata.Distance, rtdata.TestDistance, rtdata.TestPoint))
                {
                    if(Voxel->Inside(rtdata.TestPoint))
                    {
                        rtdata.Point = rtdata.TestPoint;
                        rtdata.Distance = rtdata.TestDistance;
                        rtdata.Triangle = Triangle;
                    }
                }
            }

            if(rtdata.Triangle)
            {
                rtdata.Color = rtdata.Triangle->Color;

                float NdotL = -dot(rtdata.Triangle->N, Ray);

                if(NdotL < 0.0f)
                {
                    NdotL = 0.0f;
                }

                rtdata.Color *= 0.75f * NdotL + 0.25f;

                return rtdata.Color;
            }
        }

        float min = t.x;

        if(t.y < min) min = t.y;
        if(t.z < min) min = t.z;

        if(t.x == min)
        {
           voxel.x += step.x;
           if(voxel.x == out.x) break;
           t.x += delta.x;
        }

        if(t.y == min)
        {
           voxel.y += step.y;
           if(voxel.y == out.y) break;
           t.y += delta.y;
        }

        if(t.z == min)
        {
           voxel.z += step.z;
           if(voxel.z == out.z) break;
           t.z += delta.z;
        }
    }

    return vec3(0.0f);
}

float CUniformGrid::VoxelToWorldX(float x)
{
    return x * VoxelSize + Min.x;
}

float CUniformGrid::VoxelToWorldY(float y)
{
    return y * VoxelSize + Min.y;
}

float CUniformGrid::VoxelToWorldZ(float z)
{
    return z * VoxelSize + Min.z;
}

vec3 CUniformGrid::VoxelToWorld(const vec3 &Voxel)
{
    return Voxel * VoxelSize + Min;
}

int CUniformGrid::WorldToVoxelX(float x)
{
    return (int)floor((x - Min.x) / VoxelSize);
}

int CUniformGrid::WorldToVoxelY(float y)
{
    return (int)floor((y - Min.y) / VoxelSize);
}

int CUniformGrid::WorldToVoxelZ(float z)
{
    return (int)floor((z - Min.z) / VoxelSize);
}

vec3 CUniformGrid::WorldToVoxel(const vec3 &World)
{
    vec3 voxel = (World - Min) / VoxelSize;

    voxel.x = floor(voxel.x);
    voxel.y = floor(voxel.y);
    voxel.z = floor(voxel.z);

    return voxel;
}

// ----------------------------------------------------------------------------------------------------------------------------

CRayTracer::CRayTracer()
{
    ColorBuffer = NULL;

    Triangles = NULL;
    TrianglesCount = 0;

    SuperSampling = false;
}

CRayTracer::~CRayTracer()
{
}

bool CRayTracer::Init()
{
    if(InitScene() == false)
    {
        return false;
    }

    UniformGrid.Generate(Triangles, TrianglesCount);

    return true;
}

void CRayTracer::RayTrace(int x, int y)
{
    if(ColorBuffer != NULL && Triangles != NULL && TrianglesCount > 0)
    {
        vec3 Color, Voxel = UniformGrid.WorldToVoxel(Camera.Position);

        if(!SuperSampling)
        {
            Color = UniformGrid.Traverse(Voxel, Camera.Position, normalize(*(vec3*)&(Camera.RayMatrix * vec4((float)x + 0.5f, (float)y + 0.5f, 0.0f, 1.0f))));
        }
        else
        {
            for(float yy = 0.125f; yy < 1.0f; yy += 0.25f)
            {
                for(float xx = 0.125f; xx < 1.0f; xx += 0.25f)
                {
                    Color += UniformGrid.Traverse(Voxel, Camera.Position, normalize(*(vec3*)&(Camera.RayMatrix * vec4((float)x + xx, (float)y + yy, 0.5f, 1.0f))));
                }
            }

            Color /= 16.0f;
        }

        BYTE *colorbuffer = (LineWidth * y + x) * 3 + ColorBuffer;

        colorbuffer[2] = Color.r <= 0.0f ? 0 : Color.r >= 1.0 ? 255 : (BYTE)(Color.r * 255);
        colorbuffer[1] = Color.g <= 0.0f ? 0 : Color.g >= 1.0 ? 255 : (BYTE)(Color.g * 255);
        colorbuffer[0] = Color.b <= 0.0f ? 0 : Color.b >= 1.0 ? 255 : (BYTE)(Color.b * 255);
    }
}

void CRayTracer::Resize(int Width, int Height)
{
    this->Width = Width;
    this->Height = Height;

    if(ColorBuffer != NULL)
    {
        delete [] ColorBuffer;
        ColorBuffer = NULL;
    }

    if(Width > 0 && Height > 0)
    {
        LineWidth = Width;

        int WidthMod4 = Width % 4;

        if(WidthMod4 > 0)
        {
            LineWidth += 4 - WidthMod4;
        }

        ColorBuffer = new BYTE[LineWidth * Height * 3];

        memset(&ColorBufferInfo, 0, sizeof(BITMAPINFOHEADER));
        ColorBufferInfo.bmiHeader.biSize = sizeof(BITMAPINFOHEADER);
        ColorBufferInfo.bmiHeader.biPlanes = 1;
        ColorBufferInfo.bmiHeader.biBitCount = 24;
        ColorBufferInfo.bmiHeader.biCompression = BI_RGB;
        ColorBufferInfo.bmiHeader.biWidth = LineWidth;
        ColorBufferInfo.bmiHeader.biHeight = Height;

        Camera.VPin[0] = 1.0f / (float)Width;
        Camera.VPin[5] = 1.0f / (float)Height;

        float tany = tan(45.0f / 360.0f * (float)M_PI), aspect = (float)Width / (float)Height;

        Camera.Pin[0] = tany * aspect;
        Camera.Pin[5] = tany;
        Camera.Pin[10] = 0.0f;
        Camera.Pin[14] = -1.0f;

        Camera.CalculateRayMatrix();
    }
}

void CRayTracer::Destroy()
{
    UniformGrid.Delete();

    if(Triangles != NULL)
    {
        delete [] Triangles;
        Triangles = NULL;
        TrianglesCount = 0;
    }

    if(ColorBuffer != NULL)
    {
        delete [] ColorBuffer;
        ColorBuffer = NULL;
    }
}

void CRayTracer::ClearColorBuffer()
{
    if(ColorBuffer != NULL)
    {
        memset(ColorBuffer, 0, LineWidth * Height * 3);
    }
}

void CRayTracer::SwapBuffers(HDC hDC)
{
    if(ColorBuffer != NULL)
    {
        StretchDIBits(hDC, 0, 0, Width, Height, 0, 0, Width, Height, ColorBuffer, &ColorBufferInfo, DIB_RGB_COLORS, SRCCOPY);
    }
}

// ----------------------------------------------------------------------------------------------------------------------------

bool CMyRayTracer::InitScene()
{
    bool Error = false;

    if(Error)
    {
        return false;
    }

    int BoxesCount = 256;

    TrianglesCount = 12 * BoxesCount;

    Triangles = new CTriangle[TrianglesCount];

    int t = 0;

    srand(GetTickCount());

    for(int i = 0; i < BoxesCount; i++)
    {
        vec3 m = -8.0f + 16.0f * vec3((float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX);
        vec3 s = 0.25f + 1.5f * vec3((float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX);
        vec3 color = vec3((float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX, (float)rand() / (float)RAND_MAX);

        mat3x3 RS;
        
        RS = RS * mat3x3(rotate(360.0f * (float)rand() / (float)RAND_MAX, vec3(1.0f, 0.0f, 0.0f)));
        RS = RS * mat3x3(rotate(360.0f * (float)rand() / (float)RAND_MAX, vec3(0.0f, 1.0f, 0.0f)));
        RS = RS * mat3x3(rotate(360.0f * (float)rand() / (float)RAND_MAX, vec3(0.0f, 0.0f, 1.0f)));
        
        RS = RS * mat3x3(scale(s.x, s.y, s.z));

        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f, -0.5f,  0.5f), m + RS * vec3( 0.5f, -0.5f, -0.5f), m + RS * vec3( 0.5f,  0.5f, -0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f,  0.5f, -0.5f), m + RS * vec3( 0.5f,  0.5f,  0.5f), m + RS * vec3( 0.5f, -0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f, -0.5f, -0.5f), m + RS * vec3(-0.5f, -0.5f,  0.5f), m + RS * vec3(-0.5f,  0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f,  0.5f,  0.5f), m + RS * vec3(-0.5f,  0.5f, -0.5f), m + RS * vec3(-0.5f, -0.5f, -0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f,  0.5f,  0.5f), m + RS * vec3( 0.5f,  0.5f,  0.5f), m + RS * vec3( 0.5f,  0.5f, -0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f,  0.5f, -0.5f), m + RS * vec3(-0.5f,  0.5f, -0.5f), m + RS * vec3(-0.5f,  0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f, -0.5f, -0.5f), m + RS * vec3( 0.5f, -0.5f, -0.5f), m + RS * vec3( 0.5f, -0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f, -0.5f,  0.5f), m + RS * vec3(-0.5f, -0.5f,  0.5f), m + RS * vec3(-0.5f, -0.5f, -0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f, -0.5f,  0.5f), m + RS * vec3( 0.5f, -0.5f,  0.5f), m + RS * vec3( 0.5f,  0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f,  0.5f,  0.5f), m + RS * vec3(-0.5f,  0.5f,  0.5f), m + RS * vec3(-0.5f, -0.5f,  0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3( 0.5f, -0.5f, -0.5f), m + RS * vec3(-0.5f, -0.5f, -0.5f), m + RS * vec3(-0.5f,  0.5f, -0.5f), color);
        Triangles[t++] = CTriangle(m + RS * vec3(-0.5f,  0.5f, -0.5f), m + RS * vec3( 0.5f,  0.5f, -0.5f), m + RS * vec3( 0.5f, -0.5f, -0.5f), color);
    }

    Camera.Look(vec3(0.0f, 0.0f, 10.0f), vec3(0.0f, 0.0f, 0.0f), true);

    return true;
}

// ----------------------------------------------------------------------------------------------------------------------------

CMyRayTracer RayTracer;

// ----------------------------------------------------------------------------------------------------------------------------

CWnd::CWnd()
{
}

CWnd::~CWnd()
{
}

bool CWnd::Create(HINSTANCE hInstance, char *WindowName, int Width, int Height)
{
    WNDCLASSEX WndClassEx;

    memset(&WndClassEx, 0, sizeof(WNDCLASSEX));

    WndClassEx.cbSize = sizeof(WNDCLASSEX);
    WndClassEx.style = CS_OWNDC | CS_HREDRAW | CS_VREDRAW;
    WndClassEx.lpfnWndProc = WndProc;
    WndClassEx.hInstance = hInstance;
    WndClassEx.hIcon = LoadIcon(NULL, IDI_APPLICATION);
    WndClassEx.hIconSm = LoadIcon(NULL, IDI_APPLICATION);
    WndClassEx.hCursor = LoadCursor(NULL, IDC_ARROW);
    WndClassEx.lpszClassName = "Win32CPURayTracerWindow";

    if(RegisterClassEx(&WndClassEx) == 0)
    {
        ErrorLog.Set("RegisterClassEx failed!");
        return false;
    }

    this->WindowName = WindowName;

    this->Width = Width;
    this->Height = Height;

    DWORD Style = WS_OVERLAPPEDWINDOW | WS_CLIPSIBLINGS | WS_CLIPCHILDREN;

    if((hWnd = CreateWindowEx(WS_EX_APPWINDOW, WndClassEx.lpszClassName, WindowName, Style, 0, 0, Width, Height, NULL, NULL, hInstance, NULL)) == NULL)
    {
        ErrorLog.Set("CreateWindowEx failed!");
        return false;
    }

    return RayTracer.Init();
}

void CWnd::RePaint()
{
    x = y = 0;
    InvalidateRect(hWnd, NULL, FALSE);
}

void CWnd::Show(bool Maximized)
{
    RECT dRect, wRect, cRect;

    GetWindowRect(GetDesktopWindow(), &dRect);
    GetWindowRect(hWnd, &wRect);
    GetClientRect(hWnd, &cRect);

    wRect.right += Width - cRect.right;
    wRect.bottom += Height - cRect.bottom;

    wRect.right -= wRect.left;
    wRect.bottom -= wRect.top;

    wRect.left = dRect.right / 2 - wRect.right / 2;
    wRect.top = dRect.bottom / 2 - wRect.bottom / 2;

    MoveWindow(hWnd, wRect.left, wRect.top, wRect.right, wRect.bottom, FALSE);

    ShowWindow(hWnd, Maximized ? SW_SHOWMAXIMIZED : SW_SHOWNORMAL);
}

void CWnd::MsgLoop()
{
    MSG Msg;

    while(GetMessage(&Msg, NULL, 0, 0) > 0)
    {
        TranslateMessage(&Msg);
        DispatchMessage(&Msg);
    }
}

void CWnd::Destroy()
{
    RayTracer.Destroy();

    DestroyWindow(hWnd);
}

void CWnd::OnKeyDown(UINT Key)
{
    switch(Key)
    {
        case VK_F1:
            RayTracer.SuperSampling = !RayTracer.SuperSampling;
            RePaint();
            break;
    }

    if(Camera.OnKeyDown(Key))
    {
        RePaint();
    }
}

void CWnd::OnMouseMove(int X, int Y)
{
    if(GetKeyState(VK_RBUTTON) & 0x80)
    {
        Camera.OnMouseMove(LastX - X, LastY - Y);

        LastX = X;
        LastY = Y;

        RePaint();
    }
}

void CWnd::OnMouseWheel(short zDelta)
{
    Camera.OnMouseWheel(zDelta);

    RePaint();
}

void CWnd::OnPaint()
{
    PAINTSTRUCT ps;

    HDC hDC = BeginPaint(hWnd, &ps);

    static DWORD Start;
    static bool RayTracing;

    if(x == 0 && y == 0)
    {
        RayTracer.ClearColorBuffer();

        Start = GetTickCount();

        RayTracing = true;
    }

    DWORD start = GetTickCount();

    while(GetTickCount() - start < 125 && y < Height)
    {
        int x16 = x + 16, y16 = y + 16;

        for(int yy = y; yy < y16; yy++)
        {
            if(yy < Height)
            {
                for(int xx = x; xx < x16; xx++)
                {
                    if(xx < Width)
                    {
                        RayTracer.RayTrace(xx, yy);
                    }
                    else
                    {
                        break;
                    }
                }
            }
            else
            {
                break;
            }
        }

        x = x16;

        if(x >= Width)
        {
            x = 0;
            y = y16;
        }
    }

    RayTracer.SwapBuffers(hDC);

    if(RayTracing)
    {
        if(y >= Height)
        {
            RayTracing = false;
        }

        DWORD End = GetTickCount();

        CString text = WindowName;

        text.Append(" - %dx%d", Width, Height);
        text.Append(", Time: %.03f s", (float)(End - Start) * 0.001f);

        SetWindowText(hWnd, text);

        InvalidateRect(hWnd, NULL, FALSE);
    }

    EndPaint(hWnd, &ps);
}

void CWnd::OnRButtonDown(int X, int Y)
{
    LastX = X;
    LastY = Y;
}

void CWnd::OnSize(int Width, int Height)
{
    this->Width = Width;
    this->Height = Height;

    RayTracer.Resize(Width, Height);

    RePaint();
}

// ----------------------------------------------------------------------------------------------------------------------------

CWnd Wnd;

// ----------------------------------------------------------------------------------------------------------------------------

LRESULT CALLBACK WndProc(HWND hWnd, UINT uiMsg, WPARAM wParam, LPARAM lParam)
{
    switch(uiMsg)
    {
        case WM_CLOSE:
            PostQuitMessage(0);
            break;

        case WM_MOUSEMOVE:
            Wnd.OnMouseMove(LOWORD(lParam), HIWORD(lParam));
            break;

        case 0x020A: // WM_MOUSWHEEL
            Wnd.OnMouseWheel(HIWORD(wParam));
            break;

        case WM_KEYDOWN:
            Wnd.OnKeyDown((UINT)wParam);
            break;

        case WM_PAINT:
            Wnd.OnPaint();
            break;

        case WM_RBUTTONDOWN:
            Wnd.OnRButtonDown(LOWORD(lParam), HIWORD(lParam));
            break;

        case WM_SIZE:
            Wnd.OnSize(LOWORD(lParam), HIWORD(lParam));
            break;

        default:
            return DefWindowProc(hWnd, uiMsg, wParam, lParam);
    }

    return 0;
}

int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR sCmdLine, int iShow)
{
    SetThreadPriority(GetCurrentThread(), THREAD_PRIORITY_HIGHEST);

    if(Wnd.Create(hInstance, "CPU ray tracer 02 - Uniform grid", 800, 600))
    {
        Wnd.Show();
        Wnd.MsgLoop();
    }
    else
    {
        MessageBox(NULL, ErrorLog, "Error", MB_OK | MB_ICONERROR);
    }

    Wnd.Destroy();

    return 0;
}