ITK入门教程（16）法线作为像素类型

1.概述

使用几何对象表面的点和与这些点相关的法线来表示几何对象是很常见的。使用itk::PointSet 类可以方便地对这种结构进行实例化。

用来表示表面法线和函数梯度的常见类是itk::CovariantVector。一个协变向量在仿射变换中的行为方式与向量不同，特别是在各向异性的缩放比例上。如果一个协变向量表示一个函数的梯度，那么变换后的协变向量将也是变换后函数的有效梯度，这种特性是常规向量所不具备的。

2.过程详解

下面的代码展示了如何用一个协变向量CovariantVector作为点集类PointSet的像素类型PixelType。这个例子演示了一个可形变模型如何在位势函数（potential function）的梯度的影响下移动。

为了使用协变向量类，必须包含它的头文件和点集PointSet 的头文件。

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

协变向量类是在用来表示空间坐标和空间维数的类型上进行模板化的。由于像素类型是独立于点类型的，所以我们可以自由地选择任意维数的协变向量来作为像素类型。但是，在这里我们要阐述一下可形变模型的原则，要求表示梯度的向量的维数和空间中的点的维数相同。

const unsigned int Dimension = 3;
typedef itk::CovariantVector< float, Dimension > PixelType;

然后我们使用像素类型(事实上是协变向量)来对点集类型进行实例化，进而创建一个点集对象。

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

下面的代码生成了一个圆并将梯度值分配给点。这个例子计算得到的协变向量元素用来表示圆上的法线。

PointSetType::PixelType gradient;
PointSetType::PointType point;
unsigned int pointId = 0;
const 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; // flat on the Z plane
gradient[0] = std::sin(angle);
gradient[1] = std::cos(angle);
gradient[2] = 0.0; // flat on the Z plane
pointSet->SetPoint( pointId, point );
pointSet->SetPointData( pointId, gradient );
pointId++;
}

3.代码展示

//  我们想在这里说明可变形模型的精神。然后要求表示梯度的向量与空间中的点具有相同的维数。
//  用来表示表面法线和函数梯度的常见类是itk::CovariantVector。一个协变向量在仿射变换中的行为方式与向量不同，特别是在各向异性的缩放比例上。
//  如果一个协变向量表示一个函数的梯度，那么变换后的协变向量将也是变换后函数的有效梯度，这种特性是常规向量所不具备的。
//
//  面的代码展示了如何用一个协变向量CovariantVector作为点集类PointSet的像素类型PixelType。
//  这个例子演示了一个可形变模型如何在位势函数（potential function）的梯度的影响下移动。

int
main(int, char *[])
{
  //  CovariantVector类是在用于表示空间坐标的类型和空间维度上模板化的。由于PixelType独立于PointType，我们可以自由选择协变向量的任何维度作为像素类型。然而,
  //  我们想在这里说明可变形模型的精神。然后要求表示梯度的向量与空间中的点具有相同的维数。

  constexpr unsigned int Dimension = 3;
  using PixelType = itk::CovariantVector<float, Dimension>;

  //  然后我们使用PixelType(实际上是CovariantVectors)实例化PointSet类型，并随后创建一个PointSet对象。
  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; // flat on the Z plane
    gradient[0] = std::sin(angle);
    gradient[1] = std::cos(angle);
    gradient[2] = 0.0; // flat on the Z plane
    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;
    ++pixelIterator;
    ++pointIterator;
  }
  //
  //  CovariantVector类不会重载Point的+操作符。也就是说，协变向量不能加到点上以得到新的点。
  //  此外，由于我们在示例中忽略了物理，我们还被迫手动在梯度分量和点的坐标之间做非法的加法。
  //  注意，ITK几何类中没有一些基本运算符完全是有意为之，目的是防止错误地使用它们所代表的数学概念。


  //  我们终于可以访问所有的点并打印出新的值。
  //
  pointIterator = pointSet->GetPoints()->Begin();
  pointEnd = pointSet->GetPoints()->End();
  while (pointIterator != pointEnd)
  {
    std::cout << pointIterator.Value() << std::endl;
    ++pointIterator;
  }


  return EXIT_SUCCESS;
}

4.结果展示