diff --git a/.gitignore b/.gitignore index 6c3feba..c4dbf21 100644 --- a/.gitignore +++ b/.gitignore @@ -31,4 +31,6 @@ *.out *.app -.DS_Store \ No newline at end of file +.DS_Store +backup/ +build/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..469742a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 3.15.2) +# 设置工程名称 +project(LibTIN VERSION 1.0 LANGUAGES CXX) +# 添加配置配件编写的函数 +include(CMakePackageConfigHelpers) + +message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) +# CMake默认的安装路径 Windows下为C:/Program\ Files/${Project_Name} Linux/Unix下为/usr/local +message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX}) +# CMake默认的变异类型为空 +message(STATUS "Build type: " ${CMAKE_BUILD_TYPE}) + +# 添加源文件地址 +add_subdirectory(src/) \ No newline at end of file diff --git a/LibTINConfig.cmake.in b/LibTINConfig.cmake.in new file mode 100644 index 0000000..62c9d59 --- /dev/null +++ b/LibTINConfig.cmake.in @@ -0,0 +1,15 @@ +@PACKAGE_INIT@ + +set(@PROJECT_NAME@_Version "@PROJECT_VERSION@") +set_and_check(@PROJECT_NAME@_INSTALL_PREFIX "${PACKAGE_PREFIX_DIR}") +set_and_check(@PROJECT_NAME@_INC_DIR "${PACKAGE_PREFIX_DIR}/include") +set_and_check(@PROJECT_NAME@_INCULDE_DIR "${PACKAGE_PREFIX_DIR}/include") +set_and_check(@PROJECT_NAME@_LIB_DIR "${PACKAGE_PREFIX_DIR}/lib") +set_and_check(@PROJECT_NAME@_LIBRARY_DIR "${PACKAGE_PREFIX_DIR}/lib") + +set(@PROJECT_NAME@_LIB tin) +set(@PROJECT_NAME@_LIBRARY tin) +set(@PROJECT_NAME@_FOUND 1) + +# include target information +include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") \ No newline at end of file diff --git a/topo b/data/topo similarity index 100% rename from topo rename to data/topo diff --git a/job.log b/data/topo_TIN.log similarity index 100% rename from job.log rename to data/topo_TIN.log diff --git a/topo.msh b/data/topo_TIN.msh similarity index 100% rename from topo.msh rename to data/topo_TIN.msh diff --git a/topo.neigh b/data/topo_TIN.neigh similarity index 100% rename from topo.neigh rename to data/topo_TIN.neigh diff --git a/demo.cpp b/demo.cpp deleted file mode 100644 index fe01ea5..0000000 --- a/demo.cpp +++ /dev/null @@ -1,146 +0,0 @@ -#include "tin.h" -#include "iostream" -#include "fstream" -#include "iomanip" - -#include "gctl/utility.h" - -#define GRID_FILE "grid-file|grid_file" -#define GRID_PARA "grid-para|grid_para" -#define MAXI_ERR "maxi-error|maxi_error" -#define LOG_FILE "log-file|log_file" -#define NEIGH_FILE "neighbor-file|neighbor_file|neigh-file|neigh_file" -#define SINGLE_Z "single-z" - -int main(int argc, char const *argv[]) -{ - if (argc < 2) - { - std::cerr << "Usage: ./tin \n"; - std::cerr << "Options:\n"; - std::cerr << "grid-file = \n"; - std::cerr << "grid-para = /////\n"; - std::cerr << "maxi-error = 0>\n"; - std::cerr << "log-file = (optional)\n"; - std::cerr << "neighbor-file = (optional)\n"; - std::cerr << "single-z = yes|no\n"; - return -1; - } - - gctl::getoption gopt(argv[1]); - gopt.show_options(); - - // read dem grid - double xmin, xmax, ymin, ymax, dx, dy; - gctl::parse_string_to_value(gopt.get_value(GRID_PARA), '/', xmin, dx, xmax, ymin, dy, ymax); - - int xnum = round((xmax - xmin)/dx) + 1; - int ynum = round((ymax - ymin)/dy) + 1; - std::vector topo(xnum*ynum); - - double tmp_d; - std::ifstream infile(gopt.get_value(GRID_FILE)); - if (gopt.get_value(SINGLE_Z, false) != "NullValue") - { - for (int i = 0; i < xnum*ynum; ++i) - { - infile >> topo[i]; - } - infile.close(); - } - else - { - for (int i = 0; i < xnum*ynum; ++i) - { - infile >> tmp_d >> tmp_d >> topo[i]; - } - infile.close(); - } - - double maxi_err = std::atof(gopt.get_value(MAXI_ERR).c_str()); - std::vector err_records; - std::vector tin_vert; - std::vector tin_ele; - dem2tin(topo, xmin, xmax, ymin, ymax, dx, dy, tin_vert, tin_ele, maxi_err, &err_records); - - // Write a Gmsh's .msh file - std::ofstream outfile(gopt.get_value(GRID_FILE)+".msh"); - outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<id + 1 << " " << std::setprecision(16) - << tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl; - } - outfile<<"$EndNodes"<vert[j]->id + 1; - } - outfile << std::endl; - } - outfile << "$EndElements"<< std::endl; - outfile<<"$NodeData"<id + 1 << " " << std::setprecision(16) << tin_vert[i]->elev << std::endl; - } - outfile << "$EndNodeData" << std::endl; - outfile.close(); - - // write a neighbor file - std::string neigh_name = gopt.get_value(NEIGH_FILE, false); - if (neigh_name != "NullValue") - { - outfile.open(neigh_name+".neigh"); - outfile << tin_ele.size() << std::endl; - for (int i = 0; i < tin_ele.size(); i++) - { - outfile << i + 1; - for (int j = 0; j < 3; j++) - { - if (tin_ele[i]->neigh[j] != nullptr) - { - outfile << " " << tin_ele[i]->neigh[j]->id + 1; - } - else outfile << " -1"; - } - outfile << std::endl; - } - outfile.close(); - } - - // Write a log file - std::string log_name = gopt.get_value(LOG_FILE, false); - if (log_name != "NullValue") - { - std::ofstream logfile(log_name+".log"); - logfile << "# Insertion Maxi-Error\n"; - for (int i = 0; i < err_records.size(); ++i) - { - logfile << i+1 << " " << err_records[i] << std::endl; - } - logfile.close(); - } - - // Destroy memories allocated by the dem2tin function - for (int i = 0; i < tin_vert.size(); ++i) - { - delete tin_vert[i]; - } - - for (int i = 0; i < tin_ele.size(); ++i) - { - delete tin_ele[i]; - } - return 0; -} \ No newline at end of file diff --git a/job.txt b/job.txt deleted file mode 100644 index f01fa00..0000000 --- a/job.txt +++ /dev/null @@ -1,6 +0,0 @@ -grid-file = topo -grid-para = 0/10/1000/0/10/1000 -maxi-error = 1.0 -log-file = job -neighbor-file = topo -single-z = yes \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..14b584d --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,68 @@ +# 设定源文件文件夹 +aux_source_directory(lib LIBTIN_SRC) + +# 以下部分为库的编译 +# 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库 +# 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加 +add_library(tin SHARED ${LIBTIN_SRC}) +# 首先添加静态库的生成命令 +add_library(tin_static STATIC ${LIBTIN_SRC}) +# 设置静态库的输出名称从而获得与动态库名称相同的静态库 +set_target_properties(tin_static PROPERTIES OUTPUT_NAME "tin") +# 设置输出目标属性以同时输出动态库与静态库 +set_target_properties(tin PROPERTIES CLEAN_DIRECT_OUTPUT 1) +set_target_properties(tin_static PROPERTIES CLEAN_DIRECT_OUTPUT 1) +# 设置动态库的版本号 +set_target_properties(tin PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}) +# 设置库文件的输出地址 +set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) + +# 设置编译选项 +set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3") + +set(CONFIG_FILE_PATH lib/cmake/${PROJECT_NAME}) + +configure_package_config_file(${PROJECT_SOURCE_DIR}/${PROJECT_NAME}Config.cmake.in + ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake + INSTALL_DESTINATION ${CONFIG_FILE_PATH} + NO_CHECK_REQUIRED_COMPONENTS_MACRO) + +write_basic_package_version_file(${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake + VERSION ${PROJECT_VERSION} + COMPATIBILITY SameMajorVersion) + +# 库的安装命令 +if(WIN32) + install(TARGETS tin DESTINATION lib) + install(TARGETS tin_static DESTINATION lib) +else() + install(TARGETS tin tin_static + EXPORT ${PROJECT_NAME}Targets + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib) + install(EXPORT ${PROJECT_NAME}Targets + DESTINATION ${CONFIG_FILE_PATH}) + install(FILES + ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake + ${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake + DESTINATION ${CONFIG_FILE_PATH}) +endif() +# 头文件安装命令 +install(FILES lib/tin.h DESTINATION include/tin) + +# 以下部分为例子程序的编译 +# 设置可执行文件的输出地址 +set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin) + +# 例子的编译方法 +macro(add_demo name) + # 添加可执行文件 命令行 + add_executable(${name} demo/${name}.cpp) + # 为安装文件添加动态库的搜索地址 在Windows下并没有什么用 直接忽略 + set_target_properties(${name} PROPERTIES INSTALL_RPATH ${CMAKE_INSTALL_PREFIX}/lib) + # 链接动态库 + target_link_libraries(${name} PUBLIC tin) +endmacro() + +add_demo(demo) diff --git a/src/demo/demo.cpp b/src/demo/demo.cpp new file mode 100644 index 0000000..5619ae9 --- /dev/null +++ b/src/demo/demo.cpp @@ -0,0 +1,94 @@ +#include "../lib/tin.h" +#include "iostream" +#include "fstream" +#include "iomanip" + +int main(int argc, char const *argv[]) +{ + // read dem grid + std::vector topo(10201); + std::ifstream infile("data/topo"); + for (int i = 0; i < 10201; ++i) + { + infile >> topo[i]; + } + infile.close(); + + std::vector err_records; + std::vector tin_vert; + std::vector tin_ele; + dem2tin(topo, 0, 1000, 0, 1000, 10, 10, tin_vert, tin_ele, 1.0, &err_records); + + // Write a log file + std::ofstream logfile("data/topo_TIN.log"); + logfile << "# Insertion Maxi-Error\n"; + for (int i = 0; i < err_records.size(); ++i) + { + logfile << i+1 << " " << err_records[i] << std::endl; + } + logfile.close(); + + // Write a Gmsh's .msh file + std::ofstream outfile("data/topo_TIN.msh"); + outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<id + 1 << " " << std::setprecision(16) + << tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl; + } + outfile<<"$EndNodes"<vert[j]->id + 1; + } + outfile << std::endl; + } + outfile << "$EndElements"<< std::endl; + outfile<<"$NodeData"<id + 1 << " " << std::setprecision(16) << tin_vert[i]->elev << std::endl; + } + outfile << "$EndNodeData" << std::endl; + outfile.close(); + + // write a neighbor file + outfile.open("data/topo_TIN.neigh"); + outfile << tin_ele.size() << std::endl; + for (int i = 0; i < tin_ele.size(); i++) + { + outfile << i + 1; + for (int j = 0; j < 3; j++) + { + if (tin_ele[i]->neigh[j] != nullptr) + { + outfile << " " << tin_ele[i]->neigh[j]->id + 1; + } + else outfile << " -1"; + } + outfile << std::endl; + } + outfile.close(); + + // Destroy memories allocated by the dem2tin function + for (int i = 0; i < tin_vert.size(); ++i) + { + delete tin_vert[i]; + } + + for (int i = 0; i < tin_ele.size(); ++i) + { + delete tin_ele[i]; + } + return 0; +} \ No newline at end of file diff --git a/tin.h b/src/lib/tin.cpp similarity index 80% rename from tin.h rename to src/lib/tin.cpp index 3b207dc..85e3ef3 100644 --- a/tin.h +++ b/src/lib/tin.cpp @@ -1,37 +1,56 @@ /** - * @defgroup TIN + * ___________ __ + * /_ __/ _/ | / / + * / / / // |/ / + * / / _/ // /| / + * /_/ /___/_/ |_/ + * + * C++ library of the Triangular Irregular Network (TIN) * - * @brief Generation of a Triangular Irregular Network (TIN) from a dense DEM grid. + * Copyright (c) 2021-2031 Yi Zhang (zhangyiss@icloud.com) + * All rights reserved. * - * @author Yi Zhang - * @date 2021-09-16 + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. */ -#ifndef _TIN_DELAUNAY_H -#define _TIN_DELAUNAY_H -#include "cmath" -#include "vector" -#include "algorithm" - -#define ZERO 1e-5 +#include "tin.h" // Start vertex definition -struct vertex2dc +vertex2dc::vertex2dc() { - unsigned int id; // index of the vertex - double x, y; // position of the vertex - double elev; // elevation at the vertex + x = y = elev = NAN; + id = 0; +} - vertex2dc() : x(NAN), y(NAN), elev(NAN), id(0) {} - vertex2dc(double inx, double iny, double inelev, unsigned int inid) {set(inx, iny, inelev, inid);} - void set(double inx, double iny, double inelev, unsigned int inid) - { - x = inx; y = iny; elev = inelev; id = inid; - return; - } -}; +vertex2dc::vertex2dc(double inx, double iny, double inelev, unsigned int inid) +{ + set(inx, iny, inelev, inid); +} -bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == operator for vertex2dc type +void vertex2dc::set(double inx, double iny, double inelev, unsigned int inid) +{ + x = inx; y = iny; elev = inelev; id = inid; + return; +} + +// overload the == operator for vertex2dc type +bool operator ==(const vertex2dc &a, const vertex2dc &b) { if(fabs(a.x - b.x) <= ZERO && fabs(a.y - b.y) <= ZERO) { @@ -40,7 +59,8 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper return false; } -bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test if the three points are on the same line +// Test if the three points are on the same line +bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) { // ｜(y3−y1)(x2−x1)−(y2−y1)(x3−x1)｜ if (fabs((c_ptr->y - a_ptr->y)*(b_ptr->x - a_ptr->x) - (b_ptr->y - a_ptr->y)*(c_ptr->x - a_ptr->x)) <= ZERO) @@ -50,7 +70,8 @@ bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test return false; } -void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points +// calculate the circumcircle from three points +void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) { double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x)); double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y; @@ -64,23 +85,23 @@ void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, doubl // End vertex definition // Start DEM definition -struct triangle; - -struct dem_point +dem_point::dem_point() { - double x, y; // position of the DEM location - double elev; // elevation at the DEM location - double err; // error of the TIN with respect to the elevation - triangle *host; + x = y = elev = NAN; + err = 0.0; + host = nullptr; +} - dem_point() : x(NAN), y(NAN), elev(NAN), err(0.0), host(nullptr) {} - dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);} - void set(double inx, double iny, double inelev) - { - x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr; - return; - } -}; +dem_point::dem_point(double inx, double iny, double inelev) +{ + set(inx, iny, inelev); +} + +void dem_point::set(double inx, double iny, double inelev) +{ + x = inx; y = iny; elev = inelev; err = 0.0; host = nullptr; + return; +} bool compare_dem_point(dem_point *a, dem_point *b) // determination function for std::sort { @@ -89,66 +110,59 @@ bool compare_dem_point(dem_point *a, dem_point *b) // determination function for } // End DEM definition -/* Start triangle definition - * v2 - * /\ - * / \ - * n2 / \ n1 - * / \ - * /------------\ - * v0 n0 v1 - */ -struct triangle +// Start triangle definition +triangle::triangle() { - int id; - vertex2dc *vert[3]; // vertex of the triangle - triangle *neigh[3]; // neighbors of the triangle - double cx, cy; // center of the triangle's circumcircle - double cr; // radius of the circumcircle - std::vector hosted_dem; + vert[0] = vert[1] = vert[2] = nullptr; + neigh[0] = neigh[1] = neigh[2] = nullptr; +} - triangle() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[2] = nullptr;} - triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);} - void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) - { - vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr; - neigh[0] = neigh[1] = neigh[2] = nullptr; - circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); - return; - } +triangle::triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) +{ + set(v0ptr, v1ptr, v2ptr); +} - void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr) - { - neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr; - return; - } +void triangle::set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) +{ + vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr; + neigh[0] = neigh[1] = neigh[2] = nullptr; + circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); + return; +} - bool bound_location(double inx, double iny) // Test if the location is inside the triangle - { - double l1x, l1y, l2x, l2y; - for (int i = 0; i < 3; i++) - { - l1x = vert[(i+1)%3]->x - vert[i]->x; - l1y = vert[(i+1)%3]->y - vert[i]->y; - l2x = inx - vert[i]->x; - l2y = iny - vert[i]->y; +void triangle::set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr) +{ + neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr; + return; +} - if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge - { - return false; - } - } - return true; - } +// Test if the location is inside the triangle +bool triangle::bound_location(double inx, double iny) +{ + double l1x, l1y, l2x, l2y; + for (int i = 0; i < 3; i++) + { + l1x = vert[(i+1)%3]->x - vert[i]->x; + l1y = vert[(i+1)%3]->y - vert[i]->y; + l2x = inx - vert[i]->x; + l2y = iny - vert[i]->y; - double interpolate(double inx, double iny) // Interpolate the elevation of the given location inside the triangle - { - double a1 = 0.5 * ((vert[1]->x - inx)*(vert[2]->y - iny) - (vert[1]->y - iny)*(vert[2]->x - inx)); - double a2 = 0.5 * ((vert[2]->x - inx)*(vert[0]->y - iny) - (vert[2]->y - iny)*(vert[0]->x - inx)); - double a3 = 0.5 * ((vert[0]->x - inx)*(vert[1]->y - iny) - (vert[0]->y - iny)*(vert[1]->x - inx)); - return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3); - } -}; + if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge + { + return false; + } + } + return true; +} + +// Interpolate the elevation of the given location inside the triangle +double triangle::interpolate(double inx, double iny) +{ + double a1 = 0.5 * ((vert[1]->x - inx)*(vert[2]->y - iny) - (vert[1]->y - iny)*(vert[2]->x - inx)); + double a2 = 0.5 * ((vert[2]->x - inx)*(vert[0]->y - iny) - (vert[2]->y - iny)*(vert[0]->x - inx)); + double a3 = 0.5 * ((vert[0]->x - inx)*(vert[1]->y - iny) - (vert[0]->y - iny)*(vert[1]->x - inx)); + return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3); +} /** * @brief Flip neighboring triangles and their neighbors @@ -447,9 +461,8 @@ triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4]) * @param[in] maxi_err Threshold to quit the algorithm. The default is 1e-0 * @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex. */ -void dem2tin(const std::vector &dem, double xmin, double xmax, double ymin, double ymax, - double dx, double dy, std::vector &out_verts, std::vector &out_tris, - double maxi_err = 1e-0, std::vector *err_records = nullptr) +void dem2tin(const std::vector &dem, double xmin, double xmax, double ymin, double ymax, double dx, + double dy, std::vector &out_verts, std::vector &out_tris, double maxi_err, std::vector *err_records) { if (!out_verts.empty()) out_verts.clear(); if (!out_tris.empty()) out_tris.clear(); @@ -712,5 +725,3 @@ void dem2tin(const std::vector &dem, double xmin, double xmax, double ym } return; } - -#endif // _TIN_DELAUNAY_H \ No newline at end of file diff --git a/src/lib/tin.h b/src/lib/tin.h new file mode 100644 index 0000000..b298fd8 --- /dev/null +++ b/src/lib/tin.h @@ -0,0 +1,223 @@ +/** + * ___________ __ + * /_ __/ _/ | / / + * / / / // |/ / + * / / _/ // /| / + * /_/ /___/_/ |_/ + * + * C++ library of the Triangular Irregular Network (TIN) + * + * Copyright (c) 2021-2031 Yi Zhang (zhangyiss@icloud.com) + * All rights reserved. + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +#ifndef _TIN_DELAUNAY_H +#define _TIN_DELAUNAY_H +#include "cmath" +#include "vector" +#include "algorithm" + +#define ZERO 1e-5 + +// Start vertex definition +struct vertex2dc +{ + unsigned int id; // index of the vertex + double x, y; // position of the vertex + double elev; // elevation at the vertex + + vertex2dc(); + + /** + * @brief Construct a new vertex2dc object + * + * @param inx initial x coordinate + * @param iny initial y coordinate + * @param inelev initial elevation + * @param inid initial index + */ + vertex2dc(double inx, double iny, double inelev, unsigned int inid); + + /** + * @brief Set a vertex2dc object + * + * @param inx initial x coordinate + * @param iny initial y coordinate + * @param inelev initial elevation + * @param inid initial index + */ + void set(double inx, double iny, double inelev, unsigned int inid); +}; + +/** + * @brief Compare two vertexes + * + * @param a vertex a + * @param b vertex b + * @return true the two vertexes are at the same location + * @return false the two vertexes are not at the same location + */ +bool operator ==(const vertex2dc &a, const vertex2dc &b); + +/** + * @brief Test if the three points are on the same line + * + * @param a_ptr pointer of the vertex a + * @param b_ptr pointer of the vertex b + * @param c_ptr pointer of the vertex c + * @return true the three vertexes are on the same line + * @return false the three vertexes are not on the same line + */ +bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr); + +/** + * @brief Calculate the circumcircle from three points + * + * @param v0 pointer of the vertex v0 + * @param v1 pointer of the vertex v1 + * @param v2 pointer of the vertex v2 + * @param cx x coordinate of the returned circumcircle + * @param cy y coordinate of the returned circumcircle + * @param cr squared radius of the returned circumcircle + */ +void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr); +// End vertex definition + +// Start DEM definition +struct triangle; + +struct dem_point +{ + double x, y; // position of the DEM location + double elev; // elevation at the DEM location + double err; // error of the TIN with respect to the elevation + triangle *host; + + dem_point(); + + /** + * @brief Construct a new DEM point object + * + * @param inx initial x coordinate + * @param iny initial y coordinate + * @param inelev initial elevation + */ + dem_point(double inx, double iny, double inelev); + + /** + * @brief Set a DEM point object + * + * @param inx initial x coordinate + * @param iny initial y coordinate + * @param inelev initial elevation + */ + void set(double inx, double iny, double inelev); +}; +// End DEM definition + +/* Start triangle definition + * v2 + * /\ + * / \ + * n2 / \ n1 + * / \ + * /------------\ + * v0 n0 v1 + */ +struct triangle +{ + int id; + vertex2dc *vert[3]; // vertex of the triangle + triangle *neigh[3]; // neighbors of the triangle + double cx, cy; // center of the triangle's circumcircle + double cr; // radius of the circumcircle + std::vector hosted_dem; + + triangle(); + + /** + * @brief Construct a new triangle object + * + * @param v0ptr pointer of the vertex 0 + * @param v1ptr pointer of the vertex 1 + * @param v2ptr pointer of the vertex 2 + */ + triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr); + + /** + * @brief Set a triangle object + * + * @param v0ptr pointer of the vertex 0 + * @param v1ptr pointer of the vertex 1 + * @param v2ptr pointer of the vertex 2 + */ + void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr); + + /** + * @brief Set neighbors of a triangle object + * + * @param n0ptr pointer of the neighboring vertex 0 + * @param n1ptr pointer of the neighboring vertex 1 + * @param n2ptr pointer of the neighboring vertex 2 + */ + void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr); + + /** + * @brief Test if the location is inside the triangle + * + * @param inx x coordinate of the input location + * @param iny y coordinate of the input location + * @return true the input location is inside the triangle + * @return false the input location is not inside the triangle + */ + bool bound_location(double inx, double iny); + + /** + * @brief Interpolate the elevation of the given location inside the triangle + * + * @param inx x coordinate of the input location + * @param iny y coordinate of the input location + * @return double the interpolated elevation at the input location + */ + double interpolate(double inx, double iny); +}; +// End triangle definition + +/** + * @brief Generate the TIN from the DEM grid + * + * @param[in] dem Input DEM grid (Ordered from lower left corner to the upper right corner) + * @param[in] xmin The minimal coordinate of the DEM grid on the x-axis + * @param[in] xmax The maximal coordinate of the DEM grid on the x-axis + * @param[in] ymin The minimal coordinate of the DEM grid on the y-axis + * @param[in] ymax The maximal coordinate of the DEM grid on the y-axis + * @param[in] dx Data spacing of the DEM grid on the x-axis + * @param[in] dy Data spacing of the DEM grid on the y-axis + * @param out_verts The output vector of vertex's pointers. The user need to destroy the memories allocated by the function before destroy the vector + * @param out_tris The output vector of triangle's pointers. The user need to destroy the memories allocated by the function before destroy the vector + * @param[in] maxi_err Threshold to quit the algorithm. The default is 1e-0 + * @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex. + */ +void dem2tin(const std::vector &dem, double xmin, double xmax, double ymin, double ymax, + double dx, double dy, std::vector &out_verts, std::vector &out_tris, + double maxi_err = 1e-0, std::vector *err_records = nullptr); + +#endif // _TIN_DELAUNAY_H \ No newline at end of file diff --git a/tin b/tin deleted file mode 100755 index f7ff42b..0000000 Binary files a/tin and /dev/null differ diff --git a/tin_backup.h b/tin_backup.h deleted file mode 100644 index d13f24c..0000000 --- a/tin_backup.h +++ /dev/null @@ -1,461 +0,0 @@ -/** - * @defgroup TIN - * - * @brief Generation of a Triangular Irregular Network (TIN) from a dense DEM grid. - * - * @author Yi Zhang - * @date 2021-09-16 - */ - -#ifndef _TIN_DELAUNAY_H -#define _TIN_DELAUNAY_H -#include "cmath" -#include "vector" -#include "algorithm" - -#define ZERO 1e-5 - -// Start vertex definition -struct vertex2dc -{ - unsigned int id; // index of the vertex - double x, y; // position of the vertex - double elev; // elevation at the vertex - - vertex2dc() : x(NAN), y(NAN), elev(NAN), id(0) {} - vertex2dc(double inx, double iny, double inelev, unsigned int inid) {set(inx, iny, inelev, inid);} - void set(double inx, double iny, double inelev, unsigned int inid) - { - x = inx; y = iny; elev = inelev; id = inid; - return; - } -}; - -bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == operator for vertex2dc type -{ - if(fabs(a.x - b.x) <= ZERO && fabs(a.y - b.y) <= ZERO) - { - return true; - } - return false; -} - -bool is_collinear(vertex2dc *a_ptr, vertex2dc *b_ptr, vertex2dc *c_ptr) // Test if the three points are on the same line -{ - // ｜(y3−y1)(x2−x1)−(y2−y1)(x3−x1)｜ - if (fabs((c_ptr->y - a_ptr->y)*(b_ptr->x - a_ptr->x) - (b_ptr->y - a_ptr->y)*(c_ptr->x - a_ptr->x)) <= ZERO) - { - return true; - } - return false; -} -// End vertex definition - -// Start edge definition -struct edge -{ - vertex2dc *vert[2]; // vertex of the edge - - edge() {vert[0] = vert[1] = nullptr;} - edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);} - void set(vertex2dc *v0ptr, vertex2dc *v1ptr) - { - vert[0] = v0ptr; vert[1] = v1ptr; - return; - } -}; - -bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type -{ - if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) || - (a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0])) - { - return true; - } - return false; -} -// End edge definition - -// Start triangle definition -struct dem_point; - -struct triangle -{ - vertex2dc *vert[3]; // vertex of the triangle - double cx, cy; // center of the triangle's circumcircle - double cr; // radius of the circumcircle - std::vector circum_dem; - - triangle() {vert[0] = vert[1] = vert[2] = nullptr;} - triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);} - void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) - { - vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr; - - double s = 0.5 / ((vert[1]->x - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x)); - double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y; - double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y; - - cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s; - cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s; - cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here - return; - } - - bool bound_location(double inx, double iny) // Test if the location is inside the triangle - { - double l1x, l1y, l2x, l2y; - for (int i = 0; i < 3; i++) - { - l1x = vert[(i+1)%3]->x - vert[i]->x; - l1y = vert[(i+1)%3]->y - vert[i]->y; - l2x = inx - vert[i]->x; - l2y = iny - vert[i]->y; - - if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge - { - return false; - } - } - return true; - } - - double interpolate(double inx, double iny) // Interpolate the elevation of the given location inside the triangle - { - double a1 = 0.5 * ((vert[1]->x - inx)*(vert[2]->y - iny) - (vert[1]->y - iny)*(vert[2]->x - inx)); - double a2 = 0.5 * ((vert[2]->x - inx)*(vert[0]->y - iny) - (vert[2]->y - iny)*(vert[0]->x - inx)); - double a3 = 0.5 * ((vert[0]->x - inx)*(vert[1]->y - iny) - (vert[0]->y - iny)*(vert[1]->x - inx)); - return (a1*vert[0]->elev + a2*vert[1]->elev + a3*vert[2]->elev)/(a1 + a2 + a3); - } -}; -// End triangle definition - -// Start DEM definition -struct dem_point -{ - double x, y; // position of the DEM location - double elev; // elevation at the DEM location - double err; // error of the TIN with respect to the elevation - triangle *host; // host triangle of the DEM location - std::vector circum_host; // triangles which circumcircles include the location - - dem_point() : x(NAN), y(NAN), elev(NAN), host(nullptr) {} - dem_point(double inx, double iny, double inelev) {set(inx, iny, inelev);} - void set(double inx, double iny, double inelev) - { - x = inx; y = iny; elev = inelev; host = nullptr; - return; - } -}; - -bool compare_dem_point(dem_point *a, dem_point *b) -{ - if (a->err > b->err) return true; - return false; -} -// End DEM definition - -/** - * @brief Generate the TIN from the DEM grid - * - * @param[in] dem Input DEM grid (Ordered from lower left corner to the upper right corner) - * @param[in] xmin The minimal coordinate of the DEM grid on the x-axis - * @param[in] xmax The maximal coordinate of the DEM grid on the x-axis - * @param[in] ymin The minimal coordinate of the DEM grid on the y-axis - * @param[in] ymax The maximal coordinate of the DEM grid on the y-axis - * @param[in] dx Data spacing of the DEM grid on the x-axis - * @param[in] dy Data spacing of the DEM grid on the y-axis - * @param out_verts The output vector of vertex's pointers. The user need to destroy the memories allocated by the function before destroy the vector - * @param out_tris The output vector of triangle's pointers. The user need to destroy the memories allocated by the function before destroy the vector - * @param[in] maxi_err Threshold to quit the algorithm. The default is 1e-0 - * @param[in] err_records If this pointer is not NULL, record maximal error values after each insertion of vertex. - */ -void dem2tin(const std::vector &dem, double xmin, double xmax, double ymin, double ymax, - double dx, double dy, std::vector &out_verts, std::vector &out_tris, - double maxi_err = 1e-0, std::vector *err_records = nullptr) -{ - if (!out_verts.empty()) out_verts.clear(); - if (!out_tris.empty()) out_tris.clear(); - if (err_records != nullptr && !err_records->empty()) err_records->clear(); - - if (dx <= 0.0 || dy <= 0.0 || maxi_err <= 0.0) return; - if (xmin >= xmax || ymin >= ymax || (xmin + dx) > xmax || (ymin + dy) > ymax) return; - - int xnum = round((xmax - xmin)/dx) + 1; - int ynum = round((ymax - ymin)/dy) + 1; - - if (dem.size() != xnum*ynum) return; - - // Prepare the DEM points - dem_point *tmp_dem = nullptr;; - std::vector dem_grid(xnum*ynum); - std::vector::iterator d_iter; - for (int i = 0; i < ynum; ++i) - { - for (int j = 0; j < xnum; ++j) - { - dem_grid[j + i*xnum] = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]); - } - } - - vertex2dc *tmp_vert = nullptr; - tmp_vert = new vertex2dc(xmin, ymin, dem_grid[0]->elev, out_verts.size()); // lower left corner - out_verts.push_back(tmp_vert); - - d_iter = dem_grid.begin(); - tmp_dem = *d_iter; delete tmp_dem; - dem_grid.erase(d_iter); - - tmp_vert = new vertex2dc(xmax, ymin, dem_grid[xnum-2]->elev, out_verts.size()); // lower right corner. Note the first location is already erased - out_verts.push_back(tmp_vert); - - d_iter = dem_grid.begin() + (xnum - 2); - tmp_dem = *d_iter; delete tmp_dem; - dem_grid.erase(d_iter); - - tmp_vert = new vertex2dc(xmax, ymax, dem_grid[xnum*ynum-3]->elev, out_verts.size()); // upper right corner. Note the first two locations are already erased - out_verts.push_back(tmp_vert); - - d_iter = dem_grid.begin() + (xnum*ynum - 3); - tmp_dem = *d_iter; delete tmp_dem; - dem_grid.erase(d_iter); - - tmp_vert = new vertex2dc(xmin, ymax, dem_grid[xnum*(ynum-1) - 2]->elev, out_verts.size()); // upper left corner. Note the first two locations are already erased - out_verts.push_back(tmp_vert); - - d_iter = dem_grid.begin() + (xnum*(ynum-1) - 2); - tmp_dem = *d_iter; delete tmp_dem; - dem_grid.erase(d_iter); - - triangle *tmp_tri = nullptr; - std::vector cnst_tri, new_tri; - std::vector::iterator t_iter; - - if (!is_collinear(out_verts[0], out_verts[1], out_verts[2])) // Do not create triangle if the vertexes are collinear - { - tmp_tri = new triangle(out_verts[0], out_verts[1], out_verts[2]); // order the vertex anti-clock wise - out_tris.push_back(tmp_tri); - } - - if (!is_collinear(out_verts[0], out_verts[2], out_verts[3])) - { - tmp_tri = new triangle(out_verts[0], out_verts[2], out_verts[3]); // order the vertex anti-clock wise - out_tris.push_back(tmp_tri); - } - - // Find host triangle for all DEM locations - for (int i = 0; i < dem_grid.size(); ++i) - { - for (int t = 0; t < out_tris.size(); ++t) - { - if (out_tris[t]->bound_location(dem_grid[i]->x, dem_grid[i]->y)) - { - dem_grid[i]->host = out_tris[t]; - break; // already found, no need to search more - } - } - } - - // Find circum_host triangles for all DEM locations - double dist; - for (int i = 0; i < dem_grid.size(); ++i) - { - for (int t = 0; t < out_tris.size(); ++t) - { - dist = (out_tris[t]->cx - dem_grid[i]->x) * (out_tris[t]->cx - dem_grid[i]->x) - + (out_tris[t]->cy - dem_grid[i]->y) * (out_tris[t]->cy - dem_grid[i]->y); - if ((dist - out_tris[t]->cr) <= ZERO) // Points on the circumcircle are also included - { - dem_grid[i]->circum_host.push_back(out_tris[t]); - out_tris[t]->circum_dem.push_back(dem_grid[i]); - // no beak here. There might be more than one triangle's circumcircle includes the DEM location - } - } - } - - // loop all DEM data to find the location with maximal error - for (int i = 0; i < dem_grid.size(); ++i) - { - dem_grid[i]->err = fabs(dem_grid[i]->host->interpolate(dem_grid[i]->x, dem_grid[i]->y) - dem_grid[i]->elev); - } - - // Sort dem_grid in the desceding order with respect to the error - std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point); - - bool removed; - edge tmp_edge; - std::vector cnst_edge; - std::vector::iterator e_iter; - - while (dem_grid[0]->err >= maxi_err) // quit til the threshold is meet - { - if (err_records != nullptr) - { - err_records->push_back(dem_grid[0]->err); - } - - // create a new vertex - tmp_vert = new vertex2dc(dem_grid[0]->x, dem_grid[0]->y, dem_grid[0]->elev, out_verts.size()); - out_verts.push_back(tmp_vert); - - // Move triangles which circumcircles include the new vertex to the cnst_tri and remove it from out_tris - cnst_tri.clear(); - for (int i = 0; i < dem_grid[0]->circum_host.size(); ++i) - { - cnst_tri.push_back(dem_grid[0]->circum_host[i]); - } - - for (int c = 0; c < cnst_tri.size(); ++c) - { - for (t_iter = out_tris.begin(); t_iter != out_tris.end(); ) - { - tmp_tri = *t_iter; - if (cnst_tri[c] == tmp_tri) - { - t_iter = out_tris.erase(t_iter); - break; // no need to search more - } - else t_iter++; - } - } - - // remove cnst_tri from its circumed DEM's circum triangle list - for (int c = 0; c < cnst_tri.size(); ++c) - { - for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i) - { - tmp_dem = cnst_tri[c]->circum_dem[i]; - for (t_iter = tmp_dem->circum_host.begin(); t_iter != tmp_dem->circum_host.end(); ) - { - if (cnst_tri[c] == *t_iter) - { - t_iter = tmp_dem->circum_host.erase(t_iter); - break; - } - else t_iter++; - } - } - } - - // remove dem_grid[0] from its circumed triangle's circum DEM list - for (int c = 0; c < cnst_tri.size(); ++c) - { - for (d_iter = cnst_tri[c]->circum_dem.begin(); d_iter != cnst_tri[c]->circum_dem.end(); ) - { - if (dem_grid[0] == *d_iter) - { - d_iter = cnst_tri[c]->circum_dem.erase(d_iter); - break; - } - else d_iter++; - } - } - - // clear host and circumcircle triangles for the used DEM location - d_iter = dem_grid.begin(); - tmp_dem = *d_iter; tmp_dem->circum_host.clear(); delete tmp_dem; - dem_grid.erase(d_iter); - - // loop to remove duplicate edges - cnst_edge.clear(); - for (int c = 0; c < cnst_tri.size(); ++c) - { - for (int e = 0; e < 3; ++e) - { - tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]); - - removed = false; - for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); ) - { - if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge - { - e_iter = cnst_edge.erase(e_iter); - removed = true; - break; // no need to search more - } - else e_iter++; - } - - if (!removed) // not a duplicate edge, add to the cnst_edge - { - cnst_edge.push_back(tmp_edge); - } - } - } - - // construct new triangles and add to out_tris - new_tri.clear(); - for (int c = 0; c < cnst_edge.size(); ++c) - { - if (!is_collinear(cnst_edge[c].vert[0], cnst_edge[c].vert[1], tmp_vert)) // Do not create triangle if the vertexes are collinear - { - tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], tmp_vert); // order the vertex anti-clock wise - out_tris.push_back(tmp_tri); - new_tri.push_back(tmp_tri); - } - } - - // loop all DEM data to update host triangles - for (int c = 0; c < cnst_tri.size(); ++c) - { - for (int i = 0; i < cnst_tri[c]->circum_dem.size(); ++i) - { - tmp_dem = cnst_tri[c]->circum_dem[i]; - for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new host - { - if (new_tri[n]->bound_location(tmp_dem->x, tmp_dem->y)) - { - tmp_dem->host = new_tri[n]; - tmp_dem->err = fabs(new_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev); - break; // already found, no need to search more - } - } - } - } - - // Find circum_host triangles for all DEM locations - // cnst_tri's circum area doesn't over cover new_tri's circum area - for (int i = 0; i < dem_grid.size(); ++i) - { - for (int n = 0; n < new_tri.size(); ++n) // search in newly created triangles to find new circumcircle triangles - { - dist = (new_tri[n]->cx - dem_grid[i]->x) * (new_tri[n]->cx - dem_grid[i]->x) - + (new_tri[n]->cy - dem_grid[i]->y) * (new_tri[n]->cy - dem_grid[i]->y); - if ((dist - new_tri[n]->cr) <= ZERO) // Points on the circumcircle are also included - { - new_tri[n]->circum_dem.push_back(dem_grid[i]); - dem_grid[i]->circum_host.push_back(new_tri[n]); - // no beak here. There might be more than one triangle's circumcircle includes the DEM location - } - } - } - - - // destroy memories used by cnst_edge - for (int c = 0; c < cnst_tri.size(); ++c) - { - tmp_tri = cnst_tri[c]; - tmp_tri->circum_dem.clear(); - delete tmp_tri; tmp_tri = nullptr; - } - - // Sort dem_grid in the desceding order with respect to the error - std::sort(dem_grid.begin(), dem_grid.end(), compare_dem_point); - } - - if (err_records != nullptr) - { - err_records->push_back(dem_grid[0]->err); - } - - // destroy remaining DEM data - for (int i = 0; i < dem_grid.size(); ++i) - { - tmp_dem = dem_grid[i]; - delete tmp_dem; tmp_dem = nullptr; - } - - return; -} - -#endif // _TIN_DELAUNAY_H \ No newline at end of file