Discrete active models and applications

Pattern Recognition, Vol. 30, No. 5, pp. 817-835. 1997 ~ 1997 Pattern Recognition Society. Published by Elsevier Science Lid Printed in Great Britain. All rights reserved 0031 3203197 $17.00+.00

Pergamon

PII: S0031-3203(96)00120-3

DISCRETE ACTIVE MODELS AND APPLICATIONS SANDOR FEJES and AZRIEL ROSENFELD* Computer Vision Laboratory, Center for Automation Research, University of Maryland, College Park, MD 20742-3275, U.S.A. (Received 6 February 1996; in revised form 8 July 1996; received for publication 2 August 1996) Abstract--Optimization processes based on "active models" play central roles in many areas of computational vision as well as computational geometry. Unfortunately, current models usually require highly complex and sophisticated mathematical machinery and at the same time they suffer from a number of limitations which impose restrictions on their applicability. In this paper a simple class of discrete active models, called migration processes (MPs), is presented. The processes are based on iterated averaging over neighborhoods defined by constant geodesic distance. It is demonstrated that the MP model--a system of self-organizing active particles-has a number of advantages over previous models, both parametric active models ("snakes") and implicit (contour evolution) models. Due to the generality of the MP model, the process can be applied to derive natural solutions to a variety of optimization problems,including defining (minimal) surface patches given their boundary curves; finding shortest paths joining sets of points; and decomposing objects into "primitive" parts. (t) 1997 Pattern Recognition Society. Published by Elsevier Science Ltd. Discrete active models Optimization Minimal surface Shape decomposition

Geometric diffusion processes

l. INTRODUCTION Computational vision is usually defined as reconstruction of a scene from the images obtained by a visual observer. This task immediately leads to inverse problems, which due to the incompleteness of visual information, tend to be ill-posed, i.e. the problems do not have stable or unique solutions. Physics-based modeling as a regularization framework has become a favorite tool for attacking such problems by imposing additional constraints on the allowable solutions in terms of a (physical) model and recasting the original problem in terms of variational principles. (1) Therefore, optimization processes play central roles in many areas of computational vision, as well as many areas of computational geometry. Generally, active models are used to implement these processes. Since these problems typically require solving systems of partial differential equations, highly complex and sophisticated mathematical machinery is usually involved. Parametric active models, often referred to as "snakes, ''~2/ represent energy minimization splines which are controlled by constraints. This approach has been successfully used for various reconstruction tasks; however, in general, the method has limitations if the topology of the scene is not known, and this imposes restrictions on its applicability. Implicit models (contour evolution) (3) involve geometry-driven energy minimization principles. By using costly, higher-dimensional scene representations, (4) they are able to overcome some of the topological problems. Unfortunately, neither of these approaches is suitable for application to general sets, which is important in many situations, because both * Author to whom correspondence should be addressed. E-mail: [email protected]. 817

Minimal path

approaches make explicit or implicit assumptions about the scene. We may ask whether the solution of these optimization problems really requires such relatively complex mathematical tools, which at the same time impose limitations on the general applicability of these tools. It seems clear that if we want to find a simple and easy-to-handle approach to solving such problems we need to (1) incorporate energy minimization in an implicit, e.g. geometric, way, as is done in implicit active models. Such an approach also needs to (2) be applicable to arbitrary sets, and not only to specific types such as curves, surfaces, or solids. This should also (3) allow the problem of topology to be solved by simply ignoring it. It should be noted that many physical processes in nature (e.g. diffusion) conform to these principles. We also know that such processes are primarily based on local interactions of "active" particles, whose superposition gives rise to global properties of the system. In this paper we show that it is possible to formulate such an approach and to use it to solve even highly complex optimization problems. In Section 2 we briefly introduce the general concept of a Migration Process (MP). A more detailed study can be found in companion papers. (5"6) In Section 3 we relate the MP approach to active geometric models. In Section 4 we present applications of the MP model to various optimization problems. First, we demonstrate the generation of surface patches given their boundary curves. The problem of finding minimal paths joining pairs of points or sets of points is also addressed, and the construction of minimal surfaces is also demonstrated. Based on these elementary, optimization-based processes, we present an approach to three-dimensional shape decomposition, posed as a

818

S. FEJES and A. ROSENFELD

complex optimization problem. To tackle this problem we define a sequence of migration processes involving various constraints, and we present results on simple synthetic images. 2. MIGRATION PROCESSES

In this section we very briefly present a general definition of a migration process and define some general properties of constraints on such processes. Further details can be found in two companion papers. (~'6) It is understood that all sets considered here are subsets of the m-dimensional Euclidean space R m (m - 1 , 2 , . . . ) unless we specify otherwise. A migration process (MP) is an iterative mapping of a set P in which, at each iteration, each point p of P is displaced to a new position pP, whose choice depends on a subset of P (a "neighborhood" o f p in P). It is convenient to define the mapping of a set P -~ P' represented by an MP as the composition of two simpler mappings. The first mapping is called the neighborhood operator .~U : P --~ 2 P which defines the neighborhood of each point p as a subset of P. The second one is the placement function YP : 2 P -~ R m, which determines the location of the new point p' as a function of the neighborhood of the old point p. Note that images under ~ are not necessarily points of P. Given neighborhood operator ..,V and placement function ,~, we define the displacement function Y as the composition of ~,V"and ~ as follows: : P -+ P' = { p ' : p' = ~(..4 (p)); p 6 P}, where ,~- = ~.,-. Finally, a migration process is defined as iterative application of a displacement function: i f " : (P--~ p')". Continuity of a migration process is important in order to insure that connected sets do not break up, but remain connected. It is defined in the standard way, i.e. ,~- is continuous over P, if for all e>O, there exists 8>0 such that for allp, q C P, [p - ql < 6 implies ],~-(p) - ,~-(q)l < ~. The interaction of a migrating set with its environment is controlled by means of constraints. A hypersurface constraint is a restriction on an MP in which the displacement of points that lie on a hypersurface is restricted to that hypersurface. Formally, let H be a hypersurface in ~ " ; then ~ satisfies the hypersurface constraint defined b y / 4 , if P C_ H implies ,~(p) C_ H. A fixed point constraint aims at regulating the displacement of set P by fixing some of its points. Let F be a subset of P; then , ~ satisfies the fixed point constraint defined by F i f p E F implies ~-(p) = p.

2.1. Migration by averaging In our analysis we consider MPs in which the placement function ~ is defined by averaging. It can be shown that any set being migrated by averaging cannot increase in diameter, and remains bounded by the convex hull of the initial set. Using certain geometric restrictions on the type of the neighborhood we can obtain stronger results for

migrating sets. The migration of connected sets by averaging using equigeodesic distance neighborhoods is performed by taking each point p of set Q into the centroid of its 7>0 equigeodesic neighborhood {q: qCQAlp-ql ~-7}, where ] . ] ~ denotes the intrinsic (geodesic) distance induced by set Q. In companion papers, (5'6) we analyze the migration of closed curves and extend the results to more general connected sets. We show that the process is continuous for closed curves, i.e. the curve remains closed (a Peano curve), and furthermore, its total arc length is nonincreasing. The analysis of the convergence of this MP can be reduced to that of an MP defined on a polygonal approximation of the curve. In the case of the MP for polygons the migration of the polygon can be analytically described using its (moment of) inertia as a measure and performing a frequency domain analysis. We can show that migrating curves shrink. It can also be shown that the MP for closed curves can be considered as a discrete approximation to the evolution of the curve by a geometric heat equation, which has received extensive attention in the literature recently. (7's) This MP applied to each point of a curve produces a displacement which is approximately proportional to the curvature of the curve at that point. Due to the generality of the MP model the migration of more general types of connected sets (e.g. blobs, shells, patches) is straightforward. We have shown, (5'6) however, that the analysis of an MP defined by an equigeodesic distance neighborhood applied to more general types of sets becomes much more complex. These more general types of MPs are not necessarily continuous. We can qualitatively describe the migration of connected sets by exploiting the relation between the arc length distance function of a curve and the geodesic distance function of a connected set, Q, and reduce the analysis of convergence of this MP to that of the set of all closed curves embedded in the set. This yields a bound on the diameter of a migrating connected set, which also can be shown to shrink. As an analogue to the migration of closed curves, we have shown by geometric considerations that the displacement A p of a boundary point p of a connected set can be approximated as a linear combination of (mean) curvature-dependent and constant terms: Ap ~ ~(H0 + H), where ~ and /4o are positive constants, and H denotes curvature or mean curvature. These two terms can be considered as reaction and diffusion terms, respectively. The reaction component contributes to general thinning of migrating sets, whereas the diffusion term is responsible for smooth changes in the geometry. Useful properties of the combination of such processes were studied in a recent paper. (9)

2.2. Implementation of the MP model The MP model easily lends itself to simple, fixed grid implementations. We can implement a discretized ap-

Discrete active models and applications

819

Fig. I. DMPs of a space curve: free migration (top row), migration subject to a single fixed point constraint (middle row), migration subject to two fixed point constraints (bottom row).

proach in zm where no higher precision is used than the resolution provided by the grid. This can be done by combining the averaging with the quantization operator to yield a discrete version of the MP model, a Discrete Migration Process (DMP). ~5"6) As in the continuous case, one can restrict the neighborhood operator to discrete equigeodesic distance neighborhoods so that the geometric properties of MPs also hold (approximately) for these discrete cases. A DMP needs to be "continuous," in order to assure that connected sets remain connected during migration. Unfortunately, the continuity of DMPs based on discrete equigeodesic distance neighborhoods is not maintained; therefore interpolation needs to be incorporated. We have shown ~5'6~ that there exists an interpolation process that assures connectivity of migrating sets in ym and at the same time this DMP shrinks bounded connected sets in finitely many iterations.

2.3. Experiments with DMPs In this section we illustrate DMPs for different types of migrating sets in both two and three dimensions. We also illustrate the effects of fixed point and hypersurface constraints on the migration. In the implementation, we slightly modify the concept presented in Section 2.2. Using discrete time delays, ~5'1°)~ we limit the displacement of a point in one iteration to one unit distance, and conversely, we control the rate of displacement to a repetition rate defined by the time delay. Experimental results demonstrate advantageous properties of the implementation based on this model. 2 This implementation can be regarded as a discrete interpolation of the original DMP, producing a "dense" locus of the migrating points in the set.

JThis approach also has some features in common with Cellular Automata ModelsJ Ij~ 2The analytical study of the implementation based on time delays is beyond the scope of this paper.

2.3.1. DMPs for different types of sets. Figure 1 shows DMPs for a closed space curve both as free migrations and as migrations subject to fixed point constraints. Regarding the free migrations, we see that independently of the initial shapes of the curves, they become circular before they vanish. This is due to the fact that the migration of a closed curve (and hence also of its discrete version) is closely related to the evolution of the closed curve by the heat equation, as pointed out in Section 2.2. In this latter process the curve evolves proportionally to its curvature. It has been shown ~7"j2 ~4) that both planar and embedded closed curves evolve in such a way that no singularities are generated, and their total arc length decreases as fast as possible ("curve shortening flow"). In the case of a single fixed point, the migrated curve also shrinks to a point, but this point is the fixed point, independently of the initial shape of the curve. This is due to the fact that since DMPs are continuous processes, migrated sets never break apart. If two (or more) fixed points are involved in the migration, the migrating curve never vanishes, but reaches as stationary state a shape "stretched" between the fixed points, which is independent of the shape of the initial curve, as illustrated in the bottom row of Fig. 1. An alternative explanation for the existence of the stationary state can be obtained from a force model of the migration (see Section 3). According to this model, the stationary state is reached as soon as force equilibrium is established for the entire set, including the fixed points, as a closed system. ~ When we impose a hypersurface constraint, the displacement of each point of the migrating set interacts with the surface. In this case, only the tangential components of the displacement with respect to the surface are considered and the normal components are ignored. Figure 2 illustrates such DMPs for a spatial patch.

3In the force model fixed points are considered as external forces.

820

S. FEJES and A. ROSENFELD

Fig. 2. A DMP of a spatial patch subject to hypersufface constraints.

An important property of the MP model is that the migrating sets are treated as sets of individual points; differential geometric properties (e.g. curvature) of these sets are not involved in the iteration process as they are in the case of contour evolution modelsS '8'12-17~ This fact makes the model applicable to arbitrary sets with arbitrary topologies. We saw that the migration of the points in general ("thick") connected sets can be characterized by a boundary-curvature-dependent term and a constant term, in contrast to the case of a closed curve, where only a curvature-dependent term is present. We notice that, in general, the iterations at the beginning of the migration are a combination of constant and curvature-dependent migration terms [reaction and diffusion components~8'ls)]. If the patch or blob becomes elongated the migration gradually turns into a purely curvature-dependent one (diffusion). The curvature-dependent component has smoothing effects on the migrating set, while the constant term contributes mostly to the shrinking. The constant component plays an important role in the behavior of migrating sets and contributes a beneficial effect: shrinking can take place even in concave regions, which results in a general thinning of the migrating set. This property is clearly illustrated in

Fig. 3, where the concave border of a hole in the patch propagates toward the outer boundary due to the constant component of the migration (the curvature component works in the opposite direction). As soon as the two propagating "wavefronts" (one from the hole border, the other from the outer boundary) meet, the patch becomes a closed curve and subsequently migrates as a curve. Simple examples of the DMPs of blobs are shown in Fig. 4. Notice again that the migrating sets never break; this is insured by the self-adjusting mechanism of the combination of the two migration terms. It has been shown ~15~ that in three dimensions, closed surfaces (or solids) evolved by mean curvature can generate singularities (e.g. breaks) unless they are convex. For instance, the dumb-bell shape in Fig. 4 breaks along its elongated central region. In another paper, ~9) the "mean-Gaussian" curvature, a nonlinear combination of the mean and Gaussian curvatures, is used to eliminate this problem. The MP model by its nature is free from this problem.

2.3.2. Speeding up constraint propagation. Constraints play important roles in migration processes. They provide a means for the environment to interact with the process. It is important that the enforcement of constraints, particularly fixed point

Fig. 3. DMPs of a patch with a hole: free migration (top) and migration subject to a fixed point constraint (bottom).

Fig. 4. DMPs of a blob: free migration (top) and migration subject to a fixed point constraint (bottom).

Discrete active models and applications constraints, 4 take place in as few iterations as possible. This is especially important when the fixed points are far apart. Evidently, the propagation of fixed point constraints is primarily determined by the neighborhood distance k of the displacement function. Experiments with migration processes show (5) that the stationary state of a migration process involving fixed point constraints is practically independent of the choice of this parameter provided its value is small relative to the dimensions of the migrating set. If the neighborhood distance k is changed during a migration process, then as shown in our frequency domain analysis of the MP of polygons, (5'6) changes in the value of k even within a small range result in a modification of the frequency characteristics of the migration process, which has beneficial (though relatively moderate) effects on convergence. Changing the value of k has effects similar to those of multigrid methods, (2°'21) although we do not change the size of the underlying grid (i.e. we do not resample the data). For migration processes, speedup is important primarily when multiple fixed point constraints are used, since free MPs usually require relatively few iterations to reach stationary states (i.e. to vanish).

3. MIGRATION PROCESSES AS ACTIVE GEOMETRIC MODELS

In this section we introduce a physical model of the MP which provides a better understanding of the mechanism of the MP and its relationship with active geometric models. The model also allows us to treat constrained MPs in a regularization framework. Early computational vision leads to inverse mathematical problems which tend to be ill-posed: the existence, uniqueness, and stability of their solutions cannot be maintained without additional constraints. In order to make these problems well-posed, the visual data has to be regularized. °'22) The role of regularization is to insure the existence and uniqueness of the solution of this illposed problem by using a model which imposes constraints on the allowable solutions (such as continuity, smoothness, etc.), turning the problem into a well-posed variational (optimization) problem.

821

lyzed (5'6) its differential geometrical properties, and have concluded that the displacement of a point on the curve is (approximately) proportional to the curvature at that point. It can be easily shown that this type of motion corresponds to that of an elastic string submerged in a viscous medium. Since we also concluded that the migration of an arbitrary set is closely related to that of a set of closed curves, we can generalize this model and treat any migrating set as a deformable elastic shape. The elastic force model of the migration of a set Q can be formulated as follows: Each point q in the "),-neighborhood ~ ' ~ (p) ofp exerts an "attractive force" Fpqonp which results in a "virtual displacement" ~kpq ~- pq (see Fig. 5). Applying the superposition principle, the resultant force at p is thus

Fp = f Fpq Jq~s..(p) and the resulting displacement is thus the average of the individual displacements. This force model can also be used to account for the behavior of a migrating set, e.g. during transitions between shells, blobs, patches, and curves. If, as a result of the migration, the region around p becomes "thinner" than 7, the thinning automatically slows down, and finally stops, which is due to the change in the distribution of the elementary forces (Fig. 6). The migration of a set can be subjected to a hypersurface constraint H by considering only the components of the forces Fp tangential to the hypersurface; this causes the displaced point to remain on the surface. In the case of a fixed point constraint F, the fixed points p E F cannot be moved by any force. Since the MP has to be continuous, points in the neighborhood of p, which do not necessarily belong to F, are also affected by the fact that p is fixed; each of these points can be displaced only by a limited amount, to insure that continuity is not violated. We note that there are no interactions among the points of a migrating set, i.e. the displacements are independent. Hence, the physical behavior of a migrating set can be modeled as that of an "elastic fluid." The benefits of this "elastic liquidity" will be demonstrated in Sections 3.2 and 4.

3.1. A physical model of the MP The migration of a set is specified by a displacement function ~ : p ~ p' acting on the points of the set, which we can regard as the particles of an object. Physical processes in nature that act on the particles of an object are based on the superposition of elementary ("atomic") interactions among the particles (Huygens' principle). We can define a (qualitative) physical force model for our migration processes. Consider, for instance, the migration of a closed curve. We have ana-

4The "propagation" of hypersurface constraints is not necessary, since they are applied in parallel.

Q

Q

P

P

Fig. 5. The force model of the ME Each point p in the migrating set is subject to elementary forces Fpq by the points q in its 7-neighborhood (left); each of these forces results in a virtual elementary displacement Apq (right). The average of these virtual displacements is the actual displacement Ap.

822

S. FEJES and A. ROSENFELD

QI

QU

Fig. 6. The thinning of a migrating set (left) slows down (middle) and finally stops (right) as the thickness of set Q becomes less than % We can define an energy functional for a migrating set which fits constrained migration processes into a regularization framework. Let us define an elastic energy function Eelasti o which is purely based on the geometry of the set. Since the MP decreases the inertia of the migrating set, ~5'6) it is desirable to define as the inertia: Eelastic(O) d e f j ( Q ) ~ ~Q IP

goI 2 dp.

(1)

Since the migration can also be subject to constraints (hypersurface and fixed point constraints), these also need to be taken into consideration as an additional energy term Eco, caused by the external forces of these constraints. Since our constraints are not allowed to be violated (they are "hard constraints"), we have Econ(Q) = ( 0 oo

if the constraint is fulfilled, if the constraint is violated.

We define the total energy EMp(Q) of a migrating set Q as the sum of these two terms: EMp(Q) defEcon(Q) + Eelastic(Q).

(2)

Now we can relate the energy model of a constrained migrating set to the regularization framework of ill-posed problems as follows. Let H be a hypersurface and let F C H be a set of points lying on the hypersurface. Any migration process that satisfies the hypersurface constraint H and the fixed point constraint F can be regarded as solving the problem of defining a (connected) set K C H which lies on the hypersurface and contains F. It is clear that this problem is ill-posed (under-constrained), since there is an infinite number of such sets. We can restrict the set of allowable solutions by imposing the constraint that K has to have minimal inertia. However, even this restriction does not necessarily lead to a wellposed problem, since the solution may still not be unique. 5 Instead, let us construct an approximation to set K as follows: Let Q c H be any (connected) set lying on H and containing F, and define an MP in which Q is the initial set, subject to hypersurface constraint H and fixed point constraint F. 5For example, there can exist an infinite number of geodesics between two points on a surface.

We know ~5'6) that this MP decreases (and hence minimizes) the energy function, and it clearly leads to a unique result, given Q as initial condition. Therefore, we can consider Qn as a regularized (approximate) solution to the original problem, i.e. Q,, ~ K (with n and "~ as parameters), which has been obtained by minimizing the energy functional EMp (Q) [equation (2)]. We can relate this to the variational expression of the regularization framework by regarding the hypersurface and fixed point constraints as representing the "visual data" y available to the observer. Since by definition, Q, satisfies both constraints, i.e. it "fits" the visual data, the discrepancy term, which corresponds to E~,,n(Qn), is evidently zero. On the other hand, the regularization constraint corresponds to the elastic energy Eelastic(Qn), which is to be minimized given Q as initial condition. 6 We can conclude that the MP provides a regularization approach to solving certain ill-posed problems, where the regularization amounts to minimizing the energy functional EMr, [equation (2)] of the migrating set. In the literature, such problems are usually solved by using socalled active geometric models. In Section 3.2 we review previous work on such models and compare them to the energy model of the MP (in brief, the MP model).

3.2. Comparison of the MP model to other models There are two main approaches to active geometric models. Parametric models, often called "snakes, ''(2) have been extensively studied in the literature. They are specified by a set of parameters defined initially, which do not (necessarily) change during the optimization process. This class of active models can also be regarded as performing interpolation by spline functionals. Implicit models, ~3'4'23) on the other hand, are intrinsic and based on a partial differential equation; they do not depend on any particular parameterization. These models are also often referred to as contour evolution models. The two classes of active models represent very different approaches; however, from a

6Notice the potentially confusing terminology: the (hard) constraints of the MP play the role of the data discrepancy term, and the elastic energy of the set serves as the regularization constraint.

Discrete active models and applications mathematical standpoint, they share many common features. The MP model combines several advantageous features of both the parametric and implicit models, while at the same time it eliminates some of their difficulties. We also point out limitations of the MP model. The MP model can be regarded as an implicit, geometry-based energy minimization scheme, which thus categorizes it as belonging to the class of implicit models. The MP model is completely free of topological limitations, because it does not involve interactions among the points of the migrating set. This fact greatly simplifies its implementation, since there is no need for a higher-dimensional embedding surface as in the case of other implicit models. ~4) The basic approach used in the MP model is closely related to the active particles of Szeliski, Tonnesen, and Terzopoulos; ~24) however, since the MP model is implicit, it has a significantly simpler numerical implementation than this parametric approach. We saw that the MP not only handles topological changes but also has the ability to avoid certain types of singularities such as breaks, due to its force model (Figs 5 and 6); thus it inherently solves the problem of topological instabilities discussed by other authors. (15'25~ The geometric properties of a migrating set can be characterized by a combination of curvature-dependent and constant components (Section 2). Using the concept of the reaction-diffusion scale space introduced by Kimia, Tannenbaum, and Zucker, ~8) the curvature motion of a migrating set corresponds to a trajectory in this space. The diffusion (curvature-dependent) term produces smooth boundaries, while the reaction (constant) component produces constant propagation, which acts as the "inflation force ''<23) for parametric models, or as the "constant correction" term for the level propagation model, O) making the model capable of propagating into deep cavities. Since the MP model is suitable for direct implementation (i.e. there is no need for an embedding surface), it can be designed to incorporate user interaction, just like snakes; this is a significant advantage of parametric models over the contour evolution methods. A concept similar to the hypersurface constraint is used in some other papers, (26'27) where the authors generate geodesics on single-valued surfaces. The notion of fixed points appears in the "zip lock"-snakes of Neuenschwander and Fua, (28'29) where the authors use these points as basis points to generate minimal energy surfaces. In implicit geometric models the energy minimization process is based on a geometrical structure (PDE), and like parametric models, they are sensitive to initialization, the choice of constraints, and the difficulty of finding the global minimum, i.e. the convexity of the energy functional is not guaranteed. However, by means of a special initialization technique (see Section 4), exploiting the ("elastic fluid") properties of the MP model, small local valleys and ridges in the potential field of the energy function can be bridged; this approach resembles the Graduated Non-convexity algorithm. ~3°)

823

We also mention two limitations of the MP model as compared to the other active models: The current MP model cannot handle gray scale images. However, the MP concept could be extended to this type of data by modifying the averaging function. The MP model implements only the "hard penalty" concept in its constraints but not the notion of "soft" constraints; this restricts the applicability of the current model. 4. APPLICATIONS In this section we present applications of migration processes. First we demonstrate the generation of shortest paths joining pairs of points and sets of points for both Euclidean and geodesic distances. Then we give methods of constructing surface patches and minimal surface patches based on the MP model. Finally, using these results, we present a three-dimensional shape decomposition method and illustrate it with simple examples.

4.1. The loci of migrating (discrete) closed curves If we consider the locus of a migrating closed space curve, we would expect it to be surface-patch-like. Since the migration of a curve is in fact a discrete time process, this surface patch description is meaningful only if we use a parameterized model of the migration. Exploiting the relation between the migration by averaging over equigeodesic neighborhoods and the process of evolution by heat flow, let the surface patch : S 1 x [0, T) ---, R ~

(3)

be defined by

~ ( s , t) _ o2~(s, t) Ot

Os2

,

Ct(s,0) --- ~(s),

(4)

where the initial curve Cg(s) is the boundary condition, s is the arc length parameter of the curve, and t denotes time. In generating a surface patch using the locus of a migrating curve, we take advantage of the implementation of DMPs using time delays. This approach implements displacements of points which are limited to one unit of discrete distance. The resulting locus of migrating points is thus "dense," and does not require interpolation over several unit distances, as the original model would. Figures 7 and 8 illustrate the locus of a migrating discrete space curve. This simple method of generating a surface patch from a closed contour is especially important in three (and more) dimensions. Also, the method does not make any assumptions about the initial curve (see Fig. 8 for an example involving a knotted space curve, which yields a self-intersecting surface patch). In fact, the initial set does not even have to be a curve, but can be any connected set; because of the "dense" data provided by the time-delay-based implementation of the DMP, the problem of interpolation is also solved. Experimental results suggest that the three-dimensional surface patch defined by equation (4) resembles a minimal surface (see Figs 7 and 8), i.e a surface which

824

S. FEJES and A. ROSENFELD present an improved MP-based method for path estimation which tends to avoid local minima. Finally, in Section 4.2.3, an MP-based method of estimating the minimal length of a path connecting an arbitrary set of points is presented.

4.2.1. Using the M P model to f i n d minimal paths .

4.2. Finding minimal paths

Kiryati and Sz6kely (33~ propose a method of finding shortest paths and estimating distances between points on surfaces. This method combines length estimation of three-dimensional curves with graph-theory-based route finding. The costly process of exhaustive search in a sparse graph representing the voxels on the digital surface is accelerated by applying efficient heap data structures. A completely different approach is used by Kimmel, Amir, and Bruckstein. (34~ Equal geodesic contours (i.e. geodesic "waves") are propagated from the two points generating a geodesic distance map on the surface which then yields the geodesic route between the two points. The implementation is based on the level set propagation algorithm introduced by Osher and Sethian. ~4) We have already seen that migration of a connected set with (two) fixed points leads to a stationary state where the set is "stretched" between the points. In Fig. 9 a space arc is shown with its two endpoints defined as fixed points. As already demonstrated, the arc never vanishes and in the stationary state the remaining set, which acts like a stretched elastic string, forms an (approximately) straight line segment which can evidently be considered as an estimate of the shortest path between the two points. We can take advantage of this property and solve the problem of estimating the minimal path by minimizing the energy of the migrating set. The involvement of hypersurface constraints in the migration is also straightforward (Fig. 10). However, since the energy functional of the curve in this minimization scheme may be non-convex due to the hypersurface constraint, the migrating arc may be trapped in a local energy "valley." Deterministic optimization techniques based on local operations often suffer from this weakness. Stochastic techniques [e.g. simulated annealing (35~] or methods which involve exhaustive search in the parameter space [e.g. dynamic programming (36~] can overcome this problem, but usually at a high computational cost. Exploiting the properties of migrating sets,

The problem of finding minimal paths between points is of great importance in many areas of computational geometry (e.g. shape analysis, terrain analysis, motion planning, etc.). Especially important are methods which directly make use of discrete data, the most commonly available form of data, and do not involve continuous, differential-geometry-based techniques. In this section we show how migration processes can be applied to the estimation of shortest paths by exploiting the optimization features of the MP model. In Section 4.2.1 we give a brief review of two current approaches to the problem and outline the principles of our method using the MP model. In Section 4.2.2 we

Fig. 9. A DMP of a space arc with its two end points fixed. The stationary state of the migrated arc defines a path of minimal length between the fixed points.

Fig. 7. The locus of a migrating closed space curve.

\ :~.'<4a',

....

Fig. 8. The locus of a migrating closed, knotted space curve (note that the resulting surface patch is self-intersecting).

has the smallest possible surface area given its boundary curve. It can be shown that such a surface has mean curvature H = 0 [Plateau's problem~3J'32~]. This is the surface which is formed by, e.g. a soap bubble or an elastic membrane, since this corresponds to the leastenergy formation of the surface given its boundary. However, it can be shown that a closed curve evolved under Euclidean heat flow (or alternatively, migrated by averaging over equi-arc length neighborhoods) sweeps out a surface patch which is not necessarily minimal (see Appendix A).

Discrete active models and applications

825

Fig. 11. The DMP of a planar patch with two fixed points. The near stationary state (right) of the migration defines a nearly minimal path between the fixed points.

Fig. 10. The DMP of an arc subject to both fixed point and hypersurface constraints (note that the stationary state in this case is a local, but not the global, minimum).

specifically the thinning property in the migration of "thick" sets due to the constant term, we can often overcome the problem of local minima, as discussed in Section 4.2.2. 4.2.2. Avoiding local minima in the migration . Finding geodesics on surfaces is in general not a convex optimization problem. Geodesics are locally shortest paths in the sense that any small perturbation of the curve increases its length. ~37) In many cases there is an infinite number of geodesics between two points on a surface (e.g. antipodal points of a sphere). Taking these facts into account, we can conclude that proper initialization of the migrated set, restricting the possible solutions to a subspace of the solution space, should be used. An effective shortest path finding scheme can be defined by initializing the migrating set not as an arc as in Fig. 10, which imposes too strong a restriction on the possible solution space, but rather as a "thick" set of points. In Fig. 11 a planar version of this approach is illustrated. The patch-like initial set soon becomes thin and turns into an arc-like shape due to the characteristics of the MP model. Since the MP (a DMP) has no privileged directions, the process seeks a symmetric solution with respect to the shape of the initial set until the constant migration term disappears, which happens after the migrating set becomes thin. From this point on, only the curvature-dependent term remains; the energy minimization of the migrating set (see also Section 3.1) is not "biased" and the remaining arc seeks a minimal energy (shortest) path within the (local or global) energy valley in which it is currently located. 7 This initialization method usually allows the MP to overcome the traps of local energy valleys and find an

7We emphasize that the change in the relative sizes of the two components of the migration (constant and curvaturedependent) is not abrupt, but is gradual.

estimated geodesic path at the stationary state of the migration. Figure 12 illustrates this method of finding an estimate of the minimal path on a surface. Our energy-minimization approach has some features in common with the GNC (graduated non-convexity) algorithm introduced by Blake and Zisserman. ~3°) In this algorithm the non-convex energy function of the optimization process is approximated by a convex function which is able to bridge local potential valleys. Our MPbased method using our special initialization scheme can also be considered as a convex approximation to the original energy function, where the minimum of the approximating convex energy function, s which is superimposed on the original function at the beginning of the iteration process, is located near the medial axis of the initial set. 4.2.3. Finding minimal paths between sets o f points. The problem of finding a set of paths, connecting an arbitrary set of points, which has the shortest total length, is very important in many optimization problems such as complex shape analysis or network planning. When length is measured by Euclidean distances, this problem is known as the problem of finding the Steiner tree 9 for the given set of points. It has been shown that this problem is NP-hard, and it is impractical to solve it for more than about 20-25 points. ~38) It is clear that using an MP-based method with the given points as fixed points, we can easily construct estimated solutions to this problem, with no more difficulty than in the case of a pair of points presented in Section 4.2.2. In our approach, we again exploit the constant migration term of the MP model, which is available in the cases of patch-like and blob-like sets. This term allows us to construct thin connected sets even in concave regions; this is generally necessary if three or more fixed points are involved in the migration. Figure 13 illustrates DMPs of a planar patch subject to sets of fixed point constraints. In its stationary state, the initial patch-like set becomes thin (typically 1-2 points thick). It corresponds to the minimum (elastic) energy of the initial set with the constraints, and thus can be regarded as an estimate of the system of paths of shortest total length connecting all the fixed points. Figure 14

8This is represented by the constant migration tenn. 9Note the basic difference between the Minimum Spanning Tree and the Steiner Tree; in the latter case the edges do not have to end in points of the set, but vertices can be added to the set. It follows that the total path length of the Steiner Tree is at most equal to that of the Minimal Spanning Tree.

826

S. FEJES and A. ROSENFELD

Fig. 12. A DMP of a spatial patch (top right) subject to two fixed point constraints and a hypersurface consrtaint (top left). The stationary state (middle right and bottom left) of the migrated patch defines a minimal path on the surface between the pair of fixed points.

4.3. Defining discrete minimal surfaces ~? .......

9: d' Fig. 13. DMPs of planar patches (left) subject to different numbers of fixed point constraints. The stationary states (right) of the migrated patches define paths of minimum total length. shows a similar example in three dimensions for three connected sets of fixed points. The most general and hence most important shortest path problem involves a geodesic distance function, i.e. distance measured on a surface instead of Euclidean distance. We can easily extend our approach to this case by using the surface as a hypersurface constraint for the migration. We note that this generalization of the problem is more complex due to the geodesic distance function involved. Figure 15 illustrates an estimated system of paths of minimum total length for five points on a surface. For comparison, an estimated system of shortest Euclidean paths, i.e. a Steiner tree, is given for a similar distribution of five points in a plane.

Defining minimal surfaces, i.e. surfaces of minimal surface area, is also of great importance in optimization problems. Such surfaces have the shapes of stretched elastic membranes or soap bubbles on wire frames332~ Our MP-based method which defines minimal digital paths, as introduced in Section 4.2.3, can be straightforwardly extended to the construction of digital minimal surfaces. In order to do this, the set of fixed points has to be restricted to a closed curve of fixed points, and the initial set must contain all these fixed points and have no holes. As the stationary state of the migration of this initial set is approached, the set takes on its minimum (elastic) energy, and hence can be considered as an estimate of the minimal surface having the closed curve of fixed points as its boundary curve. Figure 16 shows a DMP of a patch (shell) with its rim as a closed curve of fixed points (top row). As the stationary state of the migration of this patch is approached, it takes the shape of a planar patch which evidently corresponds to the patch of minimal surface area having the rim as its boundary curve. On the bottom row of the same figure an additional fixed point is introduced, which is not connected to the original set of fixed points. These examples clearly demonstrate that the migration of a connected set can be best modeled as a strained elastic, viscous fluid rather than a strained elastic

Discrete active models and applications

827

Fig. 14. A DMP of a blob subject to fixed point constraints. The stationary state defines paths of minimum total length.

Fig. 15. A DMP of a spatial patch subject to a number of fixed point constraints (top right) and a hypersurface constraint (top left). The stationary state (middle right and bottom left) of the migrated patch defines a system of paths of minimal overall length on the surface connecting the fixed points. A comparable planar Steiner tree is shown on the lower right.

solid, since the particles in the migrating set have no point-by-point connections to their neighbors; they can slide along each other, just as particles in a fluid do. A different initialization scheme is shown in Fig. 17, where the generation of a minimal surface defined by a space curve is illustrated. In the absence of a hypersurface constraint, the D M P of a connected set regarded as an optimization process usually has a convex energy function; hence the stationary state of the migrating set is primarily determined by the set of fixed points (see also Fig. 1). The symmetric shape of the initial set in

Fig. 17 was chosen purely for convenience in generating the synthetic data. In our analysis of the surface patch generated by the locus of a migrating curve in Section 4.1, we showed that the generated surface is not necessarily minimal. We can now give an intuitive, qualitative explanation of this result, emphasizing the principal difference between the surface in Fig. 17, which is minimal, and the surface in Fig. 7, which is not. The basic property of a minimal surface is the global equilibrium of all the local forces. In the terminology of

828

S. FEJES and A. ROSENFELD

Fig. 16. A DMP of a spatial patch subject to a fixed point constraint along its boundary (top row) and an additional single fixed point at the bottom (bottom row). Notice the behavior of the patch in regions above the single fixed point; the particles can merge and slide over one another, so that the patch behaves like an elastic, viscous fluid rather than an elastic solid.

OOew Fig. 17. A DMP of a blob subject to a fixed point constraint defined by a closed space curve (top). The stationary state of the migration generates a minimal surface patch having the closed curve as a boundary. The bottom row shows the stationary state from several viewpoints.

the MP model, this can be considered as global enforcement of all the fixed point constraints. During the generation of a surface patch by the locus of a migrating closed curve, no such global enforcement is taking place. Instead, only the "constraints" m in the local neighborhoods are taken into consideration as each point along the curve is displaced to its new location. There is no significant time to propagate these local constraints to each point of the curve in one iteration. In contrast to this, global propagation of all the constraints takes place at one time in the DMPs illustrated in Figs 16 and 17. This difference is also noticeable in the flattened shape of the surface patch in Fig. 7, which was influenced only by the local neighborhoods, while the corresponding region of the minimal surface patch in Fig. 17 is less flattened due to the global enforcement of the constraints.

4.4. Shape decomposition Recognition of complex shapes is one of the most challenging problems in computational vision. Its com-

]OWe can consider the effects of attraction forces as constraints in a general sense.

plexity can be decreased by decomposing the complex shape into parts and recognizing the individual parts.(39,4°) In this section we first give a brief review of some of the major approaches to shape decomposition, emphasizing their advantages and disadvantages. Then we formulate shape decomposition based on optimization principles and present an approach using the MP model. Finally, we demonstrate our approach by experimental results using simple, synthetic images.

4.4.1. Approaches to shape decomposition. There are many different approaches in the literature, based on widely varying principles, to the problem of shape decomposition. None of these approaches has proved to be generally applicable to the problem; they are limited in many respects and require many assumptions. One of the b e s t - k n o w n m e t h o d s is b a s e d on d e c o m posing complex objects into sets of parametrically defined primitive shapes [superquadrics, generalized cylinders~4°-45)]. This technique, which is also known as the "recognition by parts" approach, has been used in many applications due to its simplicity. An advantage of this method is that it directly gives descriptions of the

Discrete active models and applications

829

parts from their parametric models. Since the flexibility by Hoffman and Richards but they left these questions of this method is limited in describing more complex, completely open. We shall see that such contours can still e.g. articulated, shapes, it is usually used in combination be used to suggest a part decomposition. It is convenient to describe surface points based on with other techniques. ~46) The construction of morphogenetic sequences from a their curvaturest66)--specifically, to represent them in a shape attempts to describe and decompose the shape by coordinate system defined by their two principal curvameans of a differential description of these sequences as tures. [Koenderink~67)uses a similar representation of the filter outputs. One of these approaches implements cur- surface points of a shape; he introduces polar coordivature-based contour evolution,~9'~8'~9"47~9)while others nates, in which the angle represents a shape index (conuse morphological filtering.~5°) The principles of this vex or concave umbilic, saddle, etc.) and the radius approach are relatively simple; however, the interpreta- represents curvedness.] In this representation, part-like regions and part boundaries can be associated with tion and evaluation of the results is quite limited. Skeleton-based decomposition techniques have different areas, as shown in Fig. 19. The regions on a surface consisting of points which proved their usefulness in two-dimensional shape decomposition and description problems. ~5~-53)There have satisfy rule (i) will be called part boundary regions. On also been attempts to extend these techniques to three the other hand, the regions whose points satisfy rules (ii) and (iii) will be called part-like regions. In order to dimensions.~54~ An important approach to shape decomposition de- decompose a shape into parts, we have to locate the fines parts of a shape based on the curvatures of their points of the part boundary regions and construct closed boundaries or surfaces. This represents a very general cut lines which define locations for the separation of the and comprehensive approach to the decomposition pro- parts. Unfortunately, we have already seen that part boundblem, but it has remained mostly theoretical so far. Psychophysical studies have established the key role ary regions do not always define closed contours of the boundary curvatures of shapes in the recognition [Fig. 18(c)]; they may also define separated curve pieces. of shapes by humans. Koenderink and van Doom ~55'56~ It is therefore necessary to find curve pieces which first suggested the importance of surface curvature, smoothly continue one another and to join them into continuous contours. This problem of grouping curve specifically of parabolic regions, in defining the "parts" of a complex shape. They showed that parabolic regions pieces into curves has been extensively studied. ~68)Evion a surface do not intersect and always define closed dently, it does not have a unique solution, and additional curves on the surface. A shape decomposition model constraints have to be used. The most intuitively plaubased on parabolic regions as part boundaries has also sible approach reduces the problem to that of fitting been formulated. ~57~ This model has, however, many curves of least energy. ~69'7°) However, this approach limitations and often provides solutions which do not can be applied only in simple cases, and cannot handle multiple joints [Fig. 18(d)]. coincide with intuition (e.g. over-segmentation). Hoffman et al. (39'58) extended Koenderink's ideas and Our approach to shape decomposition can be summarsuggested that shape part boundaries should ideally be ized as involving the following three stages: defined by contours of negative extrema of (principal) 1. Identify points belonging to part boundary regions. curvatures. Further studies in this direction are presented 2. Construct part boundaries as smooth closed curves in references (59-63). The inference of three-dimen(cut lines) which contain these points. sional shapes based on similar principles was analyzed 3. Separate ("cut apart") the parts of the shape along in references (62,64,65). these cut lines. We can extend the approach of Hoffman and Richards by characterizing both part boundaries and the parts It is clear that these stages do not necessarily yield themselves in terms of surface curvature according to unique results; hence these problems, like many others in the following rules: 1. Local extrema of a negative principal curvature (such as hyperbolic and concave parabolic points) are associated with boundaries of parts. 2. Convex elliptic and parabolic regions correspond to parts ("protuberances"). 3. Concave elliptic and parabolic regions correspond to negative parts ("dimples"). Figure 18 illustrates simple examples of shape parts. Here (a) shows a convex elliptic region surrounded by a closed contour of hyperbolic points and (b) shows a concave elliptic region surrounded by a closed contour of hyperbolic points. Note, however, that the contours of hyperbolic points may have branches (d) and they may not even be closed (c). These problems were also noticed

Parabolic Elliptic l l

Hyperbolic

Fig. 18. The use of surface curvature to define pa~s and boundaries (see text).

830

S. FEJES and A. ROSENFELD tion (5) and equation (6) and leads to part decompositions of shapes. In this framework the energy function will be represented by the energy of the migrating set, and the constraints which are implicitly involved in shape decomposition operator EA will be defined by the part boundary and part-like regions of the shape using different fixed point and hypersurface constraints. In effect, this will break up the operator -ZA into a sequence of operators E~,I, E~,2, "=~,3. The implementation of stages (A), (B) and (C) as an initialization step and a sequence of migration processes is as follows.

Decomposition steps using the MP model /

Fig. 19. The classification of surface points of a shape according to their principal curvature.

early vision, tend to be ill-posed and need to be regularized [see reference (22)]. j j

4.4.2. A regularization framework using the MP model. In order to obtain unique results at stages (B) and (C) we need to recast them into a regularization framework, involving variational principles, which introduces restrictions on the allowable solutions (see also Section 3). We now outline a general shape decomposition framework based on such principles. Let us define a (universal) shape decomposition operator E~ : S -~ {Si},

(5)

where S - Uicl Si denotes the shape to be decomposed; Si = v A(S)i is the ith part obtained by the decomposition, so that Si A Sj = 0 (iCj), and A is a scale parameter. In order to restrict the allowable decompositions, we also define a cost (energy) function • whose minimum is the unique solution:

U ffJ(Pk) >_ Z k

~(U'a(S)i)'

(6)

i

where UkPk = S and Pt fq Pm -- ~ (lCm). In words: shape decomposition operator --a provides the set of parts such that the total energy assigned to the composition of these parts is minimal. In earlier sections we have shown how the MP model can be applied to various optimization problems. We now construct a sequence of three migration processes which implements the decomposition model defined by equa-

IJStage (A) does not require regularization. It can be performed using any surface curvature estimation algorithm. It will be demonstrated in Section 4.4.3 that due to our use of a regularization framework, it is not even necessary to provide the complete set of these points, but only a small number of them.

0. Find some of the points of each part boundary region on the surface of the shape, and initialize connected sets on the surface of the shape which contain these points but do not intersect any part-like region. 1. Migrate these sets, subject to the hypersurface constraint represented by the shape itself and the fixed points found in step (0), until a stationary state is approached. This provides a set of closed cut lines which contain the points found in step (0). These lines are smooth because of the curvature-dependent component of the migration, and they form closed (local) geodesics on the surface of the shape. 2. Migrate the set of cut lines until they vanish. The loci of this migration generate a set of surface patches. 3. Using the cut lines as sets of fixed points, migrate the set of surface patches until a stationary state is approached. This process generates a system of minimal surface patches defined by the set of closed cut lines as boundary curves. The set of parts is then obtained by using this set of minimal surface patches as part separators. These steps are illustrated in Fig. 20.12

4.4.3. Experimental results. This section gives some experimental results of our part decomposition scheme using synthetic images (Fig. 21). These preliminary examples demonstrate that our shape decomposition method using the MP model yields intuitively satisfying results. Our method combines the characteristics of the MP model, as a regularization tool, with the intuitive rules of the decomposition model. The resulting system unifies the role of boundary curvature as key to locating part boundaries with the intuitive notions of symmetry seeking and minimality of contours and surfaces (extrema). We believe that a general solution to shape decomposition, if such a solution exists, will have to incorporate these elementary principles and natural characteristics in essentially this way. We should note that our boundary-based decomposition model has limitations. Flat shapes and thick shapes having similar local surface curvature maps are difficult

12For simplicity this figure shows an example which does not involve a natural part decomposition; it is intended only to illustrate the sequence of steps.

Discrete active models and applications

831

11, f

11, f

Fig. 20. Shape decomposition using migration process: (a) placement of some fixed points at part boundary regions; (b) initialization of connected sets on the surface of the shape; (c) migration of these sets subject to hypersurface and fixed point constraints until a stationary state is approached which yields sets of cut lines; (d) generation of sets of surface patches defined by the loci of the free migration of the sets of cut lines; (e) migration of the surface patches with the sets of cut lines fixed until stationary states are approached yielding sets of minimal surface patches; (f) generation of parts by using the set of minimal surfaces as part separators. to differentiate in our approach. Concave parabolic regions, which belong to part boundary regions by definition, can also represent simple bends in flattened shapes, which do not necessarily give rise to intuitive parts, as they do in the case of thick shapes [see also the model of the "shape triangle"~47)].

In our examples, the initialization of the migrating sets in step 2 was performed manually. The problem of automatically initializing these sets has not yet been solved. The limitation of our approach discussed in the preceding paragraph suggests that initialization could take advantage of techniques which do not suffer from

832

S. FEJES and A. R O S E N F E L D

i~.

!N,

P ••~

"

~

4

! P

Fig. 21. Shape decomposition using migration processes (each row presents an example).

these limitations (e.g. skeleton-based methods). Such a hybrid approach could combine different decomposition principles in a natural way, leading to more general decomposition schemes.

5. CONCLUSIONS AND FUTURE WORK

The concept of modeling a complex, distributed system as a multitude of individual, self-organizing particles is not new. When Turing ~Tn)introduced reaction

Discrete active models and applications diffusionl3--a mathematical model incorporating an array of cells as elementary reaction sites--he was primarily motivated by the need to characterize processes in which chemicals react with and diffuse into each other. A similar principle was embodied in cellular automata m o d e l s ( U ) ~ i s c r e t e - t i m e - b a s e d systems in which the interactions of cells are defined by logic-based rules. Although in general each individual particle (cell) in such a system obeys only simple, locally applicable rules, the multiplicity of particles gives rise to a highly complex autonomous system. The incorporation of the local rules imposes restrictions not only on the local, but evidently also on the global, characteristics of the system. This is because of the iterative processes involved, which propagate the local rules (constraints) in both the time and the space domain of the system. This is a general mechanism of optimization in which, as the system becomes stationary, its state exhibits an extremum--an equilibrium. In this paper the simple, yet novel MP model was presented. In this model the interaction of each particle in the migrating system is defined by local (and in contrast to previous approaches; geometry-based) rules involving averaging. It was demonstrated that MPs, which are geometric diffusion-based iterative processes, can be successfully applied to solving a number of highly sophisticated problems involving optimization. We attempted to present a self-contained treatment of this class of processes, but several possible extensions still suggest themselves, as outlined in the following paragraphs. We have focused our attention on MPs which involve bi-valued ("segmented") data. The extension of the averaging-based model to gray level data is straightforward. Such an extension will provide more flexibility to our approach. The current MP model incorporates only fixed point and hypersurface constraints which are all "hard" constraints. Many applications could benefit from the use of "soft" constraints, i.e. constraints which have more control flexibility. Such constraints can be easily incorporated into the MP model, e.g. by introducing a scaling factor that attenuates the value of the displacement of each point. This extension can also be directly applied to gray scale (i.e. unsegmented) data. Our MP model, which involves averaging over equigeodesic neighborhoods, represents a migrating set in which each particle exerts an attraction force on its neighborhood. We can easily extend this model to include a repulsion force, also defined on a geometric basis. This extension could be highly beneficial, since it would provide a migrating set with the capability of controllable splitting. This may provide a potential solution to the problem of automatic initialization of migrating sets in the shape decomposition application described in Section 4.4. If such a repulsion force is utilized, the initialized migrating set can simply be defined as the closed surface of the shape to be decomposed. Activating the

833

repulsion component at part-like regions gives rise to holes in the set, which produce modifications to the topology of the migrating set when required. Finally, we note that since the MP model--similarly to other, active-particle-based approaches--involves simple, local operations which can be performed independently of each other, efficient, massively parallel implementation of the model is possible. This paper has demonstrated some of the potential of approaches involving systems of self-organizing active particles. However, it seems clear that much of this potential has not yet been exploited. Although diverse applications based on similar principles in the physics and computational vision literature (m'24'73"74) are known, there are still many areas where these approaches could be beneficial. We believe that very simple models like the MP model, which incorporate local interaction rules based on implicit formulations (e.g. on geometry), can be successfully applied in such areas. APPENDIX A. THE LOCUS OF A MIGRATING SPACE CURVE

Let the surface patch C~ : SI x [0, 1 ) ~

~m be defined

by O ~ ( s , t ) _ 02~(s,t) Ot Os 2 '

C~(s,0) = C~(s),

(AI)

where the initial curve ~6~(s)is the boundary condition, s is the arc length parameter of the curve, and t denotes time. Let us introduce the notation Xy for the partial derivative of X with respect to y, i.e. OX/Oy, and similarly for higher order derivatives Xyy .~ OZX/Oy 2. For simplicity, we omit the parameterization of the functions. We denote the scalar product (dot product) of vectors a and b by

(a,b). In this notation, equation (A1) becomes c~. _ ~d,,~, ~(0) = ~6~.

(A2)

Consider the second fundamental form (31) of surface patch ~(s), which describes the mean curvature: H

1 GD - 2FD' + ED" ~ (tq + ~2) = EG - F 2 '

(A3)

where E-

(g's, g's},

D = - ( U , g',.~),

G = (crV,,or"t),

D' =

F = (~s, c~t},

D"= -(N,%,),

(N, c~o~t),

(A4)

where N is the unit normal, H is the mean curvature, and eq and ~2 are the principal curvatures. Since ~, is parameterized by its arc length s, it follows (37) that (~,, ~ed = 1. Calculating the partial derivative of this expression with respect to s, we obtain 2 ¢ e s , " e , s ) = o.

J3Notice the "accidental" coincidence between this terminology and that related to wave propagation. (72)

From equation (A2) it follows that (cgs, %) = 0; hence the scalar F in equation (A4) is zero.

834

S. FEJES and A. ROSENFELD

The unit surface normal N can be expressed in terms o f the vector product o f the first partial derivatives (surface tangents) (75) as N -- -

I~'~

-

(A5)

× ~e,l

20.

Hence from equation (A2) it follows that N is perpendicular to tess, so that scalar D in equation (A4) is also zero. We thus have

<~, d.,)(N, CO'u)

O" H--

19.

(N,

cdu)

--

21.

22. 23.

where we have used the fact that E = ( ~ s , ~ s ) ¢ 0. Evidently, the numerator o f H is not necessarily zero. We conclude that the surface patch ~ ( s , t) is not necessarily a minimal surface.

24.

25. REFERENCES

26. 1. M. Bertero, T. A. Poggio and V. Torre, Ill-posed problems in early vision, Proc. IEEE 76, 869-889 (1988). 2. M. Kass, A. Witkin and D. Terzopoulos, Snakes: Active contour models, Int. J. Comput. Vision 1, 321-331 (1988). 3. V. Caselles, E Catte, T. Coll and E Dibos, A geometrical model for active contours in image processing, Numer Math. 66, 1-31 (1993). 4. S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, Z Comput. Phys. 79, 12-49 (1988). 5. S. Fejes, Migration processes: Theory and applications, PhD. Thesis, Technical Report CS-TR-3603, University of Maryland, College Park, Maryland (December 1995). 6. S. Fejes and A. Rosenfeld, Migration processes, J. Math. Imaging Vision (to appear). 7. M. A. Grayson, Shortening embedded curves, Ann Math. 129, 71-111 (1989). 8. B. B. Kimia, A. Tannenbaum and S. W. Zucker, On the evolution of curves via a function of curvature. I. The classical case, J. Math. Anal. Appl. 163, 438-458 (1992). 9. B. B. Kimia, A. R. Tannenbaum and S. W. Zucker, Shapes, shocks, and deformations. I: The components of twodimensional shape and the reaction-diffusion space, Int. J. Comput. Vision 15, 189-224 (1995). 10. S. E Thompson, Growth models for shapes, Technical Report CAR-TR-743, University of Maryland, College Park, Maryland (December 1994). 11. T. Toffoli and N. Margolus, Cellular Automata Machines. MIT Press, Cambridge, Massachusetts (1987). 12. M. Gage and R. S. Hamilton, The heat equation shrinking convex plane curves, J. Differential Geom. 23, 69-96 (1986). 13. M. A. Grayson, The heat equation shrinks embedded plane curves to round points, J. Differential Geom. 26, 285-314 (1987). 14. C. L. Epstein and M. Gage, The curve shortening flow, in Wave Motion, A.J. Chorin and A.J. Majda, eds, pp. 15-59. Springer, Berlin (1987). 15. G. Huisken, Flow by mean curvature of convex surfaces into spheres, J. Differential Geom. 20, 237-266 (1984). 16. S. Angenent, Parabolic equation for curves on surfaces, Part II. Intersections, blow-up and generalized solutions, Ann. Math 133, 171-215 (1991). 17. G. Sapiro and A. Tannenbaum, Area and length preserving geometric invariant scale-spaces, IEEE Trans. Pattern Analysis Mach. IntelL 17, 67-72 (1995). 18. B. B. Kimia, A. Tannenbaum and S. W. Zucker, Toward a computational theory of shape: An overview, Technical

27.

28.

29.

30. 31.

32. 33.

34.

35. 36. 37. 38. 39. 40. 41.

42.

Report TR-CIM-89-13, McGill University, Montreal, Canada (June 1989). P. Neskovic and B. B. Kimia, Three-dimensional shape representation from curvature dependent surface evolution, Proc. Int. Conf. Image Process., Austin, Texas, pp. 6-10 (November 1994). W. L. Briggs, A Multigrid Tutorial, 1st edn. SIAM, Providence, Rhode Island (1987). D. Terzopoulos, Image analysis using multigrid relaxation methods, IEEE Trans. Pattern Analysis Mach. lntell. 8, 129-139 (1986). A. N. Tichonov and V. J. Arsenin, Solution of Ill-Posed Problems. Winston, New York (1977). L. D. Cohen, On active contour models and balloons, Comput. Vision Graphics Image Process. Image Understanding 53, 211-218 (1991). R. Szeliski, D. Tonnesen and D. Terzopoulos, Modelling surfaces of arbitrary topology with dynamic particles, Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognition, New York, pp. 82-87 (June 1993). H. Tek and B. B. Kimia, Shock-based reaction-diffusion bubbles for image segmentation, Technical Report LEMS138, Brown University, Providence, Rhode Island (1994). V. Caselles, R. Kimmel, G. Sapiro and C. Sbert, Minimal surfaces: A three-dimensional segmentation approach, Technical Report EE 973, Technion, Haifa, Israel (June 1995). V. Caselles, R. Kimmel and G. Sapiro, Geodesic active contours, Proc. Int. Conf. Comput. Vision, Cambridge, Massachusetts, pp. 694-699 (June 1995). W. Neuenschwander, E Fua, G. Sz~kely and O. Kiibler, Making snakes converge from minimal initialization, Proc. Int. Conf. Pattern Recognition, Jerusalem, Israel, pp. 613615 (October 1994). W. Neuenschwander, E Fua, G. Szrkely and O. Kiibler, Deformable Velcro surfaces, Proc. Int. Conf. Comput. Vision, Cambridge, Massachusetts, pp. 828-833 (June 1995). A. Blake and A. Zisserman, Visual Reconstruction. MIT Press, Cambridge, Massachusetts (1987). E. A. Lord and C. B. Wilson, The Mathematical Description of Shape and Form. Ellis-Horwood, Chichester, U.K. (1984). J. L. Troutman, Variational Calculus with Elementary Convexity. Springer, Berlin (1983). N. Kiryati and Sz, 6kely, G. Estimating shortest paths and minimal distances on digitized three-dimensional surfaces, Pattern Recognition 26, 1623-1637 (1993). R. Kimmel, A. Amir and A. M. Bruckstein, Finding shortest paths on surfaces using level sets propagation, IEEE Trans. Pattern Analysis Mach. lntell. 17, 635-639 (1995). S. Krikpatrick, C. D. Gellatt and M. P. Vecchi, Optimization by simulated annealing, IBM Thomas G. Watson Research Center, Yorktown Heights, New York (1982). G. A. Korn and T. M. Korn, Mathematical Handbook for Scientists and Engineers, 2nd edn. McGraw-Hill, New York (1968). H. W. Guggenheimer, Differential Geometry. McGrawHill, New York (1963). F. P. Preparata and M. 1. Shamos, Computational Geometry: An Introduction. Springer, Berlin (1985). D. D. Hoffman and W. A. Richards, Parts of recognition, Cognition 18, 65-96 (1984). A. Pentland, Recognition by parts, Technical Report 406, SRI International, Menlo Park, California (1986). A. Gupta and R. Bajcsy, Part description and segmentation using contour surface and volumetric primitives, SP1E Sensing and Reconstruction of Three-Dimensional Objects and Scenes 1260, 203-214 (1990). E Solina and R. Bajcsy, Recovery of parametric models from range images: The case tor superquadrics with global deformations, IEEE Trans. Pattern Analysis Mach. lntell. 12, 131-147 (1990).

Discrete active models and applications 43. T. Binford, Visual perception by computer, Proc. IEEE Conf. on Systems and Control, Miami, Florida (December 1971). 44. T. Fan, G. Medioni and R. Nevatia, Recognizing 3-d objects using surface descriptions, IEEE Trans. Pattern Analysis Mach. lntell. 11, 1140-1157 (1989). 45. E Ulupinar and R. Nevatia, Shape from contour: Straight homogeneous generalized cylinders and constant cross section generalized cylinders, IEEE Trans. Pattern Analysis Mach. Intell. 17, 120-135 (1995). 46. D. Terzopoulos and D. Metaxas, Dynamic 3D models with local ~nd global deformations: deformable superquadrics, 1EEE@rans. Pattern Analysis Mach. lntell. 13, 703 714 (1991). 47. B. B. Kimia, A. Tannenbaum and S. W. Zucker, On the shape triangle, in Visual Form: Analysis and Recognition, C. Arcelli, L. P. Cordelia and G. Sanniti di Baja, eds, pp. 307-323. Plenum Press, New York (1992). 48. K. Siddiqi and B. B. Kimia, Parts of visual form: Computational aspects, Technical Report LEMS114, Brown University, Providence, Rhode Island (1994). 49. K. Siddiqi, K. Tresness and B. B. Kimia, Parts of visual form: Ecological and psychophysical aspects, Technical Report LEMS-104, Brown University, Providence, Rhode Island (1994). 50. I. Pitas and A. N. Venetsanopoulos, Morphological shape decomposition, IEEE Trans. Pattern Analysis Mach. Intell. 12, 3 8 4 5 (1990). 51. H. Blum, A transformation for extracting new descriptors of shape, in Models for the Perception of Speech and Visual Form, W. Wathen-Dunn, ed., pp. 362-380. MIT Press, Cambridge, Massachusetts (1967). 52. H. Blum and R. N. Nagel, Shape description using weighted symmetric axis features, Pattern Recognition 10, 167-180 (1978). 53. H. Rom and G. Medioni, Hierarchical decomposition and axial shape description, 1EEE Trans. Pattern Analysis Mach. Intell. 15, 973-981 (1993). 54. L. R. Nackman and S. M. Pizer, Three-dimensional shape description using the symmetric axis transform I: Theory, IEEE Trans. Pattern Analysis Mach. Intell. 7, 187-202 (1985). 55. J. J. Koenderink and A. J. van Doorn, Photometric invariants related to solid shape, Opt. Acta 27, 981-996 (1980). 56. J. J. Koenderink and A. J. van Doom, The shape of smooth objects and the way contours end, Perception 11, 129-137 (1982).

835

57. H. Rom and G. Modioni, Part decomposition and description of 3D shapes, Proc. Int. Conf. on Pattern Recognition, Jerusalem, Israel, pp. 629-632 (October 1994). 58. B. M. Bennett and D. D. Hoffman, Shape decompositions for visual recognition: the role of traversality, in Image Understanding 1985-1986, W. Richards and S. Ullman, eds, pp. 215-256. Ablex, Norwell, Massachusetts (1989). 59. W. Richards and D. D. Hoffman, Codon constraints on closed 2D shapes, Comput. Vision Graphics Image Process. 31, 265-281 (1985). 60. M. Leyton, Symmetry-curvature duality, Comput. Vision Graphics Image Process. 38, 327-341 (1987). 61. M. Leyton, A process-grammar for shape, Artificial lntell. 34, 213-247 (1988). 62. W. A. Richards, J. J. Koenderink and D. D. Hoffman, Inferring three-dimensional shapes from two-dimensional silhouettes, J. Opt. Soc. Am. 4, 1168-1175 (1987). 63. J. M. H. Beusmans, D. D. Hoffman and B. M. Bennett, Description of solid shape and its inference from occluding contours, J. Opt. Soc. Am. 4, 1155-1167 (1987). 64. J. J. Koenderink, What does the occluding contour tell us about solid shape, Perception 13, 321-330 (1984). 65. L. B. Wolff and J. Fan, Curvature segmentation from multiple illumination using a photometric invariant, Proc. ARPA Image Understanding Workshop, Monterey, California, pp. 1547-1554 (November 1994). 66. R. M. Haralick, L. T. Watson and T. J. Laffey. The topographic primal sketch, Int. J. Robot. Res. 2(1), 50-72 (1983). 67. J. J. Koenderink, Solid Shape. MIT Press, Cambridge, Massachusetts (1990). 68. D. W. Jacobs, Grouping for Recognition, AI Memo 1177. MIT, Cambridge, Massachusetts (November 1989). 69. B. K. P. Horn, The Curve of Least Energy, AI Memo 612. MIT, Cambridge, Massachusetts (January 1981). 70. S. Ullman, Filling-in the gaps: The shape of subjective contours and a model for their generation, Biol. Cybernetics 69, 1-6 (1976). 71. A. Taring, The chemical basis of morphogenesis, Phil. Trans. Royal Soc. Sen B 237, 37-72 (1952). 72. J. Smoller, Shock Waves and Reaction--Diffusion Equations. Springer, Berlin (1983). 73. A. J. Chorin, Flame advection and propagation algorithms, J. Comput. Phys. 35, 1-11 (1980). 74. P. Chambers and A. Rockwood, Visualization of solid reaction-diffusion systems, IEEE Comput. Graph. Appl. 15(5), 7-11 (1995). 75. A. Gray, Modern Differential Geometry o[ Curves and Surfaces. CRC Press, Boca Raton, Florida (1993).

About the Author - - AZRIEL ROSENFELD is a tenured Research Professor, a Distinguished University Professor, and Director of the Center for Automation Research, at the University of Maryland in College Park; he also holds Affiliate Professorships in the Departments of Computer Science and Psychology and in the College of Engineering. He holds a Ph.D. in Mathematics from Columbia University (1957), Rabbinic Ordination (1952) and a Doctor of Hebrew Literature degree (1955) from Yeshiva University, and honorary Doctor of Technology degrees from Linkoping University, Sweden (1980) and Oulu University, Finland (1994). He is a widely known researcher in the field of computer image analysis, has published over 25 books and over 500 book chapters and journal articles, is an Associate Editor of over 25 journals, and has directed over 50 Ph.D. dissertations. He is a Fellow of the IEEE (1971) and the Washington Academy of Sciences (1988); a founding Fellow of the AAAI (1990), the ACM (1994) and the IAPR (1994), and has won numerous professional society awards.

About the A u t h o r - - S A N D O R FEJES was born in Szeged, Hungary, on 26 June 1968. He received his M.S.

degree in Electrical Engineering in 1992 at the Technical University of Budapest, Hungary. As a graduate student he was with the KFKI Research Institute of Measurement and Computing Techniques of the Hungarian Academy of Sciences. In 1994 he visited the Center for Automation Research, University of Maryland at College Park as a Fulbright Scholar, where he completed his Ph.D. dissertation in 1995. Since then he has been a postdoctoral faculty member at Maryland. His research interests are in real-time image processing and computer vision, discrete geometry, artificial intelligence, and parallel architectures.