M. Grabner, W. Trobin, M. Rüther, H. Bischof, R. Wozelka
Institute for Computer Graphics and Vision, Graz University of Technology, Austria
grabner@icg.tugraz.at
S. Millington,
MR Centre for Excellence, Medical University of Vienna, Austria; Center for Applied Biomechanics, University of Virginia, Charlottesville, USA
A. Jammernegg,
Department of Surgery, Elisabethinen Hospital Graz, Austria
We present a system to measure geometric properties (e.g., volume, surface area, and contact area) of the cartilage layer in human ankle joints. The cartilage and subchondral bone surfaces are sampled with a stereophotogrammetric device. A configurable pipeline of processing steps is applied to the sampled surface data to compute the desired quantities.
The meshes are cleaned (i.e., disconnected parts removed, holes filled, and noise reduced) and aligned to each other such that they represent the cartilage layer as closely as possibly. A hierarchical stitching approach creates from two separate (opposing) meshes a single closed surface, which is used for volume calculations.
The method has been evaluated with data from 12 lower leg specimens. While our approach confirms the results of previous contact area studies for low curvature regions of the inspected cartilage and bone surfaces, it is superior to existing (mostly mechanical) methods for highly curved regions.
We also discuss biomechanical testing of articular cartilage and present a custom-built prototype for obtaining stiffness maps of biological soft tissue.
Accurate quantitative descriptions of the surface geometry of joints are essential for assessing the state of the joint (e.g., after accidents) and selecting an appropriate therapy. However, it is currently unclear if state of the art in vivo methods such as 3D MRI techniques are sufficiently accurate for this purpose. It is therefore highly desirable to develop a “gold standard” by which existing measurement approaches can be validated, even if it is an in vitro method (i.e., one that can not be applied to the patient since the joint is destroyed by the measurement procedure).
In order to characterize the mechanical properties of a diarthrodial joint, accurate measurements of the articular cartilage thickness and the variation in thickness across the surface of the joint are required. It is also essential to have baseline measurements of quantitative geometric parameters of cartilage, such as thickness and volume in healthy joints, if we are to use cartilage thickness and cartilage volume etc. as image based biomarkers in degenerative joint diseases, such as osteoarthritis.
Most previous investigations of articular cartilage thickness have dealt with the knee, however, there have been relatively few studies of the ankle and other joints with thin articular cartilage layers. A schematic view of the human ankle joint is given in Figure 1 , taken from [Wikipedia, 2007].
A variety of methods have been used to measure the thickness of cartilage in different joints both in vitro and in vivo with varying accuracy. These include anatomical sections, needle probe measurements, stereophotographic techniques, X-ray measurements, computer tomography sections, ultrasound, and magnetic resonance imaging. Many of these techniques require disarticulation of the joint and/or may alter the thickness due to deformation during contact. Additionally several of these techniques use indirect methods to make measurements of cartilage thickness, e.g., X-ray and contrast enhanced CT.
MRI is becoming more widely available and offers many benefits over others methods for measuring quantitative parameters. Using modern segmentation algorithms and 3D reconstruction techniques MRI can be used to quantify changes in thickness and volume, in vivo.
Fig. 1: Human ankle joint (right foot)
12 fresh frozen foot and ankle specimens were taken from 12 male cadavers, mean age 61.5 years (range 51-75 years). All specimens were acquired in accordance with state and federal laws. Ethical approval for the study was provided by the University of Virginia institutional review board and human usage review panel.
The lower limb complexes were stored at -25∘C. Prior to testing each specimen was allowed to thaw at room temperature for 24 hours. After thawing the joint was disarticulated and the soft tissues were removed from around the tibia-fibula complex and the talus. Specimens were then potted in custom made cups using a fast setting resin (R1 Fastcast, Goldenwest manufacturing inc. CA, USA) ensuring the articular surface was above the level of the potting material. The talus was elevated above the potting material by inserting three screws into the inferior surface of the talus, leaving part of the screw shafts projecting into the potting material so that the screws became rigidly embedded into the resin and fixed the talus in position. During preparation the specimens were kept hydrated with phosphate buffered saline containing protease inhibitor (Sigma-aldrich, USA).
For the stiffness tests, cartilage specimens were acquired from patients undergoing total knee replacement surgery. Ethical approval for this study was granted by the medical University of Vienna ethics commission.
Fig. 2: Talus with intact cartilage layer rigidly mounted on the potting cup equipped with photo targets
The potting cups incorporated a flange with photo targets fixed to it, the rigid fixation of the specimen ensured that there was no motion of the specimen relative to the photo targets (see Figure 2 ).
The cartilage surface was then captured by the Advanced TOpometric System
(ATOS II SO, Capture 3-d, CA, USA), giving a point cloud consisting of
ca. 70,000 samples. The measurement error is specified by the vendor as 
2μm.
The articular cartilage was then dissolved by submerging the specimen in a 5% sodium hypochlorite solution for 6-8 hours to reveal the intact subchondral bone. During this process the specimen did not move relative to the photo targets.
The imaging process was repeated for the subchondral bone surfaces, and both surfaces were registered in a common coordinate system by evaluating the common photo targets.
We present a highly flexible system to derive geometric properties (such as volume, surface area, and contact area) from a set of meshes describing cartilage and subchondral bone surfaces. A configurable workflow is proposed which performs a series of operations on the input meshes to derive the desired quantities. We discuss potential problems due to limitations of the 3D acquisition procedure and propose solutions.
We review previous work in several areas related to our measurement workflow. Section 2.1 gives a brief overview of other cartilage measurement methods found in the medical literature (see [Eckstein and Glaser, 2004, Millington et al., 2007a, Millington et al., 2007b] for more details). The remainder of this section covers various computer graphics techniques useful in our context.
[Ateshian et al., 1994] developed an analytical stereophotographic technique for measuring the thickness of articular cartilage and applied the technique to the knee, but not to thin articular layers such as the ankle joint. The reported technique was a non contact method which allowed precise measurements (?90 μm). Unfortunately the process was slow, labour intensive and allowed only a discrete number of points to be digitized. Furthermore, it was not possible to reconstruct the entire joint surface from one data set.
Initial research on quantitative cartilage parameters has focused on the knee joint which has the thickest cartilage in the body and is relatively easy to segment as it does not have large congruent areas. However, attempts to utilize MRI for quantitative measurements of the cartilage layers of the ankle have been limited by the achievable image resolution and the techniques for determining the cartilage boundaries. With the exception of the work of [Eckstein and Glaser, 2004] there have been few attempts to analyse the accuracy of MRI measurements of cartilage especially in thin congruent cartilage layers. Moreover, the accuracy of MRI based quantitative measurements in thin congruent cartilage layers have not been validated against an independent gold standard.
There have been several reported experimental studies of ankle joint contact characteristics with wide variations in methodology and loading conditions resulting in varying reports of the extent and location of the contact area. The most widely used technique for measuring joint contact area in the ankle has been Fuji film [Calhoun et al., 1994]. This is a thin film containing ink bubbles which burst on pressure, allowing to estimate the contact area of two opposing surfaces. Fuji film has several well established limitations. It requires dissection of the joint so the film can be inserted, it is prone to crinkling, slippage and shearing which can produce erroneous results and it has a thickness on the order 0.3 mm which can alter the joint contact characteristics in highly congruent curved joints, such as the ankle.
From a biomechanical point of view, articular cartilage can be modeled as a biphasic material consisting of a solid and a fluid phase [Mak et al., 1987]. In the context of their biphasic theory Mak et al. also described a way to determine all material parameters (i.e. the permeability and the elastic coefficients of the solid phase) from a so-called indentation test. These tests are performed by indenting the material using a probe (the indenter), either using a constant load or a constant displacement and measuring the changes in displacement or load, respectively. The indenter material is chosen to have stiffness values orders of magnitude higher than the soft tissue under test. Commonly used probe tip geometries are cylinders, either with a flat or a hemispherical end, and pyramids, because they are relatively easy to model in boundary value problems.
The stiffness properties of articular cartilage vary significantly over the joint surface, depending on the predominant stress levels [Swann and Seedhom, 1993]. Therefore multiple indentation tests have to be performed to obtain a dense stiffness map. It is important to note that the test sites should not overlap, as it takes a few minutes for articular cartilage to fully recover from a previous indentation.
The state of the art in indentation creep testing is a custom built apparatus, e.g. as described by Roemhildt et al.[Roemhildt et al., 2006]. Their device is built using a servo-motor driven linear stage as actuator and the measurements are performed with a load cell and a laser interferometer to determine load and displacement, respectively. The cartilage thickness at the test site is estimated via a needle probe test [Athanasiou et al., 1991].
Among the large number of different surface distance metrics, the Hausdorff distance is particular useful and intuitive since it describes the largest distance found between any point on the first surface and the closest point on the second surface (and vice versa). This is important in the context of mesh simplification where the deviation of the simplified mesh from the original mesh must be assessed [Klein et al., 1996].
Often we are not only interested in the Hausdorff distance as a single measure of the distance between two meshes, but also in the minimum distance of individual points on one mesh from the other mesh. Several tools for mesh distance computation have been created in this context. The Metro system [Cignoni et al., 1998] calculates the Hausdorff distance for any pair of meshes given. Moreover, the surface can be colored according the the deviation from the original data for visualization purposes. Mesh is a similar tool with improved performance and reduced memory consumption [Aspert et al., 2002]
A common approach found in medical literature to estimate cartilage thickness is to compute the distance from a sampled point on the cartilage surface to the corresponding subchondral bone surface along the normal to the cartilage surface at that sampled point [Ateshian et al., 1991, Cohen et al., 1999]. This method has been reported to be more accurate than the vertical distance and proximity methods [Heuer et al., 2001]. Note, however, that the normal distance is not necessarily the shortest distance between a point on one surface and its nearest neighbor on the other surface, in particular in the presence of noise.
Sampled surface representations (i.e., point clouds) are often incomplete due to occlusion, accessibility restrictions, specular reflections, and similar effects. Filling the resulting holes by properly extending the existing surface is mandatory for creating closed surfaces.
[Davis et al., 2002] proposed a volume-based method where the surface is represented as the zero set of a signed distance function. By volumetric diffusion, this function is extended into empty regions until a closed non-overlapping zero set can be obtained.
Another class of methods is based on the moving least squares (MLS) method, where the surface is directly defined as the interpolation of the input vertices. Holes are filled by projecting their vicinity onto a 2D plane and triangulating the resulting polygon [Wang and Oliveira, 2003], or by directly operating in 3D space [Tekumalla and Cohen, 2004]. For more details on scattered data approximation see [Wendland et al., 2004].
Meshes derived from a physical measurement process are subject to noise. For many applications it is desirable to reduce the amount of noise, while at the same time keeping important features of the object. An anisotropic denoising approach inspired by image-space methods has been presented by [Fleishman et al., 2003]. Since smoothing methods tend to shrink the input mesh, [Desbrun et al., 1999] identified volume preservation as an important issue.
Stitching together independently created parts of an object is important in many situations, including modeling, surface reconstruction, and finite elements.
Techniques to reconstruct a surface from multiple range images include [Soucy and Laurendeau, 1992, Turk and Levoy, 1994]. These methods rely on a certain amount of redundancy in the data, i.e., overlapping surface patches. A polygonal surface with topological singularities can be “repaired” (i.e., cleaned from geometric inconsistencies such as small cracks in the surface and topological inconsistencies such as duplicate mesh primitives) with the approach by [Guéziec et al., 1998]. However, none of these methods can fill larger gaps in the mesh.
Many techniques exist to create surfaces from contours. [Meyers et al., 1992] presented a sophisticated approach for planar contours, which (with some modifications) is also applicable to contours that can be projected onto a plane. The method proposed by [Payne and Toga, 1994] is suitable for non-parallel (but still planar) contours. An entirely different approach, which represents the desired surface as an iso-surface of a partial differential equation, has been proposed by [Cong and Parvin, 2000]. However, while these techniques can deal with contours of complex topology, none of them is immediately applicable to contours of arbitrary geometry.
[Field, 2001] presented a method to create transitional volumetric meshes (composed by tetrahedral and pyramidal elements) between two given triangle meshes for applications in finite element analysis.
Fig. 3: Processing pipeline configured to compute area, thickness, and volume
We employ a user-configurable pipeline of processing kernels to define the series of operations to be performed on the input data. Each kernel takes at least one input mesh and creates a new mesh as output (such as the mesh stitching procedure, see Section 3.6 ) or performs calculations on its input (such as surface area or volume, see Section 3.7 ). Both input and output meshes are kept on a stack, which allows easy combination of the tasks according to the desired application.
The input data are created by the ATOS machine and its accompanying software, the data acquisition procedure is explained in more detail in Section 1.2 . The desired output is available after application of appropriate kernels, such as depicted in Figure 3 for the case of geometric queries about the data. The remainder of this section describes the individual building blocks in detail.
The output created by the ATOS system is a set of disjoint meshes of different sizes. However, in all our experiments the mesh we are interested in is simply the largest connected component. This is identified by a region labeling method, starting at an arbitrary triangle and recursively marking all connected triangles as belonging to the first region. This procedure is repeated for all unlabeled triangles, and all regions except the largest one are discarded.
We use a variation of the hole filling approach by [Wang and Oliveira, 2003]. It has already been noted by [Wang and Oliveira, 2003] that it is important to maintain a uniform triangle density across the hole to be filled and its neighborhood. The average triangle area A is computed over all triangles adjacent to edges along the border of the hole. The maximum triangle size is set to 2A, the minimum angle to 30∘. These constraints are taken into account by the triangulation algorithm [Shewchuk, 1996], generating a triangle mesh with similar density as its neighborhood.
Fig.4: Total surface area is most influenced by deviations in normal direction 
and least influenced in tangent direction 
(i.e., ∣An - A0∣ > ∣At - A0∣).
While additive noise on the vertex locations has little influence on the volume of a closed mesh (normal-distributed perturbations tend to cancel out each other), the impact on total surface area is larger. Intuitively, the maximum change in surface area occurs for a deviation in the direction of the local surface normal (see Figure 4 , this observation can even be used as a defining property of the normal vector [Meyer et al., 2002]). However, the area error is small if the vertex displacement is small compared to the triangle edge length (i.e., vertex distance), which turned out to be the case for the data we were using (see Section 5.2 for details).
In our system, we use the technique proposed by [Fleishman et al., 2003]. We reduce the mesh shrinking effect by applying the denoising procedure to the individual cartilage surfaces before stitching them together.
As a precondition for some of the following stages, we determine for every vertex on mesh M1 its nearest neighbor on mesh M2 (and vice versa). This search is accelerated by an octree. We use these vertex-vertex distances to find corresponding surfaces (Section 3.5 ) and for our recursive mesh stitching algorithm (Section 3.6 ).
Fig. 5: The point v on triangle mesh M closest to v0 lies on an edge or triangle adjacent to v1 (the vertex in M closest to v0) unless there exists an edge or triangle in M intersecting S1 with vertices v3 and v4 (and v5 in the case of a triangle, not shown here) outside of sphere S2. We did not observe this case in any of the cartilage data sets available to us.
(a) different mesh extents
(b) diverging meshes
Fig. 6: Mesh stitching procedure. Step 1: for each vertex vi on M, the nearest neighbor v′i on M′ is found. Step 2: non-corresponding parts of the mesh (due to (a) different extents or (b) artifacts caused by other tissue) are removed. Step 3: the gap between the meshes is filled.
However, approximating the Hausdorff distance between two triangle meshes by the vertex-vertex distances is sufficiently accurate only if the distances between the vertices within each of the two meshes is small compared to the distance between the meshes (see also [Heuer et al., 2001]). This condition is not met in our case since the thickness of the cartilage layer approaches zero at the border. To determine the distance of vertex v0 to the mesh M (see Figure 5 ), we therefore also calculate the distances between v0 and all triangles and edges adjacent to v1, where v1 is the closest vertex to v0 in M. The minimum distance found by this procedure (r1 in Figure 5 ) is our final result. Note that this method might fail to report the true minimal distance as indicated in Figure 5 by the vertices v3 and v4. Both of them are outside of sphere S2, but they are corners of a triangle intersecting S1. However, we did not observe any artifacts due to this possible inaccuracy.
Fig. 7: A situation where greedy mesh stitching (i.e., selecting the shortest edge) fails: after creating edge (v3,v3′), the method “gets stuck” at v3 since (v3,v4′) is shorter than (v4,v3′), then (v3,v5′) is shorter than (v4,v4′), and so on.
Although the cartilage and bone surfaces are registered to each other (see Section 1.2 ), their borders do not necessarily match since the surfaces are scanned independently. If one of the meshes exceeds the other one as depicted in Figure 6(a) , an incorrectly large distance would be reported for the exceeding part, which has therefore to be removed. We identify these offending regions by inspecting a triangle T = (v1,v2,v3) on M and the nearest neighbors vi′ of vi (i = 1…3) on M′ (see Figure 6(a) , v3 is not shown in this 2D sketch). If all vi′, i = 1…3, lie on the boundary of M′, triangle T is discarded (step 2 in Figure 6(a) ).
There is another source for possible errors. Other tissue (adjacent to the cartilage layer) that has not been completely removed may cause artifacts, making the cartilage layer appear thicker at the boundary (see Figure 6(b) ). To alleviate this problem, we define a maximum thickness d0 and shrink both meshes M and M′ until the distance between them is not larger than d0 anywhere at the boundary. Note that this procedure has no impact on the interior of the meshes (where their distance, i.e., the thickness of the cartilage layer, is typically much larger than d0), as long as two closed boundary curves exist (one on M and one on M′) with a distance from each other not larger than d0. Only if the above condition is not met (by selecting a too small d0) the method fails and discards the entire mesh. In our experiments, d0=2 mm gave best results (see Section 5 ).
The goal of this stage is to create a single closed surface (of genus zero) describing the cartilage surface, from which the cartilage volume can be calculated in a straightforward way. Since the distance of corresponding vertices along the surface borders is typically much smaller than the average mesh thickness, the choice of a particular tiling has little influence on the computed volume. We therefore employ a simple and efficient technique.
Due to the removal of non-corresponding parts of the mesh as explained in Section 3.5 , the boundaries of the two meshes to be stitched together are almost parallel. However, a straightforward greedy stitching algorithm (i.e., traversing both boundaries and selecting the shortest edge at each step) would fail in the unlikely but possible (and actually observed) case illustrated in Figure 7 .
stitch(0, n, 0, m) function stitch(i, j, i′, j′) if j - i > j′- i′ then if j - i = 1 then if j′ = i′ then create triangle (vi,vj,vi′′) else if d(vi,vj′′) < d(vj,vi′′) then create triangle (vi,vj′′,vj) create triangle (vi,vj′′,vi′′) else create triangle (vj,vi′′,vi) create triangle (vj,vi′′,vj′′) end if end if else k := ⌊  ⌋ k′ := argminl=i′…j′d(vk,vl′) stitch(i, k, i′, k′) stitch(k, j, k′, j′) end if else (similar as above with exchanged roles of M and M′) end if end function |
Fig. 8: Recursive mesh stitching: M has the larger boundary polygon, it is subdivided at vertex vk, the boundary of M′ is subdivided at the boundary vertex vk′′ closest to vk. The procedure is repeated recursively according to Algorithm 1 .
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Instead we recursively subdivide the pair of boundary polygons (v0,…,vn-1,vn = v0)
of M, (v0′,…,vm-1,vm′ = v0′) of M′, where we assume that v0 and v0′ are the closest
pair along the boundaries (more precisely, d(v0,v0′) ≤ d(vi,vj′), 0 ≤ i < n,
0 ≤ j < m, d(⋅,⋅) denotes the Euclidean distance). This assumption can always be
met by selecting indices accordingly. Then we subdivide the larger of the two
boundary polygons (in terms of number of vertices) in the middle (i.e., at vertex vi, i = ⌊
⌋, with the additional assumption n > m), and the smaller polygon at vertex vj′ closest to vi. We proceed recursively until the larger polygon segment contains
only two vertices. In this case we create either one triangle (if the smaller segment
has only one vertex) or two triangles (if the smaller segment has two vertices).
Algorithm 1
summarizes this procedure, Figure 8
shows the first subdivision step by
an example.
The total surface area of a mesh is simply the sum of the areas of all triangles in the mesh. The volume V is calculated as
where n is the number of triangles in a two-manifold mesh without boundary, and V (Ti) is the volume of a tetrahedron with base triangle i and apex v0. Although v0 can be chosen arbitrarily, we set it to the center of the mesh to make volume computation (almost) translation-invariant even under finite-precision arithmetic.
Since the average thickness d of the cartilage layer is small compared to the size of the entire joint, it can be approximated by the expression
where V is the volume of the closed mesh, and A1, A2 are the surface areas of the input meshes before they are stitched together.
A somewhat different task can be addressed in the same framework. Knowledge of the joint contact area under physiological loads and throughout the range of motion is essential for understanding the biomechanics of the ankle joint.
To determine the contact area between two adjacent cartilage layers under a given load, we mount the joint on a machine that applies to it a load similar to natural (in vivo) conditions and rigidly attach photo targets to the two opposing bones. The relative positions of the involved bones is measured from the photo targets while the load is applied. Then the load is removed, and the relaxed cartilage surfaces are measured independently, each registered in relation to the bone it belongs to (again using the photo targets). By transforming all measurements into a single coordinate system, we yield two overlapping meshes, where the overlapping area of the relaxed meshes is approximately the same as the contact area of the cartilage surfaces under load. This approximation is valid for relatively small deformations (such as in our case) according to the Winkler elastic foundation model [Johnson, 1985].
However, determining the overlapping regions is not straightforward since the
surfaces are not closed, so we need to apply some heuristics to decide which part of
the surface should be considered in contact (i.e., overlapping) with the opposite
surface. We define two constants ε1 and ε2 and treat all triangles on mesh A with a
distance of less than ε1 from the corresponding triangle on the opposite mesh B as
part of the contact area. Moreover, overlapping regions are considered in contact for
distances up to ε2 to account for a certain amount of deformation (see numbers “1”
in Figure 9
). We used ε1=0.01 mm and ε2=1 mm in our experiments. To distinguish
overlapping and non-overlapping regions, we examine the surface normal vector 
. In
the non-overlapping case, 
points towards the opposite surface, while the opposite
surface is found in the negative direction -
in the overlapping case (see
numbers “2” in Figure 9
). Finally we disregard regions that extend beyond the
opposite surface’s border similar to Section 3.5
(see number “3” in Figure 9
).
in non-overlapping regions and -
in overlapping regions; [3] regions of the
surface extending beyond the border of the opposing surface are excluded.
In addition to morphological measurements, recent developments in magnetic resonance imaging allow to acquire biochemical (i.e., functional) properties of articular cartilage in vivo. For a proper validation of these measurements a biomechanical model is required as a “gold standard.” In order to obtain the material parameters to build such a model indentation tests, as described in Section 2.2 , have to be performed. This section briefly sketches the system we designed and built for that purpose.
The central element of the setup (as seen in Figure 10(a) ) is a robotic arm with six degrees of freedom. The robot’s hand is extended with a custom-built vice to mount bone-cartilage specimens of roughly 10 × 10 × 7mm3. Opposite the vice, attached to a rigid frame, there is a 50N load cell for force measurements during the indentation tests, two grayscale cameras (only one of them is visible in Figure 10(a) ), and a pan-tilt unit (PTU) carrying a laser source. Further in the background a cold light source with a swan neck light guide has been mounted.
Our system uses a 70mm long stainless steel cylinder with a diameter of 3mm and a hemispherical tip as an impervious probe. Depending on the cartilage geometry, the indentation trials are performed at up to nine user specified locations on the surface. In order to accurately indent the cartilage specimen at a normal to the surface, the coordinates of the indenter tip and the orientation of the indenter have to be determined.
For every indentation the approaching vice is tracked with 72fps, and simultaneously the force data is acquired with 200Hz. As soon as the robotic arm reaches the target point 0.4mm below the cartilage surface, it stops and maintains this position for 600s while gathering force data of the relaxation process.
See [Trobin et al., 2007] for more details on mechanical setup, calibration, and surface reconstruction.
In this section we focus on the technical parameters of our methods. More details about numerical results (i.e., the actual cartilage geometry measured by our method) are given in [Millington et al., 2007a].
We evaluated our method on a number of meshes from the data set obtained according to Section 1.2 . For our evaluation, we performed two types of experiments:
Input size and execution time for different tasks are summarized in table 1 . The “input” time includes reading the mesh from a local hard disc and parsing its OpenInventor ASCII representation. The most time-consuming stage is the distance computation (Section 3.4 ). The various mesh cleaning tasks (mesh of interest, Section 3.1 , and hole filling, Section 3.2 ) are subsumed under “clean&fill”. Mesh stitching (Section 3.6 ) is not applicable for contact area calculation since there is no closed surface to consider. The “geometry” timing stands for the operations discussed in Section 3.7 , where no volume is calculated for the contact area experiments, hence the shorter times there. The performance of the method described in Section 3.8 is listed in the row “contact”. Finally, the resulting meshes are written to disc (“output”).
In this section we show that the vertex location noise in our input data is sufficiently low to not severely affect the results. We apply the denoising method [Fleishman et al., 2003] to several data sets with different parameters and compare the resulting quantities (area and volume).
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
We chose the closeness parameter σc to be equal to the average vertex distance in our data sets, which is 0.18 mm. Moreover, as in [Fleishman et al., 2003], we investigate a neighborhood of size r = 2σc around each vertex for the denoising algorithm. This gives an average of 15 vertices to be considered to determine the displacement per vertex and iteration. However, since we do not have any ground truth available and the surface does not contain any region that is known to be perfectly flat, we can not apply the semi-automatic parameter estimation procedure proposed by [Fleishman et al., 2003]. Instead, we choose the smoothness parameter σs=0.03 mm. Other values in the range between the vendor-specified location error and the average vertex distance gave similar results.
Table 2 shows results for different numbers of iterations of the denoising method applied to the mesh shown in Figure 11 . As discussed in Section 3.3 , the surface area error is larger than the volume error, but still negligibly low. We therefore did not include mesh denoising in the experiments presented in Section 5.1 .
|
Finally we investigate the influence of noise on contact area estimation. Since we determine the overlapping area of two surfaces (see section 3.8 ) which intersect at a small angle, one would expect a large error due to inaccurate vertex locations. However, under the reasonable assumptions that the errors are distributed around zero and that the denoising method [Fleishman et al., 2003] finds a good approximation of the error-free surface, the overall error in the contact area is negligible, too. Table 3 compares the estimated contact area for different numbers of denoising iterations. We report the relative contact area here to compensate for the area shrinking effect observed in Table 2 . The remaining contact area reduction with increasing iteration count is due to the volume shrinking effect, which we did not compensate for. These results demonstrate that the noise present in the data has no significant impact on the accuracy of the contact area evaluation.
|
For quick assessment of the numerical results, we applied standard visualization techniques. A color-coded thickness (i.e., mesh distance) map is depicted in Figure 11 , similar to the results achieved by [Cignoni et al., 1998, Aspert et al., 2002]. The contact area between talus and fibula for different loading positions is shown in Figure 12 .
Figure 13 shows the displacement (solid plot) and load data (dashed plot) as a function of time. Initially (t = 0s) the displacement is h = -6.0mm, i.e., the cartilage surface is 6.0mm away from the indenter tip. Then the robot accelerates and approaches the cartilage surface in normal direction. It should be noted that the apparent tremor in the displacement plot is mainly caused by the deviations from the simulated path due to the discrete stepping of the manipulator’s actuators. Additional errors result from the thereby induced vibrations and measurement uncertainties.
At t = 1.63s the cartilage surface touches the indenter tip (indicated by the two dotted lines) and the load starts to ramp up. Once the final position at a displacement of h = 0.4mm is reached the robot maintains this position and the stress relaxation begins.
We presented a framework to derive useful characterizations of human ankle joints from a set of surfaces generated by a stereophotogrammetric device. Due to its flexibility the system can easily be used to implement and evaluate various mesh processing algorithms. The performance is sufficient since realtime processing is not an issue. The results obtained by our method confirm results of previous studies for low curvature regions of the cartilage surface. Moreover, our method performs equally well for highly curved regions since it operates fully in 3D (unlike other methods which use 2D projections or require a thin film to be inserted into the joint which is subject to wrinkling). We demonstrate that measurement inaccuracies (i.e., noise) is not an issue by showing that explicit noise removal has negligible impact on the results.
The next important step in this project is to compare the geometric properties derived by the system proposed in this paper with measurements from MRI data. If the MRI-based methods can be shown to be sufficiently accurate to evaluate the state of the joint (and to decide on further therapy), the benefits for the patients will be huge since the inspection of the joint can then be entirely non-invasive. We expect the method to be applicable to other joints as well without significant modifications.
Future work on the mechanical measurement procedure includes a biomechanical evaluation of the gathered data via a biphasic model consisting of an isotropic homogeneous linear elastic solid and an incompressible, inviscid fluid.
This work was supported by the Austrian Science Fund under the grant FWF P18110-B15 and by Zimmer Inc. USA.
[Aspert et al., 2002] Aspert, N., Santa-Cruz, D., and Ebrahimi, T. (2002). MESH: Measuring errors between surfaces using the Hausdorff distance. In Proceedings of the IEEE International Conference on Multimedia and Expo, volume I, pages 705–708.
[Ateshian et al., 1994] Ateshian, G. A., Kwak, S. D., Soslowsky, L. J., and Mow, V. C. (1994). A stereophotogrammetric method for determining in situ contact areas in diarthrodial joints, and a comparison with other methods. Journal of Biomechanics, 27(1):111–124.
[Ateshian et al., 1991] Ateshian, G. A., Soslowsky, L. J., and Mow, V. C. (1991). Quantification of articular surface topology and cartilage thickness in knee joints using stereophotogrammetry. Journal of Biomechanics, 24(8):761–776.
[Athanasiou et al., 1991] Athanasiou, K. A., Rosenwasser, M. P., Buckwalter, J. A., Malinin, T. I., and Mow, V. C. (1991). Interspecies comparisons of in situ intrinsic mechanical properties of distal femoral cartilage. Journal of Orthopaedic Research, 9(3):330–340.
[Calhoun et al., 1994] Calhoun, J. H., Li, F., Ledbetter, B., and Viegas, S. (1994). A comprehensive study of pressure distribution in the ankle joint with inversion and eversion. Foot and Ankle, 15(3):125–33.
[Cignoni et al., 1998] Cignoni, P., Rocchini, C., and Scopigno, R. (1998). Metro: measuring error on simplified surfaces. Computer Graphics Forum, 17(2):167–174.
[Cohen et al., 1999] Cohen, Z. A., McCarthy, D. M., Kwak, S. D., Legrand, P., Fogarasi, F., Ciaccio, E. J., and Ateshian, G. A. (1999). Knee cartilage topography, thickness, and contact areas from mri: in-vitro calibration and in-vivo measurements. Journal of the OsteoArthritis Research Society International, 7:95–109.
[Cong and Parvin, 2000] Cong, G. and Parvin, B. (2000). Surface recovery from planar sectional contours. In ICPR, pages 106–109.
[Davis et al., 2002] Davis, J., Marschner, S. R., Garr, M., and Levoy, M. (2002). Filling holes in complex surfaces using volumetric diffusion. In Cortelazzo, G. M. and Guerra, C., editors, Proceedings of the 1st International Symposium on 3D Data Processing Visualization and Transmission (3DPVT-02), pages 428–438, Los Alamitos, CA. IEEE Computer Society.
[Desbrun et al., 1999] Desbrun, M., Meyer, M., Schröder, P., and Barr, A. H. (1999). Implicit fairing of irregular meshes using diffusion and curvature flow. In Rockwood, A., editor, SIGGRAPH 99 Conference Proceedings, volume 33 of Annual Conference Series, pages 317–324, New York. ACM Press. ISBN 0-201-48560-5.
[Eckstein and Glaser, 2004] Eckstein, F. and Glaser, C. (2004). Measuring cartilage morphology with quantitative magnetic resonance imaging. Seminars in Musculoskeletal Radiology, 8(4):329–353.
[Field, 2001] Field, D. A. (2001). Automatic generation of transitional meshes. International Journal For Numerical Methods In Engineering, 50:1861–1876.
[Fleishman et al., 2003] Fleishman, S., Drori, I., and Cohen-Or, D. (2003). Bilateral mesh denoising. In Hodgins, J., editor, SIGGRAPH 2003 Conference Proceedings, Annual Conference Series, pages 950–953. ACM SIGGRAPH, ACM Press.
[Guéziec et al., 1998] Guéziec, A., Taubin, G., Lazarus, F., and Horn, W. (1998). Converting sets of polygons to manifold surfaces by cutting and stitching. In Ebert, D., Hagen, H., and Rushmeier, H., editors, Proceedings of IEEE Visualization ’98, pages 383–390. IEEE, IEEE Computer Society Press.
[Heuer et al., 2001] Heuer, F., Sommers, M., Reid, J. B., and Bottlang, M. (2001). Estimation of cartilage thickness from joint surface scans: Comparative analysis of computational methods. In Proceedings Bioengineering Conference ASME, volume 50, pages 569–570.
[Johnson, 1985] Johnson, K. L. (1985). Contact Mechanics. Cambridge University Press.
[Klein et al., 1996] Klein, R., Liebich, G., and Straßer, W. (1996). Mesh reduction with error control. In Yagel, R. and Nielson, G. M., editors, Proceedings of IEEE Visualization ’96. IEEE, IEEE Computer Society Press.
[Mak et al., 1987] Mak, A. F., Lai, W. M., and Mow, V. C. (1987). Biphasic indentation of articular cartilage – I. theoretical analysis. Journal of Biomechanics, 20(7):703–714.
[Meyer et al., 2002] Meyer, M., Desbrun, M., Schröder, P., and Barr, A. H. (2002). Discrete differential-geometry operators for triangulated 2-manifolds. In VisMath.
[Meyers et al., 1992] Meyers, D., Skinner, S., and Sloan, K. (1992). Surfaces from contours. ACM Transaction on Graphics, 11(3):228–258.
[Millington et al., 2007a] Millington, S., Grabner, M., Wozelka, R., Anderson, D., Hurwitz, S., and Crandall, J. (2007a). Quantification of ankle articular cartilage topography and thickness using a high resolution stereophotography system. Journal of Osteoarthritis and Cartilage, 15(2):205–211.
[Millington et al., 2007b] Millington, S., Li, B., Tang, J., Trattnig, S., Hurwitz, S., Crandall, J., and Acton, S. (2007b). Quantitative and topographical evaluation of ankle articular cartilage using high resolution MRI. Journal of Orthopaedic Research, 25(2):143–151.
[Payne and Toga, 1994] Payne, B. A. and Toga, A. W. (1994). Surface reconstruction by multiaxial triangulation. IEEE Computer Graphics and Applications, 14(6):28–35.
[Roemhildt et al., 2006] Roemhildt, M. L., Coughlin, K. M., Peura, G. D., and Fleming, B. C. (2006). Material properties of articular cartilage in the rabbit tibial plateau. Journal of Biomechanics, 39(12):2331–2337.
[Shewchuk, 1996] Shewchuk, J. R. (1996). Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In Lin, M. C. and Manocha, D., editors, Applied Computational Geometry: Towards Geometric Engineering, volume 1148 of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag. From the First ACM Workshop on Applied Computational Geometry.
[Soucy and Laurendeau, 1992] Soucy, M. and Laurendeau, D. (1992). Multi-resolution surface modeling from multiple range views. In Conf. on Computer Vision and Pattern Recognition (CVPR ’92), pages 348–353.
[Swann and Seedhom, 1993] Swann, A. C. and Seedhom, B. B. (1993). The stiffness of normal articular cartilage and the predominant acting stress levels: implications for the aetiology of osteoarthrosis. British Journal of Rheumatology, 32(1):16–25.
[Tekumalla and Cohen, 2004] Tekumalla, L. S. and Cohen, E. (2004). A hole-filling algorithm for triangular meshes. Technical Report UUCS-04-019, School of Computing, University of Utah.
[Trobin et al., 2007] Trobin, W., Rüther, M., Millington, S., and Bischof, H. (2007). A vision-based system for biomechanical testing of articular cartilage. In Ponweiser, W., Vincze, M., and Beleznai, C., editors, Proceedings of the 31st AAPR/OAGM, Performance Evaluation for Computer Vision, pages 113–120. Österreichische Computer Gesellschaft.
[Turk and Levoy, 1994] Turk, G. and Levoy, M. (1994). Zippered polygon meshes from range images. In Glassner, A., editor, Computer Graphics (SIGGRAPH ’94 Proceedings), volume 28 of Annual Conference Series, pages 311–318. ACM SIGGRAPH, ACM Press. ISBN 0-89791-667-0.
[Wang and Oliveira, 2003] Wang, J. and Oliveira, M. M. (2003). A hole-filling strategy for reconstruction of smooth surfaces in range images. In Proceedings of SIBGRAPI conference, pages 11–18.
[Wendland et al., 2004] Wendland, H., Ablowitz, M. J., Davis, S. H., and Hinch, E. J. (2004). Scattered Data Approximation. Cambridge University Press.
[Wikipedia, 2007] Wikipedia (2007). Ankle joint. Visited on November 9th, 2007: http://en.wikipedia.org/wiki/Ankle.
