-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy pathmath.js
More file actions
953 lines (826 loc) · 33.4 KB
/
math.js
File metadata and controls
953 lines (826 loc) · 33.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
/*
* 常用空间分析函数: hfhan
* 包含 点、线、面 之间的相互关系
* 本文主要记录一些自己开发中常用的一些分析函数,比较正宗和全面的可以看一些空间分析库,比如前端的 Turf 和 JSTS
*/
'use strict';
// 点 -------------------------------------------------------------------------
// 点到点的距离
export function dist2d(coord1, coord2) {
let dx = coord1[0] - coord2[0];
let dy = coord1[1] - coord2[1];
return Math.sqrt(dx * dx + dy * dy)
}
//dist2d([1, 1], [4, 5]) // 5
//判断两个点集是否相等
export function equals(coord1, coord2) {
let equals = true;
for (let i = coord1.length - 1; i >= 0; --i){
if (coord1[i] != coord2[i]) {
equals = false;
break
}
}
return equals
}
//equals([1, 3], [1, 3]) // true
//equals([1, 3, 1], [1, 3, 2]) // false
// 线 -------------------------------------------------------------------------
// 线段的长度
export function formatLength(coords) {
coords = coords || []
let length = 0
//通过遍历坐标计算两点之前距离,进而得到整条线的长度
for (let i = 0, leng = coords.length - 1; i < leng; i++) {
length += dist2d(coords[i], coords[i + 1]);
}
return length
}
//formatLength([[1, 0], [2, 1], [3, 0], [4, 2]]) // 5.06449510224598
// 根据偏移量偏移一条线段
export function lineOffset(coords, deltaX, deltaY){
deltaX = deltaX || 0
deltaX = isNaN(deltaX) ? 0 : deltaX
deltaY = deltaY || 0
deltaY = isNaN(deltaY) ? 0 : deltaY
if(deltaX == 0 && deltaY == 0)return coords
coords.forEach(coord => {
coord[0] += deltaX;
coord[1] += deltaY;
})
return coords
}
//lineOffset([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]], 2, 3)
// [[3, 4], [4, 6], [5, 6], [6, 5], [4, 3]]
//根据距离和角度偏移一条线段
export function offsetLine(line, d, angle){
//let line = [[1, 1], [4, 2], [5, 5], [7, 6]]
//let d = 1
//let angle = 30
//斜边长度d已知,角度angle已知
//对边长度就是y的偏移量 就是 d * sin(angle) ==> d * Math.sin(angle * Math.PI / 180)
//邻边长度就是x的偏移量 就是 d * cos(angle) ==> d * Math.cos(angle * Math.PI / 180)
let ox = d * Math.cos(angle * Math.PI / 180)
let oy = d * Math.sin(angle * Math.PI / 180)
return line.map(coords => [coords[0] + ox, coords[1] + oy])
}
//offsetLine([[1, 1], [4, 2], [5, 5], [7, 6]], 5, 30)
//[[5.330127018922194,3.4999999999999996],[8.330127018922195,4.5],[9.330127018922195,7.5],[11.330127018922195,8.5]]
// 线段上距离点P最近的一个点
export function getNearestCoord(point, lines) {
let d, res = { dist: Infinity }
if(!(Array.isArray(lines[0]) && Array.isArray(lines[0][0]))){
lines = [lines]
}
lines.forEach(function(coords){
for (let i = 0; i < coords.length; i++) {
d = dist2d(point, coords[i]);
if (d < res.dist) {
res.dist = d;
res.index = i
res.point = coords[i]
}
}
})
return res.point ? res : null;
}
//getNearestCoord([2.2, 3.1], [[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// {dist: 0.5099019513592785, index: 1, point: [2, 3]}
/*
* 点P到线段AB的最短距离
* 使用矢量算法,计算线AP在线段AB方向上的投影,当需要计算的数据量很大时,这种方式优势明显
* 特殊情况如点在线段上、点在端点、点在线段延长线上等等的情况全部适用于此公式,只是作为特殊情况出现,无需另作讨论。这也是矢量算法思想的优势所在。
* 函数返回值:point 投影坐标 dist 点P到投影距离 type 垂足位置,不为0表示垂足在线段外
*/
export function pointToSegmentDist(point, point1, point2){
let x = point[0], x1 = point1[0], x2 = point2[0]
let y = point[1], y1 = point1[1], y2 = point2[1]
//线段AB 为一个点
if(x1 == x2 && y1 == y2) return {
type: 0,
point: point1,
dist: 0
}
let cross = (x2 - x1) * (x - x1) + (y2 - y1) * (y - y1);
//let r = cross / d2
//r < 0 点P的垂足在线段AB外,且点P距离线段AB最近的点为A
//r = 0 点P的垂足和点P距离线段AB最近的点为A
if (cross <= 0) return {
type: 1,
point: point1,
dist: Math.sqrt((x - x1) * (x - x1) + (y - y1) * (y - y1))
};
let d2 = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
//r > 1 点P的垂足在线段AB外,且点P距离线段AB最近的点为B
//r = 1 点P的垂足和点P距离线段AB最近的点为B
if (cross >= d2) return {
type: 2,
point: point2,
dist: Math.sqrt((x - x2) * (x - x2) + (y - y2) * (y - y2))
};
let r = cross / d2;
let px = x1 + (x2 - x1) * r;
let py = y1 + (y2 - y1) * r;
return {
type: 0,
point: [px, py],
dist: Math.sqrt((x - px) * (x - px) + (py - y) * (py - y))
};
}
//pointToSegmentDist([1, 3], [0, 1], [3, 1]);
//{ dist: 2, point: [1, 1], type: 0 }
/*
* pointToSegmentDist变种:点P到直线AB的垂足及最短距离
*/
export function pointToFootDist(point, point1, point2){
let x = point[0], x1 = point1[0], x2 = point2[0]
let y = point[1], y1 = point1[1], y2 = point2[1]
//线段AB 为一个点
if(x1 == x2 && y1 == y2) return {
type: true,
point: point1,
dist: 0
}
let dx = x2 - x1;
let dy = y2 - y1;
//r < 0 点P的垂足在线段AB外,且点P距离线段AB最近的点为A
//r > 1 点P的垂足在线段AB外,且点P距离线段AB最近的点为B
//r = 0 点P的垂足和点P距离线段AB最近的点为A
//r = 1 点P的垂足和点P距离线段AB最近的点为B
let r = (dx * (x0 - x1) + dy * (y0 - y1)) / (dx * dx + dy * dy || 0);
let px = x1 + dx * r;
let py = y1 + dy * r;
return {
type: r >= 0 && r <= 1, //true 垂足在线段内 false 垂足在线段外
point: [px, py], //垂足
dist: Math.sqrt((x - px) * (x - px) + (py - y) * (py - y)) //点P到垂足距离
};
}
//pointToFootDist([3, 1], [1, 0], [3, 0])
// {dist: 1, point: [3, 0], type: true}
/*
* 已知三角形ABC三个顶点坐标(P1, P2, P3),底边BC上有一点D,且AD长度r已知,求D点坐标
* BC的长度为d2,设点D到B点的比例为k,则D点坐标为
* x = (x3 - x2) * k + x2
* y = (y3 - y2) * k + y2
* 根据AD距离得出公式:(x1 - x)² + (y1 - y)² = r²
* 代入x,y 得:((y3 - y2)² + (x3 - x2)²) * k² + (2 * (x3 - x2) * (x2 - x1) + 2(y3 - y2) * (y2 - y1)) * k + (x1² + x2² + y1² + y2² - 2 * x1 * x2 - 2 * y1 * y2 - r²) = 0
* 一元二次方程的解得出k:a * x² + b * x + c = 0 ==> X = (-b ± (b² - 4ac) ^ 1/2) / 2a
*/
export function getPointByDist(p1, p2, p3, r){
if(equals(p2, p3)){
return p2
}
let x1 = p1[0], x2 = p2[0], x3 = p3[0]
let y1 = p1[1], y2 = p2[1], y3 = p3[1]
let a = (y3 - y2) ** 2 + (x3 - x2) ** 2
let b = 2 * (x3 - x2) * (x2 - x1) + 2 * (y3 - y2) * (y2 - y1)
let c = x1 ** 2 + x2 ** 2 + y1 ** 2 + y2 ** 2 - 2 * x1 * x2 - 2 * y1 * y2 - r ** 2
let z = b ** 2 - 4 * a * c
if(z < 0){
return null
}
/* if(z == 0){
let k = -b / (2 * a)
let x = (x3 - x2) * k + x2
let y = (y3 - y2) * k + y2
return [x, y]
} */
let k1 = (-b + Math.sqrt(z)) / (2 * a)
let k2 = (-b - Math.sqrt(z)) / (2 * a)
let kk = Math.min(k1 < 0 ? Infinity : k1, k2 < 0 ? Infinity : k2)
if(kk > 1){
return null
}
let x = (x3 - x2) * kk + x2
let y = (y3 - y2) * kk + y2
return [x, y]
}
//getPointByDist([3, 4], [1, 1], [4, 0], 3.6) //[1.0202273020830888, 0.9932575659723037]
/*
* 计算点P到直线AB的距离
* 使用三角形面积,计算线AP在直线AB方向上的投影
* Area = |(1/2)(x1 * y2 + x2 * y3 + x3 * y1 - x2 * y1 - x3 * y2 - x1 * y3)|
* Bottom = Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
* Area = .5 * Bottom * H
* Height = Area / .5 / Bottom
*/
export function verticalDistance(point, point1, point2){
if(equals(point1, point2)){
return dist2d(point, point1)
}
let area = Math.abs(0.5 * (point1[0] * point2[1] + point2[0] * point[1] + point[0] * point1[1] - point2[0] * point1[1] - point[0] * point2[1] - point1[0] * point[1]));
let bottom = dist2d(point1, point2)
let height = area / 0.5 / bottom;
return height;
}
//verticalDistance([4, 2], [1, 0], [3, 0]) // 2
/*
* 计算点P到直线AB的距离
* 根据公式,计算线AP在直线AB方向上的投影
* d = Math.abs((A * x + B * y + C) / Math.sqrt( A * A + B * B))
* 公式介绍:https://zhuanlan.zhihu.com/p/26307123
* 直线公式:y = kx + b
* 直线公式:kx - y + b = 0
* 假设 B = -1,则 A = k,C = b
* 直线斜率:k = (y1 - y2) / (x1 - x2)
* 常数计算:b = y - kx
*/
export function verticalDistance2(point, point1, point2){
//如果point1[0] == point2[0] 说明是条竖着的线
if(point1[0] == point2[0]){
return Math.abs(point[0] - point1[0])
}else{
let k = (point1[1] - point2[1]) / (point1[0] - point2[0])
let b = point1[1] - k * point1[0]
return Math.abs((k * point[0] - point[1] + b) / Math.sqrt( k * k + 1))
}
}
//verticalDistance2([4, 2], [1, 0], [3, 0]) // 2
/*
* 计算纬度和经度指定的两点之间的距离(米)
* 在计算涉及到经纬度的计算时(长度、面积等),如果只是粗略的计算,那么我们使用坐标点计算,将计算结果按照每经纬度111km转换即可。这个结果越往两极失真越大
* 如果想要精细点,就需要用到地球长半轴参数去计算
* WGS84 坐标系下的长半轴参数 radius = 6378137 m
* WGS84(World Geodetic System 1984) 坐标系是为GPS全球定位系统使用而建立的坐标系统。是国际公认的坐标系
* 北京54坐标系下的长半轴参数 radius = 6378245 m
* 1954年北京坐标系可以认为是前苏联1942年坐标系的延伸。它的原点不在北京而是在前苏联的普尔科沃。所以误差较大,缺点较多。
* 西安80坐标系下的长半轴参数 radius = 6378140 ± 5 m
* 1980西安坐标系所采用的IAG1975椭球,其长半轴要比WGS84椭球长半轴的值大3米左右,而这可能引起地表长度误差达10倍左右
* CGCS2000(2000 国家大地坐标系) 长半轴参数 radius = 6378137 m
* 2000国家大地坐标系,是我国当前最新的国家大地坐标系,是全球地心坐标系在我国的具体体现,其原点为包括海洋和大气的整个地球的质量中心。
* 正弦曲线投影下的长半轴参数 radius = 6370997 m
* 正弦曲线投影是一种等面积的伪圆柱投影。规定纬线投影为平行直线,经线投影为对称于中央经线的正弦曲线,同一纬线上经距相等,纬距向两极缩小。主要用于小比例尺世界地图
* 公式:haversin(d / r) = haversin(φ2 - φ1) + cos(φ2) * haversin(Δλ)
* 其中:R为地球半径,取6378137; φ1, φ2 表示两点的纬度; Δλ 表示两点经度的差值
*/
export function haversineDistance(c1, c2) {
let radius = 6378137
let lat1 = toRadians(c1[1]);
let lat2 = toRadians(c2[1]);
let deltaLatBy2 = (lat2 - lat1) / 2;
let deltaLonBy2 = toRadians(c2[0] - c1[0]) / 2;
let a = Math.sin(deltaLatBy2) * Math.sin(deltaLatBy2) + Math.sin(deltaLonBy2) * Math.sin(deltaLonBy2) * Math.cos(lat1) * Math.cos(lat2);
return 2 * radius * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
}
export function toRadians(_){
//1°=π/180° 1rad=180°/π
//角度转弧度
return _ * Math.PI / 180
}
//haversineDistance([117, 34], [118, 39]) //563732.8911197125米
//haversineDistance([118, 39], [117, 34]) //563732.8911197125米
/*
* 根据中心点坐标和距离(米),计算范围
*/
export function haversineLnglat(orign, dist){
let radius = 6378137
// 求东西两侧的的范围边界。在haversin公式中令φ1 = φ2(维度相同)
let lat = toRadians(orign[1])
let dlng = 2 * Math.asin(Math.sin(dist / (2 * radius)) / Math.cos(lat));
// 弧度转换成角度
dlng = toDegrees(dlng);
// 然后求南北两侧的范围边界,在haversin公式中令 Δλ = 0
let dlat = dist / 6378137;
// 弧度转换成角度
dlat = toDegrees(dlat);
return [orign[0] - dlng, orign[1] - dlat, orign[0] + dlng, orign[1] + dlat]
}
function toDegrees(_){
//1°=π/180° 1rad=180°/π
//弧度转角度
return _ * 180 / Math.PI
}
//haversineLnglat([117, 34], 5000) //[116.94582179841282, 33.955084235794025, 117.05417820158718, 34.044915764205975]
//haversineLnglat([117, 34], 6000) //[116.93498615776213, 33.94610108295283, 117.06501384223787, 34.05389891704717]
/*
* 根据一个点切割线段
* 点不必落在线段上
* 常用于计算鼠标位置距离线段起点或终点的距离
* 简单版的是使用 getNearestCoord,而不是 pointToSegmentDist
* getNearestCoord 寻找的是线上真实存在的一点
* pointToSegmentDist 寻找的点位于线上,但不一定真实存在与线上
*/
export function getClosestPoint(coords, point){
if(!coords || coords.length < 2) return [[], []]
let squaredDistance = {dist: Infinity}, index
for(let i = 1; i < coords.length; i++){
let d = pointToSegmentDist(point, coords[i - 1], coords[i])
if(d.dist < squaredDistance.dist){
squaredDistance = d
index = i
}
}
if(index === undefined){
return [coords, []]
}
let prearr = coords.slice(0, index)
if(prearr.length && !equals(squaredDistance.point, prearr[prearr.length - 1])){
prearr.push(squaredDistance.point)
}
let nextarr = coords.slice(index)
if(nextarr.length && !equals(squaredDistance.point, nextarr[0])){
nextarr.unshift(squaredDistance.point)
}
return [prearr, nextarr]
}
//getClosestPoint([[0, 1], [2, 3], [5, 6], [3, 4], [4, 4]], [3, 5])
/*
[
[[0,1],[2,3],[3.5,4.5]],
[[3.5,4.5],[5,6],[3,4],[4,4]]
]
*/
/*
* 线段与线段的交点
* 解线性方程组, 求线段AB与线段CD的交点坐标,如果没有交点,返回null
*/
export function intersects(coords1, coords2) {
let x1 = coords1[0][0];
let y1 = coords1[0][1];
let x2 = coords1[1][0];
let y2 = coords1[1][1];
let x3 = coords2[0][0];
let y3 = coords2[0][1];
let x4 = coords2[1][0];
let y4 = coords2[1][1];
//斜率交叉相乘 k1 = (y4 - y3) / (x4 - x3) k2 = (y2 - y1) / (x2 - x1)
//k1 k2 同乘 (x4 - x3) * (x2 - x1) 并相减得到denom
let denom = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1));
// 如果分母为0 则平行或共线, 不相交
if (denom === 0) {
return null;
}
let numeA = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3));
let numeB = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3));
let uA = numeA / denom;
let uB = numeB / denom;
// 交点在线段1上,且交点也在线段2上
if (uA >= 0 && uA <= 1 && uB >= 0 && uB <= 1) {
let x = x1 + (uA * (x2 - x1));
let y = y1 + (uA * (y2 - y1));
return [x, y];
}
return null;
}
//intersects([[0,0],[1,1]], [[3,0],[2,1]]) //null
//intersects([[0,0],[1,1]], [[3,0],[0,1]]) //[0.75, 0.75]
/*
* 线段与线段的交点
* 两点式直线公式:(x - x1) / (x2 - x1) = (y - y1) / (y2 - y1)
// --> x(y2 - y1) + y(x1 - x2) + x1(y1 - y2) + y1(x2 - x1) = 0
// A = y2-y1, B = x1-x2, C = x1(y1 - y2) + y1(x2 - x1)
// (x,y) = d2/(d1+d2) * (x3,y3) + d1/(d1+d2) * (x4,y4)
// d2/(d1+d2) = |A*x4 + B*y4 + C|/(|A*x4 + B*y4 + C|+|A*x3 + B*y3 + C|)
// d1/(d1+d2) = |A*x3 + B*y3 + C|/(|A*x4 + B*y4 + C|+|A*x3 + B*y3 + C|)
*/
export function detectIntersection(coords1, coords2) {
let a = coords1[0]
let b = coords1[1]
let c = coords2[0]
let d = coords2[1]
// whether intersect ?
let intersect = (isLeft(a, b, c) ^ isLeft(a, b, d)) && (isLeft(c, d, a) ^ isLeft(c, d, b))
if (!intersect) return null;
let d2 = Math.abs((b[1] - a[1]) * d[0] + (a[0] - b[0]) * d[1] + a[0] * (a[1] - b[1]) + a[1] * (b[0] - a[0]));
let d1 = Math.abs((b[1] - a[1]) * c[0] + (a[0] - b[0]) * c[1] + a[0] * (a[1] - b[1]) + a[1] * (b[0] - a[0]));
let x = d2 / (d1 + d2) * c[0] + d1 / (d1 + d2) * d[0];
let y = d2 / (d1 + d2) * c[1] + d1 / (d1 + d2) * d[1];
return [x, y];
}
function isLeft(point, point1, point2) {
let x1 = point1[0] - point[0]
let y1 = point1[1] - point[1]
let x2 = point2[0] - point[0]
let y2 = point2[1] - point[1]
return (-y1 * x2 + x1 * y2) > 0
}
/*
* 线与线的所有交点
* 接收两条线的点集合,求取line2 和 line1的交点
* 如果规定count参数,代表返回前多少个交点,否则全部返回
* 返回交点集合,集合中的元素属性包括 index1:交点在line1的位置,index2:交点在line2的位置,coords:交点坐标
* 代码参考truf.lineIntersect:http://turfjs.org/docs/#lineIntersect
* 如果交点为线的点,则会重复返回,包括返回的数据结构,都是为下面切割面函数splitPolygon服务的
*/
export function lineIntersect(line1, line2, count){
let result = []
for(let i = 1; i < line1.length; i++){
let coords1 = [line1[i-1], line1[i]]
//求取数据的边界范围,函数放在了下面的多边形中
let bbox1 = bbox(coords1)
for(let j = 1; j < line2.length; j++){
let coords2 = [line2[j-1], line2[j]]
let bbox2 = bbox(coords2)
//判断两个边界范围的关系: bbox1 是否包含 bbox2,函数放在了下面的多边形中
if(isIntersects(bbox1, bbox2)){
let p = intersects(coords1, coords2)
if(p){
result.push({index1: i, index2: j, coords: p})
if(count && (result.length >= count)){
return result
}
}
}
}
}
return result
}
//lineIntersect([[0, 0], [0, 1], [1, 3], [3, 2], [5, 0]], [[-1, 0], [4, 3]])
//[{"index1":1,"index2":1,"coords":[0,0.6]},{"index1":3,"index2":1,"coords":[2.6363636363636367,2.1818181818181817]}]
/*
* 判断线段路径是顺时针还是逆时针
* 使用格林公式计算
* 返回值 d < 0 '顺时针' : d > 0 '逆时针'
* truf对应实现truf.booleanClockwise://turfjs.org/docs/#booleanClockwise
*/
export function isClockwise(line){
let length = line.length
let d = 0
if(length.length < 3)return d;
//如果不是闭合线,则改为闭合线
if(!equals(line[0], line[length - 1])){
length = (line = line.slice()).push(line[0]);
}
//沿着多边形的边求曲线积分,若积分为正,则是沿着边界曲线正方向(逆时针),反之为顺时针
//最后一个点是开始点,只参与末尾点(倒数第二个点)的计算,不单独计算
for(let i = 0; i < length - 1; i++){
d += -0.5 * (line[i + 1][0] - line[i][0]) * (line[i + 1][1] + line[i][1])
}
return d
}
//isClockwise([[0,0],[1,1],[1,0],[0,0]]) // -0.5 顺时针
//isClockwise([[0,0],[1,0],[1,1],[0,0]]) // 0.5 逆时针
/*
* 判断线段路径是顺时针还是逆时针
* 使用端点判断
* 为正时,p1-p2-p3 路径的走向为逆时针,
* 为负时,p1-p2-p3 走向为顺时针,
* 为零时,p1-p2-p3 所走的方向不变,亦即三点在一直线上
* 返回值 d < 0 '顺时针' : d > 0 '逆时针'
*/
export function isClockwise2(line){
let length = line.length
let d = 0
if(length.length < 3)return d;
//如果不是闭合线,则改为闭合线
if(!equals(line[0], line[length - 1])){
length = (line = line.slice()).push(line[0]);
}
//循环遍历多边形的坐标选取X或者Y值中最大或者最小的点,这个点必然是凸点
//然后取该点前后各一个点Pm-1、Pm+1,组成向量(Pm-1,Pm)、(Pm,Pm+1)。然后进行向量叉乘即可判断出顺时针或逆时针
let max = line[0][0], maxIndex = 0;
for(let i = 1; i < length - 1; i++){
if(line[i][0] > max){
maxIndex = i
max = line[i][0]
}
}
//找到端点 p2 ,然后取改点前后各一点 p1 p3
let p1 = line[maxIndex ? (maxIndex - 1) : (length - 2)]
let p2 = line[maxIndex]
let p3 = line[maxIndex + 1]
//然后根据三个点组成的两个向量乘积判断,进行向量叉乘即可判断出顺时针或逆时针
//d = (x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1)
d = (p2[0] - p1[0]) * (p3[1] - p2[1]) - (p3[0] - p2[0]) * (p2[1] - p1[1])
return d
}
//isClockwise2([[0,0],[1,1],[1,0],[0,0]]) // -1 顺时针
//isClockwise2([[0,0],[1,0],[1,1],[0,0]]) // 1 逆时针
//根据连续的三个点,用矢量判断夹角大小和方向
/*
向量点积、向量点乘,又称向量的积(不是向量积)、数量积
两个向量的数量积等于它们对应坐标的乘积的和。即:若a=(x1,y1),b=(x2,y2),则a·b=x1·x2+y1·y2
公式:a·b = xa * xb + ya * yb
公式:a·b = |a||b|·cosθ θ为向量a与向量b的夹角
几何意义:向量a在向量b方向上的投影与向量b的模的乘积。是一个标量
几何意义:数量积a·b等于a的长度|a|与b在a的方向上的投影|b|cosθ的乘积
cosθ = a·b / (|a||b|)
θ = Math.acos(a·b / (|a||b|)) 0 < θ <= 2PI
θ = θ * 180 / Math.PI
公式:|a| = Math.sqrt(xa * xa + ya * ya) 向量a的长度
公式:|b| = Math.sqrt(xb * xb + yb * yb) 向量b的长度
向量叉积、向量叉乘,又称向量积(不是向量的积):
公式:c = a x b = xa * yb - ya * xb
公式:|c| = |a x b| = |a||b|·sinθ
几何意义:c是垂直a、b所在平面,且以|b|·sinθ为高、|a|为底的平行四边形的面积。是一个矢量
向量a × 向量b(×为向量叉乘),
若结果小于0,表示向量b在向量a的顺时针方向;
若结果大于0,表示向量b在向量a的逆时针方向;
若等于0,表示向量a与向量b平行。
*/
export function getAngleBy3Point(point1, point2, point3) {
let xa = point2[0] - point1[0]
let xb = point3[0] - point2[0]
let ya = point2[1] - point1[1]
let yb = point3[1] - point2[1]
//direction 大于0 逆时针, 小于0 顺时针, 等于0 平行
let angle = 0, direction = 0
let _a = Math.sqrt(xa * xa + ya * ya)
let _b = Math.sqrt(xb * xb + yb * yb)
if(_a && _b){
let p = xa * xb + ya * yb
angle = Math.acos(p / (_a * _b))
angle = angle / Math.PI * 180
direction = xa * yb - ya * xb
direction = direction < 0 ? -1 : direction > 0 ? 1 : 0
}
return { angle, direction }
}
// 判断一个多边形是不是凹多边形
// 凸多边形两个相邻的的向量方向应该是一样的(排除平行线时的情况),所以乘积不会小于0,小于0则是凹多边形
export function IsConcavePolygon(points){
let direction = 0
for(let i = 0, l = points.length; i < l; i++){
let res = getAngleBy3Point(points[i == 0 ? (l - 1) : (i - 1)], points[i], points[i == (l - 1) ? 0 : (i + 1)])
//乘积小于0,说明方向不一致,为凹多边形
if(direction * res.direction < 0){
return true
}
direction = direction || res.direction
}
return false
}
// IsConcavePolygon([[0, 0], [0, 5], [5, 5], [2, 2], [5, 0]]) // true
// IsConcavePolygon([[0, 0], [0, 5], [5, 5], [5, 0], [0, 0]]) // false
/*
* 简化线 与 平滑线
* 简化线 与 平滑线 比较复杂,这里就不展示代码了
* 简化线 又可以叫做线的抽稀,一般使用的方法为 道格拉斯-普克,这也是我使用的方法
* 道格拉斯-普克算法(Douglas–Peucker algorithm,亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点,并减少点的数量的一种算法。其思想主要是保留关键点。
* Douglas-Peucker算法描述:
* 1、在线首尾两点A,B之间连接一条直线AB,该直线为曲线的弦;
* 2、得到曲线上离该直线段距离最大的点C,计算其与AB的距离d;
* 3、比较该距离与预先给定的阈值tolerance的大小,如果小于tolerance,则该直线段作为曲线的近似,该段曲线处理完毕。
* 4、如果距离大于阈值tolerance,则用C将曲线分为两段AC和BC,并分别对两段曲线进行1~3的处理。
* 5、当所有曲线都处理完毕时,依次连接各个分割点形成的折线,即可以作为曲线的近似。
* 除了可以使用道格拉斯-普克算法对线进行抽稀外,还可以使用Wang-Müller算法(保留关键折弯)、Zhou-Jones算法(保留加权有效面积)、Visvalingam-Whyatt算法(保留有效面积)等。
* 在实际应用中发现,现实是残酷的,如果阈值过小,简化后的控制点会很多,用户很难进行操作。而阈值变大,控制点变少了,但是经过平滑之后,线面变形严重,比较失真。真实和美观之间很难平衡,难以兼得。
* 平滑线 一般使用的方法为 贝赛尔曲线插值和B样条曲线插值
* 在应用中发现,两者在某些时候都会有些问题,所以自己对贝赛尔曲线算法做了一点改变,具体的可以看https://segmentfault.com/a/1190000031626358
*/
// 多边形 -------------------------------------------------------------------------
// 获取多边形边界范围
export function bbox(coords) {
// x/经度最小值 y/纬度最小值 x/经度最大值 y/纬度最大值
let res = [Infinity, Infinity, -Infinity, -Infinity];
coords.forEach(coord => {
if (res[0] > coord[0]) res[0] = coord[0];
if (res[2] < coord[0]) res[2] = coord[0];
if (res[1] > coord[1]) res[1] = coord[1];
if (res[3] < coord[1]) res[3] = coord[1];
})
return res;
}
//bbox([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [1, 0, 4, 3]
//判断两个边界范围的关系: a 是否包含 b
export function isContains(a, b) {
return a[0] <= b[0] &&
a[1] <= b[1] &&
b[2] <= a[2] &&
b[3] <= a[3];
}
//isContains([1, 0, 4, 3], [2, 2, 5, 5]) //false
//isContains([1, 0, 4, 3], [2, 1, 3, 2]) //true
//判断两个边界范围的关系: a 与 b 是否有交集
export function isIntersects(a, b) {
return b[0] <= a[2] &&
b[1] <= a[3] &&
b[2] >= a[0] &&
b[3] >= a[1];
}
//isIntersects([1, 0, 4, 3], [2, 2, 5, 5]) //true
//isIntersects([1, 0, 4, 3], [5, 2, 5, 5]) //false
// 获取多边形中心
export function center(coords) {
let ext = bbox(coords);
let x = (ext[0] + ext[2]) / 2;
let y = (ext[1] + ext[3]) / 2;
return [x, y];
}
//center([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.5, 1.5]
// 获取多边形重心
export function centroid(coords) {
let xSum = 0, ySum = 0, len = 0;
coords.forEach(coord => {
xSum += coord[0];
ySum += coord[1];
len++;
})
return [xSum / len, ySum / len];
}
//centroid([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.4, 1.8]
// 获取多边形质心
export function centerOfMass(coords) {
let centr = centroid(coords);
let neutralizedPoints = coords.map(function (point$$1) {
return [
point$$1[0] - centr[0],
point$$1[1] - centr[1]
];
});
let sx = 0, sy = 0, sArea = 0;
let x1, x2, y1, y2, a;
for (let i = 0, len = coords.length - 1; i < len; i++) {
x1 = neutralizedPoints[i][0];
y1 = neutralizedPoints[i][1];
x2 = neutralizedPoints[i + 1][0];
y2 = neutralizedPoints[i + 1][1];
// a 是计算有符号面积和最终坐标的公因子
a = x1 * y2 - x2 * y1;
// sArea 用来计算有符号面积的和
sArea += a;
// sx 和 sy 是用于计算最终坐标的总和
sx += (x1 + x2) * a;
sy += (y1 + y2) * a;
}
if (sArea === 0) {
return centr;
} else {
// 计算有符号面积,并因式分解: x = 1 / 6A
let area = sArea * 0.5;
let areaFactor = 1 / (6 * area);
return [
centr[0] + areaFactor * sx,
centr[1] + areaFactor * sy
]
}
}
//centerOfMass([[1, 1], [2, 3], [3, 3], [4, 2], [2, 0]])
// [2.569230769230769, 1.8735042735042735]
/*
* 获取多边形面积
* 原理:△ABC的面积就是'向量AB'和'向量AC'两个向量叉积的绝对值的一半。其正负表示三角形顶点是在右手系还是左手系。
* 所以可以把多边形以多边形内一点为顶点,拆分出一个个三角形,再求取面积
* 设 n 边形的顶点依次是 (x1,y1)(x2,y2)......(xn,yn)
* 那么:s = (x1y2-x2y1)/2 + (x2y3-x3y2)/2 +......+ (xny1-x1yn)/2
* 获取的 s 是有向面积,需要取绝对值
*/
export function getArea(polygon){
let area = 0;
let len = polygon.length
let x1 = polygon[len - 1][0]
let y1 = polygon[len - 1][1]
for(let i = 0; i < len; i++){
let x2 = polygon[i][0];
let y2 = polygon[i][1];
area += y1 * x2 - x1 * y2;
x1 = x2;
y1 = y2
}
return Math.abs(area / 2)
}
//getArea([[0,0],[0,10],[10,10],[10,0]]) //100
//getArea([[0,0],[0,3],[4,0]]) //6
/*
* 计算多边形投影到地球上时的近似面积(平方米)
* 注意,如果环是顺时针方向的,这个区域将是正的,否则它将是负的
* 注意,投影及其长半轴参数radius的介绍这里不做叙述,可以看上面haversineDistance的相关介绍
*/
export function geodesicArea(coords) {
let radius = 6378137
let area = 0, len = coords.length;
if(len <= 2){ return area }
let x1 = coords[len - 1][0];
let y1 = coords[len - 1][1];
for (let i = 0; i < len; i++) {
let x2 = coords[i][0], y2 = coords[i][1];
area += toRadians(x2 - x1) * (2 + Math.sin(toRadians(y1)) + Math.sin(toRadians(y2)));
x1 = x2;
y1 = y2
}
area = area * radius * radius / 2
return Math.abs(area)
}
//geodesicArea([[118, 39], [117, 34], [117, 33], [116, 36], [117, 40]]) //64519860945.307144平方米
/*
* 判断点是否在多边形内
* 射线法(ray casting)或者奇偶规则法(even odd rule)
* 从这个点做一条射线,计算它跟多边形边界的交点个数,如果交点个数为奇数,那么点在多边形内部,否则点在多边形外。
* 对于射线法,需要排除以下几种情况:
* 1、点在多边形的边上
* 2、点和多边形的顶点重合
* 3、射线经过多边形顶点
* 4、射线刚好经过多边形的一条边
* truf对应实现truf.booleanPointInPolygon://turfjs.org/docs/#booleanPointInPolygon
*/
export function pointInPolygon(point, polygon) {
let px = point[0], py = point[1], flag = false
let sx, sy, tx, ty
for(let i = 0, l = polygon.length, j = l - 1; i < l; j = i, i++) {
sx = polygon[i][0]
sy = polygon[i][1]
tx = polygon[j][0]
ty = polygon[j][1]
// 点与多边形顶点重合
if((sx === px && sy === py) || (tx === px && ty === py)) {
return true
}
// 判断线段两端点是否在射线两侧
if((sy < py && ty >= py) || (sy >= py && ty < py)) {
// 线段上与射线 Y 坐标相同的点的 X 坐标
let x = sx + (py - sy) * (tx - sx) / (ty - sy)
// 点在多边形的边上
if(x === px) {
return true
}
// 射线穿过多边形的边界
if(x > px) {
flag = !flag
}
}
}
// 射线穿过多边形边界的次数为奇数时点在多边形内
return flag
}
/*
* 判断点是否在多边形内
* 回转数法
* 用线段分别连接点和多边形的全部顶点,计算所有点与相邻顶点连线的夹角,计算所有夹角和,最后根据角度累加值计算回转数
* 注意每个夹角都是有方向的,所以有可能是负值。360°(2π)相当于一次回转。
* 当回转数为 0 时,点在闭合曲线外部。
*/
export function pointInPolygon2(point, polygon) {
let px = point[0], py = point[1], sum = 0
let sx, sy, tx, ty
for(let i = 0, l = polygon.length, j = l - 1; i < l; j = i, i++) {
sx = polygon[i][0]
sy = polygon[i][1]
tx = polygon[j][0]
ty = polygon[j][1]
// 点与多边形顶点重合
if((sx === px && sy === py) || (tx === px && ty === py)) {
return true
}
// 点与相邻顶点连线的夹角
let angle = Math.atan2(sy - py, sx - px) - Math.atan2(ty - py, tx - px)
// 确保夹角不超出取值范围(-π 到 π)
if(angle >= Math.PI) {
angle = angle - Math.PI * 2
} else if(angle <= -Math.PI) {
angle = angle + Math.PI * 2
}
sum += angle
}
// 计算回转数并判断点和多边形的几何关系
return Math.round(sum / Math.PI) !== 0
}
/*
* 判断点是否在多边形内
* openlayers中的方法,也是使用的射线法,但是他根据https://mentin.medium.com/which-predicate-cb608b470471中的描述
* 使用的判断依据是contains而不是covers,即点位于线上不算做包含
* 但它仍然有问题,下面两个例子要么都为false,要么都为true
* pointInPolygon3([110, 0], [[110, 0], [100, 10], [110, 30], [120, 10]]) => false
* pointInPolygon3([110, 0], [[110, 0], [110, 50], [150, 50], [150, 0]]) => true
*/
export function pointInPolygon3(point, polygon) {
let x = point[0], y = point[1]
let wn = 0;
let x1 = polygon[0][0];
let y1 = polygon[0][1];
for (let i = 0; i < polygon.length; i++) {
let x2 = polygon[i][0];
let y2 = polygon[i][1];
if (y1 <= y) {
if (y2 > y && (x2 - x1) * (y - y1) - (x - x1) * (y2 - y1) > 0){
wn++;
}
} else if (y2 <= y && (x2 - x1) * (y - y1) - (x - x1) * (y2 - y1) < 0){
wn--;
}
x1 = x2;
y1 = y2
}
return wn !== 0
}
/*
* 两个多边形之间的关系:相交、包含、相离
* 0、先判断两个路径是否有交点,是的话就是相交
* 1、遍历A多边形上的点,判断是否有坐标点在B多边形内 --- 返回结果 a
* 2、遍历B多边形上的点,判断是否有坐标点在A多边形内 --- 返回结果 b
* 如果a、b都为true,则两个多边形相交
* 如果a为true,b为false,则多边形B包含多边形A
* 如果a为false,b为true,则多边形A包含多边形B
* 如果a、b都为false,则两个多边形远离
* truf对应实现 包含 truf.booleanContains://turfjs.org/docs/#booleanContains
* truf对应实现 交叉 truf.booleanCrosses://turfjs.org/docs/#booleanCrosses
* truf对应实现 相离 truf.booleanDisjoint://turfjs.org/docs/#booleanDisjoint
*/
export function judge(coordsA, coordsB){
//先判断线有没有交点
let points = lineIntersect(coordsA, coordsB, 1)
if(points.length){
return 3 //相交
}
let boola = coordsA.some(item => {
return pointInPolygon3(item, coordsB)
}) ? 1 : 0;
let boolb = coordsB.some(item => {
return pointInPolygon3(item, coordsA)
}) ? 1 : 0;
//return ['相离', 'A包含B', 'B包含A', '相交'][boola * 2 + boolb]
return boola * 2 + boolb
}