关于医学图像处理中钙化积分的计算

前言

最近在做一个医学图像处理的项目，项目的内容就是在拉直的血管上面找到钙化的部分，并且计算出钙化积分。

效果图如下所示：  

第一章   医学背景知识介绍

1.1 冠脉动脉钙化以及钙化积分定义

      冠状动脉钙化（CAC）的定性检测根据病变大小超过某一面积或体积阈值、密度超过规定的CT值阈值即可区分冠状动脉斑块内的病变是否为钙化。冠状动脉钙化积分（CACS）中传统的积分法AS定义面积阈值为1 mm²、CT值阈值为130 Hu。Detrano等使用8.16mm³的体积作为钙化灶大小的阈值。质量积分^[10]定义面积阈值为连续的3个像素、CT值阈值为130 Hu。Nelson等^[11]所使用的面积阈值为4个紧邻的像素，CT值阈值为130 Hu。另有研究者认为^[12]不应使用某一固定的CT值作为病灶的密度阈值，因为同一密度的物体在不同CT机上的CT值可能有差异，故推荐使用一个固定的含钙浓度值（100 mg/cc）作为阈值。然而Glodny等^[13]却尝试在CCTA图像中将冠状动脉上CT值大于600 Hu的像素分离出来作为钙化灶来研究，这种方法虽有偏颇，会明显低估钙化斑块负荷，但其结果却和常规CACS有高度线性相关性，甚至认为该方法可用于替代常规CACS以避免CACS扫描导致的电离辐射，这一观点尚待进一步的考证。

     Agatston等^[2]首次使用了CACS对冠状动脉钙化斑块进行了定量分析，称之为Agatston 积分（Agatston Score, AS），
AS= Σ （CT值 × 面积 × 权重系数）。
CT值越高，权重系数越大，具体如下：130～199 Hu为1，200～299 HU为2，300～399 Hu为3，等于或大于400 Hu为4。
1998年Callister等^[14]报道了CAC的容积积分（Volume Score，VS）的定量分析。VS等于所有钙化斑块的容积之和。VS积分法所计算的是钙化的体积，而与钙化灶的密度变化没有对应关系。
2005年Nelson等^[11]使用了质量积分（Mass Score，MS）进行CAC的定量测量。MS的检测相对复杂，该方法须使用羟磷灰石钙校准体模，此体模由数个含钙浓度不同的模块组成，且浓度已知。体模与受检者须同时进行扫描。通过测量各个模块的CT值便可确定CT值与含钙浓度的计算关系，即：

体模CT值 = 斜率 × 体模含钙浓度 + 截距。

测出各个模块的CT值，可计算出斜率和截距。然后测量冠状动脉钙化斑块的CT值，根据上述关系计算出钙化斑块的含钙浓度，即：

钙化斑块的含钙浓度 = (钙化斑块平均CT值－截距) ÷ 斜率。

最后通过软件计算出冠状动脉全部钙化斑块的质量积分，即：

MS = ∑（钙化灶的含钙浓度 × 容积）。

AS是较客观的CAC量化分析方法，早己受到了影像界的认可，临床应用广泛，一直是传统的钙化斑块定量分析参数，但此方法的重复性不高，使病人间的对照及随访观察结果可信度下降。VS的重复性与AS相比有所改善，可用于高危人群随访及冠心病人治疗后复查。VS主要缺点是由于受部分容积效应的影响，另一方面VS没有考虑钙化斑块的CT值，而CT值的高低可能代表了斑块的演变阶段，因此可能存在一定的片面性，而且此方法的重复性仍有不足。由于每位患者的含钙浓度均是在体模校准的基础上获得，不受设备、扫描协议或患者个体差异的影响，其结果具有可比性，因而MS有良好的可重复性和准确性。在扫描协议标准化后，多中心的体模研究^[8]使用不同生产厂商的CT机、不同机型的CT机进行CACS测量，结果为AS、VS与MS的变异系数分别为7.9%、4.0%、2.8%， MS具有最高的重复性和准确性，这与绝大部分研究结果一致。MS是国际心脏CT标准化协会物理工作组（the Physics Task Group ofthe International Consortium on Standardization in Cardiac CT）推荐的钙化定量分析方法^[12]。

1.2窗宽、窗位、CT值

      人体组织CT值的范围为-1000到+1000共2000个分度，人眼不能分辨这样微小灰度的差别，仅能分辨16个灰阶。为了提高组织结构细节的显示， 能分辨CT值差别小的两种组织，操作人员可以根据诊断需要调节图像的对比度和亮度，这种调节技术称为窗技术-窗宽、窗位的选择。

      窗宽是指显示图像时所选用的CT值范围，在此范围内的组织结构按其密度高低从白到黑分为16个等级（灰阶）。如窗宽为160Hu，则可分辨的CT值为160/16=10Hu，即两种组织CT值的差别在10Hu以上者即可分辨出来。 由此可见窗宽的宽窄直接影响图像的对比度；窄窗宽显示的CT值范围小，每级灰阶代表的CT值幅度小，因而对比度强，可分辨密度较接近的组织或结构，如检查脑组织宜选用窄的窗宽；反之，窗宽加宽则每级灰阶代表的

CT值幅度增大，对比度差，适于分辨密度差别大的结构如肺、骨质。窗位是指窗宽上、下限CT值的平均数。因为不同组织的CT值不同，欲观察其细微结构最好选择该组织的CT值为中心进行扫描,这个中心即为窗位。窗位的高低影响图像的亮度：窗位低图像亮度高呈白色；窗位高图像亮度低呈黑色，但在实际操作中尚须兼顾其它结

构选用适当的窗位。总之，如要获得较清晰且能满足诊断要求的CT图像，必须选用合适的窗宽、窗位，否则

图像不清楚，还往往难以达到诊断要求，降低了CT扫描的诊断效能。

正常人体组织的CT值：

类别	CT值(HU)
水	0±10
脑脊液	3-8
血浆	3-14
水肿	7-17
脑白质	25-32
脑灰质	30-40
血液	13-32
血块	64-84
肝脏	50-70
脾脏	50-65
胰腺	45-55
肾脏	40-50
肌肉	40-80
胆囊	10-30
脂肪	-20～-80
钙化	80-300
空气-200HU以上	骨骼+400以上

     窗宽，窗位，CT值之间的数量关系：通俗的说：窗位是窗宽的中位数数值，窗宽是范围，比如窗宽250，窗位50，CT值范围是-75~175HU，CT值的范围的中位数必须是50，而范围则是250。换算成公式的话如下：

 CT值的下限x = L-W/2，上限y = L+W/2；

其中，L代表窗位的数值，W代表窗宽数值。

   

 第二章  算法设计

2.1待续

 

第三章  部分代码

3.1vtk to OpenCV

/*****************************************************/
/*
@function name:vtk2ipl
@function：vtk格式图片转化到opencv格式
@parameters
srcDataVtk,sliceOrder:输入，srcDataVtk是vtkImageData格式体数据，sliceOrder是要抽取的体数据的slice编号
dstDataMat :输出,dstDataMat 是opencv格式图像
@result:
@Notes:
*/
/*****************************************************/


void calScore::vtk2ipl(vtkImageData *srcDataVtk,Mat &dstDataMat,int sliceOrder)
{
    vtkImageShiftScale* shiftScale = vtkImageShiftScale::New();
    shiftScale->SetInput(srcDataVtk);
    shiftScale->SetShift(-(800 - 0.5*1680));//800level,1680windows 这里
    shiftScale->SetScale(255.0/1680);
    shiftScale->SetOutputScalarTypeToUnsignedChar();
    shiftScale->ClampOverflowOn();
    shiftScale->Update();

    if (c_mode_XYZ =="XZ")//xz平面
    {
        cv::Mat image(dims[0], dims[2], CV_8UC1);

        for (int i=0; i<dims[0]; ++i) 
        {
            for (int j=0; j<dims[2]; ++j) 
            {
                image.at<unsigned char>(cv::Point(j,i)) = 
                    *static_cast<unsigned char*>(shiftScale->GetOutput()->GetScalarPointer(i,sliceOrder,j));
            }
        }
        imwrite("D:\cvimg.bmp",image);
        image.copyTo(dstDataMat);
    }

    if (c_mode_XYZ =="YZ")//yz平面
    {
        cv::Mat image(dims[1], dims[2], CV_8UC1);

        for (int i=0; i<dims[1]; ++i) 
        {
            for (int j=0; j<dims[2]; ++j) 
            {
                image.at<unsigned char>(cv::Point(j,i)) = 
                    *static_cast<unsigned char*>(shiftScale->GetOutput()->GetScalarPointer(sliceOrder,i,j));
            }
        }
        imwrite("D:\cvimg.bmp",image);
        image.copyTo(dstDataMat);
    }

    if (c_mode_XYZ =="XY")//xy平面
    {
        cv::Mat image(dims[0], dims[1], CV_8UC1);

        for (int i=0; i<dims[0]; ++i) 
        {
            for (int j=0; j<dims[1]; ++j) 
            {
                image.at<unsigned char>(cv::Point(j,i)) = 
                *static_cast<unsigned char*>(shiftScale->GetOutput()->GetScalarPointer(i,j,sliceOrder));
                
            }
        }
        imwrite("D:\cvimg.bmp",image);
        image.copyTo(dstDataMat);
    }
    
    //free
    shiftScale->Delete();
}

　3.2 OpenCV to VTK

/*****************************************************/
     /*
     @function name:ipl2vtk
     @function：opencv格式图片转化到vtk格式
     @parameters
     _src:输入，-src 是OpenCV数据
     _dest :输出,_dest 是vtk格式图像
     @result:
     @Notes:
     */
    /*****************************************************/
void calScore::ipl2vtk(IplImage* _src, vtkImageData* _dest)
{
    vtkImageImport *importer = vtkImageImport::New();
    if ( _dest )
    {
        importer->SetOutput( _dest );
    }
    importer->SetDataSpacing( 1, 1, 1 );
    importer->SetDataOrigin( 0, 0, 0 );
    importer->SetWholeExtent(  0, _src->width-1, 0, _src->height-1, 0, _src->nChannels-1 );
    importer->SetDataExtentToWholeExtent();
    importer->SetDataScalarTypeToUnsignedChar();
    importer->SetNumberOfScalarComponents( _src->nChannels );
    importer->SetImportVoidPointer( _src->imageData );
    importer->Update();

    //free
    importer->Delete();

}

3.3 抽取slice

/*****************************************************/
     /*
     @function name:extractSlice
     @function：抽取slice
     @parameters
     _srcvtk,sliceOrder:输入，_srcvtk是vtkImageData格式体数据，sliceOrder是要抽取的体数据的slice编号
     @result:返回抽取的slice数据，数据格式是vtkImageData
     @Notes:
     */
    /*****************************************************/
vtkImageData *calScore::extractSlice(vtkImageData *_srcvtk,int sliceOrder)
{
    _srcvtk = vtkImageData::New();
    _srcvtk->SetDimensions(imageData->GetDimensions());
    _srcvtk->SetScalarTypeToDouble();
    _srcvtk->SetNumberOfScalarComponents(4);
    _srcvtk->AllocateScalars();
    _srcvtk->DeepCopy(imageData);
    _srcvtk->Update();
    

    vtkExtractVOI*extractVOI = vtkExtractVOI::New();
    //extractVOI->SetInput(_srcvtk);
    extractVOI->SetInput(imageData);

    if (c_mode_XYZ =="XZ")//xz平面
    {
        extractVOI->SetVOI(0,dims[0],sliceOrder,sliceOrder,0,dims[2]);
    }

    if (c_mode_XYZ =="YZ")//yz平面
    {
        extractVOI->SetVOI(sliceOrder,sliceOrder,0,dims[1],0,dims[2]);
    }

    if (c_mode_XYZ =="XY")//xy平面
    {
        extractVOI->SetVOI(0,dims[0],0,dims[1],sliceOrder,sliceOrder);
    }
    
    extractVOI->Update();

    //vtkImageData *imgData2D = vtkImageData::New();
    //imgData2D->DeepCopy(extractVOI->GetOutput());
    
    //return imgData2D;
    return extractVOI->GetOutput();
}

第四章 总结

4.1待续

 

 

 

 

引用：1.1节http://www.dxy.cn/bbs/topic/23603201?keywords=%20%20%20%20

         1.2节 http://blog.sina.com.cn/s/blog_4c88d09a0100k8b1.html   http://www.dxy.cn/bbs/thread/20979764#20979764