V-S 分水岭变换算法

// V-S版本的分水岭算法
bool Watershed::WatershedTransformOfVS(IplImage *src, int &basinCount)
{
#define MASK -2 /* initial value of a threshold level 标记为当前正在处理的梯度层 */
#define INIT -1 /* initial value of imgOut */
#define WSHED 0 /* value of the pixels belonging to the watersheds */
#define INIT_DIST 0 /* initial value of imgDist */
    /*
  . --input:  imgIn, decimal image;
	--output: image of the labeled watersheds;
  . Initializations:
	--Value INIT is assigned to each pixel of imgOut:
	--current_label = 0
	--current_dist: integer variable
	--imgDist:work image (of distances),initialized to 0;
    */

	if (src == NULL)
	{
		return false;
	}
	
	const int width		= src->width;
	const int height	= src->height;
	const int widthStep = src->widthStep;
	const int channels  = src->nChannels;
	uchar *imgIn = (uchar *)src->imageData;
	
	if (channels != 1) // src必须是单通道图像
	{ 
		return false;
	}

	// V-S 分水岭变换算法
	//
	// 初始化过程
	//
	BASIN_NUM *imgOut = new BASIN_NUM[height * widthStep];          // 上面标记有水坝的输出图像
	for (int i = 0; i < height * widthStep; i ++) imgOut[i] = INIT; // 初始化输出图
	
	int *imgDist = new int[height * widthStep];       // 上面标记关联距离的图像
	for (int i = 0; i < height * widthStep; i ++) 
		imgDist[i] = INIT_DIST;   // 初始化距离图
	
	int currentNum  = 0; // 盆地的编号
	int currentDist = 0; // 

	int Direct[][2] = {{0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}}; // 八个方位

	//
	// 梯度排序过程
	//
	vector<POINT_GRADIENT> vPointsGradient; // 带梯度值的所有点的有序集合
	for (int col = 0; col < height; col ++)
	{
		for (int row = 0; row < width; row ++)
		{
			if (255 == (int)imgIn[col * widthStep + row])
			{
				continue;
			}
			POINT_GRADIENT point;
			point.x = row;
			point.y = col;
			point.gradient = (int)imgIn[col * widthStep + row];
			vPointsGradient.push_back(point);
		}
	}
	sort(vPointsGradient.begin(), vPointsGradient.end(), cmp);

	//
	// 淹没过程
	//
	unsigned int nCurGradient;	                 // 当前处理的梯度级 
	vector<POINT_GRADIENT> vPointsCurGrad;       // 当前梯度级的所有点的集合

	int iPoint = vPointsGradient.size() - 1;         // 遍历点的下标
	nCurGradient = vPointsGradient[iPoint].gradient; // 最小梯度值

	for ( ; iPoint >= 0; iPoint --) // 从小到大遍历所有梯度点
	{
		if (nCurGradient == vPointsGradient[iPoint].gradient) // 得到每一层的梯度点集合
		{
			vPointsCurGrad.push_back(vPointsGradient[iPoint]);
			continue;
		}
		
		// >>>处理梯度级为nGradient 的该层上的像素点
		
		//
		// 扩展每个盆地
		//
		POINT_GRADIENT p ; // 当前点
		POINT_GRADIENT p_; // 当前点p 的领域点
		queue<POINT_GRADIENT> fifoQue;
		for (int iPoint = 0; iPoint < vPointsCurGrad.size(); iPoint ++) // 遍历当前梯度级的每个点
		{
			p = POINT_GRADIENT(vPointsCurGrad[iPoint].x, vPointsCurGrad[iPoint].y, 
				vPointsCurGrad[iPoint].gradient); // 当前点
			imgOut[p.x + p.y * widthStep] = MASK;

			bool flag = false; // 如果p 的领域点中存在盆地或者分水岭则flag为true
			for (int iDirect = 0; iDirect < 8; iDirect ++) 
			{
			    // p_ 当前点p的领域点
				p_.x = p.x + Direct[iDirect][0];
				p_.y = p.y + Direct[iDirect][1];
				if (p_.x < 0 || p_.x >= width || p_.y < 0 || p_.y >= height)
				{
					continue;
				}
				p_.gradient = (int)imgIn[p_.x + p_.y * widthStep];

				if (imgOut[p_.x + p_.y * widthStep] > 0 || imgOut[p_.x + p_.y * widthStep] == WSHED)
				{
					flag = true;
					break;
				}
			}
			if (flag)
			{ /* 此时p 是距离盆地或分水岭最近一排的点 */
				imgDist[p.x + p.y * widthStep] = 1;
				fifoQue.push(p);                     
			}
		}  
		/* 此时fifoQue里存放的是和每个盆地的边缘距离为1 的梯度值为n+1层的点，其中n为盆地边缘的梯度值 */

		POINT_GRADIENT fictPoint = POINT_GRADIENT(0, 0, 256); // fictitious_pixel 虚拟像素点，队列遍历一次结束的标记点
		currentDist = 1; 
		fifoQue.push(fictPoint);

		while (1) // repeat indefinitely
		{
			p = fifoQue.front();
			fifoQue.pop();
			if (p.gradient == fictPoint.gradient) // if p = fictPoint
			{ /* bfs 完了一圈 */
				if (fifoQue.empty()) 
					break;
				else 
				{
					fifoQue.push(fictPoint);
					currentDist += 1;        // bfs 下一圈的点，这些点的距离变大一个单位
					p = fifoQue.front();
					fifoQue.pop();
				}
			}

			//>>> 

			for (int iDirect = 0; iDirect < 8; iDirect ++)
			{
				// p_ 当前点p的领域点
				p_.x = p.x + Direct[iDirect][0];
				p_.y = p.y + Direct[iDirect][1];
				if (p_.x < 0 || p_.x >= width || p_.y < 0 || p_.y >= height)
				{ 
					continue; 
				}
				p_.gradient = (int)imgIn[p_.x + p_.y * widthStep];
			
				if (imgDist[p_.x + p_.y * widthStep] < currentDist && 
					(imgOut[p_.x + p_.y * widthStep] > 0 || imgOut[p_.x + p_.y * widthStep] == WSHED))
				{ /* 如果p_ 是一个已经被标记为盆地或者是分水岭上的点 */
					if (imgOut[p_.x + p_.y * widthStep] > 0) // 如果p_ 已经是盆地
					{
						if (imgOut[p.x + p.y * widthStep] == MASK || 
							imgOut[p.x + p.y * widthStep] == WSHED)
						{ /* p为未标记点或是分水岭 */
							imgOut[p.x + p.y * widthStep] = imgOut[p_.x + p_.y * widthStep];
						}
						else if(imgOut[p.x + p.y * widthStep] != imgOut[p_.x + p_.y * widthStep])
						{ /* p 和p_ 分别属于不同的盆地 */
							imgOut[p.x + p.y * widthStep] = WSHED; // 建水坝
						}
					}
					else // 如果p_ 是分水岭上的点
					{
						if (imgOut[p.x + p.y * widthStep] == MASK)
						{ /* 如果p 未标记 */
							imgOut[p.x + p.y * widthStep] = WSHED; // 自己盆地的分水岭
						}
					}
				}
				else if (imgOut[p_.x + p_.y * widthStep] == MASK && 
					imgDist[p_.x + p_.y * widthStep] == 0)
				{
					imgDist[p_.x + p_.y * widthStep] = currentDist + 1;
					fifoQue.push(p_);
				}
			}

			//<<<
		}

		//
		// 寻找新的极小区
		//
		for (int iPoint = 0; iPoint < vPointsCurGrad.size(); iPoint ++) // 遍历当前梯度的每个点
		{
			p = POINT_GRADIENT(vPointsCurGrad[iPoint].x, vPointsCurGrad[iPoint].y, 
				vPointsCurGrad[iPoint].gradient); // 当前点
			imgDist[p.x + p.y * widthStep] = 0; /* p 点的关联距离重置为 0 */

			if (imgOut[p.x + p.y * widthStep] == MASK)
			{
				currentNum += 1; // 新增一个盆地
				fifoQue.push(p);
				imgOut[p.x + p.y * widthStep] = currentNum;
				
				// 扩展极小区
				POINT_GRADIENT pt; // 临时点
				while (! fifoQue.empty())
				{
					pt = fifoQue.front();
					fifoQue.pop();

					POINT_GRADIENT pt_; // 临时点pt 的领域点
					for (int iDirect = 0; iDirect < 8; iDirect ++)
					{
						pt_.x = pt.x + Direct[iDirect][0];
						pt_.y = pt.y + Direct[iDirect][1];
						if (pt_.x < 0 || pt_.x >= width || pt_.y < 0 || pt_.y >= height)
						{
							continue;
						}
						pt_.gradient = (int)imgIn[pt_.x + pt_.y * widthStep];
						if (pt_.gradient != pt.gradient)
						{
							continue;
						}
						
						if (imgOut[pt_.x + pt_.y * widthStep] == MASK)
						{
							fifoQue.push(pt_);
							imgOut[pt_.x + pt_.y * widthStep] = currentNum;
						}
					}
				}
			}
		}

		// <<<处理完梯度级为nGradient 的该层上的像素点

		nCurGradient = vPointsGradient[iPoint].gradient;
		vPointsCurGrad.clear();
		vPointsCurGrad.push_back(vPointsGradient[iPoint]);
	}
	basinCount = currentNum;
	for (int col = 0; col < height; col ++)
	{
		for (int row = 0; row < width; row ++)
		{
			if (imgOut[col * widthStep + row] == 0) imgOut[col * widthStep + row] = 255;
			imgIn[col * widthStep + row] = (uchar)imgOut[col * widthStep + row];
		}
	}

	return true;
}