tincore.cpp

并行TIN生成算法, 基于DeWall算法理论实现

// 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