В этой небольшой заметке предлагается рассмотреть несколько способов компиляции модуля для WebAssembly, используя три разных подхода. Мы реализуем решение одной и той же задачки на трёх языках и скомпилируем полученный код в модуль WebAssembly. Будем использовать:

  • Emscripten для компиляции кода, написанного на c++

  • AssemblyScript для компиляции кода, написанного на, собственно, AssemblyScript

  • wasm-pack для компиляции кода, написанного на Rust

План такой:

  • Во введении мы обсудим постановку задачи и немножко поговорим о технологии WebAssembly

  • В программной части мы реализуем функциональность модуля на трёх языках: c++, AssemblyScript и Rust. Поговорим о том, какие при этом возникают сложности и как их можно обойти

  • Подведём небольшой итог. Станет видно, какая технология хорошая, а какая не очень

В конце планируется два бонуса. Первый бонус - это пример простого web-приложения, использующего один из скомпилированных нами модулей. Второй бонус - демонстрация того, как этот модуль можно использовать в программе на Python.

Введение

Предполагается, что читатель знаком хотя бы понаслышке с технологией WebAssembly и выполнил один из HelloWorld туториалов для одной из перечисленных выше технологий. Впрочем, можно было и не выполнять. Обычно эти HelloWorld-ы описывают то, как сделать модуль для вычисления чисел Фибоначчи, ну или в крайнем случае - для нахождения числа делителей данного числа. Текущая заметка предлагает сделать маленький (ну вот буквально всего ничего) шажок от HelloWorld-а и решить чуть более содержательную задачу по сравнению с вычислением чисел Фибоначчи.

Итак, задача. Предлагается написать модуль, который выполняет следующее. На вход подаётся множество точек на плоскости. Модуль строит триангуляцию Делоне с вершинами в этих точках, упаковывает построенные треугольники в BVH и реализует метод поиска треугольника, содержащего любую указанную точку.

Будем считать, что мы знаем, что такое триангуляция Делоне. Если кто забыл (или не знал) - можно вспомнить на Википедии. Есть много алгоритмов её построения. Предлагается использовать один из них вот отсюда. Может быть этот алгоритм не сильно быстрый, но по крайней мере простой и его можно реализовать довольно короткой программой. В любом случае сам алгоритм не имеет значения. Заметка ведь не об этом.

Что такое BVH (Bounding Volume Hierarhy) тоже, наверное, все знают. Ещё в школе, поди, в восьмом классе писали разные реализации. Ну а если кому не повезло или он забыл своё счастливое детство - можно тоже ознакомиться с определением на Википедии. И снова, по большому счёту, что это такое не особо важно.

И ещё пару слов о том, что такое WebAssembly. Есть определение на всё той же Википедии, есть в документации Mozill-ы. Можно по разному расставлять акценты в вопросах о том, как относиться к технологии WebAssembly и для чего её можно использовать. На мой взгляд одно из преимущество модулей WebAssembly - это их универсальность в смысле независимости от платформы. Эти модули можно (и нужно) запускать не только в браузере и NodeJS. Есть множество стационарных итерпретаторов и компиляторов в настоящие машинные коды. Вот тут большой список разных таких проектов. Помимо этого модули WebAssembly можно встраивать в программы скриптовых (читай медленных) языках программирования, и тем самым какие-то вычислительно тяжёлые куски программы выполнять гораздо быстрее. Технология Wasmer позволяет, в частности, использовать модули WebAssembly внутри программ как на Python, так и на некоторых других языках программирования.

Писать код непосредственно на WebAssembly можно. Но нужно это примерно столько же раз, сколько писать на обычном ассемблере, и даже реже. Поэтому обычно пишут на каком-нибудь человеческом языке, а потом полученную программу компилируют в WebAssembly. И вот тут вот и кроется одна из проблем. Дело в том, что довольно много языков программирования компилируются в WebAssembly. Вот тут есть список. Много там всякого разного написано. Но мало модуль скомпилировать, его ведь потом надо использовать из какого-то другого места. Например, вызвать функцию, реализованную в модуле, и получить результат. Но вот беда, язык WebAssembly понимает только числовые типы данных: целые и с плавающей точкой разной битности. И с этим нет проблем, если мы хотим вызвать функцию, которая принимает только числовые аргументы и выдаёт в качестве результата тоже число. А если мы, например, хотим посчитать сумму элементов массива. Это ведь надо как-то передать массив в модуль. Для этого его надо разместить в памяти модуля, а программа уже должна сама уметь его из этой памяти взять. Вот и проблема: откуда я как программист знаю в каком виде представлена память модуля после компиляции, как оттуда что-то брать и, прости господи, что-то туда класть? Говорю же - проблемы. Да, они решаются, но процедура использования откомпилированного модуля становится весьма нетривиальной.

Однако не всё так плохо. Есть несколько технологий, которые не только компилируют модуль WebAssembly, но и создают обвес (так называемый gluing code), который и делает все вот эти неудобные действия, заставляющие задумываться и отчаиваться. Эти технологии перечислены в самом начале. Может быть есть ещё какие-то другие, но давайте для определённости ограничимся этими тремя. Они, всё-таки, наиболее популярны.

Программная часть

Итак, приступим к программированию. Ещё раз повторю, что наша цель состоит в том, чтобы запрограммировать упаковку в BVH треугольников триангуляции Делоне. Начнём с реализации на c++, потом сделаем то же самое на AssemblyScript, а закончим на Rust.

Вот ссылка на репозиторий с кодом.

Emscripten и c++

Предлагается разыграть такой сценарий. Представим, что мы как будто бы пишем обычную программу на c++, не имея ввиду какую-то последующую компиляцию в WebAssembly. Это удобно, так как можно использовать привычные средства разработки, отладки и чего там ещё нужно. Я на c++ никогда особо ничего не писал, поэтому буду делать всё в VisualStudio, создав проект для обычного консольного приложения. Потом, когда уже всё напишем, тогда и будем думать о том, как всё это дело скомпилировать в WebAssembly.

Программа на c++

Итак, создаём пять файлов: point.h delaunay.h delaunay.cpp bvh.h и bvh.cpp

Файл point.h содержит структуру для описания точки на плоскости.

point.h
#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;
};

Файл delaunay.h содержит только описание сигнатуры функции для построения триангуляции Делоне. На вход она принимает массив точек, а возвращает массив индексов этих точек, описывающих треугольники. Предполагается, что первые три числа итогового массива - это индексы вершин первого треугольника, следующие три - второго, и так далее.

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

Файл delaunay.cpp содержит непосредственно реализацию алгоритма построения триангуляции Делоне. Там используется несколько вспомогательных функций и структур, но в целом ничего сложного. Какого бы то ни было, разбора использующегося алгоритма здесь давать не будем, так как это уведёт нас далеко в сторону.

delaunay.cpp
#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;
}

Файл bvh.h содержит описание классов, необходимых для построения самого дерева BVH. Объекты, которые складываются в bvh - это отдельная структура Triangle, предназначенная для описания треугольников на плоскости.

bvh.h
#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;
};

Ну и последний файл bvh.cpp содержит реализации методов для всех нужных классов.

bvh.cpp
#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;
	}
}

Теперь можем воспользоваться написанной программой. Создадим вот такой файл main.cpp.

main.cpp
#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;
    }
}

Откомпилируем его и запустим. Он выведет на экран правильный ответ

[2, 0, 3, 0, 2, 4, 3, 0, 1, 1, 0, 5, 0, 4, 5]
<(-7.000000, 3.000000), (-0.500000, 3.000000), (2.000000, -10.000000)>

Первая строчка - индексы треугольников триангуляции. Вторая строчка - вершины треугольника, содержащего начало координат.

Всё, программа готова, теперь можно начинать её экспорт в модуль WebAssembly. Но сначала надо определиться, что мы, собственно, хотим экспортировать.

  • Функцию для построения триангуляции

  • Класс BVHNode

  • Метод sample у класса BVHNode, который находит треугольник, содержащий данную точку

Использование Emscripten

Emscripten - это набор инструментов, позволяющих откомпилировать код на c и c++ в модуль WebAssembly. Как устанавливать Emscripten сказано в официальной документации. По сути дела надо просто клонировать репозиторий на гитхабе. Сам Emscripten содержит несколько инструментов для компиляции c++ кода в WebAssembly. Предлагается использовать embind

И опять же, для того, чтобы всё было максимально просто, постараемся делать как можно больше вещей вручную. Это позволит избежать ситуаций, когда "оно там само". На самом деле я просто по другому не умею и прячу это за демонстрацией заботы о читателе.

Создаём отдельный файл, который будет содержать описание того, что мы хотим экспортировать, а также некоторые дополнительные действия по преобразованию данных между нашей программой на c++ и теми данными, которые приходят извне, а также туда уходят. Назовём его delaunay_api.cpp

Вот его код

delaunay_api.cpp

#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);
}

Тут нужны комментарии по поводу того, что делается. Самое главное - это последний участок кода EMSCRIPTEN_BINDINGS(delaunay_module) { /*...*/ }В нём как раз и описывается что мы хотим экспортировать в модуль. Проще всего описать экспорт функций. Здесь это сделано одной строчкой

function("build_triangulation", &build_triangulation);

Смысл её совершенно понятен. Указываем имя функции, по которому она будет вызываться извне, а также ссылку на непосредственную реализацию этой функции. Здесь мы используем промежуточную функцию build_triangulation, которая принимает и возвращает аргументы типа emscripten::val. Этот тип - это обёртка над теми данными, которые передаются в модуль и забираются из него. Назначение функции build_triangulation состоит в следующем:

  • Преобразовать пришедшие данные в правильный формат. Команда convertJSArrayToNumberVector<float> возвращает std::vector<float>

  • Сформировать из полученного набора чисел массив точек

  • Вызвать для этого массива точек нашу готовую функцию triangulate

  • Преобразовать полученный массив целых чисел в то, что снаружи будет понято как Int32Array. Это делается тремя командами

    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);

В общем-то ничего сложного. Однако, как мы видим, требуется некоторая работа по преобразованию данных.

Для того, чтобы экспортировать класс BVHNode мы сделаем дополнительный класс-обёртку, который будет содержать нужные нам конструкторы и только те методы, которые мы хотим экспортировать. Итак, создаём класс BVHNodeWrapper

BVHNodeWrapper
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;
};

В нём реализуем два конструктора. Первый принимает на вход только массив с координатами точек. Мы уже умеем его преобразовывать в формат std::vector<float>. Далее строим триангуляцию, потом треугольники и упаковываем их в BVH, сохраняя ссылку на BVH во внутренней переменной.

Первый конструктор класса BVHNodeWrapper
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);
}

Другой конструктор принимает массив с координатами и массив с индексами вершин треугольников. В этом случае просто разбираем их и строим из треугольников BVH.

Второй конструктор класса BVHNodeWrapper
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);
}

Ну и, наконец, метод sample. На вход он принимает два числа. Строим из них точку, вызываем метод sample у объекта, сохранённого во внутренней переменной. Далее формируем ответ в виде массива координат вершин (или пустого массива, если точка не принадлежит никакому треугольнику).

Метод sample класса BVHNodeWrapper
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;
}

Описание экспорта класса выглядит следующим образом

class_<BVHNodeWrapper>("BVHNode")
    .constructor<emscripten::val>()
    .constructor<emscripten::val, emscripten::val>()
    .function("sample", &BVHNodeWrapper::sample);

Здесь указано, чтобы класс BVHNodeWrapper экспортировался с именем BVHNode. Потом сказано, какие у него конструкторы и какие методы (он тут один).

Теперь скомпилируем наш файл delaunay_api.cpp в модуль WebAssembly. Для этого надо использовать специальный компилятор Emscripten. Чтобы не использовать никаких дополнительных сущностей, сделаем это вручную.

Итак, сначала запускаем файл emcmdprompt.bat из директории со скаченным Emscripten-ом. Это для тех, у кого Windows. У кого Linux - те сами разберутся, что запускать. В результате откроется окно терминала с уже прописанными всеми необходимыми путями для поиска зависимостей.

Теперь в запущенном терминале переходим в директорию с проектом. Далее каждый из файлов delaunay.cpp и bvh.cpp компилируем в промежуточный формат *.o. Чтобы не загромождать директорию проекта лучше создать отдельную папку \output\libs\ и сохранять результат туда. Делается это командами

emcc .\delaunay.cpp -o output\libs\delaunay.o -c
emcc .\bvh.cpp -o output\libs\bvh.o -c

Синтаксис простой: после emcc указывается что компилировать, после ключа -o куда, а ключ -c говорит, что надо создать объектный файл. Теперь, наконец, мы готовы скомпилировать наш модуль. Для этого выполняем следующую команду

emcc .\delaunay_api.cpp output\libs\delaunay.o output\libs\bvh.o -o output\delaunay.js -lembind -s MODULARIZE -s ALLOW_MEMORY_GROWTH=1 -s EXPORT_NAME=delaunay -Oz

Структура команды такая:

  • После emcc указывается основной файл для компиляции

  • Далее указываются все ранее созданные объектные файлы delaunay.o и bvh.o. Из них Emscripten будет брать реализации необходимых функций

  • После ключа -o указывается куда сохранить итоговый модуль и с каким именем. Указывается имя файла с расширением *.js, рядом будет создан собственно модуль с таким же именем, но расширением *.wasm

  • Ключ -lembind говорит, что надо использовать embind

  • Далее идёт набор параметров, каждый из которых начинается с ключа -s

    • MODULARIZE позволяет потом использовать модуль в NodeJS, подключая его с помощью require

    • ALLOW_MEMORY_GROWTH=1 говорит, что в процессе работы модулю разрешается увеличивать объём используемой памяти. Такое требуется, когда выполняются объёмные вычисления. Без этого параметра во время вычислений может возникнуть ошибка переполнения памяти

    • EXPORT_NAME задаёт имя экспортированного модуля

  • ключ -Oz задаёт максимальную степень оптимизации с упором на размер. Без этого результирующий *.wasm файл получится уж больно большой

Вот и всё. Непосредственно файл с модулем delaunay.wasm получился 43 килобайта, а файл с обвесом delaunay.js всего ничего - аж 57 килобайт. И это уже минифицированный. Ну уж что получилось.

Вот пример использования:

var factory = require('../cpp/output/delaunay.js');

factory().then((instance) => {
	const coordinates = [
		-0.5, 3.0,
		8.0, -3.5,
		-7.0, 3.0,
		2.0, -10.0,
		7.0, 8.0,
		8.0, 5.5];
	const triangle_indices = instance.build_triangulation(coordinates);
	console.log(triangle_indices);
	
	let bvh = new instance.BVHNode(coordinates);
	let s = bvh.sample(0.0, 0.0);
	console.log(s);
});

В качестве ответа выдаёт в точности то же самое, что и чистая программа на c++. Работает правильно.

AssemblyScript

Если кто не слыхал, то AssemblyScript - это специальный язык программирования, который компилируется только в WebAssembly, и ни для чего другого не предназначен. По своему синтаксису он очень похож на TypeScript.

Установить AssemblyScript можно через npm. Удобно это будет сделать глобально. Всё-таки это не та вещь, которая меняется от проекта к проекту.

По большому счёту этого достаточно, что компилировать модули WebAssembly. Всё, что теперь требуется - это вызвать команду

asc index.ts

Здесь index.ts это файл с программой. Содержимое откомпилированного модуля в текстовом формате будет выведен в консоль.

То, что AssemblyScript предназначен исключительно для WebAssembly накладывает некоторые неудобства. Они связаны с процедурой написания и отладки более или менее комплексных программ. Ну вот как, например, наша. Тяжко каждый раз компилировать модуль, вызывать его из заранее заготовленной программки на JavaScript, смотреть что он выдаёт, чертыхаться. Поэтому предлагается использовать следующую организацию рабочей среды.

Инициализируем пустой проект.

npm init

После этого устанавливаем локально @assemblyscript/wasi-shim. Так и пишем

npm install -D @assemblyscript/wasi-shim

Это позволит использовать WASI изнутри откомпилированного модуля. В частности - выводить в консоль всякие текстовые сообщения, не совершая никаких дополнительных действий. Как раз то самое "оно там само".

Файлы с кодом программы на AssemblyScript принято складывать в папку./assembly/, а результаты компиляции в папку ./build/. Поэтому создаём в директории ./assembly/ файл index.ts. Это будет точка входа при разработке. Содержимое такое:

function main(): void {
    console.log("Hello World");
}

main();

Теперь, чтобы включить в модуль функции для выполнения команды console.log, его надо откомпилировать с дополнительными параметрами. В общем пишем

asc assembly/index.ts -o build/index.wasm --target debug --config ./node_modules/@assemblyscript/wasi-shim/asconfig.json

Наконец можно запустить наш модуль index.wasm в какой-нибудь среде выполнения, например в WasmTime. Для этого никакой обвес на JavaScript не нужен.

wasmtime ./build/index.wasm

Каждый раз писать команду для компиляции, а потом ещё одну для запуска в WasmTime хлопотно. Поэтому лучше в конфигурационном файле проекта package.json заранее сделать заготовки

"scripts": {
    "develop": "asc assembly/index.ts -o build/index.wasm --target debug --config ./node_modules/@assemblyscript/wasi-shim/asconfig.json",
    "execute": "npm run develop && wasmtime ./build/index.wasm"
},

Теперь компиляция и запуск программы происходит нажатием одной кнопки npm run execute. Можно программировать.

Программа на AssemblyScript

Общая логика ровно та же самая, что была для программы на c++. Создаём в директории ./assembly/ три файла point.ts delaunay.ts и bvh.ts

Файл point.ts содержит описание класса Point для задания координат точек на плоскости.

points.ts
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() + ")";
    }
}

Файл delaunay.ts содержит все функции для построения триангуляции. Самая главная функция это triangulate, которая принимает на вход StaticArray<Point>, содержащий точки на плоскости, а возвращает StaticArray<i32>, содержащий номера вершин треугольников триангуляции.

delaunay.ts
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);
}

В файл bvh.ts описываем класс BVHNode. В общем всё как и раньше.

bvh.ts
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;
        }
    }
}

Проверим, как работают наши функции, использовав их на проверочных данных в файле index.ts

index.ts
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();

Выполним npm run execute. В результате получим

2,0,3,0,2,4,3,0,1,1,0,5,0,4,5
<(-7.0, 3.0), (-0.5, 3.0), (2.0, -10.0)>

Это правильный ответ.

Компиляция модуля

К сожалению на настоящий момент AssemblyScript поддерживает экспорт только функций. Раньше можно было и классы, при этом для использования откомпилированного модуля использовался один универсальный загрузчик. А потом разработчики решили сместить немного акценты, отказались от универсального загрузчика в пользу генерации JavaScript обвеса, предоставляющего доступ только к указанным экспортированным сущностям. Ну и на настоящий момент этот обвес генерируется только для функций и констант.

Ограничение на использование только функций по большому счёту ни на что не влияет. Просто вместо создания экземпляра класса надо вызывать функцию, которая создаёт этот экземпляр и возвращает его (а на самом деле лишь указатель на него в памяти), а вместо вызова метода класса надо вызывать функцию, которая принимает указатель на объект, вызывает у него метод и возвращает результат.

Ещё раз определимся, что мы хотим экспортировать в модуль

  • Функцию для построения триангуляции. Она должна принимать простой числовой массив, задающий координаты точек, а возвращать массив целых чисел, описывающий треугольники

  • Две функции для построения BVH. Первая из них должна принимать массив с координатами точек и возвращать указатель на построенный объект BVH. Вторая функция - принимать два массива: с координатами точек и с номерами вершин треугольников, а возвращать точно так же указатель на объект BVH.

  • Функцию поиска треугольника, содержащего данную точку. Она должна принимать указатель на объект BVH и два числа, а возвращать либо пустой массив, если точка не принадлежит никакому треугольнику, либо массив из шести чисел, задающих вершины треугольника, содержащего точку.

Короче, создаём файл delaunay_api.ts

delaunay_api.ts
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);
    }
}

Все функции, которые мы хотим экспортировать в модуль WebAssembly, должны быть помечены словом export. Первая из нужных нам функций

function build_triangulation(coordinates: Float32Array): Int32Array

В этой функции для передачи данных используются массивы Float32Array и Int32Array. Обвес, который генерирует AssemblyScript, как раз и будет конвертировать массивы, которые передаются со стороны JavaScript вот в эти, и обратно

Следующие две нужных нам функции

function build_bvh(coordinates: Float32Array): BVHNode

function build_bvh_triangulation(coordinates: Float32Array, triangles: Int32Array): BVHNode

Они возвращают указатели на объект. Со стороны JavaScript это выглядит как обычное целое число. Поэтому напрямую в JavaScript с этими объектами работать нельзя.

Четвёртая функция

function sample_bvh(bvh: BVHNode, x: f32, y: f32): Float32Array

Компилируем

asc assembly/delaunay_api.ts -o build/delaunay.wasm --bindings esm -O

Смысл почти всех параметров команды понятен. --bindings esm говорит о том, что надо сгенерировать JavaScript обвес. -O указывает на уровень оптимизации. При компиляции итогового модуля мы не хотим, чтобы в него включались функции для поддержки WASI, поэтому никаких дополнительных конфигурационных файлов не указываем. Ну и сам модуль во время работы конечно же не должен вызывать никаких команд вроде console.log.

Можно добавить ещё один флаг компиляции --uncheckedBehavior always. В результате скорость выполнения модуля увеличивается примерно на 20%. Это достигается за счёт того, что при выполнении модуля пропускаются все проверки на корректность индексов элементов массивов. В общем, в модуле, откомпилированном с таким ключом, можно без сообщений об ошибке взять пятый элемент трёхэлементного массива. Ну а то, что в таком случае программа будет работать некорректно - так это уже дело программиста всё так спланировать, чтобы подобного не происходило. Но зато работает быстрее.

В результате получим файл delaunay.wasm размером 19 килобайт и delaunay.js размером 6 килобайт (и это не минифицированный).

Пример использования

import * as exports from "./delaunay.js";

const coordinates = [
	-0.5, 3.0,
	8.0, -3.5,
	-7.0, 3.0,
	2.0, -10.0,
	7.0, 8.0,
	8.0, 5.5];
const triangle_indices = exports.build_triangulation(coordinates);
console.log(triangle_indices);

let bvh = exports.build_bvh(coordinates);
let s = exports.sample_bvh(bvh, 0.0, 0.0);
console.log(s);

Ответ получается тоже правильный.

wasm-pack и Rust

Будем придерживаться того же подхода, что и в случае c++. Представим, что нам надо просто написать программу для построения триангуляции Делоне на Rust-е. Прошу не серчать за простецкий код, не особо использующий идеоматические вычурные конструкции Rust-а. Компилятор Rust-а настолько крутой, что выдаёт быстрые результаты вне зависимости от того, как это написано. То, что нам, неучам, и надо.

Программа на Rust

Ничтоже сумняшеся создаём пустой проект. Пока даже сторонние крэйты никакие не нужны.

cargo init

Дальше в папке ./src/ создаём три файла point.rs, delaunay.rs и bvh.rs. Содержимое этих файлов примерно такое же, что и раньше (но только на Rust-е).

point.rs
#[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}
    }
}

delaunay.rs
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
}

bvh.rs
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
        }
    }
}

Никакой специфики тут особо нет. Ну кроме той, что привносит сам Rust.

В файле main.rs проверяем, как всё работает

mod point;
mod delaunay;
mod bvh;

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

fn main() {
    let points = vec![
        Point::new(-0.5, 3.0),
        Point::new(8.0, -3.5),
        Point::new(-7.0, 3.0),
        Point::new(2.0, -10.0),
        Point::new(7.0, 8.0),
        Point::new(8.0, 5.5)
        ];
    let triangle_indices = triangulate(&points);
    println!("{:?}", triangle_indices);

    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())));
    }

    let tree = BVHNode::new(triangles);
    let s = tree.sample(&Point::new(0.0, 0.0));

    println!("{:?}", s);
}

Всё должно работать нормально.

Использование wasm-pack

Сам wasm-pack можно установить с официального сайта. Так-то у него много возможностей, но мы будем использовать лишь некоторые. Подход к компиляции модуля WebAssembly тут примерно такой же, какой был при использовании Emscripten. Нам потребуется создать некоторое количество дополнительных обёрток и функций, в которых, в частности, реализовать преобразования между данными, которые приходят извне, и которые используются в наших уже написанных программах.

Итак, создаём в директории ./src/ файл lib.rs. В нём-то и будем всё делать. Но теперь нам надо добавить в конфигурационный файл Cargo.toml разъяснения, чтобы Rust не путался. Добавляем такие строчки

[lib]
name = "delaunay"
path = "src/lib.rs"
crate-type = ["cdylib", "rlib"]

[[bin]]
name = "delaunay_app"
path = "src/main.rs"

Это говорит Rust-у, что в файле main.rs начинается бинарный крэйт, а в файле lib.rs - библиотечный. Ещё укажем две зависимости, которые нам понадобятся

[dependencies]
wasm-bindgen = "0.2.84"
js-sys = "0.3.61"

В отличие от AssemblyScript и подобно Emscripten wasm-pack позволяет экспортировать в модуль классы. Есть только ограничение, что эти классы могут иметь лишь один конструктор. Поэтому сейчас нам надо экспортировать следующее:

  • Функцию для построения триангуляции

  • Класс BVHNode с конструктором, принимающим массив координат точек на плоскости

  • Функцию, которая возвращает объект класса BVHNode, построенного по поданным на вход двум массивам, содержащих как координаты точек, так и номера вершин треугольников

  • Метод класса BVHNode, возвращающего координаты вершин треугольника, содержащего данную точку

В общем пишем такое

lib.rs
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()
            },
        }
    }
}

Основные шаги, что тут написано.

Сначала говорим

use wasm_bindgen::prelude::*;

Эта команда делает что надо.

Потом описываем функцию для построения триангуляции

#[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()
}

В этой функции стоит отметить три вещи:

  • Она помечается атрибутом #[wasm_bindgen], который означает, что функция будет экспортирована в модуль WebAssembly

  • Входной аргумент имеет тип &[f32]. Обвязка, которая генерируется wasm-pack, без проблем конвертирует числовые массивы JavaScript в этот тип

  • Тип выходных данных Array. Этот тип представляет собой что-то вроде обобщённого массива, который может содержать разные данные, и который преобразуется в массив JavaScript. Мы преобразуем массив типа Vec<usize> в Array с помощью команды

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

Возможно вместо результата Array стоит возвращать Uint32Array. Крэйт js-sys содержит много всяких дополнительных структур, предназначенных для обмена данными с JavaScript-ом. Но оставим как есть.

Для того, чтобы экспортировать класс, мы сначала объявляем обёртывающую структуру

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

Она помечается атрибутом #[wasm_bindgen], у которого указан дополнительный параметр js_name = BVHNode. Он задаёт имя класса внутри модуля.

Потом описываем как будто бы методы экспортированного класса

#[wasm_bindgen(js_class = BVHNode)]
impl BVHNodeWrapper {
}

Снова используем параметр js_class = BVHNode атрибута wasm_bindgen чтобы указать имя класса, для которого будут реализованы методы. Все функции, которые содержатся внутри области impl будут экспортированы в модуль. Каждую из них помечать отдельно не надо.

Чтобы пометить метод как конструктор класса используется атрибут #[wasm_bindgen(constructor)]

#[wasm_bindgen(constructor)]
pub fn new(coordinates: &[f32]) -> BVHNodeWrapper {
    /*...*/
}

Далее статическая функция вместо второго конструктора класса. У неё нет аргумента &self

pub fn new_trianglulation(coordinates: &[f32], triangles: &[i32]) -> BVHNodeWrapper{
    /*...*/
}

И, наконец, метод для поиска треугольника по точке. Он возвращает Array, как мы уже умеем делать

pub fn sample(&self, x: f32, y: f32) -> Array {
    /*...*/
}

Тут, возможно, кто-то снова захочет возвращать не обобщённый Array, а плоский Float32Array. Это можно. Для этого в самом начале пишем

use js_sys::Float32Array;

Потом при наличии валидного сэмпла создаём массив длины 6 (ведь треугольник описывается шестью числами, по два на координаты каждой вершины)

let to_return = Float32Array::new_with_length(6);

После этого заполняем компоненты массива координатами вершин

to_return.set_index(0, t.a.x);
to_return.set_index(1, t.a.y);
// ... and so on

Вот и всё. Теперь запускаем команду

wasm-pack build

В результате в директории ./pkg/ будут созданы в том числе два файла: delaunay.js объёмом 7 кб и delaunay_bg.wasm объёмом 31 кб. По умолчанию сразу делается релизный билд, так что никаких дополнительных параметров указывать не надо. Единственное, можно указать дополнительный параметр --target nodejs если хочется сделать модуль для использования в NodeJS, и --target web, если для использования в браузере. Хотя, конечно, у команды wasm-pack есть много других полезных параметров и возможностей.

Итого

У нас есть три модуля WebAssembly, которые делают одно и то же, но получены разными способами. Самое время их сравнить. Но при этом стоит помнить, что любое сравнение - от лукавого, и его результаты могут вызвать противоречивые чувства.

Диаграмма размеров скомпилированных *.wasm файлов. Чем меньше - тем лучше.

Мы запустили небольшой бенчмарк для проверки быстродействия модулей в среде NodeJS. Эксперимент состоял в том, чтобы нагенерировать кучу точек, а потом попросить все три модуля построить триангуляции для разного числа точек. Ну и замерялось время вызова функции построения триангуляции. Внутренние преобразования данных между JavaScript и модулем учитывались.

Вот график быстродействия для разного числа точек. Чем меньше - тем лучше.

А вот график аналогичного эксперимента, но по построению BVH для различного числа точек. Снова, чем меньше значение - тем лучше.

Выводы такие. AssemblyScript хороший, потому что генерирует самый маленький размер модуля, и работает не так чтобы уж сильно медленнее остальных. Rust хороший, потому что пусть размер его модуля и не самый маленький, зато работает быстрее всех. Emscripten тоже молодец.

Мне больше нравится AssemblyScript, так как он кажется концептуально проще по сравнению с Rust-ом.

Бонус 1: демонстрация

В этом небольшом разделе продемонстрируем как можно использовать модуль, откомпилированный из Rust-а, в маленьком вроде как web-приложении. Это web-приложение делает следующее:

  • По нажатию кнопки генерирует случайный набор точек на плоскости

  • Строит триангуляцию по построенным точкам и выводит её на экран

  • При перемещении указателя манипулятора типа "мышь" закрашивает отдельным цветом тот треугольник, над которым находится указатель

Вот маленькая гифка с демонстрацией

А вот ссылка на само приложение.

Так как приложение маленькое, то нет резону разбивать его на много файлов и использовать широкий инструментарий web-технологий. Напишем всё по-простому в одном html файле.

index.html
<!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>

Тут стоит отметить только одно - как подключать модуль при загрузке web-страницы.

Создаём

<script type="module">
</script>

Внутри него пишем

import init, { BVHNode, build_triangulation } from "./delaunay.js";

async function wasm_run() {
    await init();
    /*...*/
}

После завершения асинхронного выполнения init() можно использовать функции модуля.

Бонус 2: интеграция в Python

Обвес, который генерирует Rust, настолько читабельный, что можно попробовать портировать его на Python. Это позволит использовать WebAssembly модуль из Python.

Есть несколько модулей для Python, которые позволяют загружать модули WebAssembly. Все они предоставляют низкоуровневый доступ в том числе к функциям, который экспортируются из модуля, и к памяти. А больше нам ничего и не требуется. Будем использовать Wasmer.

pip isntall wasmer

Прежде чем переходит к написанию кода, посмотрим повнимательнее на модуль, который откомпилировал Rust. Для этого его надо любой подходящей утилитой преобразовать в текстовый формат. А можно даже и не преобразовывать, я сам это сделаю. В общем там будет много строчек, но нам надо два фрагмента: почти в самом начале и почти в самом конце:

(import "wbg" "__wbindgen_number_new" (func (;0;) (type 23)))
(import "wbg" "__wbindgen_object_drop_ref" (func (;1;) (type 5)))
(import "wbg" "__wbg_new_b525de17f44a8943" (func (;2;) (type 11)))
(import "wbg" "__wbg_newwithlength_0da6f12fbc1ab6eb" (func (;3;) (type 3)))
(import "wbg" "__wbg_set_17224bc548dd1d7b" (func (;4;) (type 2)))
(import "wbg" "__wbg_push_49c286f04dd3bf59" (func (;5;) (type 1)))
(import "wbg" "__wbindgen_throw" (func (;6;) (type 0)))

(export "build_triangulation" (func 68))
(export "wbg_bvhnode_free" (func 43))
(export "bvhnode_new" (func 52))
(export "bvhnode_new_trianglulation" (func 48))
(export "bvhnode_sample" (func 55))
(export "__wbindgen_malloc" (func 63))

После ключевого слова import написаны имена функций, которые должны быть реализованы в Python и переданы в модуль. После слова export наоборот указаны имена тех функций, которые можно вызывать из модуля.

Итак, создаём файл delaunay_wasm.py и пишем в нём заготовку для класса загрузчика

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

class Loader:
    def __init__(self, wasm_path: str):
        pass

Единственный входной аргумент для конструктора - путь до файла *.wasm.

Теперь создаём инстанс WebAssembly модуля и сохраняем его в переменную self._wasm

class Loader:
    def __init__(self, wasm_path):
        store = Store()
        module = Module(store, open(wasm_path, "rb").read())
        imports = ImportObject()
        self._wasm = Instance(module, imports)

Этот код работать не будет, так как в переменной 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()

Переменная self._uint32_memory представляет собой ссылку на один большой массив, который интерпретирует память модуля как набор 32-битных чисел. Аналогичным образом self._float32_memory интерпретирует память модуля как массив 32-битных чисел с плавающей точкой.

Второе - сохраним в отдельные переменные ссылки на все нужные нам доступные функции. Значение wasm_exports является итератором по которому можно проитерироваться и получить пары (имя функции, функция)

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

Теперь можно начинать подготовку к реализации импортов. Создаём в конструкторе нашего класса дополнительные переменные

self._heap = [None] * 128
self._heap.extend([None, None, True, False])
self._heap_next = len(self._heap)
self._wasm_vector_len = 0

После этого пишем реализацию шести функций: _add_heap_object, _get_object, _drop_object, _take_object, _pass_array_f32_to_wasm и _pass_array_32_to_wasm

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

Эти реализации - почти что дословный порт аналогичных функций из обвязки, которую генерирует Rust.

Теперь функции для импортов

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"

У аргументов этих функции обязательно должны быть аннотации типов, в противном же случае Wasmer будет ругаться. И ещё стоит отметить, что аргумент самой первой функции импорта _import_number_new должен быть типа f64. Мы указываем его явно в виде строки, так как тип float интерпретируется Wasmer-ом как f32.

Подключаем эти импорты к модулю. В конструкторе загрузчика пишем

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),
})

Теперь всё должно работать без ошибок. Ну как работать - если в основной программе написать

loader = Loader("delaunay_bg.wasm")

то всё должно выполниться без ошибок. Хотя, конечно, пока ничего, кроме загрузки модуля, и не делается

Самое время реализовать внутри класса Loader вызов функции для построения триангуляции

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]

Теперь сделаем отдельный класс BVHNode

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))

Конструктор этого класса принимает ссылку на объект загрузчика и два массива. Загрузчик используется для вызова функций модуля. При этом если указаны оба входных массива, то используется метод создания объекта по координатам точек и треугольникам, а если только один - то только по координатам точек.

Помимо этого реализован метод __del__, который автоматически вызывается питоном в случае удаления объекта. Мы в свою очередь при этом должны вызвать соответствующую функцию модуля для удаления этого объекта из памяти.

Всё. Вот итоговый файл delaunay_wasm.py

delaunay_wasm.py
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))

Модуль WebAssembly можно использовать

loader = Loader("delaunay_bg.wasm")
coordinates = [-0.5, 3.0, 8.0, -3.5, -7.0, 3.0, 2.0, -10.0, 7.0, 8.0, 8.0, 5.5]
t = loader.build_triangulation(coordinates)
print(t)

bvh = BVHNode(loader, coordinates, t)
print(bvh.sample(0.0, 0.0))

В результате будет выведено

[2, 0, 3, 0, 2, 4, 3, 0, 1, 1, 0, 5, 0, 4, 5]
[-7.0, 3.0, -0.5, 3.0, 2.0, -10.0]

Это правильный ответ.

По скорости Wasmer выполняет код примерно в 1.3 - 1.4 раза медленнее, чем NodeJS. В целом неплохо. С этой точки зрения можно воспринимать WebAssembly как платформо независимую альтернативу обычным бинарным питоновским модулям, которые обычно компилируются из того же c++ для ускорения тяжёлых вычислений.