基于形状的三维医学图像插值(二)

上篇存在的问题

这周和小组同学一起讨论，分析存在问题主要是由于以下两个原因

• 间隔1度取一个方向，再在这个方向上做插值，这个1度的间隔本身就漏掉一些点，并且会存在在边界处，越来越宽的情况
• 沿着当前方向计算时取出的点的坐标是浮点数，直接四舍侮五入造成的误差

解决方案

根据以上问题，首先解决需要解决的是对于待插值层的所有点都遍历进行插值，从而保证每个点都被计算到，避免出现空洞的情况。然后再根据待插值点与当前层中心组成向量的方向，在上下两层中沿着该方向搜索根据距离相似度找到匹配点。

实现

根据以上思路，先求待插值点和中心点所在向量与单位正方向的夹角,公式如下

向量余弦公式

反余弦函数

通过以上两步都存在判断正负的问题，容易造成较大的计算误差，考虑到求方向夹角主要是为了让上下两层都在此方向上求夹角的正余弦，因此转换思维由当前向量和单位正方向向量求正余弦

由以上公式可以看出，正余弦的正负取决于向量的正负，且单位正方向向量都可省略，解决由余弦求夹角，再由夹角转为正余弦造成的符号误差

代码

以下只列出了现在修改的相关代码，首先将求质心修改为求二值轮廓最大内接圆的圆心

def getCentriod(curSliceBoundMap):
    '''
    :param curSliceBoundMap:当前层的二值轮廓图
    :return: 当前层最大接圆圆心以及半径
    '''
    curSliceCopy = curSliceBoundMap.copy()
    curSliceCopy = curSliceCopy.astype(np.uint8)  # 将float64转为uint8才能调用findcontous
    contours, _ = cv2.findContours(curSliceCopy, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    if len(contours) == 1:
        raw_dist = np.empty(curSliceCopy.shape, dtype=float)
        # 二重循环遍历当前矩阵中的每个点
        for i in range(curSliceCopy.shape[0]):
            for j in range(curSliceCopy.shape[1]):
                raw_dist[i, j] = cv2.pointPolygonTest(contours[0], (j, i), True)  # 计算二维矩阵中的点在轮廓上还是轮廓内、外
    # 获取最大值即内接圆半径，中心点坐标
    minVal, maxVal, _, maxDistPt = cv2.minMaxLoc(raw_dist)  # 返回值是距离矩阵中的最小值、最大值、最小值的坐标、最大值的坐标
    maxVal = abs(maxVal) # 取绝对值
    drawing = curSliceBoundMap.copy()
    drawing[drawing == 1] = 128
    cv2.circle(drawing, maxDistPt, np.int16(maxVal), (255, 255, 255), 1, cv2.LINE_8, 0)
    # 找drawing中最小的行、列
    rows, cols = np.where(drawing == 255)
    minRow = min(rows)
    maxRow = max(rows)
    minCol = min(cols)
    maxCol = max(cols)
    #圆心你
    centerRow = (minRow + maxRow) / 2
    centerCol = (minCol + maxCol) / 2
    center = [centerRow,centerCol]
    return center, maxVal

根据两向量求当前方向夹角

def getTheta(baseVector, curVector):
    '''
    :param baseVector:  基准向量
    :param curVector: 当前向量
    :return: 向量之间的夹角
    '''
    # baseLength = np.sqrt(baseVector[0] ** 2 + baseVector[1] ** 2) # 单位向量模长为1,可以省略
    curLength = np.sqrt(curVector[0] ** 2 + curVector[1] ** 2)

    # 求点积
    vectorMul = curVector[1]
    # 求余弦
    cosTheta = vectorMul / (curLength + 1e-8)
    assert -1 <= cosTheta <= 1, 'cosTheta value {} is exceed'.format(cosTheta)  # 在表达式为false的时候触发异常
    theta = math.acos(cosTheta)  # 反余弦函数参数范围在[-1,1]
    # print('theta:  ', theta)
    return theta

遍历待插值层每一个坐标点，避免出现空洞的情况，考虑为了使得相邻层的中心在一条直线上，待插值层就不重新计算，而是根据上下两层计算距离加权均值坐标

def getInterArea(curSliceBoundMap, preSliceBoundMap, nextSliceBoundMap, curCentriod, preCentriod, nextCentriod, preImage, nextImage, curRadius=0, preRadius=0, nextRadius=0,lambdaValue=0.5):
    # 待插值层中心不用计算,为了保证相邻层的中心在同一直线上
    curCentriod[0] = (preCentriod[0] + nextCentriod[0]) * lambdaValue
    curCentriod[1] = (preCentriod[1] + nextCentriod[1]) * (1 -lambdaValue)
    print(curCentriod, preCentriod, nextCentriod)
    rows, cols = np.where(curSliceBoundMap == 1) # 找到所有轮廓中的点坐标

    if len(rows) == 0 or len(cols) == 0: # 说明没找到轮廓返回空
        return None

    #baseVector = [0,1]  # 用(0,1)计算出的cosTheta超范围

    curSliceBoundMapNew = np.zeros_like(curSliceBoundMap, np.float32)  # 避免修改原图
    for i in tqdm.tqdm(range(len(rows))):
        # 逐行遍历坐标点
        r = rows[i]
        c = cols[i]
        if r == curCentriod[0] and c == curCentriod[1]:
            continue
        # 求当前点和质心连接向量与基准之间夹角
        curVector = [(r - curCentriod[0]), (c - curCentriod[1])] # 终点减去起点
        # theta = getTheta(baseVector, curVector)
        #theta = getThetaAngle(baseVector, curVector)
    

        length = np.sqrt(curVector[0]**2 + curVector[1]**2) # 求模长
        # 避免除0操作
        assert 0 < length, 'length {} is exceed'.format(length)
        sinCos = np.array(curVector) / length  # 正余弦

        curDirectedPoints = getEdgePointAngle(curSliceBoundMap, sinCos, [r, c])
        preDirectedPoints = getEdgePointAngle(preSliceBoundMap, sinCos, preCentriod)
        nextDiretedPoints = getEdgePointAngle(nextSliceBoundMap, sinCos, nextCentriod)

        curRadius = getLength(curDirectedPoints[0], curDirectedPoints[1], curCentriod[0], curCentriod[1])
        preRadius = getLength(preDirectedPoints[0], preDirectedPoints[1], curCentriod[0], curCentriod[1])
        nextRadius = getLength(nextDiretedPoints[0], nextDiretedPoints[1], curCentriod[0], curCentriod[1])

        curLength = getLength(r,c,curCentriod[0],curCentriod[1])
        curRatio = curLength / curRadius  # 求当前点和当前半径的长度比,四舍五入的误差,边界点选择向上取整
        # 比值标准差
        assert 0 < curRatio <= 1, 'curRatio value {} is exceed'.format(curRatio)
        # 根据theta求上下两层对应方向对应比列的点
        preMatchPoint = getMatchPointAngle(preSliceBoundMap, theta, curRatio, preCentriod, preRadius)
        nextMatchPoint = getMatchPointAngle(nextSliceBoundMap, theta, curRatio, nextCentriod, nextRadius)
        # 根据上下两层所求匹配点当前层八邻域的距离加权平均
        preMeanValue = getNeighGreyValue(preMatchPoint, preImage, preSliceBoundMap)
        nextMeanValue = getNeighGreyValue(nextMatchPoint, nextImage, nextSliceBoundMap)
        # 上下层距离加权平均
        meanValue = lambdaValue * preMeanValue + (1 - lambdaValue) * nextMeanValue
        #print('当前灰度值:  ', meanValue)
        curSliceBoundMapNew[r][c] = meanValue
    return curSliceBoundMapNew

效果

求最大内接圆圆心

插值层

可以看出相较上一篇的效果，经过此次修改，空洞和稀疏的情况得到很大改善了，但还存在一些问题（上图中箭头所指两处），考虑此情况是由于在同方向上找边界，由于耳朵位置造成上一层和下一层的边界位置不都在耳朵位置取得，这是由于不同层耳朵轮廓不同造成的

目前遗留问题

上图中左边表示上一层，根据中间待插值层当前待插值点与中心所在方向，左右两幅图求的边界不一致，导致这一方向的点所求的距离加权平均灰度误差较大

针对以上问题，目前考虑的解决方案

• 寻找边界的优化，将二值轮廓限制在一个范围内，边界点为当前范围内最远的边界而不是第一个边界
• 直接去掉耳朵区域的影响，我们考虑在分割脑卒中病灶其病灶区域并不会出现在耳朵区域，所以可以直接去掉耳朵位置区域

下周再继续修改此问题