ITK 实例 （三）

ITK 官方网址：https://itk.org/ITK/project/project.html

（三） ITK 中的点集

（1） 创建点集

下面实例展示了如何创建点集，如何将生成的指定点插入到点集中，以及如何从点集中读取指定点。

#include "itkPointSet.h"

int main()
{
	using PointSetType = itk::PointSet<unsigned short, 3>;
	PointSetType::Pointer pointSet = PointSetType::New();

	using PointType = PointSetType::PointType;

	PointType p0;
	p0[0] = -1.0;  // x coordinate
	p0[1] = -1.0;  // y coordinate
	p0[2] = 0.0;  // z coordinate

	PointType p1;
	p1[0] = 1.0; 
	p1[1] = -1.0;
	p1[2] = 0.0;

	PointType p2;
	p2[0] = 1.0;
	p2[1] = 1.0;
	p2[2] = 0.0; 

	pointSet->SetPoint( 0, p0 );
	pointSet->SetPoint( 1, p1 );
	pointSet->SetPoint( 2, p2 );

	const unsigned int numberOfPoints = pointSet->GetNumberOfPoints();
	std::cout << "number in the pointSet is: " <<
		numberOfPoints << std::endl;

	PointType pp;
	bool pointExists = pointSet->GetPoint(0, &pp); //将点集中的 0 号点读入到 pp 中

	if (pointExists)
	{
		std::cout << "pp point is: " << pp << std::endl;
	}

	return EXIT_SUCCESS;
}

（2）访问点集中的点

由于使用 SetPoint 将点加入点集并不是一个高效的方式，因此可以通过  ITK 中的点容器首先将所有点放置到点容器中，再将这个点容易插入到点集中，将点放置到点的容器用的是 InsertElement() ，然后用 SetPoitnts 将点容器放到点集中。获取访问点集中点的权限用的是 GetPoints()。

#include "itkPointSet.h"

int main()
{
	using PointSetType = itk::PointSet<unsigned short, 3>;

	using PointsContainer = PointSetType::PointsContainer;
	PointsContainer::Pointer points = PointsContainer::New();

	using PointType = PointSetType::PointType;
	PointType p0;
	PointType p1;
	PointType p2;

	p0[0] = -1.0;
	p0[1] = 0.0;
	p0[2] = 0.0;  // point 0 ={-1,0,0}

	p1[0] = 1.0;
	p1[1] = 0.0;
	p1[2] = 0.0;  // point 1 = {1,0,0}

	p2[0] = 0.0;
	p2[1] = 1.0;
	p2[2] = 2.0;  //point 2 = {0,1,2}

	unsigned int pointId = 0;
	points->InsertElement(pointId++, p0);
	points->InsertElement(pointId++, p1);
	points->InsertElement(pointId++, p2);

	PointSetType::Pointer pointSet = PointSetType::New();
	pointSet->SetPoints(points);

	PointsContainer::Pointer points2 = pointSet->GetPoints();  // 获取点集中的点的容器
	using PointsIterator = PointsContainer::Iterator;

	PointsIterator pointIterator = points->Begin();  //迭代器指向容器的首位置
	PointsIterator end = points->End();

	while (pointIterator != end)
	{
		PointType p = pointIterator.Value();  //迭代器指向的当前数值
		std::cout << p << std::endl;  //打印出这个点
		++pointIterator;  // 迭代器移动到下一个位置
	}

	std::cout << pointSet->GetNumberOfPoints() << std::endl;  //输出点的数量
	std::cout << points2->Size() << std::endl;

	return EXIT_SUCCESS;
}

（3）访问点集中点的数据

用 GetPoints( ) 可以访问点集中的点，但并不是访问所有点的高效的方式，若想遍历点集中的所有点，最有效的方法是使用点迭代器。

#include "itkPointSet.h"

int main()
{
	using PixelType = unsigned short;
	using PointSetType = itk::PointSet<PixelType, 3>;

	PointSetType::Pointer pointSet = PointSetType::New();
	unsigned int dataId = 0;
	PixelType value = 79;
	pointSet->SetPointData(dataId++, value);

	const bool found = pointSet->GetPointData(dataId, &value);
	if (found)
	{
		std::cout << "Pixel value = " << value << std::endl;
	}

	using PointDataContainer = PointSetType::PointDataContainer;
	PointDataContainer::Pointer pointData = PointDataContainer::New();

	unsigned int pointId = 0;

	PixelType value0 = 34;
	PixelType value1 = 67;

	pointData->InsertElement(pointId++, value0);
	pointData->InsertElement(pointId++, value1);  //先将点的数据插入到点数据容器中

	pointSet->SetPointData(pointData);
	PointDataContainer::Pointer pointData2 = pointSet->GetPointData();
	
	using PointDataIterator = PointDataContainer::Iterator;
	PointDataIterator pointDataIterator = pointData2->Begin();
	PointDataIterator end = pointData2->End();
	while (pointDataIterator != end)
	{
		PixelType p = pointDataIterator.Value();   // 读取迭代器指向的点的数据
		std::cout << p << std::endl;  //打印数据
		++pointDataIterator;  // 移动迭代器到下一个位置
	}

	return EXIT_SUCCESS;
}

（4）RGB 作为像素类型 

下面的实例中，用点的位置的两倍数值来表示 RGB 各分量的值。

#include "itkRGBPixel.h"
#include "itkPointSet.h"

int main()
{
	using PixelType = itk::RGBPixel<float>;
	using PointSetType = itk::PointSet<PixelType, 3>;
	PointSetType::Pointer pointSet = PointSetType::New();

	PointSetType::PixelType pixel;
	PointSetType::PointType point;

	unsigned int pointId = 0;
	constexpr double radius = 3.0;

	for (unsigned int i = 0; i < 360; i++)
	{
		const double angle = i*itk::Math::pi / 180.0;
		point[0] = radius * std::sin(angle);
		point[1] = radius * std::cos(angle);
		point[2] = 1.0;
		pixel.SetRed(point[0] * 2.0);
		pixel.SetGreen(point[1] * 2.0);
		pixel.SetBlue(point[2] * 2.0);
		pointSet->SetPoint(pointId, point);
		pointSet->SetPointData(pointId, pixel);
		pointId++;
	}

	using PointIteratorStyle = PointSetType::PointsContainer::ConstIterator;  
    //因为不希望修改容器迭代器指向的位置的数值，所以用 constIterator
	PointIteratorStyle pointIterator = pointSet->GetPoints()->Begin();
	PointIteratorStyle pointEnd = pointSet->GetPoints()->End();

	while (pointIterator != pointEnd)
	{
		point = pointIterator.Value();
		std::cout << point << std::endl;
		++pointIterator;
	}

	using PointDataIteratorStyle = PointSetType::PointDataContainer::ConstIterator;
	PointDataIteratorStyle pixelIterator = pointSet->GetPointData()->Begin();
	PointDataIteratorStyle pixelEnd = pointSet->GetPointData()->End();

	while (pixelIterator != pixelEnd)
	{
		pixel = pixelIterator.Value();
		std::cout << pixel << std::endl;
		++pixelIterator;
	}
	
	return EXIT_SUCCESS;
}

（5）矢量作为像素类型

注意：ITK 中 PixelType 独立于 PointType。下面的代码生成一个球，将向量值赋给每个点，计算向量的分量，以表示圆的切线。

#include "itkVector.h"
#include "itkPointSet.h"

int main()
{
	constexpr unsigned int Dimension = 3;
	using PixelType = itk::Vector<float, Dimension>;

	using PointSetType = itk::PointSet<PixelType, Dimension>;
	PointSetType::Pointer pointSet = PointSetType::New();

	PointSetType::PixelType tangent;  //切线
	PointSetType::PointType point;

	unsigned int pointId = 0;
	constexpr double radius = 300.0;

	for (unsigned int i = 0; i < 360; i++)
	{
		const double angle = i* itk::Math::pi / 180.0;
		point[0] = radius * std::sin(angle);
		point[1] = radius * std::cos(angle);
		point[2] = 1.0;  
		tangent[0] = std::cos(angle);
		tangent[1] = -std::sin(angle);
		tangent[2] = 0.0;
		pointSet->SetPoint(pointId, point);
		pointSet->SetPointData(pointId, tangent);
		pointId++;
	}

	using PointDataIterator = PointSetType::PointDataContainer::ConstIterator;
	PointDataIterator pixelIterator = pointSet->GetPointData()->Begin();
	PointDataIterator pixelEnd = pointSet->GetPointData()->End();

	using PointIterator = PointSetType::PointsContainer::Iterator;
	PointIterator pointIterator = pointSet->GetPoints()->Begin();
	PointIterator pointEnd = pointSet->GetPoints()->End();

	while (pixelIterator != pixelEnd && pointIterator != pointEnd)
	{
		pointIterator.Value() = pointIterator.Value() + pixelIterator.Value();
		++pixelIterator;
		++pointIterator;
	}

	pointIterator = pointSet->GetPoints()->Begin();
	pointEnd = pointSet->GetPoints()->End();

	while (pointIterator != pointEnd)
	{
		std::cout << pointIterator.Value() << std::endl;
		++pointIterator;
	}
	return EXIT_SUCCESS;

（6）法线作为像素类型

使用几何对象表面上的点和与这些点相关联的法线来表示几何对象是很常见的。表示曲面法线和函数梯度的类是：itk::CovariantVector，它与 itk::Vector 的区别在于它在仿射变换下的行为，如果协变向量（covariant vector）表示函数的梯度，则变换后的协变向量仍然是变换函数的有效梯度，这与普通向量（vector）的性质不同。

#include "itkCovariantVector.h"
#include "itkPointSet.h"

int main()
{
	constexpr unsigned int Dimension = 3;
	using PixelType = itk::CovariantVector<float, Dimension>;
	
	using PointSetType = itk::PointSet<PixelType, Dimension>;
	PointSetType::Pointer pointSet = PointSetType::New();

	PointSetType::PixelType gradient;
	PointSetType::PointType point;

	unsigned int pointId = 0;
	constexpr double radius = 300.0;

	for (unsigned int i = 0; i < 360; i++)
	{
		const double angle = i *std::atan(1.0) / 45.0;
		point[0] = radius * std::sin(angle);
		point[1] = radius * std::cos(angle);
		point[2] = 1.0;
		gradient[0] = std::sin(angle);
		gradient[1] = std::cos(angle);
		gradient[2] = 0.0;
		pointSet->SetPoint(pointId, point);
		pointSet->SetPointData(pointId, gradient);
		pointId++;
	}

	using PointDataIterator = PointSetType::PointDataContainer::ConstIterator;
	PointDataIterator pixelIterator = pointSet->GetPointData()->Begin();
	PointDataIterator pixelEnd = pointSet->GetPointData()->End();

	using PointIterator = PointSetType::PointsContainer::Iterator;
	PointIterator pointIterator = pointSet->GetPoints()->Begin();
	PointIterator pointEnd = pointSet->GetPoints()->End();

	while (pixelIterator != pixelEnd && pointIterator != pointEnd)
	{
		point = pointIterator.Value();
		gradient = pixelIterator.Value();
		for (unsigned int i = 0; i < Dimension; i++)
		{
			point[i] += gradient[i];
		}
		pointIterator.Value() = point;
		//std::cout << pointIterator.Value() << std::endl;
		++pixelIterator;
		++pointIterator;
	}

	pointIterator = pointSet->GetPoints()->Begin();
	pointEnd = pointSet->GetPoints()->End();

 	while (pointIterator != pointEnd)
 	{
 		std::cout << pointIterator.Value() << std::endl;
 		++pointIterator;
	}

	return EXIT_SUCCESS;
}