添加 移除离群点 算法

coordinate_convert/c/studio_proj_c.c

@@ -75,17 +75,17 @@ void gauss_to_lonlat(const double central, const double gx, const double gy, dou
 
 void lonlat_to_mercator(const double lon, const double lat, double *mctx, double *mcty)
 {
-    // 经度转换：phi = lon * (M_PI / 180.0)
+    // 经度转换：phi = lon * (AO_M_PI / 180.0)
     double phi = DEG2RAD_C(lon);
 
-    // 纬度转换：lamb = lat * (M_PI / 180.0)
+    // 纬度转换：lamb = lat * (AO_M_PI / 180.0)
     double lamb = DEG2RAD_C(lat);
 
     // 墨卡托投影公式：mctx = WGS84_A_C * phi
     *mctx = WGS84_A_C * phi;
 
-    // 墨卡托投影公式：mcty = WGS84_A_C * log(tan(M_PI / 4 + lamb / 2))
-    *mcty = WGS84_A_C * log(tan(M_PI / 4 + lamb / 2));
+    // 墨卡托投影公式：mcty = WGS84_A_C * log(tan(AO_M_PI / 4 + lamb / 2))
+    *mcty = WGS84_A_C * log(tan(AO_M_PI / 4 + lamb / 2));
 }
 
 void mercator_to_lonlat(const double mctx, const double mcty, double *lon, double *lat)
@@ -94,7 +94,7 @@ void mercator_to_lonlat(const double mctx, const double mcty, double *lon, doubl
     *lon = RAD2DEG_C(mctx / WGS84_A_C);
 
     // 纬度转换：使用逆墨卡托公式
-    *lat = RAD2DEG_C(2 * atan(exp(mcty / WGS84_A_C)) - M_PI / 2);
+    *lat = RAD2DEG_C(2 * atan(exp(mcty / WGS84_A_C)) - AO_M_PI / 2);
 }
 
 // 将高斯投影坐标转换为墨卡托投影坐标

coordinate_convert/c/studio_proj_c.h

@@ -19,7 +19,6 @@
 #define RAD2DEG_C(rad) ((rad) * 180.0 / AO_M_PI)
 
 #define WGS84_A_C 6378137.0                                     // 长半轴
-#define WGS84_E2_C 0.0066943799901413165                        // 第一偏心扁率的平方
 #define INVF_C 298.257223563                                    // 扁率的倒数
 #define WGS84_B_C 6356752.3142                                  // 短半轴
 #define WGS84_F_C 1 / INVF_C                                    // 扁率

line_vacuate/studio_geo_algo_c.c

@@ -9,6 +9,15 @@
  */
 
 #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)
+{
+    const double dx = p1->x - p2->x;
+    const double dy = p1->y - p2->y;
+    return dx * dx + dy * dy;
+}
 
 double point_line_dist_square_c(const studio_point_c *p, const studio_point_c *start, const studio_point_c *end)
 {
@@ -27,13 +36,11 @@ double point_line_dist_square_c(const studio_point_c *p, const studio_point_c *s
     {
         xx = start->x;
         yy = start->y;
-    }
-    else if (param > 1)
+    } else if (param > 1)
     {
         xx = end->x;
         yy = end->y;
-    }
-    else
+    } else
     {
         xx = start->x + param * C;
         yy = start->y + param * D;
@@ -44,7 +51,8 @@ double point_line_dist_square_c(const studio_point_c *p, const studio_point_c *s
 }
 
 // Douglas-Peucker算法，用于简化多边形或曲线
-void douglas_peucker_c(const studio_line_c *points, size_t start, size_t end, double epsilon, unsigned int *indices, unsigned int *index_count)
+void douglas_peucker_c(const studio_line_c *points, size_t start, size_t end, double epsilon, unsigned int *indices,
+                       unsigned int *index_count)
 {
     // 计算最大距离
     double max_dist = 0;
@@ -96,12 +104,12 @@ void remove_duplicates(unsigned int **indices, unsigned int *size)
         return;
     }
 
-    unsigned int *unique_indices = (unsigned int *)malloc(*size * sizeof(unsigned int));
+    unsigned int *unique_indices = (unsigned int *) malloc(*size * sizeof(unsigned int));
     unsigned int unique_count = 0;
 
     // 创建一个哈希表
     unsigned int hash_table_size = *size * 4; // 哈希表大小可以稍微增大以减少冲突 这个可能有问题
-    bool *hash_table = (bool *)malloc(hash_table_size * sizeof(bool));
+    bool *hash_table = (bool *) malloc(hash_table_size * sizeof(bool));
     for (unsigned int i = 0; i < hash_table_size; i++)
     {
         hash_table[i] = false; // 初始化哈希表
@@ -154,7 +162,7 @@ bool line_vacuate_c(const studio_line_c *line, const int max_points, const doubl
     }
 
     // 分配内存空间存储线段点的索引
-    unsigned int *indices = (unsigned int *)malloc(line->size * sizeof(unsigned int));
+    unsigned int *indices = (unsigned int *) malloc(line->size * sizeof(unsigned int));
     unsigned int index_count = 0;
 
     double t_epsilon = epsilon;
@@ -167,7 +175,7 @@ bool line_vacuate_c(const studio_line_c *line, const int max_points, const doubl
         // 将indices 中的数据去重
         remove_duplicates(&indices, &index_count);
 
-        if (index_count <= (unsigned int)max_points)
+        if (index_count <= (unsigned int) max_points)
         {
             break;
         }
@@ -184,4 +192,62 @@ bool line_vacuate_c(const studio_line_c *line, const int max_points, const doubl
     free(indices);
     status = true;
     return status;
-}
+}
+
+
+bool remove_outliers_c(const studio_line_c *line, double distance, studio_line_c *outliers_line)
+{
+    bool status = false;
+    if (!line || !outliers_line || line->size <= 1) // 添加空指针和长度校验
+    {
+        printf("line or outliers_line is null or line size <= 1");
+        return status;
+    }
+    double dist_square = distance * distance;
+    // 经纬度转高斯克吕格投影
+    studio_line_c gauss_line = studio_line_c_init();
+    double central = (int) (line->data[0].x) / 3 * 3;
+    for (unsigned int i = 0; i < line->size; i++)
+    {
+        double gx, gy;
+        lonlat_to_gauss(central, line->data[i].x, line->data[i].y, &gx, &gy);
+        studio_point_c point = studio_point_init(gx, gy);
+        studio_line_c_add_point(&gauss_line, point);
+    }
+
+    // 离群点检测 ------------------------------------------------
+    if (gauss_line.size<= 1)
+    {
+        printf("gauss_line size <= 1");
+        return status;
+    }
+    studio_line_c gs_filtered_line = studio_line_c_init(); // 新建过滤后的线段
+    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]);
+        if (d_sq <= dist_square)
+        {
+            studio_line_c_add_point(&gs_filtered_line, gauss_line.data[i]);
+        }
+        else
+        {
+            int a=0;
+        }
+    }
+    studio_line_c_destroy(&gauss_line); // 释放原始数据的高斯投影
+
+    // 高斯克吕格转经纬度
+    for (unsigned int i = 0; i < gs_filtered_line.size; i++)
+    {
+        double lon, lat;
+        gauss_to_lonlat(central, gs_filtered_line.data[i].x, gs_filtered_line.data[i].y, &lon, &lat);
+        studio_point_c point = studio_point_init(lon, lat);
+        studio_line_c_add_point(outliers_line, point);
+    }
+
+    studio_line_c_destroy(&gs_filtered_line); // 释放移除离群点后的高斯投影
+
+    status = true;
+    return status;
+}

line_vacuate/studio_geo_algo_c.h

@@ -15,6 +15,12 @@
 
 /////////////////// 矢量(线)抽稀算法 ///////////////////
 
+/// 计算两点之间的距离的平方
+/// @param p1
+/// @param p2
+/// @return
+double distance_squared(const studio_point_c *p1, const studio_point_c *p2);
+
 /// 计算点到线段的距离的平方
 /// \param p 点
 /// \param start 线段起点
@@ -29,7 +35,8 @@ double point_line_dist_square_c(const studio_point_c *p, const studio_point_c *s
 /// \param epsilon 抽稀容差,若某点到当前线段的最大垂直距离超过epsilon，则保留该点,单位与适量坐标系的单位一致
 /// \param indices 存储简化后的点索引
 /// \param index_count 存储简化后的点索引数量
-void douglas_peucker_c(const studio_line_c *points, size_t start, size_t end, double epsilon, unsigned int *indices, unsigned int *index_count);
+void douglas_peucker_c(const studio_line_c *points, size_t start, size_t end, double epsilon, unsigned int *indices,
+                       unsigned int *index_count);
 
 /// 矢量线段抽稀算法
 /// \param line 线段
@@ -50,4 +57,16 @@ unsigned int hash(unsigned int value, unsigned int size);
 /// \param size 点索引数量
 void remove_duplicates(unsigned int **indices, unsigned int *size);
 
+
+/////////////////// 移除线段中的离群点 ///////////////////
+
+
+/// 移除线段中的离群点
+/// @param line 原始线段
+/// @param distance 离群点距离
+/// @param outliers_line 滤波后的线段
+/// @return
+bool remove_outliers_c(const studio_line_c *line, double distance, studio_line_c *outliers_line);
+
+
 #endif  // STUDIO_GEO_ALGO_C_H

line_vacuate/studio_geo_c.c

@@ -31,9 +31,9 @@ studio_line_c studio_line_c_init()
 {
     studio_line_c line;
     line.size = 0;
-    line.capacity = 0;
-    line.data = (studio_point_c *)malloc(line.capacity * sizeof(studio_point_c));
-    return line;  // 注意: 返回结构体副本 (需配合 C99 或更高标准)
+    line.capacity = 4;
+    line.data = (studio_point_c *) malloc(line.capacity * sizeof(studio_point_c));
+    return line; // 注意: 返回结构体副本 (需配合 C99 或更高标准)
 }
 
 // 销毁线段 (释放内存)
@@ -55,7 +55,7 @@ void studio_line_c_add_point(studio_line_c *line, studio_point_c point)
     if (line->size >= line->capacity)
     {
         line->capacity += 10;
-        studio_point_c *new_data = (studio_point_c *)realloc(line->data, line->capacity * sizeof(studio_point_c));
+        studio_point_c *new_data = (studio_point_c *) realloc(line->data, line->capacity * sizeof(studio_point_c));
         if (!new_data)
         {
             // 此处可添加错误处理 (例如: 终止程序或返回错误码)
@@ -81,7 +81,7 @@ const studio_point_c *studio_line_c_get_point(const studio_line_c *line, unsigne
     {
         return &(line->data[index]);
     }
-    return NULL;  // 越界返回空指针
+    return NULL; // 越界返回空指针
 }
 
 ////////////// 矩形 //////////////
@@ -110,7 +110,8 @@ void studio_rect_correct(studio_rect_c *rect)
 
 bool studio_rect_intersect(const studio_rect_c *r1, const studio_rect_c *r2)
 {
-    return !(r1->right_bottom.x < r2->left_top.x || r1->left_top.x > r2->right_bottom.x || r1->right_bottom.y > r2->left_top.y || r1->left_top.y < r2->right_bottom.y);
+    return !(r1->right_bottom.x < r2->left_top.x || r1->left_top.x > r2->right_bottom.x || r1->right_bottom.y > r2->
+             left_top.y || r1->left_top.y < r2->right_bottom.y);
 }
 
 ////////////// 圆 //////////////
@@ -136,7 +137,8 @@ studio_triangle_c studio_triangle_init(studio_point_c a, studio_point_c b, studi
 
 double studio_triangle_oriented_area(const studio_triangle_c *triangle)
 {
-    return (triangle->a.x * (triangle->b.y - triangle->c.y) + triangle->b.x * (triangle->c.y - triangle->a.y) + triangle->c.x * (triangle->a.y - triangle->b.y)) / 2.0;
+    return (triangle->a.x * (triangle->b.y - triangle->c.y) + triangle->b.x * (triangle->c.y - triangle->a.y) + triangle
+            ->c.x * (triangle->a.y - triangle->b.y)) / 2.0;
 }
 
 double studio_triangle_area(const studio_triangle_c *triangle)
@@ -161,7 +163,8 @@ double studio_ellipse_circumference(const studio_ellipse_c *ellipse, int Ramanuj
 {
     if (Ramanujan == 1)
     {
-        return AO_M_PI * (3 * (ellipse->rx + ellipse->ry) - sqrt((3 * ellipse->rx + ellipse->ry) * (ellipse->rx + 3 * ellipse->ry)));
+        return AO_M_PI * (3 * (ellipse->rx + ellipse->ry) - sqrt(
+                              (3 * ellipse->rx + ellipse->ry) * (ellipse->rx + 3 * ellipse->ry)));
     }
     return 2 * AO_M_PI * sqrt((ellipse->rx * ellipse->rx + ellipse->ry * ellipse->ry) / 2);
-}
+}

line_vacuate/studio_geo_c.h

@@ -45,11 +45,13 @@ bool studio_point_equal(const studio_point_c *p1, const studio_point_c *p2);
 /// 线
 /****************************************/
 
+// 该变量类型的定义为： studio_line_c temp_line = studio_line_c_init();
+// 并且 需要调用 studio_line_c_destroy(&temp_line); 释放内存
 typedef struct
 {
-    studio_point_c *data;   // 点数据指针
-    unsigned int size;      // 当前元素个数
-    unsigned int capacity;  // 当前分配的内存容量
+    studio_point_c *data; // 点数据指针
+    unsigned int size; // 当前元素个数
+    unsigned int capacity; // 当前分配的内存容量
 } studio_line_c;
 
 // 初始化线段 (默认容量为4)

line_vacuate/task_algo_c.cpp

@@ -26,8 +26,8 @@ extern "C"
 
 void studio_line_to_c_line(const studio_line &line, studio_line_c *c_line)
 {
-//    c_line->size = line.size();
-    c_line->data = (studio_point_c *)malloc(line.size() * sizeof(studio_point_c));
+    //    c_line->size = line.size();
+    c_line->data = (studio_point_c *) malloc(line.size() * sizeof(studio_point_c));
     for (size_t i = 0; i < line.size(); ++i)
     {
         studio_point_c p = studio_point_init(line[i].x, line[i].y);
@@ -44,6 +44,7 @@ void c_line_to_studio_line(const studio_line_c *c_line, studio_line &line)
     }
 }
 
+///////////////////// 移除离群点 ////////////////////////////////
 int main()
 {
     printf("\n\n===================== %s =====================\n\n", __FILE__);
@@ -59,63 +60,117 @@ int main()
     std::vector<studio_geo_coll> res_collections;
     std::vector<studio_geo_coll> collections;
     silly::geo::utils::read_geo_coll(path, collections);
-    for (auto &coll : collections)
+
+    for (auto &coll: collections)
     {
-        // ------------- 转换为高斯投影 -------------
-        studio_line gauss_line;
-        double central = static_cast<int>(coll.m_line[0].x / 3) * 3;
-        for (auto &point : coll.m_line)
-        {
-            double gx = 0.0;
-            double gy = 0.0;
-            lonlat_to_gauss(central, point.x, point.y, &gx, &gy);
-            gauss_line.push_back(studio_point(gx, gy));
-        }
-        studio_line_c gauss_line_c = studio_line_c_init();
-        studio_line_c vac_gauss_line_c = studio_line_c_init();
-        studio_line_to_c_line(gauss_line, &gauss_line_c);
-
-        // 简化线段，目标点数为28个
-        int max_points = 28;
-        double epsilon = 1.0;
-
-        bool res = line_vacuate_c(&gauss_line_c, max_points, epsilon, &vac_gauss_line_c);
-        if (!res)
-        {
-            std::cout << "Failed to simplify line." << std::endl;
-            return 1;
-        }
-
-        // 高斯投影在转回经纬度
+        studio_line_c line_c = studio_line_c_init();
+        studio_line_to_c_line(coll.m_line, &line_c);
+        std::cout << "原始数据量：" << coll.m_line.size() << std::endl;
+
+        studio_line_c outliers_line = studio_line_c_init();
+        remove_outliers_c(&line_c, 1, &outliers_line);
+
         studio_line simplified_line;
-        for (int i = 0; i < vac_gauss_line_c.size; i++)
-        {
-            double lon = 0.0;
-            double lat = 0.0;
-            const studio_point_c *p = studio_line_c_get_point(&vac_gauss_line_c, i);
-            gauss_to_lonlat(central, p->x, p->y, &lon, &lat);
-            simplified_line.push_back(studio_point(lon, lat));
-        }
+        c_line_to_studio_line(&outliers_line, simplified_line);
+        std::cout << "简化后数据量：" << simplified_line.size() << std::endl;
 
         studio_geo_coll temp;
         temp.m_type = enum_geometry_type::egtLineString;
         temp.m_line = simplified_line;
         res_collections.push_back(temp);
 
-        studio_line_c_destroy(&gauss_line_c);
-        studio_line_c_destroy(&vac_gauss_line_c);
+        studio_line_c_destroy(&line_c);
+        studio_line_c_destroy(&outliers_line);
 
         break;
     }
+
+
     std::string output_path;
 #ifdef IS_WINDOWS
-    output_path = "D:/5_file/2_readfile/geojson/multi_point/fitting_examples_res_1.geojson";
+    output_path = "D:/5_file/2_readfile/geojson/multi_point/remove_outliers.geojson";
 #else
-    output_path = "/home/wyj/myself/2_data/2_geojson/multi_point/fitting_examples_res_1_c.geojson";
+    output_path = "/home/wyj/myself/2_data/2_geojson/multi_point/remove_outliers.geojson";
 #endif
     silly::geo::utils::write_geo_coll(output_path, res_collections);
-
     silly::geo::utils::destroy_gdal_env();
-
     return 0;
 }
+
+
+// //////////////////// 路径压缩算法测试 //////////////////////
+// int main()
+// {
+//     printf("\n\n===================== %s =====================\n\n", __FILE__);
+//     silly::geo::utils::init_gdal_env();
+//     std::string path;
+// #ifdef IS_WINDOWS
+//     path = "D:/5_file/2_readfile/geojson/multi_point/fitting_examples.geojson";
+// #else
+//     path = "/home/wyj/myself/2_data/2_geojson/multi_point/fitting_examples.geojson";
+// #endif
+//
+//     std::cout << "path: " << path << std::endl;
+//     std::vector<studio_geo_coll> res_collections;
+//     std::vector<studio_geo_coll> collections;
+//     silly::geo::utils::read_geo_coll(path, collections);
+//     for (auto &coll : collections)
+//     {
+//         // ------------- 转换为高斯投影 -------------
+//         studio_line gauss_line;
+//         double central = static_cast<int>(coll.m_line[0].x / 3) * 3;
+//         for (auto &point : coll.m_line)
+//         {
+//             double gx = 0.0;
+//             double gy = 0.0;
+//             lonlat_to_gauss(central, point.x, point.y, &gx, &gy);
+//             gauss_line.push_back(studio_point(gx, gy));
+//         }
+//         studio_line_c gauss_line_c = studio_line_c_init();
+//         studio_line_c vac_gauss_line_c = studio_line_c_init();
+//         studio_line_to_c_line(gauss_line, &gauss_line_c);
+//
+//         // 简化线段，目标点数为28个
+//         int max_points = 28;
+//         double epsilon = 1.0;
+//
+//         bool res = line_vacuate_c(&gauss_line_c, max_points, epsilon, &vac_gauss_line_c);
+//         if (!res)
+//         {
+//             std::cout << "Failed to simplify line." << std::endl;
+//             return 1;
+//         }
+//
+//         // 高斯投影在转回经纬度
+//         studio_line simplified_line;
+//         for (int i = 0; i < vac_gauss_line_c.size; i++)
+//         {
+//             double lon = 0.0;
+//             double lat = 0.0;
+//             const studio_point_c *p = studio_line_c_get_point(&vac_gauss_line_c, i);
+//             gauss_to_lonlat(central, p->x, p->y, &lon, &lat);
+//             simplified_line.push_back(studio_point(lon, lat));
+//         }
+//
+//         studio_geo_coll temp;
+//         temp.m_type = enum_geometry_type::egtLineString;
+//         temp.m_line = simplified_line;
+//         res_collections.push_back(temp);
+//
+//         studio_line_c_destroy(&gauss_line_c);
+//         studio_line_c_destroy(&vac_gauss_line_c);
+//
+//         break;
+//     }
+//     std::string output_path;
+// #ifdef IS_WINDOWS
+//     output_path = "D:/5_file/2_readfile/geojson/multi_point/fitting_examples_res_1.geojson";
+// #else
+//     output_path = "/home/wyj/myself/2_data/2_geojson/multi_point/fitting_examples_res_1_c.geojson";
+// #endif
+//     silly::geo::utils::write_geo_coll(output_path, res_collections);
+//
+//     silly::geo::utils::destroy_gdal_env();
+//
+//     return 0;
+// }