tincore1.cpp

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

// TinCore.cpp
// 邓雪清, 2007/09/27

#include <math.h>
#include "stdafx.h"
#include "TinCore.h"

//#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);
	//////////////////////////////////
	long i, k, ddim = 2, onum = *pnum;

	long *ppid = new long[onum];
	if (!ppid)
		return 0;
	memset(ppid, 0, sizeof(long) * (onum));

	IDX_GRID *uidx = new IDX_GRID;
	//////////////////////////////////
	// 计算坐标边界
	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);

	long csum = uidx->num[0] * uidx->num[1];
	if (csum == 0)
	{
		delete []ppid;	delete uidx;
		return 0;
	}
	uidx->idx = new IDX_CELL[csum];
	if (!uidx->idx)
	{
		delete []ppid;	delete uidx;
		return 0;
	}
	memset(uidx->idx, 0, sizeof(IDX_CELL) * csum);

	for (i = 0; i < onum; i++)
	{
		k = CalcGridCell(uidx, (double*)(ppts + i));
		AddPid(uidx->idx + k, i);
	}
	//////////////////////////////////
	double	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, (double*)ppts, ddim, stop, dm2, pid, cid))
	{
		delete []ppid;	FreeCell(uidx->idx, csum);
		delete [](uidx->idx);	delete uidx;
		return 0;
	}
	//////////////////////////////////
	CTssList *ATL = new CTssList;

	MakeSimplex(uidx, (double*)ppts, ddim, ppid, pnum, vnum, stop, pid, cid, ATL);

	::_ftime_s(entm);

	*vnum = ATL->GetSize();

	if (*vnum > 0)
	{
		long *Tmp = *vidx = new long[*vnum * 3];
		TIN_TGON *TGon = (TIN_TGON*)(ATL->RemoveHead());
		for (i = 0; i < *vnum; i++, TGon = (TIN_TGON*)(ATL->RemoveHead()))
		{
			k = 3 * i;
			Tmp[k + 0] = TGon->pid[0];
			Tmp[k + 1] = TGon->pid[1];
			Tmp[k + 2] = TGon->pid[2];	delete TGon;
		}
	}
	else
		*vidx = 0;

	delete ATL;	delete []ppid;	FreeCell(uidx->idx, csum);
	delete [](uidx->idx);	delete uidx;

	return 1;
}

long DelaunayTriangulation2D(POINT3D *ppts, TSS_SNT32 *pnum, TSS_FLT32 ddup, TSS_SNT32 *stop, TSS_SNT32 *tnum, TSS_SNT32 **tset, TSS_TIMEB *betm, TSS_TIMEB *entm)
{
	if (!ppts || *pnum <= 0 || ddup < 0)
		return 0;

	::_ftime_s(betm);
	//////////////////////////////////
	long ddim = 3;

	long *ppid = new long[*pnum];
	if (!ppid)
		return 0;
	memset(ppid, 0, sizeof(long) * (*pnum));

	IDX_GRID *uidx = new IDX_GRID;
	//////////////////////////////////
	// 计算坐标边界
	uidx->bnd[0] = uidx->bnd[2] = ppts[0].x;
	uidx->bnd[1] = uidx->bnd[3] = ppts[0].y;
	for (long 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);

	long csum = uidx->num[0] * uidx->num[1];
	if (csum == 0)
	{
		delete []ppid;	delete uidx;
		return 0;
	}
	uidx->idx = new IDX_CELL[csum];
	if (!uidx->idx)
	{
		delete []ppid;	delete uidx;
		return 0;
	}
	memset(uidx->idx, 0, sizeof(IDX_CELL) * csum);

	long i, k;
	for (i = 0; i < *pnum; i++)
	{
		k = CalcGridCell(uidx, (double*)(ppts + i));
		AddPid(uidx->idx + k, i);
	}
	//////////////////////////////////
	*pnum = 0;
	// 消除重合点
	RemoveDupA(uidx, pnum, ddim, (double*)ppts, ppid, stop);
	// 测试
	{
/*		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;
*/	}
	//////////////////////////////////
	double	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, (double*)ppts, ddim, stop, dm2, pid, cid))
	{
		delete []ppid;	FreeCell(uidx->idx, csum);	delete [](uidx->idx);	delete uidx;
		return 0;
	}
	//////////////////////////////////
	CTssList *ATL = new CTssList;

	MakeSimplex(uidx, (double*)ppts, ddim, ppid, pnum, tnum, stop, pid, cid, ATL);

	::_ftime_s(entm);

	*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, (double*)ppts, ddim))
			{
				tgon = (TIN_TGON*)(ATL->RemoveHead());
				while (tgon)
				{
					delete tgon;
					tgon = (TIN_TGON*)(ATL->RemoveHead());
				}
				delete ATL;
				delete []ppid;	FreeCell(uidx->idx, csum);	delete [](uidx->idx);	delete uidx;
				return 0;
			}
		}*/
	}

	if (*tnum > 0)
	{
		long *tmp = new long[*tnum * 6];
		TIN_TGON *tgon = (TIN_TGON*)(ATL->RemoveHead());
		for (i = 0; i < *tnum; i++, tgon = (TIN_TGON*)(ATL->RemoveHead()))
		{
			tmp[i * 6 + 0] = tgon->pid[0];
			tmp[i * 6 + 1] = tgon->pid[1];
			tmp[i * 6 + 2] = tgon->pid[2];	delete tgon;
		}
		
		*tset = tmp;
	}
	else
		*tset = 0;

	delete ATL;

	delete []ppid;	FreeCell(uidx->idx, csum);	delete [](uidx->idx);	delete uidx;

	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)
{
#ifdef _TIN_DEBUG_
			if (pid == 350 || pid == 354)
				TSS_SNT32 test = pid;
#endif
	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))