A multiscale analysis model applied to natural surfaces

MVA '92 IAPR Workshop on Machine Vision Applications Dec. 7-9.1992, Tokyo A Multiscale Analysis Model Applied to Natural Surfaces Frkdkric FALZON Gkrard GIRAUDON Marc BERTHOD INRIA, BP.93, F-06902 SOPHIA ANTIPOLIS Cedex Abstract Smoothing (regularization), interpolation and surface reconstruction are well known subjects in computer vision. The major difficulty is to choose a model well suited for one of these goals and driven by a minimum number of parameters. Another problem also arises when we want to do one of these opemtiow adaptatively, i.c. local features are processed in keeping m'th the application domain (e.g. cartography). Our goal is to present a discrete operator driven by only one parameter, allowing both global and local prvcessing of a surface and, well suited to smoothing, interpolation and surface reconstruction. 1 Introduction In this paper we present a discrete adaptative surface model well suited to the problem of surface regularization and rewnstruction. This model is based on a well known analog network analogy and, we show that, adaptativity is given by its ability to change its behaviour without modifying its intrinsic structure. A more complete study of this model has been previously given in 141 and [5] but here we focus more on applications. Until now, a lot of works have been published in reconstruction and interpolation of surfaces [6, 2, 1, 14, 15, 13, 121. Our contribution is focussed on the framework of DEM (Digital Elevation Model) processing and based on the description of a model and a discrete operator driven by a minimum number of parameters (only one). 2 A discrete scaled operator Models baaed on analog networks have been already used in various applications [ l l , 7,8,9].Our goal is to show that, from this kind of model, we can derive an operator to solve both smoothing from sparse data and surface reconstruction from relevant 2D structures. In the following we give the expression of this operator, then we give two examples. (The complete description of the model from which we obtain our operator is given in [4,5]). Let R be a set of nodes pi distributed on a line (1D m e ) , we note V(pi) = {pk E 0 , 1 i - k 1 = 1) , the neighbourhood of the node pi ; i n t 0 = {pi E R, pi) C R), the interior of the set a ; i3R = R - intR , the boundary of R ; In , the indicator function of the set R ; Let and 4 r. E lR+ the coefficient which connect two points i and j E R N x Nthe finite matrix given by with yi = 1 + 1 . kE "b., a: This operator will appear in a linear _system V = I?;' p,where f' is a given vector. In the 2D csse V is an image stored in a vector form and V is the resulting transformed image. 2.1 First case - Smoothing operator We wnsider an isotropic distribution of a , i.e. we make the assumption : at = a; = a = Const V i , k E int(R) . In this is contractive i.e. it decreases case the discrete operator ';?I the norm of all vector it is applied to (see [S]). In a geometric point of view, this operator reduce the curvature of all discrete function (vector). Then we cao solve the problem of smoothing from sparse data. We replace the previous 1D operator by its 2D analogous, and we consider dconnectivity. Let V be a 2D discrete function. The value of the grey level of the resulting image, called V is given by An example of smoothing is given in fig(3). The data test given in fig(3)(a) are voluntary the same as [13] to compare the results. The ha1 solution in fig(3)(g) is very closed to the solution obtain by [13] with thin plate interpolant. Second case : Reconstruction operator We consider an anisotropic distribution of a,i.e. a: # a; Vi # 2.2 k . In this case, it is possible to put asymmetrical weights between two points and therefore, the evolution law of the points are not identical. Consider the following scheme for natural surface rewnstruction [I] Take D the set of all points that belong to ridge lines or stream networks (assumed to be given by an external proms); [2] Two points i and j are. comected by a: in one direction in the reverse direction. if we put a: # one and by 4 4 can influence a point more or less, according to the.value affected to the corresponding parameter. If we want a point i not to be influenced by others, but it should itself influence the others, we take a; = m and 4 = a : then the point i remains the same for all values of a ,i.e. the location and altitude are p r e s e ~ e dacross the rewnstruction. [3] When all parameters are initialized, we solve the linear system r,V = for increasing values of a. We choose to tix only ridge lines, stream networks and boundary points, all the others are left to 0 (no information). Figs (4)(a,b) show the 3D representation of the original DEM, and the points that have to be fixed to make the reconstruction. These points are obtained with the algorithm described in [3]. Figs (4)(c,d,e,f) show the reconstruction obtained for different numbers of iterations (we solve I',Vntl = Vn with V" = V ). In figs (4)(g,h) we give the contour lines of the original and reconstructed DEM. These results are obtained by considering all surface points for the initialization of the algorithm. Some points are fixed (as in the first case) and the others are initialized with the original data. This leads to adopt a less brutal approach than in the first case where we have considered that no information is equivalent to 0. We can see that valleys are narrowed and the relief is more accentuated, but spatial morphology is relatively well conserved despite the smoothing. 2.3 A remark about fixed point across scales In [4] we show that if a is the same between all points r, is the discretization of the heat equation with certain boundary Suppose we are conditions (the 1D heat equation is fi = (f,.). in the ID case, then fixing a point with a> = oo and 4 = a comes to put a 1 on the diagonal element of the line i of r, and 0 elsewhere on this line, that is 3 Conclusions We have presented a discrete operator which moults from an electrical analogy (see [4, 51). By making two different ~ s u m p t i o n s on the neighbourhood relations, we proposed two applications based on the use of the same operator. The first wncerns smoothing from sparse data and the sewnd concerns natural surface rewnstruction. The results an encouraging and show the ability of this operator to take into sccount different neighbourhood relations. Moreover parallel implementation can be envisaged, due to the discrete nature and the structure of this operator. Acknowledgments We would like to thanks N.Maman, P.Leymarie, S.Mathieu, J.M.Mal6 and B.Vaaselle for their help in some parts of this study. References [I] A. Blake et A. Zis-an Cambridge MA 1987. Visual ReeoMtruction, MIT Press [2] T.E. Boult, J.R. Kenda Viaral arrface naoMttcrction wing sparse depth data, in Proc CVPR pp 68-76, 1986. This leads to an uncoupled system and this problem can be solved piecewise. However there may exist discontinuities of the first derivative on frontier points of each subdomain because the stationnary solution tends to be linear between two boundary points (see fig (1)). Suppose that V is the domain that includes each solution (noted fD. ) computed in each separate subdomain Vi of D . Each of these fD. is obtained by solving the problem on each subdomain delimited by fixed points (without considering the natural boundary points a n ) . Then V E i n t n , V = U V i is a set of contiguous domains and we can write the global soluCardD tion as fD = fDb lo, , where each fDb is solution of the k=1 following problem : assuming that the space variable z beiongs to [O,l] in each subdomain, { the diffusion coefficient and f the initial condition, we have to solve f? = { f; with conditions t ) = fDb(0) on a? and fDh(l,t ) = f (1) on a? . f Now if we want that on frontier points the solution belongs to C1 we have to introduce new wnditions at these points. Suppose we want a null first derivative on frontier points. Then the problem has to be solved with four boundary wnditions and this is impossible because the degree of the partial differential equation is two. Then we consider the following problem where we act on the curvature of the initial condition instead of the initial condition itself : with conditions fDb(O,t) = j D ~ ( 0 )on a? , fDL(1,t) = jDb(l) on a? , !f?!k9 = 0 on B? and a? . Here and a? denotes the &t and right boundary. This last equation has been evoked in [lO](pp.ll). Numerically, the connection between the scale parameter of this model and the first one (analog network) is straightforward by putting a = C-I . The boundary conditions are inserted in the new matrix r, by putting Neumann conditions (null derivative) in the f i s t and last line, and Dirichlet conditions (fixed points) in the sewnd and last but one line. In the initial case the numerical scheme was given by I,.' = I-dlA and here ,?J = I+a-I (here A is the laplacian). One can see here the utility to consider the problem on a finite domain. In dimension two this problem is different because the constraints are 2D primitives (e.g. ridge lines), not inevitably closed and, in this case, it is not always possible to make piecewise computations. 2 [3] J. Fairfield, P. Leymarie Drainage Netcoorrb Rum Grid Digital Elevation Models Water Reaourws Rmemch, Vo1.27, No.5, Pages 709-717, May 91. [4] F. Falzon, G. Giraudon, M. Berthod Vera un ModLle pour 19AnalyseMulti~chellede Surface 3D INRIA Report No 1639, 1992. [5] F. Falwn, G. Giraudon, M. Berthod A Multiscale Analysis Model Applied to Natural Surfaces IEEE Workshop on Applications of Computer Vision. Palm Spriwp, 1992. [6] W.E.L. Grimson A Computational Study of the Human Early Visual System, MIT Press 1981. [7] B.K.P. Horn Robot Vision The MIT Press McGraw-Hill Book wmpany - 1986. [8] A. Lumsdaine, J.L. Wyatt, 1.M.Elfadel Nonlinear Analog Networks for Image Smoothing and Segmentation MIT Artificial Intelligence Laboratory. A.I. Memo No. 1280 Jan 91. - [9] M. Mahowald, C. Mead Une d i n e en ailicium Pour la Science - No 165 - Jul91. [lo] T.Poggio, V.Torre, C.Koch An Analog Model of Computation For the Ill-Posed Problem of Early Vision MIT A1 Laboratory, A.LMemo 783, C.B.1.P Paper 002, May 84 - [ll] T.Poggio, V.Torre, C.Koch Computational Vision and Regularization Theory Nature Vo1.317. 26 Sept 1985 [12] S.S. Sinha, B.G. Schunck Surface Appruzimation Using Weighted Splines, in Proc CVPR Hawai Jun 91. [13] R. Szeliski Fast surface interpolation wing hierarchical basis function, IEEE PAMI, Vol 12(6), pp 513-528, Jun 90. [14] D. Tenopoulos The Computation of visible-surface npnsentotiom IEEE Trans on PAM1 Vol lO(4) pp 417-438, 1988. [15] D. Tenopoulos, M.Vasileacu Sampling and R w o ~ t r u c t i o n with Adaptative Mesh-, in Proc CVPR Hawai 91. Figure 1: Two connected subdomains. At the connection point (in the circle) the solution w not C' (4 + (b) Figure 2: Solutions of rgu (2) : (a) Some solutions with j(z) = b ( x ) 26(z - 1) ; (b) Stationnary solution with j(z)= z . (in the 2 cases the 2 extremal points of the x domaw are fized and the derivatives on these points are 0) Figure 3: (a) Original points ;(b)(c)(d)(e)(f)(g) Intermediary interpolation steps . - Figure 4: ( a ) Original DEM; ( b ) Ridge line, stream network and boundary points ( 814 points 20% of the original DEM points); ( c ) 2 iterations; ( d ) 20 itemtions; (e) 80 iterations; (f) 150 itemtions. (we have taken a = 2). (g) Contour lines of the DEM given in (a); (h) Contour lines of the reconstructed DEM.