delone.c

// [INFO] Average time ~ 16.1 msec (10k triangles)

#include "../../common.h"

#if VIS
#include <GL/glew.h>
#include <GLFW/glfw3.h>

struct GL {
    GLuint VAO;
    GLuint VBO;
    GLFWwindow *window;
    GLuint program;
};
#endif

enum LOCALIZE_ACTION {
    LOCAL,
    SPLIT,
    DISCARD
};

const f32 PI = 3.1415926f;
const f32 ERR = 0.5f;

// NOTE: controlls how close we can get to an edge
// before snapping to it 
const f32 EDGE_TOLERANCE = 1e-6;
const f32 CYCLES_PER_MSEC_CPU = 3200000.0f;
const f32 COEFF = 0.7f;

const u32 N = 10000;
const u32 FLIPSTACK_SIZE = 100;
const u32 SLEEP_USEC = 50000;

const u32 WIDTH  = 12800000;
const u32 HEIGHT = 7200000;

#if VIS
const char *vertex_shader =
"#version 330\n"
"in vec2 position;\n"
"void main() {\n"
"    gl_Position = vec4(position.x / 12800000.0f * 1.6f - 0.8f,"
"                               position.y / 7200000.0f * 1.6f - 0.8f,"
"                               0.0f, 1.0f);\n"
"}";
const char *fragment_shader =
"#version 330\n"
"out vec4 color;\n"
"void main() {\n"
"    color = vec4(1.0f);\n"
"}";

static struct GL
init_opengl()
{
    struct GL res;
    
    glfwInit();
    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
    glfwWindowHint(GLFW_RESIZABLE, GLFW_FALSE);
    glfwWindowHint(GLFW_SAMPLES, 8);
    
    res.window = glfwCreateWindow(1280.0f, 720.0f, "MathMod Lab2", NULL, NULL);
    
    glfwMakeContextCurrent(res.window);
    glewExperimental = GL_TRUE;
    glewInit();
    
    //glPointSize(2.0f);
    //glLineWidth(2.0f);
    glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
    
    glGenBuffers(1, &res.VBO);
    glGenVertexArrays(1, &res.VAO);
    
    glBindVertexArray(res.VAO);
    glBindBuffer(GL_ARRAY_BUFFER, res.VBO);
    glEnableVertexAttribArray(0);
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, NULL);
    
    GLuint vs = glCreateShader(GL_VERTEX_SHADER);
    GLuint fs = glCreateShader(GL_FRAGMENT_SHADER);
    
    glShaderSource(vs, 1, &vertex_shader, NULL);
    glShaderSource(fs, 1, &fragment_shader, NULL);
    
    glCompileShader(vs);
    glCompileShader(fs);
    
    res.program = glCreateProgram();
    
    glAttachShader(res.program, vs);
    glAttachShader(res.program, fs);
    glLinkProgram(res.program);
    
    glBindVertexArray(res.VAO);
    
    return(res);
}
#endif

static u32
edge_index(f32 *triangles, u32 neighbour_triangle_index, u32 triangle_index, u32 edge_number)
{
    f32 edge_x1 = triangles[triangle_index + edge_number * 2 + 0];
    f32 edge_y1 = triangles[triangle_index + edge_number * 2 + 1];
    
    f32 edge_x2 = triangles[triangle_index + (edge_number * 2 + 2) % 6];
    f32 edge_y2 = triangles[triangle_index + (edge_number * 2 + 3) % 6];
    
    f32 x1, y1, x2, y2;
    
    for (u32 edge = 0; edge < 3; ++edge) {
        x1 = triangles[neighbour_triangle_index + (edge * 2 + 0) % 6];
        y1 = triangles[neighbour_triangle_index + (edge * 2 + 1) % 6];
        
        x2 = triangles[neighbour_triangle_index + (edge * 2 + 2) % 6];
        y2 = triangles[neighbour_triangle_index + (edge * 2 + 3) % 6];
        
        if (fabs(x1 - edge_x1) < ERR && fabs(y1 - edge_y1) < ERR &&
            fabs(x2 - edge_x2) < ERR && fabs(y2 - edge_y2) < ERR) {
            return(edge);
        }
        
        // NOTE: the same edge with reversed vertices
        if (fabs(x1 - edge_x2) < ERR && fabs(y1 - edge_y2) < ERR &&
            fabs(x2 - edge_x1) < ERR && fabs(y2 - edge_y1) < ERR) {
            return(edge);
        }
    }
    
    return(4); // TODO: replace somehow
}

// NOTE: perform a flip, old triangles are overwritten
static void
flip(f32 *triangles, s32 *neighbours, 
     u32 triangle_index_1, u32 edge_1, 
     u32 triangle_index_2, u32 edge_2)
{
    // triangle one
    f32 x1_1 = triangles[triangle_index_1 + 2 * edge_1 + 0];
    f32 y1_1 = triangles[triangle_index_1 + 2 * edge_1 + 1];
    
    f32 x2_1 = triangles[triangle_index_1 + (2 * edge_1 + 2) % 6];
    f32 y2_1 = triangles[triangle_index_1 + (2 * edge_1 + 3) % 6];
    
    f32 opposite_x1 = triangles[triangle_index_1 + (2 * edge_1 + 4) % 6];
    f32 opposite_y1 = triangles[triangle_index_1 + (2 * edge_1 + 5) % 6];
    
    u32 neighbour_index_1 = triangle_index_1 / 2;
    u32 triangle_number_1 = triangle_index_1 / 6;
    // =======
    
    // triangle two
    f32 opposite_x2 = triangles[triangle_index_2 + (2 * edge_2 + 4) % 6];
    f32 opposite_y2 = triangles[triangle_index_2 + (2 * edge_2 + 5) % 6];
    
    u32 neighbour_index_2 = triangle_index_2 / 2;
    u32 triangle_number_2 = triangle_index_2 / 6;
    
    f32 x1_2 = triangles[triangle_index_2 + 2 * edge_2 + 0];
    f32 y1_2 = triangles[triangle_index_2 + 2 * edge_2 + 1];
    
    // NOTE: edge_2 is the same as edge_1, not reverse of edge_1
    bool same_edge_order = (fabs(x1_2 - x1_1) < ERR) && (fabs(y1_2 - y1_1) < ERR);
    // =======
    
    s32 old_neighbours_1[] = {
        neighbours[neighbour_index_1 + 0],
        neighbours[neighbour_index_1 + 1],
        neighbours[neighbour_index_1 + 2]
    };
    
    s32 old_neighbours_2[] = {
        neighbours[neighbour_index_2 + 0],
        neighbours[neighbour_index_2 + 1],
        neighbours[neighbour_index_2 + 2]
    };
    
    {
        // overwrite triangle one
        triangles[triangle_index_1 + 0] = opposite_x1;
        triangles[triangle_index_1 + 1] = opposite_y1;
        
        triangles[triangle_index_1 + 2] = x1_1;
        triangles[triangle_index_1 + 3] = y1_1;
        
        triangles[triangle_index_1 + 4] = opposite_x2;
        triangles[triangle_index_1 + 5] = opposite_y2;
        
        neighbours[neighbour_index_1 + 0] = old_neighbours_1[(edge_1 + 2) % 3];
        neighbours[neighbour_index_1 + 1] = old_neighbours_2[(edge_2 + (same_edge_order ? 2 : 1)) % 3];
        neighbours[neighbour_index_1 + 2] = triangle_number_2;
    }
    
    {
        // overwrite triangle two
        triangles[triangle_index_2 + 0] = opposite_x1;
        triangles[triangle_index_2 + 1] = opposite_y1;
        
        triangles[triangle_index_2 + 2] = x2_1;
        triangles[triangle_index_2 + 3] = y2_1;
        
        triangles[triangle_index_2 + 4] = opposite_x2;
        triangles[triangle_index_2 + 5] = opposite_y2;
        
        neighbours[neighbour_index_2 + 0] = old_neighbours_1[(edge_1 + 1) % 3];
        neighbours[neighbour_index_2 + 1] = old_neighbours_2[(edge_2 + (same_edge_order ? 1 : 2)) % 3];
        neighbours[neighbour_index_2 + 2] = triangle_number_1;
    }
    
    // NOTE: update neighbour's neighbours, we update 4 triangles here
    
    for (u32 edge = 0; edge < 2; ++edge) {
        s32 neighbour_triangle_number = neighbours[triangle_number_1 * 3 + edge]; 
        if (neighbour_triangle_number != -1) {
            s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                            triangle_number_1 * 6, edge);
            neighbours[neighbour_triangle_number * 3 + neighbour_edge] = triangle_number_1;
        }
    }
    
    for (u32 edge = 0; edge < 2; ++edge) {
        s32 neighbour_triangle_number = neighbours[triangle_number_2 * 3 + edge]; 
        if (neighbour_triangle_number != -1) {
            s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                            triangle_number_2 * 6, edge);
            neighbours[neighbour_triangle_number * 3 + neighbour_edge] = triangle_number_2;
        }
    }
}

#if 0
// NOTE: do two line segments intersect?
// https://stackoverflow.com/a/565282
static bool
segments_intersect(f32 ax, f32 ay, f32 bx, f32 by, f32 cx, f32 cy, f32 dx, f32 dy)
{
    f32 ca_x = cx - ax;
    f32 ca_y = cy - ay;
    
    f32 ba_x = bx - ax;
    f32 ba_y = by - ay;
    
    f32 dc_x = dx - cx;
    f32 dc_y = dy - cy;
    
    f32 ca_x_ba = ca_x * ba_y - ca_y * ba_x;
    f32 ca_x_dc = ca_x * dc_y - ca_y * dc_x;
    f32 ba_x_dc = ba_x * dc_y - ba_y * dc_x;
    
    float p = 1.0f / ba_x_dc;
    float t = ca_x_dc * p;
    float u = ca_x_ba * p;
    
    return((t >= 0.0f) && (t <= 1.0f) && (u >= 0.0f) && (u <= 1.0f));
}
#endif

// NOTE: points (ax, ay) and (bx, by) are on the opposite sides
// of the line (x1, y1) - (x2, y2)
// https://math.stackexchange.com/a/162733
static bool
opposite_of_line(f32 ax, f32 ay, f32 bx, f32 by, f32 x1, f32 y1, f32 x2, f32 y2)
{
    return(((y1 - y2) * (ax - x1) + (x2 - x1) * (ay - y1)) * 
           ((y1 - y2) * (bx - x1) + (x2 - x1) * (by - y1)) < -1.0f * ERR);
}

static bool
close_to_line(f32 px, f32 py, f32 ax, f32 ay, f32 bx, f32 by)
{
    // TODO: faster!
    f32 pa_x = ax - px;
    f32 pa_y = ay - py;
    f32 pb_x = bx - px;
    f32 pb_y = by - py;
    
    f32 cos_angle = pa_x * pb_x + pa_y * pb_y; // (pa * pb)
    
    return(fabs(cos_angle) < EDGE_TOLERANCE);
}

// NOTE: is a point inside a 2D triangle?
// https://stackoverflow.com/a/14382692
static bool
inside(f32 *triangles, u32 triangle_index, f32 x, f32 y)
{
    f32 x1 = triangles[triangle_index + 0];
    f32 y1 = triangles[triangle_index + 1];
    
    f32 x2 = triangles[triangle_index + 2];
    f32 y2 = triangles[triangle_index + 3];
    
    f32 x3 = triangles[triangle_index + 4];
    f32 y3 = triangles[triangle_index + 5];
    
    f32 A = 0.5f * (-1.0f * y2 * x3 + y1 * (x3 - x2) + x1 * (y2 - y3) + x2 * y3);
    f32 a_sign = signf(A);
    
    f32 s = (y1 * x3 - x1 * y3 + (y3 - y1) * x + (x1 - x3) * y) * a_sign;
    f32 t = (x1 * y2 - y1 * x2 + (y1 - y2) * x + (x2 - x1) * y) * a_sign;
    
    //printf("%f %f %f %f\n", s, t,  (s + t), 2.0f * A * a_sign);
    
    return(s > 0.0f && t > 0.0f && (s + t) < 2.0f * A * a_sign);
}

static bool
delone_condition(f32 *triangles, u32 triangle_index, u32 opposite_vertex_offset, f32 x, f32 y)
{
    // NOTE: (x2, y2) MUST be the vertex opposite to (x, y)
    f32 x2 = triangles[triangle_index + opposite_vertex_offset + 0];
    f32 y2 = triangles[triangle_index + opposite_vertex_offset + 1];
    
    f32 x1 = triangles[triangle_index + (opposite_vertex_offset + 2) % 6];
    f32 y1 = triangles[triangle_index + (opposite_vertex_offset + 3) % 6];
    
    f32 x3 = triangles[triangle_index + (opposite_vertex_offset + 4) % 6];
    f32 y3 = triangles[triangle_index + (opposite_vertex_offset + 5) % 6];
    
#if 1
    f32 b = (x2 - x1) * (x2 - x3) + (y2 - y1) * (y2 - y3);
    f32 c = (x - x1) * (x - x3) + (y - y1) * (y - y3);
    
    if (b < -1e-4 && c < -1e-4) {
        return(false);
    } else if (b >= 1e-4 && c >= 1e-4) {
        return(true);
    }
    
    f32 a = fabsf((x - x1) * (y - y3) - (x - x3) * (y - y1));
    f32 d = fabsf((x2 - x1) * (y2 - y3) - (x2 - x3) * (y2 - y1));
    
    return(a * b + c * d >= ERR);
#else
    f32 ax = x1 - x2;
    f32 ay = y1 - y2;
    f32 bx = x3 - x2;
    f32 by = y3 - y2;
    f32 dot = ax * bx + ay * by;
    f32 cross = ax * by - ay * bx;
    f32 alpha = fabs(atan2f(cross, dot));
    
    ax = x1 - x;
    ay = y1 - y;
    bx = x3 - x;
    by = y3 - y;
    dot = ax * bx + ay * by;
    cross = ax * by - ay * bx;
    f32 beta = fabs(atan2f(cross, dot));
    
    return(alpha + beta <= PI);
#endif
}

static void
split_edge(f32 *triangles, s32 *neighbours, u32 triangle_number, u32 edge_1, u32 count, f32 x, f32 y)
{
    u32 neighbour_number = neighbours[triangle_number * 3 + edge_1];
    u32 old_triangle_index_1 = triangle_number * 6; 
    u32 old_triangle_index_2 = neighbour_number * 6;
    u32 edge_2 = edge_index(triangles, old_triangle_index_2, old_triangle_index_1, edge_1);
    u32 new_triangle_number = count;
    u32 head = count * 6;
    
    // old triangle one
    f32 x1_1 = triangles[old_triangle_index_1 + 2 * edge_1 + 0];
    f32 y1_1 = triangles[old_triangle_index_1 + 2 * edge_1 + 1];
    
    f32 x2_1 = triangles[old_triangle_index_1 + (2 * edge_1 + 2) % 6];
    f32 y2_1 = triangles[old_triangle_index_1 + (2 * edge_1 + 3) % 6];
    
    f32 opposite_x1 = triangles[old_triangle_index_1 + (2 * edge_1 + 4) % 6];
    f32 opposite_y1 = triangles[old_triangle_index_1 + (2 * edge_1 + 5) % 6];
    
    // old triangle two
    f32 opposite_x2 = triangles[old_triangle_index_2 + (2 * edge_2 + 4) % 6];
    f32 opposite_y2 = triangles[old_triangle_index_2 + (2 * edge_2 + 5) % 6];
    
    f32 x1_2 = triangles[old_triangle_index_1 + 2 * edge_2 + 0];
    f32 y1_2 = triangles[old_triangle_index_1 + 2 * edge_2 + 1];
    
    bool same_edge_order = (fabs(x1_2 - x1_1) < ERR) && (fabs(y1_2 - y1_1) < ERR);
    
    s32 old_neighbours_1[] = {
        neighbours[triangle_number * 3 + 0],
        neighbours[triangle_number * 3 + 1],
        neighbours[triangle_number * 3 + 2]
    };
    
    s32 old_neighbours_2[] = {
        neighbours[neighbour_number * 3 + 0],
        neighbours[neighbour_number * 3 + 1],
        neighbours[neighbour_number * 3 + 2]
    };
    
    {
        // triangle one (replace)
        triangles[old_triangle_index_1 + 0] = x1_1;
        triangles[old_triangle_index_1 + 1] = y1_1;
        
        triangles[old_triangle_index_1 + 2] = x;
        triangles[old_triangle_index_1 + 3] = y;
        
        triangles[old_triangle_index_1 + 4] = opposite_x1;
        triangles[old_triangle_index_1 + 5] = opposite_y1;
        
        u32 old_neighbour_index_1 = triangle_number * 3;
        neighbours[old_neighbour_index_1 + 0] = new_triangle_number;
        neighbours[old_neighbour_index_1 + 1] = neighbour_number;
        neighbours[old_neighbour_index_1 + 2] = old_neighbours_1[(edge_1 + 2) % 3];
    }
    
    {
        // triangle two (replace)
        triangles[old_triangle_index_2 + 0] = x2_1;
        triangles[old_triangle_index_2 + 1] = y2_1;
        
        triangles[old_triangle_index_2 + 2] = x;
        triangles[old_triangle_index_2 + 3] = y;
        
        triangles[old_triangle_index_2 + 4] = opposite_x1;
        triangles[old_triangle_index_2 + 5] = opposite_y1;
        
        u32 old_neighbour_index_2 = neighbour_number * 3;
        neighbours[old_neighbour_index_2 + 0] = new_triangle_number + 1;
        neighbours[old_neighbour_index_2 + 1] = triangle_number;
        neighbours[old_neighbour_index_2 + 2] = old_neighbours_1[(edge_1 + 1) % 3];
    }
    
    {
        // triangle three (new)
        triangles[head++] = x1_1;
        triangles[head++] = y1_1;
        
        triangles[head++] = x;
        triangles[head++] = y;
        
        triangles[head++] = opposite_x2;
        triangles[head++] = opposite_y2;
        
        u32 new_neighbour_index_1 = new_triangle_number * 3;
        neighbours[new_neighbour_index_1 + 0] = triangle_number;
        neighbours[new_neighbour_index_1 + 1] = new_triangle_number + 1;
        neighbours[new_neighbour_index_1 + 2] = old_neighbours_2[(edge_2 + (same_edge_order ? 2 : 1)) % 3];
    }
    
    {
        // triangle four (new)
        triangles[head++] = x2_1;
        triangles[head++] = y2_1;
        
        triangles[head++] = x;
        triangles[head++] = y;
        
        triangles[head++] = opposite_x2;
        triangles[head++] = opposite_y2;
        
        u32 new_neighbour_index_2 = (new_triangle_number + 1) * 3;
        neighbours[new_neighbour_index_2 + 0] = neighbour_number;
        neighbours[new_neighbour_index_2 + 1] = new_triangle_number;
        neighbours[new_neighbour_index_2 + 2] = old_neighbours_2[(edge_2 + (same_edge_order ? 1 : 2)) % 3];
    }
    
    // outer triangle one
    s32 neighbour_triangle_number = neighbours[triangle_number * 3 + 2]; 
    if (neighbour_triangle_number != -1) {
        s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                        triangle_number * 6, 2);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = triangle_number;
    }
    // =======
    
    // outer triangle two
    neighbour_triangle_number = neighbours[neighbour_number * 3 + 2]; 
    if (neighbour_triangle_number != -1) {
        s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                        neighbour_number * 6, 2);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = neighbour_number;
    }
    // =======
    
    // outer triangle three
    neighbour_triangle_number = neighbours[new_triangle_number * 3 + 2];
    if (neighbour_triangle_number != -1) {
        s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                        new_triangle_number * 6, 2);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = new_triangle_number;
    }
    // =======
    
    // outer triangle four
    neighbour_triangle_number = neighbours[(new_triangle_number + 1) * 3 + 2];
    if (neighbour_triangle_number != -1) {
        s32 neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                        (new_triangle_number + 1) * 6, 2);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = new_triangle_number + 1;
    }
    // =======
}

static u32
localize(f32 *triangles, s32 *neighbours, f32 x, f32 y, u32 count,
         u32 start, enum LOCALIZE_ACTION *action, bool readonly)
{
    s32 current_triangle = start;
    
    f32 opposite_x, opposite_y;
    f32 edge_x1, edge_y1, edge_x2, edge_y2;
    
    while (!inside(triangles, current_triangle * 6, x, y)) {
        ASSERT(current_triangle != -1);
        
        // NOTE: find an edge or edges to which the point is 'close enough'
        u32 total_close = 0;
        u32 last_close = 0;
        for (u32 edge = 0; edge < 3; ++edge) {
            edge_x1 = triangles[current_triangle * 6 + edge * 2 + 0];
            edge_y1 = triangles[current_triangle * 6 + edge * 2 + 1];
            
            edge_x2 = triangles[current_triangle * 6 + (edge * 2 + 2) % 6];
            edge_y2 = triangles[current_triangle * 6 + (edge * 2 + 3) % 6];
            
            if (close_to_line(x, y, edge_x1, edge_y1, edge_x2, edge_y2)) {
                ++total_close;
                last_close = edge;
            }
        }
        
        if (total_close > 1) {
            *action = DISCARD;
            if (!readonly) {
                return(0);
            } else {
                return(current_triangle);
            }
        }
        
        if (total_close == 1) {
            if (!readonly) {
                split_edge(triangles, neighbours, current_triangle, last_close, count, x, y);
                *action = SPLIT;
            }
            return(current_triangle);
        }
        
        // NOTE: find an edge for which (x, y) lies on the opposite side 
        // of this edge as opposed to the 'third' point of triangle 
        // ('third' meaning the one not included in the edge)
        for (u32 edge = 0; edge < 3; ++edge) {
            edge_x1 = triangles[current_triangle * 6 + edge * 2 + 0];
            edge_y1 = triangles[current_triangle * 6 + edge * 2 + 1];
            
            edge_x2 = triangles[current_triangle * 6 + (edge * 2 + 2) % 6];
            edge_y2 = triangles[current_triangle * 6 + (edge * 2 + 3) % 6];
            
            opposite_x = triangles[current_triangle * 6 + (edge * 2 + 4) % 6];
            opposite_y = triangles[current_triangle * 6 + (edge * 2 + 5) % 6];
            
            s32 opposite = opposite_of_line(x, y, opposite_x, opposite_y, edge_x1, edge_y1, edge_x2, edge_y2);
            if (opposite) {
                current_triangle = neighbours[current_triangle * 3 + edge];
                break;
            }
        }
    }
    
    *action = LOCAL;
    
    return(current_triangle);
}

// NOTE: split triangle in three, replacing the original one
static void
split_triangle(f32 *triangles, s32 *neighbours, u32 triangle_index, u32 head, f32 x, f32 y)
{
    u32 old_triangle_number = triangle_index / 6; // triangle-to-be-split's number
    u32 head_triangle_number = head / 6;          // number of the first added triangle (the new one)
    u32 neighbour_index = triangle_index / 2;     // where neighbours of the old triangle are written
    u32 neighbour_head = head / 2;                // where to write the new neighbours
    
    // NOTE: add new vertices
    {
        // triangle one (new)
        triangles[head++] = x;
        triangles[head++] = y;
        
        triangles[head++] = triangles[triangle_index + 0];
        triangles[head++] = triangles[triangle_index + 1];
        
        triangles[head++] = triangles[triangle_index + 2];
        triangles[head++] = triangles[triangle_index + 3];
        
        neighbours[neighbour_head++] = old_triangle_number;
        neighbours[neighbour_head++] = neighbours[neighbour_index];
        neighbours[neighbour_head++] = head_triangle_number + 1;
    }
    
    {
        // triangle two (new)
        triangles[head++] = x;
        triangles[head++] = y;
        
        triangles[head++] = triangles[triangle_index + 2];
        triangles[head++] = triangles[triangle_index + 3];
        
        triangles[head++] = triangles[triangle_index + 4];
        triangles[head++] = triangles[triangle_index + 5];
        
        neighbours[neighbour_head++] = head_triangle_number;
        neighbours[neighbour_head++] = neighbours[neighbour_index + 1];
        neighbours[neighbour_head++] = old_triangle_number;
    }
    
    {
        // triangle three (overwrite)
        triangles[triangle_index + 2] = x;
        triangles[triangle_index + 3] = y;
        
        neighbours[neighbour_index + 0] = head_triangle_number;
        neighbours[neighbour_index + 1] = head_triangle_number + 1;
        neighbours[neighbour_index + 2] = neighbours[neighbour_index + 2];
    }
    
    // NOTE: update neighbour's neighbours. We only need to update
    // 'outer' neighbours, as inner triangles are already up to date
    // [ont2]
    
    // NOTE, IMPORTNANT: outer edge indices depend on the triangle splitting
    // code above (namely on the order of vertices in three new triangles)!!!
    
    u32 our_edge;
    u32 neighbour_edge;            // edge index of outer neighbours
    s32 neighbour_triangle_number; // outer neighbour number
    
    // outer triangle one
    our_edge = 1; // <---- THIS '1' depends on the code above!
    neighbour_triangle_number = neighbours[head_triangle_number * 3 + our_edge]; 
    if (neighbour_triangle_number != -1) {
        neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                    head_triangle_number * 6, our_edge);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = head_triangle_number;
    }
    // =======
    
    // outer triangle two
    our_edge = 1; // <---- THIS '1' depends on the code above!
    neighbour_triangle_number = neighbours[(head_triangle_number + 1) * 3 + our_edge]; 
    if (neighbour_triangle_number != -1) {
        neighbour_edge =  edge_index(triangles, neighbour_triangle_number * 6,
                                     (head_triangle_number + 1) * 6, our_edge);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = head_triangle_number + 1;
    }
    // =======
    
    // outer triangle three
    our_edge = 2; // <---- THIS '2' depends on the code above!
    neighbour_triangle_number = neighbours[old_triangle_number * 3 + our_edge]; 
    if (neighbour_triangle_number != -1) {
        neighbour_edge = edge_index(triangles, neighbour_triangle_number * 6,
                                    old_triangle_number * 6, our_edge);
        neighbours[neighbour_triangle_number * 3 + neighbour_edge] = old_triangle_number;
    }
    // =======
}

static void
add_superstructure(f32 *triangles, s32 *neighbours, f32 x, f32 y)
{
    // NOTE: triangle one
    triangles[0] = -0.1f * WIDTH;
    triangles[1] = 1.1f * HEIGHT;
    
    triangles[2] = 1.1f * WIDTH;
    triangles[3] = 1.1f * HEIGHT;
    
    triangles[4] = x;
    triangles[5] = y;
    
    neighbours[1] = 1;
    neighbours[2] = 3;
    
    // NOTE: triangle two
    triangles[6] = triangles[2];
    triangles[7] = triangles[3];
    
    triangles[8] = 1.1f * WIDTH;
    triangles[9] = -0.1f * HEIGHT;
    
    triangles[10] = x;
    triangles[11] = y;
    
    neighbours[4] = 2;
    neighbours[5] = 0;
    
    // NOTE: triangle three
    triangles[12] = triangles[8];
    triangles[13] = triangles[9];
    
    triangles[14] = -0.1f * WIDTH;
    triangles[15] = -0.1f * HEIGHT;
    
    triangles[16] = x;
    triangles[17] = y;
    
    neighbours[7] = 3;
    neighbours[8] = 1;
    
    // NOTE: triangle four
    triangles[18] = triangles[14];
    triangles[19] = triangles[15];
    
    triangles[20] = triangles[0];
    triangles[21] = triangles[1];
    
    triangles[22] = x;
    triangles[23] = y;
    
    neighbours[10] = 0;
    neighbours[11] = 2;
}

static f32 *
triangulate(f32 *points, s32 *neighbours, u32 *triangle_count)
{
    ASSERT(N >= 1);
    
    u32 cell_count = COEFF * powf(N, 3.0f / 8.0f);
    u32 *STATIC_CACHE = calloc(cell_count * cell_count, sizeof(u32));
    
    f32 x, y;
    u32 found_triangle;
    enum LOCALIZE_ACTION action;
    
    u32 count = 0;
    f32 *triangles = calloc(6 * (4 + 2 * N), sizeof(f32));
    
    add_superstructure(triangles, neighbours, points[0], points[1]);
    count += 4;
    
    u32 flip_stack[FLIPSTACK_SIZE];
    u32 flipstack_head = 0;
    
    for (u32 i = 2; i < N * 2; i += 2) {
        x = points[i];
        y = points[i + 1];
        
        u32 cache_index_x = x / WIDTH * cell_count;
        u32 cache_index_y = y / HEIGHT * cell_count;
        u32 cache_index = cache_index_y * cell_count + cache_index_x;
        
        found_triangle = localize(triangles, neighbours, x, y, count, 
                                  STATIC_CACHE[cache_index], &action, false);
        
        // NOTE: creates four triangles if found_triangle < 0
        // (replaces two and creates two new ones). This means that
        // count += 2 is correct still
        if (action == LOCAL) {
            split_triangle(triangles, neighbours, found_triangle * 6, count * 6, x, y);
            STATIC_CACHE[cache_index] = found_triangle;
            count += 2;
        } else if (action == SPLIT) {
            STATIC_CACHE[cache_index] = found_triangle;
            count += 2;
        }
        
        // NOTE: push three new triangles to the stack
        ASSERT(flipstack_head == 0);
        
        flip_stack[flipstack_head++] = found_triangle;
        flip_stack[flipstack_head++] = count - 1;
        flip_stack[flipstack_head++] = count - 2;
        
        while (flipstack_head != 0) {
            u32 pop_triangle = flip_stack[--flipstack_head];
            
            for (u32 edge = 0; edge < 3; ++edge) {
                s32 pop_neighbour_number = neighbours[pop_triangle * 3 + edge];
                if (pop_neighbour_number != -1) {
                    u32 adjacent_edge = edge_index(triangles, pop_neighbour_number * 6, pop_triangle * 6, edge);
                    f32 opposite_x;
                    f32 opposite_y;
                    
                    opposite_x = triangles[pop_neighbour_number * 6 + (adjacent_edge * 2 + 4) % 6];
                    opposite_y = triangles[pop_neighbour_number * 6 + (adjacent_edge * 2 + 5) % 6];
                    
                    if (!delone_condition(triangles, pop_triangle * 6, 2 * ((edge + 2) % 3), opposite_x, opposite_y)) {
                        flip(triangles, neighbours, pop_triangle * 6, edge, pop_neighbour_number * 6,adjacent_edge);
                        flip_stack[flipstack_head++] = pop_triangle;
                        flip_stack[flipstack_head++] = pop_neighbour_number;
                    }
                }
            } // end of find edge we intersect
        }
    }
    
    *triangle_count = count;
    
    return(triangles);
}

s32
main(void)
{
    // NOTE: init random series
    u64 SEED = 1555667014;
    srand(SEED);
    xor_state = rand();
    
    // NOTE: [x0, y0, x1, y1, ...]
    f32 *points = malloc(2 * N * sizeof(f32));
    for (u32 i = 0; i < N * 2; i += 2) {
        points[i + 0] = rand_f32() * WIDTH;
        points[i + 1] = rand_f32() * HEIGHT;
    }
    
    u32 triangle_count = 0;
    s32 *neighbours;
    posix_memalign((void **) &neighbours, 8 * sizeof(void *), 3 * (4 + 2 * N) * sizeof(s32));
    
    const u32 ITERATIONS = 100;
    u64 begin = rdtscp();
    f32 *triangles;
    for (u32 i = 0; i < ITERATIONS; ++i) {
        memset(neighbours, -1, 3 * (4 + 2 * N) * sizeof(s32));
        triangle_count = 0;
        triangles = triangulate(points, neighbours, &triangle_count);
    }
    u64 end = rdtscp();
    
    u64 cycles = end - begin;
    printf("[INFO] Average time ~ %.1f msec\n", 
           (f64) cycles / CYCLES_PER_MSEC_CPU / ITERATIONS);
#if VIS
    struct GL init = init_opengl();
    
    glBindBuffer(GL_ARRAY_BUFFER, init.VBO);
    glBufferData(GL_ARRAY_BUFFER, 6 * triangle_count * sizeof(f32), triangles, GL_STATIC_DRAW);
    
    glUseProgram(init.program);
    
    while (!glfwWindowShouldClose(init.window)) {
        glfwPollEvents();
        glClear(GL_COLOR_BUFFER_BIT);
        
        glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
        glDrawArrays(GL_TRIANGLES, 0, triangle_count * 3);
        
        glfwSwapBuffers(init.window);
    }
    
    glfwTerminate();
#endif
    
    
    return(0);
}