I. Introduction
1.3 Organization of this thesis
Up to now, we have mentioned some traditional image morphing approaches. In this paper, we will introduce mathematical morphology, level set methods, their implementations, discuss experimental results, and make our conclusion.
In chapter 2, the contents include mathematical morphological fundamentals, distance transformations, and morphological interpolation methods. Dilation and erosion are two basic operations in mathematical morphology. Distance transformations play important roles in mathematical interpolations. Morphological interpolations contain Hausdroff distance interpolation, proposed by Serra [14] [15], Distance-based interpolation function, proposed by Meyer [16], Median interpolation, proposed by Beucher [19], and Shape-based morphological interpolation, proposed by Udupa [20].
In chapter 3, the contents include the theory of level set methods and the application of them. Osher and Sethian introduce level sets for computing moving wave front. Since that time, level set methods have been proved useful for curve and surface reconstruction, shape reconstruction, image processing, and geometric modeling. Especially, level sets have been proven a robust morphing approach.
Moreover, we will affix a blending (morphing) method proposed by Whitaker [32]
[33] and a radial basis function using level set ideas proposed by Mullan, et. Al [37].
In chapter 4, after introducing mathematical morphology and level sets in details, we will implement those methods step by step. There, we will present our
results with figures. In chapter 5, we will make a conclusion about this paper.
Chapter 2
Mathematical Morphology
To attain the goal that image morphing should be satisfied with the human perception, mathematical morphology and level set methods both provide good theoretical frameworks for shape morphing. Mathematical morphology is a powerful methodology for the quantitative analysis of geometrical structures. Mathematical morphology [8] [9] was initiated in the last 1960s by Matheron and Serra at the Fontainebleau School of Mines in France. Originally, it was applied to analyze images from geological or biological specimens. In chapter 2, the basic theory and the know-how of mathematical morphological interpolations will be presented.
2.1 Fundamentals
Dilation and Erosion are two basic morphologic operations. They are close related to the Minkowski set addition and subtraction.
2.1.1 Minkowski Algebra
The Minkowski subtraction of B from A , written as A B0 , is defined by
In here, B is usually called a structuring element.
Figure2-2. The upper left is the original set A , the right one is the structuring element S and show S∨ in the same time; the lower left is after dilating A , the right one is after eroding A
Another example for binary dilation and erosion with structuring element being the 3 by 3 square is shown in Figure 4.
Figure2-3. (a) Original image
for all grayscale images f in ( )F E .
2.2 Distance Transformation and Morphology
Distance transformations play major roles in the morphological interpolations.
A distance transformation [10] converts a digital binary image into a graytone image consisting of feature and non-feature pixels. A few of common distance transformations will be presented here. city-block distance does not allow traveling the distance between center and corner directly. In this case, all pixels with city-block distances from a less than or equal to some radius r form a diamond centered at a. The pixels b with D a b4( , ) 1= are the 4-neighbors of the center a, as shown in the following table.
2 1 2
Table 2-1 (a) Table 2-1 (b) chessboard distances allow a vector to go to its four corners at one step. In this case, all pixels with chessboard distances from a less than or equal to some radius r form a square centered at a. The vector b with D a b8( , ) 1= are the 8-neighbors of
where a and b are vectors in R with n a=( ,..., )a1 an and b=( ,..., )b1 bn .The Euclidean distance can easily calculated by using Pythagorus' theorem. a2+b2 = . c2 This is the usual metrics widely used in many areas. The disadvantages of it are the results are not always correct and the process is slow in some cases. Due to the complex feature geometry errors can be occurred, but they are always small.
Table2-3. Euclidean distance
◆ Chamfer distance (weighted distance)
Gunilla Borgefors [12][13] original chafer-distance method calculates distances of non-object pixels from the nearest object pixel. Borgefors reviews a number of metrics in 2 and 3 dimensions. Chamfer distance transformations are produced in two raster scans over the image. In the forward scan, the mask starts in the upper left corner of the picture, moves from left to right and from top to bottom.
In the backward scan, it starts in the lower right corner, moves from right to left and from bottom to top. The local distances, d1 and d2, in the mask pixels are added to the pixel values in the distance map and the new value of the zero pixel is the minimum of the five sums.
Table2-4. (a) Table2-4. (b)
To show in equations as:
Forward:
, min( 1, 1 2, 1, 1, 1, 1 2, , 1 1, , )
i j i j i j i j i j i j
v = v− − +d v− +d v− + +d v − +d v (2.10)
Backward:
, min( , , 1, 1, 1, 1 2, , 1 1, 1, 1 2)
i j i j i j i j i j i j
v = v v+ +d v− + +d v + +d v+ + +d (2.11)
For instance, for d1= and1 d2 = ∞ , we have the city block metric; for d1 =d2 = , 1 the chess board metric; for d1= and1 d2 = 2, Montanari's metric ; and for d1 = 2 and d2 = , Barrow's approximation; the most common is 3-4 chamfer. 3
◆Hausdorff distance
In generally, the “distance” means the shortest path, if a point a∈A is said to be at distance d to set B , we usually assume that d is the distance from a to the nearest point of b∈B. Regularly, this is called a minimin function, because the distance 0 between a single point a and an object (compact set) B 1 is defined as:
( , ) inf{ ( , ) |d a B = d a b b∈B} (2.12)
And this distance equation can also written by using the morphological notion of dilation, such as:
( , ) inf{ |d a B = λ a∈Dλ( )}B (2.13) distance is the maximum distance of a set to the nearest point in the other set. More formally, Hausdorff distance between sets A and B is defined as metric between these points. Changing the last equation to morphological notation:
( , ) inf{ρ A B = λ : A⊆Dλ( ) ,B B⊆Dλ( )}A (2.17)
2.3 Morphological Interpolation
To go into particulars, we will discuss some chief interpolating methods in the following :
• Hausdorff distance interpolation(Jean Serra)
• Distanced-based interpolation function (Frenand Meyer)
• Median interpolation (S.Beucher)
• Shape-based morphological interpolation (J. K. Udapa)
2.3.1 Hausdorff distance interpolation [14] [15] (Jean Serra)
The theoretical background of this interpolation method is the Hausdorff distance, and this method was first proposed by Serra in 1994. Hausdorff distance interpolation has one important disadvantage - the resulting interpolated objects are relatively large - even considerably larger than the input ones.The interpolator is a function which produces the interpolated objects between images in the interpolation sequence. It has three arguments: two input objects and the interpolation level. The interpolation level, denoted by α (a real number 0≤ ≤α 1) indicates a position of the interpolated object in the interpolation sequence. If α =0, the interpolated object equals the initial object, denoted by P . If α =1, the interpolated object equals the final object, denoted by Q . Then, we denote an intermedium of interpolation level 0≤ ≤α 1 by the following statement :
(1 )
( ) ( ) , 0 1
Zα =Dαρ P ID −α ρ Q ≤ ≤α (2.18)
where set Zα turns out to be the intersection of dilates of P and of Q by the structuring element of radii αρ and (1−α ρ) respectively.
Figure2-4 (a) source image Figure2-4 (b) Z0.25 Figure2-4 (c) Z0.5
Figure2-4 (d) Z0.75 Figure2-4 (e) target image
Figure2-4.(a) to (e) The example of Hausdorff distance interpolation
2.3.2 Distance-based interpolation (interpolation function)
“Distance-based interpolation” also named “interpolation function” was first proposed by Frenand Meyer [16]. This is a method of the binary morphological interpolation which makes use of interpolation function. It is based on the function which describes the relative distance between the objects and can be applied to binary and mosaic images. This method has a simple constraint: the two input objects must have non-empty intersection.
Let using P denote the source image and Q denote the target image, see Figure2-5 (a). We begin with overlapping the two images to get their intersection,
denoted by R P Q= I . And we will use the intersection R as a temporary target. It means that we will begin with the source image P and transform it into intersection R while on the other hand begin with intersection R and transform it into the target image Q .
Let P and Q be two sets (objects), (R=PIQ)⊂P. Interpolation function of P and Q is obtained from two base geodesic distance functions d and 1 d . 2 Distance function d describes a distance from the complement of P , denotes P1 ∨, in the area P R/ , and the distance function d describes a distance from 2 R in the area P R/ where P R/ means the area in P but not in R . An example is shown in the following figure.
Figure2-5. (a)
Figure2-5. (b) Figure2-5. (c)
In Figure2-5. (a) the ellipse presents P , the horizontal shape presents Q , and the coloration presents P QI ; (b) shows d ; (c) shows 1 d 2 constant, the boundaries of the interpolation sets are all arranging to lines, as shown in the Figure2-6 (f).
Figure2-6 (a) source P Figure2-6 (b) target Q
Figure2-7 (a) din source P Figure2-7 (b) Thr0.5(intPQ=d)in source P
In the interpolation step, the start image P is transformed to the target image Q . The action can be separated into two small actions. The first step begins with P shrinking to P QI and the second one begins with P QI growing to Q . Those two steps can run at the same time. The interpolation set Z at distance k k from P and (1−k) from Q can be represented by the union of those two terms:
(intP ) 1 (intQ )
k k P Q k P Q
Z =Thr I UThr− I (2.21)
Figure 2-8 (a) P Figure 2-8 (b) Z0.25 Figure 2-8 (c) Z0.5
Figure 2-8 (d) Z0.75 Figure 2-8 (e) Q
Figure 2-8.(a) to (e) Show the step by step that P morphing to Q ;
2.3.3 Median Set interpolation (J.Serra, S.Beucher)
Morphological median is one of the useful morphological tools for image interpolation. The basic theory was primarily developed for sets. Median interpolation [19] is applicable to binary, mosaic, graytone and color images.
◎ Binary image way, the influence zone defined by equation(2.22) is also represents a median set between two sets, so we can get
( , )M P Q =IZ QP( ) (2.23)
Morphological median of nested sets is defined using dilation and erosion, and it has been proved that the median set satisfies the following equation
( , )M P Q {(Q B) (P B)} size λ , both with the elementary structuring element B
More generally, we extend the simple case Q⊂ to the partial inclusion of P of the second one. It is produced by the iterative generation of new medians:
1. First iteration:
3. M5/8=M(M1/2 , M3/4) 4. M7/8=M(M3/4 ,Y) 4. …and so on…
Figure2-10 (a) Figure2-10 (b) Figure2-10 (c)
Figure2-10 (d) Figure2-10 (e) Figure2-10 (a) to (e) interpolation sequence
◎ Gray-level image
One can observe the equation that we just talked about the median element turns out to be an increasing mapping of its two operands. Hence it extends to a unique graytone mapping (where the transforms of the cross sections are the cross sections of the transforms). The implementation of this graytone operator is similar to that of the set case. For this reason, we input a pair of initial gray-level images ( , )f g and we can formula it in the following way
( , ) sup{m f g = ∀λ: inf[Dλ(inf( , )),f g Eλ(sup( , ))]}f g (2.26)
where 1, 2,...λ = are increasing values, and the dilation and erosion of object X of a given size λ are defined by, respectively:
( ) ( B( B(... B( )))) performed with non-flat (cylindrical, cone, etc.) structuring element.
We will use the last equation to construction the algorithm of median image generation. So we start from two initial input images ( , )f g and define three working images z w m as 0, 0, 0
Then we compute the iterated value as follows
1
where D and E are respectively dilation and erosion with the non-flat structuring element. The iteration will execute until idempotence and finally
( , ) i , i i 1
m f g =m∞ =m m =m+ (2.31)
where ( , )m f g is a new image of morphological median images of images f and g .
2.3.4 Shape-based interpolation
Shape-based interpolation comes from numerical analysis literature. It was first proposed by J. K. Udapa [20]. Generally, shape-based interpolation algorithms are widespread employed on binary images [21] [22] [23]. These interpolation methods consider shape features extracted from the object sets. We will consider one-to-one object interpolation at the beginning.
Shape-based interpolation converts binary images into distance maps by distance transformation functions such as chamfer or city-block distance template [21].
Another illustration is converting the segmented slice image into gray-level images.
These templates are used to efficiently approximate the shortest distance between the pixel and the contour of the object.
In the following figure, it is part from an object, we estimate the distance-from-boundary values assigned to pixels. We take b as the interval length.
The heavy line indicates the boundary of the object thresholded at 0.
Table2-5. Show the distance transformation values from the boundary, the heavy line is the boundary
The distance function Raya and Udupa used is a version of the city-block distance in the past. City-block distance can be used to calculate this distance efficiently, but it is a relatively bad approximation to Euclidean distance. For this reason, some improvement can be achieved by using chamfer distance. That is obviously better performance of shape-based interpolation by using a distance function that approximates Euclidean distance more closely. In interpolation step, we can calculate these by two consecutive chamfering processes for the inside and then for the outside. One method of interpolating between slices is linear interpolation.
Chapter 3
Level Set Method
In Chapter 2, we have discussed the theory and interpolation method of mathematical morphology completely. In this chapter, we will continue deliberating the other morphing method – level set and discussing its theory and the application by the same token.
The level set method was first introduced by Stanley Osher and James A.
Sethian in 1987 [24] [25] [26]. It has been proven to be a robust method to do morphing. The implicit surfaces are the groundwork of level set method. So it is essential to present the implicit surfaces at the beginning. Besides, there are some brief introductions to associate with this method.
3.1 Implicit functions
In n-dimensional space. There are two ways to define functions, implicitly [26]
[27] [28] and explicitly. Most of the equations we have dealt with have been explicit equations that are forms of y= f x x( ,1 2,...,xn) . In the meantime, the implicit definition of y as a function of x1 to xn is the same as a relation which expressed indirectly by an equation such as f x( ,...,1 x yn, )= , where k k is a function or a constant.
Definition
An equation of the form ( , )f x y = (where f is smooth) often defines a 0
curve in the plane. When this is the case, we say that the curve is defined implicitly (as opposed to explicitly) because we may not be able to solve for y in terms of x. Nevertheless, assuming that a local solution exists, the method of implicit
differentiation often allows us to solve for implicit differentiation ' dy
y = dx as a function of x and y at any point ( , )x y along this curve.
In general, an implicit function defined in 2-dimension space divides the space into two parts, inside region and outside region, and the border between the inside and the outside is called the interface or isocontour. In n-dimensional space, the inside and outside regions are n-dimensional object, while the interface is less than one-dimension. By the way, interface curves are usually limited to be closed.
The following figure is an example of the implicit function.
Figure3-1. The implicit function φ to identify inside region and outside
3.2 Signed distance function
In general, a distance function [26] [29] [30] ( )d xv that there are several ways to calculate distance ( )d xv
such as Euclidean distance, chessboard distance…., it can be referred to chapter 2.2 Distance transform.
A signed distance function is an implicit function φ with | ( ) |φ xv = d x( )v
Also, the opposite definition to negative on the exterior region, positive on the interior region, and zero on the boundary is allowed.
3.3 Level set Formula
After that, level set method using implicit surfaces. The original idea of level set [26] [31] can be traced back to Hamilton - Jacobi approach to numerical solutions of a time-dependent equation for moving implicit surfaces. In other words, level set
method gives an interface Γ in R of codimension one, bounding a region Ω , we n want to analyze and compute its subsequent motion under a velocity field v as the next position. More clearly, the idea behind level set methods is: instead of evolving a curve in a two dimensional plane, which requires additional parameters of the curve, evolve a 2D function (surface) in 3D . And this is a much easier problem.
For analyzing the level set method, we should introduce various basically know-how at first.
( )
Figure3-2. Show the relation of signed mean curvature κ and the appearance
◆ The one dimension Heaviside function is defined by
where H means the Heaviside function just described above.
◆ The Dirac delta function is the directional derivative of the Heaviside
1
or we can express another formula
( ) t 0
H ∇ + = (3.10) φ φ
where F is called Hamilton principle function.
After introducing the all above notations, we will carry on the main equation of level set. In order to avoid problems with instabilities, deformation of surface elements and complicated surgical procedures for topological repair of interfaces, the implicit function φ both represents the interface and the moving interface.
Level set method provides the mathematical and numerical mechanisms for computing surface (image) deformation as time-varying values of φ by solving a partial differential equation on the 3D grid. To describing the movement of the surface, we denote the point of xv∈Γ
lies on the surface and it deforms as d x dt
v . In this case, there are generally two options for representing such surface movement implicitly, those are static form and dynamic form.
★ Static :
The static formulation is represented as F x( )=k t( ), the F function is single while the k is varied.
Figure3-3. The static formulation F x y( , )=k t( )
This equation can be solved efficiently by starting with a single surface using the fast marching method. Nevertheless, the static level set method has a unfavorable problem : it’s motion strictly inward or outward. In other words, the surface can not pass back over itself overtime, and we can observe this phenomenon in the next figure.
Figure3-4 (a) Figure3-4. (b)
Figure 3-4(a) the original shape and (b) the change of the shape and notice that the motion can not pass back
★ Dynamic
On the other side, the dynamic formulation is indicated that ( , )F x t = , the k ( , )
F x t is evolved while the k is fixed.
Figure. 3-5 The dynamic formulation of F x y t( , , )= k
Then evolving φ( , )x tv with embedding time t automatically gives the
evolution. The term vx
This dynamic form provides arbitrarily motions which move forward and backward and cross back over their own paths. Therefore, we using front tracking scheme to solving this equation. Compared with static form, the dynamic form is more flexible, so the remainder of this discuss is focus on the dynamic case.
There is a widespread use that choosing the simple convection ( or advection ) equation
t V 0 φ + ⋅∇ =uv φ
(3.13) to define the evolution of the implicit function where the t subscript denotes a
temporal partial derivative in the time variable t and the Vuv
is a speed function of the velocity of each point on the implicit surface. The level set equation also can rewrite as a well-known form
where V is a constant component of velocity in the normal direction, it also wise n known as the normal velocity. Let us return to our main subject, the level set equation is a kind of Hamilton-Jacobi equation where H(∇ = ⋅∇φ) Vuv φ
or H(∇ =φ) Vn|∇ . φ|
To avoid oscillation, we can use the characteristic function to identify the shapes and boundaries. We calculate the volume integral of a function f over the interior region Ω at first : −
( )f x χ−( )x d x
∫
Ω v v v (3.15)where the region of integral is all of Ω , but the χ− prunes out the exterior region Ω automatically. We rewrite the volume integral of a function f over the interior +
region Ω using the Heaviside function: − the volume integral, the surface integral of the function f over the boundary ∂Ω is defined as :
( ) ( ) f x δ x d x
∫
Ω v v v (3.17)3.4 Deforming Method
In this section, we will talk about some image morphing using level set method.
3.4.1 Image blending
Ross T. Whitaker [32] [33] was proposed a method about the post-processing after image morphing. Once one image has been morphed, it needs a low-level post-processing to calculating the detail and color that are not captured by the coordinate transforms.
In advance, it is said that the image morphing need to input two images- source and target. To produce intermediums of morphing sequence, it was always use two input images in traditional. In the mean while, a convenient skill is proposed that using the new manufactured meddle images as the materials in place of the two original input images.
Figure3-6. (a) Figure3-6. (b)
Figure3-6 (a) (b) show the in-between images manufacture with difference kinds stuff
The two input images are regarded as two continuous functions T:ℜ → ℜ , 2 the domain ℜ denotes the image domain and the range 2 ℜ denotes the set of grayscale values of the image. Let F(x,y) and G x y represent the source and ( , )
The two input images are regarded as two continuous functions T:ℜ → ℜ , 2 the domain ℜ denotes the image domain and the range 2 ℜ denotes the set of grayscale values of the image. Let F(x,y) and G x y represent the source and ( , )