CN103908276B - 多层容积ct灌注成像源影像的自动化减影处理方法 - Google Patents

Abstract

本发明公开了一种多层容积CT灌注成像源影像的自动化减影处理方法，该方法首先根据CT灌注影像序列计算CT成像设备的球管每转一周所扫描的层数n，并确定扫描时间中每一扫描时刻容积成像数据的n幅图像，获取每一时刻所有图像的容积CT值集合，然后确定CT减影的被减时刻，由被减时刻后的容积CT值集合减去被减时刻的容积CT值集合得到减影容积CT值集合，最后根据减影容积CT值集合与其对应文件元信息集合，获得减影DICOM文件。通过该方法实现容积CTP源影像的自动化减影处理与成像，提高了减影的效率与可靠性，为形成一站式的减影CTP软件奠定了基础。

Description

多层容积CT灌注成像源影像的自动化减影处理方法

技术领域

本发明涉及容积减影成像的计算机自动化处理技术领域，具体涉及一种多层容积CT灌注成像源影像的自动化减影处理方法。

背景技术

CT灌注成像（CTPerfusion，CTP）是在静脉注射对比剂的同时，对感兴趣区层面或容积层块进行连续动态CT扫描，以获得所选层面或容积层块内每一个像素的时间─密度曲线（TimeDensityCurve，TDC），并根据此曲线通过不同的数学模型转换和计算机伪彩处理得到局部组织的血流量（BloodFlow，BF）、血容量（BloodVolume，BV）、平均通过时间（MeanTransitTime，MTT）、达峰时间（TimetoPeak，TTP）和表面通透性（permeabilitysurfacearea，PS）等血流动力学参数和灌注伪彩图。

CTP可以从形态学、功能代谢诸方面提供疾病的诊断信息，但若检查部位组织在CT图像上表现出高密度影（如骨组织、血肿、钙化、成骨性病变等），就限制了CTP技术的应用。这是因为高密度影本身，往往掩盖了对比剂血流灌注前、后局部组织密度及微循环血流动力学的变化，使得CTP难以甚至无法提供诊病信息。近年有文献报道，使用减影CTP（SubtractionCTPerfusion，SCTP）技术从软件的图像后处理角度能够有效地克服因组织本身高密度影对CTP诊断的影响[1,2]。但SCTP技术没有整合，后处理时间过长，自动化程度低，限制了其临床推广应用，需要研发一站式全自动计算机分析系统。SCTP一站式全自动计算机分析系统需要解决两个关键问题，一是CTP源影像容积数据的全自动计算机减影，生成SCTP源影像；二是找到合适数学模型，根据层面或容积层块内图像减运算后的每一个像素的TDC，计算机自动转换和生成局部组织BF、BV、MTT、TTP和PS等血流动力学参数值和灌注图像。

图像的减运算（即减影技术）在医学上已有广泛的应用，典型技术莫过于数字减影血管造影（DigitalSubtractionAngiography，DSA）[3]，近些年来CT、MRI领域也有涉猎。Hanna等[4]使用减影技术对肌肉与骨骼肿瘤MR灌注影像序列的对比剂前、后的的图像进行了减影处理，发现血管恶性组织、脂肪和出血在减影后的图像上对比度得到了增强。Terai等[5]使用动态对比剂增强减影磁共振方法诊断睾丸扭转。李征宇等[6]使用MR动态增强时间减影灌注成像技术对周围型肺癌肺动脉血供定性与定量评估研究时在工作站上将强化前后的图像行减影处理。Huang和Chamberli[7]使用计算机化的自动方法对CT灌注影像序列的对比剂增强前后的图像进行减处理，在减影后的图像上发现了对比剂注入前和注入后的CT图像上肉眼观察不到的病灶。但这些CT、MRI应用与已报道的SCTP一样，大多基于影像工作站上的减影软件[1,2,4-6]对单层图像（未涉及容积层块图像）进行半自动的减影处理，操作繁琐，费时且容易出错。

发明内容

针对现有技术中存在的缺陷，本发明的目的在于提供一种多层容积CT灌注成像源影像的自动化减影处理方法，通过该方法实现容积CTP源影像的自动化减影处理与成像。

为实现上述目的，本发明采用的技术方案如下：

一种多层容积CT灌注成像源影像的自动化减影处理方法，包括以下步骤：

（1）进行容积层CT灌注成像扫描，获取容积成像数据，根据扫描时间和容积成像数据计算CT成像设备的球管每转一周所扫描的层数n；所述的容积成像数据包括扫描得到的CT灌注影像序列中的DICOM图像；

（2）确定扫描时间中每一扫描时刻容积成像数据的n幅图像,获取每一扫描时刻n幅图像的容积CT值集合；

（3）确定CT减影的被减时刻，根据被减时刻的容积CT值集合得到其它扫描时刻的减影容积CT值集合；

（4）根据减影容积CT值集合与其对应的文件元信息集合，获取容积成像数据的减影DICOM文件。

进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（1）中，计算CT成像设备的球管每转一周所扫描的层数n的具体方式为：

（1-1）计算CT灌注影像序列中所有图像的图像平面与原点的最小距离，并存储在一个一维数组中；

（1-2）删除所述一维数组中的重复的元素，数组中剩余元素的个数为CT成像设备的球管每转一周所扫描的层数n。

进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（2）中，获取每一时刻所有图像的容积CT值集合的具体方式为：

（2-1）根据所获取的CT灌注影像序列中DICOM图像的图像序列号InstanceNunber确定某一扫描时刻timei容积成像数据的n幅图像；

（2-2）根据n幅图的图像序列号InstanceNunber的顺序依次读入对应图像的像素数值，并转换成对应图像的CT值，得到n个一维CT值数组；

（2-3）将所得到的n个一维CT值数组按照其对应图像的图像序列号InstanceNunber串接得到timei时刻的容积CT集合。

进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（3）中，所述的其它扫描时刻的减影容积CT值集合包括被减时刻之前的扫描时刻的减影容积CT值集合和被减时刻之后的扫描时刻的减影容积CT值集合。

进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（3）中，获取被减时刻之后的扫描时刻的减影容积CT值集合的具体方式为：

（3-1）确定被减时刻subtractedi，获取subtractedi时刻的容积CT集合HUset0；

（3-2）获取被减时刻subtractedi之后的某一扫描时刻t的容积容积CT集合HUset；

（3-3）计算某一扫描时刻t的减影容积CT值集合deltaHUset，计算公式为:deltaHUset=HUset-HUset0。

进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（3）中，被减时刻subtractedi之前的所有扫描时刻的减影容积CT值集合等于subtractedi+1时刻的减影容积CT集合。

再进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（3）中，如果减影容积CT值集合中某一元素的数值小于0，将该元素的数值置为0。

更进一步，如上所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，步骤（4）中，获取容积成像数据的减影DICOM文件的具体方式为：

假设需要获取t时刻的容积成像数据中第j幅图像的DICOM文件，首先获取t时刻的减影容积CT值集合deltaHUset，对集合deltaHUset中的第j个元素deltaHUset[j]进行线性变换得到第j幅图像的像素值，根据转换得到的第j幅图像的像素值与t时刻容积成像数据中第j幅图像的DICOM文件元信息得到subtractedi时刻容积成像数据第j幅图像的DICOM文件；其中，1≤j≤n。

本发明的有益效果在于：本发明所述的方法实现了容积灌注源影像的减影CTP源影像的全自动生成，为形成一站式的减影CTP软件奠定了基础。该方法全自动减影，速度快，不会因为人为的因素误操作，而生成错误的图像序列，提高了减影的效率与可靠性。

附图说明

图1为本发明一种多层容积CT灌注成像源影像的自动化减影处理方法的流程图；

图2为具体实施方式中计算某一时刻的减影容积CT值集合的流程图。

具体实施方式

下面结合说明书附图与具体实施方式对本发明做进一步的详细说明。

图1示出了本发明一种多层容积CT灌注成像源影像的自动化减影处理方法的流程图，该方法主要包括以下步骤：

步骤S11：根据CT灌注影像序列计算CT成像设备的球管每转一周所扫描的层数；

进行容积层CT灌注成像扫描，获取容积成像数据，根据扫描时间timenumber和扫描得到的容积成像数据计算CT成像设备的球管每转一周所扫描的层数n；所述的容积成像数据包括扫描得到的CT灌注影像序列中的DICOM图像。

CT灌注影像序列中的图像的存储格式为DICOM文件格式，DICOM文件一般由文件元信息（FileMetaInformation）和DICOM数据集合（DICOMDataSet）组成。数据元素（DataElement）是DICOM格式文件最基本的组成单元。它由4个部分组成：标签（Tag）、数据描述VR（ValueRepresentation）、数据长度VL（ValueLengh）以及数据域（Value）。使用数据元素标签Tag为(0020,0032)和(0020,0037)的数值就可以计算出图像平面与原点间的最小距离，用以表示图像位置。在CT灌注影像序列中，两幅图像的图像平面与原点间的最小距离相同，就认为这两幅图像处于同一个层面。计算CT灌注影像序列中所有图像的图像平面与原点间的最小距离，并存储在1个一维数组中。剔除这个数组中重复的元素，此时剩下元素的个数n就是CT成像设备的球管每转一周获得的层数。而这个影像序列的数据扫描时间timenumber是影像序列的总共拥有的图像数除以这个层数n。

步骤S12：根据每一扫描时刻容积成像数据的图像获取每一时刻所有图像的容积CT值集合；

确定扫描时间中每一扫描时刻容积成像数据的n幅图像,获取每一扫描时刻所得到的n幅图像的容积CT值集合。若多幅DICOM图像同属于一个医学图像序列，那么它们一般有一个唯一标识的SeriesInstanceUID(0020,000E)（图像序列的唯一标识符）和StudyInstanceUID(0020,000D)（研究的唯一标识符），可以根据每幅DICOM图像的图像序列号InstanceNumber(0020,0013)（标识图象的号码）确定某一扫描时刻容积成像数据是由哪几幅图像组成的。例如，进行8层CT灌注成像扫描时，即CT成像设备的球管每转一周同时扫描8层，假设数据采集（扫描时间）30s产生240幅图像，则每一扫描时刻的容积成像数据由8幅图像组成，而扫描得到的每一层的图像也可以确定，设获取到的240幅DICOM的文件名即图像序列号为IM1，IM2，……，IM239，IM240，因为一共是8层，每一层30幅图像（文件），所以每一层包含的图像（文件）情况如下：

第1层：IM1,IM9,……,IM233；

第2层：IM2,IM10,……,IM234；

第3层：IM3,IM11,……,IM235；

第4层：IM4,IM12,……,IM236；

第5层：IM5,IM13,……,IM237；

第6层：IM6,IM14,……,IM238；

第7层：IM7,IM15,……,IM239；

第8层：IM8,IM16,……,IM240；

共计240幅图像（文件）。

由步骤S11可知，CT成像设备的球管每转一周所扫描的层数n，因此每一扫描时刻的容积成像数据都有n幅图像，如第1秒获得的容积图像包括第1幅图像到第n幅图像，第1秒获得的容积图像包括第n+1幅图像到第2n幅图像。第i秒获得的容积图像包括第i*n+1幅图像到第(i+1)*n幅图像，这样，第i秒获得的容积的CT值集合是依次读取第i*n+1幅图像到第(i+1)*n幅图像的CT值，形成n个1维数组，并按图像对应的InstanceNumber从小到大的顺序串接成一个一维数组。计算容积CT值集合的具体方式如下：

例如，第timei时刻容积成像数据由n幅图像组成，其中timei≥1且timei≤timenumber。这n幅图像的InstanceNumber形成一个等差数列。这个等差数列的起始值是(timei-1)×n+1，而公差是1，一共有n项。其中n是CT成像设备的球管每转一周获得的层数。确定了某扫描时刻容积成像数据是由哪几幅图像组成后，就可以根据这些图像InstanceNumber从小到大的顺序读入对应图像的像素数值，再变换成CT值。每读入一幅图像就得到一个一维CT值数组，这样就可以得到n个一维CT值数组。把这n个一维CT值数组按照对应图像的InstanceNumber从小到大的顺序串接成一个一维数组。这个一维数组存储的就是第timei时刻容积CT值集合。

步骤S13：计算CT减影容积CT值集合；

确定CT减影的被减时刻，根据被减时刻的容积CT值集合得到其它扫描时刻的减影容积CT值集合。假设被减时刻记为subtractedi，被减时刻subtractedi之后的某一扫描时刻记为t，计算被减时刻之外的其它扫描时刻的减影容积CT值集合包括被减时刻之前的扫描时刻的减影容积CT值集合和被减时刻之后的扫描时刻的减影容积CT值集合。本实施方式中计算被减时刻之后某一扫描时刻t的减影容积CT值集合的具体方式如图2所示：

首先确定被减时刻subtractedi，然后使用步骤S12中计算容积CT值集合的方法获得被减时刻的CT值集合HUset0，再使用步骤S12中计算出的第t时刻的容积CT值集合HUset，然后用HUset减去HUset0得到t时刻的减影容积CT值集合deltaHUset，即被减时刻之后的t时刻的减影容积CT值集合deltaHUset＝HUset-HUset0。按此方法，可以计算出大于subtractedi的所有扫描时刻（被减时刻之后）的减影容积CT值集合。而小于等于被减时刻subtractedi的所有扫描时刻的减影容积CT值集合等于第subtractedi+1时刻的减影容积CT值集合。通过上述计算即可得到所有扫描时刻的减影容积CT值集合。由于deltaHUset、HUset和HUset0都是一维数组，那么：

deltaHUset[i]＝HUset[i]-HUset0[i]

其中，deltaHUset[i]、HUset[i]和HUset0[i]分别表示deltaHUset、HUset和HUset0下标为i的元素的数值。为了便于计算，deltaHUset中小于0的数值都置为0。

步骤S14：根据减影容积CT值集合与其对应的文件元信息集合，获取容积成像数据的减影DICOM文件。

根据减影容积CT值集合与其对应的文件元信息集合，获取容积成像数据的减影DICOM文件，具体方式如下：

假设某时刻t的减影容积CT值集合是deltaHUset，那么，该时刻对应的原始容积数据第j幅图像贡献的减影CT值数据是一个一维数组deltaHUset[j]，其中1≤j≤n。而deltaHUset[j]等于deltaHUset中从下标(j-1)×height×width到j×height×width之间的数值，其中，height和width分别是影像序列中任意一幅图像的高和宽，参见文献《基于DICOM标准的医学图像的显示方法研究》。再根据这个第j幅图像的DICOM文件元信息以及对deltaHUset[j]进行线性变换，就可以把deltaHUset[j]存储为一幅DICOM文件。这样就可以得到所有时刻容积成像数据对应的减影DICOM文件。

其中，对是deltaHUset[j]进行线性变换为现有技术，具体变化方法可参见文献《基于DICOM标准的医学图像的显示方法研究》。DICOM文件格式一般包括两部分：文件元信息（FileMetaInformation）和DICOM数据集合（DICOMDataSet）。文件元信息包含DICOM数据集合的识别信息，可分为3部分：文件前言（FilePreamble）、DICOM前缀（DICOMPrefix）和文件元要素（FileMetaElements），DICOM文件元信息可以通过Gdcm、DCMTK、CTN、ITK、VTK开源免费下载的DICOM开发工具包获取，也可以使用Matlab和MITK获取。将DICOM文件元信息与线性变换后的deltaHUset[j]结合得到DICOM文件也是现有技术，可以使用Gdcm、DCMTK、CTN、ITK、VTK等开源免费下载的DICOM开发工具包可以得到，也可以使用Matlab和MITK得到。例如，使用Matlab的dicomwrite函数即可办到。把DICOM文件元信息和对deltaHUset[j]线性变换的结果作为dicomwrite函数的输入参数，就可以得到DICOM文件。

本发明所述的减影处理方法与现有的半自动手工技术相比，提供了一种在PC机上对容积CT灌注（CTP）源影像全自动生成减影CTP源影像的方法，全自动减影，速度快，减影过程按照指定的程序进行，不会因为人为的因素误操作，而生成错误的图像序列。因此，减影的可靠性得到了保证。促进SCTP的临床应用以及理论研究，填补国产医学影像软件领域的一些空白点。

下面结合具体的实施例对本发明所述的方法进行进一步说明。

实施例

本实施例中CT扫描CT机进行2层CT灌注成像扫描，即球管每转一周同时扫描2层，数据采集4s产生8幅图像。每幅图像的存储格式为DICOM文件格式。8幅图像分别记为IM1、IM2…IM8，则第1秒得到的容积成像数据是IM1和IM2，第2秒得到的容积成像数据是IM3和IM4，第3秒得到的容积成像数据是IM5和IM6，第4秒得到的容积成像数据是IM7和IM8。其中，IM1、IM3、IM5、IM7是第一层随时间变化得到的数据信息，IM2、IM4、IM6、IM8是第二层随时间变化得到的数据信息。

如果令第1秒是被减时刻，那么第1秒扫描得到的两幅图像为IM1和IM2，根据2幅图的InstanceNunber的顺序依次读入IM1和IM2对应的像素数值，并转换成相应的CT值，得到2个一维CT值数组；读取IM1得到以为CT数组记为CTHU1，读取IM2得到以为CT数组记为CTHU2，将CTHU1和CTHU2串接成更大的一维数组VolumeCTHU1得到第1秒的容积CT集合。假设CTHU1有n1个元素，CTHU2有n2个元素，那么VolumeCTHU1有n1+n2个元素。VolumeCTHU1表示第1秒容积成像数据的CT值集合。同理，可以得到第2秒至第4秒的CT值集合VolumeCTHU2、VolumeCTHU3、VolumeCTHU4。在线性代数中VolumeCTHU1、VolumeCTHU2、VolumeCTHU3和VolumeCTHU4可以看做向量，用向量VolumeCTHU2减去VolumeCTHU1得到deltaVolumeCTHU2，用向量VolumeCTHU3减去VolumeCTHU1得到deltaVolumeCTHU3，用向量VolumeCTHU4减去VolumeCTHU1得到deltaVolumeCTHU4。对于deltaVolumeCTHU2前n1个元素可以组成向量deltaCTHU3，而后n2个元素可以组成向量deltaCTHU4。把deltaCTHU3和deltaCTHU4进行线性变换得到deltaX3和deltaX4。在得到VolumeCTHU2的同时，通过DICOM开发工具包可以得到IM3的文件元信息metadata3和IM4的文件元信息metadata4。使用DICOM开发工具包（例如Matlab的dicomwrite函数）把IM3的文件元信息metadata3和deltaX3作为输入参数可以得到DICOM文件。使用DICOM开发工具包（例如Matlab的dicomwrite函数）把IM4的文件元信息metadata4和deltaX4作为输入参数可以得到DICOM文件。同理，可以得到2层8幅CT灌注影像序列对应的减影CT灌注影像序列文件。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (8)

1.一种多层容积CT灌注成像源影像的自动化减影处理方法，包括以下步骤：

(1)进行容积层CT灌注成像扫描，获取容积成像数据，根据扫描时间和容积成像数据计算CT成像设备的球管每转一周所扫描的层数n；所述的容积成像数据包括扫描得到的CT灌注影像序列中的DICOM图像；

(2)确定扫描时间中每一扫描时刻容积成像数据的n幅图像,获取每一扫描时刻n幅图像的容积CT值集合；

(3)确定CT减影的被减时刻，根据被减时刻的容积CT值集合得到其它扫描时刻的减影容积CT值集合；

(4)根据减影容积CT值集合与其对应的文件元信息集合，获取容积成像数据的减影DICOM文件。

2.如权利要求1所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(1)中，计算CT成像设备的球管每转一周所扫描的层数n的具体方式为：

(1-1)计算CT灌注影像序列中所有图像的图像平面与原点的最小距离，并存储在一个一维数组中；

(1-2)删除所述一维数组中的重复的元素，数组中剩余元素的个数为CT成像设备的球管每转一周所扫描的层数n。

3.如权利要求2所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(2)中，获取每一时刻所有图像的容积CT值集合的具体方式为：

(2-1)根据所获取的CT灌注影像序列中DICOM图像的图像序列号InstanceNunber确定某一扫描时刻timei容积成像数据的n幅图像；

(2-2)根据n幅图的图像序列号InstanceNunber的顺序依次读入对应图像的像素数值，并转换成对应图像的CT值，得到n个一维CT值数组；

(2-3)将所得到的n个一维CT值数组按照其对应图像的InstanceNunber串接得到timei时刻的容积CT集合。

4.如权利要求3所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(3)中，所述的其它扫描时刻的减影容积CT值集合包括被减时刻之前的扫描时刻的减影容积CT值集合和被减时刻之后的扫描时刻的减影容积CT值集合。

5.如权利要求4所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(3)中，获取被减时刻之后的扫描时刻的减影容积CT值集合的具体方式为：

(3-1)确定被减时刻subtractedi，获取subtractedi时刻的容积CT集合HUset0；

(3-2)获取被减时刻subtractedi之后的某一扫描时刻t的容积CT集合HUset；

(3-3)计算某一扫描时刻t的减影容积CT值集合deltaHUset，计算公式为:deltaHUset＝HUset-HUset0。

6.如权利要求5所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(3)中，被减时刻subtractedi之前的所有扫描时刻的减影容积CT值集合等于subtractedi+1时刻的减影容积CT集合。

7.如权利要求6所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(3)中，如果被减时刻之后的扫描时刻的减影容积CT值集合中某一元素的数值小于0，将该元素的数值置为0。

8.如权利要求7所述的一种多层容积CT灌注成像源影像的自动化减影处理方法，其特征在于：步骤(4)中，获取容积成像数据的减影DICOM文件的具体方式为：

假设需要获取t时刻的容积成像数据中第j幅图像的DICOM文件，首先获取t时刻的减影容积CT值集合deltaHUset，对集合deltaHUset中的第j个元素deltaHUset[j]进行线性变换得到第j幅图像的像素值，根据转换得到的第j幅图像的像素值与t时刻容积成像数据中第j幅图像的DICOM文件元信息得到subtractedi时刻容积成像数据第j幅图像的DICOM文件；其中，1≤j≤n。