From 0b52bae14ed33b2c365b134cee9fb4bd055d7192 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Thu, 14 Oct 2021 21:38:45 +0800 Subject: [PATCH] update to tin library --- .gitignore | 4 +- CMakeLists.txt | 14 + LibTINConfig.cmake.in | 15 + topo => data/topo | 0 job.log => data/topo_TIN.log | 0 topo.msh => data/topo_TIN.msh | 0 topo.neigh => data/topo_TIN.neigh | 0 demo.cpp | 146 ---------- job.txt | 6 - src/CMakeLists.txt | 68 +++++ src/demo/demo.cpp | 94 ++++++ tin.h => src/lib/tin.cpp | 211 +++++++------- src/lib/tin.h | 223 +++++++++++++++ tin | Bin 103360 -> 0 bytes tin_backup.h | 461 ------------------------------ 15 files changed, 528 insertions(+), 714 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 LibTINConfig.cmake.in rename topo => data/topo (100%) rename job.log => data/topo_TIN.log (100%) rename topo.msh => data/topo_TIN.msh (100%) rename topo.neigh => data/topo_TIN.neigh (100%) delete mode 100644 demo.cpp delete mode 100644 job.txt create mode 100644 src/CMakeLists.txt create mode 100644 src/demo/demo.cpp rename tin.h => src/lib/tin.cpp (80%) create mode 100644 src/lib/tin.h delete mode 100755 tin delete mode 100644 tin_backup.h 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 f7ff42b93da4108f1abcf11f5f01cbda26b7a4a2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 103360 zcmeFa3wV^(weY{+nOtU)0MT5zB@;kR5ET@GK`@g9Xz+p-qf)g^LO>@#0J$irNB}Q` zK_yBnReK0tYG#BUJO^9w*h5fjW3jC$Ra-rE5Lze3o`N7|kTBnGy_e)oNYI}1edquG zzQ;Vzo4xjJ?X}lld+oK?-tRlVKmNCqDN1RIKL<}5Pg8%2Qr`b6Dn@Sw;ipo;{xQ)WD#WO)+J?DFy@vzIL~BLv;? zy_I2fX2eFLjKrm19e;7EwgUzzFJHW5>CD9)N!{_C*lOS#<}`{(+>Xy>;)93!pMkI3 zKYQN7*^5Z&j?eRgfiK0xCviJIbqamAR8^JFURE`G;ga$tl{35J8?oKMw}CN}KZ$q2 zXC!z0FE6iJFn@6;1H0pEdda|dp*g4$kL##?DgeI4OBPknzs1bzj_>n-Gw@ZJL`mHK z*PYh#@}=`{ubw}vyn6o43nap;O~N0UL{ z`zG@xZiq3{z(;IBX%gbj~E~)HHF!X{eJ zT0Uo#t$JSpy z{TDo=E#l1#X?GltUvFWJDniurN906EeF2kz_%RumIL2%;9CxS%Ykn>@GS?v<-oTb_?83Ta^U|`4#fG| z^LW(0{dec*JovVM&_CH7loqcqOoMRp$8ig?nu~eE+5af=}J!s~NJ+ zS2KxvQ^4~ab?b#t|0kddhCb%ndq}_$Cs+{8dg*hW> ztC%vcyKc^DBK=v{?3{-DiX27SUnrk9t|}*wI-Bz=b5>EdQmJ+g-jwhRC9fDfLxDT! zrfx0eYbjT$eC-bCD}tUp>Q;bX+V@bVfqcCWDe#a-cp84s=vU0?*n47uBY6B4XRxh2 zE_m$v_~6m25`yh_CI;JI)#{Hpwa`(g3Jv$D`YZle)m9p(&UhpA!1Q1!v?zH+!4oy@ zE8m&hUh&RctxxjY_~c}lmVCKu#fnFhDL=0|zb$3EFW5K57hSd1SMzkeuciX}UjL!+ zPxy07&U(VVqzV60;e&#Iso4v2s_3iqql&zG#;Yp-)|~IcqpESY>ipYHy6{hUl}FlN zgvaBA@0B@E<hgJor!qAEd3>wDH`JeKl{<=CiJaIrW5}lQt@8QO>jUM`a7WojkZ#=fUip z$vO{a=X?V1_ETPDX*YD{(ZBelZakpu5O^?&K1;tX9^BiF2a3LW;6c%YzM4Gf=q^7V zcrd({@volScI90w9_U4!J?^E2iQcG*~Ft_#igGL-#ElcWC(pF9$@~`KWfg@k@ zMWt>Fb%#^8TuayM>T&8yT(5CGR5YMx5!P%s!LbBFQT1iuXl7a=2P?Se8{6 zYxOOODthO&Kb(H!wJ$Ehp_kt5$gT^Vuk<4u>78Z1`Yx z=H)B??TNYVmG6kmUxv)jnKry`?j6-_mw~@%Rnvr-JbC**p61{w_BrO~eK>Rekhtt= zl2#C}>hnDnIrp~5$_m_SZ($#GxFkW9J({MX{+~_`=4$TNTIe5iuWIu{%TVZQN%Td> z!Gr!<_MWk^l1C#6s@~nc;+?4;^ zP5zg`|A3AE1K_{Ng3kYBUvw<|&(^Z`421vwM1%ip)76XaZv0>A?1BGh8vOs`Fts5H z|1aq`u^`v7=Zt`>-^6Oy+&Q(zaRc*b>V+jC`w*W z{O=F{?QL27uO9A;CLw1{@v6*4nR!3+1y`-|)vWI0YoFKG7p+1^-G`pOh4^LdC6CIw zRu6pca(~;RD}6Of#x2TuaopmZA0;i#d5rLV;?bmKIq!~No|8sAJ9%-=g;EB6mT}3_ zoOhEJ0x!(@ocN@aM= z*N(@^mSp;(+l3zLccDkeeHwaf`Q7R9pwFvK+M50WdS=u2Ad7hSX`H!cKs5|#zUz@fD zyV_fkBl<>cJ+X0Xv0XK9r)~Zuy0#!s)fbuj0T=7$6N7z^Jk~d}E&|V89<}`|*71R? zy(^vBt+UcDKfd`wRW{Jw8LyOhpU$}N231x{KOWLly6RIqeKB(;@Kawf_pbP0?nvx} zXH+obgci&_;RvdGp%L3iZ2z;M>mKaLo1rxkx?jf*yJy_ZIZ5cpN&UOp{WptEntyZ7 z+1OXblu2~^bHd;_$6cKh#x7T+J&TRj-R{@0*__yDv&Z@2$LyT5fS<vy;k8+QR&x2(FKfQ;OwHW?OWqdeKo!H@9X-OOy4Tt#|ZjXTkwClZwW@< z5@gI6|Ftqt{)^-Gb^Ve){U7g7l8UZH*ZQ$_-bhs2-QPpz`WtBxr#qBnvv>8l&ECz& z-p$3{U4=f9ylQZ`v3cvzmHv<~DDl>KwLQ7O5p)Yq(*I$c@K515jg=X8VC+Y7c zsJ$ZVtA`k}9!XMXxUkWy(N9y+Pwo}LCkoI#vM+J5FUdN*;+@AGS#_nXnYkyJ@A)d4 zy-(?L-Nl?&_xoxbjC&mOOZI6p2hV2S-_P8Y`7U!==6iAfnK>f&vd_NRl>Pe(Zq7M} z`TqJi^a6V}(PgratI+5BEjdrKU&|X;%^nN?1N*ij%Ifny(bpcwJfFqh^fl@^z0{#h zE%kha`h0)*>*sqBJdl1BIBj!+zPVGFGw3nt*Ra#{OZsH>=fB(+`Z4V5`f@tm#rzJe z#=fq;z3n51BlyMu^zrgL+cNKR1YOxG>UPK;(rxY`jeVdVXI`u{_kqLpJ!E?8X6A+D zU5A`iAjht3=Ed-?`^bg@Lneot^U6cq+DA6z8}{)ON7l*GHRWyASPiu4go5ul z>Va9{=}o`=@92xzFtg~-X?5D^^gv(yA!qx74`lhGnP)kINvt*1LwvzGtog0zxVpZ{ zb#u>PZDNf`CeM$JGe_1t^8Yi@xARRmemT~FD_H{`T(R|ut*il--)`=eBf8e5PQ~h` z2`z*LSDs&20AA&;ohWOjtjq5Htj)Ub>YM4EN>lHYqIdkVe*SmrvrbC=7TTA-i9Wjm zn_!ccZ2HfLqq`*DK)Z(SOV)j0QdXyddjDMEeza#=T!C%rR*WT z!(L(vdy2{IPbURihL*6`D6aY2MMX92J@mb5tR#N#-!DlfErq;3>_hsp59vppRO+4q zjx)iP#y%uHh!I@>mD3qI!Cs_!&`+lwV?Poc^z5_;20b&aX3*2qE+6#NwDSjTo|ZDG zVVdlr{Tt$DxbV?Lu2Ipk9$zqq4~acX8TLUma1ORT`%j4n*eh_Yce~$~hCe4-TQDc5 zvY;AYXibjVH$SK14>dWpFV4?Nim_M7tH^odq$PLA^3oTH`A z&rz)DZnZ3@!8@1qYCT=jl1NL3CYP4JrzJ~8uUX+9Sbewe(ZW7~R_z39)HP~nc&LiT zRwe{vmoSIpYNt>0s10rCq*{lsJ!Olk!Drc?tJKbBmQmx`AC{CI^cUMcnsM~cX#&hM6S_4Y$cs< z)d}jmEo0SAtAF;sT_k;rQ|*0Rs$KZXdjuW*v-d3zy8LF}T8{PHx0UoQoxUBVZvpz& zNV!`2_Gjp^_wCw#X?54oza72wEgri6be#&mV;i$yn|<0tD6$Z_#f}vfYrG-ugDg*0<_I->h#g%1Pgz`v?2>JIY&q6PhiVe3<2V_kR_av=;w9#AOF$wy0!TAW&K+T-Y?;og|}lJWM7|Tvr(*l{S^I7>i_(ip0&51 z%-T4cGn@bF{yy-Wq0dBCv&2u0{yK~C+~_pU2xR;%xeb2&&BivLax%6fdl}nboFccz zdBOFRfBbqCe3bJQ8C&bTAln*S#YzS4j0 zllM5rE$@7@F)pQ?jPWmWde%+7jq#r;Zy#f!#Tw^-B2U)W0?oCjut*(_UgX=k|20+i z74qQD^92i6ID+}&CqOn`XAy5Feuj#+pvwy5)vYJaQPFcVRaq>TvpV!z(hA?BpCZqF z*kj*@{dOJpTZ_{@yagY5#OWFyIafvHY$=ipKK$ts>b3MU)~0Y8eh_fvf;Tx9<4C9W zZ7J7+y*O}%`w6AV{XdA&Ci)4 z=dOg$l70(*0NK~S?yAWd`4eBwCz6goKpQW5qcTU%mG2?161bRj<1Do%X9(veui)o% z5>|n)IKL{VM29PLX6Db#c{QK&Xwr(ZIfvygLb$JL4?Jg_11gs}zb{bH6X*D%i&={z z=iEo zSufP`B`36TtPiYzx#&X~zu}xADtxPMcqV5)$0Q{{b^ zdG49cx#v-Q1hN(fC@1p%)nG&3uf?|LO~%iHSEa1ak8=kia1pv}@uXmM?v#RMM`QZ= z=Rt5ZK$Dzl+VOl!A3i5tcqr$vGG-!o>72!8VEed`v;OGl$XKP1pOlliva;{m)~7i4 zbuC}qrs5|DT}MvpXF-ud&VmkT^=GFsuVl`U?at^tibI4eNbq<`(S&*FR zW`*`s#?@~?UE6|$U<8|V%yHpEX6Px>e@B0_>7%{=aAfW%>EXmF!N$+J_T|W#oO|Lg zh(NR0nn}pNFVCOD`CVE627HOnnfCb;od4gePUJKHkH@L>rOf|p#fOJ3@Z!g`%E`KC zl@t4O415w^9EBH_En3jnl^1e0E_SG#ji(_a7j0k%o*tK{_O>wh3bjoA z9PCTZa0K9Bm+FFAX$(pJH%9MX2s-hS$a&GX_rXg4DF z7T}pxw$|2u5^c|*?N9Z#lk_vZ?rn>VcE>x*#M{gnxMJVHy$4;FSYxQPrimKov!AB< zqBGF*ayLQj4pmz_K{IXYscu6rTXsrmyo!2g%d$-i#;u4G>|d6@jk`~xy`3s z1#;_#2C=EH!Ir+g7klI)%9JA4EzmXoO2&vig0)s183eBw7sa^v;N#zkf6O_m@KWS{ zY&7O04 zGUV;iUiuTGjP&PmKu2KhdwvQJ?kv>V)NOP zbv<7(P?gDAAaiFR`EqaN2F@xaU-WxngZ5Rxl6>0 zqvl>-aC#pg}$~yjiXp0~p z|6^^t)!}N|mOfr8Z8PTo(``?sZTuwZ5ymlaB=$eUbNXFyC{|WY+kX5bk8&5O`YFq% ztjH0crD0Q6nl>f(7dacUY)Y{wjr)An*p!@okFUU%tT1iLYW*xo&WhwtgxIfnte?Ay zI|$`IYaVuFu?|=0_ZDBlrc|VPvN>l;<}L;IYb*MxcI$rGs!wBO72JQz!nTq-O@l;- znRl9kb5Fg~w8^~FWVfv%a}3++nSsWgrV0H9)jiS6ou=PWCLbNQe35adOl(=vV~_n+ zY>KSVgXE9Gc98o@c087CCG#hOz7{*FIoQ><_%Icnj*NHYH+Ie6C;7->NB(nN^WT(w zWThkjr(N@3Bfn|5FWSWXu7Cy!E1~yG)(h6@9@l}XGrvY}U1RFkN<(j51HP|Y2Ubv~ zH@!8Hyz4)RZQqFfk?X|o$@=eRJ@&HhyRnU8tf@`ROYDt$?@+NhvUQuo#XQ+bpET0t zj-UA0ZU&dj*Z*|`PgswOB?YBCwGUuD> zJZnE`zel^a_|*9-b4Kh?yUr`6oa`rW8DPv;x!d|xFLql!Wn@2T-J><;kT1GX*Lm5Y ztH~FgXWfyt@O0=ro1M|qoQs%t#UH$*p&gf~*6`b?U=3L}#W6nLmJy%&H zPdDebQ|5GYZZY|%opab3Vw+gDfdA&Mw!rSbU2P5j?5_DcCEsRi_-A#^e?jtXwuZl| zYkrXY?zYAhx!;SeF-bz~kg#coNLZ-nWrZe62<@_`)xO{!FZPhYCyxR{kHB9EOxpst z0gL}Aa)~c`48PQ4SNVb!j;uZDCrTc5;YZojH_%#*zd=9G8l+x47yPc2^t!Oq@CE${ z970nr=~8#4BdgVyX!yOvrm*|F#3$my##{-`A^5vm?pWNGf0pUr#P4DEYZc)`viw@d z@#U17Hfce|slJ?Yn@zgi^!GJ$A7AuyEZgvJ4Zwf(yKK8@nH|5(?>1z`WUpGE!6DkV0KO+y?eyU ze_)Om%!O3{1R{mu%EJ=AE)0;cn0At!hU+0 zL7^8Y!*`K%8h-}t2EK`W8IKzX97rbZk2wr)IK=ATC zsq{WPf|svCrOUkr!7DmN@QQ2*Ub%Bep^~`0(Uf|+)d7f7QP|in;`hU?;P2Y zhQDsiENPmaTW9(aMd3yao7K4sJs5Jue$0#ak}gG(+}PY+z{E->}|}3 z&e*TL5edtiz4`sy#)i)=YTtiBk^X+l{#-)N?dX&HnXjg*4qY~%x z6Tb&%NxH!oczvKMLFdcQ;7era;QF-93H3521vfiqq!+3Bd`+v5+!+_muT=Gs=bX{l z^G;;KcxR{&^u%_>MY#`lW@L?fr_>8lPs$cZJ<7LGS7dUs;BPI7`unNvk>^UHu2ny1 zYoR@l=Bjgfd^W7{@NQLbjf&^CwTqmKlpa@-?n@n z*!t$@jg$ZWdE`~@P7O^k~CxTyV>oPUGFaM>4pF9%<#yi{feTSC69UG|Af3R zebJvtT;_+wokh_f6CXp_$fkYX7I-53m2vgF5O}!S4S$bqmpnK8Z8_lD zDPv!H2t3rgU+^LWE{Xr%7hNN9WTE-*#6RF1ocNtW(=U9{yChHeeU!Xcxrahr_j&hYaDiqJ^D4{u$uDUWE@<~@!rPaQg9qEanv`; zyxtV=b&ayk&a#)jC+GlQ`wm`}0PVC-SSmHmRTH7XiL z=j})LiC#N^emjWn`wO})T#^VcL9VM2EcU^-CCg#oH)h`K#Ut;qJ4|b(A_Yk-ZzS;bm^Q z@QI6VQp!~q0oNMZ6TP;70&5iUyNK&$2J@~r^< zEsU=4_0pD~veMQnuhCW|@!r~MH1W&1<=rN{l5n%f))w@g)k|AfQdZi^8)LLJg?Mjm zJ!s;WwyI2c7UAx}ojiPnekJwN)`gUnwsv1=v~@1=-r8DZ;+M8cO?WC{Tve*IiZ(* zw)9!Q!-WrjXf#%q*Uz*)>6_TMk3vgN{QfuU3vc!hHTb=~hxV=c@s-}D*pr!|y(Y{> znGJaLN$KCatG_-AlT( zl{eUE>u%z`we>6VEIE3@gr6i7+4&0oEsZyLDfoUqz&YdkIHSK7ErtC8t*eo>ag2km z!#Il(JTjiEk*Cq{!${B4DaTrP+UA|I{aTr9rdz;pnS4lhie#M|o;PIqs#~sGr%9Q1BFKxNZwpR5q z+R{$bR;jJ6923r!kT&w*{kgzB$<^Pv*CqX_OfuSt5bv!&7uec}H{k@r-A5ee{xKuu zKsRmar9b=3wzNc}t>2!eEjHyHy!(`L!n;;8l>U4P+*5!4jdXo{oJJeliTBo@uS|UU z{5RpZ3FZBAY3~Ezp4xbdbZKL^GTL~8cyDcdOr9kJJ4|?|8CKGdUjg^j##+*)jm_;2 zeVt!Tytg*qF!4!0HkojP8Df7=-vZoI8w*L7HdY-s+NdVpTN~R%1>o^)xW@=K$QtBLp4#$zTvX`{x3=M#2cYwsqlr~XVJUE0!)8f{%nytlUQBF`F^ z=_Y&)p|#d}@D1D7+NFmbGcrFg=56=2R{AJ3i;nwSve{-&Yb|T$6w3CrX7&*p%r!F+ zxI5j}nmJ9%>G#9<2961zMc94KJis?9HowKmn+>@+O*U0q;E~_guq+I_;9y>s{_A- zBTj9l>_qX$xSEtV_ci4$K)1IPI=m4lHh!Va#%I5CfP1dcF#TcrA1Qpr8<0KordPbp z67DMTCb2(nS?r6Jvwu97z2i#u=W-8n=&wzE&VIn`hrT~w3!d6s`axZ(C*PKeEzqK| zxjgI@E(2fd0t4R|@X3DTcEaWXd|ReSy--|?&+q$&udEH(IJPO-dzU!e|Thf>fcb!GL%#pmk z#vI8Z-hGaUT+1AuO`iBfg!h-2@B~6RTNnE)^OugR1(!U(_?=_Mj2%WlELzIN{^EO6 z-|Miy`s;WM8|#iX*BRU;3jTug`#&x?ao*?FR17b2)>Rd>8l( zehdixS+A2`{}JJF&KyF>&)4at(bWGYdP(=2m^5~uCx>V+U*xY`GyDiotL1NOHOAT_ ze$ceqiMsxYDQ~1oEkBWk{=w&O_=Gil3`*BO3&0bJHF@(r;vd-U^-K6sk)eOa_Mm^B zHuVpOwYr!z*X?DGHgnZFEc|Wz*B2~W{8N7jXeOyaa{pSw)yQ#^Tm&v~9 zR`PZ~WA=%$m9vX1WJA{DWB+8N%iMWyyJJS=oay@hN%kzooRb-Q7TFJEhvtLFT9*?= zwjPytO_raSO}gw0y3TP`KSHNT>HC$lfqRxURbM0JkgNTCr>?|=ClShC;zMLD3e96! z3)`MAiq4R=a6yr|7G~;ep{^$=`xtdZZ^xe3qIwy9EzIn)7P_2S##)HJl6{ZIwicFA zU-;1c&yE>xeB7VX?)9q@>RtG{TUY}kmpB-U%+L)^qg-ICV}_O|d%EQMJe}^W&=f+; z=BYLPlA`4Q)xuW(kQmfzM^!P%u%zD;y8x>NiT9?nDM%vH|Cq`irPle{lY z95Hb4(==t7%WrSH4IKKNf|G>Y8Mu(LM+DbLyS)Loth29rn+c^~uAtSg^!n?pe$`<& z+$cW(j1a#oW%U22R=?8gTRDFc8d9Zy=tKQ%btT_akvKjNJzl4x9Q17nBM*A~F%^~a z(q_T4jvmho+-6((Pn&Jwv;B9yt+db=61JMS{;s0!66gGA4IyV?(;L^!XcOFTP}e%s z)i__~tgpWP&MBP#P39bM66b*>oC_9nu3N;p?yu*`ce27+ye|eG(F2i>SQCf(qJJTN z;w0-O^Y>%PgJ*$b`THh=U(&mFu8snp`2HSnQzs=a!w$=(kfs%2bFV6 zu?>ZHhe^9wzA?hQJjv6}W8`IrJ|*vHN7<~f0FuA{v`ayc3Dk*UCxl_Kbz$w9&V3q|C3eE=HtFi+QY})IFRRhCOcETN&wx`uJJr)KFN$7Dy3EZT z{mlFreB%6OoRsJMCIbH?e;fI72f>h=^w20|Oy)>_pmToxl`?X+H(0N`Pw$JI!<;4Y zU7W>9+nm1ygr*O;Hz92^55wFY*$8g^tVY`ZJ@?9`4rQCAeafga)%FLi_KEwrM^Z<8 z7?0d#7|(M7&-rG21mU?na<`#?C!a^|Hk{3K7LU{$$}@x~k4Nq{$i08L+c20%A1~Ip z<^*U>=Ds<7yQm|5k2fIc$n+>PP3~4$b8J z{t9uS&!SDz`MoneKg$e-?~N-gU7OsgYt1#p$C~kmW{i*jPD|G&*RKb+=!+K#MStt_ z9G&z_i8B`?GEc?+Cmu&z4j!vd(x&SUtIbpOyR^B*Y_s`xtIhQK#b%ouC6xL4F?aUS z?+p?c{Vs9z`;!v4>31nFZ5G_tvwp{LGURVbr~FM39LV1a36T+ftw6uuD{=IDvYww2 zYUi9G_G%*EuI`Wi9)KPnh&~^LULVZ8jvVfFB-a0RlM4RDsiJMpxKPymj>vJoA@UXD zv1{eH;4Wu>eT_JFt7C>f4;DCPd@)q3|5t~j{ym3QzkMcmidffU&T@9=6 zozSrQz6mbwe{g=iQ_Ajl4h+r3ZW21eJe!@4y&q%qEO0Tt=iryc9%*p6>SS*<3cQjp z^?qQhv!6RdU&0@&d?V#ut9;%db@+bA9C%um!Fj#Rr7-P@P36Fr_yu)be9zgzy##-` zzJ`o-@|{JIA-+o#4ZkXO+Q863o_8qQxcHCbME)9_Guj5rxb$n9*x=kaA;LP-2yeqx zN~ix*bp`F&x^>k$D&L1Ez-$!n4LVj*Uj&8sX z_+AjX((iLwbTx{ugRXFwR}TT$_L(-@)@B&a7GX)Mh4jj2Y8oMWU*Fc%e=esAh_h- zR`U;%egn@7KYD!sBxC>Z@8s7O^ZqC6w)hMJk7&BgcmGbQ_yWX!kPw>nP};v& z89a(iH+jUkUja^;kGJY!X6Q?3JbKl?qXYQ%ZHe*STb-}e5qn*DkPdGwJ}7v!z{Hz> zwS_k=G)Bjl9+ER*!QXyo5wh}%{vFkx@Vlf9-B)ztPciYAd<}kU zzm&!qqV&TB4(UfVW3yzn(T<<-sTiW_xlgFeP-JRnJ8@<^7X9b+Rp~v^UrL?c_=BF% z@#cT06Yqb*XPI|t-RK|K4gV?hPwAb0=8{2wNoPCz%yy!B*f|aw?Z~>;UH5A_#M@dogRRtz!j zOMguJ{wZeK7$faJI%w+r&Ms%Ux;~J#SM-9kbzCo(5egFu9{K*c%u%7|1wAjFw&gp| zyMpf#e~m}rb|>#u@h#@|J9jt!3cU7llW+4%Uv>WC1Cl!kA|Ki~d}XwG4LB;8YfCPh z;sg$}E7qVWDSJ9tH+hh1pE^r%s%M-GRJKJX-h9u(OqY& z!v|OkZ-Ylik&k!j!xHHuZMzshKXm$u%bhd%4!)#IS$s01Ek4K{R?$rr-DuB@#V(e5 zyN0O44X*?q7rwAD-M$3ANZokq>0@N$lc8riWl8wjR{k>UlZNi$JKlFX`ArU0FS4V2 ztj7szyYNAD{t^9~Y{;@H3xP|0(ExrB@q5Vk31ZxXN@mOh4x^07h%OJBx_DP#b?Xp( z1+~3UM|T}Iv7)o)B#wfZ*`R zZ2gS+ag2G|ea!VfGMB~QX?>g8F8@tRUEDyMtBQS5?3wBL(smDY3(eSVkN9=EyYjnX zhyyz%J>*EOojy$J_w{Y9=nJj)6>VLe)HJc;K5gsFq}qwA93@+2EnJ=1#N6S1P=|Bt z%*@(}-Pf|1)4BK9^=^2R#P}sMehz(2cIBduw*rs;4GFvvZH!yKXh?)9DeX8b0TzMf!dI(cg7(u+J_;E2biY+TV5?bibo*}^h1Wnk)qlGVG zzrBdegcE)6wrFeB>beOtS3fYJa`oB?6{`ai!imMCmyrH2=?{_qAnEHSh+gwQ6&GCj zkhXK>NiDdNZ@7UEKh*d2J;vI8oqnkO^iSqr`<<5~7nd<-rt&R+r*3oJ3@u}k{g_ij zZo1BUk&^)55D+wmO z;XCiI<41V$(?V5W%J&8=9OB=uKAAg4#$YNuF4TAr{@JyU7qk~WT1WroI|rjj)jcja zE0Y;B%4vK9@niab1$`B|o8h;g_63K0AJG9GkCt9%$<;-aQ_0K;Bfkz`SDncA7s&AD zhbFZJ`M%~5`8^EyDl$2S^3|b9K^bfL{-unK(55tZUB%g|WBp;w%c_xcDd&b4d>dMD z@8SEvQGeg>=-+A-U+a)jfyZSIR-&US&{>-1*jkC+sz7gPniIGg_zhqM{2Smsz}i-s z6LHXj>~BYwLn5!Wz?ICinP*Q54rMH=q%HatPk*b&`=SZZpvy7mO`jmg2hqR!98FQT zZiZfc{Tz<0Azxx6w=fSaz6h;NGS7d?H(oWL*Wb6Ot%dQD`BF*Um3J#$F73J+In(7% z#@Swezxgd6)~pR8cP`cvx2H(=JwHVoh3Je1Lbqp%cdhb_bssJtSE%L9-Az00c-}cr zQBm!1?zk{(ldNNcUqTmoF7k*^mw4<|ZEtVyWk|b`6p?phPvaUn1^<(0?54Z5;3niZ zz}V$-9{e=x+S<#ec<~pH-pE*sk7eUS#lauv_@XNl)TmngzE3Ol;(?X?c1oP8U-OV7 z7+}sUR?52?{pduV*YJ($rOMH|3S8npl5-DvH)|8~YZGOK*Bh~cTv~QrfBLZTV+Z&Z zd|@BM99qLS(50>2e78z|4@Yp-!k++iH6q_@&^KY$YC}(Egzo2FK(5MY-7mkBVB(Vc zVdtX*#OHM{{rC&#IhpoUrnLuZq^#T813hK#f!<{ggpc6FiP&R;FR))T&RxXD#y4WW2Rm`v0SM`rowZ|Ddy-Rc1Roy0yca z-DOXpkVQFXmj190H~9L#$=7!^_62F7J>*+^k5^2*!wIeZM&++G{hr)Z(u60af7Ckm zQL+z_IeTKkOZZAO>^ui+my$JsY*0#QwRPvV{fu_sM>p1iQ4-m?H=J=zS^0! zC!f2!xf$xP+k+1I59i_E@~w4b^>EU&fyXcV;_i1YJG|nN+Z{=D!#`Owy^uBIdF}+c z`-l!xhx56o;fk(&$3>attm(&jq)+?tDarTB4&Yxp$hv|I>t}po>o)Kdp(90yvWAQ< z1rCxn6~ADk104eFBW*vSl<_2TcB|*J9+F>5IVlt7$usd(xYX!;;1<#f$d8zKJfu~U zCU|PW6`ai81RkNEB{ zZ?o$7VlkcOsr)vK&@B8HdMzAPqMvxzxCmTB!6oA`l}E z{UN+^ktbv7UX$A<@}k%0_t&C~Q5JL6!mO#(gqDydLI1JWg8v}k7w?Wl|{f; zpw9(E#{Uy~o2VmxV&TIVCu5I&#lE1Ovr`F|eUvz65wuJVXsM!K9#AUv+-m>0%iPIz zJ}viMX{&*!l1JL&o~|wrLbu4nOzKJ*^09G@_bxi>UaBzeQiP-^uT7T{Q_^R)S6|MzlM(UXr}@?4T_Uc?zA z`L*<=nJ@r9ck{>@pXfn9X>ag|d^Gb&{R$onX9e{FdcH5ZmV7^t_^G9y;Hot1NSo5v zu=Ew!j#KoX=sVGSy1pZncEUWn&3d)a(hS_fW9bGdBY3U8C{rifaZ8yB9*bTNV-g@< z&12Ck`c>Dr=o-C!^se0bvD@;N-{tr%|AriE$g=~v#%FXV>oaz~{WrZz6@88C)FaCBLh#lRD=3WV#Jq=Rr?4 zk}k5nd3A1^wA;`_eX+%4jCH&wj-Oll68lK@l7gd>F`h~NN_alASk*txS}Hn6be@~F zN$>~A|2qA668$f2H3-d?e#Cy#^&@HQ*`kfmD`mRqM|9|E^&|STH~lDbP{Def4}IY{ z`b0l6yoORKE(l0stub|#(@pr2ydroPmLvPYIT^=pHNu1v!(A$?b#Q!Cs z*n@uNgjG-K$efV2WlmJcywc}HMyrp#qRa{QrqksdGQxPsyx<=C8HGdDtqshX&N-6N z+CKEH{9Es++wLf+>b>ht1z(si?(J>G`>{~mkY zKPLHhzW)f{z)n=7G}+@;$y%Km%Ak(e7oUd#2 z`*Lsz-9sdV_S5Z=a~ zmqWE+IFWBXQ_i^GDc>N>3h~{a4TW>Zjw|IJLJNKf_J6lVD3iw?uqz!|Ivp~OQZ}Ht z$LPe5fgA{L!_++nKO}88eR$WAu?PRg;d8t7$NEN%(WlJTO8i)x`YPQwDBq;{x^WbK zO8feK9NPRNZAzVM2xBGM-k#QD`?_Jtx%i`3Q%7W-b#ingXX(jE~p#jNv0k|BElKUVL#b{DE?3L6`pmW37xpGye8b9(cuf zmCr2HM$fH*XDhUfWqUZ$xFhv%uC zGS~bv*Xn(8HaqrV%8IV}3-bR6u!lR9G7f9ac8erfD5_gnoogwH@88-x0d8=Q~1^p~Z!tEW#b3K$WUr`*hbk zLdVQILL(?Ay6Uzmx*R#j&`H)Ep?vUKcZ9CB)sg+# zf5IQDd~bJz#MfZ_HhBM#+!6X2bzJ5hAyeKv?g*72_r@I|@xu%V{f4skdn3ubQt+ZX-y;WQv%@yOhebF$v<2Yplef34G!9(s;Yr~i^p`k&PE z(n9}2Xs7=rJv}2dn$Ye?mAN4}{Np?QobJt9z1t`HfirpXrh4)Bx!B8dcS6=7eA{iV zmyLUd73ce+PwTkTLu&}7eDj!2KjkvgI2*2iNH5dr50(93Pjuu{SNyo_3!>fWa9@0S zI#!u+5R&;y;q z7a7;ahg6&EE`u+%@a6k@TN$B=COxk6Eqc;J7k1FmdDc9+7g_@68MGMf8?+d7b<$EW zqIX(kjK-U^3@~YFEU;)v3-#0U(nDwUK+}jGXwugU^seQ559Ay4$^Jz6bd)u2ZC_Uh zpQfIp&5*MxgHLJjNo(`}d#C)QnY0~=Gv#MkCv6eZ3_pD*pT6QZAiCcJ=uT5)=;>(+ z4mmwd|Eaf|5&E4;Q@*E@rdLgx{-}o{)19=jw>X7188?224O=EX^dNOs^1g+Crf+9J z%O10qwKT|DUcjBowXAKQQpSm%n8R5h=aY5Sl$|qN3s!TsRhr-#6QF!E{oFXtT^js zKd?2WNJXWOQ{iJWGLS-FuAuF?8sF8?lJ~rVKcw8nT?qIsICrzxs&ZtHO-c7fuhdfZ z2;Qxl_WW$}BiJ;_jKk`av5Wf zDtvY6d9GC9E598H4OV%DT`D|N*7zVd1J6faMD&A!hTUzF$%we(PyUts_<3PD3luL5<-!qlRDsRo- zrOaQMlRG#cTgsX2vLVi32t6|Ad~`AIt&03!cJao?B|aCL&V~Pb_@?VT=&1R@Zts=I z;yPrpPD`$vcfJ;!gPzg96;EHQnXBne<*m)>eb2s%Jimjn=MIrEKTqA;`;dvBN?&<| zAAdu3Vr%4{el~uz{`k`d;8z=ne{B$cN$%?NK9tcYc~*HD*Pia^8}_%n_n58Q-@n5iw(ReG zW53y9e-9P=dyV)@`eT3h4}D5myInY%_QJ{L9sPkmXGYwSUd{r3U<90xhC zxDH>3e&1_l=iYV*cr0Gd5&IiHehMFT{(WNbPvW>>O-S`S06s0piZ&H zt#@{qxVmjFK0E1KZ+(e`d+*DBo!<1&ZbIox;Ln|X`MI8#7J8G=&i{RSdPZm(;n&Rx zIV%yLOqBWXCG$c28`gY~c_8vH?GE@*>Gvi||EOl%Xw3)afY`0E1#!U$a^AxD2ybOR zY^H5%K9p9e`l;RK!&l~f2vAP&fAGGV;V=2R`H&s@_x9LhLSN^6$O^skRp)%D7PUn2jXJ`MkE@?R!68H3S;!rRE3 zoibCP=cR=T3GFm~M^FF9<=;!aZ@L%#rt}ZfX-N;o5$g26-bsHH--5_kKl1GKXX$k^ zLh`OxtVGsr^IAnO2YfisbPl@mEP?la51<*$&ZT`Rr3r!}B$u_^CA zeof6751*y|#@%Yhk}CASCGT@jmG{5Ut|jmJlr88c@5fDff0J^8^9FF1GJcl+uR#A- zoTC4;LJnkC=+ot$J8W|IHNhtD@~te9_h#^JCSS(=3)Tm#j_C7;$g||#-;2D9oniQ% zGj;txP5NWXd&IOStar2y;GZ+}KfV+ECvzy+Euy8ikg@}v6+t@X&#^>>*f?QiJZ#DMotBN06`G2?G7qs~k&!Ax#;-h%PRA|~Uu+ngMP$4mcqL!j zm_QxQlJz>`Q~v|Ju*x@6J_CQQMPri45cCNxLhGlT1w?ikzJl-tv;khh83BjRho{Xk z0RMlGCi>jZ@L)%vpho)T)*PO6((+%x+UJviDV+^ikEXt6jYq)Ej&>=^? z`RBeVjUVZBq=$A93XP@z)JfxZJufZv8==WYb;Ti6&grPzC_IG%W(1fE157mu4KX{)MD-s+mbcbxfNemd{CsTs8sjrpA(dd$SV z`A2Goq)jF+b4uVxfipPQx5}P-?G8m zZ18#;ywC<;ZG%VJ;8YuY=ph?^8@$a1ueZSqZSd7Lc%%(ZwZVrTwBfhG+idW98@$j4 zUu}a&+Tc_hd}y5wzYX4IgV)>Og*NzV8$8kmr`q5{bvFDqc$*DgZ-W=w;Hz!$NE@7L zgAYAm!*7GP+2Hjyc%cox+6Irb!KpU*&{`XQ8@$a1ueZSqZSd7Lc%%(ZwZVt@u7V}2 zZ`t5&Hh8@aUTA}_w!tH9aHz6JBX2_FKUz(Jgpk-hs3Cj2RI$b{v*Sz}U4nS;RCs5&fusfSGX&%pm? z!XE%fO!$4^6r7$?_Yc4$Ojv$z^>;bp-1=k+p6f%)c{4le?J+=LeZ|H6dl z0sq#7tAP;{y*y`)S^LfOD&XTLTmhWOB|@R$`@k6{d;{<>6TTjp-+^~@yoWdFFPNvfg$Mnrtaz;zmt*9Ou54TDnw*XP6qXg&wyS?$Ob_~l%I z!ybVTdJZ9E`_74LoTT1AXPHuO<@+~GQX7X$=99xkx|&DD?SRfvQt;(*aoZ=Urg3)x zdi@R5AD_Q%lG<^JP!O3A_wpq5*#rpNH_=7Hn}uYezYBA#4T&$xkcCDMFJ}j2rM~dTs z4;5zvuPYe`OPfoE-9t)ADk-&-Skz-rG{tK_k5`AF6sj78lZkoP&Eyf^x7OxbD`RjBk2coJg?DiuE5Q?o~A-|REM{CJTDij_dSxn zKF_nGP`#Wd@XA;^26RI|z z<=IlGUOP)l9y-hQzE^EHdnE8{XS;TcSL@HoUOyh$9<J1+HeV z+H%2Yimo4-{h3$o8#!o;S9^1mtF2J&8#PR+wo&PT>n==*cva(t>5a(1g%j2nsfLRN zHN%aI(t%H0>&Y^A5(g8oz;sOpWEO=jY#J%NE&5jf_fFmyK(B2o&#HiDAh6)o6Ybl>O)MuJ# zJN{x%jMDPe1|l1@xJHMzO_O$B)*4AU=!m1}AxGQ^2cP)u#3uY$k^(0hG#8wBQvVx@gN|#Rj_aOy z!NvUsaDg3Ma2Q+%1s85qqs~GZgC{i?EjMe@+CFU(6h$1VpT((-&Xh=;dd8XlR-D?= zF&-Dp6FN6)E=s?qNr|_#NnknXNCo!>X9}Z-a%zrK+dAqDo|UrQsn%<*El%~amhy~~ z)+W8|RO=j13MIn70xegeyqBfIy=SykoSfU?DDsiYhUPQP*V$Ljhx9*^nEwjga&c9$u^?c8Sss;0JuD<0w>FcFQ>Z*lHsu#>(JU&U?vZ#8Np6Z$C znIO6ID{G{Lp1!bhQKh6WtErxU=>-#JE#u#^no0?ln@P)cn4}~{M)i5K7cE+_Nb1~H zIdAD~&!wJ`RbFX7c4rJlrF^!kh8?isn#K^vWS3CZO;$?E)?#eOyH%GrzkQx+_$sa&Ei zJpV$^NY6MmZ1VhBoyi_GY}$fZvllA>G7Naw>8ZFr)*4p4){3cUiG=q1cG&tT8}6Rh_RFQl9fhbujoYIBjMEsiB^! z^B2!vBxCKlQbuR8j82m35Vi?rmyREA#%AX)uD)Y-`4Z3k0?#?~&+*J&RyBL}ti>MT zQO(jC&*G&E7cN+|WcDnNlvB5>MKTdZ7Y*V;2wbo1ek^^Ac&#QJ%0dq|`#f zx3-|w387V`sC>0WB>*W$t)|>F+{!sGS@k=VrTWequlh{5Nol@vl~7gzj+x5go2e2~ zs#HpLm5LuxrBwbb<@C%}iGyxd&S&Prnt94~s7A$S&sWN~KqUtjsQB#*l>7aqDzWT# zztBkYc9nOV+dmRt6!QEt5%&o}dTG#Tk_ z79O))cYNLB4Rpx^1}6q}z?K;$%$!vqG!fL9^U5X1TYTjPgDckKYX2%=GT< z+1s|mRy$T`&$Gj+X5>R;iLTVz^8`*c`ICYildmFXc)|?ZmK*7Z?l8j7?li(}Y-OZ; z_rJ*9UDN+pd)FS{RFUsb(t@ZJil7!Q%0oa<)TS@m3m>h}00kCnTiE*A(=2$dzlokDTHJwvLJOf+J-T-0wu*^aUjC0k!`-NE10dQ=&l7yMoA*CeK2zTiUe{Wy-p)W>Oi zTTnW&96ZlaK5Dt#Cuq6-D1Df(O0`_+T3YTqlqopQA--4Pr3CBsp=?CirnU!VKa`cI zmy43^2h%M-Wc%UIV7{pBS0$6}hd%$euD_^uR2iwX+P+Bn>T!6kqkL4U9%rR`-m94N zu9oxeq5QY5r!sj1m1@2!i}ASk$?*AVxq2r5ID^SLL#_k0T}=_%rP||?Gr$?(3~&ZG z1DpZQ0B3+Rz!~5Sa0WO7oB_@NXMi)n8Q=_X1~>zp0nPwtfHS}u;0$mEI0Kvk&H!hC zGr$?(3~&ZG1DpZQ0B3+Rz!~5Sa0WO7oB_@NXMi)n8Q=_X1~>zp0nPwt;Qto`-apd) zsM%Yo%tKj#@(z@hC_N~BC?7$&2IVG{O(=Ju+>7!9l*dtig|ZW6;x^nLi*gvsaVRrT z=AxX9vIJ#0N+-(2DCr%;|n*$;k0gX~8hjB*6Z?%L1uIRl&l z&H!hCGr$?(3~&ZG1DpZQ0B3+Rz!~5Sa0WO7oB_@NXMi)n8Q=_X1~>zp0nPwtfHS}u z;0$mEI0Kvk&H!hCGr$?(3~&ZG1DpZQ0B3+Rz!~5Sa0WO7oB_@NXMi)n8Q=_X1~>zp z0nPwtpm!M0qS+7w3xNZF7n|hBT{WukciglFQ|AuORO%XYcvUgc-$pd{i7MQJq zueC?5%JRA^B#)h%YUy#qyl!VTu{Mjetf{Bwa@>`U1vZntRyma&o$8=!(2nm6GhD%N zep5%kZ})h83snC$(4Vq7k352rEzoV_lC+N>b=b6rZsOD-vrg zmh*SqKiKJLVd#(2RJb!)ZYrn|-0k`Xb(+FuqTC9{*@JeIE%exekYi$n*DE?Dm)l`C zS!GYHSd%Ka8g%u8z)2BdYBE=e=E^FuqT1vthml=gYb}Sv)n#ItY%iC`r&?&+RkVGhNbJ+7|ml=#klIN*(!3YTU z$-!BI#^q|Qz_J{C{-QQBLt_TLx|~XrTeO(Wl7|erRFOIag1aDW;c94gn8GER!r8LJ zt#mL8MjnhG+KEzH2G%A6xFc{+q1<##O(U)n>13svulNvgozN_k(<#9zMT^c}wWyp= zZ!eg1Np5)|3=*Z+K3&f(<~Xa}l|Z<;%Ba|^>jj%|Ef>xz%2clQ`d+kWvw_(R#q641 zFl#CbG_yf#FDgF*?j!6TS1svnQMy=0#$jU|X7bYL7;lP0*9CJuEM0=TC9HDLydOHJ z8PnlrgP!!(IJ7e@Rv0<7>3ieQ*Rn_Y&=$WHb!Kmm!-UHJ@i=s?ZlgkfM|`?6>GrAk zbnha%o$=`^(~(cFnLVQP-EgqhX_ceA%Yg}7S(EDH(pFE@8jUtCMRd+dRVIgqa0S#M z(x}PwFE{D1a=>!|mKOkGYcM)CF>2Wegmz4;;6+B*xsDKqm4=Of>|K zDqNCdjR}k+ZZXMLlF>vmz;Uz?&_Y1jLdoSIsFffKJdWm&AUQ?yHRMP-@81ui!{Je5uM3ajKU$X36C?OpC0m+Dwv$J1%JY;g>N4h+YdzR-!wr;5IV`99oNC1Hu9(O1_GzRpVr%jM4`Ds2c0aPv0FE|eZZ>3WkbHpYoEHT z?<##Md@2+7>yiea)phF5>UQf|5+CVv$Ek#IaD2vuO-X0$5zkmg^S?*jh&Vk3`Msa!--@^du?H~WZyDle5jQTU1 z|C(YK;=RA2_f;nIMcjh;AmZF7X#VhFwETX= zPQ;~aY5p$6-yy#Aa$3HA9nD{ixOF|nyAdaEpg8FYTD};u6ENZ5`y|ajg8B7`#}23E zk0Y)@JiCFG--Eb$BgOwYf|lRAiQ*R!?|6pdkt1n-<8u@*N8GlV;#R=@;ZmB=O!0e& zn_r~(1Y&Os#itQ(MtuG#s=pWUFvJ~*uSeX8cnadluTXn)5f>w_L|lh>5#k+)e~tJ9 z#2XR!dzIR2Mm!nuTZo;Ak0M@$_^*Jk=;ITPAs%%l#Yx+#y?KZyAzqHS3h`mYO^An$ z2L1DWLQ5;vzXfr>HzjCC15W=;HQxO*+9tY!|NDu#KAn4Bsez1oJbXVv@fS{=mNkr2zbSiei$be7eEJJxemKppkgTG|(z(lQnHiPFe*va6B7<`z)qxxv=+Zg;RgA@8{%NH?t9fOZC zc;Iwd#Ne$A?qKl9ep-8_3|`0JcNyIGe69W@20I!2I|d(R@PG@n_KXZ(#Ne$A zKFQ!TxR!#)kB!04Gx!Szk4w_(S1@=3gAX(Ke;ND(g9luswLhN0H!`?_!3_-F!QgfV z_v^2%Zy1A98GI9is~P+(gFj~Q-~psQpa1&3TET(qc@l#q1}|dp%MAXO!KpA+Y0fu; zpJVXn3?4I3t8Z3tp#4=0b~AV(gO?~c(BHKTeu2TSGx#urPcZl!2A>E1KJ*`ELHZfO z;0X*)V{i_GZ(?vhg9{k^O9qP!u3)f(!3!AdV{jdV>lti;ZdX96gd{_{3sMy%D6DF327Ik-H`S`dK1!HkoH1)8`3^V??BoQ=>Q}*q{kqwf%G^eaLI+IAi>9{ zrbhVkEF}2&x=PPO=AR(7LOKYE^BmwA=nzH-@M71H7Tt3KrD)fW97vIf4>w3=)Q>kv z;clO#K%rP34MFjL|Fa>eH1x3+gaw5>1{0M>WbE5+P@mJw@43U6?Cm$)VXA;fv^{-! zuD;t2V_kW$O^)}@`|KWDQirbkCOgPrT<@@h_2YVb9b}-Z_ts&Gg?nQiRLuYKAn>l5 zux*B{|GU>9(f3;!)B}c2T|wC!zX5 zor+S6Aw5bhhV>|4x2s|NL9L8ZGi(HZ9A?y@{n%B5gY@Ij4vNcluyFiEh4Nl%l-oQ7dEvfYRvT zZ3g5(bejy6l3@+TvjYHU+_-lSkk&8}>rtIIb$N|X8E50PM+j-e5r8oJ0pMt0 zCx8-wq&Y(O2GC$6j6zbwk3x?|fk-3KLPLy13ym!jEi|Y|jV2;e2_qt3 zNc@l*I25w`0BE9*Khg{k7!I-1;O#2mKMop@pI<&x3E+V{EqP9)Q1&Qe^Ua0DVxF_G z7(!=AwP+1bwZqb-pRL?Z!QG#M3A=_!ozQzih<=N?0(JUjnk9GLELDrW?xO3Is?F5d z4Hey7Gf`@*b>iG)k*!kXt(t5@qf?^(B-s=O@no<+(<=Id27}K}bjBD*I4Y&?9L%VU zMDP)#(w(SWrWT9Q3xm|}nkYxhWfG(rG)$v6=ybT@QJyUm8(g1A3L-TF!Iog!SDQC`{Vb)v3F z-OLB7_rQ1nv1wzSvm@-dA)7}i;o9x-Fh>Tm$fG_W#O-i-@(T;5l?BI0rin%p;%GG1 zz+4OueT#!9*613tMcDZ1h5%F8E~eBirbo6DuJ$yjp0Oa?$*v`< z@p@{+>UAw)D_7vv73VyVU`sk0NADMSyQcNQR_w7BK#xG@J|*- zvKJ+Ac}*l6mLMKAd`+Op>>$(hN&SrkwbmS^DrQM7%2I}efpj{HR;|1r>Doc~*ZBu< z(?ovXf%+v~`3Rw%XP2#>m{Kon%`HVw461%}uNhQZ$V@%G(Z*gh$*kjnyGHOjBRU|V zr#2l%0{!g=X84+6VuMS@sN83&lXtweumdP=JH6EFSp122+p5$NwTlsNEs?y3xa~!7 zDB`up`WW3$f&YRTu-i)i59<84SJ_K9oBB7+csv~#OAuy&&|JY^Y(r355LKH_gs13 z$(<`tSjOyby?yjknP-X~`!Z|RqM=7KALxJCdw*G0dpLKxSiSPq#j8sk$8=ZTcF&!A zR{rJH6T4PDzQr}_@2@94eXcn-vuR@b^6Q*0+AjKe*(ED?iEC_YJ{?@_8Ts{XL-voD zxbDmqu31M{NW1=P=!1@@OIBY!zvXw@co59+`6R$g)U z%-O?alK`P*|+ zPG59=#SrhJ2P^0PJbCLSzGWZgEVBxp-J^w zfBk&*HGfL&zxIP0ZqLj8hw*=Q7LFXd@p$zo&HHv_=AZoe%kk?^?k#wD>nGVa7g(nE zYe;+ByRmbyb=CM$KW%yOeZzz6j89IJ<|b^pwcn>}3lA4u@ZVGac4kJ_xMv@pIsEpH xHN&%hm@wbn^!nA~%+4p?v1SMjTOOEZn{ed4Q6;uKL*CzKzAWBl{r#Sw{|T0Q`85Cl 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