// [header]
// A very basic raytracer example.
// [/header]
// [compile]
// c++ -o raytracer -O3 -Wall raytracer.cpp
// [/compile]
// [ignore]
// Copyright (C) 2012  www.scratchapixel.com
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
// [/ignore]
#include <bits/stdc++.h>
#include <openacc.h>

float dot(float *a,float *b){
    return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);
}

void mul(float *a,float b,float *out){
    out[0] =  (a[0]*b);
    out[1] =  (a[1]*b);
    out[2] =  (a[2]*b);
}

void mul_vecs(float *a,float *b,float *out){
    out[0] =  (a[0]*b[0]);
    out[1] =  (a[1]*b[0]);
    out[2] =  (a[2]*b[0]);
}


void sub(float *a,float *b,float *out){
    out[0] = (a[0] - b[0]);
    out[1] = (a[1] - b[1]);
    out[2] = (a[2] - b[2]);
}

void add(float *a,float *b,float *out){
    out[0] = (a[0] + b[0]);
    out[1] = (a[1] + b[1]);
    out[2] = (a[2] + b[2]);
}

void crossProduct(float *a,float *b,float *out){
    out[0] = a[1] * b[2] - a[2] * b[1];
    out[1] = a[2] * b[0] - a[0] * b[2];
    out[2] = a[0] * b[1] - a[1] * b[0];
}

void normalize(float *a){
        float len2 = dot(a,a);
        if (len2 > 0) {
            float invNor = 1 / sqrt(len2);
            a[0] *= invNor; a[1] *= invNor; a[2] *= invNor;
        }
    }

//[comment]
// This variable controls the maximum recursion depth
//[/comment]
#define MAX_RAY_DEPTH 5

float mix(const float &a, const float &b, const float &mix)
{
    return b * mix + a * (1 - mix);
}

bool intersect(
    float *rayorig,
    float *raydir,
    float triangles[1][3][3],
    int i,
    float *tr)
{
    float kEpsilon = 1e-8;
    float v0[3];
    float v1[3]; 
    float v2[3];
    float orig[3] = {0,0,0};
    for( auto i=0; i <3 ; i++){
        v0[i]=triangles[0][0][i];
        v1[i]=triangles[0][1][i];
        v2[i]=triangles[0][2][i]; 
    }
    float v0v1[3],v0v2[3],pvec[3],tvec[3],qvec[3];
    sub(v1,v0,v0v1);
    sub(v2,v0,v0v2); 
     
    crossProduct(raydir,v0v2,pvec); 
    float det = dot(v0v1,pvec); 

    // ray and triangle are parallel if det is close to 0
    if (fabs(det) < kEpsilon) return false; 
    float invDet = 1 / det; 
 
    sub(orig , v0 , tvec); 
    float u = dot(pvec,tvec) * invDet; 
    if (u < 0 || u > 1) return false; 
 
    crossProduct(tvec,v0v1,qvec); 
    float v = dot(raydir,qvec) * invDet; 
    if (v < 0 || u + v > 1) return false; 
    
    *tr = dot(qvec,v0v2) * invDet;
    // printf("%f",t);
    return true; 
    }

//[comment]
// This is the main trace function. It takes a ray as argument (defined by its origin
// and direction). We test if this rarintersection point, otherwise it returns
// the background color.
//[/comment]
void trace(
    float *rayorig,
    float *raydir,
    const int &depth,
    float triangles[1][3][3],
    float *surface_color,
    float *emission_color,
    float *reflectivity,
    float *transparency, 
    int n ,
    float *pixel
    )
{
    //if (raydir.length() != 1) std::cerr << "Error " << raydir << std::endl;
    float tnear = INFINITY;
    int triag = -1;
    // find intersection of this ray with the sphere in the scene
    for (unsigned i = 0; i < n; ++i) {
        float t0 = INFINITY, t1 = INFINITY;
        if (intersect(rayorig, raydir, triangles ,i , &t0)) {
            // printf("%f" , vab );
            if (t0 < 0) t0 = t1;
            if (t0 < tnear) {
                tnear = t0;
                triag = i;
            }
        }
    }
    // if there's no intersection return black or background color
    if (triag == -1){
        pixel[0] = pixel[1] = pixel[2] = 2; 
        return;
    } 
    float surfaceColor[3] = {0,0,0}; // color of the ray/surfaceof the object intersected by the ray
    float phit[3];
    mul(raydir, tnear , raydir);
    add(rayorig , raydir ,phit); // point of intersection
    float v0[3];
    float v1[3]; 
    float v2[3];
    // float orig[3] = {0,0,0};
    for( auto i=0; i <3 ; i++){
        v0[i]=triangles[0][0][i];
        v1[i]=triangles[0][1][i];
        v2[i]=triangles[0][2][i]; 
    }
    float v0v1[3],v0v2[3];
    sub(v1,v0,v0v1);
    sub(v2,v0,v0v2);  
    float nhit[3];
    crossProduct(v0v1,v0v2,nhit); // normal at the intersection point
    normalize(nhit); // normalize normal direction
    // If the normal and the view direction are not opposite to each other
    // reverse the normal direction. That also means we are inside the sphere so set
    // the inside bool to true. Finally reverse the sign of IdotN which we want
    // positive.
    float bias = 1e-4; // add some bias to the point from which we will be tracing
    bool inside = false;
    if (dot(raydir,nhit) > 0) mul(nhit,-1,nhit), inside = true;
    if ((transparency > 0 || reflectivity > 0) && depth < MAX_RAY_DEPTH) {
        float facingratio = -1*dot(raydir,nhit);
        // change the mix value to tweak the effect
        float fresneleffect = mix(pow(1 - facingratio, 3), 1, 0.1);
        // compute reflection direction (not need to normalize because all vectors
        // are already normalized)
        float refldir[3] , tmp[3] , tmp2[3] , tmp3[3];
        mul(nhit , 2 * dot(raydir,nhit),tmp);
        sub(raydir ,tmp , refldir); 
        normalize(refldir);
        float reflection[3];
        mul(nhit , bias,tmp);
        add(phit , tmp , tmp);
        trace(tmp, refldir, depth + 1 , triangles, surface_color,emission_color,reflectivity,transparency,n,reflection);
        float refraction[3] = {};
        // if the sphere is also transparent compute refraction ray (transmission)
        if (transparency) {
            float ior = 1.1, eta = (inside) ? ior : 1 / ior; // are we inside or outside the surface?
            float cosi = -1* dot(nhit,raydir);
            float k = 1 - eta * eta * (1 - cosi * cosi);
            float refrdir[3];
            mul(raydir,eta,tmp);
            mul( nhit , (eta *  cosi - sqrt(k)) , tmp2);
            add ( tmp , tmp2 , refrdir);
            normalize(refrdir);
            mul(nhit , bias,tmp);
            sub(phit , tmp , tmp);
            trace(tmp, refrdir, depth + 1 , triangles, surface_color,emission_color,reflectivity,transparency,n,refraction);
        }
        // the result is a mix of reflection and refraction (if the sphere is transparent)
        

        mul( reflection , fresneleffect , tmp);
        mul( refraction , (1 - fresneleffect) * (*transparency) , tmp2);
        add( tmp , tmp2 , tmp3);
        mul_vecs ( tmp3 , surface_color, surfaceColor);
    }
    // else {
    //     // it's a diffuse object, no need to raytrace any further
    //     for (unsigned i = 0; i < n; ++i) {
    //         if ( emission_color[0] > 0) {
    //             // this is a light
    //             float transmission[3] = {1,1,1};
    //             float lightDirection[3];
    //             lightDirection[:] = nhit[:];
    //             normalize(lightDirection);
    //             for (unsigned j = 0; j < n; ++j) {
    //                 if (i != j) {
    //                     float t0;

    //                     if (intersect(phit + nhit * bias, lightDirection, triangles, i , &t0)) {
    //                         transmission[0] = 0;
    //                         transmission[1] = 0;
    //                         transmission[2] = 0;
    //                         break;
    //                     }
    //                 }
    //             }

    //             surfaceColor += surface_color * transmission * std::max(float(0), dot(lightDirection,nhit)) * emission_color;
    //         }
    //     }
    // }
    
    add(surfaceColor , emission_color,pixel);
}

//[comment]
// Main rendering function. We compute a camera ray for each pixel of the image
// trace it and return a color. If the ray hits a sphere, we return the color of the
// sphere at the intersection point, else we return the background color.
//[/comment]
void render( float triangles[1][3][3],float *surface_color,float *emission_color,float reflectivity,float transparency, int n)
{
    unsigned width = 640, height = 480;
    float image[width * height][3] = {};
    float invWidth = 1 / float(width), invHeight = 1 / float(height); 
    float fov = 30, aspectratio = width / float(height); 
    float angle = tan(M_PI * 0.5 * fov / 180.);
    float orig[3]= {0,0,0};
    // #pragma omp parallel for collapse(2)
    for (unsigned y = 0; y < height; ++y) {
        for (unsigned x = 0; x < width; ++x) {
            float xx = (2 * ((x + 0.5) * invWidth) - 1) * angle * aspectratio;
            float yy = (1 - 2 * ((y + 0.5) * invHeight)) * angle;
            float raydir[3] ={xx, yy, -1};
            normalize(raydir);
            trace(orig, raydir, 0 , triangles, surface_color,emission_color,&reflectivity,&transparency,n,image[ y*width + x]);
        }
    }
    // Save result to a PPM image (keep these flags if you compile under Windows)
    std::ofstream ofs("./untitled.ppm", std::ios::out | std::ios::binary);
    ofs << "P6\n" << width << " " << height << "\n255\n";
    for (unsigned i = 0; i < width * height; ++i) {
        ofs << (unsigned char)(std::min(float(1), image[i][0]) * 255) <<
               (unsigned char)(std::min(float(1), image[i][1]) * 255) <<
               (unsigned char)(std::min(float(1), image[i][2]) * 255);
    }
    ofs.close();

}

//[comment]
// In the main function, we will create the scene which is composed of 5 spheres
// and 1 light (which is also a sphere). Then, once the scene description is complete
// we render that scene, by calling the render() function.
//[/comment]
int main(int argc, char **argv)
{
    int n = 1;
    srand48(13);
    float triangles[1][3][3] = {{{-1, -1, -5} , { 1, -1, -5} , { 0,  1, -5}}};
    float surface_color[3] = {5.0,  -1, -15};
    float emission_color[3] = {1,1,1};
    float reflectivity = 0.5, transparency=0.5;


    render(triangles,surface_color,emission_color,reflectivity,transparency,n);
    
    return 0;
}