#define _Dot(a,b) ((a).x*(b).x+(a).y*(b).y+(a).z*(b).z)
#define _Len2(a) _Dot(a,a)

float f_IntersectSphere_t(Vector v_ray_org, Vector v_ray,
                          Vector v_sph_org, float f_radius)
{
    float d, b;
    Vector q;

    q.x = v_ray_org.x - v_sph_org.x;
    q.y = v_ray_org.y - v_sph_org.y;
    q.z = v_ray_org.z - v_sph_org.z;

    b = _Dot(v_ray, q);
    d = b * b - _Len2(q) + f_radius * f_radius;

    if(d < 0)
        return -1;
    // žádný průsečík

    return b - (float)sqrt(d);
}

Vector v_IntersectSphere_t(Vector v_ray_org, Vector v_ray,
                           Vector v_sph_org, float f_radius)
{
    float d, b;
    Vector q;

    q.x = v_ray_org.x - v_sph_org.x;
    q.y = v_ray_org.y - v_sph_org.y;
    q.z = v_ray_org.z - v_sph_org.z;

    b = _Dot(v_ray, q);
    d = b * b - _Len2(q) + f_radius * f_radius;

    if(d < 0)
        return v_ray_org;
    // žádný průsečík

    d = b - (float)sqrt(d);
    // vzdálenost

    v_ray_org.x += v_ray.x * d;
    v_ray_org.y += v_ray.y * d;
    v_ray_org.z += v_ray.z * d;
    // bližší průsečík

    return v_ray_org;
}

int RayCast(Vector *start, Vector *end, Face *fac)
{
    Vector mid, v1, v2, pl;
    float a, b, t, d;

    a = start.x * fac->normal.a + start.y * fac->normal.b +
        start.z * fac->normal.c + fac->normal.d;
    b = end.x * fac->normal.a + end.y * fac->normal.b +
        end.z * fac->normal.c + fac->normal.d;
    // vzdálenosti od roviny

    if(a * b > 0)
        return 0;
    // pokud jsou oba body na stejné straně

    t = -a / (b - a);
    mid.x = start->x + t * (end->x - start->x);
    mid.y = start->y + t * (end->y - start->y);
    mid.z = start->z + t * (end->z - start->z);
    // spočítá průsečík

    v1.x = -fac->normal.a;
    v1.y = -fac->normal.b;
    v1.z = -fac->normal.c;
    // první vektor (pro všechny roviny stejný)

    v2.x = fac->vertex[1]->x - fac->vertex[2]->x;
    v2.y = fac->vertex[1]->y - fac->vertex[2]->y;
    v2.z = fac->vertex[1]->z - fac->vertex[2]->z;
    // první hranový vektor

    _Cross(v1, v2, pl);
    _Normalize(pl);
    d = - (pl.x * fac->vertex[2]->x +
           pl.y * fac->vertex[2]->y +
           pl.z * fac->vertex[2]->z);
    // spočítá rovinu

    if(mid.x * pl.x + mid.y * pl.y + mid.z * pl.z + d < -1e-3)
        return 0;
    // otestuje průsečík

    v2.x = fac->vertex[0]->x - fac->vertex[1]->x;
    v2.y = fac->vertex[0]->y - fac->vertex[1]->y;
    v2.z = fac->vertex[0]->z - fac->vertex[1]->z;
    // druhý vektor

    _Cross(v1, v2, pl);
    _Normalize(pl);
    d = - (pl.x * fac->vertex[1]->x +
           pl.y * fac->vertex[1]->y +
           pl.z * fac->vertex[1]->z);
    // spočítá rovinu

    if(mid.x * pl.x + mid.y * pl.y + mid.z * pl.z + d < -1e-3)
        return 0;
    // otestuje průsečík

    v2.x = fac->vertex[2]->x - fac->vertex[0]->x;
    v2.y = fac->vertex[2]->y - fac->vertex[0]->y;
    v2.z = fac->vertex[2]->z - fac->vertex[0]->z;
    // treti vektor

    _Cross(v1, v2, pl);
    _Normalize(pl);
    d = - (pl.x * fac->vertex[0]->x +
           pl.y * fac->vertex[0]->y +
           pl.z * fac->vertex[0]->z);
    // spočítá poslední rovinu

    if(mid.x * pl.x + mid.y * pl.y + mid.z * pl.z + d < -1e-3)
        return 0;
    // otestuje průsečík

    return 1;
}

int RayCast_v2(Vector *start, Vector *end, Face *fac)
{
    Vector mid;
    Vector t[3];
    int i, p;
    float a, b, k, angle;
    float ilen, plen;

    a = start.x * fac->normal.a + start.y * fac->normal.b +
        start.z * fac->normal.c + fac->normal.d;
    b = end.x * fac->normal.a + end.y * fac->normal.b +
        end.z * fac->normal.c + fac->normal.d;
    // vzdalenosti od roviny

    if(a * b > 0)
        return 0;
    // zase zkontroluje jestli nejsou na stejné straně

    k = -a / (b - a);
    mid.x = start->x + k * (end->x - start->x);
    mid.y = start->y + k * (end->y - start->y);
    mid.z = start->z + k * (end->z - start->z);
    // spočítá průsečík

    angle = 0.0;
    p = 2;

    t[p].x = fac->vertex[p]->x - mid.x;
    t[p].y = fac->vertex[p]->y - mid.y;
    t[p].z = fac->vertex[p]->z - mid.z;
    // potřebujeme ještě předchozí vektor

    plen = _Len(t[p]);

    for(i = 0; i < 3; i ++) {
        t[i].x = fac->vertex[i]->x - mid.x;
        t[i].y = fac->vertex[i]->y - mid.y;
        t[i].z = fac->vertex[i]->z - mid.z;
        //
        ilen = _Len(t[i]);
        angle += (float)acos((DOT(t[i], t[p])) / (ilen * plen));
        plen = ilen;
        p = i;
    }
    // počítá úhly

    if(fabs(angle - D_PI) < 0.001)
        return 1;
    // pokud protne, vrátí 1

    return 0;
}

struct Face {
    int n_id, n_used;
    int n_material;
    //
    Face *n_neighbour[3];
    //
    Vertex *vertex[3];
    //
    Plane normal;
};

int b_SameVertex(Vertex *p_vertex_a, Vertex *p_vertex_b)
{
    float f, d;

    f = (p_vertex_a->x - p_vertex_b->x);
    f *= f;
    d = (p_vertex_a->y - p_vertex_b->y);
    f += d * d;
    d = (p_vertex_a->z - p_vertex_b->z);
    // vzdálenost mezi vrcholy
    
    return ((f + d * d) < 0.01f);
}

void Find_FaceNeighbours(Object *p_object)
{
    int n_edge_a, n_edge_b;
    int n_edge_c, n_edge_d;
    int i, j, k, l;

    for(i = 0; i < p_object->n_face_num; i ++) {
        p_object->face[i].n_neighbour[0] = NULL;
        p_object->face[i].n_neighbour[1] = NULL;
        p_object->face[i].n_neighbour[2] = NULL;
    }

    for(i = 0; i < p_object->n_face_num; i ++) {
        for(n_edge_a = 0; n_edge_a < 3; n_edge_a ++) {
            n_edge_b = (1 + n_edge_a) % 3;
            l = (1 + n_edge_b) % 3;
            // indexy vertexů hrany (a, b)
            for(j = i + 1; j < p_object->n_face_num; j ++) {
                for(n_edge_c = 0; n_edge_c < 3; n_edge_c ++) {
                    n_edge_d = (n_edge_c + 1) % 3;
                    // indexy vertexů druhé hrany (c, d)

                    if((b_SameVertex(p_object->face[i].vertex[n_edge_a],
                                    p_object->face[j].vertex[n_edge_d]) &&
                       b_SameVertex(p_object->face[i].vertex[n_edge_b],
                                    p_object->face[j].vertex[n_edge_c])) ||
                       (b_SameVertex(p_object->face[i].vertex[n_edge_a],
                                    p_object->face[j].vertex[n_edge_c]) &&
                       b_SameVertex(p_object->face[i].vertex[n_edge_b],
                                    p_object->face[j].vertex[n_edge_d]))) {
                        k = (n_edge_d + 1) % 3;
                        if(p_object->face[i].normal.a * p_object->face[j].vertex[k]->x +
                           p_object->face[i].normal.b * p_object->face[j].vertex[k]->y +
                           p_object->face[i].normal.c * p_object->face[j].vertex[k]->z +
                           p_object->face[i].normal.d > -.01 &&
                           p_object->face[j].normal.a * p_object->face[i].vertex[l]->x +
                           p_object->face[j].normal.b * p_object->face[i].vertex[l]->y +
                           p_object->face[j].normal.c * p_object->face[i].vertex[l]->z +
                           p_object->face[j].normal.d > -.01) {
                            p_object->face[i].n_neighbour[n_edge_a] = &p_object->face[j];
                            p_object->face[j].n_neighbour[n_edge_c] = &p_object->face[i];
                        }
                    }
                    // facy se vzájemně označí za sousedy, pokud
                    // leží odvrácené od sebe (objekt je konvexí)
                }
            }
        }
    }
    // hledá sousední facy
}