//
// LICENCE:
//  The MIT License (MIT)
//
//  Copyright (c) 2017 Karim Naaji, karim.naaji@gmail.com
//
//  Permission is hereby granted, free of charge, to any person obtaining a copy
//  of this software and associated documentation files (the "Software"), to deal
//  in the Software without restriction, including without limitation the rights
//  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
//  copies of the Software, and to permit persons to whom the Software is
//  furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included in all
//  copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
//  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
//  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
//  OUT OF OR IN CONNECTION WITH THE
//
// HOWTO:
//  #define SDF_IMPLEMENTATION
//  #define SDF_DEBUG // Only if assertions need to be checked
//  #include "sdf.h"
//

#ifndef SDF_H
#define SDF_H

// ------------------------------------------------------------------------------------------------
// SDF PUBLIC API
//

#ifndef SDF_HELPERS
#include <stdlib.h> // malloc, calloc, free
#endif

typedef struct sd_vertex {
    union {
        float v[3];
        struct {
            float x;
            float y;
            float z;
        };
    };
} sd_vertex_t;

typedef struct sd_mesh {
    sd_vertex_t* vertices;
    unsigned int* indices;
    unsigned int nindices;
    unsigned int nvertices;
} sd_mesh_t;

typedef sd_vertex_t sd_vec3_t;

unsigned int* sd_texture3d(sd_mesh_t const* mesh, int resx, int resy, int resz);
sd_mesh_t* sd_mesh_alloc(unsigned int nvertices, unsigned int nindices);
void sd_mesh_free(sd_mesh_t* mesh);

// Voxelizer Helpers, define your own if needed
#ifndef SDF_HELPERS
#define SDF_HELPERS 1
#define SD_MIN(a, b) (a > b ? b : a)
#define SD_MAX(a, b) (a > b ? a : b)
#define SD_FINDMINMAX(x0, x1, x2, min, max) \
    min = max = x0;                         \
    if (x1 < min) min = x1;                 \
    if (x1 > max) max = x1;                 \
    if (x2 < min) min = x2;                 \
    if (x2 > max) max = x2;
#define SD_CLAMP(v, lo, hi) SD_MAX(lo, SD_MIN(hi, v))
#define SD_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
#define SD_FREE(T) free(T)
#define SD_CALLOC(T, N) ((T*) calloc(N * sizeof(T), 1))
#define SD_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
#ifdef SDF_DEBUG
#define SD_ASSERT(STMT) if (!(STMT)) { *(int *)0 = 0; }
#else
#define SD_ASSERT(STMT)
#endif // SDF_DEBUG
#endif // SDF_HELPERS

//
// END SDF PUBLIC API
// ------------------------------------------------------------------------------------------------

#endif // SDF_H

#ifdef SDF_IMPLEMENTATION

#include <math.h>       // ceil, fabs & al.
#include <string.h>     // memcpy

typedef unsigned long   u64;
typedef unsigned int    u32;
typedef unsigned short  u16;
typedef unsigned char   u8;
typedef signed long     s64;
typedef signed short    s16;
typedef signed int      s32;
typedef signed char     s8;
typedef float           f32;
typedef double          f64;

#define SDF_EPSILON               (0.0000001)

typedef struct sd_aabb {
    sd_vertex_t min;
    sd_vertex_t max;
} sd_aabb_t;

typedef struct sd_triangle {
    union {
        sd_vertex_t vertices[3];
        struct {
            sd_vertex_t p1;
            sd_vertex_t p2;
            sd_vertex_t p3;
        };
    };
    sd_vertex_t centroid;
    sd_vec3_t normal;
} sd_triangle_t;

typedef struct sd_ray {
    sd_vertex_t origin;
    sd_vec3_t direction;
} sd_ray_t;

sd_vec3_t sd__vec3_sub(sd_vec3_t* a, sd_vec3_t* b) {
    sd_vec3_t v;
    v.x = a->x - b->x;
    v.y = a->y - b->y;
    v.z = a->z - b->z;
    return v;
}

sd_vec3_t sd__vec3_add(sd_vec3_t* a, sd_vec3_t* b)
{
    sd_vec3_t v;
    v.x = a->x + b->x;
    v.y = a->y + b->y;
    v.z = a->z + b->z;
    return v;
}

sd_vec3_t sd__vec3_multiply(sd_vec3_t* a, float scalar)
{
    sd_vec3_t v;
    v.x = a->x * scalar;
    v.y = a->y * scalar;
    v.z = a->z * scalar;
    return v;
}

sd_vec3_t sd__vec3_multiply(sd_vec3_t* a, sd_vec3_t* b)
{
    sd_vec3_t v;
    v.x = a->x * b->x;
    v.y = a->y * b->y;
    v.z = a->z * b->z;
    return v;
}

sd_vec3_t sd__vec3_divide(sd_vec3_t* a, float scalar)
{
    sd_vec3_t v = *a;
    v.x /= scalar;
    v.y /= scalar;
    v.z /= scalar;
    return v;
}

sd_vec3_t sd__vec3_divide(sd_vec3_t* a, sd_vec3_t* b)
{
    sd_vec3_t v;
    v.x = a->x / b->x;
    v.y = a->y / b->y;
    v.z = a->z / b->z;
    return v;
}

float sd__vec3_dot(sd_vec3_t* v1, sd_vec3_t* v2)
{
    return v1->x * v2->x + v1->y * v2->y + v1->z * v2->z;
}

sd_vec3_t sd__vec3_normalize(sd_vec3_t* v)
{
    float invLength = 1.f / sqrt(sd__vec3_dot(v, v));
    sd_vec3_t normalized = sd__vec3_multiply(v, invLength);
    return normalized;
}

void sd_mesh_free(sd_mesh_t* mesh)
{
    SD_FREE(mesh->vertices);
    mesh->vertices = NULL;
    mesh->nvertices = 0;
    SD_FREE(mesh->indices);
    mesh->indices = NULL;
    mesh->nindices = 0;
    SD_FREE(mesh);
}

sd_mesh_t* sd_mesh_alloc(int nvertices, int nindices)
{
    sd_mesh_t* mesh = SD_MALLOC(sd_mesh_t, 1);
    mesh->indices = SD_CALLOC(unsigned int, nindices);
    mesh->vertices = SD_CALLOC(sd_vertex_t, nvertices);
    mesh->nindices = nindices;
    mesh->nvertices = nvertices;
    return mesh;
}

float sd__map_to_3dgrid(float position, float voxelSize, bool min)
{
    float vox = (position + (position < 0.f ? -1.f : 1.f) * voxelSize * 0.5f) / voxelSize;
    return (min ? floor(vox) : ceil(vox)) * voxelSize;
}

sd_vec3_t sd__vec3_cross(sd_vec3_t* v1, sd_vec3_t* v2)
{
    sd_vec3_t cross;
    cross.x = v1->y * v2->z - v1->z * v2->y;
    cross.y = v1->z * v2->x - v1->x * v2->z;
    cross.z = v1->x * v2->y - v1->y * v2->x;
    return cross;
}

bool sd__vertex_equals_epsilon(sd_vertex_t* v1, sd_vertex_t* v2) {
    return fabs(v1->x - v2->x) < SDF_EPSILON &&
           fabs(v1->y - v2->y) < SDF_EPSILON &&
           fabs(v1->z - v2->z) < SDF_EPSILON;
}

void sd__aabb_init(sd_aabb_t* aabb)
{
    aabb->max.x = aabb->max.y = aabb->max.z = -INFINITY;
    aabb->min.x = aabb->min.y = aabb->min.z = INFINITY;
}

sd_aabb_t sd__triangle_aabb(sd_triangle_t* triangle)
{
    sd_aabb_t aabb;

    sd__aabb_init(&aabb);

    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            aabb.max.v[i] = SD_MAX(aabb.max.v[i], triangle->vertices[j].v[i]);
            aabb.min.v[i] = SD_MIN(aabb.min.v[i], triangle->vertices[j].v[i]);
        }
    }

    return aabb;
}

sd_vertex_t sd__aabb_center(sd_aabb_t* a)
{
    sd_vertex_t boxcenter;
    boxcenter = sd__vec3_add(&a->min, &a->max);
    boxcenter = sd__vec3_multiply(&boxcenter, 0.5f);
    return boxcenter;
}

sd_vertex_t sd__aabb_half_gc(sd_aabb_t* a)
{
    sd_vertex_t gc;

    gc.x = fabs(a->max.x - a->min.x) * 0.5f;
    gc.y = fabs(a->max.y - a->min.y) * 0.5f;
    gc.z = fabs(a->max.z - a->min.z) * 0.5f;

    return gc;
}

sd_aabb_t sd__aabb_merge(sd_aabb_t* a, sd_aabb_t* b)
{
    sd_aabb_t merge;

    merge.min.x = SD_MIN(a->min.x, b->min.x);
    merge.min.y = SD_MIN(a->min.y, b->min.y);
    merge.min.z = SD_MIN(a->min.z, b->min.z);

    merge.max.x = SD_MAX(a->max.x, b->max.x);
    merge.max.y = SD_MAX(a->max.y, b->max.y);
    merge.max.z = SD_MAX(a->max.z, b->max.z);

    return merge;
}

s8 sd__ray_triangle_intersection(sd_ray_t* ray,
    sd_triangle_t* triangle)
{
    sd_vec3_t e1 = sd__vec3_sub(&triangle->p2, &triangle->p1);
    sd_vec3_t e2 = sd__vec3_sub(&triangle->p3, &triangle->p1);
    sd_vec3_t s1 = sd__vec3_cross(&ray->direction, &e2);
    f32 divisor = sd__vec3_dot(&s1, &e1);

    if (divisor == 0.0) {
        return -1;
    }

    f32 invDivisor = 1.0f / divisor;
    sd_vec3_t d = sd__vec3_sub(&ray->origin, &triangle->p1);
    f32 b1 = sd__vec3_dot(&d, &s1) * invDivisor;
    if (b1 < 0.0f || b1 > 1.0f) {
        return -1;
    }

    sd_vec3_t s2 = sd__vec3_cross(&d, &e1);
    f32 b2 = sd__vec3_dot(&d, &s2) * invDivisor;
    if (b2 < 0.0f || b1 + b2 > 1.0f) {
        return -1;
    }

    f32 t = sd__vec3_dot(&e2, &s2) * invDivisor;
    if (t < 0.0f) {
        return -1;
    }

    return 1;
}

sd_aabb_t sd__aabb(sd_triangle_t* triangles,
    u32 nTriangles)
{
    sd_aabb_t aabb;

    sd__aabb_init(&aabb);

    for (u32 i = 0; i < nTriangles; ++i) {
        if (aabb.min.x == INFINITY) {
            aabb = sd__triangle_aabb(&triangles[i]);
        } else {
            sd_aabb_t taabb = sd__triangle_aabb(&triangles[i]);
            aabb = sd__aabb_merge(&aabb, &taabb);
        }
    }

    return aabb;
}

sd_triangle_t* sd__triangles(sd_mesh_t const* m, u32 nTriangles)
{
    sd_triangle_t* triangles = SD_MALLOC(sd_triangle_t, nTriangles);

    for (u32 i = 0, ti = 0; i < m->nindices; i += 3, ti++) {
        SD_ASSERT(m->indices[i+0] < m->nvertices);
        SD_ASSERT(m->indices[i+1] < m->nvertices);
        SD_ASSERT(m->indices[i+2] < m->nvertices);

        for (u32 index = 0; index < 3; index++) {
            u32 j = m->indices[i + index];
            triangles[ti].vertices[index] = m->vertices[j];
        }

        sd_vec3_t centroid = {0.0f, 0.0f, 0.0f};
        for (u32 index = 0; index < 3; index++) {
            centroid = sd__vec3_add(&centroid, &triangles[ti].vertices[index]);
        }
        centroid = sd__vec3_divide(&centroid, 3.0f);

        sd_vertex_t p0 = triangles[ti].vertices[0];
        sd_vertex_t p1 = triangles[ti].vertices[1];
        sd_vertex_t p2 = triangles[ti].vertices[2];

        sd_vec3_t e0 = sd__vec3_sub(&p1, &p0);
        sd_vec3_t e1 = sd__vec3_sub(&p2, &p0);

        sd_vec3_t normal;
        normal = sd__vec3_cross(&e0, &e1);
        normal = sd__vec3_normalize(&normal);

        triangles[ti].normal = normal;
        triangles[ti].centroid = centroid;
    }

    return triangles;
}

unsigned int* sd_texture3d(sd_mesh_t const* m,
    int resx,
    int resy,
    int resz)
{
    u32* texture3d = SD_MALLOC(u32, resx * resy * resz);
    u32 nTriangles = (u32)(m->nindices / 3.0f);
    sd_triangle_t* triangles = sd__triangles(m, nTriangles);
    sd_aabb_t aabb = sd__aabb(triangles, nTriangles);
    sd_vec3_t extent = sd__vec3_sub(&aabb.max, &aabb.min);
    sd_vec3_t size = {(f32)resx, (f32)resy, (f32)resz};
    sd_vec3_t step = sd__vec3_divide(&extent, &size);
    f32 minDistance = +INFINITY;
    f32 maxDistance = -INFINITY;
    f32* distances = SD_MALLOC(f32, resx * resy * resz);
    s8* signs = SD_MALLOC(s8, resx * resy * resz);

    for (f32 z = aabb.min.z, tz = 0; z < aabb.max.z && tz < resz; z += step.z, tz++) {
        for (f32 y = aabb.min.y, ty = 0; y < aabb.max.y && ty < resy; y += step.y, ty++) {
            for (f32 x = aabb.min.x, tx = 0; x < aabb.max.x && tx < resx; x += step.x, tx++) {
                sd_vertex_t p = {x + step.x * 0.5f, y + step.y * 0.5f, z + step.z * 0.5f};
                u32 index = tx + ty * resx + tz * (resx * resy);
                f32 triangleDistance = INFINITY;
                s8 sign;

                for (u32 ti = 0; ti < nTriangles; ++ti) {
                    sd_vec3_t pToTriangle = sd__vec3_sub(&triangles[ti].centroid, &p);
                    f32 distance = sqrtf(sd__vec3_dot(&pToTriangle, &pToTriangle));
                    if (distance < triangleDistance) {
                        triangleDistance = distance;
                        f32 dotProduct = sd__vec3_dot(&pToTriangle, &triangles[ti].normal);
                        sign = dotProduct > 0.0f ? -1 : +1;
                    }
                }

                signs[index] = sign;
                distances[index] = triangleDistance;

                minDistance = SD_MIN(minDistance, triangleDistance);
                maxDistance = SD_MAX(maxDistance, triangleDistance);
            }
        }
    }

    for (u32 tz = 0; tz < resz; tz++) {
        for (u32 ty = 0; ty < resy; ty++) {
            for (u32 tx = 0; tx < resx; tx++) {
                u32 index = tx + ty * resx + tz * (resx * resy);
                f32 distanceNormalized = distances[index] / fabsf(maxDistance - minDistance);
                f32 signedDistance = (signs[index] * distanceNormalized * 0.5f + 0.5f) * 255.0f;

                texture3d[index] = 0xff000000
                     | ((u8)signedDistance) << 16
                     | ((u8)signedDistance) <<  8
                     | ((u8)signedDistance) <<  0;
            }
        }
    }

    SD_FREE(distances);
    SD_FREE(triangles);

    return texture3d;
}

#undef SDF_EPSILON

#endif // SDF_IMPLEMENTATION