最近有需求要求计算六边形，我研究了一下有几种方式
- h3.js
- turf.js
实现的思路就是先使用h3的函数或者turf.hexgrid来计算出来相关的多边形，然后按照多边形和六边形的关系做一次filter
使用过后，发现以下问题
1.h3.js实现的效率比较高,但是生成的多边形是斜的,因此表现效果上不佳
2.使用turf生成的hexgrid实现的效率很低，很多多边形的时候会很卡，>4s以上，但是生成的多边形方方正正
因此便有了疑问h3.js为何效率比较高，我看了一下h3.js的源码,libh3.js明显是webasembly实现的方式,所以我决定把turf的hexgrid函数用c实现一遍,
1.turf hexgrid源码的翻译
//turf.cc
#ifndef EM_PORT_API
# if defined(__EMSCRIPTEN__)
# include<emscripten.h>
# if defined(__cplusplus)
# define EM_PORT_API(rettype) extern "C" rettype EMSCRIPTEN_KEEPALIVE
# else
# define EM_PORT_API(rettype) rettype EMSCRIPTEN_KEEPALIVE
# endif
# else
# if defined(__cplusplus)
# define EM_PORT_API(rettype) extern "C" rettype
# else
# define EM_PORT_API(rettype) rettype
# endif
# endif
#endif
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
double* hexes;//存储多边形
double min(double x, double y) {
return x < y ? x : y;
}
EM_PORT_API(void) free_space(){
if(hexes){
free(hexes);
}
hexes=nullptr;
}
double _abs(double x) {
return x > 0 ? x : -1 * x;
}
double max(double x, double y) {
return x > y ? x : y;
}
const double DEGREE2RADIANS = 0.017453292519943295;
const double EARTH_RADIUS = 6371008.8;
EM_PORT_API(double) distance(double lon1, double lat1, double lon2, double lat2) {
double radLat1 = lat1 * DEGREE2RADIANS;
double radLat2 = lat2 * DEGREE2RADIANS;
double a = radLat1 - radLat2;
double b = lon1 * DEGREE2RADIANS - lon2 * DEGREE2RADIANS;
double s = 2 * asin(sqrt(pow(sin(a / 2), 2) +
cos(radLat1) * cos(radLat2) * pow(sin(b / 2), 2)));
s = s * EARTH_RADIUS;
return s;
}
struct Point {
double x;
double y;
Point(double x, double y) :x(x), y(y) {}
Point() :x(0), y(0) {}
};
EM_PORT_API(bool) pointinpolygon(double x1, double y1, double* polygonarr, int count) {
//A point is in a polygon if a line from the point to infinity crosses the polygon an odd number of times
bool odd = false;
for (int i = 0, j = count - 1; i < count; i++) {
if (((polygonarr[2*i+1] > y1) != (polygonarr[2*j+1] > y1)) // One point needs to be above, one below our y coordinate
// ...and the edge doesn't cross our Y corrdinate before our x coordinate (but between our x coordinate and infinity)
&& (x1 < (polygonarr[2*j] - polygonarr[2*i]) * (y1 - polygonarr[2*i+1]) / (polygonarr[2*j+1] - polygonarr[2*i+1]) + polygonarr[2*i])) {
// Invert odd
// System.out.println("Point crosses edge " + (j + 1));
// totalCrosses++;
odd = !odd;
}
//else {System.out.println("Point does not cross edge " + (j + 1));}
j = i;
}
return odd;
}
const double consines[6] = { 1,
0.5000000000000001,
-0.4999999999999998,
-1,
-0.5000000000000004,
0.49999999999999933 };
const double sines[6] = {
0,
0.8660254037844386,
0.8660254037844387,
1.2246467991473532e-16,
-0.8660254037844385,
-0.866025403784439
};
double arr[14];
EM_PORT_API(double*) hexagon(double lon, double lat, double rx, double ry) {
for (int i = 0; i < 6; i++) {
double x = lon + rx * consines[i];
double y = lat + ry * sines[i];
arr[i * 2] = x;
arr[i * 2 + 1] = y;
}
arr[12] = arr[0];
arr[13] = arr[1];
return arr;
}
//六边形的算法
EM_PORT_API(double*) prehexgrid(double* bbox, double cellSide,double* mask,int maskCount,int* cc) {
double west = bbox[0];
double south = bbox[1];
double east = bbox[2];
double north = bbox[3];
double centerY = (south + north) / 2.0;
double centerX = (west + east) / 2.0;
double xFraction = (cellSide * 2) / distance(west, centerY, east, centerY);
double cellWidth = xFraction * (east - west);
double yFraction = (cellSide * 2) / distance(centerX, south, centerX, north);
double cellHeight = yFraction * (north - south);
double radius = cellWidth / 2;
double hex_width = radius * 2;
double hex_height = sqrt(3) / 2 * cellHeight;
double box_width = east - west;
double box_height = north - south;
double x_interval = (3.0 / 4) * hex_width;
double y_interval = hex_height;
double x_span = (box_width - hex_width) / (hex_width - radius / 2);
int x_count = int(floor(x_span));
double x_adjust = (x_count * x_interval - radius / 2 - box_width) / 2 -
radius / 2 +
x_interval / 2;
// adjust box_height so all hexagons will be inside the bbox
int y_count = int(floor((box_height - hex_height) / hex_height));
double y_adjust = (box_height - y_count * hex_height) / 2;
bool hasOffsetY = y_count * hex_height - box_height > hex_height / 2;
if (hasOffsetY) {
y_adjust -= hex_height / 4;
}
printf("%d,%d,%ld", x_count, y_count, (x_count + 1) * (y_count + 1) * 14 * sizeof(double));
if (hexes) free(hexes);
hexes = (double*)malloc((x_count + 1) * (y_count + 1) * 14 * sizeof(double));
int index = 0;
*cc = 0;
for (int x = 0; x <= x_count; x++) {
for (int y = 0; y <= y_count; y++) {
bool isOdd = (x % 2 == 1);
if (y == 0 && isOdd) continue;
if (y == 0 && hasOffsetY) continue;
double center_x = x * x_interval + west - x_adjust;
double center_y = y * y_interval + south + y_adjust;
if (isOdd) {
center_y -= hex_height / 2;
}
double* hex = hexagon(center_x, center_y, cellWidth / 2, cellHeight / 2);
bool inPoly = false;
for (int i = 0; i < 6; i++) {
if (pointinpolygon(hex[i * 2], hex[i * 2 + 1], mask, maskCount)) {
inPoly = true;
break;
}
}
if (inPoly) {
for (int i = 0; i < 7; i++) {
hexes[index] = hex[i * 2];
hexes[index+1] = hex[i * 2 + 1];
index += 2;
}
(*cc)++;
}
}
}
return hexes;
}
int main() {
double arr[4] = {
120.124391,
30.136342,
120.236837,
30.2376661 };
double mask[606] = {
120.221058,30.237666,120.218373,30.235087,120.213307,30.230547,120.211297,30.228716,120.208612,30.226514,120.203348,30.222039,120.196561,30.216567,120.187584,30.211086,120.182867,30.207929,120.177622,30.205625,120.165029,30.20131,120.148694,30.195716,120.142721,30.19364,120.140885,30.193007,120.140174,30.192838,120.136944,30.19197,120.134121,30.190587,120.131764,30.188743,120.129668,30.18685,120.12804,30.184778,120.126671,30.18283,120.125951,30.18111,120.125111,30.179444,120.124626,30.177352,120.124391,30.172199,120.12459,30.167261,120.125128,30.163251,120.126001,30.159457,120.127474,30.156009,120.130377,30.151164,120.133044,30.147956,120.138594,30.142077,120.146144,30.137174,120.148451,30.138966,120.151484,30.141541,120.151779,30.142047,120.152336,30.142541,120.152773,30.142747,120.153497,30.142937,120.154209,30.143183,120.154693,30.143925,120.154942,30.144073,120.155892,30.144333,120.155305,30.145231,120.154992,30.145199,120.154744,30.145958,120.155311,30.146005,120.155667,30.146189,120.156763,30.147137,120.157946,30.147951,120.158722,30.14868,120.159077,30.148889,120.159903,30.148223,120.160536,30.147474,120.160276,30.147127,120.159646,30.147089,120.159665,30.146401,120.159794,30.146044,120.160135,30.146121,120.161411,30.146764,120.162004,30.147104,120.162314,30.147387,120.162887,30.148108,120.163804,30.148916,120.164114,30.149848,120.164218,30.150019,120.164363,30.150089,120.164815,30.150059,120.165645,30.149592,120.165798,30.149784,120.166304,30.149541,120.166459,30.149609,120.167039,30.149247,120.167748,30.149061,120.168275,30.148647,120.169313,30.149234,120.169687,30.14953,120.16992,30.1499,120.170471,30.151121,120.171031,30.151651,120.171291,30.151832,120.172924,30.152478,120.173428,30.15277,120.174529,30.153646,120.175204,30.153887,120.176357,30.152938,120.1767,30.152381,120.177282,30.151913,120.177653,30.151412,120.178158,30.150979,120.178569,30.150268,120.178862,30.149422,120.180332,30.149818,120.180668,30.149326,120.18116,30.148164,120.181388,30.147826,120.181723,30.147243,120.18314,30.146317,120.184041,30.145388,120.184139,30.145363,120.186651,30.146597,120.186758,30.14656,120.187069,30.145964,120.186185,30.144433,120.186265,30.144335,120.186562,30.144447,120.187224,30.144965,120.188097,30.145145,120.188179,30.145304,120.188479,30.145566,120.188746,30.145586,120.190779,30.146556,120.191541,30.146825,120.191791,30.14716,120.19188,30.147118,120.192554,30.145647,120.19322,30.144793,120.193958,30.144826,120.195352,30.142458,120.196023,30.143079,120.196262,30.143133,120.196385,30.142888,120.196015,30.142564,120.196031,30.142421,120.196501,30.142238,120.196781,30.141628,120.196842,30.141218,120.19699,30.140893,120.197165,30.140764,120.198019,30.140516,120.199551,30.13969,120.200471,30.139436,120.201378,30.13956,120.202359,30.139871,120.202959,30.140163,120.204085,30.141057,120.204503,30.142364,120.204752,30.142711,120.205039,30.142921,120.205262,30.143025,120.205641,30.142611,120.205885,30.142215,120.207219,30.140855,120.207884,30.139947,120.208046,30.139301,120.209109,30.138965,120.209555,30.138712,120.20952,30.138274,120.208739,30.13756,120.208691,30.137251,120.208742,30.137131,120.210281,30.136504,120.210951,30.136342,120.211569,30.136376,120.213295,30.137893,120.214533,30.136915,120.214714,30.136913,120.215921,30.137828,120.217398,30.138647,120.218973,30.139728,120.21947,30.139948,120.219436,30.140051,120.218801,30.140584,120.218286,30.141561,120.216933,30.143075,120.216922,30.143299,120.21706,30.143445,120.217618,30.143524,120.217699,30.143706,120.217644,30.144917,120.21833,30.145888,120.218789,30.147601,120.218904,30.14778,120.219234,30.148027,120.220174,30.148424,120.220334,30.148576,120.220309,30.148731,120.21986,30.149268,120.219853,30.149489,120.220016,30.150049,120.219899,30.150514,120.220105,30.150769,120.220585,30.151225,120.22057,30.151896,120.221061,30.152526,120.221659,30.153759,120.223046,30.153813,120.223807,30.153504,120.224188,30.153495,120.22443,30.153591,120.224592,30.153866,120.224576,30.154195,120.224379,30.155129,120.224175,30.155523,120.223138,30.15669,120.223169,30.157257,120.223239,30.157441,120.223479,30.157708,120.225774,30.159177,120.225853,30.159582,120.226174,30.160155,120.226082,30.160374,120.225331,30.161289,120.225099,30.161799,120.225024,30.162066,120.225043,30.162656,120.224975,30.162712,120.22377,30.162598,120.223136,30.162777,120.222762,30.162807,120.2219,30.16272,120.221269,30.162432,120.220538,30.162581,120.220299,30.162559,120.220059,30.162396,120.219817,30.162012,120.219629,30.161868,120.219116,30.161748,120.218652,30.161799,120.218535,30.161576,120.217946,30.161303,120.217869,30.161315,120.217727,30.161521,120.217591,30.161542,120.216859,30.16132,120.216749,30.161433,120.216719,30.162347,120.216568,30.162382,120.21655,30.162712,120.218127,30.163232,120.218035,30.163641,120.218222,30.163788,120.218224,30.164162,120.217985,30.164506,120.217946,30.164753,120.218469,30.165197,120.218689,30.165493,120.219519,30.166093,120.21948,30.166831,120.219966,30.167029,120.219943,30.1672,120.219362,30.168161,120.219357,30.168367,120.219582,30.168605,120.219455,30.168706,120.218892,30.1688,120.218863,30.168879,120.219003,30.16894,120.219655,30.168925,120.2203,30.169686,120.220315,30.169952,120.219936,30.169766,120.219997,30.170025,120.22055,30.170523,120.220884,30.17057,120.221182,30.170495,120.221388,30.170606,120.221718,30.170632,120.22204,30.17084,120.22223,30.171062,120.224395,30.17125,120.224375,30.171648,120.224058,30.172268,120.224013,30.172775,120.223553,30.173695,120.223417,30.173849,120.223236,30.173891,120.222285,30.173643,120.221875,30.17442,120.221398,30.175512,120.222014,30.175496,120.222031,30.176295,120.222045,30.178883,120.223863,30.178899,120.224656,30.179167,120.224682,30.179393,120.224447,30.181476,120.224576,30.181885,120.233498,30.18204,120.233677,30.183505,120.23531,30.1831,120.236733,30.182979,120.236837,30.18458,120.233982,30.185901,120.234065,30.187528,120.234198,30.188784,120.234249,30.195259,120.234317,30.201283,120.234399,30.208242,120.234512,30.209543,120.23442,30.20999,120.234459,30.213544,120.234952,30.213663,120.235111,30.225426,120.235232,30.228138,120.235041,30.229126,120.234435,30.229798,120.233707,30.229988,120.231887,30.230181,120.227398,30.232514,120.226516,30.232896,120.22396,30.2363,120.221058,30.237666 };
int count = 0;
prehexgrid(arr, 50,mask,606/2,&count);
printf("\n%d", count);
bool c=pointinpolygon(120.2316,30.2130, mask, 606/2);
printf("%d", c);//in
return 0;
}
这里记一个小插曲，本来我用chatgpt来生成一个点在多边形的模板，它给我一个java模板，我一看注释写的很全，但是采用后发现还是有很多的问题
public class PolygonUtils {
/**
* 判断点是否在多边形内
*
* @param point 检测点
* @param pts 多边形的顶点
* @return 点在多边形内返回true, 否则返回false
*/
public static boolean isPointInPoly(Point2D.Double point, List<Point2D.Double> pts) {
int N = pts.size();
boolean boundOrVertex = true; //如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回true
int intersectCount = 0;//cross points count of x
double precision = 2e-10; //浮点类型计算时候与0比较时候的容差
Point2D.Double p1, p2;//neighbour bound vertices
Point2D.Double p = point; //当前点
p1 = pts.get(0);//left vertex
for (int i = 1; i <= N; ++i) {//check all rays
if (p.equals(p1)) {
return boundOrVertex;//p is an vertex
}
p2 = pts.get(i % N);//right vertex
if (p.x < Math.min(p1.x, p2.x) || p.x > Math.max(p1.x, p2.x)) {//ray is outside of our interests
p1 = p2;
continue;//next ray left point
}
if (p.x > Math.min(p1.x, p2.x) && p.x < Math.max(p1.x, p2.x)) {//ray is crossing over by the algorithm (common part of)
if (p.y <= Math.max(p1.y, p2.y)) {//x is before of ray
if (p1.x == p2.x && p.y >= Math.min(p1.y, p2.y)) {//overlies on a horizontal ray
return boundOrVertex;
}
if (p1.y == p2.y) {//ray is vertical
if (p1.y == p.y) {//overlies on a vertical ray
return boundOrVertex;
} else {//before ray
++intersectCount;
}
} else {//cross point on the left side
double xinters = (p.x - p1.x) * (p2.y - p1.y) / (p2.x - p1.x) + p1.y;//cross point of y
if (Math.abs(p.y - xinters) < precision) {//overlies on a ray
return boundOrVertex;
}
if (p.y < xinters) {//before ray
++intersectCount;
}
}
}
} else {//special case when ray is crossing through the vertex
if (p.x == p2.x && p.y <= p2.y) {//p crossing over p2
Point2D.Double p3 = pts.get((i + 1) % N); //next vertex
if (p.x >= Math.min(p1.x, p3.x) && p.x <= Math.max(p1.x, p3.x)) {//p.x lies between p1.x & p3.x
++intersectCount;
} else {
intersectCount += 2;
}
}
}
p1 = p2;//next ray left point
}
if (intersectCount % 2 == 0) {//偶数在多边形外
return false;
} else { //奇数在多边形内
return true;
}
}
}
总之，这段代码翻译成c后有点问题,后面有时间再看
//计算点是否在多边形内
EM_PORT_API(bool) pointinpolygon2(double x1, double y1, double* polygonarr, int count) {
bool boundOrVertex = true; //如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回true
int intersectCount = 0;//cross points count of x
double precision = 2e-11; //浮点类型计算时候与0比较时候的容差
struct Point p1, p2;
struct Point p(x1, y1);//当前点
p1.x = polygonarr[0];
p1.y = polygonarr[1];//leftVertex
for (int i = 1; i <= count; i++) {
if (_abs(p.x - p1.x) < precision && _abs(p.y - p1.y) < precision) {
return boundOrVertex;// p is an vertex
}
p2.x = polygonarr[i * 2];
p2.y = polygonarr[i * 2 + 1];//right vertex
if (p.x<min(p1.x, p2.x) || p.x>max(p1.x, p2.x)) {
//ray is outside of our intersets
p1.x = p2.x;
p1.y = p2.y;
continue;
}
if (p.x > min(p1.x, p2.x) && p.x < max(p1.x, p2.x)) {
if (p.y <= max(p1.y, p2.y)) {
if (p1.x == p2.x && p.y >= min(p1.y, p2.y)) {
return boundOrVertex;
}
if (p1.y == p2.y) {//ray is vertical
if (p1.y == p.y) {//overlies on a vertical ray
return boundOrVertex;
}
else {//before ray
++intersectCount;
}
}
else {//cross point on the left side
double xinters = (p.x - p1.x) * (p2.y - p1.y) / (p2.x - p1.x) + p1.y;
if (_abs(p.y - xinters) < precision) {
return boundOrVertex;
}
if (p.y < xinters) {
++intersectCount;
}
}
}
}
else {//special case when ray is crossing through the vertex
if (p.x == p2.x && p.y <= p2.y) {//p crossing over p2
Point p3(polygonarr[(i + 1) * 2], polygonarr[(i + 1) * 2 + 1]);
if (p.x >= min(p1.x, p3.x) && p.x <= max(p1.x, p3.x)) {//p.x lies between p1.x & p3.x
++intersectCount;
}
else {
intersectCount += 2;
}
}
}
p1.x = p2.x;
p1.y = p2.y;//next ray left point
}
if (intersectCount % 2 == 0) {//偶数在多边形外
return false;
}
else { //奇数在多边形内
return true;
}
}
编译
使用docker apiaryio/emcc 1.38.11镜像,之所以用这个镜像,是因为我需要导出malloc,free函数，因为传入参数中有数组，便于开辟heap空间，传递参数
docker pull apiaryio/emcc 1.38.11
docker run --rm -it -v /root/cpp/:/src apiaryio/emcc:1.38.11
这里我们设置allow_memory_growth用于在动态调整内存大小, wasm参数为0时生成wasm文件，wasm=1时生成js文件【asm.js】文件，
emcc turf.cpp -s ALLOW_MEMORY_GROWTH=1 -s WASM=0 -o libhexgrid.html
这里我们为何要生成asm.js格式文件呢？
1.经过实践，我发现放到vue项目中，使用wasm文件，需要wasmModule.onRuntimeInitialized回调后才能完备使用，实际上就是要等一段时间，在vue中嵌入起来打包等过程中困扰了我好多天，我最终放弃这个wasm的若干好处，还是参考h3.js使用方式
// @ts-nocheck libh3-browser.js
var libh3 = (
function(libh3) {
libh3 = libh3 || {};
var Module=typeof libh3!==...//asm.js代码
return libh3
}
)(typeof libh3 === 'object' ? libh3 : {});
export default libh3;
我们也这样搞个对象导出
var libhexgrid=(
libhexgrid=libhexgrid||{}
\\生成代码
return libhexgrid
)(typeof libhexgrid === 'object' ? libhexgrid : {});
export default libhexgrid;
然后直接import使用即可
这里有个优化就是这个生成代码有NODE,SHELL环境判断代码，可以删除掉，但是也不是必须
经过了上诉的总总折腾，我们提升了hexgrid的效率在5-10倍左右，实现了高效率和图形美观
最后贴一段如何使用传递参数,主要是将数据塞进数组，然后计算出结果后再组装=这个不是c语言强项，样例中使用的是高德坐标系，地图是84的，所以有点偏差
cal_data() {
let data = wasmModule._distance(120, 30, 120, 31);
console.log(data)
console.log(Mockdata)
let bbox =turf.bbox(Mockdata)
// let bbox = [119, 30, 120, 31]
let bboxPtr = wasmModule._malloc(8 * bbox.length);
for (let i = 0; i < bbox.length; i++) {
wasmModule.HEAPF64[bboxPtr / 8 + i] = bbox[i];
}
let mask=Mockdata.features[0].geometry.coordinates.flat(4)
// let mask = [
// 119, 30,
// 120, 30,
// 120, 31,
// 119, 31,
// 119, 30
// ];
console.log(JSON.stringify(mask))
let maskPtr = wasmModule._malloc(8 * mask.length);
for (let i = 0; i < mask.length; i++) {
wasmModule.HEAPF64[maskPtr / 8 + i] = mask[i];
}
let inPtr = wasmModule._malloc(4);
wasmModule.HEAP32[inPtr / 4] = 0;
let resultPtr = wasmModule._prehexgrid(bboxPtr, 50, maskPtr, mask.length/2, inPtr);
let count = wasmModule.HEAP32[inPtr >> 2];
console.log(count)
let hexes = []
for (let i = 0; i < count; i++) {
let hex = []
for (let j = 0; j < 7; j += 1) {
let x = wasmModule.HEAPF64[resultPtr / 8 + 14 * i + j * 2]
let y = wasmModule.HEAPF64[resultPtr / 8 + 14 * i + j * 2 + 1]
hex.push([x, y])
}
// hexes.push(hex)
let poly = turf.polygon([hex], { cnt: 0 })
hexes.push(poly);
// console.log(wasmModule.HEAPF64[resultPtr>>3+i])
}
console.log(JSON.stringify(turf.featureCollection(hexes)))
this.map.getSource('maine').setData(turf.featureCollection(hexes))
wasmModule._free(maskPtr);
wasmModule._free(inPtr);
wasmModule._free(bboxPtr);
wasmModule._free_space();
}
