Carm图像矫正

C-arm畸变矫正及矫正参数求解

• 前面已经讲过相机图像矫正的基本原理：
  通过求解extrinsic参数 e.g T，R 和 intrinsic参数 e.g m，s
  来对畸变图像（Carm为点光源成像（由一点发射X射线），虽然接收器会有一定程度的矫正，但是边 缘区域还是会有“比较大”的畸变 e.g. barrel）进行矫正。下图中可以看出，红色线段代表畸变的程度，越靠近边缘处，畸变越严重。
  

  在项目中，我们使用的是Imag Adapter，通过特定已知图例，对矫正参数进行求解

• 首先需要对C-arm中的图案进行采集（比如下图中bead中心坐标）：
  不需要花太多力气在这个上面，
  如果想研究具体算法的同学，可以参考 openCV cv2::findContour, 和一个imutils的开源库 (reference
  [http://www.pyimagesearch.com/2015/02/02/just-open-sourced-personal-imutils-package-series-opencv-convenience-functions/])
  

def findBeadCoordinates(img):
    img = 255-img #invert the black and white for imutils
    blurred = cv2.GaussianBlur(img, (5, 5), 0)
    thresh = cv2.threshold(blurred, 60, 255, cv2.THRESH_BINARY)[1]
    cv2.imshow("thresh.img", thresh)
    cv2.waitKey()
    cnts = cv2.findContours(thresh.copy(), cv2.RETR_EXTERNAL,
                            cv2.CHAIN_APPROX_SIMPLE)
    cnts = imutils.grab_contours(cnts)
    print("cnts = " , cnts ,"\n" ,"number of bead is = ", len(cnts))
    for c in cnts:
        # compute the center of the contour
        M = cv2.moments(c)
        cX = int(M["m10"] / M["m00"])
        cY = int(M["m01"] / M["m00"])
        # draw the contour and center of the shape on the image
        cv2.drawContours(img, [c], -1, (0, 255, 0), 2)
        cv2.circle(img, (cX, cY), 7, (255, 255, 255), -1)
        cv2.putText(img, "center", (cX - 20, cY - 20),
                    cv2.FONT_HERSHEY_SIMPLEX, 0.5, (255, 255, 255), 2)
        # show the image
        cv2.imshow("Image", img)
        cv2.waitKey(0)

"""
结果为：
[[579, 571], [1130, 563], [755, 562], [954, 553], 
[401, 399], [572, 395], [949, 390], [748, 390], 
[246, 247], [393, 243], [738, 242], [562, 237], 
[106, 104], [383, 103], [242, 103], [546, 102]]
"""

在确定bead的坐标后，可以根据bead geometric pattern进行求解图像矫正参数：
Bertelsen, A., et al. “Distortion correction and calibration of intra-operative spine X-ray images using a constrained DLT algorithm.” Computerized medical imaging and graphics 38.7 (2014): 558-568.

对extrinsic parameters的求解：

 number of bead is =  16
[[[106, 104], [242, 103], [383, 103], [546, 102]],
 [[246, 247], [393, 243], [562, 237], [738, 242]],
 [[401, 399], [572, 395], [748, 390], [949, 390]],
 [[579, 571], [755, 562], [954, 553], [1130, 563]]]

假设对应的实物坐标为
[[[10, 10], [20, 10], [30, 10], [40, 10]],
 [[20, 30], [30, 30], [40, 30], [50, 30]],
 [[30, 50], [40, 50], [50, 50], [60, 50]],
 [[40, 70], [50, 70], [60, 70], [70, 70]]]

选取三个球来求解即可，选择1，6，9号位点，：

下面是对instrinc parameters的求解

考虑到C-arm adapter和intensifier影增之间是几乎平行放置，所以图像的skewness基本可以忽略（对应K matrix的（1，2））
对应intrinsic parameter 三个variables可以用DLT方法进行求解

至此，我们可以得到C-arm图像的矫正矩阵，根据实际adapter geometric pattern
但是也可以看出，geometric pattern图列的排布，和求解extrinsic parameter坐标点的选取，也会影响最后的结果，需要预先知道C-arm image畸变的类型和分布（如图1），再进行C-arm image calibration geometric pattern的设计