|
@@ -11,8 +11,9 @@
|
|
|
#include "studio_geo_algo_c.h"
|
|
|
#include "geography/studio_proj_c.h"
|
|
|
|
|
|
+/////////////////// 基础适量算法 ///////////////////
|
|
|
|
|
|
-double distance_squared(const studio_point_c *p1, const studio_point_c *p2)
|
|
|
+double distance_squared_c(const studio_point_c *p1, const studio_point_c *p2)
|
|
|
{
|
|
|
const double dx = p1->x - p2->x;
|
|
|
const double dy = p1->y - p2->y;
|
|
@@ -50,7 +51,103 @@ double point_line_dist_square_c(const studio_point_c *p, const studio_point_c *s
|
|
|
return dist_square;
|
|
|
}
|
|
|
|
|
|
-// Douglas-Peucker算法,用于简化多边形或曲线
|
|
|
+
|
|
|
+double azimuth_angle_c(const studio_point_c *p_from, const studio_point_c *p_to)
|
|
|
+{
|
|
|
+ double delta_x = p_to->x - p_from->x;
|
|
|
+ double delta_y = p_to->y - p_from->y;
|
|
|
+
|
|
|
+ double theta_rad = atan2(delta_x, delta_y);
|
|
|
+ double theta_deg = theta_rad * (180.0 / AO_M_PI);
|
|
|
+
|
|
|
+ if (theta_deg < 0.0)
|
|
|
+ {
|
|
|
+ theta_deg += 360.0;
|
|
|
+ }
|
|
|
+ return theta_deg;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+double distance_squared_gcs_c(const studio_point_c *p1, const studio_point_c *p2)
|
|
|
+{
|
|
|
+ // 经纬度转高斯投影
|
|
|
+ // double central = ((int) (p1->x) / 3) * 3;
|
|
|
+ double central = floor((p1->x + 1.5) / 3.0) * 3.0;
|
|
|
+ double gx_1, gy_1;
|
|
|
+ lonlat_to_gauss(central, p1->x, p1->y, &gx_1, &gy_1);
|
|
|
+ studio_point_c gs_p1 = studio_point_init(gx_1, gy_1);
|
|
|
+ double gx_2, gy_2;
|
|
|
+ lonlat_to_gauss(central, p2->x, p2->y, &gx_2, &gy_2);
|
|
|
+ studio_point_c gs_p2 = studio_point_init(gx_2, gy_2);
|
|
|
+ double dist_square = distance_squared_c(&gs_p1, &gs_p2);
|
|
|
+ return dist_square;
|
|
|
+}
|
|
|
+
|
|
|
+double point_line_dist_square_gcs_c(const studio_point_c *p, const studio_point_c *start, const studio_point_c *end)
|
|
|
+{
|
|
|
+ // 经纬度转高斯投影
|
|
|
+ // double central = ((int) (p->x) / 3) * 3;
|
|
|
+ double central = floor((p->x + 1.5) / 3.0) * 3.0;
|
|
|
+ double gx_p, gy_p;
|
|
|
+ lonlat_to_gauss(central, p->x, p->y, &gx_p, &gy_p);
|
|
|
+ studio_point_c gs_p = studio_point_init(gx_p, gy_p);
|
|
|
+ double gx_start, gy_start;
|
|
|
+ lonlat_to_gauss(central, start->x, start->y, &gx_start, &gy_start);
|
|
|
+ studio_point_c gs_start = studio_point_init(gx_start, gy_start);
|
|
|
+ double gx_end, gy_end;
|
|
|
+ lonlat_to_gauss(central, end->x, end->y, &gx_end, &gy_end);
|
|
|
+ studio_point_c gs_end = studio_point_init(gx_end, gy_end);
|
|
|
+ double dist_square = point_line_dist_square_c(&gs_p, &gs_start, &gs_end);
|
|
|
+ return dist_square;
|
|
|
+}
|
|
|
+
|
|
|
+double azimuth_angle_gcs_c(const studio_point_c *p_from, const studio_point_c *p_to)
|
|
|
+{
|
|
|
+ // 经纬度转高斯投影
|
|
|
+ // double central = ((int) (p_from->x) / 3) * 3;
|
|
|
+ double central = floor((p_from->x + 1.5) / 3.0) * 3.0;
|
|
|
+ double gx_from, gy_from;
|
|
|
+ lonlat_to_gauss(central, p_from->x, p_from->y, &gx_from, &gy_from);
|
|
|
+ studio_point_c gs_from = studio_point_init(gx_from, gy_from);
|
|
|
+ double gx_to, gy_to;
|
|
|
+ lonlat_to_gauss(central, p_to->x, p_to->y, &gx_to, &gy_to);
|
|
|
+ studio_point_c gs_to = studio_point_init(gx_to, gy_to);
|
|
|
+ double theta_deg = azimuth_angle_c(&gs_from, &gs_to);
|
|
|
+ return theta_deg;
|
|
|
+}
|
|
|
+
|
|
|
+studio_point_c lonlat_move_by(const studio_point_c *pt, double delta_x, double delta_y)
|
|
|
+{
|
|
|
+ // 中央经线
|
|
|
+ double central = floor((pt->x + 1.5) / 3.0) * 3.0;
|
|
|
+
|
|
|
+ // 转换为高斯平面坐标
|
|
|
+ double gx, gy;
|
|
|
+ lonlat_to_gauss(central, pt->x, pt->y, &gx, &gy);
|
|
|
+
|
|
|
+ // 应用平面坐标偏移(东向加delta_x,北向加delta_y)
|
|
|
+ const double new_gx = gx + delta_x; // 单位:米
|
|
|
+ const double new_gy = gy + delta_y;
|
|
|
+
|
|
|
+ // 转回经纬度
|
|
|
+ double new_lon, new_lat;
|
|
|
+ gauss_to_lonlat(central, new_gx, new_gy, &new_lon, &new_lat);
|
|
|
+
|
|
|
+ // 经度范围修正(-180~180)
|
|
|
+ new_lon = fmod(new_lon + 360.0, 360.0); // 先转0~360
|
|
|
+ if (new_lon > 180.0)
|
|
|
+ {
|
|
|
+ new_lon -= 360.0; // 再转-180~180
|
|
|
+ }
|
|
|
+ studio_point_c new_point = studio_point_init(new_lon, new_lat);
|
|
|
+
|
|
|
+ return new_point;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+/////////////////// 矢量(线)抽稀算法 ///////////////////
|
|
|
+
|
|
|
+
|
|
|
void douglas_peucker_c(const studio_line_c *points, size_t start, size_t end, double epsilon, unsigned int *indices,
|
|
|
unsigned int *index_count)
|
|
|
{
|
|
@@ -204,7 +301,8 @@ bool line_vacuate_gcs_c(const studio_line_c *line, const int max_points, const d
|
|
|
}
|
|
|
// 经纬度转高斯克吕格投影
|
|
|
studio_line_c gauss_line = studio_line_c_init();
|
|
|
- double central = (int) (line->data[0].x) / 3 * 3;
|
|
|
+ // double central = ((int) (line->data[0].x) / 3) * 3;
|
|
|
+ double central = floor((line->data[0].x + 1.5) / 3.0) * 3.0;
|
|
|
for (unsigned int i = 0; i < line->size; i++)
|
|
|
{
|
|
|
double gx, gy;
|
|
@@ -245,7 +343,8 @@ bool remove_outliers_c(const studio_line_c *line, double distance, studio_line_c
|
|
|
double dist_square = distance * distance;
|
|
|
// 经纬度转高斯克吕格投影
|
|
|
studio_line_c gauss_line = studio_line_c_init();
|
|
|
- double central = (int) (line->data[0].x) / 3 * 3;
|
|
|
+ // double central = ((int) (line->data[0].x) / 3) * 3;
|
|
|
+ double central = floor((line->data[0].x + 1.5) / 3.0) * 3.0;
|
|
|
for (unsigned int i = 0; i < line->size; i++)
|
|
|
{
|
|
|
double gx, gy;
|
|
@@ -264,7 +363,7 @@ bool remove_outliers_c(const studio_line_c *line, double distance, studio_line_c
|
|
|
studio_line_c_add_point(&gs_filtered_line, gauss_line.data[0]); // 保留第一个点
|
|
|
for (unsigned int i = 1; i < gauss_line.size; i++) // 遍历相邻点对
|
|
|
{
|
|
|
- double d_sq = distance_squared(&gauss_line.data[i - 1], &gauss_line.data[i]);
|
|
|
+ double d_sq = distance_squared_c(&gauss_line.data[i - 1], &gauss_line.data[i]);
|
|
|
if (d_sq <= dist_square)
|
|
|
{
|
|
|
studio_line_c_add_point(&gs_filtered_line, gauss_line.data[i]);
|