📄 tcston2.cpp
字号:
// TcsTon2.cpp
// 邓雪清 2008/11/17
#include <math.h>
#include <assert.h>
#include "TcsTon2.h"
///////////////////////////////////////////////////////////////////
// 公共函数
// 散列函数
TCS_S32 SIDE_EQLH(void *s1, void *s2)
{
TIN_S *m1 = (TIN_S *)s1;
TIN_S *m2 = (TIN_S *)s2;
return (m1->sid[0] == m2->sid[0] && m1->eid[0] == m2->eid[0]);
}
TCS_U32 SIDE_HASH(void *s1)
{
TIN_S *m1 = (TIN_S *)s1;
return (EQHASH(m1->sid[0]) ^ EQHASH(m1->eid[0]));
}
// 计算坐标边界
inline TCS_V00 ComBnd2(TIN_G *G, TCS_F64 *P, TCS_S32 N)
{
G->bnd[0] = G->bnd[2] = P[0]; P++; // x
G->bnd[1] = G->bnd[3] = P[0]; P++; // y
for (TCS_S32 i = 1; i < N; i++)
{
if (G->bnd[0] > P[0]) G->bnd[0] = P[0];
else if (G->bnd[2] < P[0]) G->bnd[2] = P[0];
P++;
if (G->bnd[1] > P[0]) G->bnd[1] = P[0];
else if (G->bnd[3] < P[0]) G->bnd[3] = P[0];
P++;
}
}
// 计算网格矩阵
inline TCS_V00 ComSan2(TIN_G *G, TCS_S32 N, TCS_S32 av)
{
TCS_S32 *num = G->num;
TCS_F64 dup = G->dup;
TCS_F64 *bnd = G->bnd;
TCS_F64 *spn = G->spn;
TCS_F64 sum = (bnd[2] - bnd[0]) * (bnd[3] - bnd[1]);
if (sum == 0.0) // 至少存在1维共线
{ num[0] = 0; return; }
sum = sum / N * av;
spn[0] = spn[1] = sqrt(sum);
if (dup >= spn[0]) // 便于处理重合点
spn[0] = spn[1] = 2 * dup;
num[0] = TCS_S32((bnd[2] - bnd[0]) / spn[0]) + 1;
num[1] = TCS_S32((bnd[3] - bnd[1]) / spn[1]) + 1;
bnd[2] = bnd[0] + num[0] * spn[0];
bnd[3] = bnd[1] + num[1] * spn[1];
}
// 计算网格单元
inline TCS_S32 ComIdx2(TIN_G *G, TCS_F64 *p)
{
TCS_S32 c = TCS_S32((p[0] - G->bnd[0]) / G->spn[0]);
TCS_S32 r = TCS_S32((p[1] - G->bnd[1]) / G->spn[1]);
return (r * G->num[0] + c);
}
inline TCS_V00 ComIdx2(TIN_G *G, TCS_F64 *p, TCS_S32 *o)
{
o[0] = TCS_S32((p[0] - G->bnd[0]) / G->spn[0]);
o[1] = TCS_S32((p[1] - G->bnd[1]) / G->spn[1]);
}
// 单元增加点号
inline TCS_V00 AddPid(TIN_C *C, TCS_S32 pid)
{
if (C->sum < C->num)
{
C->pid[C->sum] = pid;
C->sum++; C->ocn++;
}
else
{
TCS_S32 num = C->sum + 4;
TCS_S32 *buf = new TCS_S32[num];
::memcpy(buf, C->pid, sizeof(TCS_S32) * C->sum);
buf[C->sum] = pid; C->sum++; C->ocn++;
C->num = num; delete []C->pid; C->pid = buf;
}
}
// 单元移去点号
// 基于位置
inline TCS_V00 RmxPid(TIN_C *C, TCS_S32 pos)
{
C->ocn--; C->sum--;
C->pid[pos] = C->pid[C->ocn];
}
// 基于点号
inline TCS_V00 RmvPid(TIN_C *C, TCS_S32 pid)
{
TCS_S32 *P = C->pid, ocn = C->ocn;
for (TCS_S32 n = 0; n < ocn; n++)
{
if (P[n] == pid)
{
ocn--; P[n] = P[ocn]; P[ocn] = pid;
C->ocn = ocn; return;
}
}
}
// 建立网格索引
inline TCS_V00 CnsIdx2(TIN_G *G, TCS_F64 *P, TCS_S32 N, TCS_S32 av)
{
G->idx = 0;
ComSan2(G, N, av);
TCS_S32 n = G->num[0];
if (n == 0) return;
n *= G->num[1];
G->idx = new TIN_C[n];
if (!G->idx) return;
::memset(G->idx, 0, sizeof(TIN_C) * n);
for (TCS_S32 k, i = 0; i < N; i++)
{
k = ComIdx2(G, P); P += 2;
AddPid(G->idx + k, i);
}
}
// 计算外接圆心
inline TCS_V00 ComCen2(TCS_F64 *p1, TCS_F64 *p2, TCS_F64 *p3, TCS_F64 *ce)
{
// 11(+|-), 12(*), 2(/)
TCS_F64 dx21 = p2[0] - p1[0];
TCS_F64 dy21 = p2[1] - p1[1];
TCS_F64 xy21 = dx21 * dx21 + dy21 * dy21;
TCS_F64 dx31 = p3[0] - p1[0];
TCS_F64 dy31 = p3[1] - p1[1];
TCS_F64 xy31 = dx31 * dx31 + dy31 * dy31;
TCS_F64 top1 = dy21 * xy31 - dy31 * xy21;
TCS_F64 top2 = - dx21 * xy31 + dx31 * xy21;
TCS_F64 deno = dy21 * dx31 - dy31 * dx21;
ce[0] = p1[0] + 0.5 * top1 / deno;
ce[1] = p1[1] + 0.5 * top2 / deno;
}
// 判断是否共线
// 要求点互不重合
inline TCS_S32 IsLine2(TCS_F64 *p1, TCS_F64 *p2, TCS_F64 *p3)
{
// 快速精确处理p2 == p3问题
if (p2 == p3) // 如果没有该步, 由于计算误差可能不能精确处理p2 == p3问题
return 1;
// 利用三角形面积是否等于0判断
// (p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) + p3[0] * (p1[1] - p2[1]) : 3(*), 5(+|-)
// (p1[0] - p2[0]) * (p1[1] - p3[1]) - (p1[0] - p3[0]) * (p1[1] - p2[1]) : 2(*), 5(-)
// 下面可以精确处理p1 == p3问题
if ((p1[0] - p2[0]) * (p1[1] - p3[1]) - (p1[0] - p3[0]) * (p1[1] - p2[1]) == 0.0)
return 1;
else
return 0;
}
// 计算同侧符号
// 要求点互不重合
inline TCS_S32 IsSide2(TCS_F64 *p1, TCS_F64 *p2, TCS_F64 *p3)
{
// 快速精确处理p2 == p3问题
if (p2 == p3) // 如果没有该步, 由于计算误差可能不能精确处理p2 == p3问题
return 0;
// 根据三点共线条件设计
// 三点共线条件: (p1[0] - p2[0]) * (p1[1] - p3[1]) - (p1[0] - p3[0]) * (p1[1] - p2[1]) = 0
TCS_F64 side = (p1[0] - p2[0]) * (p1[1] - p3[1]) - (p1[0] - p3[0]) * (p1[1] - p2[1]);
if (side > 0)
return 1;
else if (side < 0)
return -1;
return 0;
}
// 欧拉距离平方
inline TCS_F64 EuLen2D(TCS_F64 *p1, TCS_F64 *p2)
{
TCS_F64 dx = p1[0] - p2[0];
TCS_F64 dy = p1[1] - p2[1];
return (dx * dx + dy * dy);
}
// 共圆点弧位排序
inline TCS_V00 ArcPos2(CTcsBL *L, TCS_F64 *P, TCS_F64 *ce, TCS_F64 d2)
{
d2 = sqrt(d2);
TCS_F64 dx, dy, *pt;
TCS_S32 n = L->GetSize();
TCS_F64 *s = new TCS_F64[n];
// 计算与圆心坐标差
TCS_S32 m, k, i = 0;
TIN_X *X = (TIN_X*)L->GetHead();
while (X)
{
pt = P + 2 * X->pid;
dx = pt[0] - ce[0]; dy = pt[1] - ce[1];
// 第一象限
if (dx >= 0 && dy >= 0) s[i] = dx / d2;
// 第二象限
else if (dx >= 0 && dy < 0) s[i] = 1 - dy / d2;
// 第三象限
else if (dx < 0 && dy <= 0) s[i] = 2 - dx / d2;
// 第四象限
else s[i] = 3 - dy / d2;
i++; X = (TIN_X*)L->GetNext();
}
// 排序
m = n - 1;
TCS_V00 *I, *K;
for (i = 0; i < m; i++)
{
for (k = i + 1; k < n; k++)
{
if (s[i] > s[k])
{
// 交换
dx = s[i]; s[i] = s[k]; s[k] = dx;
I = L->GetAt(i); K = L->GetAt(k);
L->SetAt(k, I); L->SetAt(i, K);
}
}
}
delete []s;
}
// 赋值边
inline TCS_V00 IniSide(TIN_S *S, TCS_S32 sip, TCS_S32 eip, TCS_S32 sic, TCS_S32 eic, TCS_S32 pid)
{
if (sip > eip)
{
S->sid[0] = sip; S->eid[0] = eip; S->sid[1] = sic; S->eid[1] = eic;
}
else
{
S->sid[0] = eip; S->eid[0] = sip; S->sid[1] = eic; S->eid[1] = sic;
}
S->pid[0] = pid;
}
// 插入边
inline TCS_S32 InsSide(CTcsHL *ASL, CTcsHL *BSL, TIN_G *G, TIN_S *S)
{
TIN_S *T = (TIN_S*)(ASL->DetData(S));
if (T) // 已存在
{
assert(T->sid[0] == S->sid[0] && T->eid[0] == S->eid[0]);
// 减少计数器
G->pid[S->sid[0]]--; G->pid[S->eid[0]]--;
T->pid[1] = S->pid[0];
// 收集边对象
BSL->PutData(T, T);
return 0;
}
// 增加计数器
G->pid[S->sid[0]]++; G->pid[S->eid[0]]++;
// 增加边对象
ASL->PutData(S, S);
return 1;
}
// 构造三角形:最大最小角原则
inline TCS_V00 ConTen2(TIN_G *G, TCS_F64 *P, TCS_S32 sid[2], TCS_S32 eid[2], CTcsBL *APL, CTcsBL *ATL, CTcsHL *ASL, CTcsHL *BSL)
{
TIN_T *T; TIN_S *S;
TIN_X *X1 = (TIN_X*)APL->RemoveHead();
TIN_X *X2 = (TIN_X*)APL->RemoveHead();
if (!X2) // 只有1点
{
// 增加三角形
T = new TIN_T;
T->pid[0] = sid[0], T->pid[1] = eid[0], T->pid[2] = X1->pid;
ATL->AddTail(T);
// 检测活动边
S = new TIN_S;
IniSide(S, sid[0], eid[0], sid[1], eid[1], X1->pid);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, sid[0], X1->pid, sid[1], X1->cid, eid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, eid[0], X1->pid, eid[1], X1->cid, sid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
// 删除归零点
if (G->pid[sid[0]] == 0) { RmvPid(G->idx + sid[1], sid[0]); G->vpn++; }
if (G->pid[eid[0]] == 0) { RmvPid(G->idx + eid[1], eid[0]); G->vpn++; }
if (G->pid[X1->pid] == 0) { RmvPid(G->idx + X1->cid, X1->pid); G->vpn++; }
delete X1;
return;
}
// 首先调整sid与eid的顺序
TCS_F64 *sp = P + 2 * sid[0];
TCS_F64 *ep = P + 2 * eid[0];
TCS_F64 *p1 = P + 2 * X1->pid;
TCS_F64 *p2 = P + 2 * X2->pid;
if (IsSide2(sp, p1, ep) != IsSide2(sp, p1, p2))
{
// 不同侧, 交换sid与eid
TCS_S32 pi = sid[0]; sid[0] = eid[0]; eid[0] = pi;
TCS_S32 ci = sid[1]; sid[1] = eid[1]; eid[1] = ci;
sp = P + 2 * sid[0]; ep = P + 2 * eid[0];
}
// 循环构网
while (X2)
{
// ep---------sp
//
// p1
//
// p2
// 暂不局部优化
// 增加三角形
T = new TIN_T;
T->pid[0] = sid[0], T->pid[1] = eid[0], T->pid[2] = X1->pid;
ATL->AddTail(T);
// 检测活动边
S = new TIN_S;
IniSide(S, sid[0], eid[0], sid[1], eid[1], X1->pid);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, sid[0], X1->pid, sid[1], X1->cid, eid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
// 增加三角形
T = new TIN_T;
T->pid[0] = eid[0], T->pid[1] = X1->pid, T->pid[2] = X2->pid;
ATL->AddTail(T);
// 检测活动边
S = new TIN_S;
IniSide(S, X1->pid, X2->pid, X1->cid, X2->cid, eid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, X2->pid, eid[0], X2->cid, eid[1], X1->pid);
if (!InsSide(ASL, BSL, G, S))
delete S;
// 搜集公共边
S = new TIN_S;
IniSide(S, X1->pid, eid[0], X1->cid, eid[1], sid[0]);
S->pid[1] = X2->pid; BSL->PutData(S, S);
// 删除归零点
if (G->pid[sid[0]] == 0) { RmvPid(G->idx + sid[1], sid[0]); G->vpn++; }
if (G->pid[eid[0]] == 0) { RmvPid(G->idx + eid[1], eid[0]); G->vpn++; }
if (G->pid[X1->pid] == 0) { RmvPid(G->idx + X1->cid, X1->pid); G->vpn++; }
if (G->pid[X2->pid] == 0) { RmvPid(G->idx + X2->cid, X2->pid); G->vpn++; }
// X2变sid
sid[0] = X2->pid, sid[1] = X2->cid;
delete X1; delete X2;
X1 = (TIN_X*)APL->RemoveHead();
X2 = (TIN_X*)APL->RemoveHead();
}
// 最多剩1点
if (X1)
{
// 增加三角形
T = new TIN_T;
T->pid[0] = sid[0], T->pid[1] = eid[0], T->pid[2] = X1->pid;
ATL->AddTail(T);
// 检测活动边
S = new TIN_S;
IniSide(S, sid[0], eid[0], sid[1], eid[1], X1->pid);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, sid[0], X1->pid, sid[1], X1->cid, eid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
S = new TIN_S;
IniSide(S, eid[0], X1->pid, eid[1], X1->cid, sid[0]);
if (!InsSide(ASL, BSL, G, S))
delete S;
// 删除归零点
if (G->pid[sid[0]] == 0) { RmvPid(G->idx + sid[1], sid[0]); G->vpn++; }
if (G->pid[eid[0]] == 0) { RmvPid(G->idx + eid[1], eid[0]); G->vpn++; }
if (G->pid[X1->pid] == 0) { RmvPid(G->idx + X1->cid, X1->pid); G->vpn++; }
delete X1;
}
}
// 计算搜索边界
inline TCS_V00 ComBnd2(TIN_G *G, TCS_F64 *p, TCS_F64 d, TCS_S32 &c1, TCS_S32 &c2, TCS_S32 &r1, TCS_S32 &r2)
{
TCS_F64 *bnd = G->bnd, *spn = G->spn;
TCS_S32 co = G->num[0], ro = G->num[1];
// 计算向左扩展量
c1 = TCS_S32((p[0] - d - bnd[0]) / spn[0]); if (c1 < 0) c1 = 0;
// 计算向右扩展量
c2 = TCS_S32((p[0] + d - bnd[0]) / spn[0]); if (c2 > co - 1) c2 = co - 1;
// 计算向下扩展量
r1 = TCS_S32((p[1] - d - bnd[1]) / spn[1]); if (r1 < 0) r1 = 0;
// 计算向上扩展量
r2 = TCS_S32((p[1] + d - bnd[1]) / spn[1]); if (r2 > ro - 1) r2 = ro - 1;
}
// 空圆检测:有异侧判断
// 2:包含点;1:有紧缩;0:无紧缩
inline TCS_S32 RevEmp2(TIN_G *G, TCS_F64 *P, TCS_F64 *p1, TCS_F64 *p2, TCS_S32 on, TCS_F64 **p3, TCS_F64 *ce, TCS_F64 &d2, TCS_S32 ci, CTcsBL *L)
{
TCS_S32 n, pi;
TCS_F64 *p4, dn;
TIN_C *C = G->idx + ci;
// 判断是否包含已闭合点
for (n = C->ocn; n < C->sum; n++)
{
pi = C->pid[n];
p4 = P + 2 * pi;
// 异侧判断
if (on == IsSide2(p1, p2, p4))
{
dn = EuLen2D(p4, ce);
// 精确处理浮点数相等问题
if (dn < d2)
{
TIN_X *X = (TIN_X*)L->RemoveHead();
while (X)
{
delete X;
X = (TIN_X*)L->RemoveHead();
}
return 2;
}
}
}
// 判断是否包含共圆点或紧缩
for (n = 0; n < C->ocn; n++)
{
pi = C->pid[n];
p4 = P + 2 * pi;
// 异侧判断
// 下面可以精确处理p1 == p4、p2 == p4和共线问题
if (on == IsSide2(p1, p2, p4))
{
dn = EuLen2D(p4, ce);
// 精确处理浮点数相等问题
if (dn == d2) // 共圆
{
if (p4 == *p3)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -