libtin/src/lib/tin.h

/**
 *   ___________   __
 *  /_  __/  _/ | / /
 *   / /  / //  |/ / 
 *  / / _/ // /|  /  
 * /_/ /___/_/ |_/ 
 * 
 * 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