?? tincore.cpp
字號:
// TinCore.cpp
// 鄧雪清, 2007/09/27
#include <math.h>
#include "stdafx.h"
#include "TinCore.h"
// 歐拉距離平方
double EuLen2D(TSS_FLT64 *pt1, TSS_FLT64 *pt2);
// 計算外接圓
void CalcCircum(TSS_FLT64 *pt1, TSS_FLT64 *pt2, TSS_FLT64 *pt3, TSS_FLT64 *cen);
//#define _TIN_DEBUG_
// 測試函數
long HavePoint(TIN_TGON *tgon, IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim);
// 從格網數據構建TIN
// 1 隱含沒有重合點
// 2 vnum為三角形數目
// 3 vidx為頂點索引, 3個一組
long GridToTin2D(POINT2D *ppts, TSS_SNT32 *pnum, TSS_SNT32 *stop, TSS_SNT32 *vnum, TSS_SNT32 **vidx, TSS_TIMEB *betm, TSS_TIMEB *entm)
{
if (!ppts || *pnum <= 0)
return 0;
::_ftime_s(betm);
//////////////////////////////////
TSS_SNT32 i, k, csum, *ppid, ddim = 2, onum = *pnum;
ppid = new TSS_SNT32[onum];
if (!ppid)
return 0;
::memset(ppid, 0, sizeof(TSS_SNT32) * (onum));
//////////////////////////////////
IDX_GRID uidx;
// 計算坐標邊界
uidx.bnd[0] = uidx.bnd[2] = ppts[0].x;
uidx.bnd[1] = uidx.bnd[3] = ppts[0].y;
for (i = 1; i < onum; i++)
{
if (uidx.bnd[0] > ppts[i].x) uidx.bnd[0] = ppts[i].x;
else if (uidx.bnd[2] < ppts[i].x) uidx.bnd[2] = ppts[i].x;
if (uidx.bnd[1] > ppts[i].y) uidx.bnd[1] = ppts[i].y;
else if (uidx.bnd[3] < ppts[i].y) uidx.bnd[3] = ppts[i].y;
}
//////////////////////////////////
// 建立網格索引
uidx.dup = 0;
CalcGridSize(&uidx, onum, 1);
csum = uidx.num[0] * uidx.num[1];
if (csum == 0)
{ delete []ppid; return 0; }
uidx.idx = new IDX_CELL[csum];
if (!uidx.idx)
{ delete []ppid; return 0; }
::memset(uidx.idx, 0, sizeof(IDX_CELL) * csum);
for (i = 0; i < onum; i++)
{
k = CalcGridCell(&uidx, (TSS_FLT64*)(ppts + i));
AddPid(uidx.idx + k, i);
}
//////////////////////////////////
TSS_FLT64 dm2 = (uidx.bnd[0] - uidx.bnd[2]) * (uidx.bnd[0] - uidx.bnd[2]);
dm2+= (uidx.bnd[1] - uidx.bnd[3]) * (uidx.bnd[1] - uidx.bnd[3]);
// 構建初始三角形
TSS_SNT32 pid[3], cid[3];
*pnum = 0;
if (!MakeFirstSimplex(&uidx, (TSS_FLT64*)ppts, ddim, stop, dm2, pid, cid))
{
FreeCell(uidx.idx, csum);
delete [](uidx.idx);
delete []ppid;
return 0;
}
//////////////////////////////////
CTssList *ATL = new CTssList;
MakeSimplex(&uidx, (TSS_FLT64*)ppts, ddim, ppid, pnum, vnum, stop, pid, cid, ATL);
::_ftime_s(entm);
*vnum = ATL->GetSize();
if (*vnum > 0)
{
TIN_TGON *TGon = (TIN_TGON*)(ATL->RemoveHead());
TSS_SNT32 *TTmp = *vidx = new TSS_SNT32[*vnum * 3];
for (i = 0; i < *vnum; i++, TGon = (TIN_TGON*)(ATL->RemoveHead()))
{
k = 3 * i;
TTmp[k + 0] = TGon->pid[0];
TTmp[k + 1] = TGon->pid[1];
TTmp[k + 2] = TGon->pid[2]; delete TGon;
}
}
else
*vidx = 0;
delete ATL; delete []ppid;
FreeCell(uidx.idx, csum);
delete [](uidx.idx);
return 1;
}
// 測試是否存在重合點
/* {
long cnt = 0;
POINT3D *pt1, *pt2;
for (i = 0; i < num; i++)
{
if (ppid[i] == -1)
continue;
pt1 = pts + i;
for (k = i + 1; k < num; k++)
{
if (ppid[k] == -1)
continue;
pt2 = pts + k;
if (pt1->x - dup <= pt2->x && pt2->x <= pt1->x + dup &&
pt1->y - dup <= pt2->y && pt2->y <= pt1->y + dup)
{
ppid[k] = -1;
cnt++;
}
}
}
if (cnt > 0)
return 0;
}
*/
// 測試是否包含點
long HavePoint(CTssList *ATL, IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim)
{
TSS_SNT32 i, tnum = ATL->GetSize();
TIN_TGON *tgon = (TIN_TGON*)(ATL->GetHead());
for (i = 0; i < tnum; i++, tgon = (TIN_TGON*)(ATL->GetNext()))
{
if (HavePoint(tgon, uidx, ppts, ddim))
return 1;
}
return 0;
}
long DelaunayTriangulation2D(POINT3D *ppts, TSS_SNT32 *pnum, TSS_FLT32 ddup, TSS_SNT32 *stop,
TSS_SNT32 *vnum, TSS_SNT32 **vidx, TSS_TIMEB *betm, TSS_TIMEB *entm)
{
if (!ppts || *pnum <= 0 || ddup < 0)
return 0;
::_ftime_s(betm);
//////////////////////////////////
TSS_SNT32 i, k, csum, ddim = 3;
TSS_SNT32 *ppid = new TSS_SNT32[*pnum];
if (!ppid)
return 0;
::memset(ppid, 0, sizeof(TSS_SNT32) * (*pnum));
//////////////////////////////////
IDX_GRID uidx;
// 計算坐標邊界
uidx.bnd[0] = uidx.bnd[2] = ppts[0].x;
uidx.bnd[1] = uidx.bnd[3] = ppts[0].y;
for (i = 1; i < *pnum; i++)
{
if (uidx.bnd[0] > ppts[i].x) uidx.bnd[0] = ppts[i].x;
else if (uidx.bnd[2] < ppts[i].x) uidx.bnd[2] = ppts[i].x;
if (uidx.bnd[1] > ppts[i].y) uidx.bnd[1] = ppts[i].y;
else if (uidx.bnd[3] < ppts[i].y) uidx.bnd[3] = ppts[i].y;
}
//////////////////////////////////
// 建立網格索引
uidx.dup = ddup;
CalcGridSize(&uidx, *pnum, 1);
csum = uidx.num[0] * uidx.num[1];
if (csum == 0)
{ delete []ppid; return 0; }
uidx.idx = new IDX_CELL[csum];
if (!uidx.idx)
{ delete []ppid; return 0; }
::memset(uidx.idx, 0, sizeof(IDX_CELL) * csum);
for (i = 0; i < *pnum; i++)
{
k = CalcGridCell(&uidx, (TSS_FLT64*)(ppts + i));
AddPid(uidx.idx + k, i);
}
//////////////////////////////////
*pnum = 0;
// 消除重合點
RemoveDupA(&uidx, pnum, ddim, (TSS_FLT64*)ppts, ppid, stop);
//////////////////////////////////
TSS_FLT64 dm2 = (uidx.bnd[0] - uidx.bnd[2]) * (uidx.bnd[0] - uidx.bnd[2]);
dm2+= (uidx.bnd[1] - uidx.bnd[3]) * (uidx.bnd[1] - uidx.bnd[3]);
// 構建初始三角形
TSS_SNT32 pid[3], cid[3];
if (!MakeFirstSimplex(&uidx, (TSS_FLT64*)ppts, ddim, stop, dm2, pid, cid))
{
FreeCell(uidx.idx, csum);
delete [](uidx.idx);
delete []ppid;
return 0;
}
//////////////////////////////////
CTssList *ATL = new CTssList;
MakeSimplex(&uidx, (TSS_FLT64*)ppts, ddim, ppid, pnum, vnum, stop, pid, cid, ATL);
::_ftime_s(entm);
if (HavePoint(ATL, &uidx, (TSS_FLT64*)ppts, ddim))
return 0;
*vnum = ATL->GetSize();
if (*vnum > 0)
{
TIN_TGON *TGon = (TIN_TGON*)(ATL->RemoveHead());
TSS_SNT32 *TTmp = *vidx = new TSS_SNT32[*vnum * 6];
for (i = 0; i < *vnum; i++, TGon = (TIN_TGON*)(ATL->RemoveHead()))
{
k = i * 6;
TTmp[k + 0] = TGon->pid[0];
TTmp[k + 1] = TGon->pid[1];
TTmp[k + 2] = TGon->pid[2]; delete TGon;
}
}
else
*vidx = 0;
delete ATL; delete []ppid;
FreeCell(uidx.idx, csum);
delete [](uidx.idx);
return 1;
}
// 測試函數
long HavePoint(TIN_TGON *tgon, IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim)
{
TSS_SNT32 c[3], r[3], i, k, n, cid, nid, cnt, sum;
TSS_SNT32 col = uidx->num[0], row = uidx->num[1];
TSS_FLT64 *pt1 = ppts + tgon->pid[0] * ddim;
TSS_FLT64 *pt2 = ppts + tgon->pid[1] * ddim;
TSS_FLT64 *pt3 = ppts + tgon->pid[2] * ddim;
TSS_FLT64 cen[2], *pt4, rad;
IDX_CELL *cell, *cidx = uidx->idx;
CalcCircum(pt1, pt2, pt3, cen); rad = EuLen2D(pt1, cen);
CalcSearchBound(uidx->bnd, uidx->spn, col, row, cen, ::sqrt(rad), c[1], c[2], r[1], r[2]);
for (i = r[1]; i <= r[2]; i++)
{
cid = i * col;
for (k = c[1]; k <= c[2]; k++)
{
cell = cidx + cid + k; cnt = cell->cnt; sum = cell->sum;
for (n = cnt; n < sum; n++)
{
nid = cell->pid[n];
pt4 = ppts + ddim * nid;
if (pt4 == pt1 || pt4 == pt2 || pt4 == pt3)
continue;
if (EuLen2D(pt4, cen) < rad)
return 1;
}
}
}
return 0;
}
///////////////////////////////////////////////////////////////////
// 基本函數
TSS_SNT32 FACE_EQLH(void const *f1, void const *f2)
{
TIN_FACE const *m1 = (TIN_FACE const *)f1;
TIN_FACE const *m2 = (TIN_FACE const *)f2;
return (m1->sid[0] == m2->sid[0] && m1->eid[0] == m2->eid[0]);
}
TSS_UNT32 FACE_HASH(void const *f1)
{
TIN_FACE const *m1 = (TIN_FACE const *)f1;
return (EQHASH(m1->sid[0]) ^ EQHASH(m1->eid[0]));
}
void CalcGridSize(IDX_GRID *uidx, TSS_SNT32 num, TSS_SNT32 avg)
{
TSS_SNT32 *mat = uidx->num;
TSS_FLT32 dup = uidx->dup;
TSS_FLT64 *bnd = uidx->bnd;
TSS_FLT64 *spn = uidx->spn;
TSS_FLT64 sum = (bnd[2] - bnd[0]) * (bnd[3] - bnd[1]);
if (sum == 0.0) // 至少存在1維共線
{ mat[0] = 0; return; }
sum = sum / num * avg;
spn[0] = spn[1] = sqrt(sum);
if (dup >= spn[0]) // 便于處理重合點
spn[0] = spn[1] = 2 * dup;
mat[0] = TSS_SNT32((bnd[2] - bnd[0]) / spn[0]);
if (mat[0] * spn[0] <= bnd[2] - bnd[0])
mat[0] += 1;
mat[1] = TSS_SNT32((bnd[3] - bnd[1]) / spn[1]);
if (mat[1] * spn[1] <= bnd[3] - bnd[1])
mat[1] += 1;
}
inline long CalcGridCell(IDX_GRID *uidx, TSS_FLT64 *pnt)
{
TSS_SNT32 c = TSS_SNT32((pnt[0] - uidx->bnd[0]) / uidx->spn[0]);
TSS_SNT32 r = TSS_SNT32((pnt[1] - uidx->bnd[1]) / uidx->spn[1]);
return (r * uidx->num[0] + c);
}
inline void CalcGridCell(IDX_GRID *uidx, TSS_FLT64 *pnt, TSS_SNT32 *pos)
{
pos[0] = TSS_SNT32((pnt[0] - uidx->bnd[0]) / uidx->spn[0]);
pos[1] = TSS_SNT32((pnt[1] - uidx->bnd[1]) / uidx->spn[1]);
}
void AddPid(IDX_CELL *cell, TSS_SNT32 pid)
{
if (cell->cnt < cell->num)
{
cell->pid[cell->cnt] = pid;
cell->cnt++; cell->sum++;
}
else
{
TSS_SNT32 num = cell->cnt + 4;
TSS_SNT32 *buf = new TSS_SNT32[num];
::memcpy(buf, cell->pid, 4 * cell->cnt);
buf[cell->cnt] = pid; cell->cnt++; cell->sum++;
cell->num = num; delete []cell->pid; cell->pid = buf;
}
}
// 用于刪除重合點
inline void DelPid(IDX_CELL *cell, TSS_SNT32 pos)
{
cell->cnt--; cell->sum--;
cell->pid[pos] = cell->pid[cell->cnt];
}
// 用于刪除閉合點
inline void DelPidEx(IDX_CELL *cell, TSS_SNT32 pid)
{
TSS_SNT32 *tmp = cell->pid, cnt = cell->cnt;
for (TSS_SNT32 n = cnt - 1; n >= 0; n--)
{
if (tmp[n] == pid)
{
cnt--; tmp[n] = tmp[cnt]; tmp[cnt] = pid;
cell->cnt = cnt; return;
}
}
}
inline void SetFace(TIN_FACE *face, TSS_SNT32 sid, TSS_SNT32 eid, TSS_SNT32 pid, TSS_SNT32 tid, TSS_SNT32 sic, TSS_SNT32 eic)
{
if (sid > eid)
{
face->sid[0] = sid; face->eid[0] = eid; face->sid[1] = sic; face->eid[1] = eic;
}
else
{
face->sid[0] = eid; face->eid[0] = sid; face->sid[1] = eic; face->eid[1] = sic;
}
face->pid = pid; face->tid = tid;
}
// 檢查用2次邊
inline void SetFace(TIN_FACE *face, TSS_SNT32 sid, TSS_SNT32 eid)
{
if (sid > eid)
{
face->sid[0] = sid; face->eid[0] = eid;
}
else
{
face->sid[0] = eid; face->eid[0] = sid;
}
}
// 歐拉距離平方
inline double EuLen2D(TSS_FLT64 *p1, TSS_FLT64 *p2)
{
TSS_FLT64 dx = p1[0] - p2[0];
TSS_FLT64 dy = p1[1] - p2[1];
return (dx * dx + dy * dy);
}
// 消除重合點
void RemoveDupA(IDX_GRID *uidx, TSS_SNT32 *pnum, TSS_SNT32 ddim, TSS_FLT64 *ppts, TSS_SNT32 *ppid, TSS_SNT32 *stop)
{
// 處理四周單元, 不包括角單元
TSS_SNT32 col = uidx->num[0], row = uidx->num[1];
if (row > 2)
RemoveDupC(uidx, pnum, ddim, ppts, ppid, 0, 1, row - 2);
if (col > 1 && row > 2)
RemoveDupC(uidx, pnum, ddim, ppts, ppid, col - 1, 1, row - 2);
if (col > 2)
RemoveDupR(uidx, pnum, ddim, ppts, ppid, 0, 1, col - 2);
if (row > 1 && col > 2)
RemoveDupR(uidx, pnum, ddim, ppts, ppid, row - 1, 1, col - 2);
// 檢測停止標志
if (!(*stop))
return;
// 處理邊界單元
RemoveDupS(uidx, pnum, ddim, ppts, ppid);
// 檢測停止標志
if (!(*stop))
return;
// 處理中間單元
if (col > 2 && row > 2)
RemoveDupM(uidx, pnum, ddim, ppts, ppid, stop, 1, col - 2, 1, row - 2);
}
// 按列處理重合點
// 要求r1 >= 1, r2 <= row - 2
void RemoveDupC(IDX_GRID *uidx, TSS_SNT32 *pnum, TSS_SNT32 ddim, TSS_FLT64 *ppts, TSS_SNT32 *ppid, TSS_SNT32 c, TSS_SNT32 r1, TSS_SNT32 r2)
{
TSS_FLT32 dup = uidx->dup;
TSS_FLT64 *bnd = uidx->bnd;
TSS_FLT64 *spn = uidx->spn;
TSS_SNT32 *num = uidx->num;
IDX_CELL *cidx = uidx->idx;
TSS_SNT32 col = num[0], row = num[1];
TSS_FLT64 *pt, ty;
IDX_CELL *cell, *neig;
TSS_SNT32 r, n, cid, *pid;
ty = bnd[1] + r1 * spn[1];
for (r = r1; r <= r2; r++)
{
cid = r * col + c; ty += spn[1];
cell = cidx + cid; pid = cell->pid; neig = cell + col;
for (n = cell->cnt - 1; n >= 0; n--)
{
pt = ppts + ddim * pid[n];
if (IsHaveDup(cell, ppts, ddim, dup, pt, n - 1))
{
ppid[pid[n]] = -1; DelPid(cell, n); pnum++;
}
else
{
// 上單元
if (pt[1] + dup >= ty && IsHaveDup(neig, ppts, ddim, dup, pt, neig->cnt - 1))
{
ppid[pid[n]] = -1; DelPid(cell, n); pnum++;
}
}
}
}
}
// 按行處理重合點
// 要求c1 >= 1, c2 <= col - 2
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -