519 lines
13 KiB
C++
519 lines
13 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _GCTL_POINT3S_H
|
|
#define _GCTL_POINT3S_H
|
|
|
|
#include "../core/macro.h"
|
|
#include "../core/exceptions.h"
|
|
#include "point3c.h"
|
|
|
|
#include "iostream"
|
|
#include "string"
|
|
#include "cmath"
|
|
#include "iomanip"
|
|
#include "regex"
|
|
|
|
namespace gctl
|
|
{
|
|
template <typename T> struct point3c;
|
|
template <typename T> struct point3s;
|
|
|
|
typedef point3s<double> point3ds;
|
|
typedef point3s<float> point3fs;
|
|
|
|
#ifndef IO_PSN
|
|
#define IO_PSN
|
|
// static variable for controlling the IO process
|
|
static int io_psn = 6;
|
|
#endif // IO_PSN
|
|
|
|
/**
|
|
* @brief 球坐标系中的一个点（经纬度表示的）
|
|
*/
|
|
template <typename T>
|
|
struct point3s
|
|
{
|
|
T rad, lon, lat;
|
|
// 构造函数与拷贝构造函数
|
|
point3s();
|
|
point3s(T r, T ln, T lt);
|
|
point3s(const point3s<T> &b);
|
|
// 析构函数
|
|
virtual ~point3s(){}
|
|
// 结构体方法
|
|
bool valid() const;
|
|
void set(T r, T ln, T lt);
|
|
void set(const point3s<T> &b);
|
|
double module() const;
|
|
point3c<T> s2c() const; // 转换为球坐标点
|
|
void str(std::string str);
|
|
/**
|
|
* @brief Sets the i/o precision of the type.
|
|
*
|
|
* @param[in] psn The desired precision
|
|
*/
|
|
void set_io_precision(int psn);
|
|
/**
|
|
* @brief 输出位置
|
|
*
|
|
* @param os 输出流
|
|
*/
|
|
void out_loc(std::ostream &os, char deli) const;
|
|
/**
|
|
* @brief 输入位置
|
|
*
|
|
* @param os 输入流
|
|
*/
|
|
void in_loc(std::istream &os);
|
|
/**
|
|
* @brief 返回位置
|
|
*
|
|
* @return 位置
|
|
*/
|
|
point3s<T> get_loc() const;
|
|
/**
|
|
* @brief 返回正常球坐标中的Phi值
|
|
*
|
|
* @return double
|
|
*/
|
|
T get_phi() const;
|
|
/**
|
|
* @brief 返回正常球坐标中的Theta值
|
|
*
|
|
* @return double
|
|
*/
|
|
T get_theta() const;
|
|
};
|
|
|
|
template <typename T>
|
|
gctl::point3s<T>::point3s()
|
|
{
|
|
rad = lon = lat = NAN;
|
|
}
|
|
|
|
template <typename T>
|
|
gctl::point3s<T>::point3s(T r, T ln, T lt)
|
|
{
|
|
set(r, ln, lt);
|
|
}
|
|
|
|
template <typename T>
|
|
gctl::point3s<T>::point3s(const point3s<T> &b)
|
|
{
|
|
set(b);
|
|
}
|
|
|
|
template <typename T>
|
|
bool gctl::point3s<T>::valid() const
|
|
{
|
|
if (std::isnan(rad) || std::isnan(lon) || std::isnan(lat)) return false;
|
|
if (std::isinf(rad) || std::isinf(lon) || std::isinf(lat)) return false;
|
|
if (rad < 0 || lon > 180 || lon < -180 || lat > 90 || lat < -90) return false;
|
|
return true;
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::set(T r, T ln, T lt)
|
|
{
|
|
if (std::isnan(r) || std::isnan(ln) || std::isnan(lt) ||
|
|
std::isinf(r) || std::isinf(ln) || std::isinf(lt))
|
|
{
|
|
throw invalid_argument("Invalid value detected. From point3ds::set(...)");
|
|
}
|
|
|
|
if (r < 0 || ln > 180 || ln < -180 || lt > 90 || lt < -90)
|
|
{
|
|
throw out_of_range("Invalid values. From point3s::set(...)");
|
|
}
|
|
|
|
rad = r; lon = ln; lat = lt;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::set(const point3s<T> &b)
|
|
{
|
|
if (std::isnan(b.rad) || std::isnan(b.lon) || std::isnan(b.lat) ||
|
|
std::isinf(b.rad) || std::isinf(b.lon) || std::isinf(b.lat))
|
|
{
|
|
throw invalid_argument("Invalid value detected. From point3s::set(...)");
|
|
}
|
|
|
|
if (b.rad < 0 || b.lon > 180 || b.lon < -180 || b.lat > 90 || b.lat < -90)
|
|
{
|
|
throw out_of_range("Invalid radius. From point3s::set(...)");
|
|
}
|
|
|
|
rad = b.rad; lon = b.lon; lat = b.lat;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
double gctl::point3s<T>::module() const
|
|
{
|
|
return rad;
|
|
}
|
|
|
|
template <typename T>
|
|
gctl::point3c<T> gctl::point3s<T>::s2c() const
|
|
{
|
|
/*
|
|
point3c<T> outc;
|
|
outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
|
|
outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
|
|
outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
|
|
return outc;
|
|
*/
|
|
|
|
point3c<T> v;
|
|
v.x = rad*cos(GCTL_Pi*lat/180.0)*cos(GCTL_Pi*lon/180.0);
|
|
v.y = rad*cos(GCTL_Pi*lat/180.0)*sin(GCTL_Pi*lon/180.0);
|
|
v.z = rad*sin(GCTL_Pi*lat/180.0);
|
|
return v;
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::str(std::string str)
|
|
{
|
|
std::smatch ret;
|
|
std::regex pattern("\\((-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?)\\)");
|
|
|
|
if (regex_search(str, ret, pattern))
|
|
{
|
|
rad = atof(std::string(ret[1]).c_str());
|
|
lon = atof(std::string(ret[2]).c_str());
|
|
lat = atof(std::string(ret[3]).c_str());
|
|
return;
|
|
}
|
|
|
|
throw runtime_error("Fail to parse the input string: " + str + ". From point3s::str(...)");
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::set_io_precision(int psn)
|
|
{
|
|
if (psn < 0)
|
|
{
|
|
throw invalid_argument("Invalid precision. From point3s::set_io_precision(...)");
|
|
}
|
|
|
|
io_psn = psn;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::out_loc(std::ostream &os, char deli) const
|
|
{
|
|
os << std::setprecision(io_psn) << rad << deli << lon << deli << lat;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void gctl::point3s<T>::in_loc(std::istream &os)
|
|
{
|
|
os >> rad >> lon >> lat;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
gctl::point3s<T> gctl::point3s<T>::get_loc() const
|
|
{
|
|
return point3s<T>(rad, lon, lat);
|
|
}
|
|
|
|
template <typename T>
|
|
T gctl::point3s<T>::get_phi() const
|
|
{
|
|
if (lon < .0) return GCTL_Pi*(2.0 - lon/180.0);
|
|
return GCTL_Pi*lon/180.0;
|
|
}
|
|
|
|
template <typename T>
|
|
T gctl::point3s<T>::get_theta() const
|
|
{
|
|
return GCTL_Pi*(0.5 - lat/180.0);
|
|
}
|
|
|
|
template <typename T>
|
|
bool operator ==(const point3s<T> &a, const point3s<T> &b)
|
|
{
|
|
if( fabs(a.lon - b.lon) < GCTL_ZERO &&
|
|
fabs(a.lat - b.lat) < GCTL_ZERO &&
|
|
fabs(a.rad - b.rad) < GCTL_ZERO)
|
|
{
|
|
return true;
|
|
}
|
|
else return false;
|
|
}
|
|
|
|
template <typename T>
|
|
bool operator !=(const point3s<T> &a, const point3s<T> &b)
|
|
{
|
|
if( fabs(a.lon - b.lon) >= GCTL_ZERO ||
|
|
fabs(a.lat - b.lat) >= GCTL_ZERO ||
|
|
fabs(a.rad - b.rad) >= GCTL_ZERO)
|
|
{
|
|
return true;
|
|
}
|
|
else return false;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator -(const point3s<T> &a, const point3s<T> &b)
|
|
{
|
|
point3s<T> m;
|
|
m.rad = a.rad - b.rad;
|
|
m.lon = a.lon - b.lon;
|
|
m.lat = a.lat - b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator +(const point3s<T> &a, const point3s<T> &b)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = a.rad + b.rad;
|
|
m.lon = a.lon + b.lon;
|
|
m.lat = a.lat + b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(int sign, const point3s<T> &b)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(float sign, const point3s<T> &b)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(double sign, const point3s<T> &b)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(const point3s<T> &b, int sign)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(const point3s<T> &b, float sign)
|
|
{
|
|
gctl::point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
point3s<T> operator *(const point3s<T> &b, double sign)
|
|
{
|
|
point3s<T> m;
|
|
m.rad = sign*b.rad;
|
|
m.lon = sign*b.lon;
|
|
m.lat = sign*b.lat;
|
|
|
|
while (m.lon > 180) m.lon -= 360;
|
|
while (m.lon < -180) m.lon += 360;
|
|
while (m.lat > 90) m.lat = 180 - m.lat;
|
|
while (m.lat < -90) m.lat = -180 - m.lat;
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
std::ostream &operator <<(std::ostream & os, const point3s<T> &a)
|
|
{
|
|
os << std::setprecision(io_psn) << a.rad << " " << a.lon << " " << a.lat;
|
|
return os;
|
|
}
|
|
|
|
template <typename T>
|
|
std::istream &operator >>(std::istream & os, point3s<T> &a)
|
|
{
|
|
os >> a.rad >> a.lon >> a.lat;
|
|
return os;
|
|
}
|
|
|
|
/**
|
|
* @brief 函数判断两个point3ds类型是否等于
|
|
*
|
|
* @param[in] a 三维空间内的一个向量或实点a。
|
|
* @param[in] b 三维空间内的一个向量或实点b。
|
|
* @param[in] cut_off 截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
|
|
* 注意比较的精度为 cut_off 减一位，比如 cut_off 为1e-10则比较的精度为小数点后9位。
|
|
*
|
|
* @return 是否相等。
|
|
*/
|
|
template <typename T>
|
|
bool isequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
|
|
{
|
|
if( fabs(a.rad - b.rad) < cut_off &&
|
|
fabs(a.lon - b.lon) < cut_off &&
|
|
fabs(a.lat - b.lat) < cut_off)
|
|
{
|
|
return true;
|
|
}
|
|
else return false;
|
|
}
|
|
|
|
/**
|
|
* @brief 函数判断两个point3ds类型是否不等于
|
|
*
|
|
* @param[in] a 三维空间内的一个向量或实点a。
|
|
* @param[in] b 三维空间内的一个向量或实点b。
|
|
* @param[in] cut_off 截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
|
|
*
|
|
* @return 是否不相等。
|
|
*/
|
|
template <typename T>
|
|
bool notequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
|
|
{
|
|
if( fabs(a.rad - b.rad) >= cut_off ||
|
|
fabs(a.lon - b.lon) >= cut_off ||
|
|
fabs(a.lat - b.lat) >= cut_off)
|
|
{
|
|
return true;
|
|
}
|
|
else return false;
|
|
}
|
|
|
|
/**
|
|
* @brief 生成直球坐标系的点数组
|
|
*
|
|
* @param out_ps 返回的点数组
|
|
* @param[in] lonmin lon最小值
|
|
* @param[in] lonmax lon最大值
|
|
* @param[in] latmin lat最小值
|
|
* @param[in] latmax lat最大值
|
|
* @param[in] dlon lon间隔
|
|
* @param[in] dlat lat间隔
|
|
* @param[in] rad 半径值
|
|
* @param[in] cor 网格起始位置
|
|
*/
|
|
template <typename T>
|
|
void grid_points_2d(array<point3s<T>> &out_ps, T lonmin, T lonmax, T latmin,
|
|
T latmax, T dlon, T dlat, T rad, matrix_corner_e cor = BtmLeft)
|
|
{
|
|
if (lonmin >= lonmax || latmin >= latmax ||
|
|
lonmin + dlon > lonmax || latmin + dlat > latmax ||
|
|
dlon <= 0 || dlat <= 0)
|
|
{
|
|
throw invalid_argument("[gctl::grid_points_2d] Invalid parameters.");
|
|
}
|
|
|
|
int lonnum = round((lonmax-lonmin)/dlon) + 1;
|
|
int latnum = round((latmax-latmin)/dlat) + 1;
|
|
|
|
out_ps.resize(lonnum*latnum);
|
|
if (cor == BtmLeft)
|
|
{
|
|
for (int j = 0; j < latnum; j++)
|
|
{
|
|
for (int i = 0; i < lonnum; i++)
|
|
{
|
|
out_ps.at(j*lonnum+i).lon = lonmin + dlon*i;
|
|
out_ps.at(j*lonnum+i).lat = latmin + dlat*j;
|
|
out_ps.at(j*lonnum+i).rad = rad;
|
|
}
|
|
}
|
|
}
|
|
else // cor == TopLeft
|
|
{
|
|
for (int j = 0; j < latnum; j++)
|
|
{
|
|
for (int i = 0; i < lonnum; i++)
|
|
{
|
|
out_ps.at(j*lonnum+i).lon = lonmin + dlon*i;
|
|
out_ps.at(j*lonnum+i).lat = latmax - dlat*j;
|
|
out_ps.at(j*lonnum+i).rad = rad;
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
}
|
|
|
|
#endif // _GCTL_POINT3S_H
|