Automatic Distortion Correction of Endoscopic ... - Semantic Scholar

3 downloads 142700 Views 4MB Size Report
We propose a fully automatic Hough-entropy-based calibration algorithm, which ... obtained from the IEEE by sending an email to [email protected].
This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

1

Automatic Distortion Correction of Endoscopic Images Captured with Wide-Angle Zoom Lens Tung-Ying Lee, Tzu-Shan Chang, Chen-Hao Wei, Shang-Hong Lai, Kai-Che Liu, Hurng-Sheng Wu

Abstract—Operation in minimally invasive surgery is more difficult since the surgeons perform operations without haptic feedback or depth perception. Moreover, the field of view perceived by the surgeons through endoscopy is usually quite limited. The goal of this paper is to allow surgeons to see wideangle images from endoscopy without the drawback of lens distortion. The proposed distortion correction process consists of lens calibration and real-time image warping. The calibration step is to estimate the parameters in the lens distortion model. We propose a fully automatic Hough-entropy-based calibration algorithm, which provides calibration results comparable to the previous manual calibration method. To achieve real-time correction, we use graphics processing unit to warp the image in parallel. In addition, surgeons may adjust the focal length of a lens during the operation. Real-time distortion correction of a zoomable lens is impossible by using traditional calibration methods because the tedious calibration process has to repeat again if focal length is changed. We derive a formula to describe the relationship between the distortion parameter, focal length, and image boundary. Hence, we can estimate the focal length for a zoomable lens from endoscopic images on-line and achieve real-time lens distortion correction. Index Terms—Endoscopic Image, Optical Distortion, WideAngle Distortion Correction, Zoom Lens Calibration

I. I NTRODUCTION

M

INIMALLY invasive surgery (MIS) is one of the main evolutions of surgical techniques and it provides great benefits to the patient. MIS uses an endoscopic camera and instruments manipulated through small incision on the skin to avoid open surgery. This technique reduces patient scars, infection risk, post-operative morbidity and patient recovery. Although a considerable amount of intervention is still performed with standard open surgery techniques, a large amount of research work is performed every year to develop new instrumentations that replace open surgery procedure by MIS procedure and to make the MIS interventions easier. Indeed, MIS procedures are still more difficult to perform Manuscript received November 3, 2012; revised March 19, 2013. This work was partially supported by National Science Council, Taiwan, R.O.C., under the grant 99-2622-E-007-019-CC2 and the grant 98-2221-E-007-089-MY3. T.-Y. Lee, C.-H. Wei, and S.-H. Lai are with Department of Computer Science, National Tsing Hua University, Hsinchu 300, Taiwan (email: [email protected], [email protected], [email protected], phone: 886-35742958). T.-S. Chang was with Institute of Information Systems and Applications, National Tsing Hua University, Hsinchu 300, Taiwan (email: [email protected]). K.-C. Liu and H.-S. Wu are with Asian Institute of TeleSurgery, Chang Bing Show Chwan Memorial Hospital, Changhua 505, Taiwan. Copyright (c) 2013 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending an email to [email protected].

than standard surgery for many reasons. Firstly, surgeons lose depth perception since the camera is monocular and the visual resolution is reduced since the video sensor resolution is still not able to match human eye resolution. Secondly, they lose the haptic feedback since they cannot touch the organs with their hands. Instead, they can obtain haptic feedback with some particular instruments. However, feeling an artery position from blood pulse, for example, is currently not feasible with instruments. Thirdly, the camera field of view is generally quite limited compared to human perception (70 ◦ against 160 ◦ for human eyes) for almost all abdominal MIS procedures. The main reason is that the wide-angle lens introduces image distortion: straight lines become curved in the distorted image, thus seriously decreasing the hand eye coordination of the surgeon and degrading the capability of size and shape determination. This small field of view for an endoscopic camera has several drawbacks. When the surgeon introduces the endoscopic camera, a long exploration phase is necessary to localize each organ position and status. And even after this step, the surgeon can never have a broad view of the patient’s organs that allows him to make faster movement. A second drawback is the instrument occultation when they are not in the camera field of view. For some procedures, more than 5 instruments are used at the same time (and hold by assistants) and they cannot be visible at the same time. Hidden instruments should not be moved since they can damage the closest organ. This means that the intervention is slowed down because all the instruments cannot be moved at the same time. Moreover, since the camera holder usually does not manipulate instruments, a very good coordination is necessary among the whole surgical team so that the camera is always focused on what is important during the surgical procedure. The goal of this paper is to allow surgeons to use an endoscopic wide-angle lens without the drawback of the lens distortion. Since the current cameras are digital, removing distortion can be realized via image processing. In this context, the goal of the paper is to develop a system which can efficiently compensate the wide-angle lens distortion under different zooming conditions. This system is divided into two parts, calibration part and image correction part. In the calibration part, the system captures several images which contain the calibration pattern, and then the distortion parameters of a lens are estimated. These distortion parameters are used to undistort image in the correction part. Although there are existing calibration methods and interpolation methods, there are several challenges to apply them to the practical problem. In practice, current endoscopy provides 60Hz frame rate for visual comfort

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

at HD resolution. It is almost impossible to achieve real-time image correction by using only central processing unit (CPU) for the computation due to the large amount of computation involved in the image correction. Secondly, most calibration methods are manual or semi-automatic, and they need some user intervention. Furthermore, traditional methods are still not feasible for a zoomable lens even though the lens can be calibrated before operation. It is because the distortion parameters will be changed when surgeons change the focal length of a lens. In order to make the system practical, a reliable, fast, and automatic calibration is required. Wide-angle cameras have been widely used in video surveillance, automotive applications, and endoscopic imaging [1]– [12]. A large field of view is obviously needed in these applications. For the application of 3D reconstruction [2], [3] and mosaicing [11], the reconstruction and mosaicing algorithms are based on the undistorted images. Undistorted images are also needed when the system uses several different modalities of images to achieve virtual guidance surgery [4]. The impact of distortion for automatic classification of celiac disease was discussed in [5], [6]. For a lens without zooming, the distortion parameters of an endoscope could be estimated offline. There exist some systems of real-time image distortion correction [4], [7], [9] for endoscopy without zooming. In the system proposed by Shahidi et al. [9], the endoscope was calibrated by using Tsai’s method with some 2D-to-3D point correspondences determined manually. Some systems [8], [10] used dot patterns for calibrating a wide-angle lens, they require careful setting of endoscopy during calibration or assume dot extraction is under controlled environment such that these dots could be organized in a vertical or horizontal fashion. Moreover, the parameter-free distortion correction methods [13], [14] require point correspondences between the input image and calibration pattern, which usually requires user interaction. Thus, an automatic and robust image calibration method is strongly demanded. A zoomable lens further complicates the problem. The calibration and correction of images acquired with a zoomable lens involves two difficulties. Firstly, it requires a mechanism to calibrate the camera for all possible zoom settings, which either is quite labor intensive [15] or needs special hardware to obtain focal length and radial distortion [16]. These methods are not suitable for calibrating a zoomable endoscope. The second difficulty is to estimate zoom setting of a lens online even if a zoomable lens is fully calibrated. In this paper, we propose a single-image Hough-based autocalibration for calibrating a wide-angle endoscope with a certain focal length. It is robust against noise and occlusion. Unlike the method by Shahidi et al. [9] that uses small line segments to estimate radial distortion, our method automatically estimates the radial distortion from several distorted straight lines. If the pattern consists of parallel lines, our framework can also utilize the prior knowledge of parallel line to improve the calibration accuracy. We propose a twostage Hough transform based method for the estimation of the distortion parameters. Furthermore, we extend the proposed system to calibrate the distortion parameters for a zoomable

2

lens. To describe the relationship between the wide-angle distortion and focal length, we derive a formula that simplifies the process of zoom calibration, and just a few zoom settings need to be calibrated. Hence, we can estimate focal length online and achieve real-time distortion correction. In this paper, we demonstrate the proposed algorithm can achieve real-time correction of wide-angle lens distortion for a zoomable lens. For the rest of this paper, we present the single-view Houghbased auto-calibration method in Section II. The relationship between the distortion parameter and focal length is derived in Section III. The flow of the whole system is described in Section IV. In Section V, we will the accuracy of the proposed auto-calibration method and the performance of the system. A preliminary version of this work, given in parts of Section II, was reported in [17]. In this paper, we propose an extended algorithm based on the two-stage Hough transform to take advantage of the parallel-line calibration patterns and develop the calibration and correction framework for a zoomable lens. II. S INGLE -I MAGE H OUGH -BASED AUTOCALIBATION Camera calibration usually uses two kinds of markers, corners and straight lines. A popular calibration toolbox [18] uses corners of chessboards for establishing point correspondences. Although some existing camera calibration methods utilize automatic chessboard detection [19], their performance usually depends on the correctness of chessboard detector, which is easily influenced by illumination. The line-based camera calibration is another popular approach [20]–[23]. In structured environments or straight-line patterns, using line features could be a more suitable option. The projected lines without distortion should be straight; hence, this principle could be used in removing optical distortion [22], [23]. Hough transform [24] is a powerful tool for extracting straight lines. Several previous works have used Hough transform to find the perspective transform [25] or vanishing points [26], [27]. From the distorted images, Habib and Morgan [28] detected clusters of distorted lines and Xu et al. [29] just used Hough transform to find the distortion center given all other distortion parameters. The most related works to our calibration method are the HTRDC method proposed by Cucchiara et al. [30] and 1-D angular Hough method [31]. These two methods are based on polynomial distortion model and Harris model. The HTRDC method focuses on using a single straight line, and tries to adjust the distortion parameter to maximize the maximal vote assuming the distortion center is known. The 1-D angular Hough method requires less space. This method assumes the orientations of lines are sparsely distributed in the angular Hough space. However, this could be violated due to the perspective effect. A plane with several parallel lines could be served as a pattern for auto-calibration (we denoted it by parallel-line-pattern), but it is a good example that the orientations of lines are not sparsely distributed if the plane of calibration pattern is not parallel to the image plane. We propose a Hough entropy based method considering all distorted straight lines, and it estimates the distortion

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

3

Fig. 1. Degenerate problem: Five Hough maps are associated to different distortion parameters (λ): −3 × 10−5 , −2.5 × 10−5 , −2 × 10−5 , −1.5 × 10−5 , −1 × 10−5 . For each Hough map, we show zero counters in whole map (green border) and its central region (red border). The entropy measure is highly influenced by the number of the non-zero region of Hough map.Although the third map is the densest among the five maps in central region,the fifth map has the smallest entropy because its map has a large number of zero cells.

parameter and the distortion center. The entropy in the Hough space is employed to measure straightness of the undistorted curves. In a naive way, all edge points are undistorted by possible parameters and then the best parameter setting is selected based on the entropy in the Hough space. However, the results will be degraded by multiple lines occurred in the image, because the entropy will be lower when the image is distorted to a smaller size. To overcome this problem, we only count the votes with most possible orientation for each edge point. In addition, the image gradient computed from the distorted image is also distorted. Hence, we derive the undistorted gradient for more compact Hough voting, and it also solves the aforementioned degenerate problem. In the following subsections, we will introduce some preliminaries, including distortion model and Hough transform, and then the Hough entropy method and gradient estimation are presented 3subsequently. A plane with parallel lines (Parallel Line Pattern) could be served as the calibration pattern in our system. We use the two-stage Hough transform to incorporate the prior of this type of calibration patterns. The first layer Hough map collects votes for straight lines, and the second layer Hough map is used for accumulating pencils of lines. The edge points produced by noise or unrelated curves will be suppressed in the second map. It can make the calibration more accurate. A. Pin-Hole Camera Model and Division Model A popular and simple pin-hole camera model contains focal length f and assumes that camera is located at the origin and faces to negative Z-axis. An object location with respect to the camera coordinate is denoted by (X, Y, Z). The corresponding image point is inversely proportional to the depth Z, i.e. x = −f X/Z, y = −f Y /Z. If the camera coordinate and world coordinate are different, there exists a rotation matrix R and translation vector T for the coordinate transformation. Optical distortion can be easily incorporated into the model. Two kinds of optical distortions, radial distortion and tangential distortion were modeled in [32]. The radial distortion at a point is represented by a displacement along the direction of the point to the camera center. The tangential direction is perpendicular to the radial direction. However, low tangential distortion could be compensated by estimating the distortion

center [33]. The relationship between the distorted radius of an image point and the undistorted radius is approximated by a low-degree polynomial. Many researchers have found that the higher-order terms could be ignored [34]. In addition, Fitzgibbon [35] suggested a division model which uses a loworder polynomial in the denominator instead of the numerator. The corresponding undistorted coordinate of an image coordinate xd = (xd , yd ) is defined by xu = (xu , yu ). The division model [35] used to model the optical distortion relates these two coordinates as follows: rd ru = , ru = kxu − ck, rd = kxd − ck (1) 1 + λrd2 where c = (cx , cy ) is the image center. Equation (1) gives the relationship of distorted radius rd and undistorted radius ru. For wide-angle distortion correction by using the division model, we just need to determine the image center (cx , cy ) and distortion parameter λ. B. Hough Transform Traditional Hough transform is to transform image features in an image domain to a line parameter space to vote for objects, i.e. lines. Any one line passing through (x, y) can be represented in polar form (r, θ), and the corresponding line equation is r = x cos θ + y sin θ, and all possible lines can be collected by Scont (x, y) = {(r, θ)|r = x cos θ + y sin θ}. For each edge point (x, y) in the image, all lines passing through (x, y) will receive one vote individually. Finally, the Hough map stores the total number of votes for each line (r, θ). In the computer program, the polar coefficients are discretized and a Hough map is a 2D accumulation array used for recording votes. The possible lines passing through (x, y) can be written by a transfer function S. S(x, y) = {(xcosθi + ysinθi , θi )|∀θi }

(2)

where θi is the sampled angle in [0, π]. The voting process for the edge point (x, y) is denoted by H(r, θ) ← H(r, θ) + vote(S, (x, y), inc)

(3)

where H is an accumulation array and inc is the increment (In traditional Hough transform, inc is equal to 1).

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

Eps/This pdf

is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816

Fig:SPLP eps/pdfON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE. TBME-01659-2012.R1. IEEE TRANSACTIONS ωi

V(ωi)

θ

r L

L

ρi

P

ρi

M

Q

θ

+H1(r, θ)

H1

Original

ωi

H2

SPLP

We propose the Hough entropy method with image gradient estimation for distortion parameter estimation. We found that sometimes the result is degraded because of the votes in the Hough space which are generated from some orientation not similar to the edge direction. Hence, we derive the gradient estimation formula and just vote the most possible orientation for each edge point. C. Hough Entropy Method First, the edge points which could belong to the distorted straight lines are extracted by using Canny edge detector. A typical range of distortion parameters and the possible regions of image center are selected. A uniform sampling is performed in this 3D parameter space. For each candidate parameter triplet (λ, cx , cy ), all edge points are undistorted as follows:   yd − cy x d − cx + cx , + cy undist(xd , yd ; λ, cx , cy ) = 1 + λrd2 1 + λrd2 (4) Based on these undistorted points, a Hough map H(r, θ) is obtained and normalized by dividing the total number of all votes. The entropy of the normalized Hough map is calculated with the following entropy formula. X p(r, θ) = H(r, θ)/ H(i, j) (5) ∀i,j

X

−p(r, θ) ln p(r, θ)

Edge map

H1 map

H2 map

Corrected

Fig. 3. An example of distortion correction and the corresponding H1 and H2 maps.

Fig. 2. Given a line L in the first Hough map H1 , the H2 and line M , QV (ω1 ), the corresponding cell H2 (ω, ρ) in the second Hough map is voted by the weight of line L.

Entropy(H) =

4

(6)

∀(r,θ),s.t.p(r,θ)6=0

Finally, the parameter triplet with the minimum entropy is selected as the best distortion parameter and image center. D. Hough Entropy with Gradient Estimation (HEwG) Fig. 1 illustrates the degenerate problem occurred in Hough entropy method. When the distortion parameter is large, the undistorted image will decrease its size. Hence, the zero counters in Hough map will increase and entropy decrease. The modified version of Hough entropy method is to use image center and distortion parameter to estimate the orientation of an edge point in undistorted image. Undistorting an image costs too much computation power, because our method examines many candidates in the parameter space. Hence, we derive how to use the image gradient in the division model to simplify the computation in the Hough transform. For estimating the orientation, we first calculate image gradient on distorted image by using the Sobel masks. Let

us focus on the edge point (xd , yd ). The normalized gradient on the distorted image is denoted by (gxd , gyd ). A neighboring point (x0d , yd0 ) along a curved line direction (perpendicular to gradient) is (xd , yd ) + (−gyd , gxd ), where  is a small number. By using the division model, we can obtain undistorted positions, (x0u , yu0 ), of these two points by using the undistortion formula. The gradient on the undistorted image, (gxu , gyu ), can be estimated by lim→0 (−(yu0 − yu ), x0u − xu ). By letting  approach to zero, we can obtain the gradient estimation formula by l’Hospital’s rule. (gxu , gyu ) ∼ (−(λA(yd − cy ) + gxd β), λA(xd − cx ) − gyd β), (7)   where A = −2 (yd − cy )gxd − (xd − cx )gyd and β = (1 + λrd2 ). For evaluating a parameter triplet, edge pixels (xu , yu ) and the associated image gradients (gxu , gyu ) are undistorted by equation (4) and (7), respectively. For each edge point, only the vote that corresponds to the estimated gradient is applied to the Hough space. SGE (x, y, gx , gy ) = {(r, θ)|r = x cos φ + y sin φ}   q φ = cos−1 gx / gx2 + gy2

(8) (9)

The best parameter setting is selected by finding the minimal entropy measure. E. Incorporation of Hough Framework and Parallel-LinePatterns (PLPs) If the calibration pattern is PLP, we employ the two-stage Hough transform to incorporate the prior of this type of calibration patterns. The edge points produced by noise or unrelated curves will be suppressed in the second map. After performing the first Hough transform in Section II-D, the output Hough map H1 is used as the input to the second Hough transform. The proposed algorithm is to evaluate the entropy of the second Hough map H2 . The projected version (projected on the image sensor) of a set of 3D parallel lines can be represented by a vanishing point V in the image space. However, the range of vanishing point is very large, even unbounded. To overcome this problem, we introduce two auxiliary points P and Q , i.e., (xP , yP ) and (xQ , yQ ), in the image for the new representation of the second Hough transform. In practice, we select the left-upper image corner as P and the image center as Q. Without loss of generality, the vanishing point V can be represented by the orientations of the two lines, P V and QV , which are denoted by (ω, ρ). We use this representation in the second Hough map

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

H2 . Note that the map is bounded because orientation angles are limited within 360 (or 180) degrees. Given the first Hough map H1 , the voting process is given as follows. A cell in H1 at position (r, θ) represents a line with supporting votes H1 (r, θ) as the line L shown in Fig. 2. All possible vanishing points of lines parallel to L are on the line L. We can vote them by using two-orientation representation. For each sample ωi , we can construct the line M passing through P with orientation ωi , ρ = xP cos ωi + yP sin ωi . The corresponding vanishing point V (ωi ) is the intersection of the line L and M . The orientation ρi of the line QV (ωi ) can be computed. The location (ωi , ρi ) in H2 corresponds to H1 (r, θ) votes. Based on the theory of projective geometry, the intersection of two lines and the line crossing two points are easily calculated by cross-product in homogeneous coordinates. Thus, we define a transfer function SP LP with the following equations. SP LP (r, θ) = {(ωi , ρi )|ρi = ρ(ωi )} −1

ρ(ω) = tan

(n1 (ω)/n2 (ω))

(n1 (ω), n2 (ω), n3 (ω)) = (xQ , yQ , 1) × V (ω) V (ω) = (cos ωi , sin ωi , −rω (ω)) × (cos θ, sin θ, r) rω (ω) = xP cos ωi + yP sin ωi

(10) (11) (12) (13) (14)

Based on the transfer function SP LP , the line L will contribute H1 (r, θ) votes at all locations in SP LP . Finally, we calculate the entropy of H2 instead of that of H1 . In Fig. 3, we depict an example of distortion correction on a chessboard image and the corresponding H1 and H2 maps. In the cascaded Hough transform [27], several Hough transforms are also applied cascadedly. However, our framework is different from the traditional cascaded Hough transform. In the first Hough map, we use radial parameterization for line, not the slope form in [27]. In the second Hough space, we introduce two auxiliary points and use two relative angles instead of the vanishing point. Hence, these two Hough spaces are limited. We do not need to fold the Hough space as that used in [27]. F. Optimization Our camera distortion calibration is an optimization process to find the best distortion parameters based on the proposed measures. Our measures are evaluated by entropy of Hough map, and it is not easy to apply continuous optimization to solve the problem. The coarse-to-fine search is a possible strategy [30]. We modified this hierarchical search as follows. First, we uniformly sample over a bounded parameter space. The sample with the minimal entropy is selected, and then we reduce the parameter search space to the neighborhood of this sample. The next round of sampling, evaluation, and reduction of parameter search space is repeated until the size of the parameter search space is less than a user-defined accuracy. The Hough entropy method with gradient estimation and the version of incorporation of PLP prior are summarized in Algorithm 1, 2, and 3 (shown in Fig. 4, 5, and 6). In our implementation, the typical range of distortion param2 eter is given by [−1/rmax , −10−7 ] where rmax is the distorted

5

Algorithm 1: Hough Entropy Framework Input: the initial range of distortion parameter [λmin , λmax ], the initial range of distortion center in x coordinate [xmin , xmax ] and in y coordinate [ymin , ymax ], image I, sample number S, user-defined tolerance  Output: the best distortion parameter λ, the best distortion center (cx , cy ) 1: Extract the set of edge points E from I, and the associated gradient Ix and Iy . 2: for l = 1 to 20 do 3: Sample S samples over [λmin , λmax ] 4: Sample S samples over [xmin , xmax ] 5: Sample S samples over [ymin , ymax ] 6: Generate S 3 possible parameter configurations from previous samples. 7: for each configuration (λi , cix , ciy ) do 8: if method is PLP extension then 9: Run Algorithm 3 10: else 11: Run Algorithm 2 12: end if 13: end for 14: Find the best configuration minimizing Hough entropy he 15: For distortion parameter and distortion center use two nearest samples to update λmax , λmin , xmax , xmin , ymax , ymin 16: if (λmax − λmin ) <  then 17: break 18: end if 19: end for Fig. 4. Hough Entropy Framework

Algorithm 2: Hough Entropy with Gradient Estimation Input: the configuration (λi , cix , ciy ) Output: the associated Hough map H1 and Hough entropy 1: Use this configuration to undistort the edge point and calculate undistorted gradient by using (4) and (7). 2: Vote for each undistorted point by using (10) in H1 . 3: Estimate Hough entropy of H1 . Fig. 5. Hough Entropy with Gradient Estimation

radius of visible regions, the possible range of distortion center is limited in a 0.1 × W by 0.1 × H region which is located in image center. When the subrange of distortion parameter λ is less than , the optimization is finished. In our implementation,  is set to −10−9 , and the number of iterations required to meet this convergence criterion is usually no more than 7. III. E XTENSION TO A Z OOMABLE L ENS A. Quadratic Prediction of Distortion Parameters We extend our algorithm to a zoomable lens. Consider the pin-hole camera model. Assume the focal length is changed to rf and the center of camera is stationary. The new projected point (xr , yr ) satisfies the following equations, xr = −rf X/Z

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

Algorithm 3: Extension of Parallel-Line Patterns Input: the configuration (λi , cix , ciy ) Output: the associated Hough entropy with incorporation of parallel-line patterns 1: (H1 , Entropy(H1 )) = Run Algorithm 2 2: for each line (r, θ) in H1 do 3: Find out all corresponding cells in H2 by using SP LP , and increase these cells by H1 (r, θ) 4: end for 5: Estimate Hough entropy of H2

Fig:zoomlens eps/pdf

Fig. 6. Extension of Parallel-Line Patterns

Fisheye lens

Zoomable lens

Radial distortion

Perspective projection

E1.F1

E1.F2

6

E2.F3

E2.F4

Fig. 8. Real endoscopic images with different focal lengths. TABLE I R EAL ENDOSCOPIC IMAGES AND PREDICTION OF DISTORTION PARAMETERS : ONE OF IMAGES CAPTURED FROM AN ENDOSCOPE IS USED TO PREDICT DISTORTION PARAMETER OF THE OTHER . Image

λ

E1.F1 E1.F2 E2.F3 E2.F4

-2.9644e-5 -1.2742e-5 -7.5446e-6 -1.6638e-6

Diameter of Circle 210 (px) 312 (px) 432 (px) ∼ 712 (px)

Linear Prediction -1.8931e-5 -1.9953e-5 -3.0577e-6 -4.1052e-6

Quadratic Prediction -2.8126e-5 -1.3430e-5 -5.6195e-6 -2.2338e-6

Fig. 7. A zoomable fisheye lens.

and yr = −rf Y /Z. Most lenses for photography follow this rule. Inspired by Kannala’s formulation [36], we divide the projection part into radial distortion and perspective projection as shown in Fig. 7. It means the incoming rays to the camera center are refracted and then the refracted rays are perspectively projected on the sensor. Originally, the division model is unrelated to focal length. Now we will derive the relationship between the change of distortion parameters and the change of focal length. Consider focal length f and rf . Based on our assumption, the image J with focal length rf could be approximated by the linearly zoomed version of image I with focal length f . A 3D point X is projected to xfd and xrf d in the image I and J, respectively, and these projected points and their undistorted versions can be well approximated with the following equarf f f f rf tions, xrf d = rxd , xu = rxu , and rd = rrd . We obtain the relationship between the change of distortion parameters and the change of focal length. The distortion parameter is inverse-quadratically proportional to focal length. xrf u =

xrf d 1 + (λf /r2 )(rdrf )2

⇒ λrf = λf /r2

(15)

We observe real endoscopic images with different focal lengths. In Fig. 8, we have two kinds of endoscopes E1 and E2 and four focal lengths F1, F2, F3, and F4. We use a circle to fit the boundary of the endoscope in the image. Given distortion parameter of one image in an endoscope, the ratio of two diameters of circles is used to predict the distortion parameter of the other image. We compare the linear prediction model with the derived quadratic prediction model in Eq. (15), and found the quadratic prediction is much better than naive linear prediction. The detailed information is summarized in Table I. B. Practical Issues Based on the above quadratic prediction model, we can use finite samples to predict distortion parameters of arbitrary continuous focal length.

In the previous examples, we extract the boundary of an endoscope to estimate the focal length. There are two practical issues about extracted boundary, completeness and existence of boundary. First, sometimes extracted boundary of the circle is not complete due to cropping by image sensor. For incomplete boundary of circle, the area of visible region is another good measure instead of diameter of circle and linear prediction can be applied. If the boundary is not available, a simple modification for endoscopes can realize the estimation step. For example, a small marker could be put behind the fisheye lens for online recognition. The misalignment is another issue. The center of zoomable device and the distortion center of a frontend fisheye lens may not be perfectly aligned. The distorted centers are located in a line at different focal lengths. To overcome this problem, we can linearly interpolate distortion center (cx , cy ) by using two sets of calibration parameters with two nearest focal lengths. IV. T HE S YSTEM O F R EAL -T IME Z OOMABLE W IDE A NGLE L ENS C ORRECTION In this section, we describe the whole distortion correction system. First, the distortion correction system with fixed focal length is introduced. Then, the system is extended to a zoomable lens. Finally, we briefly describe how we achieve real-time image distortion correction by using GPU.

A. Fixed Focal Length The system consists of two parts, calibration and image correction. The calibration could be operated offline. We just need take a single image which contains the calibration pattern, and then perform the proposed Hough-based autocalibration procedure. The distortion parameters are used for correcting lens distortion. The pattern could be a chessboard or parallel lines containing at least 8 lines and at most 20 lines, and the captured image is expected to contain curved lines near the boundary of visible regions.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

Eps/ pdf

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816

Fig:flowchart TBME-01659-2012.R1. IEEE TRANSACTIONS eps/pdfON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE. Calibration patterns

Focal length Estimation

Calibration

Distortion parameter table

Offline

Original Original streaming video

7

Focal length Estimation Correction

Online Corrected streaming video

Fig. 9. The flowchart of the proposed distortion correction system for a zoomable fisheye lens.

B. Calibration for A Zoomable Lens For extending the system to a zoomable lens, two components, focal length estimation component and distortion parameter table, are added into the system as shown in Fig. 9. The input image is binarized by an intensity threshold, and we use the structural analysis [37] for finding the outermost contour on the binary image. First, a raster scan is applied for marking some starting points. Then, the boundary is traced from starting points. An open source implementation for boundary tracing in OpenCV is used in this work. The process can achieve 535 fps at 720-by-540 image size. For images of very high resolutions, regular image downsampling by nearest interpolation can be applied first and the boundary detection is performed on the smaller image to reduce the computational cost. The outmost image boundary can be used to calculate the diameter of the circle and the area of visible regions. The focal length is related to the diameter of the circle for the image boundary as given in Eq. (15). In order to construct the distortion parameter table, we first calibrate several images with different focal lengths. The distortion parameters, distortion centers, associated diameters, and associated area values are recorded. For an incoming frame of endoscopy, we first perform boundary detection to find the distortion parameter and distortion center with the closest diameter or with two closest areas in the table. If we select parameters with two closest areas, the distortion parameter is linearly interpolated. Elsewise, the ratio of diameter is served as the ratio of the focal lengths and the distortion parameter is calculated by quadratic prediction in Section III. C. Correction Part Given distortion parameters, the image has to be corrected in real time. The operation involved in the distortion correction algorithm is parallelizable. It consists of computing the corrected pixel coordinates for all pixels and bilinear interpolation. However, the large amount of computation makes real-time processing on CPU impossible. We use GPU to speed up the computation. First, the corresponding coordinates are computed on GPU in parallel. Secondly, we use texture memory to store distorted image. Because the distorted image is just used for read-only and all threads can efficiently access the texture memory, the bilinear interpolation can also be accomplished in parallel. The traditional lookup table (LUT)

HEwG

HT1D

HEPLP

Fig. 10. Corrected results of the dataset s4c by using three methods to estimate the distortion center and parameter. It also shows the degraded result of HT1D due to perspective effect.

implementation is efficient but it requires a large amount of memory. For a zoomable lens, many lookup tables have to be stored for the LUT approach. In this work, we do not store these coordinates in advance. Instead, we compute the corrected image coordinates with GPU on the fly. V. E XPERIMENTAL R ESULTS A. Quantitative Measure and Related Methods for Calibration We implement Wang’s line-based method [23], Hough Transform-based method for Radial Lens Distortion Correction (HTRDC) [30], the measure in 1-D angular Hough method (HT1D) [31], CHT measure [27] for experimental comparison. The single image calibration method EasyCamCalib proposed in [7] needs to extract the corners on the chessboard and cannot handle fronto-parallel configuration. In our dataset, about half of the cases (47.62%) were failed due to noise, specular highlights, low illumination, or orientation of plane of calibration pattern. Hence, we mainly compare our method with line-based methods and Hough-based methods. We use synthetic datasets and real images captured from fisheye lenses and wide-angle endoscopes in our experiments. Wang’s method requires the user to select lines. In HTRDC, the user needs to select a region of interest that contains a single curve. The performance of these two methods depends on the user selection of lines. Our automatic method is based on Hough transform to aggregate the undistorted edge points and uses the undistorted image gradient to discount many error-prone votes. The results show the accuracy and efficiency of the Hough entropy method is considerably improved by properly using the image gradient information. Even if the lines are carefully selected in the other two methods, our method can provide comparable results. For evaluating different methods, we select several curved lines in the distorted images. After the correction, the undistorted points located in the same line are fitted with a straight line by using the total least squares. The root mean square of the distance from the undistorted points to the regression line is served as the quantitative measure RMSE. Our Hough entropy method with gradient estimation is denoted by HEwG, and a version without gradient estimation (HEwoG) is served as a baseline method for comparison. The version that incorporates prior of PLPs is denoted by HEPLP. For a fair comparison of of different distortion correction methods based on HT1D, HEwoG, CHT, HEwG and HEPLP, all the corresponding measures are based on the division model and optimized in the same hierarchical framework.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

8

TABLE II E XPERIMENTS O N R EAL DATASETS G IVEN D ISTORTION C ENTERS

Eps/ pdf

Original 5.00

Fig:syn

HEwGS15

HEwG S7 HEwG S21 HTRDC Avg

6.00

HEwG S15 HTRDC Best

Syn1c

2.00 1.00

HEwGS15

HEwG S7 HEwG S21 HTRDC Avg

6.00 5.00

4.00 3.00

Original RMSE (pixels)

RMSE (pixels)

Fig:syn

HEwG S15 HTRDC Best

4.00

Syn2c

3.00 2.00 1.00

0.00

method Wang(L9) Wang(L13) Wang(L17) HTRDC(Avg) HTRDC(Best) HT1D HEwoG CHT HEwG HEPLP

s1c 1.3534 1.0233 1.0021 2.4697 1.4904 1.2035 1.8230 1.8060 1.1546 1.0105

s2c 1.6120 1.1114 1.0844 2.1852 1.3037 1.6059 2.2818 1.7315 1.1633 1.0672

s4c 1.9891 1.6121 1.6180 4.2826 2.4082 1.8313 1.8443 1.8443 1.7143 1.3652

s3c

s5c

(L4) 0.9099

(L9) 1.2322

(L8) 0.9167

(L13) 1.1250

(L10) 0.9153

(L15) 1.1048

1.7105 0.7425 0.8859 0.9554 0.9554 0.6256 0.6786

3.6953 1.7959 1.0916 4.0931 7.8821 1.2421 1.0191

0.00 Δ=-20

Δ=-10

Δ=0

Δ=10

Δ=20

Δ=-20

Δ=-10

Δ=0

Δ=10

Δ=20

TABLE III E XPERIMENTS O N R EAL DATASETS

Fig. 11. RMSE on synthetic datasets. The first row is syn1c and syn2c and the undistorted results of HEwG S15.

B. Synthetic Datasets We also generate two synthetic datasets, syn1c (ground truth λ = −5.25e − 6) and syn2c (ground truth λ = −8.85e − 6). Because HTRDC needs to know the distortion center in advance, we compare the RMSE for different distortion centers. The distortion center is given by (x + ∆, y + ∆) where (x, y) is correct distortion center. The RMSE of HEwG and HTRDC are shown in Fig. 11. In our method, we try to use different sizes of Sobel operator to reduce the gradient estimation error. Our experiments show that the proposed method in conjunction with 21 × 21 Sobel filters gives the best results.

method HT1D HEwoG HEwG HEPLP

s1c 1.4764 1.9114 1.2112 1.0586

s2c 1.6280 2.3045 1.3168 1.0978

s4c 3.8058 2.1124 1.7860 1.5891

s3c 1.0195 0.9998 0.8164 0.7060

s5c 0.9451 4.0970 1.1475 0.8680

TABLE IV E FFICIENCY E VALUATION O F T HE C ALIBRATION A ND I MAGE C ORRECTION P ROCEDURES I N T HE P ROPOSED D ISTORTION C ORRECTION S YSTEM F OR D IFFERENT C OMPUTING P LATFORMS Number of points Calibration 15299 7758 7029 6100 1-core CPU 6.302 (sec) 3.182 (sec) 2.387 (sec) 1.653 (sec) 4-core CPU 1.997 (sec) 0.89 (sec) 0.671 (sec) 0.531 (sec) Image Correction 1637x1091 2847x1900 6000x6000 1-core CPU 5.097 fps 2.218 fps 0.389 fps 4-core CPU 16.006 fps 8.217 fps 1.714 fps CPU+GPU 64.103 fps 32.083 fps 10.581 fps

C. Real Datasets We have used five different wide-angle lenses to capture the dataset; namely s1c, s2c, s3c, s4c, and s5c. There are 21 images in total, and results for an image in s5c are shown in Fig. 13. All averaged RMSE are shown in Table III and Table II. In Table II, correct distorted centers are given for a fair comparison, because the HTRDC method assumes that the distortion center is given. We compare our method with Wang’s method, HTRDC, HT1D, and CHT. The results show the line selection is very critical to the results of Wang’s and HTRDC methods. However, our method does not require any user interaction, since we use all the edge points. The experimental results show our automatic methods HEwG and HEPLP produced similar results compared to those of Wang17, which requires the user manually select 17 lines on the image. As shown in Fig. 10, the performance of HT1D is degraded if the perspective effect appears in undistorted images, and both the proposed HEwG and HEPLP methods can provide pretty good correction results. The patterns used here are chessboards, and the HEPLP method benefits from the assumption of the parallel-line pattern. In Table III, we compare four methods for determining the radial distortion parameters in the division model. The results show our methods HEwG and HEPLP are very stable and provide the most accurate results when the distortion center is unknown.

TABLE V S UBJECTIVE E VALUATION O N S IX V IDEO C LIPS

Sharpness Lightness Blurring Resolution Correction

Overall

Original videos are better 13% 0% 17% 10% 0%

Corrected videos are better 27% 28% 23% 37% 88%

Not useful at all

Helpful

5%

60%

No significant difference 60% 72% 58% 53% 12% Helpful under some situations 35%

D. Efficiency of Calibration and Image Correction We have tested our distortion correction system on the computer with Intel CPU (i7-860) and Nvidia GPU (C2050). Because the number of edge point will influence the time of calibration, we test on images with different numbers of edge points. The execution time of the proposed calibration algorithm and the image correction procedure for different computing platforms is summarized in Table IV. By combining the GPU and CPU in our implementation of the correction process, the image correction can be accomplished in real time for the undistorted image size of 2847-by-1900.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

Fig. 12. Some snapshots of pairs of original (left) and corrected (right) images, including in-vivo evaluation, calibration, and training.

E. In-vivo Evaluation We also perform in-vivo evaluation for our system as shown in Fig. 12 and Fig. 14. We provide evaluation by ten endoscopy-related professionals, including surgeons and biomedical engineers. Original video streams and undistorted streams are displayed side by side for subjective comparison. The participants compare them and evaluate the sharpness, lightness, resolution, blurring, and distortion of the video streams before and after applying the proposed distortion correction procedure. Finally, the participants give an overall opinion about the distortion correction. We provide six video clips with different types of camera operations, which contain zoom-in, zoom-out, and camera movements. The results of evaluation are summarized in Table V. The evaluation results show the proposed distortion correction can effectively minimize distortion and preserve the image quality compared to the original videos. Most participants think the correction is helpful for improving the overall image quality. The original videos and corrected videos used in the evaluation are available through the link to http://undistort.erufa.com/. VI. C ONCLUSION We proposed a Hough entropy method combined with gradient estimation to correct radial distortion. Our method is automatic and the accuracy is comparable to or better than other line-based methods which require manual selection of lines. The state-of-the-art Hough-transform-based distortion correction method focuses on a region containing a single line. In order to consider all curved lines, we use the entropy in the Hough space as the straightness measure. However, naive extension is not practical because the entropy is easily influenced by the presence of multiple lines. We solve this degenerate problem by using the image gradient estimation. The estimated gradient is derived from the division model and is used to remove the error-prone votes. When the

9

calibration pattern is chessboard or parallel-line-pattern, the proposed method that incorporated prior will further improve the accuracy. Experimental results show the proposed method provides accurate distortion correction results in comparison with the state-of-the-art distortion correction methods. The Hough entropy method is used in our distortion correction system to calibrate a wide-angle lens, and then online correct the distorted streaming video. Furthermore, we extend the system to a zoomable lens. We derive the relationship between the focal length and distortion parameter for the division model. Based on this relationship, we sample some different focal lengths (typically 3 or 4) and use quadratic prediction to achieve good approximation. In this paper, we have shown convincing results and demonstrated it is possible to correct wide-angle optical distortion from images acquired from a zoomable lens in real time. Realtime distortion correction has more applications than offline correction. The accuracy of the calibration and efficiency of image correction in the proposed system are demonstrated on a CPU+GPU computing platform. The main limitation of our method is that the endoscope boundary must be present in the images. In the future, we would like to incorporate other focal length estimators to relax the limitation and make the system more flexible. The proposed distortion correction algorithm takes advantage of the parallel computing power in GPU to achieve real-time performance. It needs to be modified for other computing platforms with different computational architectures for realtime applications. R EFERENCES [1] W. Li, S. Nie, M. Soto-Thompson, C.-I. Chen, and Y. I. A-Rahim, “Robust distortion correction of endoscope,” in Proc. SPIE 6918, Medical Imaging, 2008, pp. 691 812–691 812. [2] T. Stehle, D. Truhn, T. Aach, C. Trautwein, J. Tischendorf et al., “Camera calibration for fish-eye lenses in endoscopywith an application to 3d reconstruction,” in Proc. IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI’07), 2007, pp. 1176–1179. [3] D. Hong, W. Tavanapong, J. Wong, J. Oh, and P. Groen, “3d reconstruction of colon segments from colonoscopy images,” in Proc. IEEE International Conference on Bioinformatics and BioEngineering (BIBE’09), 2009, pp. 53–60. [4] J. P. Helferty, C. Zhang, G. McLennan, and W. E. Higgins, “Videoendoscopic distortion correction and its application to virtual guidance of endoscopy,” IEEE Trans. Med. Imag., vol. 20, no. 7, pp. 605–617, 2001. [5] M. Gschwandtner, M. Liedlgruber, A. Uhl, and A. V´ecsei, “Experimental study on the impact of endoscope distortion correction on computerassisted celiac disease diagnosis,” in Proc. IEEE International Conference on Information Technology and Applications in Biomedicine (ITAB’10), 2010, pp. 1–6. [6] M. Liedlgruber, A. Uhl, and A. V´ecsei, “Statistical analysis of the impact of distortion (correction) on an automated classification of celiac disease,” in Proc. IEEE International Conference on Digital Signal Processing (DSP’11), 2011, pp. 1–6. [7] R. Melo, J. P. Barreto, and G. Falcao, “A new solution for camera calibration and real-time image distortion correction in medical endoscopy– initial technical evaluation,” IEEE Trans. Biomed. Eng., vol. 59, no. 3, pp. 634–644, 2012. [8] K. Vijayan Asari, S. Kumar, and D. Radhakrishnan, “A new approach for nonlinear distortion correction in endoscopic images based on least squares estimation,” IEEE Trans. Med. Imag., vol. 18, no. 4, pp. 345– 354, 1999. [9] R. Shahidi, M. R. Bax, C. R. Maurer Jr, J. A. Johnson, E. P. Wilkinson, B. Wang, J. B. West, M. J. Citardi, K. H. Manwaring, and R. Khadem, “Implementation, calibration and accuracy testing of an image-enhanced endoscopy system,” IEEE Trans. Med. Imag., vol. 21, no. 12, pp. 1524– 1535, 2002.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

Original

Wang

HTRDC

HEwoG

HEwG

10

HT1D

HEPLP

Fig. 13. Comparison of distortion correction results on the s5c image by using different methods.

Fig. 14. Several endoscopic images are shown in the first row, and they are captured with different endoscopies and with different focal lengths. The corrected images are shown in the second row.

[10] D. Gonzalez-Aguilera, J. Gomez-Lahoz, and P. Rodr´ıguez-Gonz´alvez, “An automatic approach for radial lens distortion correction from a single image,” IEEE Sensors J., vol. 11, no. 4, pp. 956–965, 2011. [11] R. Miranda-Luna, C. Daul, W. C. Blondel, Y. Hernandez-Mier, D. Wolf, and F. Guillemin, “Mosaicing of bladder endoscopic image sequences: Distortion calibration and registration algorithm,” IEEE Trans. Biomed. Eng., vol. 55, no. 2, pp. 541–553, 2008. [12] S.-L. Chen, H.-Y. Huang, and C.-H. Luo, “Time multiplexed vlsi architecture for real-time barrel distortion correction in video-endoscopic images,” IEEE Trans. Circuits Syst. Video Technol., vol. 21, no. 11, pp. 1612–1621, 2011. [13] S.-H. Lee, S.-K. Lee, and J.-S. Choi, “Correction of radial distortion using a planar checkerboard pattern and its image,” IEEE Trans. Consum. Electron., vol. 55, no. 1, pp. 27–33, 2009. [14] R. Hartley and S. B. Kang, “Parameter-free radial distortion correction with center of distortion estimation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 29, no. 8, pp. 1309–1321, 2007. [15] J. Oh and K. Sohn, “Semiautomatic zoom lens calibration based on the camera’s rotation,” J. Electron. Imag., vol. 20, no. 2, pp. 023 006– 023 006, 2011. [16] D. Kim, H. Shin, J. Oh, and K. Sohn, “Automatic radial distortion correction in zoom lens video camera,” J. Electron. Imag., vol. 19, no. 4, pp. 043 010–043 010, 2010. [17] T.-Y. Lee, T.-S. Chang, S.-H. Lai, K.-C. Liu, and H.-S. Wu, “Wide-angle distortion correction by hough transform and gradient estimation,” in Proc. IEEE Visual Communications and Image Processing (VCIP’11), 2011, pp. 1–4. [18] J.-Y. Bouguet. (2004) Camera calibration toolbox for matlab. [Online]. Available: www.vision.caltech.edu/bouguetj/calib doc/ [19] M. Rufli, D. Scaramuzza, and R. Siegwart, “Automatic detection of checkerboards on blurred and distorted images,” in Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS’08), 2008, pp. 3121–3126. [20] S.-Y. Chen and W.-H. Tsai, “A systematic approach to analytic determination of camera parameters by line features,” Pattern Recogn., vol. 23, no. 8, pp. 859–877, 1990. [21] B. Prescott and G. McLean, “Line-based correction of radial lens distortion,” Graphical Models and Image Process., vol. 59, no. 1, pp. 39–47, 1997. [22] F. Devernay and O. D. Faugeras, “Automatic calibration and removal of distortion from scenes of structured environments,” in Proc. SPIE International Symposium on Optical Science, Engineering, and Instrumentation, 1995, pp. 62–72. [23] A. Wang, T. Qiu, and L. Shao, “A simple method of radial distortion

[24] [25] [26]

[27]

[28] [29] [30]

[31] [32] [33] [34] [35] [36] [37]

correction with centre of distortion estimation,” J. Math. Imag. and Vis., vol. 35, no. 3, pp. 165–172, 2009. P. V. Hough, “Method and means for recognizing complex patterns,” U.S. Patent 3,069,654, Dec. 18, 1962. K. Ehrenfried, “Processing calibration-grid images using the hough transformation,” Meas. Sci. Technol., vol. 13, no. 7, p. 975, 2002. R. T. Collins and J. R. Beveridge, “Matching perspective views of coplanar structures using projective unwarping and similarity matching,” in Proc. IEEE Conference on Computer Vision and Pattern Recognition (CVPR’93), 1993, pp. 240–245. T. Tuytelaars, L. Van Gool, M. Proesmans, and T. Moons, “The cascaded hough transform as an aid in aerial image interpretation,” in Proc. IEEE International Conference on Computer Vision (ICCV’98), 1998, pp. 67– 72. A. F. Habib and M. F. Morgan, “Automatic calibration of low-cost digital cameras,” Opt. Eng., vol. 42, no. 4, pp. 948–955, 2003. D. Xu, Y. F. Li, and M. Tan, “Method for calibrating cameras with large lens distortion,” Opt. Eng., vol. 45, no. 4, pp. 043 602–043 602, 2006. R. Cucchiara, C. Grana, A. Prati, and R. Vezzani, “A hough transformbased method for radial lens distortion correction,” in Proc. IEEE International Conference on Image Analysis and Processing, 2003, pp. 182–187. E. Rosten and R. Loveland, “Camera distortion self-calibration using the plumb-line constraint and minimal hough entropy,” Mach. Vision Appl., vol. 22, no. 1, pp. 77–85, 2011. D. Brown, “Decentric distortion of lenses,” Photogramm. Eng. Remote Sens., vol. 32, no. 3, pp. 444–462, 1966. G. P. Stein, “Internal camera calibration using rotation and geometric shapes,” Master’s thesis, Massachusetts Institute of Technology, Cambridge, 1993. R. Tsai, “A versatile camera calibration technique for high-accuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses,” IEEE J. Robot. Autom., vol. 3, no. 4, pp. 323–344, 1987. A. W. Fitzgibbon, “Simultaneous linear estimation of multiple view geometry and lens distortion,” in Proc. IEEE Conference on Computer Vision and Pattern Recognition (CVPR’01), vol. 1, 2001, pp. I–125. J. Kannala and S. S. Brandt, “A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 28, no. 8, pp. 1335–1340, 2006. S. Suzuki and K. Abe, “Topological structural analysis of digitized binary images by border following,” Computer Vision, Graphics, and Image Process., vol. 30, no. 1, pp. 32–46, 1985.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2013.2261816 TBME-01659-2012.R1. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. IEEE EARLY ACCESS ARTICLE.

Tung-Ying Lee (S’08-) received the B.S. and M.S. degrees from the Department of Computer Science, National Tsing Hua University, Taiwan, in 2005 and 2006, respectively. From 2005 to 2006, he was a member of Multimedia and Knowledge Engineering Laboratory, National Tsing Hua University. He had a research visit to Nara Institute of Science and Technology, Japan, during the summer of 2009. He is currently a Ph.D. candidate in Department of Computer Science at National Tsing Hua University and a member of Computer Vision Laboratory, National Tsing Hua University. His current research interests include medical image restortion, natural image restoration, stereoscopic image processing, computer vision and machine learning.

Tzu-Shan Chang received the B.S. degree from the Department of Information Management, National Chi Nan University, Taiwan, in 2010 and the M.S. degree from the Institute of Information Systems and Applications, National Tsing Hua University, Taiwan, in 2012. In 2012, she joined ASUSTeK Computer Inc., Taiwan.

Chen-Hao Wei was born in Kaohsiung, Taiwan in 1989. He received the B.S. degree from the Department of Computer Science, National Tsing Hua University, Taiwan, in 2011. He is currently a M.S candidate in Department of Computer Science at National Tsing Hua University. His research interests include computer vision and image processing.

11

Kai-Che Liu After receiving B.S., M.S., and Ph.D. from National Cheng Kung University in 2001, 2002 and 2005, respectively, he had been working at Industrial Technology Research Institute(ITRI), in the Electronic and Optoelectronics Research Laboratories, first as an engineer and then the assistant manager in the 3D system application division for the research of parallel processing on various processor architectures, and image-processing algorithms for 2D-to-3D conversion, image-based rendering and multi-view imaging. He had also worked on medical image processing and computer aided diagnosis. He has received the Outstanding Research Award four times, the Excellent Technical Paper Award twice and the Technology Promotion Award once in ITRI. After working in ITRI for several years, he joins the medical image research lab at Asian Institute of TeleSurgery (AITS/IRCAD-Taiwan) as the director in 2010. He has also served as deputy general manager of Ming Shi Co., Ltd., a biomedical company of ShowChwan healthcare system. Moreover, he had been the professor in IRCAD in 2011 and is the faculty in COmputational Surgery International NEtwork since 2012. Currently research interesting includes computer assisted surgery, computational surgery, 3D medical system, augmented reality application and surgical robot platform. He was the visiting scholar in CMU and Cornell University.

Hurng-Sheng Wu graduated from National Defense Medical Center, Taiwan in 1982. He had done his clinical research fellow in UCSF USA during 1987– 1988. He graduated from Tulane University with Master Degree of Public Health in 2006. Prof. Wu is the pioneer of robotic surgery in South-East Asia since 2002 and master in NOTES since 2007 in Taiwan. He was the President of the Taiwan Association of Endoscopic Surgery during 2008–2010 and current Dean of Asia Institute TeleSurgery since 2008. In addition, he has continued his research leading both clinical and Research & Development in advanced minimally invasive surgery and other aspects of advanced surgical science.

Shang-Hong Lai (M’95-) received the B.S. and M.S. degrees in electrical engineering from National Tsing Hua University, Hsinchu, Taiwan, and the Ph.D. degree in electrical and computer engineering from University of Florida, Gainesville, in 1986, 1988 and 1995, respectively. He joined Siemens Corporate Research in Princeton, New Jersey, as a member of technical staff in 1995. Since 1999, he returned to Taiwan as a faculty member in the Department of Computer Science, National Tsing Hua University. He is currently a professor in the same department and the director of the Computer and Communication Center in the university. In 2004, he was a visiting scholar with Princeton University. Dr. Lai’s research interests include computer vision, visual computing, pattern recognition, medical imaging, and multimedia signal processing. He has authored more than 200 papers published in the related international journals and conferences. Dr. Lai has been a member of program committee of several international conferences, including CVPR, ICCV, ECCV, ACCV, ICPR, PSIVT and ICME. He has been an associate editor for Journal of Signal Processing Systems since 2010. Moreover, he also served as a guest editor for special issues in Journal of Visual Communication and Image Representation as well as Journal of Signal Processing Systems.

Copyright (c) 2013 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing [email protected].