Download Estimating Curvature of Space Curves in Grey-Value Images: Rieger and van Vliet's Method and more Summaries Mathematics in PDF only on Docsity!
738 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 11, NO. 7, JULY 2002
Curvature of n-Dimensional Space Curves
in Grey-Value Images
Bernd Rieger and Lucas J. van Vliet
Abstract— Local curvature represents an important shape pa- rameter of space curves which are well described by differential geometry. We have developed an estimator for local curvature of space curves embedded in -dimensional ( -D) grey-value images. There is neither a segmentation of the curve needed nor a para- metric model assumed. Our estimator works on the orientation field of the space curve. This orientation field and a description of local structure is obtained by the gradient structure tensor. The ori- entation field has discontinuities; walking around a closed contour yields two such discontinuities in orientation. This field is mapped via the Knutsson mapping to a continuous representation; from a -D vector to a symmetric^2 -D tensor field. The curvature of a space curve, a coordinate invariant property, is computed in this tensor field representation. An extensive evaluation shows that our curvature estimation is unbiased even in the presence of noise, in- dependent of the scale of the object and furthermore the relative error stays small.
Index Terms— Curvature, gradient structure tensor, Knutsson mapping, space curves in -D.
I. INTRODUCTION
I
N THIS PAPER, we present a method suitable for curvature estimation of space curves, implicitly represented by grey- level isophotes (level-sets), in -dimensional ( -D) images. The curve is embedded in the image by a grey-level difference with respect to the background. Our method works directly on the grey-value information of the image; neither a segmentation is needed to detect the curve nor a parametric fit is done at any time during the analysis. The method exploits the differential structure of images. Curvature in two-dimensional (2-D) images has been well- studied, both in segmented and in grey-level images [1]–[5]. Curvature is the first-order shape descriptor of an object and, therefore, an important feature. In 2-D, it totally determines the shape of a curve. Isophote curvature [6]–[8] can successfully be applied to edges in 2-D and three-dimensional (3-D) grey-value images, but it fails when applied to lines (space curves) [9]. To over- come the problems associated with isophote curvature, one transforms the grey-value image into an orientation map from which the curvature can be derived after solving the disconti- nuity problem [1]. In 2-D, the use of the double angle method
Manuscript received July 24, 2001; revised March 21, 2002. This work was supported in part by the Netherlands Organization for Scientific Research under Grant 612-012-003. The associate editor coordinating the review of this manu- script and approving it for publication was Prof. Ioannis Pitas. The authors are with the Pattern Recognition Group, Department of Applied Physics, Delft University of Technology, 2628 Delft, The Netherlands (e-mail: bernd@ph.tn.tudelft.nl; lucas@ph.tn.tudelft.nl). Publisher Item Identifier 10.1109/TIP.2002.800885.
is well-known [1], [10], but in 3-D it remained an obstacle that prevented the computation of curvature in 3-D. Therefore, traditional 3-D methods are applied to segmented images, or even on curves represented by ordered points, which enables one to fit a parametric model to the curve [11]–[14]. These methods rely heavily on the preceding segmentation, labeling, and orderings steps, which may fail due to noise or the presence of a bundle of space curves comparable to a lock of hair.
A. Mathematics First we present the mathematics of differential geometry that describe space curves. We have adopted a formulation that can be applied to grey-value images (Section II). In 2-D, the curvature of a curve in every point describes the shape of this curve completely. In 3-D, a second parameter, the torsion, is needed to give a full description, in 4-D another, and so forth. For a -dimensional curve, we know from the central theorem of space curves, that for given curvatures there exists a curve with these and any two such curves differ only by a translation followed by a rotation [element of SO( )] [15], [16]. The curvatures therefore totally determine the shape of a space curve but do not tell anything about its position. This makes these parameters well suited as curve descriptors. The curvature is a first order feature of a -D curve. Let be an interval, then a -mapping is called a parameter curve and is called a space curve. Fur- ther, let be the parameter of then the tangent exists for all curves. The arc length is. If , i.e., is a regular curve then and are valid parameter transformations. In the following, will be the arc length and the derivative with respect to it. In this parameterization, we have the favorable properties and. A local orthonormal basis can be constructed iteratively for a curve if are linearly indepen- dent [15]. In this basis, the Frenet-equations can be formulated [16]. For a parameter curve , unambiguous numbers exist
for which holds , where. Here, the is the th curvature of the curve. From , we see
1057-7149/02$17.00 © 2002 IEEE
RIEGER AND VAN VLIET: CURVATURE OF -DIMENSIONAL SPACE CURVES IN GREY-VALUE IMAGES 739
immediately how to compute the curvature as and from the iterative construction of the basis we obtain
From (2), we see that the curvature is always greater or equal to zero. Indeed for space curves it does not make sense to speak of a signed curvature in a coordinate independent description. In contrary to closed surfaces there is no border separating two distinct parts of space. By choosing an origin one can speak of signed curvature also for space curves. Example: The circular helix, a parameterization is
where is the pitch, the radius, and the helix is winding around the -axis. With the parameterization of the curve and the arc length factor , can be computed using the previous formulas
II. THEORY OF^ SPACE^ CURVES IN^ GREY^ -V^ ALUE^ I^ MAGES
The local orthonormal basis can be obtained from the grey-level images itself by local orientation analysis.
A. Orientation Field: A Local Orthonormal Basis
In order to obtain the orientation field along a space curve embedded in a -D image we use the gradient structure tensor (GST) [3], [9], [17]–[20]. For any image , we can always com- pute the GST: , where the lower index denotes a partial spatial derivative. The overhead bar de- notes smoothing which is done per element of , where each element is a -D image. The GST can be expanded in terms of the eigenvalues and eigenvectors as , with. The eigenvectors of the GST con- tain information about the local structure in the image. We can compute the largest and smallest eigenvalue and the as- sociated eigenvector for any dimension of the image by using the power method [21]. For the smallest eigenvalue, has to be inverted, which becomes time consuming for large. In the 2-D or 3-D case analytic solutions are possible and much faster [9], [19]. For line like structures the tangent orientation is given by the “smallest” eigenvector. A normalized line detector is the ratio [9], [19]. All derivatives are implemented as convolutions with Gaussian derivatives. The scale denotes the resolution at which the Gaussian derivatives are computed. The size of the tensor smoothing defines how local the image structure is computed. The local set of eigenvectors of the GST consists of the same set of vectors as the local orthonormal basis. The ordering, however, is different.
B. Discontinuity of the Orientation Field Unfortunately the calculated orientation field contains a discontinuity mod [18], i.e., the direction of the line is un- defined. Computation of partial spatial derivatives of the orien- tation field are not possible without some preparation. In gen- eral, a mapping to a higher dimensional space is needed to solve the discontinuity problem. For example, in 2-D, the phase jump can be resolved by doubling the angle of the gradient vector [1], [10].
C. Knutsson Mapping Removing the phase jump in a -D field is not a trivial task. Knutsson has introduced a mapping
that removes that discontinuity [18] while being distance preserving. The mapping satisfies the following three require- ments: let
- uniqueness: , this removes the phase jump;
- uniformity: for , lo- cally preserves the angle metric;
- polar separability: , information car- ried by the magnitude of normally does not depend on the angle. For a comprehensive review on the problem of orientation rep- resentation, the connection between the Knutsson mapping, and the gradient structure tensor, see [22].
D. Curvature in -D Grey-Value Images Our goal is to compute the curvature via (2). We start with the tangent orientation obtained with the GST. Here, we already have the first derivative. The discon- tinuity problem (discussed previously) prevents direct computa- tion of the derivative of along. This is solved by mapping the tangent orientation via the quadratic mapping to a con- tinuous representation. The elements of form a new -D vector. The ordering of the elements in this vector does not influence our curvature estimation, because we only evaluate a norm, as shown in (8), which is independent of a permutation of the elements of. Now, we calculate the derivative of in the direction of the tangent [15], which is again a^ -D vector
where is the functional matrix
From the uniform stretch requirement of we know how to scale the norm of a variation vector. Thus, starting from (2)
RIEGER AND VAN VLIET: CURVATURE OF -DIMENSIONAL SPACE CURVES IN GREY-VALUE IMAGES 741
Fig. 2. Isosurface plot of a grey-level helix.
Fig. 3. Cross section through a test object.
B. Ring
To start off in 2-D, we generate smooth rings with different radii and signal strength. To measure the robustness of the estimator we add different levels of Gaussian noise. We use the definition SNR , where is the variance of the Gaussian noise. The results are depicted in Fig. 4. The different noise levels are always calculated at the same radii, but slightly shifted in the figure for a better display. The error bars indicate the standard deviation over 40 runs. A tensor smoothing and a gradient smoothing was used. For SNR dB
Fig. 4. Average curvature estimation on a 2-D ring for different noise levels.
Fig. 5. Average curvature estimation on a torus for different noise levels.
the average relative error is smaller than 10%. The estimator is unbiased.
C. Torus The torus image has a grey-value range from zero to one; therefore, SNR. Again, we use and
. The performance of the estimator is tested over 20 runs for different noise levels (see Fig. 5). The error bars indicate the standard deviation. The different noise levels are again shifted in the figure for better display. The error bars include even for the high noise level (3 dB) the true value.
D. Helix The curvature (4) of the helix depends on the two parameters and which scale the helix. For increasing size of the helix radius the curvature first rises, being at its max- imum at and then decreasing (Fig. 12). In order to make a scale invariant statement about the performance of our algo- rithm we sample the scale space , generate the according test
742 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 11, NO. 7, JULY 2002
Fig. 6. Curvature of the helix, scale invariant.
Fig. 7. Relative error, scale invariant.
images and compute the curvature of the center line. Therefore, we rewrite (4) to
and (13)
Now, we deal with dimensionless quantities and , which are suitable to show that our estimation works fine over a wide range of scales. In Fig. 6, we plot the theoretical prediction and our calculations, in which the different symbols indicate helices having either the same radius or pitch. Because our estimation stays so close to the true value over a wide range of scales, we conclude that for all values in between our sampled grid points the estimation works as well. The relative error of the estimation shown in Fig. 7 is smaller than 5% over a range from . 1) Influence of Noise:
- Noise for different helix scales In the same scale invariant manner as for Fig. 6, we investigate the performance under noisy circumstances. In
Fig. 8. Curvature estimation with added noise.
Fig. 9. Relative error in curvature estimation with added noise.
Fig. 8, the results are shown for 13, 19, and 25 dB, where the error bars indicate the standard deviation for 20 runs. The estimation is unbiased, since the error bars always intersect the true curvature. The relative error is shown in Fig. 9.
- Studies along a cross section of the helix profile For one helix ( ), we add noise (SNR dB) to the image and plot the computed curvature (mean over 20 runs) and the standard deviation in Fig. 10 along a cross section of the helix (see Fig. 3). Again we see that the mean remains around the true value, and the vari- ation stays approximately constant. It should be pointed out that we see here not only the isophote line at the point where the mathematical helix would lie but a cross sec- tion through the 15 pixel diameter of the helix. Due to the regularization effect of the GST, we can estimate the true curvature even if we are not at the exact position. If we choose a small tensor smoothing , then the relative error at the exact point becomes smaller but for the surrounding values it becomes larger.
744 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 11, NO. 7, JULY 2002
Fig. 14. Curvature along a circle contour.
V. APPLICATIONS
In fluorescent confocal microscopy, especially in biological life time applications, 3-D time series are acquired [24], [25]. Typically, these images contain moving bright spots. Here, it could be of interest to compute the acceleration of the spots. The acceleration is related to the curvature of the spatiotemporal space curve formed by the moving spots. The acceleration can be computed from the 4-D image with the notation as in Section I-A as follows [26]:
This is the classical expression of the tangential and the normal components of the acceleration. The velocity can be computed by, e.g., optical flow [24], [17] or via the gradient structure tensor. In a spatiotemporal 4-D image , the orientation of the space curve, formed by a moving spot, is a measure for the velocity of the spot [27], [28]. It should be noted that the eigenvalue analysis of the gradient structure tensor which is used to measure the 4-D orientation does not give direction information. We can, however, retrieve the velocity vector if we shift the discontinuity of the orientation field to the time dimension. This is reasonable as we know there is a causal connection between the time frames, i.e.,. We can retrieve the components of the velocity vector as fol- lows:
(16)
where are the eigenvectors of the gradient structure tensor.
VI. CONCLUSION
We have shown that the curvature of space curves embedded in -D grey-value images can be estimated using the formulas given by differential geometry adjusted to the higher dimen- sional space mapped by the Knutsson mapping. Our new estima- tion formula (8) reduces in 2-D to the known expression, which
clearly indicates that our work is consistent with older work. Furthermore, the estimation is unbiased, which even holds in the presence of noise. The curvature calculation is clearly indepen- dent of the scale of the objects as shown by our computations. Our grey-value based approach is clearly superior to a discrete curvature estimation.
ACKNOWLEDGMENT
The calculations were performed with the MATLAB toolbox DIPimage [29].
REFERENCES
[1] M. van Ginkel, J. van de Weijer, P. W. Verbeek, and L. J. van Vliet, “Curvature estimation from orientation fields,” in Proc. 11th Scandinavian Conf. Image Analysis , B. K. Ersboll and P. Johansen, Eds. Kangerlussuaq, Greenland, June 7–11, 1999, pp. 545–551. [2] K. Nordberg, H. Knutsson, and G. Granlund, “Signal representation using operators,” Linköping Univ., Linköping, Sweden, Tech. Rep. LiTH-ISY-I-1342, 1992. [3] B. Jähne, Digital Image Processing , 4th ed. New York: Springer, 1997. [4] M. Worring and A. W. M. Smeulders, “Digital curvature estimation,” CVGIP: Image Understand. , vol. 58, no. 3, pp. 336–382, 1993. [5] , “The accuracy and precision of curvature estimation methods,” in Proc. 11th IAPR Int. Conf. Pattern Recognition , The Hague, The Nether- lands, 1992, pp. 139–142. [6] P. W. Verbeek, “A class of sampling-error free measures in oversampled band-limited images,” Pattern Recognit. Lett. , vol. 3, pp. 287–292, 1985. [7] P. Breton and S. W. Zucker, “Shadows and shading flow fields,” in Proc. IEEE Conf. Computer Vision and Pattern Recognition , San Francisco, CA, 1996, pp. 782–789. [8] L. J. van Vliet and P. W. Verbeek, “Curvature and bending energy in digitized 2D and 3D images,” in Proc. 8th Scandinavian Conf. Image Analysis , Tromsø, Norway, 1993, pp. 1403–1410. [9] J. van de Weijer, L. J. van Vliet, P. W. Verbeek, and M. van Ginkel, “Cur- vature estimation in oriented patterns using curvilinear models applied to gradient vector fields,” IEEE Trans. Pattern Anal. Machine Intell. , vol. 23, pp. 1035–1042, Sept. 2001. [10] H. Knutsson, “Producing a continuous and distance preserving 5-D vector representation of 3-D orientation,” in IEEE Comput. Soc. Workshop on Computer Architecture for Pattern Analysis and Image Database Management , Miami Beach, FL, Nov. 18–20, 1985, pp. 175–182. [11] M. Boutin, “Numerically invariant signature curves,” Int. J. Comput. Vis. , vol. 40, no. 3, pp. 235–248, 2000. [12] N. Kehtarnavaz and R. J. P. de Figueiredo, “A 3-d contour segmentation scheme based on curvature and torsion,” IEEE Trans. Pattern Anal. Ma- chine Intell. , vol. 10, no. 5, pp. 707–713, 1988. [13] K. V. Mardia, R. J. Morris, A. N. Walder, and J. J. Koenderink, “Esti- mation of torsion,” J. Appl. Statist. , vol. 26, pp. 373–381, 1999. [14] D. Coeurjolly, S. Miguet, and L. Tougne, “Discrete curvature based on osculating circle estimation,” in Visual Form 2001, 4th Int. Workshop on the Visual Form , C. Arcelli, L. P. Cordella, and G. Sanniti fi Baja, Eds. New York: Springer, May 2001, vol. 2059, pp. 303–312. [15] I. N. Bronstein, K. A. Semendjajew, G. Musiol, and H. Mühlig, Taschenbuch der Mathematik , 4th ed. Frankfurt, Germany: Verlag Harri Deutsch, 1999. [16] M. Spivak, A Comprehensive Introduction to Differential Geom- etry. Berkeley, CA: Publish or Perish, 1979, vol. 2. [17] M. Kass and A. Witkin, “Analyzing oriented patterns,” Comput. Vis. Graph. Image Process. , vol. 37, pp. 362–385, 1987. [18] H. Knutsson, “Representing local structure using tensors,” in Proc. 6th Scandinavian Conf. Image Analysis , Oulu, Finland, June 19–22, 1989, pp. 244–251. [19] G. M. P. van Kempen, N. van den Brink, L. J. van Vliet, M. van Ginkel, and P. W. Verbeek, “The application of a local dimensionality estimator to the analysis of 3D microscopic network structures,” in SCIA’99, Proc. 11th Scandinavian Conf. Image Analysis , Kangerlussuaq, Greenland, June 7–11, 1999, pp. 447–455. [20] S. di Zeno, “A note on the gradient of a multi-image,” Comput. Vis., Graph., Image Process. , vol. 33, pp. 166–125, 1986. [21] G. H. Golub and C. F. van Loan, Matrix Computations , 3 ed. Baltimore, MD: John Hopkins Univ. Press, 1996.
RIEGER AND VAN VLIET: CURVATURE OF -DIMENSIONAL SPACE CURVES IN GREY-VALUE IMAGES 745
[22] A systematic approach to n D orientation representation , submitted for publication. [23] L. J. van Vliet and P. W. Verbeek, “Better geometric measurements based on photometric information,” in Proc. IEEE Instrumentation and Mea- surement Technol. Conf., IMTC94 , 1994, pp. 1357–1360. [24] C. B. J. Bergsma, G. J. Streekstra, A. W. M. Smeulders, and E. M. M. Manders, “Velocity estimation of spots in 3d confocal images sequences of living cells,” Cytometry , vol. 43, no. 4, pp. 261–272, 2001. [25] E. M. M. Manders, A. E. Visser, A. Koppen, W. C. de Leeuw, R. van Liere, G. J. Brakenhoff, and R. van Driel, “Chromatin dynamics during the formation of the interphase nucleus,” J. Cell Sci. , to be published. [26] T. Frankel, The Geometry of Physics. Cambridge, U.K.: Cambridge Univ. Press, 1997. [27] J. Karlholm, “Local signal models for image sequence analysis,” Ph.D. dissertation, Linköping Univ., Linköping, Sweden, 1998. [28] B. Jähne, Spatio–Temporal Image Processing. Berlin, Germany: Springer-Verlag, 1993, vol. 751 in Lecture Notes in Computer Science. [29] C. L. Luengo Hendriks, L. J. van Vliet, B. Rieger, and M. van Ginkel. (1999) DIPimage: A scientific image processing toolbox for MATLAB. Pattern Recognit. Group, Dept. Appl. Phys., Delft Univ. Technology, Delft, The Netherlands. [Online]. Available: http://www.ph.tn.tudelft.nl/DIPlib.
Bernd Rieger received the M.Sc. degree in physics from the Munich University of Technology, Munich, Germany, in 1999. He is currently pursuing the Ph.D. degree in image processing and analysis at Delft Uni- versity of Technology, Delft, The Netherlands. His research interests include motion estimation for 3-D image sequences and development of grey- value based estimators for analysis of structure in im- ages.
Lucas J. van Vliet was born in 1965. He received the M.Sc. degree in applied physics in 1988 and the Ph.D. degree (cum laude) in 1993. His dissertation entitled “Grey-scale measurements in multi-dimen- sional digitized images” presents novel methods for sampling-error free measurements of geometric ob- ject features. He is a Full Professor in multidimensional data analysis with the Faculty of Applied Sciences, Delft University of Technology, Delft, The Netherlands. He has worked on various sensor, restoration, and measurement problems in quantitative microscopy. His current research interests include segmentation and analysis of objects, textures, and structures in multidimensional digitized images from a variety of imaging modalities. Dr. van Vliet was awarded a fellowship from the Royal Netherlands Academy of Arts and Sciences (KNAW) in 1996.
