A simple raytrace program

imprisoned-toroid.cpp

// Imprisoned Toroid -- a simple raytrace program
// written in 2012, shared for educational purposes only.
// A (more than) slightly modified version of this program produces the following:
// https://www.youtube.com/watch?v=oYdYEcjz7_4

#include <math.h>
#include <stdlib.h>

#if defined(WIN32) || defined(_WIN32) || defined(__WIN32__)
  #include <windows.h>     
#endif 

#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>     

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <omp.h>

#define REAL_MAX 1.7e308

#define PHOTON_MAP 1
#define TONE_MAPPING 1
#define MAX_RECURSION 6
#define ROTATE 0

typedef double Real;

struct Vector {
    Real x, y, z;

    Vector() {
        x = y = z = 0;
    }

    Vector(Real x0, Real y0, Real z0) {
        x = x0;
        y = y0;
        z = z0;
    }

    Vector operator*(Real a) const {
        return Vector(x * a, y * a, z * a);
    }

    Vector& operator*=(Real a) {
        x *= a;
        y *= a;
        z *= a;
        return *this;
    }

    Vector operator+(const Vector& v) const {
        return Vector(x + v.x, y + v.y, z + v.z);
    }

    Vector& operator+=(const Vector& v) {
        x += v.x;
        y += v.y;
        z += v.z;
        return *this;
    }

    Vector operator-(const Vector& v) const {
        return Vector(x - v.x, y - v.y, z - v.z);
    }

    Vector operator/(Real d) const {
        return Vector(x / d, y / d, z / d);
    }

    Vector& operator/=(Real d) {
        *this = *this / d;
        return *this;
    }

    Real operator*(const Vector& v) const { // dot product
        return (x * v.x + y * v.y + z * v.z);
    }

    Vector operator%(const Vector& v) const { // cross product
        return Vector(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
    }

    Real length() const {
        return sqrt(x * x + y * y + z * z);
    }

    Vector& normalize() {
        Real l = length();
        x /= l;
        y /= l;
        z /= l;

        return *this;
    }
};

struct Color {
   float r, g, b;

    Color() {
        r = g = b = 0;
    }

    Color(float r0, float g0, float b0) {
        r = r0;
        g = g0;
        b = b0;
    }

    Color operator*(float a) const {
        return Color(r * a, g * a, b * a);
    }

    Color operator*(const Color& c) const {
        return Color(r * c.r, g * c.g, b * c.b);
    }

    Color& operator*=(float a) {
        *this = *this * a;
        return *this;
    }

    Color operator+(const Color& c) const {
        return Color(r + c.r, g + c.g, b + c.b);
    }

    Color& operator+=(const Color& c) {
        *this = *this + c;
        return *this;
    }

    Color& operator+=(float a) {
        r += a;
        g += a;
        b += a;
        return *this;
    }

    Color operator+(float a) const {
        return Color(r + a, g + a, b + a);
    }

    Color operator-(float a) const {
        return Color(r - a, g - a, b - a);
    }

    Color operator/(const Color& c) {
        return Color(r / c.r, g / c.g, b / c.g);
    }
};

const int screenWidth = 600;    // alkalmazás ablak felbontása
const int screenHeight = 600;

struct Ray
{
    Vector point;
    Vector dir;

    void setTargetPoint(const Vector& target) {
        dir = target - point;
        dir.normalize();
    }

#define RND (2.0*(double)rand()/RAND_MAX-1.0)

    void randDir() {
        Real x;
        Real y;
        Real z;

        do {
            x = RND;
            y = RND;
            z = RND;
        } while (x*x + y*y + z*z > 1 || y > 0);

        dir = Vector(x, y, z);
        dir.normalize();
    }

#undef RND
};

#define MAT_REFLECTIVE 2
#define MAT_REFRACTIVE 4
#define MAT_ROUGH 8

#define IS_REFLECTIVE(material) (material & MAT_REFLECTIVE)
#define IS_REFRACTIVE(material) (material & MAT_REFRACTIVE)
#define IS_ROUGH(material) (material & MAT_ROUGH)

struct Object;

struct RayHit
{
    Object* object;
    Vector surfacePoint;
    Vector surfaceNorm;
    Real distance;
};

struct Object
{
    Color ka;
    Color kd;
    Color ks;

    Real shine;

    Color n;
    Color k;

    int material;

    Object(int material = 0)
        :material(material)
    {}

    void setK(Color ka_, Color kd_, Color ks_, Real shine_) {
        ka = ka_;
        kd = kd_;
        ks = ks_;
        shine = shine_;
    }

    virtual RayHit intersect(const Ray& ray) = 0;

    Color fresnel(Real angle) const {
        const Color f0 = ((n-1)*(n-1) + k*k) / ((n+1)*(n+1) + k*k);
        return f0 + (f0 * -1 + 1) * pow(1-cos(angle), 5);
    }

    virtual void addPower(const Vector& point, Real power) {}

    virtual Color getPower(const Vector& point) const {
        return Color(0, 0, 0);
    }
};

struct Sphere :public Object
{
    Vector center;
    Real radius;

    Sphere(const Vector& center, Real radius)
        :center(center),
         radius(radius)
    {}

    Real quadraticPolinomMinPosRoot(Real a, Real b, Real c) const {
        Real determinant = b*b - 4*a*c;

        if (determinant < 0.01) {
            return -1;
        }

        Real root = (-1 * b - sqrt(determinant)) / (2*a);

        if (root >= 0) {
            return root;
        } else {
            return (-1 * b + sqrt(determinant)) / (2*a);
        }
    }

    RayHit intersect(const Ray& ray) {
        Real firstIntersect = quadraticPolinomMinPosRoot(
            ray.dir * ray.dir,
            (ray.point - center) * ray.dir * 2,
            ((ray.point - center) * (ray.point - center)) - (radius * radius)
        );

        RayHit details;
        details.distance = -1;

        if (firstIntersect > 0.01) {
            details.distance = firstIntersect;
            details.object = this;
            details.surfacePoint = ray.point + ray.dir * firstIntersect;
            details.surfaceNorm = normAtPoint(details.surfacePoint);
        }

        return details;
    }

    virtual Vector normAtPoint(const Vector& point) const {
        return (point - center).normalize();
    }
};

struct InsideOutSphere :public Sphere
{
    InsideOutSphere(const Vector& center, Real radius)
        :Sphere(center, radius)
    {}

    Vector normAtPoint(const Vector& point) const {
        return (point - center).normalize() * -1;
    }
};

struct Plane :public Object
{
    Vector point;
    Vector norm;

    struct PhotonMapEntry {
        Vector point;
        Real power;
    };

    #define PHOTON_MAP_SIZE 30000
    PhotonMapEntry photonMap[PHOTON_MAP_SIZE];
    int photonMapSize;

    Plane(const Vector& point, const Vector& norm)
        :Object(MAT_ROUGH),
         point(point),
         norm(norm),
         photonMapSize(0)
    {}

    RayHit intersect(const Ray& ray) {
        Real den = ray.dir * norm;

        RayHit details;

        if (den == 0) {
            details.distance = -1;
        } else {
            Real num = (point - ray.point) * norm;

            details.distance = num / den;
            details.object = this;
            details.surfacePoint = ray.point + ray.dir * details.distance;
            details.surfaceNorm = norm;
        }

        return details;
    }

    void addPower(const Vector& point, Real power) {
        if (photonMapSize == PHOTON_MAP_SIZE) {
            return;
        }

        photonMap[photonMapSize].point = point;
        photonMap[photonMapSize].power = power;

        photonMapSize++;
    }

    Color getPower(const Vector& point) const {
        int photonCount = 0;
        Real radius = 32;
        Real totalPower = 0;

        while (photonCount < photonMapSize / 8) {
            radius *= 2;

            for (int i = 0; i < photonMapSize; i++) {
                if ((point - photonMap[i].point).length() < radius) {
                    totalPower += (photonMap[i].power / (radius * radius * M_PI));
                    photonCount++;
                }
            }
        }

        Color scaledPower = kd * totalPower;

        return scaledPower;
    }

    #undef PHOTON_MAP_SIZE
};

struct Triangle :public Object
{
    Vector va;
    Vector vb;
    Vector vc;

    Vector norm;

    Triangle() {}

    Triangle(const Vector& va, const Vector& vb, const Vector& vc)
        :va(va),
         vb(vb),
         vc(vc)
    {
        norm = ((vb - va) % (vc - va)).normalize();
    }

    RayHit intersect(const Ray& ray) {
        RayHit details;

        Real num = (va - ray.point) * norm;
        Real den = ray.dir * norm;

        if (den == 0 || num/den < 0) {
            details.distance = -1;
        } else {
            Real distance = num/den;
            Vector hitpoint = ray.point + ray.dir * distance;

            Vector areaA = (vb - va) % (hitpoint - va);
            Vector areaB = (vc - vb) % (hitpoint - vb);
            Vector areaC = (va - vc) % (hitpoint - vc);

            if (
                areaA * norm > 0
            &&  areaB * norm > 0
            &&  areaC * norm > 0
            ) {
                Vector weight(areaB.length(), areaC.length(), areaA.length());

                details.distance = distance;
                details.object = this;
                details.surfaceNorm = getNorm(weight);
                details.surfacePoint = hitpoint;
            } else {
                details.distance = -1;
            }
        }

        return details;
    }

    virtual Vector getNorm(const Vector& weight) const {
        return norm;
    }
};

struct CurvedTriangle :public Triangle
{
    Vector normA;
    Vector normB;
    Vector normC;

    CurvedTriangle() {}

    CurvedTriangle(
        const Vector& va, const Vector& vb, const Vector& vc,
        const Vector& normA, const Vector& normB, const Vector& normC
    )
        :Triangle(va, vb, vc),
         normA(normA),
         normB(normB),
         normC(normC)
    {}

    virtual Vector getNorm(const Vector& weight) const {
        return (normA * weight.x + normB * weight.y + normC * weight.z).normalize();
    }
};

struct Toroid :public Object
{
    #define TRIANGLE_COUNT 100
    CurvedTriangle triangles[TRIANGLE_COUNT];
    int triangleIndex;

    Sphere bound;

    Toroid(const Vector& center, Real R, Real r)
        :bound(center, R+r)
    {
        triangleIndex = 0;
        Vector vertices[12];
        Vector norms[12];
        int vertexIndex = 0;

        for (Real u = 0; u < 2 * M_PI + 0.001; u += M_PI / 4) {
            for (Real v = 0; v < 2 * M_PI - 0.001; v += M_PI / 3) {
                Vector vertex(
                    (R + r*cos(v)) * cos(u),
                    (R + r*cos(v)) * sin(u),
                    r * sin(v)
                );

                vertex += center;

                Vector norm =
                    Vector(-1*sin(u) * (r * cos(v) + R), cos(u) * (r * cos(v) + R), 0)
                %   Vector(-1*r*cos(u)*sin(v), -1*r*sin(u)*sin(v), r * cos(v))
                    ;
                norm.normalize();

                vertices[vertexIndex] = vertex;
                norms[vertexIndex] = norm;

                vertexIndex++;
            }

            if (u > 0) {
                for (int i = 0; i < 5; i++) {
                    triangles[triangleIndex++] = CurvedTriangle(
                        vertices[i+0],
                        vertices[i+6],
                        vertices[i+1],
                        norms[i+0],
                        norms[i+6],
                        norms[i+1]
                    );

                    triangles[triangleIndex++] = CurvedTriangle(
                        vertices[i+1],
                        vertices[i+6],
                        vertices[i+7],
                        norms[i+1],
                        norms[i+6],
                        norms[i+7]
                    );
                }

                triangles[triangleIndex++] = CurvedTriangle(
                    vertices[6],
                    vertices[0],
                    vertices[11],
                    norms[6],
                    norms[0],
                    norms[11]
                );

                triangles[triangleIndex++] = CurvedTriangle(
                    vertices[11],
                    vertices[0],
                    vertices[5],
                    norms[11],
                    norms[0],
                    norms[5]
                );
            }

            if (vertexIndex == 12) {
                for (int i = 0; i < 6; i++) {
                    vertices[i] = vertices[i+6];
                    norms[i] = norms[i+6];
                }
                vertexIndex = 6;
            }
        }
    }

    void setMaterial(int material) {
        for (int i = 0; i < triangleIndex; i++) {
            triangles[i].material = material;
        }
    }

    void setNK(const Color& n_, const Color& k_) {
        for (int i = 0; i < triangleIndex; i++) {
            triangles[i].n = n_;
            triangles[i].k = k_;
        }
    }

    void setK(Color ka_, Color kd_, Color ks_, Real shine_) {
        for (int i = 0; i < triangleIndex; i++) {
            triangles[i].setK(ka_, kd_, ks_, shine_);
        }
    }

    RayHit intersect(const Ray& ray) {
        RayHit firstIntersect;

        RayHit boundIntersect = bound.intersect(ray);
        if (boundIntersect.distance < 0) {
            return boundIntersect;
        }

        firstIntersect.distance = REAL_MAX;

        for (int i = 0; i < TRIANGLE_COUNT; i++) {
            RayHit intersect = triangles[i].intersect(ray);

            if (intersect.distance > 0 && intersect.distance < firstIntersect.distance) {
                firstIntersect = intersect;
            }
        }

        if (firstIntersect.distance == REAL_MAX) {
            firstIntersect.distance = -1;
        }

        return firstIntersect;
    }

    #undef TRIANGLE_COUNT
};

struct BubbleBox :public Object
{
    Triangle triangles[12];

    InsideOutSphere bubble;

    Sphere bound;

    BubbleBox(
        const Vector& frontBottomLeft, const Vector& frontBottomRight,
        const Vector& rearBottomRight, const Vector& rearTopRight
    )
        :bubble(
             frontBottomLeft + ((rearTopRight - frontBottomLeft) / 2),
             (frontBottomLeft - frontBottomRight).length() * 0.4
         ),
         bound(
             frontBottomLeft + ((rearTopRight - frontBottomLeft) / 2),
             (frontBottomLeft - frontBottomRight).length() * (sqrt(3) + 0.01)
         )
    {
        const Vector rearTopLeft = rearTopRight + (frontBottomLeft - frontBottomRight);
        const Vector rearBottomLeft = rearTopLeft + (rearBottomRight - rearTopRight);
        const Vector frontTopLeft = frontBottomLeft + (rearTopRight - rearBottomRight);
        const Vector frontTopRight = frontBottomRight + (rearTopRight - rearBottomRight);

        triangles[0] = Triangle(frontBottomLeft, rearBottomLeft, frontTopLeft);
        triangles[1] = Triangle(rearBottomLeft, rearTopLeft, frontTopLeft);
        triangles[2] = Triangle(rearBottomLeft, rearTopRight, rearTopLeft);
        triangles[3] = Triangle(rearBottomLeft, rearBottomRight, rearTopRight);
        triangles[4] = Triangle(rearBottomRight, frontBottomRight, rearTopRight);
        triangles[5] = Triangle(frontBottomRight, frontTopRight, rearTopRight);
        triangles[6] = Triangle(frontBottomLeft, frontTopLeft, frontTopRight);
        triangles[7] = Triangle(frontBottomLeft, frontTopRight, frontBottomRight);
        triangles[8] = Triangle(frontTopLeft, rearTopLeft, rearTopRight);
        triangles[9] = Triangle(frontTopLeft, rearTopRight, frontTopRight);
        triangles[10] = Triangle(frontBottomLeft, rearBottomRight, rearBottomLeft);
        triangles[11] = Triangle(frontBottomLeft, frontBottomRight, rearBottomRight);
    }

    RayHit intersect(const Ray& ray) {
        RayHit firstIntersect;

        RayHit boundIntersect = bound.intersect(ray);
        if (boundIntersect.distance < 0) {
            return boundIntersect;
        }

        firstIntersect.distance = REAL_MAX;

        for (int i = 0; i < 12; i++) {
            RayHit intersect = triangles[i].intersect(ray);

            if (intersect.distance > 0 && intersect.distance < firstIntersect.distance) {
                firstIntersect = intersect;
            }
        }

        RayHit bubbleIntersect = bubble.intersect(ray);

        if (bubbleIntersect.distance > 0 && bubbleIntersect.distance < firstIntersect.distance) {
            firstIntersect = bubbleIntersect;
        }

        if (firstIntersect.distance == REAL_MAX) {
            firstIntersect.distance = -1;
        }

        return firstIntersect;
    }

    void setMaterial(int material) {
        for (int i = 0; i < 12; i++) {
            triangles[i].material = material;
        }

        bubble.material = material;
    }

    void setNK(const Color& n_, const Color& k_) {
        for (int i = 0; i < 12; i++) {
            triangles[i].n = n_;
            triangles[i].k = k_;
        }

        bubble.n = n_;
        bubble.k = k_;
    }

    void setK(Color ka_, Color kd_, Color ks_, Real shine_) {
        for (int i = 0; i < 12; i++) {
            triangles[i].setK(ka_, kd_, ks_, shine_);
        }

        bubble.setK(ka_, kd_, ks_, shine_);
    }
};

struct Lightsource
{
    Vector position;
    Real intensity;

    Lightsource(const Vector& position, Real intensity)
        :position(position),
         intensity(intensity)
    {}
};

struct Scene
{
public:
    Color image[screenWidth*screenHeight];

private:
    #define MAX_OBJ_COUNT 50
    #define MAX_LIGHT_COUNT 10

    Object* objects[MAX_OBJ_COUNT];
    Lightsource* lights[MAX_LIGHT_COUNT];

    int objectCount;
    int lightCount;

    Color ambientLight;

    Vector verticalRotate(Real angle, const Vector& point, const Vector& origo) {
        Vector moved = point - origo;
        Vector rotated(
            moved.x * cos(angle) - moved.z * sin(angle),
            moved.y,
            moved.x * sin(angle) + moved.z * cos(angle)
        );

        return rotated + origo;
    }

    Vector horizontalRotate(Real angle, const Vector& point, const Vector& origo) {
        Vector moved = point - origo;
        Vector rotated(
            moved.x,
            moved.y * cos(angle) - moved.z * sin(angle),
            moved.y * sin(angle) + moved.z * cos(angle)
        );

        return rotated + origo;
    }

    Ray getRayTo(int x, int y) {
        Vector eye(300, 300, 0);
        Vector lookat(300, 300, 600);
        Vector up(0, screenHeight/2, 0);
        Vector right(screenWidth/2, 0, 0);

        Vector target = lookat + right * ((Real)2*x/screenWidth - 1) + up * ((Real)2*y/screenHeight - 1);

        Vector origo(500, 210, 900);

        Vector offset(0, 0, -250);

        eye += offset;
        target += offset;

        eye = horizontalRotate(xRotation, eye, origo);
        target = horizontalRotate(xRotation, target, origo);

        eye = verticalRotate(yRotation, eye, origo);
        target = verticalRotate(yRotation, target, origo);

        Ray ray;
        ray.point = eye;
        ray.setTargetPoint(target);

        return ray;
    }

    Ray reflectedRay(const Ray& in, const RayHit& hit) {
        Ray reflectRay;
        reflectRay.dir = (in.dir - (hit.surfaceNorm * (hit.surfaceNorm * in.dir) * 2)).normalize();
        reflectRay.point = hit.surfacePoint + reflectRay.dir * 0.01;

        return reflectRay;
    }

    bool refractedRay(const Ray& in, const RayHit& hit, Ray& refracted) {
        Real n = hit.object->n.r;
        Real cosalpha = hit.surfaceNorm * in.dir * -1;

        Vector effectiveNorm = hit.surfaceNorm;

        if (cosalpha < 0) {
            effectiveNorm *= -1;
            cosalpha *= -1;
            n = 1/n;
        }

        double disc = 1.0f - ((1.0f - cosalpha * cosalpha) / n / n);

        if (disc > 0.01) {
            refracted.dir =
                ((in.dir / n) +
                effectiveNorm * (
                    cosalpha / n
                -   sqrt(disc)
                )).normalize()
                ;

            refracted.point = hit.surfacePoint + refracted.dir * 0.1;

            return true;
        }

        return false;
    }

    void shootPhoton(Real power, const Ray& ray, int depth = 0) {
        if (depth > MAX_RECURSION) {
            return;
        }

        RayHit hit = firstIntersect(ray);

        if (hit.distance > 0) {
            Color fresnel = hit.object->fresnel(ray.dir * -1 * hit.surfaceNorm);

            if (IS_REFRACTIVE(hit.object->material)) {
                Ray refracted;
                if (refractedRay(ray, hit, refracted)) {
                    Real scaledPower = (1 + fresnel.r * -1) * power;
                    shootPhoton(scaledPower, refracted, depth + 1);
                } else {
                    fresnel = Color(1, 1, 1);
                }
            }

            if (IS_REFLECTIVE(hit.object->material)) {
                Ray reflected = reflectedRay(ray, hit);
                shootPhoton(fresnel.r * power, reflected, depth + 1);
            }

            if (depth > 0 && IS_ROUGH(hit.object->material)) {
                hit.object->addPower(hit.surfacePoint, power);
            }
        }
    }

    Color trace(const Ray& ray, int depth = 0) {
        if (depth > MAX_RECURSION) {
            return ambientLight;
        }

        Color color;

        RayHit hit = firstIntersect(ray);

        if (hit.distance > 0) {
            Color fresnel = hit.object->fresnel(ray.dir * -1 * hit.surfaceNorm);

            if (IS_REFRACTIVE(hit.object->material)) {
                Ray refracted;

                if (refractedRay(ray, hit, refracted)) {
                    Color refractedColor = trace(refracted, depth + 1) * (fresnel * -1 + 1);
                    color += refractedColor;
                } else {
                    fresnel = Color(1, 1, 1);
                }
            }

            if (IS_REFLECTIVE(hit.object->material)) {
                Ray reflected = reflectedRay(ray, hit);
                Color reflectedColor = trace(reflected, depth + 1) * fresnel;
                color += reflectedColor;
            }

            if (IS_ROUGH(hit.object->material)) {
                color += directLight(ray, hit);
#if PHOTON_MAP
                color += hit.object->getPower(hit.surfacePoint);
#endif
            }
        } else {
            color = ambientLight;
        }

        return color;
    }

    RayHit firstIntersect(const Ray& ray) {
        RayHit firstHit;
        firstHit.distance = REAL_MAX;

        for (int i = 0; i < objectCount; i++) {
            RayHit hit = objects[i]->intersect(ray);

            if (hit.distance > 0 && hit.distance < firstHit.distance) {
                firstHit = hit;
            }
        }

        if (firstHit.distance < REAL_MAX) {
            return firstHit;
        } else {
            firstHit.distance = -1;
            return firstHit;
        }
    }

    Color directLight(const Ray& ray, const RayHit& hit) {
        Color color = ambientLight * hit.object->ka;

        for (int i = 0; i < lightCount; i++) {
            Ray shadowray;
            shadowray.point = hit.surfacePoint + hit.surfaceNorm * 0.01;
            shadowray.setTargetPoint(lights[i]->position);

            RayHit barrierHit = firstIntersect(shadowray);

            Vector surfaceToLight = (lights[i]->position - hit.surfacePoint);
            Real lightDistance = surfaceToLight.length();

            if (barrierHit.distance < 0 || barrierHit.distance > lightDistance) {
                Vector surfaceToLightDir = surfaceToLight;
                surfaceToLightDir.normalize();

                Real cosTheta = surfaceToLightDir * hit.surfaceNorm;

                if (cosTheta < 0) {
                    cosTheta = 0;
                }

                Vector hitpointToEyeDir = ray.dir * -1;
                Vector halfDir = ((hitpointToEyeDir + surfaceToLightDir) / 2).normalize();
                Real cosDelta = halfDir * hit.surfaceNorm;

                if (cosDelta < 0) {
                    cosDelta = 0;
                }

                color +=
                    hit.object->fresnel(hitpointToEyeDir * hit.surfaceNorm) *
                    lights[i]->intensity *
                    (hit.object->kd * cosTheta + hit.object->ks * pow(cosDelta, hit.object->shine))
                    ;
            }
        }

        return color;
    }

    Real intensity(const Color& c) {
        return 0.21 * c.r + 0.72 * c.g + 0.07 * c.b;
    }

    void toneMap() {
        Real maxIntensity = intensity(image[0]);
        Real minIntensity = intensity(image[0]);

        for (int i = 0; i < screenWidth * screenHeight; i++) {
            Real intensityI = intensity(image[i]);

            if (intensityI < minIntensity) {
                minIntensity = intensityI;
            } else if (intensityI > maxIntensity) {
                maxIntensity = intensityI;
            }
        }

        Real alpha = 0.8;
        Real avgI = (maxIntensity + minIntensity) / 2;

        for (int i = 0; i < screenWidth * screenHeight; i++) {
            Real mult = (alpha) / (avgI + alpha * intensity(image[i]));
            image[i] *= mult;
        }
    }

public:
    Real yRotation;
    Real xRotation;

    Scene()
        :objectCount(0),
         lightCount(0),
         ambientLight(500, 500, 500),
         yRotation((Real)-1/4 * M_PI),
         xRotation(0.6f)
    {}

    void addObject(Object* object) {
        if (objectCount + 1 == MAX_OBJ_COUNT) {
            return;
        }

        objects[objectCount] = object;
        objectCount++;
    }

    void addLight(Lightsource* light) {
        if (lightCount + 1 == MAX_LIGHT_COUNT) {
            return;
        }

        lights[lightCount] = light;
        lightCount++;
    }

    void render() {
#if PHOTON_MAP
        for (int i = 0; i < lightCount; i++) {
            Ray ray;
            ray.point = lights[i]->position;

            int photonCount = 10000;

            for (int j = 0; j < photonCount; j++) {
                ray.randDir();
                shootPhoton(700000, ray);
            }
        }

        puts("photon shoot done");
#endif

#pragma omp for
        for (int y = 0; y < screenHeight; y++) {
            for (int x = 0; x < screenWidth; x++) {
                Ray ray = getRayTo(x, y);

                image[y*screenWidth + x] = trace(ray);
            }

            printf(" %d row done\n", y);
        }

#if TONE_MAPPING
        toneMap();
#endif
    }

    #undef MAX_OBJ_COUNT
    #undef MAX_LIGHT_COUNT
};


Scene scene;

Lightsource lightA(
    Vector(500, 550, 600),
    30000
);

Plane planeBottom(
    Vector(0,0,0),
    Vector(0, 1, 0)
);

Toroid toroid(
    Vector(300, 210, 1000),
    70,
    20
);

BubbleBox box(
    Vector(100, 20, 800),
    Vector(500, 20, 800),
    Vector(500, 20, 1200),
    Vector(500, 420, 1200)
);

void onInitialization( ) { 
    glViewport(0, 0, screenWidth, screenHeight);

    planeBottom.material = MAT_ROUGH;

    planeBottom.n = Color(1.5, 1.5, 1.5);
    planeBottom.k = Color(3.0, 3.0, 3.0);

    planeBottom.setK(
        Color(0.0, 0.0, 0.0),
        Color(0.682353, 0.682353, 0.682353),
        Color(0, 0, 0),
        12.8
    );

    box.setMaterial(MAT_REFRACTIVE | MAT_REFLECTIVE);

    box.setNK(
        Color(1.5, 1.5, 1.5),
        Color(0, 0, 0)
    );

    box.setK(
        Color(0.0, 0.0, 0.0),
        Color(0.588235, 0.588235, 0.588235),
        Color(0.9, 0.9, 0.9),
        96
    );

    toroid.setMaterial(MAT_REFLECTIVE);

    toroid.setNK(
        Color(0.17, 0.35, 1.5),
        Color(3.1, 2.7, 1.9)
    );

    toroid.setK(
        Color(0.24725, 0.1995, 0.0745),
        Color(0.75164, 0.60648, 0.22648),
        Color(0.628281, 0.555802, 0.366065),
        51.2
    );

    scene.addObject(&planeBottom);
    scene.addObject(&toroid);
    scene.addObject(&box);

    scene.addLight(&lightA);

    scene.render();
}

void onDisplay() {
    glClear(GL_COLOR_BUFFER_BIT);
    glDrawPixels(screenWidth, screenHeight, GL_RGB, GL_FLOAT, scene.image);

    glutSwapBuffers();
}

int main(int argc, char **argv) {
    glutInit(&argc, argv);             
    glutInitWindowSize(screenWidth, screenHeight);         
    glutInitWindowPosition(1200, 100);         
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);  

    glutCreateWindow("Imprisoned toroid");       

    glMatrixMode(GL_MODELVIEW);            
    glLoadIdentity();
    glMatrixMode(GL_PROJECTION);         
    glLoadIdentity();

    onInitialization();                 

    FILE *f = fopen("ray.ppm", "w");
    fprintf(f, "P3\n%d %d\n%d\n ", screenWidth, screenHeight, 255);
    for (int i = 0; i < screenWidth; i++) {
        for (int j = 0; j < screenHeight; j++) {
            fprintf(
                f,
                "%d %d %d ",
                (int)(scene.image[i*screenWidth + j].r * 255),
                (int)(scene.image[i*screenWidth + j].g * 255),
                (int)(scene.image[i*screenWidth + j].b * 255)
            );
        }
    }
    fclose(f);

    glutMainLoop();
    
    return 0;
}