ITK入门教程（7）ITK定义原点和间距

主要内容

1.概述
2.代码实现
3.结果
参考目录

1.概述

尽管ITK可以用来执行通用图像处理的任务，但是这个工具包的主要目的是处理医学图像数据。因此关于图像的额外信息是强制性要求的。尤其是与像素间的物理间距相关联的信息、图像在一些世界坐标中的空间位置信息极其重要。

图像原点、三维空间的体素（voxel）方向，和间距是很多应用的基础。例如配准就是在物理坐标中执行的。错误的定义间距、方向和原点在这样的处理中将导致不一致的结果。没有空间信息的医学图像不能被用于医学诊断、图像分析、特征提取、放射辅助治疗和图像指导手术。换句话说，缺乏空间信息的医学图像不仅仅是无用的，而且是很危险的。

上图说明了与itk::Image 相关的主要几何概念。在这个图中，圆圈表示像素的中心。像素值假定为像素中心的单位脉冲函数。像素间距是像素中心之间的距离，在各个坐标方向上可以是不同的。图像原点是图像中第一个像素的坐标。在这个简化示例中，体素点阵与物理空间方向完美对齐，因此图像方向是恒等映射。如果体素点阵采样相对于物理空间旋转，则图像方向将包含一个旋转矩阵。

一个像素就是含有数据值的像素中心周围的矩形区域，它可以被认为是图像网格中的Voronoi区域，如上图中右图所示。图像值的线性内插是在Delaunay区域中执行的。

图像间距用FixedArray表示，它的大小与图像维数相匹配。为了手动设定图像的间距，必须创建一个对应类型的数组。然后将这个数组的元素用相邻像素中心的间距来进行初始化。下面的程序代码阐述了图像类中处理原点和间距的方法。

ImageType::SpacingType spacing;
// Units (e.g., mm, inches, etc.) are defined by the application.
spacing[0] = 0.33; // spacing along X
spacing[1] = 0.33; // spacing along Y
spacing[2] = 1.20; // spacing along Z

使用SetSpacing( )方法来把数组分配给图像。image->SetSpacing( spacing );
使用GetSpacing( )方法可以从图像中提取间距信息。该方法返回一个FixedArray 的引用，它可以用来读取数组的内容。注意使用关键字const 表示数组是不可修改的。

const ImageType::SpacingType& sp = image->GetSpacing();
std::cout <<"Spacing = ";
std::cout <<sp[0] <<","<< sp[1] <<","<< sp[2] << std::endl;

图像原点的处理方式和间距是相似的。必须首先分配合适维数的Point。任何一个位置都可以指定为原点坐标。这些坐标和图像中第一个像素的位置相一致，图像可以用物理空间中的任意参考系统，用户必须确保同一个应用下的所有图像都是使用一致的参考系统。这在图像配准应用中极其重要。

下面的代码阐述了用于初始化图像原点的变量的创建和分配。

//coordinates of the center of the first pixel in N-D
ImageType::PointType newOrigin;
newOrigin.Fill(0.0);
image->SetOrigin( newOrigin );

使用GetOrigin( )方法也可以从图像中提取原点。这将返回Point的一个引用。该引用可以用来读取数组中的内容。再次提醒注意关键字const 表示数组内容是不可修改的。

const ImageType::PointType& origin = image->GetOrigin( );
std::cout << "Origin = ";
std::cout << origin[0] << ", "<< origin[1] << ", " << origin[2] << std::endl;

图像方向矩阵代表了图像采样和物理空间坐标系统之间的方向关系。图像方向矩阵是一个标准正交矩阵,是一个N*N矩阵，其中N是图像维数。

下面的代码阐述了用于初始化恒等图像方向的变量的创建和分配。

//coordinates of the center of the first pixel in N-D
ImageType::DirectionType direction;
direction.SetIdentity();
image->SetDirection(direction);

使用GetDirection( )方法也可以从图像中提取方向。这将返回Matrix的一个引用。该引用可以用来读取数组中的内容。注意关键字const 表示数组内容是不可修改的。

const ImageType::DirectionType&direct = image->GetDirection( );
std::cout << "Direction = "<< std::endl;
std::cout <<direct<< std::endl;

图像采样的间距、原点和方向一经初始化，图像就会正确映射像素索引到物理空间坐标，反之亦然。

下面的代码阐述了一个物理空间点是如何映射到图像索引的，以读取其邻近像素的内容。

首先，必须声明一个itk::Point 类型。这个Point模板类用来表示坐标的类型和空间维数。在这种特殊的情况下，Point 的维数必须和图像的维数相匹配。typedef itk::Point< double, ImageType::ImageDimension > PointType;像itk::Index一样，itk::Point类是一个相对较小、较简单的对象。这意味着这里不会使用itk::SmartPointer，而是像其它C++类一样简单地声明为对象实例。一旦声明Point，就可以使用传统的数组符号来访问它的元素。特别地，允许使用[]操作。出于效率原因，访问特定Point元素时不执行数组边界检查。用户必须确认数组索引在有效范围{0，Dimension−1}以内。

PointType point;
point[0] = 1.45; // x coordinate
point[1] = 7.21; // y coordinate
point[2] = 9.28; // z coordinate

图像将使用当前的原点和间距把point 映射到index。必须提供一个index 对象来接收映射结果。可以用图像类型中定义的IndexType来对index对象进行实例化：ImageType::IndexType pixelIndex;
图像类的TransformPhysicalPointToIndex( )方法可以计算出与提供的point 最接近的像素index。这个方法将检查这个index 是否被包含在当前的缓冲像素数据里。这个方法返回一个布尔值，表示这个index 结果是否在这个缓冲区间内。当返回值是false时，表示输出的index不能使用。
下面的代码阐述了point 到index 的映射，以及使用像素index访问图像像素数据的用法：

ImageType::IndexType pixelIndex;
const bool isInside = image->TransformPhysicalPointToIndex( point, pixelIndex );
if ( isInside )
{ImageType::PixelType pixelValue = image->GetPixel( pixelIndex );pixelValue += 5;image->SetPixel( pixelIndex, pixelValue );
}

记住GetPixel( )和SetPixel( )方法都是非常低效的访问像素数据的方法，当需要大量访问像素数据时应该使用图像迭代器。

下面的示例说明了对于一个指定图像，其图像索引位置和它对应的物理点表示之间的数学关系。
让我们想象存在一个图形用户接口，终端用户用鼠标手动选择了物体的左眼的体素索引位置。我们需要将那个索引位置转换为物理位置，以执行精确的激光引导手术。此时可以使用TransformIndexToPhysicalpoint方法。

const Image::IndexType LeftEyeIndex = GetIndexFromMouseClick();
ImageType::PointType LeftEyePoint;
Image->TransformIndexToPhysicalpoint(LeftEyeIndex, LeftEyePoint);

对于一个指定索引I3X1，物理位置P3X1计算如下：
P3X1=O3X1+D3X3∗diag(S3X1)3x3∗I3X1P3X1 = O3X1 +D3X3 ∗diag(S3X1)3x3 ∗I3X1P3X1=O3X1+D3X3∗diag(S3X1)3x3∗I3X1
其中D是一个标准正交方向余弦矩阵，S是图像空间对角矩阵。
相应的C/C++代码如下：

typedef itk::Matrix<double, Dimension, Dimension> MatrixType;
MatrixType SpacingMatrix;
SpacingMatrix.Fill( 0.0F );
const ImageType::SpacingType & ImageSpacing = image->GetSpacing();
SpacingMatrix( 0,0 ) = ImageSpacing[0];
SpacingMatrix( 1,1 ) = ImageSpacing[1];
SpacingMatrix( 2,2 ) = ImageSpacing[2];
const ImageType::DirectionType & ImageDirectionCosines =
image->GetDirection();
const ImageType::PointType &ImageOrigin = image->GetOrigin();
typedef itk::Vector< double, Dimension > VectorType;
VectorType LeftEyeIndexVector;
LeftEyeIndexVector[0]= LeftEyeIndex[0];
LeftEyeIndexVector[1]= LeftEyeIndex[1];
LeftEyeIndexVector[2]= LeftEyeIndex[2];
ImageType::PointType LeftEyePointByHand =
ImageOrigin + ImageDirectionCosines * SpacingMatrix * LeftEyeIndexVector;

2.代码实现

#include "itkImage.h"// 模拟鼠标点击获取图像像素值的函数
static itk::Image< unsigned short, 3 >::IndexType GetIndexFromMouseClick()
{itk::Image< unsigned short, 3 >::IndexType LeftEyeIndex;LeftEyeIndex[0]=30;LeftEyeIndex[1]=12;LeftEyeIndex[2]=9;return LeftEyeIndex;
}
int main(int, char *[])
{const unsigned int Dimension=3;typedef itk::Image< unsigned short, Dimension > ImageType;ImageType::Pointer image = ImageType::New();const ImageType::SizeType  size  = {{ 100, 100, 100}}; //Size along {X,Y,Z}const ImageType::IndexType start = {{ 0, 0, 0 }}; // First index on {X,Y,Z}ImageType::RegionType region;region.SetSize( size );region.SetIndex( start );image->SetRegions( region );image->Allocate(true); // 内存分配//图像类中处理原点和间距的方法//创建一个和图像数据类型相一致的数列spacingImageType::SpacingType spacing;//设定X、Y、Z方向像素间的间距spacing[0] = 0.33; // X方向上相邻两像素中心的间距spacing[1] = 0.33; // Y方向上相邻两像素中心的间距spacing[2] = 1.20; // Z方向上相邻两像素中心的间距//使用 SetSpacing( ) 方法指向数列spacingimage->SetSpacing( spacing );//使用 GetSpacing( ) 方法可以从图像中得到间距信息, const 表示数列是不可修改的const ImageType::SpacingType& sp = image->GetSpacing();//输出读取到的图像X、Y、Z方向像素间的间距信息std::cout << "Spacing = ";std::cout << sp[0] << ", " << sp[1] << ", " << sp[2] << std::endl;//初始化图像原点的变量newOrigin的创建和分配ImageType::PointType newOrigin;newOrigin.Fill(0.0);image->SetOrigin( newOrigin );// GetOrigin( ) 方法可以从图像中读取原点const ImageType::PointType & origin = image->GetOrigin();//输出读取到图像的原点坐标std::cout << "Origin = ";std::cout << origin[0] << ", "<< origin[1] << ", "<< origin[2] << std::endl;ImageType::DirectionType direction;direction.SetIdentity();  //设置为恒等映射image->SetDirection( direction );const ImageType::DirectionType& direct = image->GetDirection();std::cout << "Direction = " << std::endl;std::cout << direct << std::endl;//将物理空间映射到读取最近像素内容的图像index中//声明一个 itk::Point 类型。这个Point模板类用来表示坐标的类型和空间维数。typedef itk::Point< double, ImageType::ImageDimension > PointType;PointType point;point[0] = 1.45;    // x coordinatepoint[1] = 7.21;    // y coordinatepoint[2] = 9.28;    // z coordinate//图像类型中定义的 IndexType 来对 index 对象进行实例化ImageType::IndexType pixelIndex;// Point 到 index 的映射和访问图像像素数据的像素 index 的用法const bool isInside =image->TransformPhysicalPointToIndex( point, pixelIndex );if ( isInside ){ImageType::PixelType pixelValue = image->GetPixel( pixelIndex );pixelValue += 5;image->SetPixel( pixelIndex, pixelValue );}//将图像索引转为物理空间坐标const ImageType::IndexType LeftEyeIndex = GetIndexFromMouseClick();ImageType::PointType LeftEyePoint;//(1)调用函数实现将图像索引转为物理空间坐标image->TransformIndexToPhysicalPoint(LeftEyeIndex,LeftEyePoint);std::cout << "===========================================" << std::endl;std::cout << "The Left Eye Location is " << LeftEyePoint << std::endl;//(2)手动实现将图像索引转为物理空间坐标typedef itk::Matrix<double, Dimension, Dimension> MatrixType;MatrixType SpacingMatrix;SpacingMatrix.Fill( 0.0F );const ImageType::SpacingType & ImageSpacing = image->GetSpacing();SpacingMatrix( 0,0 ) = ImageSpacing[0];SpacingMatrix( 1,1 ) = ImageSpacing[1];SpacingMatrix( 2,2 ) = ImageSpacing[2];const ImageType::DirectionType & ImageDirectionCosines =image->GetDirection();const ImageType::PointType &ImageOrigin = image->GetOrigin();typedef itk::Vector< double, Dimension > VectorType;VectorType LeftEyeIndexVector;LeftEyeIndexVector[0]= LeftEyeIndex[0];LeftEyeIndexVector[1]= LeftEyeIndex[1];LeftEyeIndexVector[2]= LeftEyeIndex[2];ImageType::PointType LeftEyePointByHand =ImageOrigin + ImageDirectionCosines * SpacingMatrix * LeftEyeIndexVector;// Software Guide : EndCodeSnippetstd::cout << "===========================================" << std::endl;std::cout << "Spacing:: " << std::endl << SpacingMatrix << std::endl;std::cout << "===========================================" << std::endl;std::cout << "DirectionCosines:: " << std::endl << ImageDirectionCosines << std::endl;std::cout << "===========================================" << std::endl;std::cout << "Origin:: " << std::endl << ImageOrigin << std::endl;std::cout << "===========================================" << std::endl;std::cout << "The Left Eye Location is " << LeftEyePointByHand << std::endl;if ( (LeftEyePointByHand - LeftEyePoint).GetNorm() < 0.01F ){std::cout << "===========================================" << std::endl;std::cout << "Two results are identical as expected!" << std::endl;std::cout << "The Left Eye from TransformIndexToPhysicalPoint is " << LeftEyePoint << std::endl;std::cout << "The Left Eye from Math is " << LeftEyePointByHand << std::endl;}return EXIT_SUCCESS;
}

3.结果

参考目录

https://blog.csdn.net/qq_32809093/article/details/116886173