Content

Problem motivation

2. B Spline Representation

3. Data reduction and curve smoothing

• Smoothing strategy
• Point interpolation
• Data reduction algorithm

4. Results

1. Problem Motivation

In many computer aided applications, numerical models of actual geometry are obtaining using slice-and-reconstruct strategy. This procedure involves several characteristic steps:

1. Slicing of the given geometry,

2. Identification of the boundary for each slice

3. Digitalization of the slice contours,

4. Surface mesh reconstruction,

5. 3-D solid mesh generation.

Particular procedure for obtaining the numerical model has to satisfy several criterions:

a. To be efficient, i.e. to be automated as much as possible

2. To give accurate representation of the actual geometry
3. Numerical model has to be 'useful' i.e. analysis performed on it has to be done in real time for the application purposes.

In this work, obtaining a numerical model that is to be used in computer aided medical surgery is discussed. Examples of computer aided surgeries are hip or knee joint replacement. Computer aid in these cases includes visualization of knee/hip joint of the patient, performing numerical analysis regarding replacement surgeries regarding biomedical and biomechanical issues of the bone/tissue complex. Such analysis can be done for pre, during and postoperative purposes and should be done in a real time, in this case the time of surgery. Numerical model, thus, needs to represent the actual geometry accurately.

Usually, necessary data are obtained using CT scans. Such scans usually carry a lot of amount of data in the form of three-dimensional cells-voxels. Accurate geometry representation at the same time means a lot of data, making model to slow for the later analysis purposes. After contour identification step, amount of data per contour could be measured by thousands. Each slice of CT scan of femoral bone, for example, carries between 400 and 1000 contour points and there are between 100 and 200 scans per bone. Thus, to represent only the surface of the bone, one needs in average 100 000 points. This is unacceptable for real time analysis. Data reduction must be imposed on the raw CT scan data.

This reduction technique has to fulfill the following:

1. Has to reduce the amount of data significantly to the degree that will allow for later analysis purpose to be fast enough;
2. Has to be fast;
3. Reduced data has to represent original geometry accurately.

Particular motivation problem for the illustration purposes is Total Hip Replacement surgery. In such intervention, damaged hip joint is replaced with artificial implants. (Figures 1 and 2).

Figure 1: Damaged hip joint

Figure 2: Hip joint after Total Hip Replacement surgery

In order to perform such complicated operation, surgeon needs accurate data concerning the bone and joint stage. In addition, biomechanical analysis can provide the data useful for determining the best size, shape and position of the implants.

Contours obtained from extraction from CT scans carry usually a huge number of contour points. However, inspite the large amount of data, actual contour representation has carries and approximation error.

This is due to inaccuracy to the extracting methods used to represent boundaries. As a result, misrepresentations like sharp corners in the bone surface might occur. This is why, in or data reduction strategy, we turn not to actual, but approximate representation (still accurate enough). Also, smoothing is second important characteristic in this strategy. Main idea is to replace the original set of points, with much shorter one. Contour is replaced from many small linear segments, to small number of spline segments. Because of locality property, natural choice for spline representation is B-splines. They are natural for use in fitting algorithms. Continuation requirements are up to second order derivative and thus our choice is cubic B splines. Data reduction strategy consists of the following steps: 1. Starting interpolation - creation of cubic B-spline that will interpolate given contour points 2. Knot removal technique - remove as many knots from the starting vector such the difference between original spline and new splines are less than set error value. Error can be measured in different norms, but preferable in such norms that will give local error measure rather than an average meaning of error.

2.1. B-spline curve representation

Our goal is to represent a given set of points wit the parametric B-spline curve. B-spline curve f is of degree k with n control points is defined on knot vector t=(t₀, t₁, …, t_n+k). Written in terms of basis B-spline functions we have:

Where d_i stands for i-th control point and B_i,k is basis spline function of degree k starting at knot value i.

In our particular problem, we will consider nonuniform knot vector t, which knot values correspond to the actual distance between the contour points.

There are certain situations where it is useful to increase the degree of freedom of the spline by adding l additional knots to the already existing ones. Suppose m=n+l and new knot vector is r=(r₀, … , r_m+k) such that earlier vector t is subset of the new vector r.

Starting spline f can be written as a linear combination of this extended basis with unique coefficients c_j:

The problem is given k, n, m, t and r, compute coefficients c_j. There are several ways to do that. We can choose m new knot values e₀, …, e_m-1 and solve the interpolation problem:

 If t_j < e_j <t_k+k+1, j=0,1,…,m-1, then this system of linear equations has a unique solution for c coefficients.

Relation between spline coefficients in original and extended form can be written in matrix form:

or

Matrix B is called insertion matrix of order k+1 from knot vector t to knot vector r.

More details on this can be found in [1].

3. Data Reduction Strategy

In this section, a strategy for the 2D contour smoothing is described. We start with set of points representing a 2-D contour line. Goal is to obtain much shorter set of points that will approximate this contour accurately enough.

In order to start, form of smoothing function has to be chosen. Our choice is B-splines because of their characteristics described above. To complete this choice additional parameters should be determined.

For our purposes, we are using cubic splines (degree=3). There is no need for higher order splines. Cubic splines assure second order derivatives continuity which is enough in this case.

In the literature, several spline fitting algorithms have been described. Simplest way is to choose uniform partition on unit interval. Another way is to use arc length approximation. For our case, we choose chord length approximation. Starting from 0 eac subsequent knot has added value of Euclidean distance bettween points. Main problem with this is tat incorrect parametrization migt occur when the distance between neighboring points is too large.

Method for data reduction that is implemented in this work was proposed by Lyche [2]. In first stage, starting set of points is interpolated by closed cubic B-spline. Then knot removal algorithm is applied sequentially until process converges, i.e. until no more knot values can be removed for the given tolerance approximation. Main idea is that from initial set of contour points, only few are needed to describe the underlying geometry. We also have to keep in mind that even the starting data set carries errors in contour representation.

In the rest of this section description of data reduction algorithm implemented in this work is described.

1. Interpolation

We are given the set of points: P₀, …, P_n-1. To make it clear that curve is closed, extension is made as P_n = P_0.

As described above, knot vector corresponding to this data set is obtained as:

For numerical convenience normalization of the above vector is done as:

In order to produce B-spline basis that spans all the splines on the knot vector, extended knot vector is created using technique described in [4] (for the case of cubic splines):

Spline curve can be written as:

To calculate unknown coefficients d, we set interpolation equations:

To obtain full system of equations, we need 3 more conditions on coefficients d. Because the curve is closed and restrictions on continuity, we have;

Because of the locality property of B-splines and because degree of curve is 3 system of equation (8) has tridiagonal form and can be solved in linear time using Thomas algorithm.

Knot vector defined in previous section consists of three different parts.

where

Middle part of the knot vector t₂ consists of so called internal knots. The goal of knot removal procedure is to eliminate as many as possible internal knots.

The algorithm we are using consists of 4 main components. We refer to them as Weight, Rank, Remove and Approximation. In the first step, each internal knot value is assigned a weight that is in some sense measure how important is the knot in specifying the spline shape. The second step is to rank these knots. After ranking them, we remove the least important of them. Finally, given the reduced knot vector, we compute new spline which minimizes the ‘distance’ from the original spline.

When a reduced knot vector is computed, we may be able to continue reduction until no more knot values can be eliminate. At the end, if the reduced vector is equal to starting, tat means every knot value is important. If reduced vector do not have internal knots, that means all of them are removed. In practice, result will lie somewhere between these two extremes.

Algorithm for knot removal is:

1. Initialization step:

Determine starting spline f and knot vector t.

Set r⁰=t, g⁰=t

2. Removal procedure:

for l = 0,1,2,….

1. Compute weights for internal knots
2. Rank knots by their weights
3. Binary search on s (number of knots to remove)
1. Initialize s
2. r = r^l –s knots with minimum wights
3. new spline g = G(t;r)f (new spline is best approximation)
4. error = ||f-g|| (in specified norm)
5. if error < epsilon increase s
6. else decrease s
4. if s=0 stop
5. r^l+1=r, g^l+1=g

This algorithm generates a sequence of shrinking knot vectors and a sequence of approximations to the original interpolating spline. In practise te majority of knots are removed in first few iterations.

We now present important details of the procedure.

Let w_j denote the weight corresponding to knot r^l_j and g_j,r denote the subspace of spline functions on the knot vector r’ = r^l- r^l_j =( r^l_j,…, r^l_j-1, r^l_j+1,…). The idea is to let w_j be an approximation to the distance between starting spline and the new one based on reduced knot vector.

The weight of knot j is computed as:

w_j = ||g-g_j,r||_infinity,r where g_j,r=min||g-g_r,j|| for all g_j,r from subspace of splines on reduced knot vector.

Above minimization problem can be read as: Remove knot value t_j. For the new vector, pick the spline that will be the best approximation to the starting one in the specified norm (here we used infinity or maximum norm). For such spline, compute the error in approximating the original spline. Tat error will be the knot weight.

We now give short description on how to calculate knot weights.

Extending to the space of starting knot vector:

where

and

So if

We want to solve te problem for unknown vector .

owever, poblem has infinitely many solutions. We can obtain a solution by setting:

For c₁ we need to solve problem: find min||d¹-Bc¹|| where and matrix B is defined as:

Matrix B has k+2 rows and k+1 columns. A specific method for solving above minimization problem is presented in [1].

After the weights are computed knots need to be ranked. For each knot t_j we compute ranking number which is integer as follows:

We sort all v_j values in nondecreasing order.

We need to remove s nodes from the knot vector t. We remove all the knots whose rank is strictly less then v_s. Because many knots will have the same rank it is possible that the number of knots removed id less than s. To choose remaining knots we implement such technique that will eliminate knots uniformly across the subscripts.

If we need to remove p such nodes, and u₀<u₁<…<u_q are the knots whose ranking number is equal to v_s listed in the order in which they occur in knot vector, we remove

where

This problem is similar to one in weight calculation step but more general since we are removing more tan one knot. We want to compute spline g given starting spline f degree spline k-1 and knot vectors t and r such that distance between splines g and f are minimal in specified norm.

Spline f is given as:

Spline g is described as:

such that: ||f-g||_L2norm is in minimum.

The coefficients q_i which define spline g on starting knot vector are computed with the discrete splines insertion matrix B which is computed by the Oslo algorithm (Lyche & Morken).

Data reduction strategy described in previous section is implemented on several contour data obtained from CT scans done in UPMC Shadyside hospital. Data are provided by the courtesy of the Center for Orthopaedic Research.

CT scans performed are for the femoral bone.

Figure 3 presents procedure for contour extracting. This is done using image processing package AVS. Contour extraction is based primarily on color specifications and rate of change for the particular portions of the image. However, manual interventions are needed especially in the regions where clear color difference does not exist.

On Figure 4, results for contour smoothing for several different contours are presented. results show that on average 40 points per contour are enough to accurately represent underlying geometry.

On Figure 5 influence of choosing the norm in which knot weights are to be computed has shown. We can see that norm which carries local information in it gives better results. L₂ norm gives good approximation in average sense but locally, there might be a significant inaccuracy in curve representation.

Figure 6a shows complete CT scan data in the order they are done in reality. On figure 6b approximate smooth contours are shown. Differences in the shape are negligible but there is enormous saving on stored data.

[1] T. Lyche and K.Morken "A Data Reduction Strategy for Splines" Research report no. 107, Institute of informatics, University of Oslo, 1987

[2] T.Lyche and K.Morken "Knot Removal for Parametric B-spline Curves and Surfaces" research report No. 109, Institute of informatics, University of Oslo, 1987

[3] T.Lyche and Ricard Riensenfeld "Discrete B-Splines and Subdivision Techniques in Computer-Aided Geometric Design and Computer Graphics", Computer Grapics and Image Processing, Vol. 14, 87-111 1980

[4]Fujio Yamaguci "Curves and Surfaces in Computer Aided Geometric Design", Springer-Verlag, Berlin, 1988

Figure 3: Contour extraction using AVS

Figure 4: Results of curve smoothing for several different contours

 

Figure 6: Contour smoothing for femoral bone