{
    "configurations": [
        {
            "name": "Linux",
            "includePath": [
                "${workspaceFolder}/**",
                "/usr/include",
                "/usr/local/include",
                "/usr/include/eigen3",
                "/usr/include/python3.12"
            ],
            "defines": [],
            "compilerPath": "/usr/bin/gcc",
            "cStandard": "c17",
            "cppStandard": "c++17",
            "intelliSenseMode": "linux-gcc-arm64"
        }
    ],
    "version": 4
}

.
|-- CMakeLists.txt

cmake_minimum_required(VERSION 3.12)

project(example VERSION 1.0.0 LANGUAGES C CXX)

find_package(GEOS 3.12 REQUIRED)
find_package(Python 3.12 COMPONENTS Interpreter Development REQUIRED)
find_package(pybind11 CONFIG REQUIRED)
find_package (GeographicLib REQUIRED)
find_package(Eigen3 REQUIRED NO_MODULE)

message("GeographicLib_INCLUDE_DIRS:" ${GeographicLib_INCLUDE_DIRS})
message("GeographicLib_LIBRARIES:" ${GeographicLib_LIBRARIES})
message("GEOS_LIBRARIES:" ${GEOS_LIBRARIES})
message("GEOS_INCLUDE_DIR:" ${GEOS_INCLUDE_DIR})

file(GLOB SRC_LIST1 "*.cpp")
file(GLOB SRC_LIST2 "src/*.cpp")

include_directories(${GEOS_INCLUDE_DIR}) 
include_directories(${GeographicLib_INCLUDE_DIRS})

include_directories(
    /usr/local/include
    /usr/include/eigen3
    ${CMAKE_CURRENT_SOURCE_DIR}/include
)

link_directories(
    /usr/local/lib
)

pybind11_add_module(example ${SRC_LIST1} ${SRC_LIST2})


target_link_libraries(example 
    PRIVATE GEOS::geos
    ${GeographicLib_LIBRARIES}
    )

add_subdirectory(app)

|-- apt.install

wget https://github.com/geographiclib/geographiclib/archive/refs/tags/v2.4.tar.gz
cd geographiclib
mkdir build && cd build
cmake .. #无报错
make #无报错
make install #无报错

===================================================================================================================

apt install libeigen3-dev

===================================================================================================================


===================================================================================================================

|-- build
|-- app
| |-- CMakeLists.txt

# cmake_minimum_required(VERSION 3.12)

# project(example VERSION 1.0.0 LANGUAGES C CXX)

# find_package(GEOS 3.12 REQUIRED)
# find_package(pybind11 REQUIRED)
# find_package (GeographicLib REQUIRED)
# find_package(Eigen3 REQUIRED NO_MODULE)

# message("GeographicLib_INCLUDE_DIRS:" ${GeographicLib_INCLUDE_DIRS})
# message("GeographicLib_LIBRARIES:" ${GeographicLib_LIBRARIES})
# message("GEOS_LIBRARIES:" ${GEOS_LIBRARIES})
# message("GEOS_INCLUDE_DIR:" ${GEOS_INCLUDE_DIR})


# include_directories(${GEOS_INCLUDE_DIR}) 
# include_directories(${GeographicLib_INCLUDE_DIRS})

# include_directories(
#     /usr/local/include
#     /usr/include/eigen3
#     ${CMAKE_CURRENT_SOURCE_DIR}/../include
# )

# link_directories(
#     /usr/local/lib
# )

# # 将链接模式设为静态
# # set(CMAKE_EXE_LINKER_FLAGS "-static")
# # set(CMAKE_EXE_LINKER_FLAGS "-static -static-libgcc -static-libstdc++")


add_executable(example_exe main.cpp)

target_link_libraries(example_exe 
    PRIVATE GEOS::geos
    ${GeographicLib_LIBRARIES}
)

| |-- build
| `-- main.cpp

#include <iostream>
#include <geoutil.hpp>

int main()
{
    Eigen::Vector2d ba(1, 1);
    Eigen::Vector2d bc(1,0);
    float angle = GeoUtil::getVectorAngle(ba, bc);
    std::cout << angle << std::endl;
}

|-- include
| |-- a.h

#ifndef A_H_
#define A_H_

#include <iostream>
#include <string>

int add(int i, int j);
std::string lineInterpolatePoint(const std::string wkt, float index);
std::string getIntersect();

#endif

| |-- duckdb_udf.h

#pragma once

#include <string>

namespace DuckdbUDF
{
    float geomZDiff(const std::string &geomWkt);
    float lineStringLength2d(const std::string &geomWkt);
    float line3PointMinAngle(const std::string &geomWkt);
}

| `-- geoutil.hpp

#pragma once
#include <iostream>
#include <string>
#include <vector>
#include <geos/io/WKTReader.h>
#include <GeographicLib/Geodesic.hpp>
#include <Eigen/Dense>

namespace GeoUtil
{
    std::unique_ptr<geos::geom::LineString> getLineStringFromWkt(const std::string &wkt)
    {
        geos::io::WKTReader reader;
        std::unique_ptr<geos::geom::Geometry> geom = nullptr;
        try
        {
            geom = reader.read(wkt);
        }
        catch (const geos::util::GEOSException &e)
        {
            std::cerr << "Failed to read WKT geometry: " << e.what() << std::endl;
            return nullptr;
        }

        if (!geom || geom->getGeometryTypeId() != geos::geom::GEOS_LINESTRING)
        {
            std::cerr << "The WKT geometry is not a LineString." << std::endl;
            return nullptr;
        }

        return std::unique_ptr<geos::geom::LineString>(dynamic_cast<geos::geom::LineString *>(geom.release()));
    }
    std::vector<geos::geom::Coordinate> getCoordsFromWkt(const std::string &wkt)
    {
        std::vector<geos::geom::Coordinate> coords;
        geos::io::WKTReader reader;
        std::unique_ptr<geos::geom::Geometry> geom = nullptr;
        try
        {
            geom = reader.read(wkt);
        }
        catch (const geos::util::GEOSException &e)
        {
            std::cerr << "Failed to read WKT geometry: " << e.what() << std::endl;
            return coords;
        }
        if (geom)
        {
            std::unique_ptr<geos::geom::CoordinateSequence> coordSeq;
            if (geom->getGeometryTypeId() == geos::geom::GEOS_POINT)
            {
                std::unique_ptr<geos::geom::Point> point(dynamic_cast<geos::geom::Point *>(geom.release()));
                coordSeq = point->getCoordinates();
            }
            else if (geom->getGeometryTypeId() == geos::geom::GEOS_LINESTRING)
            {
                std::unique_ptr<geos::geom::LineString> lineString(dynamic_cast<geos::geom::LineString *>(geom.release()));
                coordSeq = lineString->getCoordinates();
            }
            else if (geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON)
            {
                std::unique_ptr<geos::geom::Polygon> polygon(dynamic_cast<geos::geom::Polygon *>(geom.release()));
                coordSeq = polygon->getExteriorRing()->getCoordinates();
            }
            for (std::size_t i = 0; i < coordSeq->getSize(); ++i)
            {
                coords.push_back(coordSeq->getAt(i));
            }
        }
        return coords;
    }
    float getDistance2d(float lon1, float lat1, float lon2, float lat2)
    {
        const GeographicLib::Geodesic &geod = GeographicLib::Geodesic::WGS84();
        double distance;
        geod.Inverse(lat1, lon1, lat2, lon2, distance);
        return distance;
    }

    double getVectorAngle(const Eigen::Vector2d &ba, const Eigen::Vector2d &bc)
    {
        // 计算点积
        double dotVal = ba.dot(bc);

        // 计算向量的范数
        double normBa = ba.norm();
        double normBc = bc.norm();

        // 计算 cos(theta)
        double cosTheta = dotVal / (normBa * normBc);

        // 将 cos(theta) 限制在 [-1, 1] 范围内
        cosTheta = std::max(-1.0, std::min(1.0, cosTheta));

        // 计算角度（弧度）
        double angle = std::acos(cosTheta);

        // 将角度转换为度
        return angle * (180.0 / M_PI);
    }

}

|-- pybind.cpp

#include <pybind11/pybind11.h>
#include "a.h"
#include "duckdb_udf.h"

PYBIND11_MODULE(example, m)
{
      m.doc() = "pybind11 example plugin"; // optional module docstring

      m.def("add", &add, "A function which adds two numbers");
      m.def("add_lambda", [](int a, int b)
            { return a + b; }, "A function which adds two numbers");
      m.def("getIntersect", &getIntersect, "getIntersect");
      m.def("lineInterpolatePoint", &lineInterpolatePoint, "lineInterpolatePoint",
            pybind11::arg("line_wkt"), pybind11::arg("index") = 0.6);
      m.def("geomZDiff", &DuckdbUDF::geomZDiff, "geom_z_diff");
      m.def("lineStringLength2d", &DuckdbUDF::lineStringLength2d, "lineStringLength2d");
      m.def("line3PointMinAngle", &DuckdbUDF::line3PointMinAngle, "line3PointMinAngle");
}

`-- src
|-- a.cpp

#include <iostream>
#include <geos/io/WKTReader.h>
#include <geos/io/WKTWriter.h>
#include <geos/linearref/LengthIndexOfPoint.h>
#include <geos/linearref/LengthIndexedLine.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Coordinate.h>
#include <string>
#include "a.h"

int add(int i, int j)
{
    return i + j;
}

std::string getIntersect()
{
    geos::geom::GeometryFactory::Ptr factory = geos::geom::GeometryFactory::create();
    geos::io::WKTReader reader(*factory);

    std::string wkt_a("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))");
    std::string wkt_b("POLYGON((5 5, 15 5, 15 15, 5 15, 5 5))");

    std::unique_ptr<geos::geom::Geometry> geom_a(reader.read(wkt_a));
    std::unique_ptr<geos::geom::Geometry> geom_b(reader.read(wkt_b));

    std::unique_ptr<geos::geom::Geometry> inter = geom_a->intersection(geom_b.get());

    geos::io::WKTWriter writer;
    writer.setTrim(true);
    std::string inter_wkt = writer.write(inter.get());
    return inter_wkt;
}
std::string lineInterpolatePoint(const std::string line_wkt, float index)
{
    geos::io::WKTReader wktReader;
    std::unique_ptr<geos::geom::Geometry> geom(wktReader.read(line_wkt));
    if (!geom)
    {
        std::cout << "wkt有错误" << std::endl;
        return "";
    }
    std::string geomType = geom->getGeometryType();
    if (geomType != "LineString")
    {
        std::cout << "wkt需要是linestring类型" << std::endl;
        return "";
    }
    std::unique_ptr<geos::geom::LineString> line(static_cast<geos::geom::LineString *>(geom.release()));
    geos::linearref::LengthIndexedLine lengthIndexedLine(line.release());
    geos::geom::Coordinate c = lengthIndexedLine.extractPoint(0.5);
    geos::geom::GeometryFactory::Ptr geomFactory = geos::geom::GeometryFactory::create();
    std::unique_ptr<geos::geom::Point> point = geomFactory->createPoint(c);
    // std::cout << point->getX() << '\t' << point->getY() << std::endl;
    geos::io::WKTWriter wktWriter;
    std::string pointWkt = wktWriter.write(point.get());
    return pointWkt;
}

`-- duckdb_udf.cpp

#include <iostream>
#include <string>
#include <geos/io/WKTReader.h>
#include "duckdb_udf.h"
#include "geoutil.hpp"

float DuckdbUDF::geomZDiff(const std::string &geomWkt)
{
    try
    {
        std::unique_ptr<geos::geom::LineString> lineString = GeoUtil::getLineStringFromWkt(geomWkt);
        if (lineString)
        {
            if (lineString->hasZ())
            {
                std::unique_ptr<geos::geom::CoordinateSequence> coordSeq = lineString->getCoordinates();
                const geos::geom::Coordinate &firstCoord = coordSeq->getAt(0);
                double maxZ = firstCoord.z;
                double minZ = firstCoord.z;
                for (std::size_t i = 1; i < coordSeq->getSize(); ++i)
                {
                    const geos::geom::Coordinate &coord = coordSeq->getAt(i);
                    if (coord.z > maxZ)
                    {
                        maxZ = coord.z;
                    }
                    if (coord.z < minZ)
                    {
                        minZ = coord.z;
                    }
                }
                return maxZ - minZ;
            }
        }
    }
    catch (const geos::util::GEOSException &e)
    {
        std::cerr << "GEOS Exception: " << e.what() << std::endl;
    }
    catch (const std::exception &e)
    {
        std::cerr << "Standard Exception: " << e.what() << std::endl;
    }

    return 0;
}

float DuckdbUDF::lineStringLength2d(const std::string &geomWkt)
{
    try
    {
        std::vector<geos::geom::Coordinate> coords = GeoUtil::getCoordsFromWkt(geomWkt);
        if (coords.size() > 0)
        {
            float distance = 0;
            for (int i = 1; i < coords.size(); ++i)
            {
                geos::geom::Coordinate coord0 = coords[i - 1];
                geos::geom::Coordinate coord1 = coords[i];
                distance += GeoUtil::getDistance2d(coord0.x, coord0.y, coord1.x, coord1.y);
            }
            return distance;
        }
    }
    catch (const std::exception &e)
    {
        std::cerr << "Standard Exception: " << e.what() << std::endl;
    }

    return 0;
}

float DuckdbUDF::line3PointMinAngle(const std::string &geomWkt)
{
    try
    {
        std::vector<geos::geom::Coordinate> coords = GeoUtil::getCoordsFromWkt(geomWkt);
        std::cout << geomWkt << std::endl;
        std::cout << coords.size() << std::endl;
        if (coords.size() > 2)
        {
            float minAngle = 180;
            for (int i = 0; i < coords.size() - 2; ++i)
            {
                geos::geom::Coordinate coord0 = coords[i];
                geos::geom::Coordinate coord1 = coords[i + 1];
                geos::geom::Coordinate coord2 = coords[i + 2];
                Eigen::Vector2d ba(coord0.x - coord1.x, coord0.y - coord1.y);
                Eigen::Vector2d bc(coord2.x - coord1.x, coord2.y - coord1.y);
                float angle = GeoUtil::getVectorAngle(ba, bc);
                if (angle < minAngle)
                {
                    minAngle = angle;
                }
            }
            return minAngle;
        }
    }
    catch (const std::exception &e)
    {
        std::cerr << "Standard Exception: " << e.what() << std::endl;
    }

    return 0;
}