Surface Reconstruction from Medical Imaging for Use in a Computer-AidedDesign (CAD) Environment
A thesis by Gregory M. Sturgeon M.S.E. Presented January 2005
In partial fulfillment of the requirements for the Masters of Science in Engineering degree

Padnos School of Engineering Grand Valley State University Advisor: Dr. Jeff Ray

Grand Valley State University Padnos School of Engineering

Certificate of Examination
Advisor ________________________ Dr. J. Ray Examining Board ________________________ Dr. H. Jack ________________________ Dr. J. Ray ________________________ Dr. J. Farris The EGR 693 thesis by Gregory M. Sturgeon entitled Surface Reconstruction from Medical Imaging for Use in a Computer-Aided-Design (CAD) Environment is accepted in partial fulfillment of the requirements for the degree of Masters of Science in Engineering

________________________ Date

________________________ Dr. H. Jack, Graduate Chair

ii

Abstract
A method was developed to create Computer-Aided-Design (CAD) models for bones of the human body utilizing medical imaging data. The human hand was chosen as the subject of the research. Computed Tomography (CT) imaging was chosen to provide a volumetric data set. This data set was visualized through an isosurfacing technique utilizing the marching cubes algorithm. The original CT data set contained slices that were not aligned with the natural orientation or long axis of the bones. Transformation matrices and linear interpolations were used to generate a data set of slices oriented along the natural axis of the bones. Contours were created on these slices through an edge-tracking method. B-Spline curves were then constructed utilizing the contour’s vertices as knot points. A consistent starting location was found on each closed B-Spline curve relative to its centroid. Points on the closed B-Spline curves were then selected to define open non-uniform B-Spline curves in Pro/Engineer, a common CAD software package. Pro/Engineer was then utilized to trim the B-Spline curves to obtain their uniform portions. Cross-curves were developed through groupings of parallel B-Spline curves in order to define a closed boundary for a Boundary Blend surface patch. These surface patches were joined to adjacent surface patches and maintained C2 curvature continuity. The method presented applied common visualization techniques to a data set from CT imaging. This provided vertices from which to construct curves and surfaces in a CAD environment resulting in the ability to create detailed anatomical CAD models.

iii

Table of Contents
LIST OF FIGURES ..........................................................................................................vii Chapter 1: Introduction and Background Information ..............................................1 1.1 Introduction ...........................................................................................................1 1.2 Medical Imaging Technologies ..............................................................................1 1.2.1 Computed Tomography ..................................................................................1 1.2.2 Magnetic Resonance Imaging .........................................................................1 1.2.3 Decision to Use CT Data.................................................................................2 1.3 CAD Background ..................................................................................................2 1.4 B-Spline Curve Background ..................................................................................2 1.4.1 B-Spline Curve Definition...............................................................................3 1.4.2 B-Spline Example ...........................................................................................4 1.4.3 Local Control of Change .................................................................................7 1.4.4 Uniform B-Spline ...........................................................................................8 1.4.5 Closed B-Spline ............................................................................................10 1.5 B-Spline Surfaces ................................................................................................10 1.5.1 Practical Consideration for Working With Surfaces ......................................12 1.6 Overview of Visualization Techniques.................................................................13 1.6.1 Color Mapping..............................................................................................13 1.6.2 Contouring ....................................................................................................14 1.6.3 Contouring Methods .....................................................................................14 1.6.3.1 Edge Tracking ........................................................................................15 1.6.3.2 Marching Squares11 ................................................................................15 1.6.4 Extending Contouring to 3-D Isosurfacing ....................................................17 1.6.5 Isosurfacing with the Marching Cubes11 ........................................................17 Chapter 2: Alternative Approaches............................................................................20 2.1 SURFdriver .........................................................................................................20 2.1.1 Using SURFdriver ........................................................................................21 2.1.2 Difficulty of Sketching Contours...................................................................21 2.1.3 Resulting Surface from SURFdriver..............................................................22 2.1.4 Importing into Pro/Engineer..........................................................................22 2.1.5 Tracing the Iges Curves.................................................................................23 2.1.6 Surfaces from Pro/Engineer ..........................................................................23 2.1.7 Discussion of Results ....................................................................................24 2.2 Isosurfacing in MATLAB12 .................................................................................24 2.2.1 Importing the CT data ...................................................................................24 2.2.2 Isosurface Function.......................................................................................26 2.2.3 Reducing the Noise .......................................................................................27 2.2.4 Discussion of Results ....................................................................................28 2.3 Working with the Isosurfaces...............................................................................28 2.3.1 Method of Regrouping the Faces and Vertices ..............................................29 2.3.2 Discovery of Joined Bones............................................................................29 2.3.3 Slice Plane Through Joints ............................................................................29

iv

2.3.4 Isosurfacing at Various Values......................................................................30 2.3.5 An Attempt to Separate the Bones.................................................................31 2.3.6 Discussion of Results ....................................................................................32 Chapter 3: Contours on Slice Planes to B-Spline Surfaces ........................................33 3.1 Reorienting the Slices of the CT Data ..................................................................33 3.1.1 Creating Points to Resample the CT Data......................................................34 3.1.2 Finding the Long Axis ..................................................................................35 3.1.3 Reorienting the Long Axis ............................................................................36 3.1.4 Creating the Uniform Grid ............................................................................38 3.1.5 Resampling the Uniform Grid .......................................................................39 3.2 Generating Contours and Fitting B-Spline Curves................................................39 3.2.1 Generating the Contours from MATLAB..........................................................39 3.2.2 Fitting B-Spline Curves to the Contour Vertices............................................39 3.2.3 Fitting B-Spline Curves to a Resampled Contour ..........................................41 3.2.4 Consistency in Starting Points.......................................................................44 3.2.5 Tracking the Contour Across Slices ..............................................................45 3.3 Utilizing Pro/Engineer to Create Surfaces from the B-Spline Curves ...................45 3.3.1 Importing to Pro/Engineer from an Ibl File Format .......................................46 3.3.1.1 Non-uniform B-Spline from an Ibl File ..................................................46 3.3.1.2 Trimming Curves to Maintain C2 Continuity ..........................................47 3.3.2 Exporting the B-Splines from MATLAB.......................................................48 3.3.3 Creating Surfaces in Pro/E ............................................................................49 3.3.4 Resulting Surfaces ........................................................................................52 3.3.5 Capping the Bones ........................................................................................52 3.3.5.1 Style Curves to Cap the Bones................................................................52 3.3.5.2 Capping the Finger .................................................................................54 Chapter 4: Summary and Conclusions.......................................................................58 Cited References:.........................................................................................................60 Appendix A: Non-uniform B-Spline Example............................................................62 Appendix B: import_ct_data.m ..................................................................................70 Appendix C: separateEntities.m .................................................................................71 Appendix D: separateBones.m...................................................................................73 Appendix E: TransformGridResample ......................................................................78 Appendix F: From Knot Points to Control Points .....................................................79 Appendix G: Fit a B-Spline Curve to Contour...........................................................81 G.1: fitBspline ...........................................................................................................81

v

G.2: knot_to_cntrl .....................................................................................................82 G.4: y_mod ...............................................................................................................83 G.5: N_array..............................................................................................................84 G.6: P_x_u_segments ................................................................................................85 G.7: P_y_u_segments ................................................................................................85 Appendix H: Resample the Contour at Equal Spaced Points....................................86 H.1: resampleContour................................................................................................86 H.2: contour_x...........................................................................................................87 H.3: contour_y...........................................................................................................87 H.4: contour_N_array................................................................................................88 H.5: contour_P_x_u_segments ..................................................................................89 H.6: contour_P_y_u_segments ..................................................................................89 H.7: contour_per_spaced ...........................................................................................90 H.7.1: contour_perimeter .......................................................................................92 H.8: contour_per_spaced_pnts...................................................................................93 Appendix I: Finding the Centroid of the Closed B-Spline Curve..............................94 I.1: centroid_x ...........................................................................................................94 I.2: centroid_y ...........................................................................................................95 I.3: findArea ..............................................................................................................96 Appendix J: Derivation of the Area, Centroid, and Perimeter for a Closed B-Spline ......................................................................................................................................97 Appendix K: P_x_cen_x ............................................................................................ 104 Appendix L: top_P_x_cen_x ..................................................................................... 106 Appendix M: Separating the Contours .................................................................... 107 M.1 contoursSeperate .............................................................................................. 107 M.2: per_spaced ...................................................................................................... 110 M.2.1: perimeter .................................................................................................. 112 M.3: per_spaced_pnts.............................................................................................. 113 Appendix N: Sample Ibl File..................................................................................... 114 Appendix O: Numbering Scheme of the Ibl Files..................................................... 115 Appendix P: Exporting to Excel ............................................................................... 116 P.1: ExportToExcel.m ............................................................................................. 116 P.2: xlscell16 ............................................................................................................ 118

vi

LIST OF FIGURES Chapter 1: Introduction and Background Information
Figure 1.1: Basis functions for a fourth order B-Spline with K = 7. …………………… 5 Figure 1.2: Basis functions (a) for a sample fourth order B-Spline (b) with n = 6. ……. 6 Figure 1.3: The local control of change varies depending on how many segments a given control point influences. …………………………………………………………………. 7 Figure 1.4: Uniform B-Spline basis functions with K = 4 shown over four intervals. …. 8 Figure 1.5: Comparison of uniform and non-uniform B-Spline curves. ……………….. 9 Figure 1.6: B-Spline surface with the control mesh and the numbering scheme of the control points. ………………………………………………………………………….. 11 Figure 1.7: The non-uniform B-Spline surface (a) interpolates the corner points of the control mesh where the uniform B-Spline surface (b) does not. ………………………. 12 Figure 1.8: Color mapped image of a typical slice of the CT data. …………………… 13 Figure 1.9: The color-mapped image (a) shows the region of interest of the CT data, 225x200 cells. The scalar values (b) are a 5x5 section as shown in the red on the colormapped image. …………………………………………………………………………. 14 Figure 1.10: Contour lines created from the vertices seen in Figure 1.9b. ……………. 15 Figure 1.11: Marching squares cases. …………………………………………………. 16 Figure 1.12: Solving the ambiguity by interpolating to find the above/below state of the center of the cell. ………………………………………………………………………. 16 Figure 1.13: The 15 marching cubes cases. …………………………………………… 17 Figure 1.14: Ambiguous cases of the marching cubes shown with their complementary cases. …………………………………………………………………………………… 18 Figure 1.15: The correct pairing of two adjoining marching cubes cases (a) will result in an isosurface with no holes. When the complementary case is incorrectly chosen (b) a hole will be introduced in the isosurface. ……………………………………………… 18

Chapter 2: Alternative Approaches
Figure 2.1: Typical images of two successive slices from the CT data set. …………... 20 Figure 2.2: A sketched contour shown with an imaged slice from the CT data set at two different magnifications.………………………………………………………………... 21 Figure 2.3: Surface reconstruction of the M-II bone from SURFdriver. ……………… 22 Figure 2.4: Boundary blend surfaces from Pro/Engineer shown with closed B-Spline curves sketched on parallel datum planes. ……………………………………………... 23 vii

Figure 2.5: Illustration of the multidimensional array storing the slices of the CT data. 25 Figure 2.6: The region of interest, enclosed in the box, was selected from the entire dataset. …………………………………………………………………………………. 26 Figure 2.7: Isosurface image corresponding to the boundary between bone and soft tissue. …………………………………………………………………………………... 27 Figure 2.8: The isosurface from the original dataset (a) shown with the isosurface from the smoothed dataset (b). ………………………………………………………………. 28 Figure 2.9: An enlarged view of the joint between two bones is shown in an edge display with its location in the shaded isosurface of the hand. ………………………………… 29 Figure 2.10: A slice plane oriented through the thumb displayed as a bitmap (a). A color map was applied to the slice plane (b). ………………………………………………… 30 Figure 2.11: Isosurfaces of four different values to investigate the separation of the bones at the joints. …………………………………………………………………………….. 31 Figure 2.12: Resulting separate entities, from separateBones, are shown in alternating shades. ………………………………………………………………………………….. 32

Chapter 3: Contours on Slice Planes to B-Spline Surfaces
Figure 3.1: Two representative slices through the index finger in the original orientation (a) and in the desired “natural” orientation (b). ………………………………………... 33 Figure 3.2: Uniformly spaced grid of points to resample the CT data. ………………. 34 Figure 3.3: The central axis shown with the three orthogonal views of the index finger used to determine its location. …………………………………………………………. 35 Figure 3.4: Translation of point P0 to the origin. ……………………………………… 36 Figure 3.5: Rotation of line P’ about the x-axis onto the x-y plane. …………………... 37 Figure 3.6: Rotation of line P’’ about the z-axis in line with the x-axis. ……………… 37 Figure 3.7: Points of the original dataset are illustrated as x’s and the resampled data points are illustrated as o’s. ……………………………………………………………. 38 Figure 3.8: Sample contour and B-Spline fit to the contour vertices. ………………… 41 Figure 3.9: Close-up of a self-intersecting loop and its location on the contour. ……... 42 Figure 3.10: Vertices of the resampled and original contours. ………………………... 43 Figure 3.11: B-Spline obtained from the resampled contour, with equally spaced points, shown with the original contour. ………………………………………………………. 44 Figure 3.12: Two separate non-uniform B-Spline curves, which share one end point. .. 46 Figure 3.13: Curvature plot for the two non-uniform B-Spline curves. The discontinuity is enlarged in the circle. ………………………………………………………………... 46

viii

Figure 3.14: The curves of Figure 3.12 can be defined as uniform B-Spline curves through trimming longer curves. ………………………………………………………. 47 Figure 3.15: C2 curvature continuity maintained through the joint of two uniform BSpline curves. The B-Spline curves were obtained through the trimming operations of Figure 3.14. …………………………………………………………………………….. 48 Figure 3.16: Non-uniform B-Spline curves (a) were trimmed to the uniform B-Spline curves (b). ……………………………………………………………………………… 49 Figure 3.17: Second direction cross curves fit through the end-points of the trimmed curves. ………………………………………………………………………………….. 50 Figure 3.18: The boundary blend surface (a) was trimmed to obtain the uniform portion (b). ……………………………………………………………………………………… 51 Figure 3.19: Two consecutive surfaces where the surfaces overlap prior to trimming the second surface (a) and the trimmed surfaces with C2 continuity (b). ………………….. 51 Figure 3.20: Resulting surfaces with C2 continuity. …………………………………... 52 Figure 3.21: Image of a CT slice to be imported into Pro/Engineer. …………………. 53 Figure 3.22: The image of Figure 3.21 was imported into Pro/Engineer and is shown offset from the imported curves used to surface the finger. …………………………… 53 Figure 3.23: B-Spline curves created through opposite groups of trim-points. ………. 54 Figure 3.24: End-cap curves intersecting at the centroid from the last curve projected onto the following slice. ……………………………………………………………….. 55 Figure 3.25: Boundary blend surface from the end-cap curves with their intersection point and the parallel curves from the slices. ………………………………………….. 56 Figure 3.26: The surfaces forming the end-cap of the bone shown with their boundary curves. …………………………………………………………………………………. 57

Chapter 4: Summary and Conclusions
Figure 4.1: Comparison of the CAD models created utilizing CT data and the surface of an actual bone. …………………………………………………………………………. 59

ix

Chapter 1: Introduction and Background Information
1.1 Introduction The focus of this research was developing a method to create Computer-Aided-Design (CAD) models for bones of the human body utilizing medical imaging. The human hand was chosen as the subject of the research due to past interest in the area. In this research medical imaging provided the data set from which curves and surfaces were developed. This allows the geometry of the surface models created in the CAD environment to accurately represent an individual’s form and structure. 1.2 Medical Imaging Technologies Medical imaging provides a data set from which curves and surfaces will be developed. This allows the geometry of the surface models in a CAD environment to accurately represent an individual’s form and structure. 1.2.1 Computed Tomography Computed tomography (CT) utilizes the transmission of X-rays through a material from many viewing angles. A single X-ray tube is rotated through a full rotation while recording projections at every 0.5° to 1°. This X-ray beam originates as a point source and after passing through the body is collimated to a single slice on an arc of detectors. These detectors measure the variation in transmission of the X-rays based on attenuation due to tissue opacity, which correlates to density. The projections of the X-rays are processed together to form an image of the densities on the slice. Each slice is a 2-D image with a defined thickness. The pixels, or cells, of the image can be thought of as voxels or a volumetric cube from which the density was measured. Sequential CT sections are imaged separately, by changing the position of the X-ray source and strip of detectors along a path perpendicular to the image, to construct a 3-D volume image1, 2. 1.2.2 Magnetic Resonance Imaging Magnetic Resonance Imaging (MRI) technology measures electromagnetic energy from resonating hydrogen protons. A magnetic field is generated to align the hydrogen protons in the subject’s body along the bore of the scanner. Gradient magnets are used to generate radio frequency (RF) pulses, which excite the unbalanced protons. These protons absorb that energy and spin or precess in a different direction. Once the RF excitation is stopped, these protons precess back to their alignment in the magnetic field. In doing so, they emit electromagnetic energy in the

1

form of an RF signal at the characteristic Larmour frequency. The proton densities and relaxation times are reconstructed in slices similar to the CT images1, 2. 1.2.3 Decision to Use CT Data The decision to use CT data over MRI was based upon the material of interest. MRI has advantages in distinguishing between types of soft tissue. Where CT imaging provides greater sensitivity between hard and soft tissues1, which aids in the identification of bone. Medical imaging was made available through the National Library of Medicine’s (NLM) Visible Human Project3. Both CT and MRI data are available upon approval from the NLM. 1.3 CAD Background The ultimate goal of this research is to accurately model bones of a human hand for use in CAD. It is important to understand how information is stored in a CAD environment. The industry standard for CAD software is parametric 3-D solid modeling software. In a CAD environment, parametric equations are favored over single-valued functions, such as y = f(x). A two dimensional parametric curve is expressed as two functions of a parameter u: x = x(u), y = y(u). This parametric form simplifies the modification of the curve. For example a line can be defined as y = mx + b, where m is the slope of the line and b is the y intercept. The parametric equation of a line, over the interval 0 ≤ u < 1, is P(u) = (1-u)P0 + u P1 where P0 and P1 are the start and end points4. The characteristics of the line can easily be changed in the parametric form by selecting new start or end points. Modifying the line of the explicit function would not be as easy. The advantages of parametric equations result in the ease of modifying designs. Curves and surfaces allow for the modeling of geometric entities with complex curvature in multiple coordinate directions. B-Spline curves and surfaces have become the industry standard for this complex modeling. They enable greater control over curvature and are often used to produce “stylized” and ergonomic consumer products. 1.4 B-Spline Curve Background B-Spline curves typically consist of multiple curve segments. Each segment is defined by a few control points, which influence the shape of the curve. B-Spline curves have several advantages over other curve types including Hermite and Bezier curves. In the Hermite and Bezier curves a change in shape of the curve affects the entire curve’s characteristics. This is known as a global propagation of local change4, 5. However, in the definition of a B-Spline curve there is an on/off condition that allows a set number of

2

control points to influence each curve segment. This leads to a curve in which a local change affects a limited portion of the curve. Additionally, the degree of the curve is independent of the number of control points. This means that numerous control points can define a B-Spline curve while maintaining the order of the curve. Typically, in modeling curves and surfaces it is desired to maintain curvature continuity or C2 continuity. This ensures the transition from one section to another remains smooth at the d2/du2 order. Specifically, the tangency of the curve is continuous through this transition. A B-Spline must be at least fourth order (K = 4) to achieve curvature continuity across curve segments. 1.4.1 B-Spline Curve Definition A B-Spline curve is defined as: p( u)

∑ i n

p i⋅ Ni , K( u )

=0

(1.1)

There are K control points that define each segment of a B-Spline curve. Where the degree of the curve is (K-1). Basis functions define the influence each control point has on the curve’s shape over the parameter u. These basis functions are defined recursively as shown in Equation 1.2.

Ni , 1( u ) = 1 Ni , 1( u ) = 0

if

t i ≤ u ≤ t i+ 1

otherwise

and
Ni , k( u ) =

(1.2)

(u − ti) Ni , k−1( u) t i+ k− 1 − ti

+

(ti+ k − u) Ni+1 , k−1(u) t i+ k − t i + 1

for integer values of k, such that 2 ≤ k ≤ K. The variables ti are the knot values that make up a knot vector, which control the effective intervals of u for the basis functions. For an open curve, the tj are calculated using K as shown in Equations 1.3.

3

tj

0

if

jn

(1.3)

for integer values of k, such that 0 ≤ j ≤ n+K.

1.4.2 B-Spline Example An example of calculating the basis functions for a fourth order curve with 7 control points (n = 6) is shown in Appendix A. Each of the seven basis functions seen below contain a Ni,1(u) term, which controls the effective range of the basis function.
N0 , 4( u) = ( 1 − u ) N3 , 1( u ) N1 , 4( u) = N2 , 4( u) = N3 , 4( u) = N4 , 4( u) = N5 , 4( u) = 1 4
3

(7u 3 − 18u2 + 12u) N3 , 1( u ) +
12

1 4

( 2 − u) N4 , 1( u) 1 6 ( 3 − u) N5 , 1( u ) 1 6 ( 4 − u) N6 , 1( u)
3 3

3

1 12 1 6 1 6 1 4

(−11u3 + 18u 2) N3 , 1(u ) + 1 (7u 3 − 36u2 + 54u − 18) N4 , 1(u ) +
3

u N3 , 1( u) +
3

1 6

(−3u3 + 12u2 − 12u + 4) N4 , 1( u) + 1 (3u 3 − 24u 2 + 60u − 44) N5 , 1( u ) +
6 1 12 1 4 12

(1.4)

( u − 1) N4 , 1( u) + ( u − 2) N5 , 1( u) +
3 3

(−7u 3 + 48u 2 − 102u + 70) N5 , 1( u) + 1 (11u 3 − 114u 2 + 384u − 416) N6 , 1( u )

(−7u 3 + 66u2 − 204u + 208) N6 , 1( u)

N6 , 4( u) = ( u − 3) N6 , 1( u )

The Ni,1(u), seen below, act as a switch turning on/off portions of the basis functions for a range of u values.
N3,1(u) = 1 =0 N4,1(u) = 1 =0 N5,1(u) = 1 =0 N6,1(u) = 1 =0 for for for for 0≤ u