三维血管中轴线特征描述

原理

• 1 凸性：
  (1) 基于观察，对于血管上的一点，以该点为中心，在与血管方向垂直的一维横截线上的灰度分布，符合抛物线特征。我们使用抛物线方程的二次项系数作为描述血管点的特征量，并将其称为该点在该方向上的凸性。
  (2) 由于血管具有方向性，所以对于某一点，需要考虑横街线的方向。所以在计算所有方向上的凸性值后，选取最大的N个方向的乘积，作为凸性值。
  (3) 由于血管粗细不同，所以对于某一点，需要考虑横截线的半径。所以在计算某一半径范围内的所有凸性值后，选取最大的凸性值，作为该点的凸性值。
  (4) 设半径范围为[Rmin,Rmax], 而方向为θi,N为选取计算最终凸性特征值的方向数。
  则最终的凸性特征值为：
  
  Convexity=r∗∏i=1NMax{Cr,r∈[Rmin,Rmax]}(1)
  
  注意：对于血管背景点灰度比血管高的情况下(如二维)，我们计算的是凸性的最大值，对于血管背景点灰度比血管低的情况下(如三维)，我们计算的是凸性的最小值（因为此时二次项系数是负的）
• 2 对称性
  (1) 对称性基于凸性中的横截线的概念来计算，在凸性的计算中，我们可以得到凸性最大值的半径，利用该半径，计算以该点为中心的两侧中心对称的点的灰度差值来计算对称性。对称差值越小，则说明对称性越好。
  (2) 设在[Rmin,Rmax]中，在半径Rmc处取得最大的凸性值。
  
  Symmetry=R∑R∈[0,Rmc]|L+(p,R,θ)−L−(p,R,θ)|(2)

实现

• 1公式推导
  对于图像上的一点 p , 半径为 r 的横截线，根据其灰度值，求得其二次拟合的二次项系数。
  (1) 用到方法：a. 最小二乘法 b. 线性方程组的求解 c.矩阵的LU分解
  即对于二维数组[x,y]进行二次线性拟合。x=[012,..2∗r],y为对应点上的灰度值。
  对于拟合方程
  
  P(xi)=a2x2i+a1xi+a0(3)
  
  使用最小二乘法计算系数a0,a1,a2。即:
  
  Q=∑i=1me2i=∑i=1m(Yi−P(xi))2,m=2r(4)
  
  求Q的最小值，即求多元函数的最小值。利用求偏导的方法进行计算，最终可以转化为线性方程组的求解，即：
  
  ⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢m+1∑i=0mxi∑i=0mx2i∑i=0mxi∑i=0mx2i∑i=0mx3i∑i=0mx2i∑i=0mx3i∑i=0mx4i⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢a0a1a2⎤⎦⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∑i=0myi∑i=0mxiyi∑i=0mx2iyi⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥(5)
  
  对系数矩阵进行LU分解，对于方程组
  
  Ax=b(6)
  
  首先进行对矩阵A进行LU分解，设A中的元素为aki ,
  
  ⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪lik=ukiukkuki=aki−∑q=1k−1uqkuqiuqq(8)
  
  计算得到LU矩阵后，由公式(8), 计算得到X:
  
  {LY=bUX=Y(7)
• 2 代码实现
  step1 对左矩阵进行LU分解。由于对于不同的点，我们的横截线的 X，是相同的，即x=[012,..2∗r]。所以可以事先用一个数组(pLM)计算好每个半径的LU分解好的矩阵。

//计算pLM（经过LU分解的系数矩阵）int i=0;int j=0;//LM means Left Matrix, it's the A in Eq.(6)//[fs,f]为半径范围const int LM_SIZE = 3;//pLM为每一个半径建立一个A系数矩阵double  *pLM = new double[LM_SIZE * LM_SIZE * (f - fs + 1)]();//nfSize为横截线的长度，fm为临时数组，用来事先计算好Eq.(5)中的左矩阵中的相关数据int nfSize  = 2 * f + 1;double *fm  = new double[nfSize * LM_SIZE];for(i = 0;i < LM_SIZE;i++){    for(j = 0;j < nfSize;j++)    {        fm[j + i * nfSize] = pow(double(j + 1), double(i));    }}//计算pLM，即每个半径的A系数矩阵int lf=0;for(lf = fs;lf <= f;lf++)  //遍历[fs,f]的每个半径长度{    for(i = 0;i < LM_SIZE;i++) //3*3系数矩阵    {        for(j = 0;j < LM_SIZE;j++) //根据Eq(5)计算        {            for(k = 0;k < (2 * lf + 1);k++)//系数矩阵中每一个元素都是累加和                pLM[j + i * LM_SIZE + (lf - fs) * LM_SIZE * LM_SIZE] += fm[k + i * nfSize] * fm[k + j * nfSize];        }    }//for i    LUDecomLeftM33(pLM + (lf - fs) * LM_SIZE * LM_SIZE);}delete []fm;

//计算二次多项式拟合的左矩阵的LU分解矩阵//LUDecomLeftM33 意思是LU分解3*3的左矩阵，分解后的LU还存储在原3*3矩阵中void LUDecomLeftM33(double p[3 * 3]){    int di, dj;    for(di = 0;di < 3;di++) {        for(dj = di;dj < 3;dj++)            p[dj + di * 3] = GetUkj(p, di, dj);        for(dj = di + 1;dj < 3;dj++)            p[di + dj * 3] = GetLik(p, dj, di);    }}//计算Ukj，利用Eq(8)inline double GetUki(double p[3 * 3], int k, int i){    int di;    double ftmp = 0;    for(di = 0;di <= k - 1;di++)        ftmp += p[di + k * 3] * p[i + di * 3];    ftmp = p[j + k * 3] - ftmp;    return ftmp;}//计算Lik，利用Eq(8)inline double GetLik(double p[3 * 3], int i, int k) {    return (p[i + k * 3] / p[k + k * 3]);}

step2 根据Eq(7)计算 x .

#define GETLM(i, j) (LeftM[(j) + (i) * LM_SIZE]) //define GETLM(i, j)是用来便于根据i，j来拿取LeftM中的值的//Y[3]为Eq(7)中的b, LeftM为L矩阵和U矩阵，Rx[2]即为x//GetPn1Coef的意思就是活得拟合多项式的二次项系数，这里把系数存在了Rx中void GetPn1Coef(double Y[3], const double LeftM[3 * 3], double RX[2]){    Y[1] = Y[1] - GETLM(1, 0) * Y[0];    Y[2] = Y[2] - GETLM(2, 0) * Y[0] - GETLM(2, 1) * Y[1];    RX[0] = Y[2] / GETLM(2, 2);    RX[1] = (Y[1] - RX[0] * GETLM(1, 2)) / GETLM(1, 1);}

step3: 利用step1和step2的方法计算出每一个半径的凸性值，从而得到最大的凸性值，作为最终的凸性值，再根据此时的半径，来计算出对称性

/**Function:MinimumPath_3DVascular*Description:compute the Convexity and Symmetry*Input (16):*   pInVol: origin 3D data*   bp: the point which needs to compute the cost, or the center point of the 1-D profile*   pLM: Left Matrix which has LU decomposed*   height: data height*   OMSize: width*height*   f: max radius*   fs: min radius*   Val_Thrd: Intensity Range*Output(8):*   pBalWeight: Convexity*   pSym:Symmetry*Return: void*/bool GetBalWeightij(const double    *pInVol,                    const CMyPoint  &bp,                    const double    *pLM,                    double          *pBalWeight,                    double          *pSym,                    int              height,                    int              OMSize,                    int              f,                    int              fs,                    const double    *Val_Thrd                    ) {    static double RY[LM_SIZE];//矩阵的右边    static double RX[2];//系数矩阵    static double pTmpRX0[40];//用来记录不同半径长度对应的凸性值    static double pCache[81];//用来存放横截线上的点对应的灰度值    int i, j, k, l;    int bwcount = 0;//用来为方向计数    int lf, fr, fb;    double dtmp, dlim, dsum, dWeight, dlfVal, drfVal;    //在26个方向，然后找到一个最大（或者最小的a系数的）fr    //fb记录具有最大凸性值的半径    for(k = -1;k <= 1;k++)    {        for(i = -1;i <= 1;i++)         {            for(j = -1;j <= 1;j++)            {                dlim = MatInf;//存放最小值                fb = f;//初始为最大半径                if(abs(k) + abs(i) + abs(j) == 0) //26个邻域点，13个方向                    return true;//故遍历到中心点时，返回结束                for(lf = fs;lf <= f;lf+=F_INC) //F_INC为步长，目的是减少计算量                {                    if((!InInterval(pInVol[bp.i + lf * i + (bp.j + lf * j) * height + (bp.k + lf * k) * OMSize], Val_Thrd)  ||                        !InInterval(pInVol[bp.i - lf * i + (bp.j - lf * j) * height + (bp.k - lf * k) * OMSize], Val_Thrd)) &&                        fb == f                        ) //背景跟目标的边界也有可能被认为是血管，去除这个边界                    {                        fb = lf;                    }                    //读取1D横截线灰度值                    //0，lf，2*lf,lf的位置是中心点                    pCache[lf] = pInVol[bp.i + bp.j * height + bp.k * OMSize];//中心点                    //若中心点的灰度超出范围，定义为边界值                    if(pCache[lf] < Val_Thrd[0])                        pCache[lf] = Val_Thrd[0];                    if(pCache[lf] > Val_Thrd[1])                        pCache[lf] = Val_Thrd[1];                    //读取灰度值                         dlfVal = pCache[lf];                    drfVal = pCache[lf];                    for(l = 1;l <= lf;l++)                    {                        dtmp = pInVol[bp.i + l * i + (bp.j + l * j) * height + (bp.k + l * k) * OMSize];                        if(InInterval(dtmp, Val_Thrd))                            dlfVal = dtmp;                        pCache[lf + l] = dlfVal;                        dtmp = pInVol[bp.i - l * i + (bp.j - l * j) * height + (bp.k - l * k) * OMSize];                        if(InInterval(dtmp, Val_Thrd))                            drfVal = dtmp;                        pCache[lf - l] = drfVal;                    }                    for(l = 0;l < LM_SIZE;l++)                        RY[l] = 0.0;                    //计算Eq(5)中的右矩阵                    for(l = 0;l <= 2 * lf;l++)                    {                        dtmp = pCache[l];                        RY[0] += dtmp;                        RY[1] += dtmp * (l + 1);                        RY[2] += dtmp * (l + 1) * (l + 1);                    }                    //根据右矩阵，系数矩阵(LU分解后的)，计算得出X,即二次拟合函数的系数                    GetPn1Coef(RY, (pLM + (lf - fs) * LM_SIZE * LM_SIZE), RX);                    //二次项系数存在RX[0]                    dtmp = RX[0];                    pTmpRX0[lf] = dtmp;                    if((dtmp*lf) < dlim)//找到最小的dtmp时的fr                    {                        dlim = dtmp*lf;                        fr = lf;                    }                }                //取最大凸性值时的半径和其相邻的半径的凸性值做平均，以平均值作为最终的半径                if(fr + 2 > f)                {                    dWeight = pTmpRX0[fr] * fr + pTmpRX0[fr - 2] * (fr - 2);                }                else                {                    dWeight = pTmpRX0[fr] * fr + pTmpRX0[fr + 2] * (fr + 2);                }                pBalWeight[bwcount] = -dWeight / 2;//凸性                //计算对称性，严格来讲，此时计算的是非对称性。在后面计算时，需要取倒数                dsum = 0.0;                for(l = 1;l <= fr;l++)                 {                    dtmp = pCache[f - l] - pCache[f + l];                    dsum += abs(dtmp);                }                pSym[bwcount] = dsum / (fr + 1);                bwcount++;            }        }    }    return false;}