#pragma once

#include <string>

struct Point
{
	Point() : x(0.0f), y(0.0f) { }
	Point(float _x, float _y) : x(_x), y(_y) { }

	std::string to_string() const
	{
		return "(" + std::to_string(x) + ", " + std::to_string(y) + ")";
	}

	float x;
	float y;
};

std::vector<int> triangulate(std::vector<Point>& points);

#include "point.h"

#include <vector>
#include <iostream>
#include <algorithm>
#include <float.h>

struct TriangleCircle
{
    // triangle point indices
    int i;
    int j;
    int k;

    // circle center
    float x;
    float y;

    // the square of the circle radius
    float radius_sq;

    std::string to_string() const
    {
        return "<" + std::to_string(i) + ", " + std::to_string(j) + ", " + std::to_string(k) + "|" + std::to_string(x) + ", " + std::to_string(y) + "|" + std::to_string(radius_sq) + ">";
    }
};

std::vector<Point> build_supertriangle(const std::vector<Point> &points)
{
    float x_min = FLT_MAX;
    float y_min = FLT_MAX;
    float x_max = FLT_MIN;
    float y_max = FLT_MIN;

    for (size_t i = 0; i < points.size(); i++)
    {
        const Point p = points[i];
        if (p.x < x_min) { x_min = p.x; }
        if (p.x > x_max) { x_max = p.x; }
        if (p.y < y_min) { y_min = p.y; }
        if (p.y > y_max) { y_max = p.y; }
    }

    float dx = x_max - x_min;
    float dy = y_max - y_min;

    float d_max = std::max(dx, dy);
    float x_mid = x_min + dx * 0.5f;
    float y_mid = y_min + dy * 0.5f;

    return { 
        Point(x_mid - 20.0f * d_max, y_mid - d_max),
        Point(x_mid, y_mid + 20.0f * d_max),
        Point(x_mid + 20.0f * d_max, y_mid - d_max) };
}

TriangleCircle circumcircle(const std::vector<Point> &points, int i, int j, int k)
{
    const Point point_i = points[i];
    const Point point_j = points[j];
    const Point point_k = points[k];

    float x_1 = point_i.x;
    float y_1 = point_i.y;

    float x_2 = point_j.x;
    float y_2 = point_j.y;

    float x_3 = point_k.x;
    float y_3 = point_k.y;

    float y1_y2 = std::abs(y_1 - y_2);
    float y2_y3 = std::abs(y_2 - y_3);

    float center_x = 0.0f;
    float center_y = 0.0f;
    float m_1 = 0.0f;
    float m_2 = 0.0f;
    float mx_1 = 0.0f;
    float mx_2 = 0.0f;
    float my_1 = 0.0f;
    float my_2 = 0.0f;

    float EPSILON = 0.00001f;

    if (y1_y2 < EPSILON)
    {
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_2 = (x_2 + x_3) / 2.0f;
        my_2 = (y_2 + y_3) / 2.0f;
        center_x = (x_2 + x_1) / 2.0f;
        center_y = m_2 * (center_x - mx_2) + my_2;
    }
    else if (y2_y3 < EPSILON)
    {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        mx_1 = (x_1 + x_2) / 2.0f;
        my_1 = (y_1 + y_2) / 2.0f;
        center_x = (x_3 + x_2) / 2.0f;
        center_y = m_1 * (center_x - mx_1) + my_1;
    }
    else
    {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_1 = (x_1 + x_2) / 2.0f;
        mx_2 = (x_2 + x_3) / 2.0f;
        my_1 = (y_1 + y_2) / 2.0f;
        my_2 = (y_2 + y_3) / 2.0f;
        center_x = (m_1 * mx_1 - m_2 * mx_2 + my_2 - my_1) / (m_1 - m_2);
        center_y = y1_y2 > y2_y3 ? (m_1 * (center_x - mx_1) + my_1) : (m_2 * (center_x - mx_2) + my_2);
    }

    float dx = x_2 - center_x;
    float dy = y_2 - center_y;

    return TriangleCircle{ i, j, k, center_x, center_y, dx*dx + dy*dy};
}

void remove_duplicates(std::vector<int> &edges)
{
    int j = edges.size();
    while (j >= 2)
    {
        int b = edges[--j];
        int a = edges[--j];

        int i = j;
        while (i >= 2)
        {
            int n = edges[--i];
            int m = edges[--i];

            if ((a == m && b == n) || (a == n && b == m))
            {
                edges.erase(std::next(edges.begin(), j), std::next(edges.begin(), j + 2));
                edges.erase(std::next(edges.begin(), i), std::next(edges.begin(), i + 2));

                j -= 2;
                break;
            }
        }
    }
}

std::vector<int> triangulate(std::vector<Point> &points)
{
    size_t points_count = points.size();

    // if the nmber is less than 3, nothing to do
    if (points_count < 3)
    {
        return std::vector<int>(0);
    }

    // sort points by x coordinate
    // create array with indices 0, 1, 2, ...
    std::vector<int> indices(points_count);
    for (size_t  i = 0; i < points_count; i++)
    {
        indices[i] = i;
    }

    // sort it by using input points
    std::sort(indices.begin(), indices.end(), [&points](int a, int b) { return points[a].x < points[b].x; });

    const std::vector<Point> st = build_supertriangle(points);
    // add points of the super triangle into input array
    points.insert(points.end(), st.begin(), st.end());

    // create buffers
    std::vector<TriangleCircle> open_list;
    std::vector<TriangleCircle> closed_list;
    std::vector<int> edges_list;

    // init open buffer by super triangle
    open_list.push_back(circumcircle(points, points_count, points_count + 1, points_count + 2));
    
    for (size_t i = 0; i < points_count; i++)
    {
        int c = indices[i];
        Point p = points[c];
        edges_list.clear();

        for (int j = open_list.size() - 1; j >= 0; j--)
        {
            TriangleCircle t = open_list[j];
            
            float dx = p.x - t.x;
            if (dx > 0.0f && dx*dx > t.radius_sq)
            {
                open_list.erase(open_list.begin() + j);
                closed_list.push_back(t);
                continue;
            }

            float dy = p.y - t.y;
            if (dx*dx + dy*dy - t.radius_sq > 0.00001f)
            {
                continue;
            }

            open_list.erase(open_list.begin() + j);
            edges_list.insert(edges_list.end(), { t.i, t.j, t.j, t.k, t.k, t.i });
        }

        // delete edges with two triangles
        remove_duplicates(edges_list);

        // add triangles for the remain edges to the open list
        int k = edges_list.size();
        while (k >= 2)
        {
            int b = edges_list[--k];
            int a = edges_list[--k];
            open_list.push_back(circumcircle(points, a, b, c));
        }
    }

    // add to the closest list all triangles from the remain open list
    for (size_t i = 0; i < open_list.size(); i++)
    {
        closed_list.push_back(open_list[i]);
    }

    open_list.clear();
    edges_list.clear();

    // form the output array
    std::vector<int> triangles;
    for (size_t i = 0; i < closed_list.size(); i++)
    {
        TriangleCircle t = closed_list[i];
        if (t.i < points_count && t.j < points_count && t.k < points_count)
        {
            triangles.insert(triangles.end(), { t.i, t.j, t.k });
        }
    }

    closed_list.clear();

    return triangles;
}

#pragma once
#include "point.h"
#include <vector>

struct AABB
{
	float x_min;
	float y_min;
	float x_max;
	float y_max;

	std::string to_string()
	{
		return "|" + std::to_string(x_min) + ", " + std::to_string(y_min) + ", " + std::to_string(x_max) + ", " + std::to_string(y_max) + "|";
	}
};

class Trinangle
{
public:
	Trinangle(const std::vector<Point>& vertices);
	~Trinangle() {};
	AABB get_aabb();
	Point get_center();
	bool is_point_inside(const Point &point);
	std::string to_string();

	Point get_a();
	Point get_b();
	Point get_c();

private:
	Point a;
	Point b;
	Point c;

	AABB aabb;
	Point center;
};

class BVHNode
{
public:
	BVHNode(const std::vector<Trinangle*> &triangles);
	~BVHNode();

	AABB get_aabb();
	bool is_inside_aabb(const Point &point);
	Trinangle* sample(const Point &point);

private:
	Trinangle* triangle;
	BVHNode* left_node;
	BVHNode* right_node;
	AABB aabb;
};

#include <iostream>
#include <float.h>

#include "bvh.h"

Trinangle::Trinangle(const std::vector<Point> &vertices)
{
	a = vertices[0];
	b = vertices[1];
	c = vertices[2];

	float x_min = FLT_MAX;
	float y_min = FLT_MAX;
	float x_max = FLT_MIN;
	float y_max = FLT_MIN;

	float x_accum = 0.0f;
	float y_accum = 0.0f;

	for (size_t i = 0; i < vertices.size(); i++)
	{
		Point v = vertices[i];
		if (v.x < x_min) { x_min = v.x; }
		if (v.x > x_max) { x_max = v.x; }
		if (v.y < y_min) { y_min = v.y; }
		if (v.y > y_max) { y_max = v.y; }

		x_accum += v.x;
		y_accum += v.y;
	}

	aabb = AABB{ x_min, y_min, x_max, y_max };
	center = Point(x_accum / vertices.size(), y_accum / vertices.size());
}

AABB Trinangle::get_aabb()
{
	return aabb;
}

Point Trinangle::get_center()
{
	return center;
}

bool Trinangle::is_point_inside(const Point& point)
{
	float as_x = point.x - a.x;
	float as_y = point.y - a.y;

	bool s_ab = ((b.x - a.x) * as_y - (b.y - a.y) * as_x) > 0;

	if (((c.x - a.x) * as_y - (c.y - a.y) * as_x > 0) == s_ab)
	{
		return false;
	}
	if (((c.x - b.x) * (point.y - b.y) - (c.y - b.y) * (point.x - b.x) > 0) != s_ab)
	{
		return false;
	}
	return true;
}

std::string Trinangle::to_string()
{
	return "<" + a.to_string() + ", " + b.to_string() + ", " + c.to_string() + ">";
}

Point Trinangle::get_a()
{
	return a;
}

Point Trinangle::get_b()
{
	return b;
}
Point Trinangle::get_c()
{
	return c;
}

AABB union_aabb(AABB &x, AABB &y)
{
	return AABB {
		std::min(x.x_min, y.x_min), std::min(x.y_min, y.y_min), 
		std::max(x.x_max, y.x_max), std::max(x.y_max, y.y_max) };
}

BVHNode::BVHNode(const std::vector<Trinangle*>& triangles)
{
	triangle = NULL;
	left_node = NULL;
	right_node = NULL;

	int objects_count = triangles.size();
	if (objects_count == 1)
	{
		triangle = triangles[0];
		aabb = triangle->get_aabb();
	}
	else
	{
		float x_median = 0.0f;
		float y_median = 0.0f;

		float x_min = FLT_MAX;
		float x_max = FLT_MIN;
		float y_min = FLT_MAX;
		float y_max = FLT_MIN;

		for (size_t i = 0; i < objects_count; i++)
		{
			Trinangle* t = triangles[i];
			Point c = t->get_center();
			x_median += c.x;
			y_median += c.y;

			if (c.x < x_min) { x_min = c.x; }
			if (c.x > x_max) { x_max = c.x; }
			if (c.y < y_min) { y_min = c.y; }
			if (c.y > y_max) { y_max = c.y; }
		}

		int split_axis = ((x_max - x_min) > (y_max - y_min)) ? 0 : 1;
		float median = split_axis == 0 ? (x_median / objects_count) : (y_median / objects_count);

		std::vector<Trinangle*> left;
		std::vector<Trinangle*> right;

		for (size_t i = 0; i < objects_count; i++)
		{
			Trinangle* t = triangles[i];
			Point c = t->get_center();
			float v = split_axis == 0 ? c.x : c.y;
			if (v < median)
			{
				left.push_back(t);
			}
			else
			{
				right.push_back(t);
			}
		}

		if (left.size() == 0)
		{
			left.push_back(right[right.size() - 1]);
			right.pop_back();
		}

		if (right.size() == 0)
		{
			right.push_back(left[left.size() - 1]);
			left.pop_back();
		}

		left_node = new BVHNode(left);
		right_node = new BVHNode(right);

		AABB left_aabb = left_node->get_aabb();
		AABB right_aabb = right_node->get_aabb();

		aabb = union_aabb(left_aabb, right_aabb);
	}
}

BVHNode::~BVHNode()
{
	if (triangle) { delete triangle; }
	if (left_node) { delete left_node; }
	if (right_node) { delete right_node; }
}

AABB BVHNode::get_aabb()
{
	return aabb;
}

bool BVHNode::is_inside_aabb(const Point& point)
{
	return aabb.x_min < point.x && aabb.y_min < point.y && aabb.x_max > point.x && aabb.y_max > point.y;
}

Trinangle* BVHNode::sample(const Point& point)
{
	if (is_inside_aabb(point))
	{
		if (!triangle && left_node && right_node)
		{
			Trinangle* left_sample = left_node->sample(point);
			Trinangle* right_sample = right_node->sample(point);

			if (!left_sample)
			{
				return right_sample;
			}
			else
			{
				if (!right_sample)
				{
					return left_sample;
				}
				else
				{
					Point l_c = left_sample->get_center();
					float l_dist = (l_c.x - point.x) * (l_c.x - point.x) + (l_c.y - point.y) * (l_c.y - point.y);

					Point r_c = right_sample->get_center();
					float r_dist = (r_c.x - point.x) * (r_c.x - point.x) + (r_c.y - point.y) * (r_c.y - point.y);

					if (l_dist < r_dist)
					{
						return left_sample;
					}
					else
					{
						return right_sample;
					}
				}
			}
		}
		else
		{
			if (triangle && triangle->is_point_inside(point))
			{
				return triangle;
			}
			else
			{
				return NULL;
			}
		}
	}
	else
	{
		return NULL;
	}
}

#include "delaunay.h"
#include "bvh.h"

#include <iostream>

void console_log(const std::vector<int>& array)
{
    std::cout << "[";
    for (size_t i = 0; i < array.size(); i++)
    {
        std::cout << array[i] << (i == array.size() - 1 ? "" : ", ");
    }

    std::cout << "]" << std::endl;
}

int main()
{
    std::vector<Point> points = { 
        Point(-0.5, 3.0), 
        Point(8.0, -3.5), 
        Point(-7.0, 3.0), 
        Point(2.0, -10.0), 
        Point(7.0, 8.0), 
        Point(8.0, 5.5) 
    };
    std::vector<int> triangle_indices = triangulate(points);

    console_log(triangle_indices);

    std::vector<Trinangle*> triangles(triangle_indices.size() / 3);
    for (size_t i = 0; i < triangles.size(); i++)
    {
        int a = triangle_indices[3 * i];
        int b = triangle_indices[3 * i + 1];
        int c = triangle_indices[3 * i + 2];
        std::vector<Point> vertices = { points[a], points[b], points[c] };
        triangles[i] = new Trinangle(vertices);
    }

    BVHNode* tree = new BVHNode(triangles);
    Point p = Point(0.0f, 0.0f);
    Trinangle* sample_result = tree->sample(p);

    if (sample_result)
    {
        std::cout << sample_result->to_string() << std::endl;
    }
    else
    {
        std::cout << "Empty sample" << std::endl;
    }
}

#include <emscripten/bind.h>
#include <emscripten/val.h>

#include "point.h"
#include "delaunay.h"
#include "bvh.h"

using namespace emscripten;

emscripten::val build_triangulation(const emscripten::val &in_coordinates)
{
	// prepare input
	const std::vector<float> coordinates = convertJSArrayToNumberVector<float>(in_coordinates);

	int points_count = coordinates.size() / 2;
	std::vector<Point> points(points_count);
	for (size_t i = 0; i < points_count; i++)
	{
		points[i] = Point(coordinates[2*i], coordinates[2*i + 1]);
	}

	std::vector<int> triangles = triangulate(points);

	// prepare output
	emscripten::val view{ emscripten::typed_memory_view(triangles.size(), triangles.data()) };
	emscripten::val result = emscripten::val::global("Int32Array").new_(triangles.size());
	result.call<void>("set", view);

	return result;
}

class BVHNodeWrapper
{
public:
	BVHNodeWrapper(const emscripten::val &in_coordinates);
	BVHNodeWrapper(const emscripten::val &in_coordinates, const emscripten::val &in_triangles);
	~BVHNodeWrapper();

	emscripten::val sample(float x, float y);

private:
	BVHNode* bvh;
};

BVHNodeWrapper::BVHNodeWrapper(const emscripten::val &in_coordinates)
{
	const std::vector<float> coordinates = convertJSArrayToNumberVector<float>(in_coordinates);

	int points_count = coordinates.size() / 2;
	std::vector<Point> points(points_count);
	for (size_t i = 0; i < points_count; i++)
	{
		points[i] = Point(coordinates[2*i], coordinates[2*i + 1]);
	}

	std::vector<int> trinagle_indices = triangulate(points);
	int trinagles_count = trinagle_indices.size() / 3;
	std::vector<Trinangle*> trinagles(trinagles_count);
	for (size_t i = 0; i < trinagles_count; i++)
	{
		int a = trinagle_indices[3 * i];
		int b = trinagle_indices[3 * i + 1];
		int c = trinagle_indices[3 * i + 2];

		std::vector<Point> vertices = { points[a], points[b], points[c] };
		trinagles[i] = new Trinangle(vertices);
	}

	bvh = new BVHNode(trinagles);
}

BVHNodeWrapper::BVHNodeWrapper(const emscripten::val &in_coordinates, const emscripten::val &in_triangles)
{
	const std::vector<float> coordinates = convertJSArrayToNumberVector<float>(in_coordinates);
	const std::vector<int> triangles = convertJSArrayToNumberVector<int>(in_triangles);

	int triangles_count = triangles.size() / 3;
	std::vector<Trinangle*> triangles_array(triangles_count);
	for (size_t i = 0; i < triangles_count; i++)
	{
		int a = triangles[3 * i];
		int b = triangles[3 * i + 1];
		int c = triangles[3 * i + 2];

		std::vector<Point> vertices = {
			Point(coordinates[2 * a], coordinates[2 * a + 1]),
			Point(coordinates[2 * b], coordinates[2 * b + 1]),
			Point(coordinates[2 * c], coordinates[2 * c + 1]) };
		triangles_array[i] = new Trinangle(vertices);
	}

	bvh = new BVHNode(triangles_array);
}

BVHNodeWrapper::~BVHNodeWrapper()
{
	delete bvh;
}

emscripten::val BVHNodeWrapper::sample(float x, float y)
{
	Point p = Point(x, y);
	Trinangle* s = bvh->sample(p);
	std::vector<float> to_return(0);
	if (s)
	{
		Point a = s->get_a();
		Point b = s->get_b();
		Point c = s->get_c();

		to_return.insert(to_return.end(), { a.x, a.y, b.x, b.y, c.x, c.y });
	}

	emscripten::val view{ emscripten::typed_memory_view(to_return.size(), to_return.data()) };
	emscripten::val result = emscripten::val::global("Float32Array").new_(to_return.size());
	result.call<void>("set", view);

	return result;
}

EMSCRIPTEN_BINDINGS(delaunay_module)
{
	class_<BVHNodeWrapper>("BVHNode")
		.constructor<emscripten::val>()
		.constructor<emscripten::val, emscripten::val>()
		.function("sample", &BVHNodeWrapper::sample);
	function("build_triangulation", &build_triangulation);
}

class BVHNodeWrapper
{
public:
	BVHNodeWrapper(const emscripten::val &in_coordinates);
	BVHNodeWrapper(const emscripten::val &in_coordinates, const emscripten::val &in_triangles);
	~BVHNodeWrapper();

	emscripten::val sample(float x, float y);

private:
	BVHNode* bvh;
};

BVHNodeWrapper::BVHNodeWrapper(const emscripten::val &in_coordinates)
{
	const std::vector<float> coordinates = convertJSArrayToNumberVector<float>(in_coordinates);

	int points_count = coordinates.size() / 2;
	std::vector<Point> points(points_count);
	for (size_t i = 0; i < points_count; i++)
	{
		points[i] = Point(coordinates[2*i], coordinates[2*i + 1]);
	}

	std::vector<int> trinagle_indices = triangulate(points);
	int trinagles_count = trinagle_indices.size() / 3;
	std::vector<Trinangle*> trinagles(trinagles_count);
	for (size_t i = 0; i < trinagles_count; i++)
	{
		int a = trinagle_indices[3 * i];
		int b = trinagle_indices[3 * i + 1];
		int c = trinagle_indices[3 * i + 2];

		std::vector<Point> vertices = { points[a], points[b], points[c] };
		trinagles[i] = new Trinangle(vertices);
	}

	bvh = new BVHNode(trinagles);
}

BVHNodeWrapper::BVHNodeWrapper(const emscripten::val &in_coordinates, const emscripten::val &in_triangles)
{
	const std::vector<float> coordinates = convertJSArrayToNumberVector<float>(in_coordinates);
	const std::vector<int> triangles = convertJSArrayToNumberVector<int>(in_triangles);

	int triangles_count = triangles.size() / 3;
	std::vector<Trinangle*> triangles_array(triangles_count);
	for (size_t i = 0; i < triangles_count; i++)
	{
		int a = triangles[3 * i];
		int b = triangles[3 * i + 1];
		int c = triangles[3 * i + 2];

		std::vector<Point> vertices = {
			Point(coordinates[2 * a], coordinates[2 * a + 1]),
			Point(coordinates[2 * b], coordinates[2 * b + 1]),
			Point(coordinates[2 * c], coordinates[2 * c + 1]) };
		triangles_array[i] = new Trinangle(vertices);
	}

	bvh = new BVHNode(triangles_array);
}

emscripten::val BVHNodeWrapper::sample(float x, float y)
{
	Point p = Point(x, y);
	Trinangle* s = bvh->sample(p);
	std::vector<float> to_return(0);
	if (s)
	{
		Point a = s->get_a();
		Point b = s->get_b();
		Point c = s->get_c();

		to_return.insert(to_return.end(), { a.x, a.y, b.x, b.y, c.x, c.y });
	}

	emscripten::val view{ emscripten::typed_memory_view(to_return.size(), to_return.data()) };
	emscripten::val result = emscripten::val::global("Float32Array").new_(to_return.size());
	result.call<void>("set", view);

	return result;
}

export class Point {
    private m_x: f32;
    private m_y: f32;

    constructor(in_x: f32 = 0.0, in_y: f32 = 0.0) {
        this.m_x = in_x;
        this.m_y = in_y;
    }

    @inline
    x(): f32 {
        return this.m_x;
    }

    @inline
    y(): f32 {
        return this.m_y;
    }

    toString(): string {
        return "(" + this.m_x.toString() + ", " + this.m_y.toString() + ")";
    }
}

import { Point } from "./point";

const EPSILON: f32 = 0.0001;

class PointExt extends Point {
    private m_index: i32;

    constructor(in_x: f32, in_y: f32, in_index: i32) {
        super(in_x, in_y);

        this.m_index = in_index;
    }

    index(): i32 {
        return this.m_index;
    }
}

class TriangleCircle {
    private m_i: i32;
    private m_j: i32;
    private m_k: i32;

    private m_x: f32;
    private m_y: f32;

    private m_radius_sq: f32;

    constructor(in_i: i32, in_j: i32, in_k: i32, in_x: f32, in_y: f32, in_r: f32) {
        this.m_i = in_i;
        this.m_j = in_j;
        this.m_k = in_k;

        this.m_x = in_x;
        this.m_y = in_y;

        this.m_radius_sq = in_r;
    }

    i(): i32 { return this.m_i; }
    j(): i32 { return this.m_j; }
    k(): i32 { return this.m_k; }

    x(): f32 { return this.m_x; }
    y(): f32 { return this.m_y; }

    radius_sq(): f32 { return this.m_radius_sq; }

    toString(): string {
        return "<" + this.m_i.toString() + ", " + this.m_j.toString() + ", " + this.m_k.toString() + "|" + this.m_x.toString() + ", " + this.m_y.toString() + "|" + this.m_radius_sq.toString() + ">";
    }
}

function maximum(a: f32, b: f32): f32 {
    if(a > b) {
        return a;
    } else {
        return b;
    }
}

function absolute(a: f32, b: f32): f32 {
    if(a > b) {
        return a - b;
    } else {
        return b - a;
    }
}

function build_supertriangle(points: StaticArray<Point>): StaticArray<Point> {
    let x_min: f32 = f32.MAX_VALUE;
    let y_min: f32 = f32.MAX_VALUE;
    let x_max: f32 = f32.MIN_VALUE;
    let y_max: f32 = f32.MIN_VALUE;

    for(let i = 0, len = points.length; i < len; i++) {
        const p = points[i];
        if (p.x() < x_min) { x_min = p.x(); }
        if (p.x() > x_max) { x_max = p.x(); }
        if (p.y() < y_min) { y_min = p.y(); }
        if (p.y() > y_max) { y_max = p.y(); }
    }

    const dx: f32 = x_max - x_min;
    const dy: f32 = y_max - y_min;

    const d_max: f32 = maximum(dx, dy);
    const x_mid: f32 = x_min + dx * 0.5;
    const y_mid: f32 = y_min + dy * 0.5;

    const vertices = new StaticArray<Point>(3);
    vertices[0] = new Point(x_mid - 20.0 * d_max, y_mid - d_max);
    vertices[1] = new Point(x_mid, y_mid + 20.0 * d_max);
    vertices[2] = new Point(x_mid + 20.0 * d_max, y_mid - d_max);
    return vertices;
}

function circumcircle(points: StaticArray<Point>, i: i32, j: i32, k: i32): TriangleCircle {
    const point_i = points[i];
    const point_j = points[j];
    const point_k = points[k];

    const x_1 = point_i.x();
    const y_1 = point_i.y();

    const x_2 = point_j.x();
    const y_2 = point_j.y();

    const x_3 = point_k.x();
    const y_3 = point_k.y();

    const y1_y2 = absolute(y_1, y_2);
    const y2_y3 = absolute(y_2, y_3);

    let center_x: f32 = 0.0;
    let center_y: f32 = 0.0;
    let m_1: f32 = 0.0;
    let m_2: f32 = 0.0;
    let mx_1: f32 = 0.0;
    let mx_2: f32 = 0.0;
    let my_1: f32 = 0.0;
    let my_2: f32 = 0.0;

    if (y1_y2 < EPSILON) {
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_2 = (x_2 + x_3) / 2.0;
        my_2 = (y_2 + y_3) / 2.0;
        center_x = (x_2 + x_1) / 2.0;
        center_y = m_2 * (center_x - mx_2) + my_2;
    } else if (y2_y3 < EPSILON)
    {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        mx_1 = (x_1 + x_2) / 2.0;
        my_1 = (y_1 + y_2) / 2.0;
        center_x = (x_3 + x_2) / 2.0;
        center_y = m_1 * (center_x - mx_1) + my_1;
    } else
    {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_1 = (x_1 + x_2) / 2.0;
        mx_2 = (x_2 + x_3) / 2.0;
        my_1 = (y_1 + y_2) / 2.0;
        my_2 = (y_2 + y_3) / 2.0;
        center_x = (m_1 * mx_1 - m_2 * mx_2 + my_2 - my_1) / (m_1 - m_2);
        center_y = y1_y2 > y2_y3 ? (m_1 * (center_x - mx_1) + my_1) : (m_2 * (center_x - mx_2) + my_2);
    }

    const dx = x_2 - center_x;
    const dy = y_2 - center_y;

    return new TriangleCircle(i, j, k, center_x, center_y, dx*dx + dy*dy);
}

// retun new size of the edges array
function remove_duplicates(edges: StaticArray<i32>, in_length: i32): i32 {
    let j = in_length;
    while(j >= 2) {
        const b = edges[--j];
        const a = edges[--j];

        let i = j;
        while(i >= 2) {
            const n = edges[--i];
            const m = edges[--i];

            if((a == m && b == n) || (a == n && b == m)) {
                edges[j] = -1;
                edges[j + 1] = -1;

                edges[i] = -1;
                edges[i + 1] = -1;

                j -= 2;
                break;
            }
        }
    }

    // next remove all -1 values in the array
    let i = -1;  // actual new index
    j = 0;  // pointer to the value in the array
    while(j < in_length) {
        if(edges[j] == -1) {
            j += 1;
        } else {
            edges[++i] = edges[j++];
        }
    }

    return i + 1;
}

function remove_from_array<T>(array: StaticArray<T>, in_length: i32, remove_index: i32): i32 {
    for(let i = remove_index + 1; i < in_length; i++) {
        array[i - 1] = array[i];
    }

    return in_length - 1;
}

export function triangulate(in_points: StaticArray<Point>): StaticArray<i32> {
    const points_count: i32 = in_points.length;
    if(points_count < 3) {
        return new StaticArray<i32>(0);
    }

    // AS does not support closures, so, create array with extended points
    let ext_points = new StaticArray<PointExt>(points_count);
    for(let i = 0; i < points_count; i++) {
        const p = in_points[i];
        ext_points[i] = new PointExt(p.x(), p.y(), i);  // assign in dex to each point
    }

    // sort by using custom comparator
    ext_points.sort((a: PointExt, b: PointExt): i32 => {
        if(a.x() < b.x()) {
            return -1;
        } else if(a.x() > b.x())
        {
            return 1;
        } else {
            return 0;
        }
    });

    // extract indices from sorted array
    let indices = new StaticArray<i32>(points_count);
    for(let i = 0; i < points_count; i++) {
        indices[i] = ext_points[i].index();
    }

    const st = build_supertriangle(in_points);

    // create copy of the input points array and extend it by points from the super triangle
    const points = new StaticArray<Point>(points_count + 3);
    for(let i = 0; i < points_count; i++) {
        points[i] = in_points[i];
    }
    points[points_count] = st[0];
    points[points_count + 1] = st[1];
    points[points_count + 2] = st[2];

    // create buffers
    const buffers_max_size = 6 * (points_count + 3);
    let open_list = new StaticArray<TriangleCircle>(buffers_max_size);
    let open_list_length = 0;  // manually count the used length

    let closed_list = new StaticArray<TriangleCircle>(buffers_max_size);
    let closed_list_length = 0;

    let edges_list = new StaticArray<i32>(6 * buffers_max_size);  // here we place 6 put for each triangle
    let edges_list_length = 0;

    open_list[open_list_length++] = circumcircle(points, points_count, points_count + 1, points_count + 2);

    for(let i = 0; i < points_count; i++) {
        const c: i32 = indices[i];
        const p: Point = points[c];

        edges_list_length = 0;
        for(let j = open_list_length - 1; j >= 0; j--) {
            const t: TriangleCircle = open_list[j];
            const dx: f32 = p.x() - t.x();

            if(dx > EPSILON && dx*dx > t.radius_sq()) {
                open_list_length = remove_from_array(open_list, open_list_length, j);
                closed_list[closed_list_length++] = t;
                continue;
            }

            const dy: f32 = p.y() - t.y();
            if(dx*dx + dy*dy - t.radius_sq() > EPSILON) {
                continue;
            }

            open_list_length = remove_from_array(open_list, open_list_length, j);
            edges_list[edges_list_length++] = t.i();
            edges_list[edges_list_length++] = t.j();
            edges_list[edges_list_length++] = t.j();
            edges_list[edges_list_length++] = t.k();
            edges_list[edges_list_length++] = t.k();
            edges_list[edges_list_length++] = t.i();
        }
        
        edges_list_length = remove_duplicates(edges_list, edges_list_length);

        let k: i32 = edges_list_length;
        while(k >= 2) {
            const b = edges_list[--k];
            const a = edges_list[--k];

            open_list[open_list_length++] = circumcircle(points, a, b, c);
        }
    }

    for(let i = 0; i < open_list_length; i++) {
        closed_list[closed_list_length++] = open_list[i];
    }

    const triangles = new StaticArray<i32>(closed_list_length * 3);
    let t_index: i32 = 0;
    for(let i = 0; i < closed_list_length; i++) {
        const t = closed_list[i];
        if(t.i() < points_count && t.j() < points_count && t.k() < points_count) {
            triangles[3*t_index] = t.i();
            triangles[3*t_index + 1] = t.j();
            triangles[3*t_index + 2] = t.k();

            t_index += 1;
        }
    }
    return triangles.slice<StaticArray<i32>>(0, 3*t_index);
}

import { Point } from "./point";

@inline
function maximum(a: f32, b: f32): f32 {
    if(a > b) { return a;
    } else { return b; }
}

@inline
function minimum(a: f32, b: f32): f32 {
    if(a > b) { return b;
    } else { return a; }
}

@inline
function square(value: f32): f32 {
    return value * value;
}

class AABB {
    private m_x_min: f32;
    private m_y_min: f32;

    private m_x_max: f32;
    private m_y_max: f32;

    constructor(in_x_min: f32 = 0.0, in_y_min: f32 = 0.0, in_x_max: f32 = 0.0, in_y_max: f32 = 0.0) {
        this.m_x_min = in_x_min;
        this.m_y_min = in_y_min;

        this.m_x_max = in_x_max;
        this.m_y_max = in_y_max;
    }

    @inline
    x_min(): f32 { return this.m_x_min; }

    @inline
    y_min(): f32 { return this.m_y_min; }

    @inline
    x_max(): f32 { return this.m_x_max; }

    @inline
    y_max(): f32 { return this.m_y_max; }

    toString(): string {
        return "|" + this.m_x_min.toString() + ", " + this.m_y_min.toString() + ", " + this.m_x_max.toString() + ", " + this.m_y_max.toString() + "|";
    }
}

export class Triangle {
    private m_a: Point;
    private m_b: Point;
    private m_c: Point;

    private m_aabb: AABB;
    private m_center: Point;

    constructor(vertices: Array<Point>) {
        this.m_a = vertices[0];
        this.m_b = vertices[1];
        this.m_c = vertices[2];

        let x_min = f32.MAX_VALUE;
        let y_min = f32.MAX_VALUE;
        let x_max = f32.MIN_VALUE;
        let y_max = f32.MIN_VALUE;

        let x_accum: f32 = 0.0;
        let y_accum: f32 = 0.0;

        for (let i = 0, len = vertices.length; i < len; i++)
        {
            const v = vertices[i];
            if (v.x() < x_min) { x_min = v.x(); }
            if (v.x() > x_max) { x_max = v.x(); }
            if (v.y() < y_min) { y_min = v.y(); }
            if (v.y() > y_max) { y_max = v.y(); }

            x_accum += v.x();
            y_accum += v.y();
        }

        this.m_aabb = new AABB(x_min, y_min, x_max, y_max);
        this.m_center = new Point(x_accum / <f32>vertices.length, y_accum / <f32>vertices.length);
    }

    @inline
    a(): Point { return this.m_a; }

    @inline
    b(): Point { return this.m_b; }

    @inline
    c(): Point { return this.m_c; }

    @inline
    aabb(): AABB {
        return this.m_aabb;
    }

    @inline
    center(): Point {
        return this.m_center;
    }

    is_point_inside(point: Point): bool {
        const as_x = point.x() - this.m_a.x();
        const as_y = point.y() - this.m_a.y();

        const s_ab: bool = ((this.m_b.x() - this.m_a.x()) * as_y - (this.m_b.y() - this.m_a.y()) * as_x) > 0;

        if (((this.m_c.x() - this.m_a.x()) * as_y - (this.m_c.y() - this.m_a.y()) * as_x > 0) == s_ab)
        {
            return false;
        }
        if (((this.m_c.x() - this.m_b.x()) * (point.y() - this.m_b.y()) - (this.m_c.y() - this.m_b.y()) * (point.x() - this.m_b.x()) > 0) != s_ab)
        {
            return false;
        }
        return true;
    }

    toString(): string {
        return "<" + this.m_a.toString() + ", " + this.m_b.toString() + ", " + this.m_c.toString() + ">";
    }
}

function union_aabb(x: AABB, y: AABB): AABB {
    return new AABB(minimum(x.x_min(), y.x_min()), 
    minimum(x.y_min(), y.y_min()), 
    maximum(x.x_max(), y.x_max()), 
    maximum(x.y_max(), y.y_max()));
}

export class BVHNode {
    private m_triangle: Triangle | null = null;
    private m_left_node: BVHNode | null = null;
    private m_right_node: BVHNode | null = null;
    private m_aabb: AABB = new AABB();

    constructor(triangles: StaticArray<Triangle>) {
        const objects_count = triangles.length;

        if(objects_count == 1) {
            this.m_triangle = triangles[0];
            this.m_aabb = triangles[0].aabb();
        } else {
            let x_median: f32 = 0.0;
            let y_median: f32 = 0.0;

            let x_min = f32.MAX_VALUE;
            let x_max = f32.MIN_VALUE;
            let y_min = f32.MAX_VALUE;
            let y_max = f32.MIN_VALUE;

            for (let i = 0; i < objects_count; i++) {
                const t = triangles[i];
                const c = t.center();
                x_median += c.x();
                y_median += c.y();

                if (c.x() < x_min) { x_min = c.x(); }
                if (c.x() > x_max) { x_max = c.x(); }
                if (c.y() < y_min) { y_min = c.y(); }
                if (c.y() > y_max) { y_max = c.y(); }
            }

            const split_axis: i32 = ((x_max - x_min) > (y_max - y_min)) ? 0 : 1;
            const median: f32 = split_axis == 0 ? (x_median / <f32>objects_count) : (y_median / <f32>objects_count);

            const left = new StaticArray<Triangle>(objects_count);
            let left_length: i32 = 0;
            const right = new StaticArray<Triangle>(objects_count);
            let right_length: i32 = 0;

            for (let i = 0; i < objects_count; i++) {
                const t = triangles[i];
                const c = t.center();
                const v: f32 = split_axis == 0 ? c.x() : c.y();
                if (v < median) {
                    left[left_length++] = t;
                } else {
                    right[right_length++] = t;
                }
            }

            if (left_length == 0)
            {
                left[left_length++] = right[right_length - 1];
                right_length -= 1;
            }

            if (right_length == 0)
            {
                right[right_length++] = left[left_length - 1];
                left_length -= 1;
            }

            let left_node = new BVHNode(left.slice<StaticArray<Triangle>>(0, left_length));
            let right_node = new BVHNode(right.slice<StaticArray<Triangle>>(0, right_length));

            const left_aabb = left_node.aabb();
            const right_aabb = right_node.aabb();

            this.m_aabb = union_aabb(left_aabb, right_aabb);
            this.m_left_node = left_node;
            this.m_right_node = right_node;
        }
    }

    @inline
    aabb(): AABB {
        return this.m_aabb;
    }

    @inline
    is_inside_aabb(point: Point): bool {
        return this.m_aabb.x_min() < point.x() && this.m_aabb.y_min() < point.y() && this.m_aabb.x_max() > point.x() && this.m_aabb.y_max() > point.y();
    }

    sample(point: Point): Triangle | null {
        if(this.is_inside_aabb(point)) {
            let triangle = this.m_triangle;
            let left_node = this.m_left_node;
            let right_node = this.m_right_node;

            if(!triangle && left_node && right_node) {
                let left_sample = left_node.sample(point);
                let right_sample = right_node.sample(point);

                if(!left_sample) {
                    return right_sample;
                } else {
                    if(!right_sample) {
                        return left_sample;
                    } else {
                        const l_c = left_sample.center();
                        const l_dist = square(l_c.x() - point.x()) + square(l_c.y() - point.y());

                        const r_c = right_sample.center();
                        const r_dist = square(r_c.x() - point.x()) + square(r_c.y() - point.y());

                        if (l_dist < r_dist) {
                            return left_sample;
                        } else {
                            return right_sample;
                        }
                    }
                }
            } else {
                if(triangle && triangle.is_point_inside(point)) {
                    return triangle;
                } else {
                    return null;
                }
            }
        } else {
            return null;
        }
    }
}

import { BVHNode, Triangle } from "./bvh";
import { triangulate } from "./delaunay";
import { Point } from "./point";

function main(): void {
    let points = new StaticArray<Point>(6);
    points[0] = new Point(-0.5, 3.0);
    points[1] = new Point(8.0, -3.5);
    points[2] = new Point(-7.0, 3.0);
    points[3] = new Point(2.0, -10.0);
    points[4] = new Point(7.0, 8.0);
    points[5] = new Point(8.0, 5.5);

    let triangle_indices = triangulate(points);
    console.log(triangle_indices.toString());

    const trinagles = new StaticArray<Triangle>(triangle_indices.length / 3);
    for(let i = 0; i < trinagles.length; i++) {
        const a: i32 = triangle_indices[3*i];
        const b: i32 = triangle_indices[3*i + 1];
        const c: i32 = triangle_indices[3*i + 2];
        trinagles[i] = new Triangle([points[a], points[b], points[c]]);
    }

    const tree = new BVHNode(trinagles);
    const p = new Point(0.0, 0.0);
    const sample_result = tree.sample(p);
    if(sample_result) {
        console.log(sample_result.toString());
    } else {
        console.log("empty sample");
    }
}

main();

import { BVHNode, Triangle } from "./bvh";
import { triangulate } from "./delaunay";
import { Point } from "./point";

export function build_triangulation(coordinates: Float32Array): Int32Array {
    const points_count = coordinates.length / 2;
    const points = new StaticArray<Point>(points_count);
    for(let i = 0; i < points_count; i++) {
        points[i] = new Point(coordinates[2*i], coordinates[2*i + 1]);
    }

    const triangle_indices = triangulate(points);

    // convert to Int32Array
    const to_return = new Int32Array(triangle_indices.length);
    for(let i = 0, len = triangle_indices.length; i < len; i++) {
        to_return[i] = triangle_indices[i];
    }

    return to_return;
}

export function build_bvh(coordinates: Float32Array): BVHNode {
    const points_count = coordinates.length / 2;
    const points = new StaticArray<Point>(points_count);
    for(let i = 0; i < points_count; i++) {
        points[i] = new Point(coordinates[2*i], coordinates[2*i + 1]);
    }

    const triangle_indices = triangulate(points);

    const trinagles = new StaticArray<Triangle>(triangle_indices.length / 3);
    for(let i = 0; i < trinagles.length; i++) {
        const a: i32 = triangle_indices[3*i];
        const b: i32 = triangle_indices[3*i + 1];
        const c: i32 = triangle_indices[3*i + 2];
        trinagles[i] = new Triangle([points[a], points[b], points[c]]);
    }

    return new BVHNode(trinagles);
}

export function build_bvh_triangulation(coordinates: Float32Array, triangles: Int32Array): BVHNode {
    const trinagles_array = new StaticArray<Triangle>(triangles.length / 3);
    for(let i = 0, len = trinagles_array.length; i < len; i++) {
        const a: i32 = triangles[3*i];
        const b: i32 = triangles[3*i + 1];
        const c: i32 = triangles[3*i + 2];
        trinagles_array[i] = new Triangle([
            new Point(coordinates[2*a], coordinates[2*a + 1]),
            new Point(coordinates[2*b], coordinates[2*b + 1]),
            new Point(coordinates[2*c], coordinates[2*c + 1])]);
    }

    return new BVHNode(trinagles_array);
}

export function sample_bvh(bvh: BVHNode, x: f32, y: f32): Float32Array {
    const sample_result = bvh.sample(new Point(x, y));
    if(sample_result) {
        const to_return = new Float32Array(6);
        const a = sample_result.a();
        const b = sample_result.b();
        const c = sample_result.c();
        to_return[0] = a.x(); to_return[1] = a.y();
        to_return[2] = b.x(); to_return[3] = b.y();
        to_return[4] = c.x(); to_return[5] = c.y();

        return to_return;
    } else {
        return new Float32Array(0);
    }
}

#[derive(Debug, Clone)]
pub struct Point {
    pub x: f32,
    pub y: f32
}

impl Point {
    pub fn new(in_x: f32, in_y: f32) -> Point {
        Point{ x: in_x, y: in_y}
    }
}

use crate::point::Point;

const EPSILON: f32 = 0.00001;

#[derive(Debug, Clone)]
struct TriangleCircle {
    i: usize,
    j: usize,
    k: usize,

    x: f32,
    y: f32,

    radius_sq: f32
}

fn maximum(a: f32, b: f32) -> f32 {
    if a > b { a } else { b }
}

fn absolute(a: f32, b: f32) -> f32 {
    if a > b { a - b } else { b - a }
}

fn build_supertriangle(points: &Vec<Point>) -> Vec<Point> {
    let mut x_min: f32 = f32::MAX;
    let mut y_min: f32 = f32::MAX;
    let mut x_max: f32 = f32::MIN;
    let mut y_max: f32 = f32::MIN;

    let points_count = points.len();

    for i in 0..points_count {
        let p = &points[i];
        if p.x < x_min { x_min = p.x; }
        if p.x > x_max { x_max = p.x; }
        if p.y < y_min { y_min = p.y; }
        if p.y > y_max { y_max = p.y; }
    }

    let dx = x_max - x_min;
    let dy = y_max - y_min;

    let d_max = maximum(dx, dy);
    let x_mid = x_min + dx * 0.5;
    let y_mid = y_min + dy * 0.5;

    vec![Point::new(x_mid - 20.0 * d_max, y_mid - d_max),
         Point::new(x_mid, y_mid + 20.0 * d_max),
         Point::new(x_mid + 20.0 * d_max, y_mid - d_max)]
}

fn circumcircle(points: &Vec<Point>, i: usize, j: usize, k: usize) -> TriangleCircle {
    let point_i = &points[i];
    let point_j = &points[j];
    let point_k = &points[k];

    let x_1: f32 = point_i.x;
    let y_1: f32 = point_i.y;

    let x_2: f32 = point_j.x;
    let y_2: f32 = point_j.y;    

    let x_3: f32 = point_k.x;
    let y_3: f32 = point_k.y;

    let y1_y2: f32 = absolute(y_1, y_2);
    let y2_y3: f32 = absolute(y_2, y_3);

    let center_x: f32;
    let center_y: f32;
    let m_1: f32;
    let m_2: f32;
    let mx_1: f32;
    let mx_2: f32;
    let my_1: f32;
    let my_2: f32;

    if y1_y2 < EPSILON {
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_2 = (x_2 + x_3) / 2.0;
        my_2 = (y_2 + y_3) / 2.0;
        center_x = (x_2 + x_1) / 2.0;
        center_y = m_2 * (center_x - mx_2) + my_2;
    } else if y2_y3 < EPSILON {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        mx_1 = (x_1 + x_2) / 2.0;
        my_1 = (y_1 + y_2) / 2.0;
        center_x = (x_3 + x_2) / 2.0;
        center_y = m_1 * (center_x - mx_1) + my_1;
    } else {
        m_1 = -(x_2 - x_1) / (y_2 - y_1);
        m_2 = -(x_3 - x_2) / (y_3 - y_2);
        mx_1 = (x_1 + x_2) / 2.0;
        mx_2 = (x_2 + x_3) / 2.0;
        my_1 = (y_1 + y_2) / 2.0;
        my_2 = (y_2 + y_3) / 2.0;
        center_x = (m_1 * mx_1 - m_2 * mx_2 + my_2 - my_1) / (m_1 - m_2);
        center_y = if y1_y2 > y2_y3 { m_1 * (center_x - mx_1) + my_1 } else { m_2 * (center_x - mx_2) + my_2 };
    }

    let dx = x_2 - center_x;
    let dy = y_2 - center_y;

    TriangleCircle { i: i, j: j, k: k, x: center_x, y: center_y, radius_sq: dx*dx + dy*dy }
}

fn remove_duplicates(edges: &mut Vec<usize>) {
    let mut j: usize = edges.len();
    while j >= 2 {
        j = j - 1;
        let b = edges[j];

        j = j - 1;
        let a = edges[j];

        let mut i: usize = j;
        while i >= 2 {
            i = i - 1;
            let n = edges[i];
            i = i - 1;
            let m = edges[i];

            if (a == m && b == n) || (a == n && b == m) {
                edges.remove(j + 1); edges.remove(j);
                edges.remove(i + 1); edges.remove(i);
                j = j - 2;
                break;
            }
        }
    }
}

pub fn triangulate(in_points: &Vec<Point>) -> Vec<usize> {
    let points_count = in_points.len();
    if points_count < 3 {
        return Vec::new();
    }

    let mut indices: Vec<usize> = Vec::with_capacity(points_count);
    for i in 0..points_count {
        indices.push(i);
    }

    indices.sort_by(|a, b| { in_points[*a].x.partial_cmp(&in_points[*b].x).unwrap() });
    
    // create copy of the in_points array
    let mut points: Vec<Point> = in_points.to_vec();
    let mut st = build_supertriangle(&points);
    points.append(&mut st);

    let mut open_list: Vec<TriangleCircle> = Vec::new();
    let mut closed_list: Vec<TriangleCircle> = Vec::new();
    let mut edges_list: Vec<usize> = Vec::new();

    open_list.push(circumcircle(&points, points_count, points_count + 1, points_count + 2));

    for i in 0..points_count {
        let c = indices[i];
        let p = &points[c];
        edges_list.clear();

        for j in (0..open_list.len()).rev() {
            let t = &open_list[j];

            let dx = p.x - t.x;
            if dx > 0.0 && dx*dx > t.radius_sq {
                closed_list.push(t.clone());
                open_list.remove(j);
                continue;
            }

            let dy = p.y - t.y;
            if dx*dx + dy*dy - t.radius_sq > EPSILON {
                continue;
            }
            
            edges_list.push(t.i); edges_list.push(t.j);
            edges_list.push(t.j); edges_list.push(t.k);
            edges_list.push(t.k); edges_list.push(t.i);

            open_list.remove(j);
        }

        remove_duplicates(&mut edges_list);

        let mut k: usize = edges_list.len();
        while k >= 2 {
            k = k - 1;
            let b = edges_list[k];
            k = k - 1;
            let a = edges_list[k];
            open_list.push(circumcircle(&points, a, b, c));
        }
    }

    for i in 0..open_list.len() {
        closed_list.push(open_list[i].clone());
    }

    open_list.clear();
    edges_list.clear();

    let mut triangles: Vec<usize> = Vec::new();
    for i in 0..closed_list.len() {
        let t = &closed_list[i];
        if t.i < points_count && t.j < points_count && t.k < points_count {
            triangles.push(t.i);
            triangles.push(t.j);
            triangles.push(t.k);
        }
    }

    closed_list.clear();

    triangles
}

use crate::point::Point;

fn minimum(a: f32, b: f32) -> f32 {
    if a < b { return a; } else { return b; }
}

fn maximum(a: f32, b: f32) -> f32 {
    if a < b { return b; } else { return a; }
}

#[derive(Debug, Clone)]
struct AABB {
    x_min: f32,
    y_min: f32,
    x_max: f32,
    y_max: f32,
}

#[derive(Debug, Clone)]
pub struct Triangle {
    pub a: Point,
    pub b: Point,
    pub c: Point,

    aabb: AABB,
    center: Point
}

impl Triangle {
    pub fn new(mut vertices: Vec<Point>) -> Triangle {
        let mut x_min: f32 = f32::MAX;
        let mut y_min: f32 = f32::MAX;
        let mut x_max: f32 = f32::MIN;
        let mut y_max: f32 = f32::MIN;

        let mut x_accum: f32 = 0.0;
        let mut y_accum: f32 = 0.0;

        for i in 0..vertices.len() {
            let v = &vertices[i];
            if v.x < x_min { x_min = v.x; }
            if v.x > x_max { x_max = v.x; }
            if v.y < y_min { y_min = v.y; }
            if v.y > y_max { y_max = v.y; }

            x_accum += v.x;
            y_accum += v.y;
        }

        let c = vertices.pop().unwrap();
        let b = vertices.pop().unwrap();
        let a = vertices.pop().unwrap();

        Triangle {
            a: a,
            b: b,
            c: c,
            aabb: AABB { x_min, y_min, x_max, y_max },
            center: Point::new(x_accum / 3.0, y_accum / 3.0)}
    }

    fn get_aabb(&self) -> AABB {
        self.aabb.clone()
    }

    fn get_center(&self) -> Point {
        self.center.clone()
    }

    fn is_point_inside(&self, point: &Point) -> bool {
        let as_x = point.x - self.a.x;
        let as_y = point.y - self.a.y;

        let s_ab = ((self.b.x - self.a.x) * as_y - (self.b.y - self.a.y) * as_x) > 0.0;

        if ((self.c.x - self.a.x) * as_y - (self.c.y - self.a.y) * as_x > 0.0) == s_ab {
            return false;
        }
        if ((self.c.x - self.b.x) * (point.y - self.b.y) - (self.c.y - self.b.y) * (point.x - self.b.x) > 0.0) != s_ab {
            return false;
        }
        return true;
    }
}

fn union_aabb(x: AABB, y: AABB) -> AABB {
    AABB {
        x_min: minimum(x.x_min, y.x_min),
        y_min: minimum(x.y_min, y.y_min),
        x_max: maximum(x.x_max, y.x_max),
        y_max: maximum(x.y_max, y.y_max),
    }
}

#[derive(Debug)]
pub struct BVHNode {
    triangle: Option<Triangle>,
    left_node: Option<Box<BVHNode>>,
    right_node: Option<Box<BVHNode>>,
    aabb: AABB
}

impl BVHNode {
    pub fn new(mut triangles: Vec<Triangle>) -> BVHNode {
        let objects_count = triangles.len();
        if objects_count == 1 {
            let t: Triangle = triangles.pop().unwrap();
            let aabb = t.get_aabb();
            BVHNode {
                triangle: Some(t),
                left_node: None,
                right_node: None,
                aabb: aabb
            }
        } else {
            let mut x_median: f32 = 0.0;
            let mut y_median: f32 = 0.0;

            let mut x_min = f32::MAX;
            let mut x_max = f32::MIN;
            let mut y_min = f32::MAX;
            let mut y_max = f32::MIN;

            for i in 0..objects_count {
                let t = &triangles[i];
                let c = t.get_center();
                x_median += c.x;
                y_median += c.y;

                if c.x < x_min { x_min = c.x; }
                if c.x > x_max { x_max = c.x; }
                if c.y < y_min { y_min = c.y; }
                if c.y > y_max { y_max = c.y; }
            }

            let split_axis = if (x_max - x_min) > (y_max - y_min) { 0 } else { 1 };
            let median: f32 = if split_axis == 0 { x_median / objects_count as f32 } else { y_median / objects_count as f32 };

            let mut left: Vec<Triangle> = Vec::new();
            let mut right: Vec<Triangle> = Vec::new();

            for _i in 0..objects_count {
                let t = triangles.pop().unwrap();
                let c = t.get_center();
                let v = if split_axis == 0 { c.x } else { c.y };
                if v < median {
                    left.push(t);
                } else {
                    right.push(t);
                }
            }

            if left.len() == 0 {
                left.push(right.pop().unwrap());
            }

            if right.len() == 0 {
                right.push(left.pop().unwrap());
            }

            let left_node = BVHNode::new(left);
            let right_node = BVHNode::new(right);

            let left_aabb = left_node.get_aabb();
            let right_aabb = right_node.get_aabb();

            BVHNode {
                left_node: Some(Box::new(left_node)),
                right_node: Some(Box::new(right_node)),
                triangle: None,
                aabb: union_aabb(left_aabb, right_aabb)
            }
        }
    }

    fn get_aabb(&self) -> AABB {
        self.aabb.clone()
    }

    fn is_inside_aabb(&self, point: &Point) -> bool {
        self.aabb.x_min < point.x && self.aabb.y_min < point.y && self.aabb.x_max > point.x && self.aabb.y_max > point.y
    }

    pub fn sample(&self, point: &Point) -> Option<Triangle> {
        if self.is_inside_aabb(point) {
           match &self.triangle {
                Some(t) => {
                    if t.is_point_inside(point) {
                        Some(t.clone())
                    } else {
                        None
                    }
                },
                None => {
                    match &self.left_node {
                        Some(node) => {
                            // left node exists
                            match &self.right_node {
                                Some(other) => {
                                    // left and right are exists
                                    let left_sample = node.sample(point);
                                    let right_sample = other.sample(point);

                                    match left_sample {
                                        Some(t) => {
                                            // left sample valid
                                            // check right sample
                                            match right_sample {
                                                Some(u) => {
                                                    // both right and left are valid
                                                    // it should not happens, but we select the closest triangle
                                                    let l_c = t.get_center();
                                                    let r_c = u.get_center();

                                                    let l_dist = (l_c.x - point.x)*(l_c.x - point.x) + (l_c.y - point.y)*(l_c.y - point.y);
                                                    let r_dist = (r_c.x - point.x)*(r_c.x - point.x) + (r_c.y - point.y)*(r_c.y - point.y);

                                                    if l_dist < r_dist {
                                                        Some(t)
                                                    } else {
                                                        Some(u)
                                                    }
                                                },
                                                None => {
                                                    // left valid, right empty
                                                    Some(t)
                                                },
                                            }
                                        },
                                        None => {
                                            // left sample empty
                                            right_sample
                                        },
                                    }
                                },
                                None => {
                                    // only left
                                    node.sample(point)
                                },
                            }
                        },
                        None => {
                            // left node is empty
                            match &self.right_node {
                                Some(node) => {
                                    // only right
                                    node.sample(point)
                                },
                                None => {
                                    // both nodes are empty
                                    None
                                }
                            }
                        }
                    }
                }
            }
        } else {
            None
        }
    }
}

mod point;
mod delaunay;
mod bvh;

use crate::point::Point;
use crate::delaunay::triangulate;
use crate::bvh::{BVHNode, Triangle};

use wasm_bindgen::prelude::*;
use js_sys::Array;

#[wasm_bindgen]
pub fn build_triangulation(coordinates: &[f32]) -> Array {
    let points_count = coordinates.len() / 2;
    let mut points: Vec<Point> = Vec::with_capacity(points_count);
    for i in 0..points_count {
        points.push(Point::new(coordinates[2*i], coordinates[2*i + 1]));
    }

    let triangle_indices = triangulate(&points);

    triangle_indices.into_iter().map(JsValue::from).collect()
}

#[derive(Debug)]
#[wasm_bindgen(js_name = BVHNode)]
pub struct BVHNodeWrapper {
    bvh: BVHNode
}

#[wasm_bindgen(js_class = BVHNode)]
impl BVHNodeWrapper {
    #[wasm_bindgen(constructor)]
    pub fn new(coordinates: &[f32]) -> BVHNodeWrapper {
        let points_count = coordinates.len() / 2;
        let mut points: Vec<Point> = Vec::with_capacity(points_count);
        for i in 0..points_count {
            points.push(Point::new(coordinates[2*i], coordinates[2*i + 1]));
        }

        let triangle_indices = triangulate(&points);

        let mut triangles: Vec<Triangle> = Vec::with_capacity(triangle_indices.len() / 3);
        for i in 0..(triangle_indices.len() / 3) {
            let a = triangle_indices[3*i];
            let b = triangle_indices[3*i + 1];
            let c = triangle_indices[3*i + 2];
            triangles.push(Triangle::new(vec!(points[a].clone(), points[b].clone(), points[c].clone())));
        }

        BVHNodeWrapper {
            bvh: BVHNode::new(triangles)
        }
    }

    pub fn new_trianglulation(coordinates: &[f32], triangles: &[i32]) -> BVHNodeWrapper {
        let triangles_count = triangles.len() / 3;
        let mut triangles_array: Vec<Triangle> = Vec::with_capacity(triangles_count);
        for i in 0..triangles_count {
            let a = triangles[3*i] as usize;
            let b = triangles[3*i + 1] as usize;
            let c = triangles[3*i + 2] as usize;

            triangles_array.push(Triangle::new(vec!(Point::new(coordinates[2*a], coordinates[2*a + 1]), Point::new(coordinates[2*b], coordinates[2*b + 1]), Point::new(coordinates[2*c], coordinates[2*c + 1]))));
        }

        BVHNodeWrapper {
            bvh: BVHNode::new(triangles_array)
        }
    }

    pub fn sample(&self, x: f32, y: f32) -> Array {
        let point = Point::new(x, y);
        let s = self.bvh.sample(&point);

        match s {
            Some(t) => {
                let to_return = Array::new_with_length(6);
                to_return.set(0, JsValue::from(t.a.x)); to_return.set(1, JsValue::from(t.a.y));
                to_return.set(2, JsValue::from(t.b.x)); to_return.set(3, JsValue::from(t.b.y));
                to_return.set(4, JsValue::from(t.c.x)); to_return.set(5, JsValue::from(t.c.y));
                to_return
            },
            None => {
                Array::new()
            },
        }
    }
}

<!DOCTYPE html>
<html>
    <head>
        <meta charset="UTF-8">
        <title>Delaunay Demo</title>
        <style>
            body, html {
                height: 100%;
                width: 100%;
                padding: 0px;
                margin: 0px;
                background-color: rgb(40, 40, 40);
            }
            
            div#central_div {
                position: relative;
                width: 100%;
                height: 100%;
            }
            
            .points_button {
                width: 80px;
                height: 35px;
                border: none;
                border-radius: 10px;
                font-size: 14px;
                font-weight: bold;
                position: absolute;
                top: 10px;
            }
            
            .points_button#btn_0 {
                left: 10px;
            }
            
            .points_button#btn_1 {
                left: 100px;
            }
            
            .points_button#btn_2 {
                left: 190px;
            }
            
            .points_button#btn_3 {
                left: 280px;
            }
        </style>
    </head>
    <body>
        <div id="central_div">
            <button class="points_button" id="btn_0" onclick="generate(10)">10</button>
            <button class="points_button" id="btn_1" onclick="generate(100)">100</button>
            <button class="points_button" id="btn_2" onclick="generate(1000)">1000</button>
            <button class="points_button" id="btn_3"onclick="generate(12000)">12000</button>
            <canvas id="canvas"></canvas>
        </div>
        <script type="module">
            import init, { BVHNode, build_triangulation } from "./delaunay.js";

            async function wasm_run() {
                await init();

                const canvas = document.getElementById("canvas");
                const ctx = canvas.getContext("2d");

                var canvas_width = document.body.clientWidth;
                var canvas_height = document.body.clientHeight - 3;

                var points = [];
                var triangles = [];
                var bvh = undefined;
                var cursor_position = [0, 0];
                var sample_coordinates = [];

                function draw_canvas() {
                    ctx.clearRect(0, 0, canvas_width, canvas_height);
                    ctx.save();
                    ctx.fillStyle = "rgb(40, 40, 40)";
                    ctx.fillRect(0, 0, canvas_width, canvas_height);
                    ctx.restore();

                    if(sample_coordinates.length == 6) {
                        ctx.save();
                        ctx.beginPath();
                        ctx.fillStyle = "rgb(143, 186, 173)";
                        ctx.strokeStyle = "rgb(143, 186, 173)";
                        ctx.moveTo(sample_coordinates[0], sample_coordinates[1]);
                        ctx.lineTo(sample_coordinates[2], sample_coordinates[3]);
                        ctx.lineTo(sample_coordinates[4], sample_coordinates[5]);
                        ctx.lineTo(sample_coordinates[0], sample_coordinates[1]);
                        ctx.fill();
                        ctx.restore();
                    }

                    ctx.save();
                    ctx.fillStyle = "rgb(129, 129, 129)";
                    const points_count = points.length / 2;
                    for(let i = 0; i < points_count; i++) {
                        const x = points[2*i];
                        const y = points[2*i + 1];
                        ctx.beginPath();
                        ctx.arc(x, y, 2, 0, 2 * Math.PI);
                        ctx.fill();
                    }
                    ctx.restore();

                    ctx.save();
                    ctx.strokeStyle = "rgb(129, 129, 129)";
                    const triangles_count = triangles.length / 3;
                    for(let i = 0; i < triangles_count; i++) {
                        ctx.beginPath();
                        const a = triangles[3*i]; const a_x = points[2*a]; const a_y = points[2*a + 1];
                        const b = triangles[3*i + 1]; const b_x = points[2*b]; const b_y = points[2*b + 1];
                        const c = triangles[3*i + 2]; const c_x = points[2*c]; const c_y = points[2*c + 1];
                        ctx.moveTo(a_x, a_y);
                        ctx.lineTo(b_x, b_y);
                        ctx.lineTo(c_x, c_y);
                        ctx.lineTo(a_x, a_y);
                        ctx.stroke();
                    }
                    ctx.restore();
                }

                function update_canvas() {
                    canvas_width = document.body.clientWidth;
                    canvas_height = document.body.clientHeight - 3;

                    canvas.width = canvas_width;
                    canvas.height = canvas_height;

                    draw_canvas();
                }

                onresize = (event) => {
                    update_canvas();
                };
                onmousemove = (event) => {
                    var rect = canvas.getBoundingClientRect();
                    cursor_position[0] = event.clientX - rect.left;
                    cursor_position[1] = event.clientY - rect.top;

                    if(bvh) {
                        sample_coordinates = bvh.sample(cursor_position[0], cursor_position[1]);
                    }

                    draw_canvas();
                };

                update_canvas();

                generate(10);

                function generate(n) {
                    const x_min = 10;
                    const y_min = 10;
                    const x_max = canvas_width - 10;
                    const y_max = canvas_height - 10;

                    points = []

                    for(let i = 0; i < n; i++) {
                        const x = Math.random() * (x_max - x_min) + x_min;
                        const y = Math.random() * (y_max - y_min) + y_min;

                        points.push(x, y);
                    }

                    triangles = build_triangulation(points);
                    bvh = new BVHNode(points, triangles);

                    draw_canvas();
                }

                window.generate = generate;
            };

            wasm_run();
        </script>
    </body>
</html>

from wasmer import Store, Module, Instance, ImportObject, Function

class Loader:
    def __init__(self, wasm_path):
        self._heap = [None] * 128
        self._heap.extend([None, None, True, False])

        self._heap_next = len(self._heap)

        self._wasm_vector_len = 0

        store = Store()
        module = Module(store, open(wasm_path, "rb").read())

        imports = ImportObject()
        imports.register("wbg", {
            "__wbindgen_number_new": Function(store, self._import_number_new),
            "__wbindgen_object_drop_ref": Function(store, self._import_drop_ref),
            "__wbg_new_b525de17f44a8943": Function(store, self._import_new),
            "__wbg_newwithlength_0da6f12fbc1ab6eb": Function(store, self._import_new_with_length),
            "__wbg_set_17224bc548dd1d7b": Function(store, self._import_set),
            "__wbg_push_49c286f04dd3bf59": Function(store, self._import_push),
            "__wbindgen_throw": Function(store, self._import_throw),
        })

        self._wasm = Instance(module, imports)
        wasm_exports = self._wasm.exports
        wasm_memory = wasm_exports.memory
        self._float32_memory = wasm_memory.float32_view()
        self._uint32_memory = wasm_memory.uint32_view()

        for e in wasm_exports:
            export_name = e[0]
            export_function = e[1]
            if export_name == "build_triangulation":
                self._build_triangulation = export_function
            elif export_name == "bvhnode_new":
                self._new_bvh = export_function
            elif export_name == "bvhnode_new_trianglulation":
                self._new_bvh_triangulation = export_function
            elif export_name == "bvhnode_sample":
                self._bvh_sample = export_function
            if export_name == "__wbindgen_malloc":
                self._malloc = export_function
            elif export_name == "__wbg_bvhnode_free":
                self._bvh_free = export_function

    def _import_number_new(self, arg: "f64") -> int:
        return self._add_heap_object(arg)

    def _import_drop_ref(self, arg: int):
        self._take_object(arg)

    def _import_new(self) -> int:
        return self._add_heap_object([])

    def _import_new_with_length(self, arg: int) -> int:
        return self._add_heap_object([None] * arg)

    def _import_set(self, arg0: int, arg1: int, arg2: int):
        self._get_object(arg0)[arg1] = self._take_object(arg2)

    def _import_push(self, arg0: int, arg1: int) -> int:
        self._get_object(arg0).append(self._get_object(arg1))
        return len(self._heap)

    def _import_throw(self, arg0: int, arg1: int):
        return "something wrong"

    def _add_heap_object(self, obj: int) -> int:
        if self._heap_next == len(self._heap):
            self._heap.append(len(self._heap) + 1)

        idx = self._heap_next
        self._heap_next = self._heap[idx]
        self._heap[idx] = obj

        return idx

    def _get_object(self, idx: int) -> int:
        return self._heap[idx]

    def _drop_object(self, idx: int):
        if idx < 132:
            return None

        self._heap[idx] = self._heap_next
        self._heap_next = idx

    def _take_object(self, idx: int):
        ret: int = self._get_object(idx)
        self._drop_object(idx)
        return ret

    def _pass_array_f32_to_wasm(self, arg) -> int:
        ptr = self._malloc(len(arg) * 4)
        for i in range(len(arg)):
            self._float32_memory[ptr // 4 + i] = arg[i]
        self._wasm_vector_len = len(arg)
        return ptr

    def _pass_array_32_to_wasm(self, arg) -> int:
        ptr: int = self._malloc(len(arg) * 4)
        for i in range(len(arg)):
            self._uint32_memory[ptr // 4 + i] = arg[i]
        self._wasm_vector_len = len(arg)
        return ptr

    def build_triangulation(self, coordinates: list[float]) -> list[int]:
        ptr = self._pass_array_f32_to_wasm(coordinates)
        length = self._wasm_vector_len
        ret = self._build_triangulation(ptr, length)
        array = self._take_object(ret)
        return [int(v) for v in array]

class BVHNode:
    def __init__(self, loader, coordinates: list[float], triangles=None):
        self._loader = loader
        if triangles is None:
            ptr = loader._pass_array_f32_to_wasm(coordinates)
            length = loader._wasm_vector_len
            self._ptr = loader._new_bvh(ptr, length)
        else:
            ptr_0 = loader._pass_array_f32_to_wasm(coordinates)
            length_0 = loader._wasm_vector_len
            ptr_1 = loader._pass_array_32_to_wasm(triangles)
            length_1 = loader._wasm_vector_len
            self._ptr = loader._new_bvh_triangulation(ptr_0, length_0, ptr_1, length_1)

    def __del__(self):
        self._loader._bvh_free(self._ptr)

    def sample(self, x:float, y: float) -> list[float]:
        return self._loader._take_object(self._loader._bvh_sample(self._ptr, x, y))