zhangyiss/libtin
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update to tin library

Browse Source
This commit is contained in:
张壹
2021-10-14 21:38:45 +08:00
parent f06d009807
commit 0b52bae14e
15 changed files with 528 additions and 714 deletions
Download Patch File Download Diff File Expand all files Collapse all files

68
src/CMakeLists.txt Normal file
View File

@@ -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)

94
src/demo/demo.cpp Normal file
View File

@@ -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<double> topo(10201);
std::ifstream infile("data/topo");
for (int i = 0; i < 10201; ++i)
{
infile >> topo[i];
}
infile.close();
std::vector<double> err_records;
std::vector<vertex2dc*> tin_vert;
std::vector<triangle*> 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 "<<std::endl;
outfile << "$Nodes" << std::endl << tin_vert.size() << std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->id + 1 << " " << std::setprecision(16)
<< tin_vert[i]->x << " " << tin_vert[i]->y << " " << tin_vert[i]->elev << std::endl;
}
outfile<<"$EndNodes"<<std::endl;
outfile << "$Elements" << std::endl << tin_ele.size() <<std::endl;
for (int i = 0; i < tin_ele.size(); i++)
{
outfile << i + 1 << " 2 0";
for (int j = 0; j < 3; j++)
{
outfile << " " << tin_ele[i]->vert[j]->id + 1;
}
outfile << std::endl;
}
outfile << "$EndElements"<< std::endl;
outfile<<"$NodeData"<<std::endl;
outfile<<1<<std::endl
<<"\"Topography (m)\"" <<std::endl
<< 1 <<std::endl<< 0.0 <<std::endl
<< 3 <<std::endl<< 0<<std::endl
<< 1 <<std::endl<< tin_vert.size() <<std::endl;
for (int i = 0; i < tin_vert.size(); i++)
{
outfile << tin_vert[i]->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;
}

727
src/lib/tin.cpp Normal file
View File

@@ -0,0 +1,727 @@
/**
* ___________ __
* /_ __/ _/ | / /
* / / / // |/ /
* / / _/ // /| /
* /_/ /___/_/ |_/
*
* 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.
*/
#include "tin.h"
// Start vertex definition
vertex2dc::vertex2dc()
{
x = y = elev = NAN;
id = 0;
}
vertex2dc::vertex2dc(double inx, double iny, double inelev, unsigned int inid)
{
set(inx, iny, inelev, inid);
}
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)
{
return true;
}
return false;
}
// 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)
{
return true;
}
return false;
}
// 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;
double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y;
cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s;
cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s;
cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here
return;
}
// End vertex definition
// Start DEM definition
dem_point::dem_point()
{
x = y = elev = NAN;
err = 0.0;
host = nullptr;
}
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
{
if (a->err > b->err) return true;
return false;
}
// End DEM definition
// Start triangle definition
triangle::triangle()
{
vert[0] = vert[1] = vert[2] = nullptr;
neigh[0] = neigh[1] = neigh[2] = nullptr;
}
triangle::triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
{
set(v0ptr, v1ptr, v2ptr);
}
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;
}
void triangle::set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
{
neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
return;
}
// 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;
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
*
* original
*
* /\
* / \
* / \
* / t \
* t_id-------\ t_id (0, 1 or 2)
* \--------/
* \ /
* \ n /
* \ /
* \/
* n_id (0, 1 or 2)
*
* flipped
*
* /|\
* / | \
* / | \
* / | \
* t_id | \ t_id (0, 1 or 2)
* \ t | n /
* \ | /
* \ | /
* \ | /
* \|/
* n_id (0, 1 or 2)
*
* @param t target triangle
* @param n neighboring triangle
* @param t_vid reference index of the target triangle
* @param n_vid reference index of the neighboring triangle
*/
void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id)
{
t->vert[(t_id+1)%3] = n->vert[n_id]; // flip t
circumcircle(t->vert[0], t->vert[1], t->vert[2], t->cx, t->cy, t->cr); // update circumcircle
n->vert[(n_id+2)%3] = t->vert[(t_id+2)%3]; // flip n
circumcircle(n->vert[0], n->vert[1], n->vert[2], n->cx, n->cy, n->cr); // update circumcircle
// set side neighbors
t->neigh[t_id] = n->neigh[(n_id+2)%3];
n->neigh[(n_id+1)%3] = t->neigh[(t_id+1)%3];
// set opposite neighbors
t->neigh[(t_id+1)%3] = n;
n->neigh[(n_id+2)%3] = t;
// set oppsite neighbors
if (t->neigh[t_id] != nullptr)
{
for (int i = 0; i < 3; i++)
{
if (t->neigh[t_id]->neigh[i] == n)
{
t->neigh[t_id]->neigh[i] = t;
break;
}
}
}
if (n->neigh[(n_id+1)%3] != nullptr)
{
for (int i = 0; i < 3; i++)
{
if (n->neigh[(n_id+1)%3]->neigh[i] == t)
{
n->neigh[(n_id+1)%3]->neigh[i] = n;
break;
}
}
}
// move hosted DEM points
dem_point *tmp_dem;
std::vector<dem_point*>::iterator d_iter;
for (d_iter = t->hosted_dem.begin(); d_iter != t->hosted_dem.end(); )
{
tmp_dem = *d_iter;
if (n->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = n;
n->hosted_dem.push_back(tmp_dem);
d_iter = t->hosted_dem.erase(d_iter);
}
else d_iter++;
}
for (d_iter = n->hosted_dem.begin(); d_iter != n->hosted_dem.end(); )
{
tmp_dem = *d_iter;
if (t->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = t;
t->hosted_dem.push_back(tmp_dem);
d_iter = n->hosted_dem.erase(d_iter);
}
else d_iter++;
}
// update errors for hosted DEM data
for (int i = 0; i < n->hosted_dem.size(); i++)
{
tmp_dem = n->hosted_dem[i];
tmp_dem->err = fabs(n->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
}
for (int i = 0; i < t->hosted_dem.size(); i++)
{
tmp_dem = t->hosted_dem[i];
tmp_dem->err = fabs(t->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
}
// Sort maximal errors for triangles t and n
std::sort(t->hosted_dem.begin(), t->hosted_dem.end(), compare_dem_point);
std::sort(n->hosted_dem.begin(), n->hosted_dem.end(), compare_dem_point);
return;
}
/**
* @brief Make sure that the input triangle meets the empty circumcircle condition
*
* @param t Input triangle
*/
void make_delaunay(triangle *t)
{
double dist;
vertex2dc *n_vert;
triangle *n_tri;
dem_point *tmp_dem;
for (int n = 0; n < 3; n++)
{
if (t->neigh[n] != nullptr) // must has neighbor on this side
{
n_tri = t->neigh[n];
for (int v = 0; v < 3; v++)
{
n_vert = n_tri->vert[v];
if (n_vert != t->vert[n] && n_vert != t->vert[(n+1)%3]) // find the opposite vertex
{
dist = (t->cx - n_vert->x) * (t->cx - n_vert->x) +
(t->cy - n_vert->y) * (t->cy - n_vert->y);
if ((dist - t->cr) < -1.0*ZERO) // A very restrict condition. The testing point must be really inside the circumcircle
{
flip_neighboring_triangles(t, n_tri, n, v);
// Make sure the triangles also meet the empty circumcircle condition after flipping
make_delaunay(t);
make_delaunay(n_tri);
return; // Neighborhood changed. The current loop is not valid any more.
}
break; // no need to search more
}
}
}
}
return;
}
triangle *split_triangle(vertex2dc *v, triangle *t, triangle *new_t[4])
{
vertex2dc *tmp_vert;
triangle *tmp_tri;
new_t[0] = new_t[1] = new_t[2] = new_t[3] = nullptr;
// Check for collinear
for (int i = 0; i < 3; i++)
{
if (is_collinear(t->vert[i], t->vert[(i+1)%3], v)) // the new vertex is on edge
{
if (t->neigh[i] == nullptr) // no neighboring triangle. create two new triangles
{
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
new_t[0]->set_neighbor(nullptr, new_t[1], t->neigh[(i+2)%3]);
new_t[1]->set_neighbor(new_t[0], nullptr, t->neigh[(i+1)%3]);
for (int n = 0; n < 2; n++)
{
if (new_t[n]->neigh[2] != nullptr)
{
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
{
if (new_t[n]->neigh[2]->neigh[k] == t)
{
new_t[n]->neigh[2]->neigh[k] = new_t[n];
break;
}
}
}
}
return nullptr;
}
// has a neighboring triangle. create four new triangles
for (int k = 0; k < 3; k++)
{
tmp_vert = t->neigh[i]->vert[k];
if (tmp_vert != t->vert[i] && tmp_vert != t->vert[(i+1)%3])
{
tmp_tri = new triangle(t->vert[i], v, t->vert[(i+2)%3]); new_t[0] = tmp_tri;
tmp_tri = new triangle(t->vert[(i+2)%3], v, t->vert[(i+1)%3]); new_t[1] = tmp_tri;
tmp_tri = new triangle(tmp_vert, v, t->vert[i]); new_t[2] = tmp_tri;
tmp_tri = new triangle(t->vert[(i+1)%3], v, tmp_vert); new_t[3] = tmp_tri;
new_t[0]->set_neighbor(new_t[2], new_t[1], t->neigh[(i+2)%3]);
new_t[1]->set_neighbor(new_t[0], new_t[3], t->neigh[(i+1)%3]);
new_t[2]->set_neighbor(new_t[3], new_t[0], t->neigh[i]->neigh[(k+2)%3]);
new_t[3]->set_neighbor(new_t[1], new_t[2], t->neigh[i]->neigh[k]);
for (int n = 0; n < 2; n++)
{
if (new_t[n]->neigh[2] != nullptr)
{
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
{
if (new_t[n]->neigh[2]->neigh[k] == t)
{
new_t[n]->neigh[2]->neigh[k] = new_t[n];
break;
}
}
}
}
for (int n = 2; n < 4; n++)
{
if (new_t[n]->neigh[2] != nullptr)
{
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
{
if (new_t[n]->neigh[2]->neigh[k] == t->neigh[i])
{
new_t[n]->neigh[2]->neigh[k] = new_t[n];
break;
}
}
}
}
break;
}
}
return t->neigh[i]; // Return the neighboring tiangle to be deleted
}
}
// The new vertex is inside the triangle. create three new triangles
for (int n = 0; n < 3; ++n)
{
tmp_tri = new triangle(t->vert[n], t->vert[(n+1)%3], v);
new_t[n] = tmp_tri;
}
// sort neighbors for new triangles
for (int n = 0; n < 3; ++n)
{
if (t->neigh[n] == nullptr)
{
new_t[n]->set_neighbor(nullptr, new_t[(n+1)%3], new_t[(n+2)%3]);
}
else
{
new_t[n]->set_neighbor(t->neigh[n], new_t[(n+1)%3], new_t[(n+2)%3]);
for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle
{
if (t->neigh[n]->neigh[k] == t)
{
t->neigh[n]->neigh[k] = new_t[n];
break;
}
}
}
}
return nullptr;
}
// 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<double> &dem, double xmin, double xmax, double ymin, double ymax, double dx,
double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris, double maxi_err, std::vector<double> *err_records)
{
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_point*> dem_tri;
std::vector<dem_point*>::iterator d_iter;
vertex2dc *tmp_vert = nullptr;
tmp_vert = new vertex2dc(xmin, ymin, dem[0], out_verts.size()); // lower left corner
out_verts.push_back(tmp_vert);
tmp_vert = new vertex2dc(xmax, ymin, dem[xnum-1], out_verts.size()); // lower right corner
out_verts.push_back(tmp_vert);
tmp_vert = new vertex2dc(xmax, ymax, dem[xnum*ynum-1], out_verts.size()); // upper right corner
out_verts.push_back(tmp_vert);
tmp_vert = new vertex2dc(xmin, ymax, dem[xnum*(ynum-1)], out_verts.size()); // upper left corner
out_verts.push_back(tmp_vert);
triangle *old_tri = nullptr, *tmp_tri = nullptr;
triangle *cnst_tri[4];
std::vector<triangle*>::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); tmp_tri = nullptr;
}
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); tmp_tri = nullptr;
}
if (out_tris.size() != 2) return;
out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);
// Find host triangle for all DEM locations
int tmp_id;
for (int i = 0; i < ynum; ++i)
{
for (int j = 0; j < xnum; ++j)
{
tmp_id = j + i*xnum;
if (tmp_id != 0 && tmp_id != (xnum-1) && tmp_id != (xnum*ynum-1) && tmp_id != (xnum*(ynum-1))) // the four corners are already used
{
tmp_dem = new dem_point(xmin + dx*j, ymin + dy*i, dem[j + i*xnum]);
for (int t = 0; t < out_tris.size(); ++t)
{
if (out_tris[t]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = out_tris[t];
tmp_dem->err = fabs(out_tris[t]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
out_tris[t]->hosted_dem.push_back(tmp_dem);
break; // already found, no need to search more
}
}
}
}
}
// Sort hosted_dem in the desceding order with respect to the error. Add maximal zeros to dem_tri
for (int t = 0; t < out_tris.size(); ++t)
{
std::sort(out_tris[t]->hosted_dem.begin(), out_tris[t]->hosted_dem.end(), compare_dem_point);
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
}
// Sort dem_tri
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
while (dem_tri[0]->err >= maxi_err) // quit til the threshold is meet
{
if (err_records != nullptr)
{
err_records->push_back(dem_tri[0]->err);
}
// find the triangle that includes dem_tri[0] and remove it from out_tris
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
old_tri = *t_iter;
if (old_tri == dem_tri[0]->host)
{
t_iter = out_tris.erase(t_iter);
break;
}
else t_iter++;
}
// remove dem_tri[0] from its host triangle's hosted DEM list
for (d_iter = old_tri->hosted_dem.begin(); d_iter != old_tri->hosted_dem.end(); )
{
if (dem_tri[0] == *d_iter)
{
d_iter = old_tri->hosted_dem.erase(d_iter);
break;
}
else d_iter++;
}
// create a new vertex
tmp_vert = new vertex2dc(dem_tri[0]->x, dem_tri[0]->y, dem_tri[0]->elev, out_verts.size());
out_verts.push_back(tmp_vert);
// Delete dem_tri[0]
tmp_dem = dem_tri[0]; delete tmp_dem;
// build new triangles
tmp_tri = split_triangle(tmp_vert, old_tri, cnst_tri);
for (int n = 0; n < 4; ++n)
{
if (cnst_tri[n] != nullptr)
{
out_tris.push_back(cnst_tri[n]);
}
}
if (tmp_tri != nullptr)
{
for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
{
if (tmp_tri == *t_iter)
{
t_iter = out_tris.erase(t_iter);
break;
}
else t_iter++;
}
// build hosted dem for the new triangles
for (int d = 0; d < old_tri->hosted_dem.size(); d++)
{
tmp_dem = old_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem);
break;
}
}
}
for (int d = 0; d < tmp_tri->hosted_dem.size(); d++)
{
tmp_dem = tmp_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem);
break;
}
}
}
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr)
{
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
}
}
// delete the old triangle
old_tri->hosted_dem.clear();
delete old_tri; old_tri = nullptr;
tmp_tri->hosted_dem.clear();
delete tmp_tri; tmp_tri = nullptr;
}
else
{
// build hosted dem for the new triangles
for (int d = 0; d < old_tri->hosted_dem.size(); d++)
{
tmp_dem = old_tri->hosted_dem[d];
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr && cnst_tri[n]->bound_location(tmp_dem->x, tmp_dem->y))
{
tmp_dem->host = cnst_tri[n];
tmp_dem->err = fabs(cnst_tri[n]->interpolate(tmp_dem->x, tmp_dem->y) - tmp_dem->elev);
cnst_tri[n]->hosted_dem.push_back(tmp_dem);
break;
}
}
}
for (int n = 0; n < 4; n++)
{
if (cnst_tri[n] != nullptr)
{
std::sort(cnst_tri[n]->hosted_dem.begin(), cnst_tri[n]->hosted_dem.end(), compare_dem_point);
}
}
// delete the old triangle
old_tri->hosted_dem.clear();
delete old_tri; old_tri = nullptr;
}
// Make sure cnst_tri meet the empty circumcircle condition
for (int n = 0; n < 4; ++n)
{
if (cnst_tri[n] != nullptr)
{
make_delaunay(cnst_tri[n]);
}
}
// get maximal errors from out_tris and sort dem_tri
dem_tri.clear(); dem_tri.reserve(out_tris.size());
for (int t = 0; t < out_tris.size(); t++)
{
if (!out_tris[t]->hosted_dem.empty())
{
dem_tri.push_back(out_tris[t]->hosted_dem[0]);
}
}
std::sort(dem_tri.begin(), dem_tri.end(), compare_dem_point);
}
if (err_records != nullptr)
{
err_records->push_back(dem_tri[0]->err);
}
// assign triangles index
for (int i = 0; i < out_tris.size(); i++)
{
out_tris[i]->id = i;
// destroy remaining DEM data
for (int d = 0; d < out_tris[i]->hosted_dem.size(); d++)
{
tmp_dem = out_tris[i]->hosted_dem[d];
delete tmp_dem; tmp_dem = nullptr;
}
}
return;
}

223
src/lib/tin.h Normal file
View File

@@ -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<dem_point*> 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<double> &dem, double xmin, double xmax, double ymin, double ymax,
double dx, double dy, std::vector<vertex2dc*> &out_verts, std::vector<triangle*> &out_tris,
double maxi_err = 1e-0, std::vector<double> *err_records = nullptr);
#endif // _TIN_DELAUNAY_H