/** ****************************************************************************** * @file : studio_geo_algo_c.c * @author : wangyingjie * @brief : None * @attention : None * @date : 2025/5/11 ****************************************************************************** */ #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) { double A = p->x - start->x; double B = p->y - start->y; double C = end->x - start->x; double D = end->y - start->y; double dot = A * C + B * D; double len_sq = C * C + D * D; double param = len_sq == 0 ? -1 : dot / len_sq; double xx, yy; if (param < 0) { xx = start->x; yy = start->y; } else if (param > 1) { xx = end->x; yy = end->y; } else { xx = start->x + param * C; yy = start->y + param * D; } double dist_square = (xx - p->x) * (xx - p->x) + (yy - p->y) * (yy - p->y); return dist_square; } // 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) { // 计算最大距离 double max_dist = 0; size_t index = start; double epsilon_square = epsilon * epsilon; // 遍历从start+1到end-1的点 for (size_t i = start + 1; i < end; ++i) { // 计算点到线段的距离的平方 double dist_square = point_line_dist_square_c(&points->data[i], &points->data[start], &points->data[end]); // 如果距离大于最大距离,则更新最大距离和索引 if (dist_square > max_dist) { index = i; max_dist = dist_square; } } // 如果最大距离大于epsilon的平方,并且索引不是start和end,则递归调用Douglas-Peucker算法 if (max_dist > epsilon_square && index != start && index != end) { douglas_peucker_c(points, start, index, epsilon, indices, index_count); douglas_peucker_c(points, index, end, epsilon, indices, index_count); } // 否则,将start和end的索引添加到indices数组中 else { indices[*index_count] = start; (*index_count)++; if (start != end) { indices[*index_count] = end; (*index_count)++; } } } unsigned int hash(unsigned int value, unsigned int size) { return value % size; } void remove_duplicates(unsigned int **indices, unsigned int *size) { if (*size == 0) { return; } 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)); for (unsigned int i = 0; i < hash_table_size; i++) { hash_table[i] = false; // 初始化哈希表 } for (unsigned int i = 0; i < *size; i++) { unsigned int current_value = (*indices)[i]; unsigned int hash_index = hash(current_value, hash_table_size); // 检查该元素是否已经在哈希表中 if (!hash_table[hash_index]) { hash_table[hash_index] = true; unique_indices[unique_count++] = current_value; } } // 释放原有的 indices 和哈希表 free(*indices); free(hash_table); // 更新 indices 指针和大小 *indices = unique_indices; *size = unique_count; } bool line_vacuate_c(const studio_line_c *line, const int max_points, const double epsilon, studio_line_c *vacuate_line) { // 初始化状态变量 bool status = false; // 如果输入的线段为空,则打印错误信息并返回false if (line->size == 0) { printf("line is empty\n"); return false; } // 如果输入的线段点数小于等于最大点数,则直接将线段点添加到vacuate_line中,并返回true if (line->size <= max_points) { for (size_t i = 0; i < line->size; ++i) { studio_line_c_add_point(vacuate_line, line->data[i]); } status = true; return status; } // 分配内存空间存储线段点的索引 unsigned int *indices = (unsigned int *) malloc(line->size * sizeof(unsigned int)); unsigned int index_count = 0; double t_epsilon = epsilon; // 循环执行Douglas-Peucker算法,直到索引点数小于等于最大点数 while (1) { index_count = 0; douglas_peucker_c(line, 0, line->size - 1, t_epsilon, indices, &index_count); // 将indices 中的数据去重 remove_duplicates(&indices, &index_count); if (index_count <= (unsigned int) max_points) { break; } t_epsilon *= 1.1; } // 将索引点添加到vacuate_line中 for (unsigned int i = 0; i < index_count; ++i) { studio_line_c_add_point(vacuate_line, line->data[indices[i]]); } // 释放内存空间 free(indices); status = true; return status; } bool line_vacuate_gcs_c(const studio_line_c *line, const int max_points, const double epsilon, studio_line_c *vacuate_line) { bool status = false; if (!line || !vacuate_line || line->size <= 1) // 添加空指针和长度校验 { printf("line or vacuate_line is null or line size <= 1"); return status; } // 经纬度转高斯克吕格投影 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); } studio_line_c vac_gauss_line = studio_line_c_init(); bool res = line_vacuate_c(&gauss_line, max_points, epsilon, &vac_gauss_line); if (!res) { printf("Failed to simplify line."); return status; } studio_line_c_destroy(&gauss_line); // 高斯克吕格转经纬度 for (unsigned int i = 0; i < vac_gauss_line.size; i++) { double lon, lat; gauss_to_lonlat(central, vac_gauss_line.data[i].x, vac_gauss_line.data[i].y, &lon, &lat); studio_point_c point = studio_point_init(lon, lat); studio_line_c_add_point(vacuate_line, point); } studio_line_c_destroy(&vac_gauss_line); 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; }