⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 tincore.cpp

📁 并行TIN生成算法, 基于DeWall算法理论实现
💻 CPP
📖 第 1 页 / 共 4 页
字号:
🔍
1234 
	}
	// 检测停止标志
	if (!(*stop))
		return 0;
	// 3.采用空圆法则搜索点pid[2]
	{
		TSS_FLT64 *pt1 = ppts + ddim * pid[0];
		TSS_FLT64 *pt2 = ppts + ddim * pid[1];
		TSS_FLT64  cen[2] = {(pt1[0] + pt2[0]) / 2, (pt1[1] + pt2[1]) / 2};

		// 计算初始搜索范围
		TSS_FLT64  rad = sqrt(EuLen2D(pt1, cen));
		CalcSearchBound(bnd, spn, col, row, cen, rad, c[1], c[2], r[1], r[2]);

		for (i = r[1]; i <= r[2]; i++)
		{
			cin = i * col;
			for (k = c[1]; k <= c[2]; k++)
			{
				cell = cidx + (cin + k);
				if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin + k, pid[2], cid[2], cen))
					return 1;
			}
		}

		// 继续搜索
		c[3] = c[1];	c[4] = c[2];		r[3] = r[1];	r[4] = r[2];
		c[5] = 0;		c[6] = col - 1;		r[5] = 0;		r[6] = row - 1;
		while (1)
		{
			c[1]--;	if (c[1] < c[5])	c[1] = c[5];
			c[2]++;	if (c[2] > c[6])	c[2] = c[6];
			r[1]--;	if (r[1] < r[5])	r[1] = r[5];
			r[2]++;	if (r[2] > r[6])	r[2] = r[6];
			if (SearchCellEMR(uidx, ppts, ddim, pt1, pt2, pid[2], cid[2], cen, col, c, r))
				return 1;	// 已经找到
			
			c[3] = c[1], c[4] = c[2], r[3] = r[1], r[4] = r[2];
			if (c[3] == c[5] && c[4] == c[6] &&
				r[3] == r[5] && r[4] == r[6])
				return 0;
		}
	}

	return 0;
}

// 单元搜索第3点(空圆法则)
inline long SearchCellEMR(IDX_GRID *uidx, IDX_CELL *cell, TSS_FLT64 *ppts, TSS_SNT32 ddim,
						  TSS_FLT64 *pt1, TSS_FLT64 *pt2, TSS_SNT32 cin, TSS_SNT32 &pid,
						  TSS_SNT32 &cid, TSS_FLT64 *cen)
{
	TSS_SNT32 n, nid;
	TSS_FLT64 *pt3, rad;

	for (n = 0; n < cell->cnt; n++)
	{
		nid = cell->pid[n];
		pt3 = ppts + ddim * nid;

		if (!IsLine(pt1, pt2, pt3))	// 不共线
		{
			CalcCircum(pt1, pt2, pt3, cen);		rad = EuLen2D(pt1, cen);	pid = nid;	cid = cin;
			SearchCellEMR(uidx, ppts, ddim, pt1, pt2, pt3, rad, pid, cid, cen);
			return 1;
		}
	}

	return 0;
}

// 判断是否共线
// 注:算法要求点集不存在重合点
inline long IsLine(TSS_FLT64 *p1, TSS_FLT64 *p2, TSS_FLT64 *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 void SearchCellEMR(IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_FLT64 *pt1,
						  TSS_FLT64 *pt2, TSS_FLT64 *pt3, TSS_FLT64 rad, TSS_SNT32 &pid,
						  TSS_SNT32 &cid, TSS_FLT64 *cen)
{
	TSS_FLT64 *pt4;
	IDX_CELL *cell, *cidx = uidx->idx;
	TSS_SNT32 i, k, n, c[4], r[4], tmp, nid;
	TSS_SNT32 col = uidx->num[0], row = uidx->num[1];

	// 计算搜索边界
	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++)
	{
		tmp = i * col; 
		for (k = c[1]; k <= c[2]; k++)
		{
			cell = cidx + (tmp + k);
			for (n = 0; n < cell->cnt; n++)
			{
				nid = cell->pid[n];
				pt4 = ppts + ddim * nid;
				// 注:rad = EuLen2D(pt1, cen)可以精确处理pt4 == pt1问题
				if (pt4 == pt3 || pt4 == pt2)// || pt4 == pt1)
					continue;

				if (EuLen2D(pt4, cen) < rad) // 有点
				{
					CalcCircum(pt1, pt2, pt4, cen);		rad = EuLen2D(pt1, cen);	pid = nid;	cid = tmp + k;
					SearchCellEMR(uidx, ppts, ddim, pt1, pt2, pt4, rad, pid, cid, cen);
					return;
				}
			}
		}
	}
}

// 批量单元搜索
inline long SearchCellEMR(IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_FLT64 *pt1,
						  TSS_FLT64 *pt2, TSS_SNT32 &pid, TSS_SNT32 &cid, TSS_FLT64 *cen,
						  TSS_SNT32 col, TSS_SNT32 *c, TSS_SNT32 *r)
{
	TSS_SNT32 i, k, cin;
	IDX_CELL *cell, *cidx = uidx->idx;

	if (r[3] != r[1])	// 即r[3] != r[5]
	{
		for (k = c[1] + 1; k < c[2]; k++)
		{
			cin = r[1] * col + k;		cell = cidx + cin;
			if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
				return 1;
		}
	}
	if (r[4] != r[2])	// 即r[4] != r[6]
	{
		for (k = c[1] + 1; k < c[2]; k++)
		{
			cin = r[2] * col + k;		cell = cidx + cin;
			if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
				return 1;
		}
	}
	if (c[3] != c[1])	// 即c[3] != c[5]
	{
		for (i = r[1] + 1; i < r[2]; i++)
		{
			cin = i * col + c[1];		cell = cidx + cin;
			if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
				return 1;
		}
	}
	if (c[4] != c[2])	// 即c[4] != c[6]
	{
		for (i = r[1] + 1; i < r[2]; i++)
		{
			cin = i * col + c[2];		cell = cidx + cin;
			if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
				return 1;
		}
	}
	// 处理4个角单元
	if (c[3] != c[5] || r[3] != r[5])
	{
		cin = r[1] * col + c[1];	cell = cidx + cin;
		if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
			return 1;
	}
	if (c[3] != c[5] || r[4] != r[6])
	{
		cin = r[2] * col + c[1];	cell = cidx + cin;
		if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
			return 1;
	}
	if (c[4] != c[6] || r[3] != r[5])
	{
		cin = r[1] * col + c[2];	cell = cidx + cin;
		if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
			return 1;
	}
	if (c[4] != c[6] || r[4] != r[6])
	{
		cin = r[2] * col + c[2];	cell = cidx + cin;
		if (SearchCellEMR(uidx, cell, ppts, ddim, pt1, pt2, cin, pid, cid, cen))
			return 1;
	}

	return 0;
}

// 计算外接圆
inline void CalcCircum(TSS_FLT64 *p1, TSS_FLT64 *p2, TSS_FLT64 *p3, TSS_FLT64 *pt)
{
	// 11(+|-), 12(*), 2(/)
	TSS_FLT64 dx21 = p2[0] - p1[0];
	TSS_FLT64 dy21 = p2[1] - p1[1];
	TSS_FLT64 xy21 = dx21 * dx21 + dy21 * dy21;

	TSS_FLT64 dx31 = p3[0] - p1[0];
	TSS_FLT64 dy31 = p3[1] - p1[1];
	TSS_FLT64 xy31 = dx31 * dx31 + dy31 * dy31;

	TSS_FLT64 top1 =   dy21 * xy31 - dy31 * xy21;
	TSS_FLT64 top2 = - dx21 * xy31 + dx31 * xy21;
	TSS_FLT64 deno =   dy21 * dx31 - dy31 * dx21;

	pt[0] = p1[0] + 0.5 * top1 / deno;
	pt[1] = p1[1] + 0.5 * top2 / deno;
}

// 构建三角形网络
long MakeSimplex(IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_SNT32 *ppid, TSS_SNT32 *pnum,
				 TSS_SNT32 *vnum, TSS_SNT32 *stop, TSS_SNT32 *pid, TSS_SNT32 *cid, CTssList *ATL)
{
	CTssHash *AFH = new CTssHash(128, FACE_EQLH, FACE_HASH);	// 工作链表
	CTssHash *BFH = new CTssHash(128, FACE_EQLH, FACE_HASH);	// 跟踪链表

	IDX_CELL *cidx = uidx->idx;

	TIN_FACE *face, *temp;
	// 构建初始三角形墙
	face = new TIN_FACE;
	SetFace(face, pid[0], pid[1], pid[2], 0, cid[0], cid[1]);
	AFH->PutData(face, face);	ppid[pid[0]]++;	ppid[pid[1]]++;

	face = new TIN_FACE;
	SetFace(face, pid[1], pid[2], pid[0], 0, cid[1], cid[2]);
	AFH->PutData(face, face);	ppid[pid[1]]++;	ppid[pid[2]]++;

	face = new TIN_FACE;
	SetFace(face, pid[2], pid[0], pid[1], 0, cid[2], cid[0]);
	AFH->PutData(face, face);	ppid[pid[2]]++;	ppid[pid[0]]++;

	TIN_TGON *tgon = new TIN_TGON;
	tgon->pid[0] = pid[0];	tgon->pid[1] = pid[1];	tgon->pid[2] = pid[2];
	tgon->tid[0] = tgon->tid[1] = tgon->tid[2] = -1;	ATL->AddTail(tgon);

	// 循环构建
	*vnum = 1;	face = (TIN_FACE*)(AFH->RemoveHead());
	while (face)
	{
		temp = 0;
		// 搜索第三点
		pid[0] = face->sid[0];	pid[1] = face->eid[0];	pid[2] = face->pid;
		cid[0] = face->sid[1];	cid[1] = face->eid[1];
		if (SearchCell(BFH, uidx, ppts, ddim, stop, pid, cid[2]))
		{
			// 加入三角形
			tgon = new TIN_TGON;
			tgon->pid[0] = pid[1];	tgon->pid[1] = pid[0];	tgon->pid[2] = pid[2];
			tgon->tid[0] = face->tid;	tgon->tid[1] = tgon->tid[2] = -1;	ATL->AddTail(tgon);
			// 减少计数器
			ppid[pid[0]]--;	ppid[pid[1]]--;
			// 收集边对象
			BFH->PutData(face, face);

			// 更新链表
			face = new TIN_FACE;
			SetFace(face, pid[0], pid[2], pid[1], *vnum, cid[0], cid[2]);
			if (InsertFace(AFH, BFH, ppid, face))
			{
				temp = face;
				face = 0;
			}

			if (!face)
				face = new TIN_FACE;
			SetFace(face, pid[2], pid[1], pid[0], *vnum, cid[2], cid[1]);
			if (InsertFace(AFH, BFH, ppid, face))
			{
				temp = face;
				face = 0;
			}

			if (face)
				delete face;
			// 删除归零点
			if (ppid[pid[0]] == 0)	{	DelPidEx(cidx + cid[0], pid[0]);	*pnum += 1;	}
			if (ppid[pid[1]] == 0)	{	DelPidEx(cidx + cid[1], pid[1]);	*pnum += 1;	}
			if (ppid[pid[2]] == 0)	{	DelPidEx(cidx + cid[2], pid[2]);	*pnum += 1;	}

			*vnum += 1;
		}
		else
		{
			// 减少计数器
			ppid[pid[0]]--;	ppid[pid[1]]--;
			// 收集边对象
			BFH->PutData(face, face);
			// 删除归零点
			if (ppid[pid[0]] == 0)	{	DelPidEx(cidx + cid[0], pid[0]);	*pnum += 1;	}
			if (ppid[pid[1]] == 0)	{	DelPidEx(cidx + cid[1], pid[1]);	*pnum += 1;	}
		}

		// 检测停止标志
		if (!(*stop))
			break;

		if (temp)
		{
			face = temp;	AFH->DelData(face);
		}
		else
			face = (TIN_FACE*)(AFH->RemoveHead());
	}

	// 释放内存
	face = (TIN_FACE*)(AFH->RemoveHead());
	while (face)
	{
		delete face;
		face = (TIN_FACE*)(AFH->RemoveNext());
	}
	face = (TIN_FACE*)(BFH->RemoveHead());
	while (face)
	{
		delete face;
		face = (TIN_FACE*)(BFH->RemoveNext());
	}
	delete AFH;	delete BFH;

	return 1;
}

// 插入边
inline long InsertFace(CTssHash *AFH, CTssHash *BFH, TSS_SNT32 *ppid, TIN_FACE *face)
{
	TIN_FACE *temp = (TIN_FACE*)(AFH->DelData(face));
	if (temp)	// 已存在
	{
		_ASSERT(temp->sid[0] == face->sid[0] && temp->eid[0] == face->eid[0]);
		// 减少计数器
		ppid[face->sid[0]]--;	ppid[face->eid[0]]--;
		// 收集边对象
		BFH->PutData(temp, temp);
		return 0;
	}

	// 增加计数器
	ppid[face->sid[0]]++;	ppid[face->eid[0]]++;
	// 增加边对象
	AFH->PutData(face, face);

	return 1;
}

// 计算同侧符号
// 注:要求点集中不存在重合点
inline long IsSide(TSS_FLT64 *p1, TSS_FLT64 *p2, TSS_FLT64 *p3)
{
	// 快速精确处理p2 == p3问题
	if (p2 == p3)	// 如果没有该步, 由于计算误差可能不能精确处理p2 == p3问题
		return 0;
	TSS_FLT64 dx21 = p2[0] - p1[0];
//	TSS_FLT64 side = (p3[1] - p1[1]) * (p2[0] - p1[0]) - (p2[1] - p1[1]) * (p3[0] - p1[0]);
	if (dx21 == 0)	// p1, p2的x坐标相等
	{
		TSS_FLT64 side = p3[0] - p1[0];
		if		(side > 0)
			return  1;
		else if (side < 0)
			return -1;
	}
	else	// side / dx21
	{
		TSS_FLT64 side = (p3[1] - p1[1]) * (p2[0] - p1[0]) - (p2[1] - p1[1]) * (p3[0] - p1[0]);// / dx21;
		if		(side > 0)
			return  1;
		else if (side < 0)
			return -1;
	}

	return 0;
}

// 定义外扩算法优化开关
#define OPTI_OU1

// 搜索第3点:空圆法则
inline long SearchCell(CTssHash *BFH, IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim,
					   TSS_SNT32 *stop, TSS_SNT32 *pid, TSS_SNT32 &cid)
{
	TSS_SNT32 cin;
	TSS_SNT32 c[8], r[8], i, k;
	IDX_CELL *cell, *cidx = uidx->idx;
	TSS_SNT32 col = uidx->num[0], row = uidx->num[1];

	TSS_FLT64 *pt1 = ppts + ddim * pid[0];
	TSS_FLT64 *pt2 = ppts + ddim * pid[1];
	TSS_FLT64 *opt = ppts + ddim * pid[2]; 

	TSS_FLT64 cen[2] = {(pt1[0] + pt2[0]) / 2, (pt1[1] + pt2[1]) / 2};
	// 计算初始搜索范围
	TSS_FLT64 rad = ::sqrt(EuLen2D(pt1, cen));
	CalcSearchBound(uidx->bnd, uidx->spn, col, row, cen, rad, c[5], c[6], r[5], r[6]);
	// 计算搜索方向
	TSS_SNT32 orn = -IsSide(pt1, pt2, opt); // 反向

	pid[2] = -1;

#ifdef OPTI_OU1

	// 采用外扩算法优化
	// 从距离圆心最近的单元开始
	TSS_SNT32 pos[2];
	CalcGridCell(uidx, cen, pos);

	if		(pos[0] < 0)		pos[0] = 0;
	else if (pos[0] > col - 1)	pos[0] = col - 1;
	if		(pos[1] < 0)		pos[1] = 0;
	else if (pos[1] > row - 1)	pos[1] = row - 1;

	k = pos[0], i = pos[1];
	cin = i * col + k;	cell = cidx + cin;
	if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin, pid, cid, cen))
		return 1;

	if (c[6] - c[5] > 0 || r[6] - r[5] > 0)
	{
		c[1] = c[2] = c[3] = c[4] = k;
		r[1] = r[2] = r[3] = r[4] = i;
		while (1)
		{
			c[1]--;	if (c[1] < c[5])	c[1] = c[5];
			c[2]++;	if (c[2] > c[6])	c[2] = c[6];
			r[1]--;	if (r[1] < r[5])	r[1] = r[5];
			r[2]++;	if (r[2] > r[6])	r[2] = r[6];
			/////////////////////////////////////////////////////////////
			if (r[3] != r[1])	// 即r[3] != r[5]
			{
				cin = r[1] * col;
				for (k = c[1] + 1; k < c[2]; k++)
				{
					cell = cidx + cin + k;
					if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin + k, pid, cid, cen))
						return 1;
				}
			}
			if (r[4] != r[2])	// 即r[4] != r[6]
			{
				cin = r[2] * col;
				for (k = c[1] + 1; k < c[2]; k++)
				{
					cell = cidx + cin + k;
					if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin + k, pid, cid, cen))
						return 1;
				}
			}
			if (c[3] != c[1])	// 即c[3] != c[5]
			{
				for (i = r[1] + 1; i < r[2]; i++)
				{
					cin = i * col + c[1];	cell = cidx + cin;
					if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin, pid, cid, cen))
						return 1;
				}
			}
			if (c[4] != c[2])	// 即c[4] != c[6]
			{
				for (i = r[1] + 1; i < r[2]; i++)
				{
					cin = i * col + c[2];	cell = cidx + cin;
					if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin, pid, cid, cen))
						return 1;
				}
			}
			// 处理4个角单元
			if (c[3] != c[5] || r[3] != r[5])
			{
				cin = r[1] * col + c[1];	cell = cidx + cin;
				if (SearchCellEMR(BFH, uidx, cell, ppts, ddim, pt1, pt2, opt, orn, cin, pid, cid, cen))
					return 1;
1234

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -