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的设计