«
»

一种简单的GPU三维图像分割算法

1 Star2 Stars3 Stars4 Stars5 Stars
Loading ... Loading ...

图像分割的定义是将整幅图像区域分割成各自互不交叉的小区域。图像分割的算法有很多，教科书上讲到一些经典的分割算法，包括阈值法、区域生长法、分裂与合并算法、边界跟踪与连接算法、分水岭算法，这些算法适合提取出灰度比较一致的各个区域；此外还有纹理分割算法，严格说来，与其称之为纹理分割还不如称之为纹理检测与分类，通过提取纹理特征，如：Laws特征、梯度共生矩阵特征、分形特征，得到每一个点或小邻域的特征向量，然后利用特征向量将每个点分到不同的纹理类中去。
有些扯远了，还是回到三维图像分割上来，由于三维图像的数据量是巨大的，暂不考虑复杂的纹理特征，我们仅仅希望有在GPU上实现一种并行的快速算法，能够自动标记出三维数据中所有灰度比较一致的各个连通区域。我们先看看教科书上的方法，阈值法不能提取多个灰度不一致的区域，区域生长不容易做到每个区域种子点自动选取，分裂与合并的逻辑稍显复杂且分割的区域边缘会出现锯齿，边界跟踪与连接在三维上根本没办法实现，分水岭又是串行递归的-_-!!!
假设已经是有了一幅二值图，我们就可以用连通区标记的算法，为二值图中的每一个小的区域标上不同的标号，这就是常用的二值图CCL（Connected Component Labeling）算法，可以扫描填充法、FloodFill等方式来具体实现这个算法。算法的原理是，判断当前点的邻域内有没有和当前点取值相同的点，如果有，就把他们标记为相同的标号。利用这个思想，如果在判断的时候，不是判断当前点的邻域内有没有和当前点取值相同的点，而是判断当前点的邻域内有没有和当前点取值相近的点——有点类似于梯度区域生长算法，就可以设计出适用于灰度图像的灰度连通区标记算法灰度图CCL算法，如果对每一个点进行如上判断，就可以预先生成一幅具有二值图，且对每个点的判断不依赖于其它点判断的结果，可以做到完全并行。这一段代码很简单，假设体数据已经做了去噪处理，然后用For循环遍历每一个点就能搞定（对于Cuda实现，就是把最里层的内容写到kernel里面），贴代码：

float fCurrentVal = 0.f;
float fThreshold = 0.1f;
int nIndex = 0;
unsigned char* puchMaskPtr = NULL;
for (i = 1; i < nDepth - 1; i ++)
  {
    for (j = 1; j < nHeight - 1; j ++)
	{
	   for (k = 1; k < nWidth - 1; k ++)
		{
		  nIndex = i * nWidth * nHeight + j * nWidth + k;
		  puchMaskPtr = puchMask + nIndex;
		  fCurrentVal = pfData[nIndex];
		   if (fCurrentVal < 0.01)
		  {
		 // 不处理被认为背景的像素
		 continue;
		  }// 1
	if (fabs(pfData[(i - 1) * nWidth * nHeight + (j - 1) * nWidth + k - 1] - fCurrentVal) < fThreshold)
	{
	   (*puchMaskPtr) ++;
	}
	// 2
	if (fabs(pfData[(i - 1) * nWidth * nHeight + (j - 1) * nWidth + k] - fCurrentVal) < fThreshold)
	{
		(*puchMaskPtr) ++;
	}
	// 3
	if (fabs(pfData[(i - 1) * nWidth * nHeight + (j - 1) * nWidth + k + 1] - fCurrentVal) < fThreshold)
	{
	 (*puchMaskPtr) ++;
	}
	// 对邻域内其他的点也做同样处理，这里略
	if (*puchMaskPtr < 8)
	{
	//只有当当前点跟周围的点比较接近的时候，才认为是目标点
	*puchMaskPtr = 0;
	}
	else
	{
	*puchMaskPtr = 255;
	       }
	    }
	  }
	}

这里考虑了每个像素的26邻域，可以按需要改造成根据18邻域（除去立方体角上的像素）或者6邻域（只考虑前后左右上下6个像素）判断；如果是Z方向的尺度于XY方向不一致，可以修改对Z方向邻域像素判断时所用的阈值。
下一步，就是三维的二值图CCL算法了，暂时没有想好怎样在GPU上实现这段代码，如果采用FloodFill这个很串行的算法，就没有放到GPU上实现的必要了。直接上三维FloodFill代码：

	stack stackIndex;
 
	int nPlanePixelNumber = nHeight * nWidth;
 
	// 只考虑6邻域
	int nIdxForeDev	=	-nPlanePixelNumber;
	int nIdxNextDev	=	nPlanePixelNumber;
	int nIdxLeftDev	=	-1;
	int nIdxRightDev=	1;
	int nIdxUpDev	=	-nWidth;
	int nIdxDownDev	=	nWidth;
 
	int nIdxCur = 0;
	Pos ptCur = posSeed;
	int nIndex = ptCur.z * nPlanePixelNumber + ptCur.y * nWidth + ptCur.x;
	unsigned short uiOldVal = puiLabel[nIndex];
	unsigned int unVolume = 0;
	puiLabel[nIndex] = uiNewVal;
	// 种子先入栈
	stackIndex.push(nIndex);
	unVolume = 1;
 
	int nMaxSize = 0;
 
	// mark 添加坐标的最大最小值
	int nIndexMin = nIndex;
	int nIndexMax = nIndex;
 
	while(stackIndex.empty() == false)
	{
		nIndex = stackIndex.top();
		stackIndex.pop();
 
		nIdxCur = nIndex + nIdxForeDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
 
		nIdxCur = nIndex + nIdxNextDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
 
		nIdxCur = nIndex + nIdxLeftDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
 
		nIdxCur	= nIndex + nIdxRightDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
 
		nIdxCur	= nIndex + nIdxUpDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
 
		nIdxCur = nIndex + nIdxDownDev;
		if (puiLabel[nIdxCur] != 0 &amp;&amp; puiLabel[nIdxCur] != uiNewVal)
		{
			puiLabel[nIdxCur] = uiNewVal;
			stackIndex.push(nIdxCur);
			unVolume ++;
		}
	}
 
	return 0;
 
}

这里考虑到为了限制堆栈的大小，只处理6邻域。

以上只是一个用于灰度图像CCL算法最简单的实现，将算法的过程分为并行和串行两部分，分别放在GPU上和CPU上运行，但算法还存在一些问题，如：没办法对灰度渐变的区域或粘连区域进行精确分割，只能分割出一整块区域。其实CCL算法作为机器视觉中对目标物进行提取的最基本算法，要想完全在GPU上实现不是这么轻松，nVidia专门为此举行了一次编程大赛来解决这个问题，参见下面的网址：

http://www.prnewswire.com/news-releases/topcoder-and-nvidia-announce-winners-of-cuda-superhero-challenge-at-nvidia-gpu-technology-conference-63121422.html

Tags: , , , | Print Print | 485 views

Leave a Reply