tincore.cpp

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

void RemoveDupR(IDX_GRID *uidx, TSS_SNT32 *pnum, TSS_SNT32 ddim, TSS_FLT64 *ppts, TSS_SNT32 *ppid, TSS_SNT32 r, TSS_SNT32 c1, TSS_SNT32 c2)
{
	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, rx;
	IDX_CELL *cell, *neig;
	TSS_SNT32 c, n, cid, *pid;

	rx = bnd[0] + c1 * spn[0];
	for (c = c1; c <= c2; c++)
	{
		cid = r * col + c;		rx += spn[0];
		cell = cidx + cid;		pid = cell->pid;	neig = cell + 1;
		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[0] + dup >= rx && IsHaveDup(neig, ppts, ddim, dup, pt, neig->cnt - 1))
				{
					ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
				}
			}
		}
	}
}

// 处理边界单元
void RemoveDupS(IDX_GRID *uidx, TSS_SNT32 *pnum, TSS_SNT32 ddim, TSS_FLT64 *ppts, TSS_SNT32 *ppid)
{
	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, rx, ty;
	TSS_SNT32 c, r, n, cid, *pid, hav;
	IDX_CELL *cell, *nei1, *nei2, *nei3;

	// 处理四个角单元
	rx = bnd[0] + spn[0];	ty = bnd[1] + spn[1];
	// 左下单元
	cell = cidx;	nei1 = cell + 1;	nei2 = cell + col;	nei3 = nei2 + 1;	pid = cell->pid;
	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
		{
			hav = 0;
			// 右单元
			if (col > 1 && pt[0] + dup >= rx)
			{
				if (IsHaveDup(nei1, ppts, ddim, dup, pt, nei1->cnt - 1))
					hav = 1;
				// 右上单元
				if (!hav && row > 1 && pt[1] + dup >= ty && IsHaveDup(nei3, ppts, ddim, dup, pt, nei3->cnt - 1))
					hav = 1;
			}
			// 上单元
			if (!hav && row > 1 && pt[1] + dup >= ty && IsHaveDup(nei2, ppts, ddim, dup, pt, nei2->cnt - 1))
				hav = 1;
			if (hav)
			{
				ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
			}
		}
	}
	// 右下单元
	cell = cidx + col - 1;	nei2 = cell + col;	pid = cell->pid;
	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 (row > 1 && pt[1] + dup >= ty && IsHaveDup(nei2, ppts, ddim, dup, pt, nei2->cnt - 1))
			{
				ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
			}
		}
	}
	// 左上单元
	cell = cidx + (row - 1) * col;	nei1 = cell + 1;	pid = cell->pid;
	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 (col > 1 && pt[0] + dup >= rx && IsHaveDup(nei1, ppts, ddim, dup, pt, nei1->cnt - 1))
			{
				ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
			}
		}
	}
	// 右上单元
	cell = cidx + row * col - 1;	pid = cell->pid;
	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++;
		}
	}
	
	// 处理下边单元边界点
	// 要求r = 0, c >= 1, c <= col - 2所有单元上边界点
	if (row > 1)
	{
		rx = bnd[0] + spn[0];	ty = bnd[1] + spn[1];
		for (c = 1; c <= col - 2; c++)
		{
			rx += spn[0];
			cell = cidx + c;	pid = cell->pid;
			nei2 = cell + col;	nei3 = nei2 + 1;
			for (n = cell->cnt - 1; n >= 0; n--)
			{
				pt = ppts + ddim * pid[n];
				// 上单元
				if (pt[1] + dup >= ty)
				{
					if (IsHaveDup(nei2, ppts, ddim, dup, pt, nei2->cnt - 1))
					{
						ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
					}
					else if (pt[0] + dup >= rx && IsHaveDup(nei3, ppts, ddim, dup, pt, nei3->cnt - 1))
					{
						ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
					}
				}
			}
		}
	}
	// 要求c = 0, r >= 1, r <= row - 2所有单元上边界点
	if (col > 1)
	{
		rx = bnd[0] + spn[0];	ty = bnd[1] + spn[1];
		for (r = 1; r <= row - 2; r++)
		{
			cid = r * col;	ty += spn[1];
			cell = cidx + cid;	pid = cell->pid;
			nei1 = cell + 1;	nei3 = nei1 + col;
			for (n = cell->cnt - 1; n >= 0; n--)
			{
				pt = ppts + ddim * pid[n];
				// 右单元
				if (pt[0] + dup >= rx)
				{
					if (IsHaveDup(nei1, ppts, ddim, dup, pt, nei1->cnt - 1))
					{
						ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
					}
					else if (pt[1] + dup >= ty && IsHaveDup(nei3, ppts, ddim, dup, pt, nei3->cnt - 1))
					{
						ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
					}
				}
			}
		}
	}
}

// 按矩阵处理重合点
// 要求c1 >= 1, c2 <= col - 2, r1 >= 1, r2 <= row - 2
void RemoveDupM(IDX_GRID *uidx, TSS_SNT32 *pnum, TSS_SNT32 ddim, TSS_FLT64 *ppts, TSS_SNT32 *ppid, TSS_SNT32 *stop, TSS_SNT32 c1, TSS_SNT32 c2, 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, rx, ty;
	TSS_SNT32 c, r, n, cid, *pid, hav;
	IDX_CELL *cell, *nei1, *nei2, *nei3;

	ty = bnd[1] + r1 * spn[1];
	for (r = r1; r <= r2; r++)
	{
		cid = r * col;		ty += spn[1];		rx = bnd[0] + c1 * spn[0];
		for (c = c1; c <= c2; c++)
		{
			rx += spn[0];
			cell = cidx + cid + c;		pid = cell->pid;
			// 右单元, 上单元, 右上单元
			nei1 = cell + 1;	nei2 = cell + col;	nei3 = nei2 + 1;
			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
				{
					// 必要时处理右单元, 上单元, 右上单元
					hav = 0;
					// 右单元
					if (pt[0] + dup >= rx)
					{
						if (IsHaveDup(nei1, ppts, ddim, dup, pt, nei1->cnt - 1))
							hav = 1;
						// 右上单元
						if (!hav && pt[1] + dup >= ty && IsHaveDup(nei3, ppts, ddim, dup, pt, nei3->cnt - 1))
							hav = 1;
					}
					// 上单元
					if (!hav && pt[1] + dup >= ty && IsHaveDup(nei2, ppts, ddim, dup, pt, nei2->cnt - 1))
						hav = 1;
					// 是否有重合点
					if (hav)
					{
						ppid[pid[n]] = -1;	DelPid(cell, n);	pnum++;
					}
				}
			}
		}
		// 检测停止标志
		if (!(*stop))
			return;
	}
}

// 判断单元中是否有重合点
inline long IsHaveDup(IDX_CELL *cell, TSS_FLT64 *ppts, TSS_SNT32 ddim,
					  TSS_FLT32 dup, TSS_FLT64 *pt, TSS_SNT32 s)
{
	TSS_FLT64 *xy;
	TSS_SNT32 *pid = cell->pid;

	for (TSS_SNT32 n = s; n >= 0; n--)
	{
		xy = ppts + ddim * pid[n];
		if (xy[0] - dup <= pt[0] && pt[0] <= xy[0] + dup &&
			xy[1] - dup <= pt[1] && pt[1] <= xy[1] + dup)
			return 1;
	}

	return 0;
}

// 释放CELL索引
inline void FreeCell(IDX_CELL *cidx, TSS_SNT32 num)
{
	for (TSS_SNT32 i = 0; i < num; i++)
	{
		if (cidx[i].num > 0)
			delete []cidx[i].pid;
	}
}

// 计算搜索边界
// 要求pnt在某个单元中
inline void CalcSearchBound(TSS_FLT64 *bnd, TSS_FLT64 *spn, TSS_SNT32 col, TSS_SNT32 row,
							TSS_FLT64 *pnt, TSS_FLT64 rad, TSS_SNT32 &c1, TSS_SNT32 &c2,
							TSS_SNT32 &r1, TSS_SNT32 &r2)
{
	// 计算向左扩展量
	c1 = TSS_SNT32((pnt[0] - rad - bnd[0]) / spn[0]);	if (c1 < 0)			c1 = 0;
	// 计算向右扩展量
	c2 = TSS_SNT32((pnt[0] + rad - bnd[0]) / spn[0]);	if (c2 > col - 1)	c2 = col - 1;
	// 计算向下扩展量
	r1 = TSS_SNT32((pnt[1] - rad - bnd[1]) / spn[1]);	if (r1 < 0)			r1 = 0;
	// 计算向上扩展量
	r2 = TSS_SNT32((pnt[1] + rad - bnd[1]) / spn[1]);	if (r2 > row - 1)	r2 = row - 1;
}

// 单元搜索最近点
// 同单元搜索
inline void SearchCellNPS(IDX_CELL *cell, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_FLT64 *pnt,
						  TSS_SNT32 cin, TSS_SNT32 &pid, TSS_SNT32 &cid, TSS_FLT64 &dnp)
{
	TSS_SNT32 n, nid;
	TSS_FLT64 *tmp, dd2;

	for (n = 0; n < cell->cnt; n++)
	{
		nid = cell->pid[n];
		tmp = (ppts + ddim * nid);
		if (tmp == pnt)	// 指针地址相同
			continue;
		dd2 = EuLen2D(tmp, pnt);
		if (dd2 < dnp)
		{	pid = nid;	cid = cin;	dnp = dd2;	}
	}
}

// 异单元搜索
inline void SearchCellNPD(IDX_CELL *cell, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_FLT64 *pnt,
						  TSS_SNT32 cin, TSS_SNT32 &pid, TSS_SNT32 &cid, TSS_FLT64 &dnp)
{
	TSS_SNT32 n, nid;
	TSS_FLT64 *tmp, dd2;

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

		dd2 = EuLen2D(tmp, pnt);
		if (dd2 < dnp)
		{	pid = nid;	cid = cin;	dnp = dd2;	}
	}
}

// 批量异单元搜索
inline void SearchCellNPD(IDX_CELL *cidx, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_FLT64 *pnt,
						  TSS_SNT32 &pid, TSS_SNT32 &cid, TSS_FLT64 &dnp, TSS_SNT32 col,
						  TSS_SNT32 *c, TSS_SNT32 *r)
{
	IDX_CELL *cell;
	TSS_SNT32 i, k, cin;

	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;
			SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
		}
	}
	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;
			SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
		}
	}
	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;
			SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
		}
	}
	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;
			SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
		}
	}
	// 处理4个角单元
	if (c[3] != c[5] || r[3] != r[5])
	{
		cin = r[1] * col + c[1];	cell = cidx + cin;
		SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
	}
	if (c[3] != c[5] || r[4] != r[6])
	{
		cin = r[2] * col + c[1];	cell = cidx + cin;
		SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
	}
	if (c[4] != c[6] || r[3] != r[5])
	{
		cin = r[1] * col + c[2];	cell = cidx + cin;
		SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
	}
	if (c[4] != c[6] || r[4] != r[6])
	{
		cin = r[2] * col + c[2];	cell = cidx + cin;
		SearchCellNPD(cell, ppts, ddim, pnt, cin, pid, cid, dnp);
	}
}

// 构建初始三角形
long MakeFirstSimplex(IDX_GRID *uidx, TSS_FLT64 *ppts, TSS_SNT32 ddim, TSS_SNT32 *stop,
					  TSS_FLT64 dm2, TSS_SNT32 *pid, TSS_SNT32 *cid)
{
	TSS_FLT64 *bnd = uidx->bnd;
	TSS_FLT64 *spn = uidx->spn;
	TSS_SNT32 *num = uidx->num;
	IDX_CELL *cell, *cidx = uidx->idx;
	TSS_SNT32 col = num[0], row = num[1];

	TSS_FLT64 dnp;	// 欧拉距离平方
	TSS_SNT32 cin, pos[2];
	TSS_SNT32 c[8], r[8], i, k;

	pid[0] = pid[1] = pid[2] = -1;
	// 1.搜索点pid[0]:遍历单元
	{
		for (i = 0; i < row; i++)
		{
			cin = i * col;
			for (k = 0; k < col; k++)
			{
				cell = cidx + cin + k;
				if (cell->cnt > 0)
				{
					pid[0] = cell->pid[0];	cid[0] = cin + k;
					break;
				}
			}
			if (pid[0] != -1)
				break;
		}
		if (pid[0] == -1)	// 不能找到第1点
			return 0;
	}
	// 2.搜索到点pid[0]最近的点pid[1]
	{
		TSS_FLT64 *pnt = ppts + pid[0] * ddim;
		// 计算点pid[0]所处单元
		CalcGridCell(uidx, pnt, &pos[0]);

		dnp = dm2;
		k = pos[0], i = pos[1];
		cin = i * col + k;	cell = cidx + cin;
		SearchCellNPS(cell, ppts, ddim, pnt, cin, pid[1], cid[1], dnp);

		c[1] = c[2] = k;	r[1] = r[2] = i;
		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;

		if (pid[1] == -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];
				SearchCellNPD(cidx, ppts, ddim, pnt, pid[1], cid[1], dnp, col, c, r);

				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])
					break;
				if (pid[1] != -1)	// 已经找到候选点
					break;
			}
		}

		if (pid[1] == -1)	// 不能找到第2点
			return 0;
		// 计算最后搜索边界
		TSS_FLT64 rad = sqrt(EuLen2D(pnt, ppts + ddim * pid[1]));
		CalcSearchBound(bnd, spn, col, row, pnt, rad, c[5], c[6], r[5], r[6]);
		// 根据搜索边界判断是否需要继续搜索
		while (1)
		{
			c[1]--;	if (c[1] < c[5])	c[1] = c[3] = c[5];	// 说明前一列已经到达搜索边界
			c[2]++;	if (c[2] > c[6])	c[2] = c[4]	= c[6];
			r[1]--;	if (r[1] < r[5])	r[1] = r[3] = r[5];	// 说明前一行已经到达搜索边界
			r[2]++;	if (r[2] > r[6])	r[2] = r[4] = r[6];
			SearchCellNPD(cidx, ppts, ddim, pnt, pid[1], cid[1], dnp, col, c, r);

			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])
				break;
		}