医学影像与坐标系转换 （以DICOM为例）_dicom矩阵与图像位置转换

医学影像与坐标系转换

研究DICOM坐标系转换，查阅网上资料时，发现大多自用文章存在术语不统一，概念混淆的情况，这里进行简单的总结，如有错误望请指出

一、常用医学影像格式

下面我将简单介绍常用的医学影像文件格式。

医学数字成像和通信(DICOM)

DICOM（Digital Imaging and Communications in Medicine）是一种用于医学图像处理的国际标准。它定义了一种格式，用于存储、传输、共享和打印医学图像信息。DICOM不仅包含图像数据，还包含了与图像相关的患者信息、设备信息和诊断信息。

主要特点

• 多模态支持：DICOM支持多种医学成像技术，包括CT（计算机断层扫描）、MRI（磁共振成像）、PET（正电子发射断层扫描）、超声等。
• 包含丰富的元数据：DICOM文件不仅包含图像数据，还包含了大量的元数据，如患者信息（姓名、性别、出生日期等）、图像采集信息（设备类型、扫描参数等）和诊断信息（报告、注释等）。
• 网络通信：DICOM定义了一套网络协议，用于在医疗设备和系统之间传输DICOM文件。这使得医生可以在不同的设备和地点查看和分析图像。
• 数据安全：DICOM支持数据加密和用户认证，以保护患者的隐私和数据的安全。
• 扩展性：DICOM是一个开放的标准，可以根据需要进行扩展。例如，可以添加新的图像类型，或者定义新的元数据字段。

数据结构

• 文件头：包含了DICOM的文件标识符和元数据，如患者信息、设备信息、图像采集参数等。
• 数据集：包含了实际的图像数据和与图像相关的元数据。数据集由一系列的数据元素组成，每个数据元素都有一个唯一的标签（Tag），用于标识数据元素的类型，如患者姓名、图像日期、图像时间等。
  • 数据元素(dataElement)：数据元素是DICOM数据结构的基本单位，每个数据元素都包含一个标签(Tag)、一个值长度和一个或多个值。标签用于标识数据元素的类型，值长度用于指示值的长度，值则是数据元素的实际内容。
    • 序列：特殊的数据元素，它的值是一个或多个数据集。序列可以用于表示复杂的数据结构，如图像序列、报告结构等。
    • 像素数据：特殊的数据元素，它包含了图像的实际像素值。像素数据通常是数据集中最大的部分。

这种数据结构使得DICOM可以表示各种复杂的医学图像数据，包括2D图像、3D体积数据、时间序列数据等。同时，通过使用元数据，DICOM还可以包含与图像相关的丰富的上下文信息，如患者信息、设备信息、诊断信息等。

数据元素实例

通过 Group & Element 的组合，确定了唯一的Tag，我们可以通过该Tag获取该类型数据元素中的信息。举例的Tag,后面我们会详细介绍，将用于坐标系的转换。

神经影像学信息技术倡议(NIfTI)

NIfTI（Neuroimaging Informatics Technology Initiative）格式是一种用于存储医学图像数据的文件格式，特别是用于神经影像学的数据。它是由NIH（美国国立卫生研究院）的NIfTI项目组开发的，目的是提供一种简单、灵活且方便的方式来存储和共享大量的神经影像学数据。

主要特点

• 支持多维数据：NIfTI格式支持存储多维数据，包括2D图像、3D体积数据和4D时间序列数据。这使得NIfTI格式非常适合用于存储MRI、fMRI和PET等类型的数据。
• 包含丰富的元数据：NIfTI文件包含了大量的元数据，如图像维度、像素间距、时间间隔、坐标系信息等。这些元数据使得NIfTI文件可以自包含，不需要额外的信息就可以正确地解析和显示图像。
• 支持数据压缩：NIfTI格式支持gzip压缩，可以显著减小文件大小，方便数据的存储和传输。
• 兼容Analyze 7.5格式：NIfTI格式是基于Analyze 7.5格式开发的，因此它与Analyze 7.5格式兼容。这意味着大多数可以读取Analyze 7.5格式的软件也可以读取NIfTI格式。

NIfTI格式通常有两种文件扩展名：.nii和.nii.gz。.nii文件是未压缩的NIfTI文件，而.nii.gz文件是gzip压缩后的NIfTI文件。

数据结构

数据主要由两部分组成：头部（header）和数据部分（data）。

1. 头部（Header）：头部包含了图像的元数据，如图像维度、像素间距、时间间隔、坐标系信息等。头部的大小固定为348字节，这使得读取和解析头部信息变得非常简单。
2. 数据部分（Data）：数据部分包含了图像的实际像素值。数据可以是二维的（例如CT或MRI图像），三维的（例如一个时间点的3D体积数据），或者四维的（例如随时间变化的3D体积数据，如fMRI数据）。

NIfTI格式还支持存储多个数据集在一个文件中，这通过在头部中指定维度信息来实现。例如，一个四维的NIfTI文件可以存储一个随时间变化的3D体积数据序列。

头文件数据实例

// NIFTI-1 头文件存储信息（顺序存放）
/*
NIFTI 第一版的数据格式的头文件大小必须为348字节，这是为了兼容新，和旧的analyze格式保持一致。就是sizeof_hdr保存的数据大小。
*/
//                       //OFFSET    SIZE    Description
int     sizeof_hdr;      //0B        4B      Size of the header.Must be 348 (bytes).
char    data_type[10];   //4B        10B     Not used; compatibility with analyze.
char    db_name[18];     //14B       18B     Not used; compatibility with analyze.
int     extents;         //32B       4B      Not used; compatibility with analyze.
short   session_error;   //36B       2B      Not used; compatibility with analyze.
char    regular;         //38B       1B      Not used; compatibility with analyze.
char    dim_info;        //39B       1B      Encoding directions(phase, frequency, slice).
short   dim[8];          //40B       16B     Data array dimensions.
float   intent_p1;       //56B       4B      1st intent parameter.
float   intent_p2;       //60B       4B      2nd intent parameter.
float   intent_p3;       //64B       4B      3rd intent parameter.
short   intent_code;     //68B       2B      nifti intent.
short   datatype;        //70B       2B      Data type.
short   bitpix;          //72B       2B      Number of bits per voxel.
short   slice_start;     //74B       2B      First slice index.
float   pixdim[8];       //76B       32B     Grid spacings(unit per dimension).
float   vox_offset;      //108B      4B      Offset into a.nii file! start point.
float   scl_slope;       //112B      4B      Data scaling, slope.
float   scl_inter;       //116B      4B      Data scaling, offset.
short   slice_end;       //120B      2B      Last slice index.
char    slice_code;      //122B      1B      Slice timing order.
char    xyzt_units;      //123B      1B      Units of pixdim[1..4].
float   cal_max;         //124B      4B      Maximum display intensity.
float   cal_min;         //128B      4B      Minimum display intensity.
float   slice_duration;  //132B      4B      Time for one slice.
float   toffset;         //136B      4B      Time axis shift.
int     glmax;           //140B      4B      Not used; compatibility with analyze.
int     glmin;           //144B      4B      Not used; compatibility with analyze.
char    descrip[80];     //148B      80B     Any text.
char    aux_file[24];    //228B      24B     Auxiliary filename.
short   qform_code;      //252B      2B      Use the quaternion fields.
short   sform_code;      //254B      2B      Use of the affine fields.
float   quatern_b;       //256B      4B      Quaternion b parameter.
float   quatern_c;       //260B      4B      Quaternion c parameter.
float   quatern_d;       //264B      4B      Quaternion d parameter.
float   qoffset_x;       //268B      4B      Quaternion x shift.
float   qoffset_y;       //272B      4B      Quaternion y shift.
float   qoffset_z;       //276B      4B      Quaternion z shift.
float   srow_x[4];       //280B      16B     1st row affine transform
float   srow_y[4];       //296B      16B     2nd row affine transform.
float   srow_z[4];       //312B      16B     3rd row affine transform.
char    intent_name[16]; //328B      16B     Name or meaning of the data.
char    magic[4];        //344B      4B      Magic string.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

上面介绍的是NIfTI-1，NIfTI-2是NIfTI-1的拓展版本，不做过多介绍。
以下NIfTI-2的头文件参考，如果magic字段的值为"n+2\0"，则文件是NIfTI-2格式。
// NIfTI-2 header
typedef struct {
    int sizeof_hdr;     // 0B   8B  Size of the header. Must be 540 (bytes).
    char magic[8];      // 8B   8B  Magic string.
    int16_t datatype;   // 16B  2B  Data type.
    int16_t bitpix;     // 18B  2B  Number of bits per voxel.
    int64_t dim[8];     // 20B  64B Data array dimensions.
    double intent_p1;   // 84B  8B  1st intent parameter.
    double intent_p2;   // 92B  8B  2nd intent parameter.
    double intent_p3;   // 100B 8B  3rd intent parameter.
    double pixdim[8];   // 108B 64B Grid spacings (unit per dimension).
    double vox_offset;  // 172B 8B  Offset into .nii file (start point of data).
    double scl_slope;   // 180B 8B  Data scaling, slope.
    double scl_inter;   // 188B 8B  Data scaling, offset.
    double cal_max;     // 196B 8B  Maximum display intensity.
    double cal_min;     // 204B 8B  Minimum display intensity.
    double slice_duration; // 212B 8B  Time for one slice.
    double toffset;     // 220B 8B  Time axis shift.
    int64_t slice_start; // 228B 8B  First slice index.
    int64_t slice_end;  // 236B 8B  Last slice index.
    char descrip[80];   // 244B 80B Any text.
    char aux_file[24];  // 324B 24B Auxiliary filename.
    int qform_code;     // 348B 4B  Use the quaternion fields.
    int sform_code;     // 352B 4B  Use of the affine fields.
    double quatern_b;   // 356B 8B  Quaternion b parameter.
    double quatern_c;   // 364B 8B  Quaternion c parameter.
    double quatern_d;   // 372B 8B  Quaternion d parameter.
    double qoffset_x;   // 380B 8B  Quaternion x shift.
    double qoffset_y;   // 388B 8B  Quaternion y shift.
    double qoffset_z;   // 396B 8B  Quaternion z shift.
    double srow_x[4];   // 404B 32B 1st row affine transform.
    double srow_y[4];   // 436B 32B 2nd row affine transform.
    double srow_z[4];   // 468B 32B 3rd row affine transform.
    int slice_code;     // 500B 4B  Slice timing order.
    int xyzt_units;     // 504B 4B  Units of pixdim[1..4].
    int intent_code;    // 508B 4B  nifti intent.
    char intent_name[16]; // 512B 16B Name or meaning of the data.
    char dim_info;      // 528B 1B  Encoding directions (phase, frequency, slice).
    char unused_str[15]; // 529B 15B Unused, filled with zeroes.
} nifti_2_header;      // Total: 544B
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

NIfTI与DICOM格式的区别

NIfTI和DICOM都是医学影像领域常用的文件格式，但它们有一些重要的区别：

• 数据结构：DICOM是一种复杂的文件格式，它不仅包含医学影像数据，还包含了大量的元数据，如患者信息、设备信息、扫描参数等。这使得DICOM能够存储丰富的信息，但也使得它的结构相对复杂。相比之下，NIfTI是一种更简单的文件格式，它主要用于存储三维或四维的医学影像数据，其元数据主要包括图像的维度、像素尺寸、数据类型等。
• 使用场景：DICOM是医院和诊所中最常用的文件格式，它被广泛用于存储和传输CT、MRI等医学影像数据。NIfTI则主要用于科研和软件开发中，它被广泛用于神经影像学和医学影像处理软件中。
• 坐标系统：DICOM使用的是LPS（左-后-上）坐标系统，而NIfTI使用的是RAS（右-前-上）坐标系统。这意味着在将DICOM数据转换为NIfTI格式时，我们需要进行坐标转换。
• 文件数量：一个DICOM序列通常由多个DICOM文件组成，每个文件包含一张图像。相比之下，一个NIfTI文件可以包含一个完整的三维或四维图像数据集，这使得NIfTI文件更易于管理和处理。
• 数据类型支持：NIfTI支持更多的数据类型，包括复数和RGB/RGBA数据，这在DICOM中是不支持的。

近乎原始光栅数据(Nrrd)

近乎原始光栅数据（Nearly Raw Raster Data，NRRD）是一种用于表示n维光栅数据的文件格式。它由Gordon Kindlmann在开发Teem库时设计，用于存储和处理医学影像数据和体渲染数据，被广泛用于医学影像处理和体渲染软件中，如3D Slicer、MITK等。

主要特点

• 简单性：NRRD格式设计简单，易于理解和实现。它的头部信息是纯文本，可以直接用文本编辑器查看。
• 灵活性：NRRD格式可以表示各种类型和维度的光栅数据，包括标量数据、向量数据、张量数据等。它支持多种数据类型，包括无符号整数、有符号整数、浮点数等。
• 元数据支持：NRRD格式可以存储各种元数据，如像素尺寸、坐标轴信息、数据单位等。
• 数据压缩：NRRD格式支持gzip压缩，可以有效地减小文件大小。

数据结构

NRRD文件通常包括一个头文件和一个数据文件。头文件是纯文本文件，包含了元数据和数据文件的路径。数据文件包含了实际的光栅数据，可以是二进制格式，也可以是文本格式。此外，NRRD也支持将头文件和数据文件合并为一个文件，这种情况下，头部信息和光栅数据之间由一个空行分隔。

• 头部信息：头部信息是纯文本，包含了元数据和数据文件的路径。元数据包括：
  • dimension：数据的维度，例如2D图像的维度为2，3D图像的维度为3。
  • sizes：每个维度的大小，例如一个256x256的2D图像的sizes为256 256。
  • type：数据类型，例如unsigned char、short、float等。
  • endian：数据的字节序，例如little或big。
  • encoding：数据的编码方式，例如raw、ascii、gzip等。
  • space：空间坐标系，例如left-posterior-superior (LPS)或right-anterior-superior (RAS)。
  • space directions：每个维度的方向向量。
  • space origin：空间坐标系的原点。
• 光栅数据：光栅数据是实际的图像数据，可以是二进制格式，也可以是文本格式。数据的排列方式由头部信息中的dimension和sizes字段决定。

二、参照系

在处理医学图像数据时，通常会首先将数据从像素坐标系转换到患者坐标系，然后再转换到世界坐标系。这样可以更方便地描述和分析患者的解剖结构，同时也可以方便地与其他设备或系统交互。

世界坐标系

存在一个场景，当我们需要描述其中绝对位置和方向时，约定一个以点O为世界原点，于O点彼此正交的三个单位轴形成的笛卡尔坐标系为这个空间的世界坐标系。世界坐标系只是无数可能的坐标系之一，但是我们需要一个坐标系来测量默认情况下的所有事物，所以我们创建了它并给它起了“世界坐标系”这个特殊名称。

解剖学坐标系

医学影像以患者为主体，因此解剖学坐标系/患者坐标系极其重要。

首先我们需要了解医学上的六个指向和三个观察位：

六个指向：

• Anteriort：前
• Posterior：后
• Left：左
• Right：右
• Superior/Head：上/头
• Inferior/Feet：下/脚

三个观察位：

• 冠状位（Coronal）：体前/体后观察人体（A-P），紫色面称为冠状面
• 矢状位（Sagittal）：体侧观察人体（ L-R），粉色面称为矢状面
• 横断位（Transverse）：头顶俯视/足底仰视人体（S-I），绿色面称为横断面或水平面

所谓的解剖坐标系，这是一个相对于特定患者的坐标系，通常用于描述患者身体内部结构的位置和方向。在患者坐标系中，原点通常被定义为患者身体的某个特定点，方向通常取上述六个方向其三。如常见的LPS坐标系，RAS坐标系。

像素坐标系

这是一个相对于特定图像的坐标系，通常用于描述图像中像素的位置。在像素坐标系中，原点通常被定义为图像的左上角，X轴从左到右，Y轴从上到下。（很多文章将像素坐标系与图像坐标系混淆，请注意）

三、坐标系的转换

为什么进行坐标系转换

• 不同的设备和软件使用不同的坐标系：在医学图像处理中，我们经常需要处理来自不同设备和软件的数据，这些设备和软件可能使用不同的坐标系，为了在这些设备和软件之间共享和处理图像数据，我们需要进行坐标系转换。例如，MRI、CT、PET等医学成像设备通常使用DICOM格式来存储图像数据，而DICOM通常采用LPS（左后上）坐标系。另一方面，一些常用的图像文件格式，如NIfTI和Nrrd，通常使用RAS（右前上）坐标系。。
• 方便数据处理和分析：在某些情况下，将图像数据转换到特定的坐标系可以简化数据处理和分析的过程。例如，将图像数据转换到与人体解剖结构对应的坐标系，可以方便我们理解和解释图像数据。
• 配准和融合：在进行图像配准（将两个或多个图像对齐）或图像融合（将两个或多个图像合并）时，我们需要将所有的图像转换到同一坐标系。
• 可视化：在进行图像可视化时，我们可能需要将图像数据转换到与显示设备对应的坐标系，例如屏幕坐标系或世界坐标系

医学图像处理软件及默认坐标系

1. 3D Slicer: 默认使用RAS（右前上）坐标系。这是一种常见的神经影像坐标系，其中X轴指向右，Y轴指向前，Z轴指向上。

2. MITK (Medical Imaging Interaction Toolkit): MITK默认使用LPS（左后上）坐标系，这是DICOM标准中定义的坐标系，其中X轴指向左，Y轴指向后，Z轴指向上。

3. FSL (FMRIB Software Library): FSL默认使用RAS（右前上）坐标系。

4. SPM (Statistical Parametric Mapping): SPM默认使用RAS（右前上）坐标系。

5. OsiriX: OsiriX默认使用LPS（左后上）坐标系。

请注意，虽然这些软件有默认的坐标系，但大多数都支持在不同的坐标系之间转换，以便处理来自不同设备和软件的图像数据。

以DICOM为例

在医学图像处理中，我们通常会设定一个世界坐标系来描述物体的绝对位置和方向。这个世界坐标系通常会选择一个与患者身体结构相关的坐标系，如LPS（左后上）或RAS（右前上）坐标系，因为这样可以更直观地理解和解释图像数据。因此，我们可以说，在医学图像处理的上下文中，解剖学坐标系（也就是LPS或RAS坐标系）实际上就是一种世界坐标系。而在DICOM中，被称为参考坐标系（Reference Coordinate System）的世界坐标系以LPS为标准。

当我们进行坐标系的变换时，所涉及的旋转、缩放、翻转、平移等操作，都称之为仿射变换，变换使用的矩阵称为仿射矩阵。在DICOM中，我们可以通过以下步骤构建从像素坐标系转到物理坐标系所要用的仿射矩阵：

1. 获取DICOM标签值：我们需要以下三个标签的值：
   • Image Position (Patient) (0020,0032)：这个标签包含三个值，分别表示图像平面左上角像素（也就是原点）在物理空间中的x、y、z坐标。这些值是相对于患者身体的，通常使用LPS（左后上）坐标系。
   • Image Orientation (Patient) (0020,0037)：这个标签包含六个值，前三个值表示图像行（row）的方向，后三个值表示图像列（column）的方向。这些值是相对于患者身体的，通常使用LPS（左后上）坐标系。每个方向都是一个单位向量，表示该方向与x、y、z轴的夹角的余弦值。
   • Pixel Spacing (0028,0030)：这个标签包含两个值，第一个值表示行间距（即垂直方向上相邻两像素的物理距离），第二个值表示列间距（即水平方向上相邻两像素的物理距离）。这些值通常以毫米为单位。
2. 计算行和列向量：我们将Image Orientation (Patient)标签的前三个值和后三个值分别乘以Pixel Spacing标签的值，得到行和列向量。
3. 计算切片向量：我们计算行和列向量的叉积，得到切片向量。
4. 构建仿射矩阵：我们将行向量、列向量和切片向量作为仿射矩阵的前三列，将Image Position (Patient)标签的值作为仿射矩阵的最后一列。最后一行是[0, 0, 0, 1]。

以下是最终的矩阵：

这个仿射矩阵可以将像素坐标系中的坐标转换到物理坐标系中。例如，如果我们有一个像素坐标P(x, y, z)，我们可以通过将P乘以这个仿射矩阵，得到物理坐标系中的坐标。

在C++中，我们可以使用VTK库与ITK库来处理DICOM文件并计算像素坐标系到物理坐标系的转换。以下是一个示例代码：

#include <itkImageFileReader.h>
#include <itkGDCMImageIO.h>
#include <itkSpatialOrientationAdapter.h>
#include <vtkMatrix4x4.h>
#include <vtkMath.h>

using ImageType = itk::Image<unsigned short, 3>;
using ReaderType = itk::ImageFileReader<ImageType>;

// 创建变换矩阵的函数
vtkSmartPointer<vtkMatrix4x4> CreateTransformMatrix(const ReaderType::Pointer &reader);
// 获取图像方向的函数
itk::SpatialOrientation::ValidCoordinateOrientationFlags GetImageOrientation(const ReaderType::Pointer &reader);
// 如果方向不是LPS，转换矩阵的函数
vtkSmartPointer<vtkMatrix4x4> TransformMatrixIfNotLPS(vtkSmartPointer<vtkMatrix4x4> matrix, itk::SpatialOrientation::ValidCoordinateOrientationFlags orientationFlag);

int main(int argc, char* argv[])
{
    // 检查命令行参数
    if (argc < 2)
    {
        std::cerr << "Usage: " << argv[0] << " <DICOM file>" << std::endl;
        return 1;
    }
    
           // 创建ITK读取器并读取图像
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName(argv[1]);
    reader->Update();
    
           // 根据源数据获取变换矩阵
    vtkSmartPointer<vtkMatrix4x4> transformMatrix = CreateTransformMatrix(reader);
    
           // 获取图像坐标系
    itk::SpatialOrientation::ValidCoordinateOrientationFlags orientationFlag 
        = GetImageOrientation(reader);
    
           // 如果坐标系不是LPS，进行矩阵变换计算
    if (orientationFlag != itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS)
    {
        transformMatrix = TransformMatrixIfNotLPS(transformMatrix, orientationFlag);
    }
    
           // 打印变换矩阵
    transformMatrix->Print(std::cout);
    
    // 使用变换矩阵从像素坐标系转到LPS坐标系
    // double pixelPoint[4] = {row, column, slice, 1};
    // double lpsPoint[4];
    // transformMatrix->MultiplyPoint(pixelPoint, lpsPoint);
    
    return 0;
}

vtkSmartPointer<vtkMatrix4x4> CreateTransformMatrix(const ReaderType::Pointer &reader)
{
    if (!reader)
    {
        std::cout << "Invalid reader pointer!" << std::endl;
        return nullptr;
    }
    // 获取图像位置、方向和间距
    ImageType::PointType imagePosition = reader->GetOutput()->GetOrigin();
    ImageType::DirectionType imageOrientation = reader->GetOutput()->GetDirection();
    ImageType::SpacingType pixelSpacing = reader->GetOutput()->GetSpacing();
    
           // 创建变换矩阵
    vtkSmartPointer<vtkMatrix4x4> transformMatrix =
        vtkSmartPointer<vtkMatrix4x4>::New();
    
           // 计算行向量和列向量
    double rowVector[3], columnVector[3];
    for (int i = 0; i < 3; ++i)
    {
        rowVector[i] = imageOrientation[i][0] * pixelSpacing[0];
        columnVector[i] = imageOrientation[i][1] * pixelSpacing[1];
    }
    
           // 在矩阵中设置行向量、列向量和图像位置
    for (int i = 0; i < 3; ++i)
    {
        transformMatrix->SetElement(i, 0, rowVector[i]);
        transformMatrix->SetElement(i, 1, columnVector[i]);
        transformMatrix->SetElement(i, 3, imagePosition[i]);
    }
    
           // 计算行向量和列向量的叉积以获取切片方向
    double sliceDirection[3];
    vtkMath::Cross(rowVector, columnVector, sliceDirection);
    
           // 在矩阵中设置切片方向
    for (int i = 0; i < 3; ++i)
    {
        transformMatrix->SetElement(i, 2, sliceDirection[i]);
    }
    
           // 设置矩阵的最后一个元素
    transformMatrix->SetElement(3, 0, 0);
    transformMatrix->SetElement(3, 1, 0);
    transformMatrix->SetElement(3, 2, 0);
    transformMatrix->SetElement(3, 3, 1);
    
    return transformMatrix;
}

itk::SpatialOrientation::ValidCoordinateOrientationFlags GetImageOrientation(ReaderType::Pointer reader)
{
    if (!reader)
    {
        std::cout << "Invalid reader pointer!" << std::endl;
        return itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_INVALID;
    }
    // 获取图像方向
    ImageType::DirectionType direction = reader->GetOutput()->GetDirection();
    
           // 将图像方向转换为itk::SpatialOrientation::ValidCoordinateOrientationFlags值
    itk::SpatialOrientationAdapter adapter;
    itk::SpatialOrientation::ValidCoordinateOrientationFlags orientationFlag = adapter.FromDirectionCosines(direction);
    
    return orientationFlag;
}

vtkSmartPointer<vtkMatrix4x4> TransformMatrixIfNotLPS(vtkSmartPointer<vtkMatrix4x4> matrix, itk::SpatialOrientation::ValidCoordinateOrientationFlags orientationFlag)
{
    if (!matrix)
    {
        std::cout << "Invalid matrix pointer!" << std::endl;
        return nullptr;
    }
    // 使用itk::SpatialOrientationAdapter将orientationFlag转换为方向余弦矩阵
    itk::SpatialOrientationAdapter adapter;
    itk::SpatialOrientation::ValidCoordinateOrientationFlags lpsFlag = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS;
    itk::SpatialOrientationAdapter::DirectionType directionCosines = adapter.ToDirectionCosines(lpsFlag);
    
           // 直接修改原始矩阵的元素，将其转换为LPS坐标系
    for (int i = 0; i < 3; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            matrix->SetElement(i, j, directionCosines[i][j]);
        }
    }
    
    return matrix;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

请注意，我们不是总能从DICOM获取到想要的Tag，当有些Tag设备无法提供的时候，我们可以根据其他Tag进行origin、spacing、orientation的计算。这需要我们对于具体数据具体分析。

因为我的软件都是将NIfTI 与Nrrd转换为DICOM再使用，关于这两种格式文件的解析如果以后有时间我再补充。

最后，因为个人时间与水平有限，部分知识、代码未进行完整的验证与测试，请大家酌情参考使用。