IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 12, DECEMBER 2010

2861

Reversible Projection Technique for Colon Unfolding Jianhua Yao*, Ananda S. Chowdhury, Member, IEEE, Javed Aman, and Ronald M. Summers

Abstract—Colon unfolding provides an efficient way to navigate the colon in computed tomographic colonography (CTC). Most existing unfolding techniques only compute forward projections. When radiologists find abnormalities or conduct measurements on the unfolded view (which is often quicker and easier), it is difficult to locate the corresponding region on the 3-D view for further examination (which is more accurate and reliable). To address this, we propose a reversible projection technique for colon unfolding. The method makes use of advanced algorithms including rotation-minimizing frames, recursive ring sets, mesh skinning, and cylindrical projection. Both forward and reverse mapping can be computed for points on the colon surface. Therefore, it allows for detecting and measuring polyps on the unfolded view and mapping them back to the 3-D surface. We generated realistic colon simulation data incorporating most colon characteristics, such as curved centerline, variable distention, haustral folds, teniae coli, and colonic polyps. Our method was tested on both simulated data and data from 110 clinical CTC studies. The results showed submillimeter accuracy in simulated data and –0.23 ± 1.67 mm in the polyp measurement using clinical CTC data. The major contributions of our technique are: 1) the use of a recursive ring set method to solve the centerline and surface correspondence problem; 2) reverse transformation from the unfolded view to the 3-D view; and 3) quantitative validation using a realistic colon simulation and clinical CTC polyp measurement. Index Terms—Biomedical image processing.

I. INTRODUCTION HE COLON is a convoluted tubular structure. The traditional way to navigate a colon surface along its centerline is not efficient due to limited view angle and the frequent obstructed view of polyps behind haustral folds. A more efficient way is to open and unfold the colon and view the flattened surface from above. Although this is not possible in optical colonoscopy, it can be done virtually in computed tomographic colonography (CTC) via computer graphics techniques. Several investigations [1]–[3] had showed that colon unfolding can in-

T

Manuscript received November 13, 2009; revised January 30, 2010, March 20, 2010, and May 10, 2010; accepted May 12, 2010. Date of publication June 10, 2010; date of current version November 17, 2010. This work was supported in part by the Intramural Research Program of the National Institutes of Health, Clinical Center. Asterisk indicates corresponding author. *J. Yao is with the National Institutes of Health, Bethesda, MD 20892 USA (e-mail: [email protected]). A. S. Chowdhury was with the National Institutes of Health, Bethesda, MD 20892 USA. He is now with the Department of Electronics and Telecommunication Engineering, Jadavpur University, Calcutta, West Bengal 70032, India (e-mail: [email protected]). J. Aman and R. M. Summers are with the National Institutes of Health, Bethesda, MD 20892 USA. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TBME.2010.2052255

crease the field of view and improve the sensitivity of polyp detection. Many techniques have been proposed to unfold the colon. They can be put in three categories. The first category is based on ray casting [4]. Vilanova and Groller [5] first locally unfolded the colon using a ray projection, and then, globally mapped the colon using a suitable parameterization. They divided the colon unfolding into three main steps: nonlinear ray casting to handle the multiple appearances of polyps, nonlinear 2-D scaling to reduce the distortion, and resampling to avoid missing polyps. Wang et al. [6] proposed a technique: using the electrical field of a charged centerline to transform the colon into a straight cylinder-like shape, they uniformly sample the planar cross sections orthogonal to the centerline, and then, applied a ray-tracing technique to generate the unfolded view. Sudarsky et al. [7] presented an efficient method based on skeletal subspace deformation and uniform ray casting orthogonal to the central path. The colon unfolding was conducted in two steps. First, the colon was straightened using mesh skinning, then the dissection was computed using the ray-casting technique. However, as some papers in this category [7] point out, it is difficult to choose the proper sampling rate in ray-casting-based techniques and structures behind the first object encountered by the rays may be obstructed. The second category is based on conformal-texture mapping [8]. Haker et al. [9] used an angle-preserving conformal mapping to map the entire colon surface onto a flat plane. They first topologically modeled the colon as an open-ended cylinder, and then, applied a finite-element method to construct a conformal map, which unfolded the colon surface onto a rectangular region. Hong et al. [10] proposed a conformal mapping based on minimizing harmonic energy to achieve angle preserving and minimum distortion. They first removed the topology noise of the surface and applied a holomorphic 1-form curve to cut open the surface. The surface was then conformal-mapped to a 2-D rectangular map. Huang et al. [11] adopted a similar conformal-mapping approach in their effort to visualize teniae coli on the unfolded view, where they demonstrated that the three-teniae coli are parallel and evenly distributed. This observation cannot be easily appreciated on the original 3-D colon view. The conformal-mapping-based techniques have to deal with texture distortion and surface parameterization. They often require high-quality surfaces and are computationally expensive. The third category is based on local projection. Paik et al. [2] proposed various map projection techniques, including cylindrical and planar projection for the fly-through of virtual colonoscopy. These increased the view angle up to 360◦ , and in some sense essentially computed local unfolded views. Vos et al. [12] projected six orthogonal images onto an unfolded cube to render the complete field of view. This was also a local

0018-9294/$26.00 © 2010 IEEE

2862

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 12, DECEMBER 2010

unfolding, but introduced extra distortion on the edge of the cube. Hoppe et al. [13] proposed a virtual colon dissection by placing a virtual camera inside the colon and perpendicular to the centerline. The camera displayed the inner colon surface with 90◦ field of view and rotated in 45◦ increments around the centerline. This technique had trouble-seeing polyps behind a fold, and introduced discontinuities when different view panels were stitched together. Most of the aforementioned colon-unfolding techniques only support forward mapping, i.e., only the projection from the 3-D colon surface to the 2-D unfolded plane is provided, the reverse is not. In the ray-casting-based techniques, multiple points may lie in the path of the same casting ray, causing ambiguity in the reverse projection. In most ray-casting techniques, a ray stops when it encounters the first object, making projection and backprojection unavailable for objects behind the first object. It is possible to circumvent this by extending the ray casting and using an additional data structure to store multiple hits by the same ray. The conformal-mapping technique requires a genus-one surface, which is not satisfied by most colon surfaces. Hong et al. [10] proposed a novel surface denoising and double covering method to convert a colon surface to a genusone surface. Although the conformal mapping is reversible if a bijective-mapping function is used, the denoising process is not reversible (or hard to reverse) due to the merging, breaking, and retriangulation of the surface. Since the denoising process changes the surface, the reverse mapping is imprecise. Without the reverse mapping, when radiologists find abnormalities or conduct measurement on the unfolded view (which is often quicker and easier than on the 3-D view), it is difficult to locate the corresponding region in the 3-D view for confirmation and further examination (which is more accurate and reliable). Since the polyps are somewhat distorted in the unfolded view, it is important to examine and measure the detections in both the 3-D and unfolded views. Since polyp detection and measurement is a major rationale for the use of CTC, the development of a quick, accurate means of matching the findings on unfolded and 3-D views would do much to improve the efficiency of interpretation of CTC. Another important issue is the precision of the reverse mapping. Although many commercially available software platforms for CTC interpretation provide an unfolded view, only a few can map points back and forth between the multiplanar reconstruction (MPR) or fly-through, and the unfolded view. Furthermore, most of the mappings are not precise and only project to the closest centerline point. Users have to navigate again in a small vicinity of a polyp to obtain an optimal view. Since polyps visualized in an unfolded view often appear distorted, the inability to verify and characterize them precisely in the fly-through view limits the use of the colon-unfolding views by radiologists in clinical practice. To address this problem, we developed an efficient colonunfolding technique, which also supports precise reverse projection, and assessed its performance on both simulated data and clinical CTC data. The results showed submillimeter accuracy in simulated data and in the polyp measurement using clinical CTC data. The novelties of our technique and the major contri-

butions are: 1) the use of a recursive ring set method to solve the centerline and surface correspondence problem (Section II-B); 2) reverse transformation from the unfolded view to the 3-D view (Section II-D); and 3) quantitative validation using both realistic colon simulation and clinical CTC polyp measurement (Section II-F and G). There are two ways to present the unfolded view of a colon surface. One is to generate a 2-D texture rendering of the unfolded view [7]–[10]. The other is to display a flat 3-D rendering of the mucosal surface [1]–[3], [5], [6], [12], [13]. The 3-D rendering techniques have the advantage to allow viewing the unfolded colon from different angles. Our technique belongs to the second category. II. METHODS Our method consists of the following sequence of steps. Given a CTC dataset, the 3-D colon surface is first segmented using intensity thresholding, region growing, and level sets [14]. The centerline of the colon is then extracted based on fast marching level set and topographical thinning [15], and rotationminimizing frames (RMF) are established along the centerline. After this, a recursive ring set technique is applied to map vertices on the colon surface to their corresponding centerline points. Next, mesh skinning is employed to straighten the colon. After this, cylindrical projection is applied to unfold the colon. Finally, reverse transformation is computed for every vertex. A. RMF The colon centerline is a 3-D spatial curve in the center of a colon, and is represented as a sequence of 3-D points from rectum to cecum. We first fit a fifth degree B-spline curve to the centerline. The B-spline curve utilizes the original centerline as control points and has explicit forms to compute tangent, normal, and binormal for every point along the curve [16]. We then discretize the centerline at 0.5 mm resolution (which is approximately the same resolution as the CT dataset). The Frenet–Serret frame (FSF) is a common tool to study the differential geometry of spatial curves [17]. The FSF of a given point c(u) on the curve is defined by its tangent t(u), normal n(u), and binormal b(u) vectors c(u) ˙ × c¨(u) ˙ t(u) = c(u) , b(u) = , n(u) = b(u) × t(u) c(u) ˙ c(u) ˙ × c¨(u) (1) where u is the curve variable, c(u) is one point on the curve, c(u) ˙ is the first degree derivative, and c¨(u) is the second degree derivative. However, FSF is not the best choice for our application, since it may exhibit strong rotational transition around the tangent vector [see Fig. 1(a)]. A so-called RMF is more appropriate [17]. RMF is formed by the tangent vector t(u) and two unit vectors f1 (u) and f2 (u) spanning the normal plane at c(u). f1 (u) and f2 (u) rotate as little as possible around t(u) between neighboring frames [see Fig. 1(b)]. They can be derived from the FSF f1 (u) = sin Ω(u)b(u) + cos Ω(u)n(u) f2 (u) = cos Ω(u)b(u) − sin Ω(u)n(u)

(2)

YAO et al.: REVERSIBLE PROJECTION TECHNIQUE FOR COLON UNFOLDING

2863

Fig. 1. FSF and RMF. (a) FSF; blue: b(u) and green: n (u). (b) RMF; blue: f1 (u) and green: f2 (u). Abrupt rotational transitions of normal (blue vector) and binormal (green vector) along the centerline are apparent in FSF. The transitions are smooth in RMF. Zoomed parts of the frames are shown.

where u is the curve variable and Ω is the angle difference between two neighboring frames  u τ (t) c(t) ˙ dt (3) Ω(u) = − 0

where τ = det{c, ˙ c¨, c¨˙}/c˙ × c¨2 is the torsion of the curve and t is the integral variable. The proof that RMF M (u) = {t(u), f1 (u), f2 (u)} is rotation minimizing along the centerline can be found in [17]. Fig. 1 illustrates the FSF and RMF of one centerline. In the following description, we use cj to represent the jth point on the centerline and Mj as its RMF. B. Recursive Ring Sets Once the RMF is established, every vertex on the colon surface must then be associated with its topologically corresponding centerline point and a local RMF. The methods in [5] and [7] rely on a distance map to establish the correspondence. A distance map stores the distances of vertices to their nearest centerline points and is often computed by an Euclidean distance transformation [18]. However, in 3-D colon, the distance map approach may fail to find the right correspondence in highly curved regions. For instance, in Fig. 2(a), c1 is the topological corresponding centerline point of v1 . However, c2 is closer to v1 than c1 in terms of Euclidean distance; therefore c2 is associated with v1 in terms of distance map. Since c1 and c2 have different RMFs, the miscorrespondence will cause incorrect transformation in the unfolding process. Another type of technique proposed in [19] is based on the colon cross section at each point along the centerline. They first applied a surface thinning to compute the skeletonization of the colon, and then, used an editing tool to prune the branches and loops to obtain the centerline. After this, cross sections orthogonal to the centerline were used to establish the centerline-to-surface correspondence. As pointed out in the paper, due to the tortuous nature of the colon, nearby cross sections may conflict with each other, resulting in missing or duplicate correspondences. Local iterative refinement [6] and sometimes manual intervention are needed to resolve the conflicts. To efficiently handle the correspondence problem, we propose a recursive ring set technique that can establish a topology-preserving correspondence. The concept of ring sets was also adopted in the centerline computa-

Fig. 2. Recursive ring sets. (a) Method scheme (solid line: colon wall, thick dashed line: centerline, and thin dashed line: ring sets). (b) Illustration of iterative ring set update. Thick dashed line is the centerline. Both red and green vertices are in the initial ring set. Red vertices remain in the initial ring set, while green vertices are put in the refinement set and eventually assigned to their corresponding ring sets. (c) Ring sets associated with the centerline (color strips). (d) Closeup view of a highly curved region.

tion in [20], where two parallel planes perpendicular to the local centerline segment were used to generate the ring sets. However, the method in [20] is computational expensive, since each vertex needs to be tested for every ring set. Our ring set technique is based on recursive refinement to obtain topology-preserving correspondence between the colon surface and its centerline. The recursive ring set technique is a four-step procedure. In the first step, for every surface vertex vi , the closest centerline point cj based on Euclidean distance is recorded as its initial corresponding centerline point, i.e., CL(vi ) = arg min (vi − cj )

(4)

cj

where CL(vi ) represents the corresponding centerline point for vi . A binary space partition (BSP) tree [21] of the centerline is created to accelerate the search. The searching time is therefore reduced from O(N ) to O(logN ), where N is the number of points on the centerline. In the second step, a ring set Rj is created for each centerline point cj . The initial set for Rj includes all vertices {vi }, whose closest centerline point is cj , i.e., ∀ vi , vi ∈ Rj , s.t. CL(vi ) = cj .

(5)

In the third step, the ring sets are recursively updated along the centerline. A connected component analysis is conducted for each ring set Rj . Two vertices are considered connected if there is an edge joining them. The largest connected surface patch remains in the ring set, and the rest are removed from Rj and put in a refinement set Q. In the fourth step, for every point vk in set Q, its closest neighbor vm not in set Q is retrieved, i.e., vm = arg min (vk − vn ) vn , vn ∈ /Q

(6)

2864

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 12, DECEMBER 2010

and vm ’s corresponding centerline point is assigned to vk , i.e., CL(vk ) = CL(vm ). The process is repeated until Q is empty. Steps 3 and 4 are repeated until the ring set for every centerline point has one and only one connected component, i.e., ∀ vi , vk ∈ Rj , there is a path connecting vi and vk in Rj . The scheme is illustrated in Fig. 2(a). Here v1 is in the ring set of c1 , so CL(v1 ) = c1 . Fig. 2(b) illustrates the iterative ring set update (steps 3 and 4). Fig. 2(c) shows the recursive ring sets associated with the centerline (each color strip indicates one ring set), and Fig. 2(d) shows a closeup view in a highly curved region. After the ring sets are defined, the radius of a ring set is defined as the average distance between every point in the ring set and its corresponding centerline point.

Fig. 3. Colon straightening using mesh skinning. (a) Original 3-D colon surface. (b) Straightened colon. (c) Closeup view of the straightened colon. The straightening method is based on mesh skinning and RMF, which transform the colon centerline into a straight line.

C. Mesh Skinning Mesh skinning [22] is a skeleton-driven deformation technique widely used in computer animation. We apply this technique to straighten a colon using its centerline as the skeleton and its surface as the skin. Here, we use a similar technique as the one in [7]. This is a two-step procedure. In the first step, the centerline is deformed to form a straight line along the z-axis. The straightened centerline is represented as follows: cj .x = cj .y = 0 cj .z = cj −1 .z + cj − cj −1 

(7)

where cj and cj are the coordinates in the original and straightened centerlines, respectively, and (x, y, z) are the 3-D coordinates. In the second step, the vertices on the colon surface are transformed using the RMF associated with their corresponding centerline point CL(vi ). The transformation is written as follows: vi = cj + Mj−1 (vi − cj )

(8)

where cj = CL(vi ) is the corresponding centerline point of vi , and Mj is the RMF of cj . In order to achieve a smooth transformation, we adopt the strategy in [7] to include a neighborhood of cj in the calculation. The new transformation is then updated as follows: vi =

k =j +δ 

  wi,k ck + Mk−1 (vi − ck )

(9)

k =j −δ

where δ is the neighborhood range, wi,k is the weight, which is defined to be inversely proportional to the distances di,k between vertex vi and centerline point ck , i.e., d−1 i,k wi,k = j =k +δ j =k −δ

d−1 i,k

di,k = vi − ck .

(10)

The neighborhood range is set empirically at 3. We have noted that the bigger the neighborhood, the smoother the transition and the straightened colon. However, a larger neighborhood may also create bigger distortion. Fig. 3 shows the straightened

colon using mesh skinning and one closeup view in the middle of the colon. D. Reversible Projection and Colon Unfolding Cylindrical projection is widely used in geographical maps to unfold the earth onto a 2-D map [23]. In the cylindrical map projection, the meridians are mapped to equally spaced vertical lines and the circles of latitude are mapped to horizontal lines. One trait of the cylindrical projection is its ability to preserve distance on the surface. We can in a similar fashion unfold the straightened colon using cylindrical projection and the colon centerline as the analog of the earth’s axis of rotation in the geographical map. The cylindrical projection maps every point on the straightened colon onto a plane at y = 0 to get the x, z coordinates. We preserve the distance-to-centerline function for every vertex by mapping it as the height elevation above the plane (y coordinate). The projection and following unfolding are formulated as follows: v˜i .z = vi .z  v˜i .y = (vi .x)2 + (vi .y)2 v˜i .x = arctan(vi .y, vi .x) × v˜i .y

(11)

where v˜i is the coordinate of vertex vi on the unfolded colon, v˜i .z is the distance along the centerline, v˜i .y is the distance to the centerline, and v˜i .x is the arc length around the centerline. We cut open the colon at 0◦ meridian during the unfolding process. An edge (vi , vj ) is deleted from the unfolded colon surface if it crosses the 0◦ meridian, i.e., arctan(vi .y, vi .x) × arctan(vj .y, vj .x) < 0. All other edges are retained. Therefore, the unfolded colon is essentially a cut-opened 3-D surface. Since there is a one-to-one correspondence between the original 3-D surface and the unfolded one, we preserve not only the topology of the original surface, but also the vertex resolution. Fig. 4 shows the unfolded colon and two closeup views, where the colon is totally collapsed, there is no colon surface and the ring sets will be empty. Our recursive ring set technique can handle this situation to get a correct unfolded surface. The projection from the 3-D colon to the unfolded colon is reversible. The reverse projection of (11) can be computed

YAO et al.: REVERSIBLE PROJECTION TECHNIQUE FOR COLON UNFOLDING

2865

Fig. 4. Unfolded colon. (a) Entire unfolded colon. (b) and (c) Closeup views. The unfolding method is based on the recursive ring sets and cylindrical projection. One section near the rectum is collapsed (arrow). At the collapsed colon segments, the ring set associated with the centerline is empty and the distention value is 0.

in two steps [see (12) and (13)]. Given a coordinate v˜i in the unfolded colon, the coordinate in the straightened colon vi can be obtained as follows:  v˜i .x vi .x = v˜i .y × cos v˜i .y  v˜i .x  vi .y = v˜i .y × sin v˜i .y vi .z = v˜i .z.

(12)

We can compute the reverse transformation from the straightened colon to the original 3-D colon using the following equation: vi = cj + Mj (vi − cj ).

(13)

Combining (12) and (13), a point v˜i on the unfolded colon can be reversely projected back to a point vi on the original 3-D colon. Since (9) smoothes the unfolded colon, the reverse projection is not exact. The relationship between the neighborhood range in (9) (controlling the smoothness) and the accuracy of the reverse projection will be reported in Section III. To get an exact reverse projection, (8) instead of (9) should be used. In our implementation, we keep the coordinates for both the unsmoothed unfolding result (8) and smoothed unfolding result (9). This reverse transformation allows a user to locate and measure lesions on the unfolded colon, and then, map them back to the 3-D colon for confirmation and further characterization. This provides an easier and quicker way than performing the measurement task directly on the 3-D colon. E. Algorithm Complexity Our technique is computationally efficient. We conducted a complexity analysis by breaking down the three main parts of the algorithm: recursive ring sets, mesh skinning, and reversible projection. In the following analysis, N is the number of centerline points and V is the number of colon vertices, where V is much greater than N . In the recursive ring sets stage, step 1 requires every colon vertex to search for its closest centerline point, which makes its

Fig. 5. Colon simulation. (a) Digital phantom of a cylindrical tube with the haustral folds and teniae coli. (b) Surface of simulated colon. (c) Surface of realcolon CTC. The digital phantom (see Fig. 6b) was generated using the centerline and distention values obtained from the real colon (see Fig. 6c). The simulated and real colons highly resemble each other. The colon simulation incorporates most colon characteristics in one model, including curved centerline, variable distention, haustral folds, teniae coli, and colonic polyps.

time complexity O(V logN ). Step 2 needs to locate the ring set for every vertex (O(V )). Step 3 involves a connected component analysis for every ring set. The complexity is O(N × K logK), where K is the average number of vertices in a ring set and K = V /N . The complexity in step 3 can be further simplified as O(N × K logK) = O(N × V /N log V /N ) = O(V log (V /N )). Step 4 is the rearrangement of vertices in set Q (O(Q)), where Q  V . Therefore, the complexity for the recursive ring sets is O(V logN ) + O(V ) + O(V log(V /N ) + O(Q) = O(V logN ). In the mesh-skinning stage, step 1 is the straightening of the colon centerline (O(N )). Step 2 is the mesh skinning of colon vertices (O(V )). The complexity for mesh skinning is then O(N ) + O(V ) = O(V ). In the reversible projection stage, every vertex on the colon is projected to a flat plane (O(V )). Therefore, the overall complexity of the entire algorithm is as follows: O(V log N ) + O(V ) + O(V ) = O(V log N ).

(14)

F. Colon Simulation and CTC Data Our colon unfolding technique was evaluated using both simulated colon data and clinical CTC data. We created a realistic colon simulation in three steps. First, a digital phantom of a cylindrical tube model with haustral folds and teniae coli was built [see Fig. 5(a)]. Haustral folds were added in the form of indentations onto the tube. The folds were modeled as a series of trisected torus with elliptic cross sections. The trisections were partitioned equally in circumferential locations 0, 120, and 240◦ on the cross-section plane. Spaces were placed at the trisection to simulate the teniae coli. Second, the centerline from a clinical CTC was computed and employed. Every point on the centerline was equipped with a RMF (Section II-A) and its associated distention value computed from the CTC data. The

2866

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 12, DECEMBER 2010

radius of the ring set was used as the distention value at the corresponding point on the centerline. Third, inverse mesh skinning (Section II-C) was applied to warp the tube phantom around the centerline to generate a realistic simulation of clinical CTC data. Simulated colonic polyps, modeled as semiellipsoids, were also added to the phantom. Fig. 5(b) shows the surface of a colon simulation. Fig. 5(c) shows the surface of a real colon, whose centerline and distention data was employed to generate the simulation in Fig. 5(b). The simulated colon highly resembles the real one. To our best knowledge, this is the first colon simulation that incorporates most colon characteristics in one model, including curved centerline, variable distention, haustral folds, teniae coli, and colonic polyps. The clinical CTC data were acquired at three institutions (National Naval Medical Center Bethesda, Naval Medical Center San Diego, and Walter Reed Army Medical Center) as follows. Patients underwent 24-h colonic preparation that consisted of oral administration of a cathartic and barium- and iodinecontaining contrast agents. Each patient was scanned supine and prone using a four- or eight-channel CT scanner (GE LightSpeed or LightSpeed Ultra). CT scanning parameters included 1.25- to 2.5-mm section collimation, 15 mm/s table speed, 1-mm reconstruction interval, 100 mA·s, and 120 kVp. Each dataset had more than 400 slices on average. To evaluate the quality of the unfolded colon, an experienced user visually examined every case. An unfolded colon was judged as satisfactory if the colon features, such as haustral folds and polyps were easily identified. We also recorded the computation time for the unfolding process.

Fig. 6. Gallery of polyps from 3-D and unfolded views. (a) Sessile adenoma polyp of 6 mm in sigmoid colon. (b) Sessile adenoma polyp of 9 mm in transverse colon. (c) Round adenoma of 6 mm in descending colon. (d) Round hyperplastic polyp of 7 mm in sigmoid colon. Pictures on left: from the 3-D views and pictures on right: from the unfolded views. Arrows point to the polyps.

G. Polyp Size Measurement on Unfolded Colon Estimating the size of a polyp is essential in CTC diagnosis, since large polyps are more likely to become malignant and cancerous. To explore the reverse transformation of our technique, we compared the polyp measurements made on both the 3-D and the unfolded views. Here, the polyp size was defined as the largest diameter across the polyp determined by the operator. In this experiment, the operator first conducted the measurement on the unfolded surface and recorded the two endpoints. The two endpoints were then mapped to the 3-D surface for additional measurement. Fig. 6 illustrates a gallery of polyps on both the 3-D view and the unfolded view, including polyps with different sizes and in different segments of the colon.

Fig. 7. Unfolded view of the simulated colon and ground truth of the unfolded view. The colon simulation is shown in Fig. 7(b). Haustral folds, teniae coli, and polyps are simulated. (a) Unfolded view of the simulated colon. (b) Ground truth of the unfolded view. TABLE I FOLD COMPARISON BETWEEN UNFOLDED VIEW (UV) AND GROUND TRUTH (GT) FOR THE SIMULATED COLONS (N = 49)

III. RESULTS We ran tests on the simulated colons to validate the capability of our unfolding technique in preserving the location, orientation, and size information. Haustral folds were used as the landmarks. Fig. 7 shows the unfolded view of the colon simulation and the ground truth of the unfolded view. We located the two endpoints of the middle folds on both the unfolded view and the ground truth view. We then compared the difference of the fold location (midpoint of two endpoints), orientation (direction of two endpoints), and size (distance between two endpoints). Table I summarizes the statistics of the comparison

for 49 folds. We were able to achieve subvoxel accuracy (voxel size is 1 mm3 ) Our method was also tested on 110 clinical CTC studies. These cases were randomly selected from a cohort of 1186 patients screened at the three medical institutes. We only select the cases with successfully extracted centerlines. The result of colon unfolding was visually evaluated. An unfolding was determined successful if the following three criteria were satisfied:

YAO et al.: REVERSIBLE PROJECTION TECHNIQUE FOR COLON UNFOLDING

2867

Fig. 10. Reverse projection accuracy. The default neighborhood range is 3, which gives 0.30 ± 0.17 mm reverse-projection error. Fig. 8. Polyp linear size measurement on 3-D and unfolded views. The two sets of measurements show high degree of correlation. The measurements on unfolded view tend to be larger.

Fig. 9. Bland–Altman plot of agreement between polyp measurements on 3-D and unfolded views. The mean percentage difference is 3.6% and the 95% limits of agreement is [−46.8%, 53.9%].

1) the haustral folds could be identified and were oriented roughly perpendicular to the colon centerline; 2) the teniae coli could be identified; and 3) polyps indentified in the 3-D view could be identified in the unfolded view. The evaluation was conducted by an experienced research scientist and verified by an experienced radiologist. All unfolded colons were determined visually acceptable. The average computation time on a 3.0 GHz Pentium processor with 4 GB RAM was 14.9 ± 3.8 s. Currently, our algorithm is not optimized for speed. We conducted linear size measurements of 103 colonic polyps on both the 3-D colon surface and the unfolded surface. Fig. 8 shows the scatter plot of the two sets of measurements, which show a high degree of correlation. The equation for the linear regression is y = 1.14x – 0.54, which indicates the unfolded view tends to give larger measurements when the colon is straightened. Fig. 9 shows the agreement between the two sets of measurements using Bland–Altman plots [24]. The mean percentage difference is –3.6% and the 95% limits of agreement is [−46.8%, 53.9%]. Our reversible projection provides an efficient and reliable polyp measurement tool. The operator can conduct the measurement on the unfolded view (which is often

quicker and easier to do), and then, reversely project the endpoints to the 3-D view for verification and generation of the 3-D measurement. In the experiment to evaluate the accuracy of the reverse projection, we first projected the linear size measurement (two endpoints) on the original surface to the unfolded surface, and then, reverse projected the two endpoints back to the original surface. The distances between the original endpoints and the reverse-projected endpoints were used to evaluate the accuracy of the reverse projection. The accuracy was affected by the neighborhood range (δ) used in (9) to smooth the unfolded view. When δ is 0 [i.e., (8)], we will get an exact reverse projection. The larger δ is, the smoother the unfolded view and the less accurate the reverse projection. This relationship is plotted in Fig. 10. The default δ is 3, which gives 0.30 ± 0.17 mm reverseprojection error. IV. DISCUSSION AND CONCLUSION Both our method and the one described in [7] used the meshskinning technique to straighten the colon and were much more efficient than other prevalent methods. However, our method includes several improvements over the one discussed in [7]. First and foremost, we developed a novel recursive ring set technique to find the topology-preserving correspondence between the centerline and vertices on the colon surface. It is more robust and does not rely on the distance map computed from the original CT data. Second, we do not use ray casting. We use a cylindrical projection instead to map the vertices to the RMF along the centerline. This preserves the resolution and topology of the 3-D surface and avoids the obstruction problem. Third, our method supports reverse transformation from the unfolded colon to the original 3-D colon. This allows us to conduct the further characterization and measurement of polyps detected on the unfolded view. Finally, our technique maps each vertex onto the unfolded plane instead of sampling the centerline and surface. This avoids the undersampling (missing polyps) and oversampling (double polyps) problems that occur in some other methods [5]. Though Paik et al. [2] adapted a similar cylindrical projection in their virtual colon fly-through to increase the view angle; the technique was limited to the local

2868

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 12, DECEMBER 2010

view field and their projection was not reversible. Another approach to achieve the reverse projection is to store the reverse coordinates for vertices on the unfolded surface and interpolate the coordinates for all other points. This approach still needs to go through all the difficult steps (ring sets, mesh skinning, and cylindrical projection) to obtain the reverse coordinates. As for our method, the reverse projection has an explicit formula, so it only takes an insignificant amount of the running time in the entire process. Compared to the “storage” approach, it has almost the same running cost, but with the following benefits: 1) no need for interpolation, therefore more precise; 2) explicit formula for every point on the surface, not just for the vertices; and 3) less storage requirement. The capability to precisely transfer from the fly-through to the unfolded view and vice versa can add important practical value and convenience for the application of CTC. Another desirable function is to locate the polyp on the unfolded surface by clicking on the MPR or fly-through views. Our reversible colonunfolding technique provides a way to obtain exact reverse projection between the unfolded and original representations of the colon surface. This technique could potentially prove to be useful in clinical practice. Furthermore, exact reverse mapping provides a more convenient way for polyp measurement. Since it is easier to pick endpoints in the unfolded view, the operator can conduct approximate measurement in the unfolded view, and then, reversely map the endpoints to 3-D view for exact measurement. It has been shown in [25] that the inter and intraoperator variability can be big for linear size measurements using 3-D views. One possible reason is that different view angles during 3-D navigation give different perceptions of the polyp, which leads the operator to make different measurements. The unfolded view provides a more consistent and larger view of the polyp and its surrounding region. Therefore, it alleviates the variation in measurement. In conclusion, we presented a computationally efficient reversible projection technique for colon unfolding based on recursive ring set formulation, mesh skinning, and cylindrical projection. The earlier techniques have their origins in diverse areas like computer animation and geographical mapping. Among all these techniques, the recursive ring set is a totally new method to solve the correspondence problem. RMF, mesh skinnings, and cylindrical projections had been employed in prior work. However, we are the first to incorporate these into a pipeline and derive an explicit formulation to obtain the reverse projection. The proposed method is quite robust and works for all colon segments including poorly distended ones. The method has been tested on realistic colon simulation and 110 clinical datasets. It achieved submillimeter level accuracy in the simulation data and generated satisfying results in the clinical data. The technique supports polyp measurement on the unfolded view, an essential feature for virtual colonoscopy. For the future, we are currently investigating effective ways to integrate our technique into colon navigation and computeraided detection, and planning to apply the realistic colon simulation in other research areas, such as supine and prone colon registration, teniae coli detection, and colon segmentation.

ACKNOWLEDGMENT The authors would like to thank P. J. Pickhardt, W. R. Schindler, and R. Choi for providing CTC and supporting data. They would also like to thank Dr. A. Dwyer for providing the insightful comments and proofreading. REFERENCES [1] K. T. Johnson, C. D. Johnson, J. G. Fletcher, R. L. MacCarty, and R. M. Summers, “CT colonography using 360-degree virtual dissection: A feasibility study,” AJR Amer. J. Roentgenol., vol. 186, pp. 90–95, 2006. [2] D. S. Paik, C. Beaulieu, R. B. Jeffrey, C. Karadi, and S. Napel, “Visualization modes for CT colonography using cylindrical and planar map projections,” J. Comput. Assist. Tomogr., vol. 24, pp. 179–188, 2000. [3] A. C. Silva, C. V. Wellnitz, and A. K. Hara, “Three-dimensional virtual dissection at CT colonography: Unraveling the colon to search for lesions,” Radiographics, vol. 26, pp. 1669–1686, 2006. [4] H. Ray, D. Silver, H. Pfister, and T. Cook, “Ray casting architectures for volume visualization,” IEEE Trans. Vis. Comput. Graph., vol. 5, no. 3, pp. 210–223, Jul. 1999. [5] A. Vilanova and E. Groller, “Geometric modeling for virtual colon unfolding,” in Geometric Modeling for Scientific Visualization, G. Brunnett, B. Hamann, H. Muller, and L. Lisen, Eds. Basel, Switzerland: Birkhauser, 2004, pp. 453–488. [6] G. Wang, G. McFarland, B. P. Brown, and M. W. Vannier, “GI tract unraveling with curved cross sections,” IEEE Trans. Med. Imag., vol. 17, no. 2, pp. 318–322, Apr. 1998. [7] S. Sudarsky, B. Geiger, C. Chefd’hotel, and L. Guendel, “Colon unfolding via skeletal subspace deformation,” in Proc. 11th Int. Conf. Med. Image Comput. Comput.-Assist. Interv., Part II (Lecture Notes in Computer Science, vol. 5242), D. Metaxas, L. Axel, G. Fichtinger, and G. Sz´ekely, Eds. New York, Sep. 6–10, 2008, pp. 205–212. [8] S. Haker, S. Angenent, A. Tannenbaurn, R. Kikinis, G. Sapiro, and M. Halle, “Conformal surface parameterization for texture mapping,” IEEE Trans. Vis. Comput. Graph., vol. 6, no. 2, pp. 181–189, Apr.–Jun. 2000. [9] S. Haker, S. Angenent, A. Tannenbaurn, and R. Kikinis, “Nondistorting flattening maps and the 3-D visualization of colon CT images,” IEEE Trans. Med. Imag., vol. 19, no. 7, pp. 665–670, Jul. 2000. [10] W. Hong, X. Gu, F. Qiu, M. Jin, and A. Kaufman, “Conformal virtual colon flattening,” in Proc. SPM, 2006, pp. 85–93. [11] A. Huang, D. Roy, and R. M. Summers, “Teniae coli guided navigation and registration for virtual colonoscopy,” presented at the IEEE Vis., Baltimore, MD, 2005. [12] F. M. Vos, I. Serlie, R. van Gelder, F. H. Post, R. Truyen, F. A. Gerritsen, J. Stoker, and A. Vossepoel, “A new visualization method for virtual colonoscopy,” in Proc. 4th Int. Conf. Med. Image Comput. Comput.-Assist. Interv. (Lecture Notes in Computer Science, vol. 2208), W. J. Niessen and M. A. Viergever, Eds. London, U.K.: Springer-Verlag, Oct. 14–17, 2001, pp. 645–654. [13] H. Hoppe, C. Quattropani, A. Spreng, J. Mattich, P. Netzer, and H.-P. Dinkel, “Virtual colon dissection with CT colonography compared with axial interpretation and conventional colonoscopy: Preliminary results,” AJR Amer. J. Roentgenol., vol. 182, pp. 1151–1158, 2004. [14] M. Franaszek, R. M. Summers, P. J. Pickhardt, and J. R. Choi, “Hybrid segmentation of colon filled with air and opacified fluid for CT colonography,” IEEE Trans. Med. Imag., vol. 25, no. 3, pp. 358–368, Mar. 2006. [15] R. L. Van Uitert and R. M. Summers, “Automatic correction of level set based subvoxel precise centerlines for virtual colonoscopy using the colon outer wall,” IEEE Trans. Med. Imag., vol. 26, no. 8, pp. 1069–1078, Aug. 2007. [16] E. W. Weisstein, CRC Concise Encyclopedia of Mathematics, 2nd ed. London, U.K./Boca Raton, FL: Chapman and Hall/CRC Press, 1999. [17] C. Maurer and B. Juttler, “Rational approximation of rotation minimizing frames using Pythagorean-hodograph cubics,” J. Geom. Graph., vol. 3, pp. 141–159, 1999. [18] O. Cuisenaire and B. Macq, “Fast euclidean distance transformation by propagation using multiple neighborhoods,” Comput. Vis. Image Understanding, vol. 76, pp. 163–172, 1999. [19] E. Sorantin, G. Werkgartner, E. Balogh, A. Bartroli, K. Palagyi, L. G. Nyul, and L. Rusko, “Virtual dissection and automated polyp detection of the

YAO et al.: REVERSIBLE PROJECTION TECHNIQUE FOR COLON UNFOLDING

[20] [21] [22] [23] [24] [25]

colon based on spiral CT—Techniques and preliminary experience on a cadaveric phantom,” Eur. Surg., vol. 34, pp. 143–149, 2002. G. Iodanescu and R. M. Summers, “Automated centerline for computed tomography colonography,” Acad. Radiol., vol. 10, pp. 1291–1301, 2003. H. Fuchs, Z. M. Kedem, and B. F. Naylor, “On visible surface generation by a priori tree structures,” Comut. Graph., vol. 14, pp. 124–133, 1980. J. P. Lewis, M. Cordner, and N. Fong, “Pose space deformation: A unified approach to shape interpolation and skeleton-driven deformation,” in Proc. ACM SIGGRAPH, 2000, pp. 165–172. J. P. Snyder, Map Projections—A Working Manual. Reston, VA: USGS, 1987. J. M. Bland and D. G. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” Lancet, vol. i, pp. 307–310, 1986. S. Yeshwant, R. M. Summers, J. Yao, D. Brickman, J. R. Choi, and P. J. Pickhardt, “Polyps: Linear and volumetric measurement at CT colonography,” Radiology, vol. 241, pp. 802–811, 2006.

Jianhua Yao received the B.S. degree from Tianjin University, Tianjin, China, the M.S. degree from Tsinghua University, Beijing, China, and the Ph.D. degree from Johns Hopkins University, Baltimore, MD, all in computer science. In 2002, he joined as a Staff Scientist the National Institutes of Health, Bethesda, MD, where he directs a Clinical Image-Processing Laboratory and also engaged with the Imaging Biomarker and ComputerAided Diagnosis (CAD) Laboratory. He has authored or coauthored more than 100 papers in journals and conference proceedings. He holds two patents in colon cancer CAD technique. His current research interests include clinical image processing, deformable model, nonrigid registration, CAD, and computed tomographic colonography.

Ananda S. Chowdhury (M’02) received the B.Sc. (Hons.) in physics from Presidency College, Calcutta, India, in 1996, the B.Tech. degree in electronics engineering from the Institute of Radiophysics and Electronics, Calcutta, in 1999, the M.E. degree in computer engineering from Jadavpur University, Calcutta, in 2001, and the Ph.D. degree in computer science from the University of Georgia, Athens, GA, in 2007. From August 2007 to December 2008, he was a Postdoctoral Fellow in the Department of Radiology and Imaging Sciences, National Institutes of Health, Bethesda, MD. He is currently a Reader in the Department of Electronics and Telecommunication Engineering, Jadavpur University, Calcutta. His research interests include computer vision, pattern recognition, and image processing. Dr. Chowdhury is a member of the IEEE Computer Society.

2869

Javed Aman received the B.S. degree in biomedical engineering from Johns Hopkins University, Baltimore, MD. He is currently a Postbaccalaureate Intramural Research Training Award Fellow at the National Institutes of Health, Bethesda, MD. His research interests include topics related to computer vision (contentbased image retrieval and computer-aided detection).

Ronald M. Summers received the B.A. degree in physics, and the M.D. and Ph.D. degrees in medicine/anatomy and cell biology from the University of Pennsylvania, Philadelphia, PA. He was with the Presbyterian, University of Pennsylvania Hospital, Philadelphia, PA, for medical internship and also with the University of Michigan, Ann Arbor, MI, for radiology residency. He was an MRI Fellow at Duke University, Durham, NC. In 1994, he joined the Diagnostic Radiology Department, National Institutes of Health Clinical Center, Bethesda, MD, where he is currently a Tenured Senior Investigator and Staff Radiologist. He is currently the Chief of the Clinical Image-Processing Service and directs the Imaging Biomarkers and Computer-Aided Diagnosis (CAD) Laboratory. He has authored or coauthored more than 200 journals, review, and conference proceedings articles. He holds nine patents. His research interests include virtual colonoscopy, CAD, multiorgan multiatlas registration, development of large radiologic image databases, thoracic and gastrointestinal radiology, and body cross-sectional imaging. Dr. Summers was the recipient of the Presidential Early Career Award for Scientists and Engineers in 2000. He is a Co-Chair of the CAD Program of the Annual International Society for Optical Engineering Medical Imaging conference.