一种基于应用的DICOM医学图像处理方法

摘要

本发明提供了一种基于应用的DICOM医学图像处理方法，该方法包括：步骤一，输入两幅待配准的DICOM图像，依次对其进行格式转换、图像去噪处理和图像分割，从分割后的图像中抽取出特征信息组成特征空间；步骤二，给定初始变换，采用线性插值方法对特征空间进行线性变换；步骤三，选择优化策略，在特征空间中寻找最优变换；步骤四，根据预定义相似度函数来测量两幅待配准的DICOM图像的相似度，如果当前变换满足预设相似度阈值，则输出该变换，否则重复步骤二和步骤三。本发明提出了改进的图像配准方法，有效提高了DICOM医学图像配准的速度和精度。

说明书

技术领域

本发明涉及医学图像处理，特别涉及一种基于应用的DICOM医学图像处理 方法。

背景技术

随着计算机科学的发展，医学影像信息学已成为国际上发展迅速的一个领 域。DICOM作为医学图像数据标准的使用使不同厂家，不同设备之间医学图像数 据传输，交换和存储成为可能。DICOM标准强大的兼容性使它越来越广泛地应用 到医学领域，逐步成为一种标准的医学图像数据格式。DICOM医学图像为临床诊 断医学、计算机辅助手术、病灶监控等提供了高质量的信息支持，但是，DICOM 图像也带来了一些问题，如时序差异、位置差异、获得设备差异等导致影像版 本信息的不一致性。怎样消除各种差异通过图像配准获取最大有效信息成为亟 待解决的问题。医学图像的配准关键在于界定医学图像之间的相似性。在临床 方面，这种方法常用于匹配一段时间内来自不同或同一患者的医学图像。如果 将普通图像配准直接应用于结构复杂的DICOM医学图像，则不仅计算量非常大， 效率低；而且容易产生误差导致结果不准确。

因此，针对相关技术中所存在的上述问题，目前尚未提出有效的解决方案。

发明内容

为解决上述现有技术所存在的问题，本发明提出了一种基于应用的DICOM 医学图像处理方法，包括：

步骤一，输入两幅待配准的DICOM图像，依次对其进行格式转换、图像去 噪处理和图像分割，从分割后的图像中抽取出特征信息组成特征空间；

步骤二，给定初始变换，采用线性插值方法对特征空间进行线性变换；

步骤三，选择优化策略，在特征空间中寻找最优变换；

步骤四，根据预定义相似度函数来测量两幅待配准的DICOM图像的相似度， 如果当前变换满足预设相似度阈值，则输出该变换，否则重复步骤二和步骤三。

优选地，所述线性插值进一步包括：

利用最近邻域点灰度值的加权平均值来表示该点灰度值，通过线性插值计 算直接修改直方图的项，将每对像素对联合直方图的贡献分散到联合直方图上 相邻的各个像素对，联合直方图上各点的值以小于1的小数增加，从而得到一 个平滑的配准函数。

优选地，所述图像去噪进一步包括：

对于每个种子像素点I(k，l，n)，选定最大搜索区域为M×M，其中k和l 代表像素坐标，n代表帧数；

种子像素点的8邻域像素表示为I(i,j，n)，i,j∈Ω，根据以下准则为每 个像素确定一个权值：

$|I(k,l,n)-I(i,j,n)|<\sqrt{2}{\sigma }_N$

其中的σ_N为图像噪声标准差的估计值；

通过加权平均得到该种子像素点的估计灰度值：

$I_e(k,l,n)=\frac{\underset{i,j\in \Omega }{\Sigma }p(i,j)I(i,j,n)}{\underset{i,j\in \Omega }{\Sigma }p(i,j)}$

其中$p(i,j)=\left(\begin{array}{cc}1& |I(k,l,n)-I(i,j,n)|<\sqrt{2}{\sigma }_N\\ 0& I(k,l,n)-I(i,j,n)|\ge \sqrt{2}{\sigma }_N\end{array}\right),$Ω表示8邻域像素坐标 空间；

然后在选定的搜索区域M×M内，以I_e(k，l，n)为种子像素点构建邻域： 将邻域内的像素组成自适应邻域，得到一个像素 点集合；

设置空域的高斯核函数为：

$g_s(i,j)=\exp (-\frac{D(I(k,l,n),I(i,j,n))}{2{\mathrm{bn}}^2})$

其中D代表种子像素点I(k,l,n)和自适应邻域像素点集合中任意一点I(i,j,n) 之间的欧氏距离，bn的值取最大搜索区域尺寸M的一半；

设置灰度域的高斯核函数为：

$g_g(i,j)=\exp (-\frac{{(I_e(k,l,n),I(i,j,n))}^2}{2{\mathrm{bn}}^2})$

其中σ_N仍为图像噪声标准差的估计值，通过空间域和灰度域的高斯核函 数，按下式进行滤波：

$\hat{I}(k,l,n)=\frac{\underset{i,j\in V}{\Sigma }\left[g_s(i,j)g_g(i,j)I(i,j,n)\right]}{\underset{i,j\in V}{\Sigma }\left[g_s(i,j)g_g(i,j)\right]}$

其中V表示自适应邻域像素点集合的坐标空间。

优选地，对于待处理的第n帧，建立一个图像缓冲器，用来暂存第(n-1)、 n、(n+1)三帧图像，首先计算以两个对应像素点为中心，大小为3×3的前后两 帧图像块之间灰度差的绝对值和SAD_prev和SAD_next：

${\mathrm{SAD}}_{\mathrm{prev}}=\underset{i,j\in \Omega }{\Sigma }|\hat{I}(i,j,n)-\hat{I}(i,j,n-1)|$

${\mathrm{SAD}}_{\mathrm{next}}=\underset{i,j\in \Omega }{\Sigma }|\hat{I}(i,j,n)-\hat{I}(i,j,n+1)|$

其中表示待处理的帧，它与其前后两帧图像都经过的空域和灰度域 滤波，Ω表示3×3图像块的像素坐标空间；

设置时域的高斯核函数为：

$g_t(\mathrm{SAD},T)=\exp \left(\frac{{\mathrm{SAD}}^2}{2T^2}\right)$

其中T是由噪声引起的误差决定的常量参数，最终的时域滤波按下式进行：

${\hat{I}}^{\prime }(k,l,n)=\frac{\hat{I}(k,l,n)+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)\hat{I}(k,l,n-1)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)\hat{I}(k,l,n+1)}{1+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)}.$

本发明相比现有技术，具有以下优点：

提出了改进的图像配准方法，有效提高了DICOM医学图像配准的速度和精 度。

附图说明

图1是根据本发明实施例的基于应用的DICOM医学图像处理方法的流程图。

具体实施方式

下文与图示本发明原理的附图一起提供对本发明一个或者多个实施例的详 细描述。结合这样的实施例描述本发明，但是本发明不限于任何实施例。本发 明的范围仅由权利要求书限定，并且本发明涵盖诸多替代、修改和等同物。在 下文描述中阐述诸多具体细节以便提供对本发明的透彻理解。出于示例的目的 而提供这些细节，并且无这些具体细节中的一些或者所有细节也可以根据权利 要求书实现本发明。

本发明提出一种DICOM医学图像配准方法。先对DICOM医学图像进行预处 理，然后测量待配准图像间的相似度，最后优化图像。图1是根据本发明实施 例的基于应用的DICOM医学图像处理方法流程图。首先抽取出图像的特征信息 组成特征空间；然后根据所抽取的特征空间确定出一个空间变换，使一幅图像 经过该变换后能够达到所定义的相似性测度。在确定变换的过程中若采取一定 的优化措施或搜索策略，则可使相似性测度更好、更快地达到最优值。

本发明寻找最佳空间和灰度变换函数，便于图像使用相似性测度进行匹配， 其过程就是寻找最优变换，其中该最优变换将图像B进行变换，以使其尽可能 多的包含图像A的信息。用互信息表示相似性测度对两图像进行精确配准，通 过对完整的医学图像寻找一个仿射变换即旋转参数和平移参数，该仿射变换对 应的两个图像的互信息最大，此时图像实现精确配准互信息用熵来表示：

对于概率分布函数为p(a)的随机变量集A，其熵H(A)定义如下：

H(A)＝-∑p(a)logp(a)a∈A

对两个图像，它们的离散的随机变量集分别是A和B，它们的边缘概率分布 函数分别是p(a)和p(b)，联合概率分布函数是P_AB(a，b)，则随机变量A和B的 联合熵定义如下：

H(A)＝-∑p(a)logp(a)logp(b)a∈Ab∈B

如果H(A/B)表示已知B是A的条件熵，那么H(A)与H(A/B)的差值，就代表 了在B中包含的A的信息，即互信息；因此两个随机变量集间的互信息定义为：

I(A，B)＝H(A)+H(B)-H(A，B)

＝H(A)-H(A/B)

＝H(B)-H(B/A)

在图像配准中，当两幅图像的空间位置达到完全一致时，其中一幅图像A 表达的关于另一幅图像B的信息，也就是对应像素灰度的互信息应为最大；用 联合概率分布和完全独立时的概率分布间的广义距离来估计互信息：

$I(A,B)=\mathrm{\Sigma p}(a,b)\log \frac{p(a,b)}{p\left(a\right)p\left(b\right)}$

对于待配准的两幅图像，可以认为它们是关于图像灰度的两个随机变量集 参考图像和浮动图像，a和b是两幅图像中相关的像素灰度值，a和b通过图像 之间的坐标变换相联系，对于离散的数字图像，联合概率分布p_AB(a，b)可以用 归一化的联合直方图表示，则有：

$I(A,B)=\underset{i,j}{\Sigma }{p}_{\mathrm{AB}}(i,j)\log \frac{p_{\mathrm{AB}}(i,j)}{p_A\left(i\right)p_B\left(j\right)}$

接下来寻找一个合适的线性变换，使得一幅试验图像经过此变换后和参考 图像的互信息最大。

本发明的配准方法总体步骤为：

1、输入两幅待配准的DICOM图像，即参考图像和试验图像，对其进行格式 转换；然后进行图像去噪和图像分割，从处理后的图像中抽取出特征信息组成 特征空间。

2、给定初始变换，使用线性插值方法进行线性变换。

3、选择优化策略，在特征空间中寻找最优变换。

4、测量参考图像和试验图像的相似度，如果当前最优变换满足预设相似度 阈值，则输出该变换，否则重复步骤2和3。

在图像预处理的图像格式转换过程中，首先在DIOCM文件中找到数据域的 起始位置，用pdata表示，根据DICOM的组标号，单位标号，获取数据的帧数 和数据长度。判断当前字节排序方式，如果排序方式是BIG_ENDIAN，则分别交 换组标号和单位标号的2个单元。若组标号为0x0002，单元标号是0x0010，则 获取UID，并根据UID判断数据的字节顺序；若组标号为0x0002，获取单元标 号，根据所获取的单元标号的值确定数据帧数，每帧中的行数和列数；若组标 号＝0x7fe0，则设置位图的文件头、信息头。判断是否需要设定颜色表，如果需 要，则获取单元标号，根据所获取的单元标号的值确定数据的存储长度，并获 取数据起始位置。设置颜色表并按照逆序方式存储pdata数据到位图文件，如 果不需要，则按照行逆序和RGB逆序的方式存储pdata数据到位图文件，最后 对位图文件进行显示。

DICOM文件可能封装了多幅图像，对于这样的DICOM文件，需要找到每幅图 像的偏移量，根据偏移量分次转换。在DICOM图像文件转换为位图文件只需把 DICOM中的像素数据输出到BMP文件即可，但在转换前需考虑窗宽和窗位。窗宽 是指图像数据显示的范围，窗位是指图像数据显示的中心值，两指标的值根据 不同的生理组织结构设定不同的值，当把DICOM文件的全部像素读入内存后， 就要考虑如何把DICOM显示到屏幕上，由于DICOM图像像素的灰阶可达到2¹⁶， 而系统的调色板只能显示2⁸种颜色，所以为了较逼真的显示医学图像的全部信 息，必须对图像数据进行窗宽、窗位的调整，把窗口区域图像数据线形的转换 到显示器的最大显示范围内，这样就可以动态调整窗宽和窗位，观察医学图像 的全部信息。可对图像的每个数据进行转换，把数据压缩到256个灰度等级。

其中：y为显示的位图灰度值，x为图像的数据，w为窗宽，c为窗位

$y\left(x\right)=\left(\begin{array}{cc}0& c-w/2>x\\ [x-(c-w/2)]255/w& c-w/2<x<c+w/2\\ 255& c+w/2<x\end{array}\right)$

由于DICOM图像和位图所规定的图像数据的存储顺序不同，所以在存放位 图数据时，其顺序是先从最后一行开始，依次向上，直到第一行。

位图文件有4部分数据项组成：位图文件头、位图信息头、颜色表和字节 阵列。初始化位图文件实现了对位图文件数据项的赋值和对颜色表的设置。设 置彩色表是将数据的灰度值映射为较少级数的*’+图像下的灰度值。位图文件 头主要包括文件类型，文件大小，字节偏移量等。每个数据项的值有DICOM格 式的行数、列数、帧数等确定。

要实现位图数据的写入，应首先确定DICOM文件中数据域的起始地址值， 该起始地址是通过DICOM数据预处理得到的，用pdata来记录该值。位图数据 的写入分两种情况，一种是DICOM格式下真彩色存储的情况，另一种是DICOM格 式下非真彩色存储的情况。两种情况在由DICOM格式转存为位图格式时都将采 取逆序方式存储原来数据域的值，不同点在于非真彩色存储时只需将DICOM数 据按行逆序存储即可转换为相应位图，而真彩色除了按行逆序存储外还需对RGB 值也逆序存储。将DICOM格式的图像转换成非真彩色存储时，每个数据单元的 值不会发生变化，但是存储位置都将发生变化，主要表现在逆序存储每一行的 数据值，即把DICOM数据中最后一行的值存储在位图非真彩色存储的第一行上， 第一行的数据存储在最后一行上。真彩色格式像素的存储比非真彩色要复杂得 多，除了采取行逆序存储外，还要对每个像素中RGB值逆序存储。

在图像预处理的图像去噪过程中，医学图像通常存在大量的斑点噪声，使 DICOM图像的质量明显下降，甚至掩盖了某些组织细节，这对图像特征的提取、 识别与分析带来了极大困难。通常采用滤波的方法来抑制噪声频率成分，从而 提高图像的信噪比，改善图像质量。本发明在现有技术的基础上，相继在空域、 灰度域、时域中进行滤波，有效抑制图像中的噪声，并保留边缘细节信息。提 高了算法的鲁棒性。有利于特征提取，识别和分析。具体的去噪过程为：

对于一个种子像素点，I(k，l，n)，选定最大搜索区域为M×M，n代表帧 数。先考虑种子像素点的8邻域像素，I(i,j，n)，并根据以下准则为每个像素 确定一个权值：$|I(k,l,n)-I(i,j,n)|<\sqrt{2}{\sigma }_N$

其中的σ_N为图像噪声标准差的估计值。然后，通过加权平均得到一个种子 像素点的估计灰度值：

$I_e(k,l,n)=\frac{\underset{i,j\in \Omega }{\Sigma }p(i,j)I(i,j,n)}{\underset{i,j\in \Omega }{\Sigma }p(i,j)}$

其中$p(i,j)=\left(\begin{array}{cc}1& |I(k,l,n)-I(i,j,n)|<\sqrt{2}{\sigma }_N\\ 0& I(k,l,n)-I(i,j,n)|\ge \sqrt{2}{\sigma }_N\end{array}\right),$Ω表示8邻域像素坐标空间。 这一步是为了减小种子像素点本身噪声较大时对自适应邻域构建的影响。

然后在选定的搜索区域M×M内，以I_e(k，l，n)为种子像素点构建邻域： 将满足上式的像素，将其纳入最终的自适应 邻域，从而得到一个像素点集合。在这一步中，阈值放大到了配合高斯 加权滤波，使自适应邻域的构建不会对这个阈值过分依赖。

设置空域的高斯核函数为：

$g_s(i,j)=\exp (-\frac{D(I(k,l,n),I(i,j,n))}{2{\mathrm{bn}}^2})$

其中D代表种子像素点I(k,l,n)和自适应邻域像素点集合中任意一点I(i,j,n) 之间的欧氏距离，bn的值取最大搜索区域尺寸M的一半。

灰度域的高斯核函数为：

$g_g(i,j)=\exp (-\frac{{(I_e(k,l,n),I(i,j,n))}^2}{2{\mathrm{bn}}^2})$

其中σ_N仍为图像噪声标准差的估计值。结合空间域和灰度域的高斯核函 数，最终的滤波按下式进行：

$\hat{I}(k,l,n)=\frac{\underset{i,j\in V}{\Sigma }\left[g_s(i,j)g_g(i,j)I(i,j,n)\right]}{\underset{i,j\in V}{\Sigma }\left[g_s(i,j)g_g(i,j)\right]}$

其中V表示自适应邻域像素点集合的坐标空间。

为了避免ghost和提高噪声抑制能力，对于待处理的第n帧，其前后两帧 与它的相关性最大。为此本发明建立一个图像缓冲器，用来暂存第(n-1)、n、(n+1) 三帧图像。首先，计算以两个对应像素点为中心，大小为3×3的图像块之间灰 度差的绝对值和SAD_prev和SAD_next：

${\mathrm{SAD}}_{\mathrm{prev}}=\underset{i,j\in \Omega }{\Sigma }|\hat{I}(i,j,n)-\hat{I}(i,j,n-1)|$

${\mathrm{SAD}}_{\mathrm{next}}=\underset{i,j\in \Omega }{\Sigma }|\hat{I}(i,j,n)-\hat{I}(i,j,n+1)|$

其中表示待处理的帧，它与其前后两帧图像都已经过了上一步的空 域和灰度域滤波，Ω表示3×3图像块的像素坐标空间。时域的高斯核函数为：

$g_t(\mathrm{SAD},T)=\exp \left(\frac{{\mathrm{SAD}}^2}{2T^2}\right)$

其中T是由噪声引起的误差决定的常量参数。最终的时域滤波按下式进行：

${\hat{I}}^{\prime }(k,l,n)=\frac{\hat{I}(k,l,n)+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)\hat{I}(k,l,n-1)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)\hat{I}(k,l,n+1)}{1+g_t({\mathrm{SAD}}_{\mathrm{prev}},T)+g_t({\mathrm{SAD}}_{\mathrm{next}},T)}.$

进一步地，在步骤2中，本发明采用线性变换和刚性配准。在医学图像配 准中，称刚体变换和仿射变换为刚性配准。当试验图像变换后坐标位置不在网 格点，通常采用插值技术对其灰度值进行估计。插值方法影响图像的平滑性、 优化性能和总体计算时间。为了提高优化效率和节省总体时间，本发明采用了 线性插值和为建立联合灰度直方输入DICOM图像而进行部分体积插值。利用最 近邻域点灰度值的加权平均值来表示该点灰度值。部分体积插值使用线性插值 权值，但并不通过邻域点加权来计算该点灰度值和更新一个简单直方图条目， 而是通过插值计算直接修改几个直方图的项，将每对像素对联合直方图的贡献 分散到联合直方图上与之相邻的各个像素对，联合直方图上各点的值以小于1 的小数增加，从而得到一个平滑的配准函数，有利于优化搜索。

步骤3中，配准方法的核心是一个含有能搜索到最小化测度函数的变换模 型参数的优化框架。在基于特征配准方法中，例如一系列点相关的点匹配，使 用直接优化技术，如最小范数的解决方案。在较为复杂的情况下，使用迭代优 化方法：

步骤A1、给定允许误差ε>0，选取初始点X₀(0，0，0)和一组线性无关的搜 索方向向量e⁽ⁱ⁾(i＝1，2，3)，其中e⁽¹⁾＝(1，0，0)，e⁽²⁾＝(0，1，0)，e⁽³⁾＝(0，0， 1)，即令e⁽ⁱ⁾分别为沿坐标轴的方向；设置迭代计数器K＝1；

步骤A2、从X₀出发沿方向e⁽ⁱ⁾(i＝1，2，3)依次进行一维搜索，得到

X⁽ⁱ⁾＝X_(i-1)+λ_i×e⁽ⁱ⁾

f(X_(i))＝min(f(X_(i-1)+λ_i×e⁽ⁱ⁾))

其中λ_i为预设常数权值，在完成了这3次一维搜索后，得到X₍₃₎；

步骤A3、检验是否已满足终止准则：若|f(X_(i))-f(X₍₀₎)|≤ε或者迭代次数 大于等于U，则停止迭代，退出；其中，ε为预先设定的阈值，U是预设的最大 迭代次数；

步骤A4、按照下式计算最速下降方向上的函数的变化量：

Δf＝max|f(X(i))-f(X(0))|(1<i<N)，并记最快下降方向为e_max；

步骤A5、取e＝X₍₃₎-X₍₀₎，Pcs＝2X₍₃₎-X₍₀₎，计算f(e)＝f(X₍₃₎-X₍₀₎)， f(Pcs)＝f(2X₍₃₎-X₍₀₎)，若f(e)≥f(x₍₀₎)或f(Pcs)<f(x₍₀₎)，且 2(f(x₍₀₎)-2f(x_(n))+f(e))[(f(x₍₀₎)-f(x_(n)))-Δf]²≥(f(x₍₀₎)-f(e))²Δf，则将X₍₃₎赋值给X₍₀₎作为新的起点，沿上面的一组旧方向向量e⁽ⁱ⁾(i＝1，2，3)重复步骤A2 至A5；若上述条件均不满足，则沿方向e＝X₍₃₎-X₍₀₎，以X₍₃₎作为起始点进行搜索 得到目标函数在此方向上的极小值点P_min；将原来的下降方向最大的e_max去掉， 而保留其它的2个方向，再加上方向e，仍得到3个方向：e⁽¹⁾，e⁽²⁾，e⁽³⁾，把此 时的P_min赋给X₍₀₎作为新起点，重复步骤A2至A5。

本发明采用的迭代优化方法进行配准参数的寻优，不必对目标函数求导， 具有二次收敛性，收敛速度快，可加快搜索最大值的速度，且便于编程实现。

综上所述，本发明提出了改进的图像配准方法，有效提高了DICOM医学图 像配准的速度和精度。

显然，本领域的技术人员应该理解，上述的本发明的各模块或各步骤可以 用通用的计算系统来实现，它们可以集中在单个的计算系统上，或者分布在多 个计算系统所组成的网络上，可选地，它们可以用计算系统可执行的程序代码 来实现，从而，可以将它们存储在存储系统中由计算系统来执行。这样，本发 明不限制于任何特定的硬件和软件结合。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本 发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围 的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围 之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或 者这种范围和边界的等同形式内的全部变化和修改例。