DICOM世界观·第一章 坐标系统·番外篇

题记：

DICOM世界观·第一章 坐标系统完成后，总感觉缺了点什么，大概有两个原因：第一，没有从基础概念说起，来形象的介绍坐标系间的各种变换；第二，没有深入到DICOM数据本身，来进行实例演示。这两方面的介绍都停留在半山腰，让读者似懂非懂或一知半解。为此近期重新翻阅了一下经典著作《Introduction to Linear Algebra, Fifth Edition(2016)》，以及相关维基百科资料。虽不能完全将书中理念呈现于此篇博文，但还是希望能够对上文有一点补充，能够更清晰、形象的给大家一个展示，为后续DICOM的相关数据操作打下坚持的基础。所以本篇记做DICOM世界观 第一章 坐标系统·番外篇。

配图选自《Introduction to Linear Algebra, 4th edition》封面

第一章 坐标系统·番外篇

1.4 变换

日常生活中我们切身体会的“变换”有平移、缩放、旋转、挤压（这个可以想象一个正方形的边框，用力挤压任意一个顶点后会变成平行四边形的场景）、扭曲等等。通过这些变换我们可以将一个物体（备注：这里不一定指的是一个真真切切的实体，可能是一个图形或图像）变成我们想要的样子。更复杂的情况参见下图演示的Photoshop中的几种变换（快捷键：Ctrl+T）：

1.4.1 线性变换

在日常的变换中，线性变换最被我们所熟知。用数学的语言来描述线性变换（在wiki百科中，线性映射linear map、线性变换linear transformation、线性函数linear function在某些情况下所表达的意思相同）拥有以下特性：
- 加法闭合：
- f(u+v)=f(u)+f(v),可扩展为f(c1v1+c2v2+c3v3+...+cnvn)=c1f(v1)+c2f(v2)+c3f(v3)+...+cnf(vn)
-“数”乘闭合：f(cu)=c(fu)

上述两个定律满足闭合的意思指的是“上述两种变换的结果与变换之前的值同属于一个坐标系空间”。言外之意可以简单的理解为：变换前后直线依然是直线，直线彼此间的比例不变，且相对的坐标原点固定不变。(详情参考知乎马同学的介绍如何通俗的解释仿射变换)
用矩阵形式（这里仅用二维空间的变换来演示）来表示线性变化大致有以下几种情况（参见wiki百科linear map）：
(1)旋转：rotation

cosθsinθ−sinθcosθ

上述矩阵表示旋转任意角度θ,例如旋转90°矩阵为：

01−10

用matlab进行演示：

%旋转角度为30度，即pi/6r1 =    0.8660   -0.5000    0.5000    0.8660%旋转角度(默认正方向是逆时针)为90度，即pi/2r2 =     0    -1     1     0

(2)缩放：scaling

a00b

两个坐标轴的缩放比例分别为x为缩放比例等于a，y缩放比例等于b。例如下图一个是等比例扩大两倍，一个是X与Y轴扩大不同倍数。

ss =     2     0     0     2ss2 =     3     0     0     2

(3)剪切：shear

10m1

或

1m01

两种情况分别表示裁切的方向沿着x轴(因为左乘运算时保持y不变)和y轴（因为左乘运算时保持x不变），例如m=2时，

s1 =     1     2     0     1s2 =     1     0     2     1

(4)挤压：squeeze

k001/k

或

1/k00k

当k=2时

(5)投影：projection

0001

或

1000

分别表示投影到y轴和投影到x轴，例如：

(6)反射：reflection

−1001

或

100−1

分别表示沿着X轴方向(Y轴不变)反射和沿着Y轴方向（X轴不变）反射，例如：

1.4.2 仿射变换

在上述线性变换中可以发现，整个系统的坐标系没有动过（言外之意是指没有对任意点进行移动），因此可知线性变换是无法表达点的平移的（这一点有一些反常识大家细细体会一下）。为了包含平移这一日常最最常见的操作，我们需要引入“仿射变换”的概念，即Affine transformation。那么相较于线性变换，仿射变换的特性有：
- p+0=p
- (p+a)+b=p+(a+b)
- q∈A,存在a∈A,使得q=p+a

通俗的来讲，仿射变换后直线依然还是直线，直线间的比例保持不变，但是原点可能会发生变化，即仿射空间可以看作是没有原点的坐标系。(详情参考知乎马同学的介绍如何通俗的解释仿射变换)

上一节的线性变换用数学表达式表示为：

y=Ax,其中y,x是向量,A是线性变换矩阵

而仿射变换的表达式为：

y=Ax+b,其中y,x,b是向量,A是线性变换矩阵

这其中多出了一个b平移向量，所以没有办法保证f(x+y)=f(x)+f(y)，——说明仿射变换不满足线性结合律的，也正以为此才称之为仿射变换，是线性变换的扩展泛化。（参考阅读博文仿射空间与仿射变换）
为了将平移向量b包含到转换矩阵A中，这里需要用到“齐次坐标”的思路，即提升一个维度，将二维的平移用三维的线性变换来实现，言外之意高维空间的线性变换可以表示低维空间的仿射变换（备注：关于vector向量、matrix矩阵、space空间之间的关系可以仔细阅读经典著作《Introduction to Linear Algebra, Fifth Edition(2016)》）
即

y=Ax+b可以写作[y1]=[A0b1][x1]

用wiki百科的一幅图演示效果如下：

即先通过升维将二维的平面图形变成三维的体，然后对三维的体进行线性变换，最后进行降维处理得到二维平面图形的平移。

1.5 旋转

这里再一次对之前DICOM世界观·第一章 坐标系统博文介绍的各类描述旋转的方式以及本章节1.4 变换小节中的旋转进行单独的补充介绍。旋转的特殊性在于根据欧拉旋转定理，如果找到物体旋转中心（即不变的点），将该点移动到坐标原点，那么除却中心点到原点的平移操作，可以用旋转来表示空间中的任意变化。因此旋转有必要细细钻研一番。三维旋转，Rotation formalisms in three dimensions，有着诸多应用，
- 几何学中，存在多种用数据变换来表示三维空间旋转的方式
- 物理学中，经典力学中用于定量描述真实世界动能的物体旋转运动

【备注】：这里要注意区别一下旋转和方向这两个概念，方向表示物体某一时刻在参考坐标系下的方位，这个方位可以看做是由想象中的某个参考位置经过旋转运动（变换）而得到的，但并不表示这个物体的方向是由某个参考位置旋转而得来的。

三维旋转轴如果过物体的质心，那么叫做自旋转(spin)，如果旋转轴不过物体的质心，叫做(revolution)，旋转轴叫做极（pole）:
- 数学上认为旋转是刚性物体的移动，相较于平移（translation）不同的是至少保持某个点固定（a point fix）。
- 所有刚性物体的移动（movements）包括旋转（rotations）、平移（translations），和两者的组合三种形式。
- 绕坐标轴（x/y/z）的旋转称之为基础旋转，任何旋转都可以分解为基础旋转的组合。
- 在飞行动力学里面，基础旋转用yaw、pitch、roll表示（即Tait-Bryan Angles）（之前DICOM世界观·第一章 坐标系统博文有详细介绍）

上一篇博文中提到了表示旋转的几种常用方式，主要有欧拉角、余弦向量（矩阵）、四元数等。这里我们用MATLAB来演示一下各种方式之间的关系和换算方式。

1.5.1 欧拉角与余弦矩阵(EA 2 DCM)

EA=Euler Angle表示欧拉角，即之前提到的各种方式下绕各个轴的旋转角度，DCM=Direction Cosine Matrices表示余弦矩阵
根据1.4小节中的旋转（这里扩展到三维空间），用MATLAB先展示一下一个正方体分别绕Z旋转0.7854弧度，绕Y旋转0.1弧度、绕X轴旋转0弧度的效果图：

【备注】：选取的0.7854和0.1等比较奇特的弧度是因为直接从Matlab官方说明文档中截取的，可以通过help angle2dcm等指令打开查看，如下图所示：

z=[cos(0.7854),sin(0.7854),0;-sin(0.7854),cos(0.7854),0;0,0,1]z =    0.7071    0.7071         0   -0.7071    0.7071         0         0         0    1.0000y=[cos(0.1),0,-sin(0.1);0,1,0;sin(0.1),0,cos(0.1)]y =    0.9950         0    -0.0998         0    1.0000          0   0.0998          0     0.9950x=[1,0,0,;0,cos(0),sin(0);0,-sin(0),cos(0)]x =     1     0     0     0     1     0     0     0     1%即yaw=0.7854,pitch=0.1,roll=0

按照之前DICOM世界观·第一章 坐标系统博文介绍的内旋方式（即按照旋转后的各自坐标轴进行旋转）我们看一下上述yaw、pitch、roll旋转后的具体效果如何，即将上面效果图的三种旋转进行串联，效果图如下：

• 第一列是原始立方体，第二列是整体旋转后的示意图，第三列是旋转结果在XY平面的投影；第四列是旋转结果在XZ平面的投影；第五列是旋转结果在YZ平面的投影。
• 第一行是只绕着Z轴旋转的结果（可以从第一行第三列的XY投影的正方形可以看出，由于截图的原因使得正方形看起来像平行四边形）即旋转矩阵为z；第二行是在第一行基础上绕着Y轴旋转的结果（从第二行第三列投影图可以看出），即旋转矩阵为y∗z；第三行是在第二行基础上绕着X轴旋转的结果，即旋转矩阵为x∗y∗z。

即整体的上述旋转矩阵A=x∗y∗z（注意矩阵右乘的顺序），

%我们验证一下，变换矩阵A=x*y*zA =    0.7036    0.7036   -0.0998   -0.7071    0.7071         0    0.0706    0.0706    0.9950

【备注】：注意旋转矩阵A的每个列向量A=(a1,a2,a3)分别对应三个余弦向量，分别表示旋转后的新坐标轴X‘，Y′，Z‘与原始坐标轴X,Y,Z中各个轴夹角的余弦。
为了验证一下我们的结果，这里直接利用angle2dcm函数计算如下：

yaw = 0.7854;  pitch = 0.1; roll = 0;dcm = angle2dcm( yaw, pitch, roll )dcm =    0.7036    0.7036   -0.0998   -0.7071    0.7071         0    0.0706    0.0706    0.9950

两个结果是一致的，至此我们清楚了欧拉旋转角与余弦矩阵之间的关系和转换方式。

1.5.2 欧拉角与四元数(EA 2 Quat)

EA=Euler Angle表示欧拉角，即之前提到的各种方式下绕各个轴的旋转角度，Quat=Quaternion表示四元数。四元数表示空间旋转的优点在上一篇博文中介绍过了，主要是节省空间并能够避免万向锁。两者之间的转换关系可以简单的根据上述“欧拉角到余弦矩阵DCM”的转换来求解出来，参考资料给出两者的换算关系如下：

由四元数计算余弦矩阵

由余弦矩阵DCM反推四元数

下面借助MATLAB计算验证一下：

z=[cos(0.7854),sin(0.7854),0;-sin(0.7854),cos(0.7854),0;0,0,1]z =    0.7071    0.7071         0   -0.7071    0.7071         0         0         0    1.0000y=[cos(0.1),0,-sin(0.1);0,1,0;sin(0.1),0,cos(0.1)]y =    0.9950         0    -0.0998         0    1.0000         0   0.0998          0     0.9950x=[1,0,0,;0,cos(0),sin(0);0,-sin(0),cos(0)]x =     1     0     0     0     1     0     0     0     1

仿照上一小节，手动计算一下：
第一步，根据欧拉角计算DCM余弦矩阵

%根据欧拉角（yaw、pitch、roll）计算余弦矩阵DCMDCM = x*y*z    0.7036    0.7036   -0.0998   -0.7071    0.7071         0    0.0706    0.0706    0.9950

第二步，根据上文公式，手动计算四元数各个参数

%公式如下：%a = 0.5  * sqrt(1+R11+R22+R33)    (not stored)在NIFTI文件格式中通过a=sqrt(1-b*b-c*c-d*d)计算而来%b = 0.25 * (R32-R23) / a       => quatern_b%c = 0.25 * (R13-R31) / a       => quatern_c%d = 0.25 * (R21-R12) / a       => quatern_da=sqrt(1+0.7036+0.7071+0.9950)/2a = 0.9227b=0.25*(-0.0706)/0.9227b = -0.0191c=0.25*(0.0998+0.0706)/0.9227c = 0.0462

第三步，用matlab自带angle2quat计算

q=angle2quat(yaw,pitch,roll)q = 0.9227   -0.0191    0.0462    0.3822

由此可以看出两者转换结果一致。

总结：

利用matlab工具，通过示例的演示，详细的介绍了欧拉角、余弦向量（矩阵）、四元数表示旋转的具体含义，以及简单介绍了彼此之间的转换，如果读者想了解更详细的信息可参考matlab中这一类的工具，包括但不仅限于angle2dcm、dcm2angle、quat2dcm、dcm2quat、dcm2quat、quat2dcm，另外可以参考NIfTI官方文档对于qform、sform的介绍（【备注】：NIfTI是在fMRI领域应用很广泛的一种医学图像格式，后续也会详细介绍，敬请期待）。

作者：[e-mail]
时间：2017-06-04