algorithm/line_vacuate/studio_geo_algo_c.c

/**
 ******************************************************************************
 * @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 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;
}