439 lines
9.0 KiB
C++
439 lines
9.0 KiB
C++
#include "dl_math.h"
|
||
|
||
namespace dl::Math
|
||
{
|
||
|
||
std::vector<dl::PointU> GetCircle(const Size& size, const PointU& xy, unsigned rad)
|
||
{
|
||
if (rad == 0)
|
||
return { xy };
|
||
|
||
std::vector<PointU> ret;
|
||
|
||
PointU r = { rad, rad };
|
||
PointU lt = Math::SubU(xy, r);
|
||
PointU rb = Each::Limit(xy + r, {}, size);
|
||
|
||
for (unsigned j = lt[1]; j <= rb[1]; ++j)
|
||
{
|
||
for (unsigned i = lt[0]; i <= rb[0]; ++i)
|
||
{
|
||
if (Math::DistancePow2U({i, j}, xy) <= rad * rad)
|
||
ret.push_back(PointU{ i, j });
|
||
}
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::PointU> RingX::GetCircle(const Size& size, const PointU& xy, unsigned rad)
|
||
{
|
||
if (rad == 0)
|
||
return { xy };
|
||
|
||
std::vector<PointU> ret;
|
||
ret.reserve(GetCircleNumMax(rad));
|
||
int r = (int)rad;
|
||
|
||
Point xy_i = toPoint(xy);
|
||
for (int y = -r; y <= r; ++y)
|
||
{
|
||
// 计算该行最大x偏移
|
||
int maxX = (int)std::sqrt(r * r - y * y);
|
||
for (int x = -maxX; x <= maxX; ++x)
|
||
{
|
||
PointU p;
|
||
if (Clamp(size, xy_i + Point{ x, y }, p))
|
||
ret.push_back(p);
|
||
}
|
||
}
|
||
assert(ret.size() <= GetCircleNumMax(rad));
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::PointU> RingX::Diffusion(const Size& size, const PointU& xy, unsigned count, FnDiffusionCheck func_valid, FnDiffusionCb func_cb)
|
||
{
|
||
MapLimit map;
|
||
map.SetSize(size);
|
||
|
||
// 返回值
|
||
std::vector<PointU> ret;
|
||
// 是否已经遍历过
|
||
std::unordered_map<unsigned, bool> map_index;
|
||
|
||
// 当前需要扩散的格子
|
||
std::list<PointU> list_open;
|
||
// 加入首节点
|
||
unsigned index = xy.x + xy.y * size.w;
|
||
map_index[index] = true;
|
||
list_open.push_back(xy);
|
||
|
||
for (unsigned i = 0; i < count; ++i)
|
||
{
|
||
std::list<PointU> list_temp; // 此轮遍历到的格子
|
||
for (auto& iter : list_open)
|
||
{
|
||
if (func_valid(iter))
|
||
{
|
||
// 添加到返回列表
|
||
ret.push_back(iter);
|
||
|
||
PointU xy_near;
|
||
xy_near = iter;
|
||
if (map.Left(xy_near))
|
||
{
|
||
unsigned index = xy_near.x + xy_near.y * size.w;
|
||
if (!map_index.contains(index))
|
||
{
|
||
map_index[index] = true;
|
||
list_temp.push_back(xy_near);
|
||
}
|
||
}
|
||
|
||
xy_near = iter;
|
||
if (map.Right(xy_near))
|
||
{
|
||
unsigned index = xy_near.x + xy_near.y * size.w;
|
||
if (!map_index.contains(index))
|
||
{
|
||
map_index[index] = true;
|
||
list_temp.push_back(xy_near);
|
||
}
|
||
}
|
||
|
||
xy_near = iter;
|
||
if (map.Bottom(xy_near))
|
||
{
|
||
unsigned index = xy_near.x + xy_near.y * size.w;
|
||
if (!map_index.contains(index))
|
||
{
|
||
map_index[index] = true;
|
||
list_temp.push_back(xy_near);
|
||
}
|
||
}
|
||
|
||
xy_near = iter;
|
||
if (map.Top(xy_near))
|
||
{
|
||
unsigned index = xy_near.x + xy_near.y * size.w;
|
||
if (!map_index.contains(index))
|
||
{
|
||
map_index[index] = true;
|
||
list_temp.push_back(xy_near);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
func_cb(i, list_temp);
|
||
list_open = list_temp;
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::Point> GetCircle(const Point& xy, unsigned rad)
|
||
{
|
||
if (rad == 0)
|
||
return { xy };
|
||
|
||
std::vector<Point> ret;
|
||
|
||
int r = rad;
|
||
for (int j = -r; j <= r; ++j)
|
||
{
|
||
for (int i = -r; i <= r; ++i)
|
||
{
|
||
if (i * i + j * j <= r * r)
|
||
ret.push_back(Point{ i, j } + xy);
|
||
}
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
dl::Position4 Slerp(const Position4& q0, const Position4& q1, const Float& t)
|
||
{
|
||
Float magnitude = sqrt(DistancePow2(q0) * DistancePow2(q1));
|
||
assert(magnitude > 0);
|
||
Float product = DotProduct(q0, q1) / magnitude;
|
||
Float product_abs = std::abs(product);
|
||
if (product_abs < 1_f - std::numeric_limits<Float>::epsilon())
|
||
{
|
||
Float theta = acos(product_abs);
|
||
Float d = sin(theta);
|
||
assert(d > 0);
|
||
|
||
Float sign = product < 0 ? -1_f : 1_f;
|
||
Float s0 = sin((1_f - t) * theta) / d;
|
||
Float s1 = sin(sign * t * theta) / d;
|
||
|
||
return q0 * s0 + q1 * s1;
|
||
}
|
||
else
|
||
return q0;
|
||
}
|
||
|
||
std::vector<dl::Point> DrawLineLerp(const Point& p0, const Point& p1)
|
||
{
|
||
std::vector<Point> ret;
|
||
int N = DistanceDiagonal(p0, p1);
|
||
for (int step = 0; step <= N; ++step)
|
||
{
|
||
double t = N ? ((double)step / N) : 0.0;
|
||
ret.push_back(Each::Lerp(p0, p1, t, std::round));
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::Point> DrawLineBresenham(const Point& p0, const Point& p1)
|
||
{
|
||
std::vector<Point> ret;
|
||
|
||
const int& x0 = p0[0];
|
||
const int& y0 = p0[1];
|
||
const int& x1 = p1[0];
|
||
const int& y1 = p1[1];
|
||
|
||
int dx = x1 - x0;
|
||
int dy = y1 - y0;
|
||
int ux = dx > 0 ? 1 : -1;
|
||
int uy = dy > 0 ? 1 : -1;
|
||
dx = abs(dx);
|
||
dy = abs(dy);
|
||
int dx_2 = dx * 2;
|
||
int dy_2 = dy * 2;
|
||
if (dx > dy)
|
||
{
|
||
int e = -dx;
|
||
int x = x0;
|
||
int y = y0;
|
||
|
||
for (x = x0; x != x1; x += ux)
|
||
{
|
||
ret.push_back({ x, y });
|
||
e = e + dy_2;
|
||
if (e > 0)
|
||
{
|
||
y += uy;
|
||
e = e - dx_2;
|
||
}
|
||
}
|
||
ret.push_back({ x, y });
|
||
}
|
||
else
|
||
{
|
||
int e = -dy;
|
||
int x = x0;
|
||
int y = y0;
|
||
for (y = y0; y != y1; y += uy)
|
||
{
|
||
ret.push_back({ x, y });
|
||
e = e + dx_2;
|
||
if (e > 0)
|
||
{
|
||
x += ux;
|
||
e = e - dy_2;
|
||
}
|
||
}
|
||
ret.push_back({ x, y });
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::Point> DrawLineOrth(const Point& p0, const Point& p1)
|
||
{
|
||
std::vector<Point> ret;
|
||
|
||
int dx = p1[0] - p0[0];
|
||
int dy = p1[1] - p0[1];
|
||
int nx = abs(dx);
|
||
int ny = abs(dy);
|
||
int sign_x = dx > 0 ? 1 : -1;
|
||
int sign_y = dy > 0 ? 1 : -1;
|
||
|
||
Point p = p0;
|
||
ret.push_back(p);
|
||
int ix = 0;
|
||
int iy = 0;
|
||
for (; ix < nx || iy < ny;)
|
||
{
|
||
int decision = (1 + 2 * ix) * ny - (1 + 2 * iy) * nx;
|
||
if (decision < 0)
|
||
{//水平
|
||
p[0] += sign_x;
|
||
ix++;
|
||
}
|
||
else
|
||
{//垂直
|
||
p[1] += sign_y;
|
||
iy++;
|
||
}
|
||
ret.push_back(p);
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
std::vector<dl::Point> DrawLineCover(const Point& p0, const Point& p1)
|
||
{
|
||
std::vector<Point> ret;
|
||
|
||
int dx = p1[0] - p0[0];
|
||
int dy = p1[1] - p0[1];
|
||
int nx = abs(dx);
|
||
int ny = abs(dy);
|
||
int sign_x = dx > 0 ? 1 : -1;
|
||
int sign_y = dy > 0 ? 1 : -1;
|
||
|
||
Point p = p0;
|
||
ret.push_back(p);
|
||
int ix = 0;
|
||
int iy = 0;
|
||
for (; ix < nx || iy < ny;)
|
||
{
|
||
int decision = (1 + 2 * ix) * ny - (1 + 2 * iy) * nx;
|
||
if (decision == 0)
|
||
{//对角线
|
||
p[0] += sign_x;
|
||
p[1] += sign_y;
|
||
ix++;
|
||
iy++;
|
||
}
|
||
else if (decision < 0)
|
||
{//水平
|
||
p[0] += sign_x;
|
||
ix++;
|
||
}
|
||
else
|
||
{//垂直
|
||
p[1] += sign_y;
|
||
iy++;
|
||
}
|
||
ret.push_back(p);
|
||
}
|
||
return ret;
|
||
}
|
||
|
||
void MergePath(std::list<Point>& list_pos)
|
||
{
|
||
size_t size = list_pos.size();
|
||
if (list_pos.size() < 3)
|
||
return;//少于3个节点不用处理
|
||
//初始向量
|
||
auto iter_0 = list_pos.begin();//前一个节点
|
||
auto iter_1 = iter_0;//后一个节点
|
||
++iter_1;
|
||
Point v = *iter_1 - *iter_0;//向量指向前一个节点
|
||
++iter_0;
|
||
++iter_1;
|
||
|
||
//直到后一个到end
|
||
for (; iter_1 != list_pos.end();)
|
||
{
|
||
Point v1 = *iter_1 - *iter_0;
|
||
if (v1 == v)
|
||
{//方向相同，删除前一个节点
|
||
list_pos.erase(iter_0);
|
||
iter_0 = iter_1++;//iter_0失效、指向iter_1，iter_1指向下一位置
|
||
}
|
||
else
|
||
{
|
||
v = v1;
|
||
++iter_0;
|
||
++iter_1;
|
||
}
|
||
}
|
||
}
|
||
|
||
double GetSolarDeclination(int day)
|
||
{
|
||
// 近似公式：赤纬 = 23.44° × sin(360/365 × (284 + N))
|
||
// 这是地理学中常用的经验公式
|
||
double angle = (360.0 / 365.0) * (284 + day);
|
||
double declination = 23.44 * std::sin(DegreesToRadians(angle));
|
||
return declination;
|
||
}
|
||
|
||
double GetSolarElevation(double latitude, double declination)
|
||
{
|
||
// H = 90° - |纬度 - 赤纬|
|
||
double diff = fabs(latitude - declination);
|
||
double elevation = 90.0 - diff;
|
||
return elevation;
|
||
}
|
||
|
||
int GetSeason(int day, double latitude)
|
||
{
|
||
bool is_northern = (latitude >= 0);
|
||
|
||
// 定义节气的大致日期（北半球）
|
||
// 春分：3月20日（第79天）
|
||
// 夏至：6月21日（第172天）
|
||
// 秋分：9月23日（第266天）
|
||
// 冬至：12月22日（第356天）
|
||
|
||
int spring_equinox = 79;
|
||
int summer_solstice = 172;
|
||
int autumn_equinox = 266;
|
||
int winter_solstice = 356;
|
||
|
||
if (!is_northern)
|
||
{
|
||
// 南半球：季节相反，交换春分/秋分，夏至/冬至
|
||
std::swap(spring_equinox, autumn_equinox);
|
||
std::swap(summer_solstice, winter_solstice);
|
||
}
|
||
|
||
// 判断季节（使用环形日期）
|
||
// 为了处理跨年，将日期映射到0-364的环上
|
||
int day_ring = day % 365;
|
||
|
||
// 简化判断：将四个节点之间的区间分配给四个季节
|
||
// 这里用区间中点划分更准确
|
||
if (day_ring >= spring_equinox && day_ring < summer_solstice)
|
||
{
|
||
return 0; // 春季
|
||
}
|
||
else if (day_ring >= summer_solstice && day_ring < autumn_equinox)
|
||
{
|
||
return 1; // 夏季
|
||
}
|
||
else if (day_ring >= autumn_equinox && day_ring < winter_solstice)
|
||
{
|
||
return 2; // 秋季
|
||
}
|
||
else
|
||
{
|
||
return 3; // 冬季（包含跨年部分）
|
||
}
|
||
}
|
||
|
||
double GetBaseTemperature(double latitude, int day)
|
||
{
|
||
// 基础纬度降温效应（年均温基准）
|
||
double base_temp = 26.0 - (std::fabs(latitude) * 0.6);
|
||
|
||
// 太阳高度角
|
||
double solar_elevation = GetSolarElevation(latitude, GetSolarDeclination(day));
|
||
|
||
// 太阳高度角修正系数
|
||
// 高度角为90°时热力系数为1.0，为0°时系数为0.0
|
||
double heat_factor = std::sin(DegreesToRadians(solar_elevation));
|
||
|
||
// 确保热力系数在合理范围
|
||
heat_factor = (std::max)(0.1, (std::min)(1.0, heat_factor));
|
||
|
||
// 季节修正（太阳高度角已经反映了季节变化，但这里做额外微调）
|
||
double seasonal_offset = 0.0;
|
||
if (GetSeason(day, latitude) == 1)
|
||
{
|
||
seasonal_offset = 2.0; // 夏季额外增温
|
||
}
|
||
else if (GetSeason(day, latitude) == 3)
|
||
{
|
||
seasonal_offset = -2.0; // 冬季额外降温
|
||
}
|
||
|
||
// 4. 组合计算：基准温度 × 热力系数 + 季节偏移
|
||
double estimated_temp = base_temp * heat_factor + seasonal_offset;
|
||
|
||
return estimated_temp;
|
||
}
|
||
|
||
} |