tcstan2.cpp

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

// TcsTan2.cpp
// 邓雪清 2008/11/05

#include <time.h>
#include <math.h>
#include "TcsTan2.h"

// 分区结构
typedef struct tagTIN_Z {
	TCS_S32	 sco[2];	// 左下角
	TCS_S32	 eco[2];	// 右上角
	TCS_S32	 vpn;
	TIN_G	*G;
} TIN_Z;

// C0:本单元;C1:右单元;C2:上单元;C3:右上单元
TCS_V00 RmvDup2(TIN_C *C0, TIN_C *C1, TIN_C *C2, TIN_C *C3, TCS_F64 *P, TCS_F64 dup, TCS_S32 *pid, TCS_S32 &vpn)
{
	TCS_S32 i, n;
	TCS_F64 *p1, *p2;

	// 本单元
	{
		// 注意C0->sum可能在运行中改变
		for (i = 0; i < C0->sum - 1; i++)
		{
			p1 = P + 2 * C0->pid[i];
			for (n = i + 1; n < C0->sum; n++)
			{
				p2 = P + 2 * C0->pid[n];
				if (p1[0] - dup <= p2[0] && p2[0] <= p1[0] + dup &&
					p1[1] - dup <= p2[1] && p2[1] <= p1[1] + dup)
				{
					// 建立重合点映射关系
					pid[C0->pid[n]] = pid[C0->pid[i]];
					// 消除重合点
					RmxPid(C0, n);
					// 点数减1
					vpn--;
					// 回位
					n--;
				}
			}
		}
	}
	// 右单元
	if (C1)
	{
		for (i = 0; i < C0->sum; i++)
		{
			p1 = P + 2 * C0->pid[i];
			for (n = 0; n < C1->sum; n++)
			{
				p2 = P + 2 * C1->pid[n];
				if (p1[0] - dup <= p2[0] && p2[0] <= p1[0] + dup &&
					p1[1] - dup <= p2[1] && p2[1] <= p1[1] + dup)
				{
					// 建立重合点映射关系
					pid[C1->pid[n]] = pid[C0->pid[i]];
					// 消除重合点
					RmxPid(C1, n);
					// 点数减1
					vpn--;
					// 回位
					n--;
				}
			}
		}
	}
	// 上单元
	if (C2)
	{
		for (i = 0; i < C0->sum; i++)
		{
			p1 = P + 2 * C0->pid[i];
			for (n = 0; n < C2->sum; n++)
			{
				p2 = P + 2 * C2->pid[n];
				if (p1[0] - dup <= p2[0] && p2[0] <= p1[0] + dup &&
					p1[1] - dup <= p2[1] && p2[1] <= p1[1] + dup)
				{
					// 建立重合点映射关系
					pid[C2->pid[n]] = pid[C0->pid[i]];
					// 消除重合点
					RmxPid(C2, n);
					// 点数减1
					vpn--;
					// 回位
					n--;
				}
			}
		}
	}
	// 右上单元
	if (C3)
	{
		for (i = 0; i < C0->sum; i++)
		{
			p1 = P + 2 * C0->pid[i];
			for (n = 0; n < C3->sum; n++)
			{
				p2 = P + 2 * C3->pid[n];
				if (p1[0] - dup <= p2[0] && p2[0] <= p1[0] + dup &&
					p1[1] - dup <= p2[1] && p2[1] <= p1[1] + dup)
				{
					// 建立重合点映射关系
					pid[C3->pid[n]] = pid[C0->pid[i]];
					// 消除重合点
					RmxPid(C3, n);
					// 点数减1
					vpn--;
					// 回位
					n--;
				}
			}
		}
	}
}

// 工作线程
typedef struct tagTUN_PAR {
	TCS_S32	*exit;		// 退出标志
	TCS_V00	*thrd;		// 线程句柄
	TCS_V00	*nven[2];	// 工作事件
	TCS_F64	*P;
	TIN_Z	*Z;
	TIN_R	*R;
} TUN_PAR;

TCS_V00 RmvDup2(TUN_PAR *TUN)
{
	TIN_Z *Z = TUN->Z;
	TCS_F64 *P = TUN->P;

	TCS_S32 c, r;
	TIN_C *C0, *C1, *C2, *C3;

	TIN_G *G = Z->G;
	TCS_S32 co = G->num[0], ro = G->num[1];

	for (r = Z->sco[1]; r <= Z->eco[1]; r++)
	{
		C0 = G->idx + r * co + Z->sco[0];
		for (c = Z->sco[0]; c <= Z->eco[0]; c++)
		{
			C1 = C0 + 1;
			C2 = C0 + co;
			C3 = C2 + 1;
			RmvDup2(C0, C1, C2, C3, P, G->dup, G->pid, Z->vpn);
			// 指针推进1列
			C0 = C0 + 1;
		}
	}
}

TCS_U32 WINAPI TunThread(LPVOID lpParam)
{
	TUN_PAR *TUN = (TUN_PAR*)lpParam;

	while (1)
	{
		::WaitForSingleObject(TUN->nven[0], INFINITE);
		if (*(TUN->exit) == 1)
			break;

		// 消除重合
		if (TUN->Z)
			RmvDup2(TUN);
		TUN->Z = 0;
		// 构网分区
		if (TUN->R && TUN->R->S)
		{
			TIN_R *R = TUN->R;
			if (R->AXI == 0)
				ConTxn2(R->G, R->S, TUN->P, R->AXV, R->ATL, R->ASL, R->BSL);
			else
				ConTyn2(R->G, R->S, TUN->P, R->AXV, R->ATL, R->ASL, R->BSL);
		}
		TUN->R = 0;

		::SetEvent(TUN->nven[1]);
		if (*(TUN->exit) == 1)
			break;
	}

	::SetEvent(TUN->nven[1]);
	TUN->thrd = 0;
	TUN->Z = 0;
	TUN->R = 0;

	return 0;
}

// 调度线程
typedef struct tagTAN_PAR {
	TCS_S32	*exit;		// 退出标志
	TCS_V00	*thrd;		// 线程句柄
	TCS_S32	 unnu;		// 线程数目
	TCS_V00	*unpa;		// 线程参数
	TIN_G	*G;
	TCS_F64	*P;
	TCS_S32	*vpn;
	TCS_S32	*nzn;
	TCS_TMB	*btm;
	TCS_TMB	*etm;
	CTcsBL	*ATL;		// 三角形链表
	CTcsHL	*BSL;		// 静止边链表
	CTcsBL	*ARL;		// 隔离区链表
} TAN_PAR;

// 1分4
TCS_V00 To4(TIN_Z *Z1, TIN_Z *Z2, TIN_Z *Z3, TIN_Z *Z4, TIN_G *G, TCS_F64 *P)
{
	TCS_S32 c = (Z1->sco[0] + Z1->eco[0]) / 2;
	TCS_S32 r = (Z1->sco[1] + Z1->eco[1]) / 2;

	// 消除隔断重合
	TCS_S32 i, k;
	TIN_C *C0, *C1, *C2, *C3;

	TCS_S32 co = G->num[0];
	TCS_S32 ro = G->num[1];
	// 处理c列网格单元
	C0 = G->idx + Z1->sco[1] * co + c;
	for (i = Z1->sco[1]; i <= Z1->eco[1]; i++)
	{
		C1 = C0 + 1;
		C2 = C0 + co;
		C3 = C2 + 1;
		RmvDup2(C0, C1, C2, C3, P, G->dup, G->pid, Z1->vpn);
		// 指针推进1行
		C0 = C0 + co;
	}
	// 处理r行网格单元
	C0 = G->idx + r * co + Z1->sco[0];
	for (k = Z1->sco[0]; k <= Z1->eco[0]; k++)
	{
		C1 = C0 + 1;
		C2 = C0 + co;
		C3 = C2 + 1;
		RmvDup2(C0, C1, C2, C3, P, G->dup, G->pid, Z1->vpn);
		// 指针推进1列
		C0 = C0 + 1;
	}

	// 构成区间
	// 右下区间
	Z2->sco[0] = c + 1;			Z2->eco[0] = Z1->eco[0];
	Z2->sco[1] = Z1->sco[1];	Z2->eco[1] = r - 1;
	Z2->vpn = 0;	Z2->G = G;
	if (Z2->sco[0] >= co - 1)	Z2->sco[0] = co - 2;
	if (Z2->eco[1] < 0)			Z2->eco[1] = 0;

	// 左上区间
	Z3->sco[0] = Z1->sco[0];	Z3->eco[0] = c - 1;
	Z3->sco[1] = r + 1;			Z3->eco[1] = Z1->eco[1];
	Z3->vpn = 0;	Z3->G = G;
	if (Z3->eco[0] < 0)			Z3->eco[0] = 0;
	if (Z3->sco[1] >= ro - 1)	Z3->sco[1] = ro - 2;

	// 右上区间
	Z4->sco[0] = c + 1;			Z4->eco[0] = Z1->eco[0];
	Z4->sco[1] = r + 1;			Z4->eco[1] = Z1->eco[1];
	Z4->vpn = 0;	Z4->G = G;
	if (Z4->sco[0] >= co - 1)	Z4->sco[0] = co - 2;
	if (Z4->sco[1] >= ro - 1)	Z4->sco[1] = ro - 2;

	// 左下区间
	Z1->eco[0] = c - 1;
	Z1->eco[1] = r - 1;
	if (Z1->eco[0] < 0)	Z1->eco[0] = 0;
	if (Z1->eco[1] < 0)	Z1->eco[1] = 0;
}

// 消除重合点
TCS_V00 RmvDup2(TAN_PAR *TAN)
{
	TIN_G *G = TAN->G;
	TCS_F64 *P = TAN->P;
	TCS_S32 *vpn = TAN->vpn;
	TCS_S32 *nzn = TAN->nzn;

	TCS_S32 c, r;
	TIN_C *C0, *C1, *C2, *C3;

	TCS_S32 co = G->num[0];
	TCS_S32 ro = G->num[1];
	// 处理最右边网格单元
	C1 = C3 = 0;
	C0 = G->idx + co - 1;
	for (r = 0; r < ro; r++)
	{
		if (r == ro - 1)
			C2 = 0;
		else
			C2 = C0 + co;
		RmvDup2(C0, C1, C2, C3, P, G->dup, G->pid, *vpn);
		// 指针推进1行
		C0 = C0 + co;
	}
	// 处理最上边网格单元
	C2 = C3 = 0;
	C0 = G->idx + (ro - 1) * co;
	for (c = 0; c < co; c++)
	{
		if (c == co - 1)
			C1 = 0;
		else
			C1 = C0 + 1;
		RmvDup2(C0, C1, C2, C3, P, G->dup, G->pid, *vpn);
		// 指针推进1列
		C0 = C0 + 1;
	}
	// 处理可能情况
	if (co == 1 || ro == 1)
		return;
	// 按点数进行分区
	TCS_S32 zn = *vpn / *nzn;
	if (zn == 0)	zn = 1;
	// 采用4分法, 即1分4
	TCS_S32 cz = 0;
	CTcsSL *S1 = new CTcsSL;
	CTcsSL *S3, *S2 = new CTcsSL;

	TIN_Z *Z1 = new TIN_Z;
	Z1->sco[0] = 0;	Z1->eco[0] = co - 2;
	Z1->sco[1] = 0;	Z1->eco[1] = ro - 2;
	Z1->vpn = 0;	Z1->G = G;	S1->PutData(Z1);
	while (S1->GetSize() + cz < zn)
	{
		Z1 = (TIN_Z*)S1->PopData();
		while (Z1)
		{
			// 1分4
			TIN_Z *Z2 = new TIN_Z;
			TIN_Z *Z3 = new TIN_Z;
			TIN_Z *Z4 = new TIN_Z;
			To4(Z1, Z2, Z3, Z4, G, P);
			if (Z1->eco[0] >= Z1->sco[0] &&	Z1->eco[1] >= Z1->sco[1])
				S2->PutData(Z1);
			else
			{	delete Z1;	*vpn += Z1->vpn;	cz++;	}
			if (Z2->eco[0] >= Z2->sco[0] &&	Z2->eco[1] >= Z2->sco[1])
				S2->PutData(Z2);
			else
			{	delete Z2;	*vpn += Z2->vpn;	cz++;	}
			if (Z3->eco[0] >= Z3->sco[0] &&	Z3->eco[1] >= Z3->sco[1])
				S2->PutData(Z3);
			else
			{	delete Z3;	*vpn += Z3->vpn;	cz++;	}
			if (Z4->eco[0] >= Z4->sco[0] &&	Z4->eco[1] >= Z4->sco[1])
				S2->PutData(Z4);
			else
			{	delete Z4;	*vpn += Z4->vpn;	cz++;	}

			// 下一个
			Z1 = (TIN_Z*)S1->PopData();
		}
		// 交换链表
		S3 = S1;	S1 = S2;	S2 = S3;
	}

	//调度线程
	TCS_V00 **nven = new TCS_V00*[TAN->unnu];
	TUN_PAR *unpa = (TUN_PAR*)TAN->unpa;
	for (long i = 0; i < TAN->unnu; i++)
		nven[i] = unpa[i].nven[1];
	Z1 = (TIN_Z*)S1->PopData();
	while (Z1)
	{
		// 等待线程空闲
		TCS_U32 idx = ::WaitForMultipleObjects(TAN->unnu, nven, 0, INFINITE);
		if (idx >= WAIT_OBJECT_0 && idx < WAIT_OBJECT_0 + TAN->unnu)
		{
			idx = idx - WAIT_OBJECT_0;
			unpa[idx].P = P;
			unpa[idx].Z = Z1;
			::SetEvent(unpa[idx].nven[0]);
		}
		S2->PutData(Z1);
		// 如果中途终止
		if (*(TAN->exit) == 1)
			break;

		Z1 = (TIN_Z*)S1->PopData();
	}
	delete []nven;
	// 等待任务完成
	for (long i = 0; i < TAN->unnu; i++)
	{
		while (1)
		{
			if (unpa[i].Z == 0 || unpa[i].thrd == 0)
				break;
			else
				::Sleep(2);
		}
	}

	// 清理内存
	Z1 = (TIN_Z*)S1->PopData();
	while (Z1)
	{
		*vpn += Z1->vpn;	delete Z1;
		Z1 = (TIN_Z*)S1->PopData();
	}
	Z1 = (TIN_Z*)S2->PopData();
	while (Z1)
	{
		*vpn += Z1->vpn;	delete Z1;
		Z1 = (TIN_Z*)S2->PopData();
	}
	delete S1;	delete S2;
}

// 搜索最邻近点
TCS_V00 FndNep2(TIN_G *G, TIN_C *C, TCS_F64 *P, TCS_F64 *p, TCS_S32 ci, TCS_F64 &dnp, TCS_S32 &pid, TCS_S32 &cid)
{
	TCS_S32 n, pi;
	TCS_F64 *p2, d2;

	for (n = 0; n < C->ocn; n++)
	{
		pi = C->pid[n];
		p2 = P + 2 * pi;

		d2 = EuLen2D(p, p2);
		if (d2 < dnp)
		{	pid = pi;	cid = ci;	dnp = d2;	}
	}
}

TCS_V00 FndNep2(TIN_G *G, TCS_F64 *P, TCS_F64 *p, TCS_S32 *c, TCS_S32 *r, TCS_S32 co, TCS_F64 &dnp, TCS_S32 &pid, TCS_S32 &cid)
{
	TIN_C *C;
	TCS_S32 i, k, ci;

	if (r[2] != r[0])	// 即r[2] != r[4]
	{
		for (k = c[0] + 1; k < c[1]; k++)
		{
			ci = r[0] * co + k;	C = G->idx + ci;
			FndNep2(G, C, P, p, ci, dnp, pid, cid);
		}
	}
	if (r[3] != r[1])	// 即r[3] != r[5]
	{
		for (k = c[0] + 1; k < c[1]; k++)
		{
			ci = r[1] * co + k;	C = G->idx + ci;
			FndNep2(G, C, P, p, ci, dnp, pid, cid);
		}
	}
	if (c[2] != c[0])	// 即c[2] != c[4]
	{
		for (i = r[0] + 1; i < r[1]; i++)
		{
			ci = i * co + c[0];	C = G->idx + ci;
			FndNep2(G, C, P, p, ci, dnp, pid, cid);
		}
	}
	if (c[3] != c[1])	// 即c[3] != c[5]
	{
		for (i = r[0] + 1; i < r[1]; i++)
		{
			ci = i * co + c[1];	C = G->idx + ci;
			FndNep2(G, C, P, p, ci, dnp, pid, cid);
		}
	}
	// 处理4个角单元
	if (c[2] != c[4] || r[2] != r[4])
	{
		ci = r[0] * co + c[0];	C = G->idx + ci;
		FndNep2(G, C, P, p, ci, dnp, pid, cid);
	}
	if (c[2] != c[4] || r[3] != r[5])
	{
		ci = r[1] * co + c[0];	C = G->idx + ci;
		FndNep2(G, C, P, p, ci, dnp, pid, cid);
	}
	if (c[3] != c[5] || r[2] != r[4])
	{
		ci = r[0] * co + c[1];	C = G->idx + ci;
		FndNep2(G, C, P, p, ci, dnp, pid, cid);
	}
	if (c[3] != c[5] || r[3] != r[5])
	{
		ci = r[1] * co + c[1];	C = G->idx + ci;
		FndNep2(G, C, P, p, ci, dnp, pid, cid);
	}
}

TCS_V00 FndNep2(TIN_G *G, TCS_F64 *P, TCS_F64 *p, TCS_S32 &pid, TCS_S32 &cid)
{
	pid = cid = -1;

	// 计算单元位置
	TCS_S32 o[2];
	ComIdx2(G, p, o);

	// 设置计算变量
	TCS_F64 *p2, d2;
	TCS_S32 n, pi, ci;
	TCS_S32 co = G->num[0];
	TCS_S32 ro = G->num[1];
	TCS_F64 dx = G->bnd[2] - G->bnd[0];
	TCS_F64 dy = G->bnd[3] - G->bnd[1];
	TCS_F64 dnp = (dx * dx + dy * dy) / 4;
	// 在自身单元搜索
	ci = o[1] * co + o[0];
	TIN_C *C = G->idx + ci;
	for (n = 0; n < C->ocn; n++)
	{
		pi = C->pid[n];
		p2 = P + 2 * pi;
		if (p2 == p) // 同点
			continue;

		d2 = EuLen2D(p, p2);
		if (d2 < dnp)
		{	pid = pi;	cid = ci;	dnp = d2;	}
	}

	// 扩展搜索算法
	TCS_S32 c[6], r[6];

	c[0] = c[1] = o[0];					r[0] = r[1] = o[1];
	c[2] = c[0];	c[3] = c[1];		r[2] = r[0];	r[3] = r[1];
	c[4] = 0;		c[5] = co - 1;		r[4] = 0;		r[5] = ro - 1;
	// 在其他单元搜索
	if (pid == -1)
	{
		while (1)
		{
			c[0]--;	if (c[0] < c[4])	c[0] = c[4];
			c[1]++;	if (c[1] > c[5])	c[1] = c[5];
			r[0]--;	if (r[0] < r[4])	r[0] = r[4];
			r[1]++;	if (r[1] > r[5])	r[1] = r[5];
			FndNep2(G, P, p, c, r, co, dnp, pid, cid);

			c[2] = c[0], c[3] = c[1], r[2] = r[0], r[3] = r[1];
			if (c[2] == c[4] && c[3] == c[5] &&
				r[2] == r[4] && r[3] == r[5])
				break;
			if (pid != -1)	// 已经找到候选点
				break;
		}
	}
	// 精化搜索
	if (pid != -1)
	{
		// 计算最后搜索范围
		ComBnd2(G, p, sqrt(dnp), c[4], c[5], r[4], r[5]);
		// 根据搜索边界判断是否需要继续搜索
		while (1)
		{
			c[0]--;	if (c[0] < c[4])	c[0] = c[2] = c[4];	// 说明前一列已经到达搜索边界
			c[1]++;	if (c[1] > c[5])	c[1] = c[3]	= c[5];
			r[0]--;	if (r[0] < r[4])	r[0] = r[2] = r[4];	// 说明前一行已经到达搜索边界
			r[1]++;	if (r[1] > r[5])	r[1] = r[3] = r[5];
			FndNep2(G, P, p, c, r, co, dnp, pid, cid);

			c[2] = c[0], c[3] = c[1], r[2] = r[0], r[3] = r[1];
			if (c[2] == c[4] && c[3] == c[5] &&
				r[2] == r[4] && r[3] == r[5])
				break;
		}
	}
}

// 空圆检测和紧缩搜索
TCS_V00 RevEmp2(TIN_G *G, TCS_F64 *P, TCS_F64 *p1, TCS_F64 *p2, TCS_F64 *p3, TCS_F64 *ce, TCS_F64 &d2, CTcsBL *L)
{
	TIN_C *C;
	TCS_F64 *p4, dn;
	TCS_S32 n, ci, pi, c[2], r[2];
	TCS_S32 i, k, co = G->num[0], ro = G->num[1];

	// 计算搜索边界
	ComBnd2(G, ce, sqrt(d2), c[0], c[1], r[0], r[1]);

	for (i = r[0]; i <= r[1]; i++)
	{
		ci = i * co; 
		for (k = c[0]; k <= c[1]; k++)
		{
			C = G->idx + ci + k;
			for (n = 0; n < C->ocn; n++)
			{
				pi = C->pid[n];
				p4 = P + 2 * pi;

				dn = EuLen2D(p4, ce);
				// 精确处理浮点数相等问题
				if		(dn == d2) // 共圆
				{
					if (p4 == p1 || p4 == p2 || p4 == p3)
						continue;
					TIN_X *X = new TIN_X;
					X->cid = ci + k;
					X->pid = pi;
					L->AddTail(X);
				}
				else if (dn <  d2 && p4 != p2 && p4 != p3) // 有点
				{
					// 紧缩搜索
					TIN_X *X = (TIN_X*)L->GetHead();
					X->cid = ci + k;
					X->pid = pi;
					// 清除其它候选
					X = (TIN_X*)L->RemoveAt(1);
					while (X)
					{
						delete X;
						X = (TIN_X*)L->RemoveAt(1);
					}
					// 计算外接圆圆心
					ComCen2(p1, p2, p4, ce);
					// 计算外接圆半径平方
					d2 = EuLen2D(p1, ce);
					// 空圆检测和紧缩搜索
					RevEmp2(G, P, p1, p2, p4, ce, d2, L);

					return;
				}
			}
		}
	}
}

// 搜索第3点(空圆法则)
TCS_V00 FndEmp2(TIN_G *G, TCS_F64 *P, TCS_F64 *p1, TCS_F64 *p2, TCS_S32 ci, TCS_F64 *ce, TCS_F64 &d2, CTcsBL *L)
{
	TCS_F64 *p3;
	TCS_S32 n, pi;
	TIN_C *C = G->idx + ci;

	for (n = 0; n < C->ocn; n++)
	{
		pi = C->pid[n];
		p3 = P + 2 * pi;

		if (!IsLine2(p1, p2, p3))	// 不共线
		{
			TIN_X *X = new TIN_X;
			X->cid = ci;
			X->pid = pi;
			L->AddTail(X);
			// 计算外接圆圆心
			ComCen2(p1, p2, p3, ce);
			// 计算外接圆半径平方
			d2 = EuLen2D(p1, ce);
			// 空圆检测和紧缩搜索
			RevEmp2(G, P, p1, p2, p3, ce, d2, L);