【ITK库学习】使用itk库进行图像分割（四）：水平集分割

1、水平集

水平集是跟踪轮廓和表面运动的一种数字化方法。基于图像的亮度均值、梯度、边缘特征的微分计算，进行水平集分割。在itk中，所有基于分割滤波器水平集必须用浮点精度来进行操作，输出的数字化类型默认为浮点型，也可以被转换为双精度，但是不能使用整型或者无符号数据类型。

2、itkFastMarchingImageFilter 快速步进分割

该类使用快速行进求解Eikonal方程，其中速度始终为非负且仅取决于位置。

使用熵满足方案执行更新，其中仅使用“逆风”邻域。 Fast Marching 的此实现使用 std::priority_queue 来定位下一个要更新的正确网格位置。快速行进以 (N log N) 步扫描 N 个网格点，以获得锋面传播通过网格时的到达时间值。

当图像梯度很高时前面的传播速率很慢而梯度很低的区域传播速率快。此方法可是轮廓不断传播直到图像中的解剖结构边缘在那写边缘前将速率降低下来。

此类基于水平集图像类型和速度图像类型进行模板化，初始前沿由两个容器指定：一个包含已知点，另一个包含试验点。 活动点是那些已经是对象一部分的点，并且考虑包含试验点。 为了使过滤器不断发展，至少必须指定一些试验点。 例如，可以将它们指定为活动点周围的像素层。

速度函数可以指定为速度图像或速度常数，速度图像使用 SetInput() 方法设置。 如果速度图像为nullptr，则使用恒速函数，并使用SetSpeedConstant()方法指定。

如果速度函数恒定且值为1，则快速行进会产生距初始活动点的近似距离函数。FastMarchingImageFilter 在 ReinitializeLevelSetImageFilter 对象中使用，以创建距零水平集的有符号距离函数。

通过设置适当的停止值可以提前终止算法，当处理时间大于停止值时，算法终止。

FastMarchingImageFilter输出的是一个时间交叉图，对每一个像素，表达了到达这个像素位置之前所花费的时间。有两种方法可以指定输出图像信息（LargestPossibleRegion、Spacing、Origin）：(a) 直接从输入速度图像复制 (b) 由用户指定，如果用户未指定所有信息，则使用默认值。

输出信息计算如下。 如果速度图像为 nullptr 或 OverrideOutputInformation 设置为 true，则根据用户指定的参数设置输出信息。 这些参数可以使用 SetOutputRegion()、SetOutputSpacing()、SetOutputDirection() 和 SetOutputOrigin() 方法指定。 否则，如果速度图像不为nullptr，则从输入速度图像复制输出信息。

可能的改进：在当前的实现中，std::priority_queue只允许从前面取出节点并从后面放入节点，要更新堆上已有的值，需要将新节点添加到堆中，失效的旧节点留在堆上，当它从顶部移除时，它将被识别为无效并且不被使用，未来的实现可以以不同的方式实现堆，从而允许更新值，这通常需要一些上移和下移函数以及从图像到堆的反向指针图像，以便找到要更新的节点。

常用的成员函数：

• Set/GetStoppingValue()：设置/获取快速行进算法停止值，当最小尝试点的值大于停止值时，算法终止
• Set/GetSpeedConstant()：设置/获取速度常数，如果速度图像为nullptr，则SpeedConstant值将用于整个级别集，默认情况下，SpeedConstant`设置为 1.0
• Set/GetTrialPoints()：设置/获取代表初始前沿的试验点容器，试验点表示为 LevelSetNodes的VectorContainer`
• SetOutsidePoints()：设置不打算评估的点的容器
• Set/GetOutputRegion/Size/Spacing/Origin/Direction()：设置/获取输出的最大可能值、大小、间距、原点、方向计算，如果速度图像为nullptr或OverrideOutputInformation为true，则根据用户指定的参数设置输出信息，如果速度图像不为nullptr`，则从输入速度图像复制输出信息
• Set/GetOverrideOutputInformation()：同上
• Set/GetCollectPoints()：设置/获取收集点标志，检测算法以收集其访问过的所有节点的容器，对于为支持窄带的水平集算法创建窄带很有用
• Set/GetAlivePoints()：设置/获取代表初始前沿的活动点容器，活动点表示为LevelSetNodes的VectorContainer`
• Set/GetNormalizationFactor()：设置/获取速度图像的标准化因子，速度图像中的值除以该因子，这允许使用具有整数像素类型的图像来表示速度
• GetProcessedPoints()：获取已处理点的容器，如果设置了 CollectPoints 标志，算法将收集所有已处理节点的容器，这对于定义支持窄带的水平集算法创建窄带非常有用
• GetLabelImage()：获取点型标签图像
• GetLargeValue()：获得巨大价值，该值用于表示分配给尚未访问的像素的时间无穷大的概念，该值默认设置为用于表示时间交叉图的像素类型max()的一半

示例代码：

#include "itkImage.h"
#include "itkFastMarchingImageFilter.h"
typedef itk::Image<float, 3> FloatImageType;

bool fastMarchingImageFilter(FloatImageType* image, FloatImageType* outImage)
{
    typedef itk::FastMarchingImageFilter<FloatImageType, FloatImageType> FastMarchingFilterType;
	typename FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
	using NodeContainer = FastMarchingFilterType::NodeContainer;
	using NodeType = FastMarchingFilterType::NodeType;

	auto seeds = NodeContainer::New();

	FloatImageType::IndexType seedPosition;
	seedPosition[0] = 256;
	seedPosition[1] = 256;
	seedPosition[2] = 100;

	constexpr double initialDistance = 1.0;

	const double seedValue = -initialDistance;

	NodeType node;                   //创建节点作为堆栈变量，并初始化
	node.SetValue(seedValue);
	node.SetIndex(seedPosition);     
	seeds->Initialize();             //初始化
	seeds->InsertElement(0, node);   //插入每个节点

    FloatImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
    //double stopTime = 100;
	fastMarching->SetInput(image);
	fastMarching->SetTrialPoints(seeds);
	fastMarching->SetSpeedConstant(1.0);
	fastMarching->SetOutputSize(size);
	//fastMarching->SetStoppingValue(stopTime);
	try
	{
		fastMarching->Update();
	}
	catch (itk::ExceptionObject& ex)
	{
		//读取过程发生错误
		std::cerr << "Error: " << ex << std::endl;
		return false;
	}

	outImage = fastMarching->GetOutput();
	return true;
}

3、itkShapeDetectionLevelSetImageFilter 快速步进分割

该类根据用户提供的边缘势图分割图像中的结构。

此类是水平集方法分段过滤器，初始轮廓向外（或向内）传播，直到它“粘住”形状边界，通过使用基于用户提供的边缘势图的水平集速度函数来完成的。

该过滤器需要两个输入：
第一个输入是初始水平集，初始水平集是包含初始轮廓/表面作为零水平集的真实图像，例如，通常使用距初始轮廓/表面的带符号距离函数；请注意，对于此算法，初始轮廓必须完全位于要分割的结构内部（或完全外部）；
第二个输入是特征图像，即边缘势图，它的一般特征是：在边缘附近的区域中具有接近于零的值，并且在形状本身内部具有接近于1的值， 通常边缘势图是根据图像梯度计算的，例如：

其中: I 是图像强度，(∇*G) 是高斯算子的导数。

参数
PropagationScaling 参数可用于在向外传播（正缩放参数）与向内传播（负缩放参数）之间切换。
可以使用PropagationScaling和CurvatureScaling参数的组合来调整生成的轮廓/表面的平滑度。CurvatureScaling参数越大，生成的轮廓越平滑，为了使该算法正确运行，CurvatureScaling 参数应为非负数。 为了遵循 Malladi 等人论文中的实现，将 PropagationScaling 设置为 ±1.0，将 CurvatureScaling设置为 ϵ。

请注意，此过滤器没有平流项，设置平流缩放不会产生任何效果。

输出
滤波器输出单个标量实值图像，输出图像中的负值表示分割区域的内部，图像中的正值表示分割区域的外部，图像的零交叉点对应于传播前沿的位置。
有关详细信息，请参阅 SparseFieldLevelSetImageFilter 和 SegmentationLevelSetImageFilter。

常用的成员函数：

• Set/GetInput()：设置初始水平集
• Set/GetFeatureImage()：设置/获取用于水平集方程速度函数的特征图像
• Set/GetPropagationScaling()：设置/获取传播速度的缩放比例
• Set/GetCurvatureScaling()：设置/获取曲率的缩放，使用此参数可增加曲率对曲面运动的影响，相对于平流和传播值而言，较高的值将提供更平滑的表面
• Set/GetMaximumRMSError()：设置/获取水平集函数的最大均方根
• Set/GetNumberOfIterations()：设置/获取迭代次数

示例代码：

#include "itkImage.h"
#include "itkShapeDetectionLevelSetImageFilter.h"

typedef itk::Image<float, 3> FloatImageType;

bool shapeDetectionLevelSetImageFilter(FloatImageType* image, FloatImageType* featureImage, FloatImageType* outImage)
{
    const double propagationScaling = 1.0;         //传播缩放
	const double curvatureScaling = 0.5;           //曲率缩放

    typedef itk::ShapeDetectionLevelSetImageFilter<FloatImageType, FloatImageType> ShapeDetectionLevelSetFilterType;
	typename ShapeDetectionLevelSetFilterType::Pointer shapeDetection = ShapeDetectionLevelSetFilterType::New();
	shapeDetection->SetInput(image);                         //初始水平集
	shapeDetection->SetFeatureImage(featureImage);           //特征图像，一个边缘潜在图像
	shapeDetection->SetMaximumRMSError(0.02);
	shapeDetection->SetNumberOfIterations(800);
	shapeDetection->SetPropagationScaling(propagationScaling);
	shapeDetection->SetCurvatureScaling(curvatureScaling);
	try
	{
		shapeDetection->Update();
	}
	catch (itk::ExceptionObject& ex)
	{
		//读取过程发生错误
		std::cerr << "Error: " << ex << std::endl;
		return false;
	}

	outImage = shapeDetection->GetOutput();
	return true;
}