MEASURING AGGLOMERATION OF AGGLOMERATED PARTICLES PICTURES

arXiv:1302.5150v4 [] 30 Jun 2013

SHIGEKI MATSUTANI, YOSHIYUKI SHIMOSAKO
Abstract. In this article, we introduce a novel geometrical index δagg , which is associated with the Euler number and is
obtained by an image processing procedure for a given digital picture of aggregated particles such that δagg exhibits the
degree of the agglomerations of the particles. In the previous work (Matsutani, Shimosako, Wang, Appl.Math.Modeling 37
(2013), 4007-4022), we proposed an algorithm to construct a picture of agglomerated particles as a Monte-Carlo simulation
whose agglomeration degree is controlled by γagg ∈ (0, 1). By applying the image processing procedure to the pictures of the
agglomeration particles constructed following the algorithm, we show that δagg statistically reproduces the agglomeration
parameter γagg . agglomeration, digital image processing procedure, Euler number

1. Introduction

In the article [11], in order to find the agglomeration effect in the electric conductivity of the nano-composite material, we study the electric conductivity in an agglomerated
continuum percolation model and show that the agglomeration of particles affects the macro-material properties. The
purpose of this article is to evaluate the agglomeration in
the binary digital images of agglomerated particles, e.g., of
electron-microscopes.
For the same purpose, so many evaluation methods and
definitions of the agglomeration are proposed to evaluate the
agglomeration. In spatial point analysis, the distribution of
the nearest distance particles, Clark-Evans index and so on
are considered [2, 7]. Further Miles considered the problem
in [12] and showed the two-dimensional overlapping ratio
of the random configurations in three dimensional space.
These investigations on the agglomeration have been done
in the framework of the statistical analysis for a point pattern R = {pi ∈ R2 } which are given as statistical configuS
rations of (finite) points. In the analysis, Rr := p∈R Ur,p
is investigated, which is a configuration of disks whose centers are R, where Uε,p := {q ∈ R2 | |q − p| < ε }. The
Euler number, the area and the perimeter of Rr for several
point processes R’s are computed as morphological indices
or the Minkowski characterization [7]. When R is given by
the point process of the Poisson type, Stoyan, Kendall and
Mecke studied their behaviors based on the study of Miles
[12] and found that
(1)
1
e(x) = (1 − x)e−x , a(x) = (1 − e−x ), ℓ(x) = e−x ,
x

Nano-composite materials have a promising future from
industrial viewpoints, since in the materials, geometrical
properties in micro-scale play crucial roles and generate
novel and various macro-material properties. By controlling the geometrical properties or shapes, we can design the
macro-material properties drastically. Following Kelvin’s
philosophy of science [6]1, it is quite important to evaluate
such geometrical properties or shapes if one needs to control
them using the scientific knowledge.
On the other hand, the smaller the particle is, the larger
the effect of the surface energy is. It means that small particles are apt to aggregate or agglomerate in general because
the agglomeration and the aggregation of the particles decrease the total surface energy and contribute to the stability of the system. When we handle materials consisting of
nano-particles, the agglomeration and the aggregation are
ones of the most important shapes since they sometimes
have an effect on the generations of the macro-materials
properties. In this article, we focus on them. It is said
that the aggregation is due to chemical effects whereas the
agglomeration comes from physical effects. Since in a computational model, there is no difference between them, we
call both agglomeration in this article, though in spatial
point analysis [7], the aggregation is chosen in general.

1

Kelvin wrote his philosophy of science, “In physical science the
first essential step in the direction of learning any subject is to find
principles of numerical reckoning and practicable methods for measuring some quality connected with it. I often say that when you can
measure what you are speaking about, and express it in numbers, you
know something about it; but when you cannot measure it, when you
cannot express it in numbers, your knowledge is of a meagre and unsatisfactory kind; it may be the beginning of knowledge, but you have
scarcely in your thoughts advanced to the state of Science, whatever
the matter may be.”[6]

where e(x), a(x), and ℓ(x) are the normalized versions of
the Euler number, the area and the perimeter of Rr , and x
is a normalized radius r [15, 9]. Mecke and Stoyan studied
the difference among point patterns given by different processes in terms of these behaviors [9]. Further Tscheschel
1

2

and Stoyan also studied the statistical reconstruction of random point patterns [16].
However in the nano-composite materials consisting of
nano-particles, the particles themselves sometimes have complicated shapes, such as ellipsoids and rods, as we investigated in [10]. In other words, the configuration is not given
by a point pattern with radius r in general and further its
parts are overlapped like (d) in Figures 2 and 3. Hence it is
basically an ill-posed problem to define the center points of
actual agglomerated particles in a given picture, e.g., of an
electron-microscopes.
Thus it is natural to consider geometrical properties of
the binary picture as a general geometrical object embedded
in R2 .
Recently MacPherson and Schweinhart proposed a novel
method which evaluates the complicatedness of the complicated geometric objects embedded in a plane R2 in terms
of the persistent homology [8]. The persistent homology
gives the homological quantities of the persistent modules
with real parameter [3, 17]. It could be regarded as a generalization of homotopical approach in traditional algebraic
topology [1], though the deformation does not preserve homotopical properties. For a geometrical object M ⊂ R2 , we
consider a family of objects with a real parameter t ∈ [0, 1],
i.e., {Mt |t ∈ [0, 1]}.
By considering
S union of the ε-neighborhood of each point
in M, Mε := p∈M Uε,p , induced from the standard Euclidean topology, MacPherson and Schweinhart evaluated
the complexity of the geometrical objects. The persistent
homology shows the distributions of topological changes
generated by the persistent modules (vector spaces) induced
from Mt′ ⊂ Mt for t′ < t.
As in [8], to investigate the effect from the standard
topology of Euclidean space and to evaluate the complexity,
we use the one-parameter family of a deformed geometrical
object, and propose a digital image processing procedure
which characterizes the shapes in pictures of the electric
microscope in this article. (In Section 3, we give the list of
assumed geometrical features of the pictures which we deal
with.) For an appropriate geometrical object M in R2 with
a characteristic lengths ℓ1 and S
ℓ2 , we also handle the family
of geometrical objects {Mt = p∈M Ut,p |t ∈ [ℓ1 , ℓ2 ]}. We
define the cumulus of the absolute differential Euler number
(CADE) by,
Z ℓ2
dχ(Mt )
dt,
(2)
E(M; ℓ2 , ℓ1 ) :=
dt
ℓ1

E(M; ℓ2 , ℓ1 ) as Ê(M; ℓ2 , ℓ1 ) for any binary pictures as an
image processing procedure, which is shown in Section 3
more precisely. Further in the nano-materials, there are
several scales and one of them is the size of the particles and
the resolution of the digital picture is given by the pixel size.
We fix ℓ1 and ℓ2 by the pixel size and the (average) radius
ρ of the particle respectively to evaluate the agglomeration
and propose an agglomeration index,
α
(Eˆp(M) − Ê(M; ρ, a)),
(3)
δagg (M) :=
Êp(M)
where p(M) is the volume fraction of M in the region W
(M ⊂ W ⊂ R2 ), Êp is the average of a “standard pattern
of volume fraction p” as mentioned in Section 3, and α is a
normalized factor 1.2.
In order to estimate our agglomeration parameter δagg ,
we performed Monte-Carlo simulations for the binary agglomeration configurations of particles whose degree of the
agglomeration is parameterized by γagg , since in the article [11], we proposed a statistical model which numerically
generates the agglomeration of particles controlled by the
parameter γagg in order to investigate the properties of the
agglomerated continuum percolation models. In Section 2,
we review the algorithm following the article [11]. For a
given parameter γagg ∈ [0, 1], we can statistically construct
the infinitely many configurations with the same level of the
agglomerations. As the continuum percolation model is the
same as a germ-grain model in the study for the point process [7], there are several other algorithms to construct aggregational germ-grain models, such as Neyman-Scott processes [7], though they are different from ours; ours is for the
actual pictures of electron-microscopes of the agglomerated
nano-composite materials with the radius ρ as mentioned in
[11]. We apply the index δagg to evaluate the agglomeration
of the agglomerated configuration which is generated by the
forward method in [11]. Then the relevancy between δagg
and γagg is shown in Section 4, i.e., in Figure 7 and Table
3. We also mention the relation between our δagg and the
well-established Clark-Evans index in Section 4.
2. Agglomerate configuration

In order to explain what is the agglomeration that we are
concerned with, we show the agglomeration configurations
in computer science, which we handled in [11]. In the article
[11], we proposed a construction of the agglomerated continuum percolation models which apparently recover geometric
properties of real nano-particles, though there are several
where χ(X) is the Euler number of X. E(M; ℓ2 , ℓ1 ) evalu- other agglomerated percolation models such as Neymanates how many topology changes occur for the deformation Scott processes [7, 15, 16]. Since our method has a single
[ℓ1 , ℓ2 ].
parameter γagg besides a typical length ρ whereas others
As we are concerned with the image processing procedure are given as point processes with several parameters, we
for images of the electron-microscopes, we will customize believe that ours is a natural model for the actual pictures

3

of electron-microscopes of the agglomerated nano-composite
materials with the radius ρ. As shown in Figures 2 and 3,
we have the agglomerated configuration of particles depending upon a agglomeration parameter γagg ∈ [0, 1]. In this
section, we show the geometrical setting of agglomerated
continuum percolation model in [11], which is modeled by
the agglomerated clusters in nature.
We set particles parameterized by their center positions
(x, y) into a box-region W := [0, L] × [0, L] at random and
get a configuration Mn as a model of continuum percolation. The particle corresponds to a disk with the same
radius ρ, Bxi ,yi := {(x, y) ∈ W | |(x, y)S− (xi , yi )| ≤ ρ}. The
n
configuration Mn is given by Mn := i=1 Bxi ,yi .

∅, we abandon the position and go on to take another uniformly random position (x, y) in W until we find the position
which supplies a connected particle Bx,y with Mn .
In other words, for the case γ ≤ γagg , the added particle must be connected with the previous configuration Mn .
Thus, γagg stands for the agglomeration of the particle system.
By monitoring the total volume fraction which is a function of Mn and is denoted by p(Mn ), we go on to put the
particles as long as p(Mn ) ≤ p for a given volume fraction p. We find the step n(p) such that p(Mn(p)−1 ) ≤ p
and p(Mn(p) ) > p. Since we assume that the difference
between p(Mn(p)−1 ) and p(Mn(p) ) is sufficiently small, we
regard p(Mn(p) ) as the volume fraction p itself hereafter
under this accuracy.
Since in the Monte-Carlo method, we use the pseudorandomness to simulate the random configuration Mn(p)
for given p and γagg , the configuration Mn(p) depends upon
the seed iS of the pseudo-randomness which we choose. We
let it be denoted by Mγagg ,p,iS or its statistical quantity by
Mγagg ,p .

(a)

(b)

Figure 1. The flowchart of the algorithm
which constructs the agglomerated configurations.
The flowchart in Figure 1 illustrates the algorithm. As
the initial state, the configuration M0 has no particle. As
the first step, for a uniform random position (x, y) ∈ W, we
set a particle Bx,y whose center is (x, y) and the radius is
ρ, i.e., M1 := Bx,y .
(c)
(d)
For the (n + 1)-th step, we take a position (x, y) at uniform random in W, and another random parameter γ at
Figure 2. The agglomerated configurauniform random
tions of p = 0.2: These (a), (b), (c), and (d)
S in [0, 1]. If γ is greater than γagg , we set
Mn+1 := Mn Bx,y . We now allow the particles to overlap
show the configurations with the agglomereach other.
ation parameter γagg = 0.0, 0.3, 0.6 and 0.9
For the case γ ≤ γagg , we first check whether the disk
respectively.
Bx,y is connected withTthe previous configuration Mn or
For sufficiently large L and L′ (ρ ≪ L′ < L) and for a
not. For the case Mn S Bx,y 6= ∅, we employ theTposition
′
:= [x, x + L′ ] × [y, y + L′ ] ⊂ W, the volume
and set Mn+1 := Mn Bx,y . Otherwise or Mn Bx,y = window W(x,y)

4

3. Evaluation of Agglomeration for a
configuration Mp,γagg ⊂ W
As we are concerned with the evaluation method as a
digital image processing procedure [13], in this section, we
illustrate our algorithm for a picture which only has binary
values. It is natural that we assume the configuration M
(implicitly Mp,γagg ,iS and a picture of nano-composite material in an electron-microscope) has the following structures:
(a)

(b)

(c)

(d)

Figure 3. The agglomerated configurations of p = 0.4: These (a), (b), (c), and (d)
show the configurations with the agglomeration parameter γagg = 0.0, 0.3, 0.6 and 0.9
respectively.

′
fraction p′ (x, y) in W(x,y)
is proportional to the number of
′
′
particles N (x, y) in W(x,y)
, i.e., there is a constant num′
ber λ such that p (x, y) = λN ′ (x, y). Further Mγagg ,p,iS is
isotropic and independently scattered and thus the asymptotic behavior of Mγagg ,p,iS is a kind of Poisson process [7,
p.66] though we have the condition L′ ≫ ρ.
Further though our algorithm for γagg > 0 is not Markov
process because each step depends on the previous configuration, we should note that it preserves a hierarchical structure for fixed iS and p > p′ ,

(1) L is sufficiently larger than ρ so that the particles
of M are a representative of sufficiently randomized configurations; we could assume the Euclidean
invariance (translation, rotation and inversion) statistically; after averaging them, the physical and geometrical quantities are invariant for any Euclidean
action E(2) up to the statistical deviation. If the
deviation is not small, we could consider the series
of {Mp,γagg ,iS | iS }. (It means that for the case of
the pictures of the electron-microscopes, we could
assume that the researchers prepare the series of
pictures of a material or materials which are produced in the same conditions.)
(2) It is assumed that the volume fraction is less than
the percolation threshold of two dimensional continuum percolation models. (For the case of nanocomposite material which is based upon the percolation theory, the volume fraction around the percolation threshold 0.2 ∼ 0.3 in three dimensional
percolation models is concerned, which is far less
than the percolation threshold of two dimensional
case 0.5 ∼ 0.7.)
(3) There are three sizes of the system or the picture
Mp,γagg (and the digital image of nano-composite
material of an electron-microscope);
(a) the (average) size of particles, which is given
by ρ,
(b) the analyzed size of the system, which is, now,
given by L as mentioned above, and
(c) the pixel size a, which is also controlled so that
we can discriminate the particles in concerned
resolution.

Under these assumptions, we consider geometry
of M.
T It
S
W,
U
is
known
that
the
ε-neighborhood,
M
=
ε,p
ε
⊂ Mγagg ,p,iS .
M
p∈M
can be realized by the so-called level set method in comAs we will show the assumed geometric properties of pic- puter science [14]. Let d : W → R be the signed distance
tures in Section 3, the pictures which we will deal with are from the boundary ∂M so that the outer side is assigned
illustrated in Figures 2 and 3, which are obtained following to the positive distance and the inner side is to the negaour algorithm. Our more concrete aim of this article is to tive one, and then the geometrical object in the level set
as Lt = d−1 (t). Mt of (t > 0) is
recover the parameter γagg for a given configuration Mγagg ,p method can be regarded
S
−1
d ([0, t]) M and Lt = ∂Mt . For t < 0 case,
in a statistical meaning. More precisely, our study is to find equal to S
statistical monotone functions of γagg and to show that one Lt = ∂( p∈∂M Ut,p ) \ L−t . Hence by means of the level
of them is δagg in (3).
set method, we can compute the more precise geometrical
γagg ,p′ ,iS

5

properties beyond the pixel size resolution even on the image defined over a subset of Z2 .
However in the digital image processing procedure, we
investigate the geometrical object up to the pixel size resolution in general. Further we must pay our attentions on
the computational cost if we apply our method to real problems in industry, though level set function method requires
higher computational cost than a simple digital image processing procedure. Hence in this article, we use the thickening scheme in the image processing procedure [13] instead
of the level set function. Though the ordinary thickening
scheme has anisotropic behavior, it does not have a serious
effect on the result because the configuration itself is isotopic or rotational invariant. We use the modified thickening scheme, which improves the anisotropic behavior shown
in Section 4. Let M(i) be the i-th thickening of M in W.
We modify the CADE (2) as an image processing procedure
by
n2
X

Ê(M; n2 a; n1 a) :=

|χ(M(i) ) − χ(M(i−1) )|.

i=n1 +1

(n ,n )

where p(M) is the volume fraction of M, Êp 2 1 is the
average of the standard patterns of volume fraction p, and
α is a normalized factor 1.2, which is chosen as a result of the
comparison with γagg (see Table 3). The standard pattern
means the pattern of γagg = 0 with the same radius in the
(n ,n )
same window W. Then δagg2 1 (Mγagg ,p ) characterizes how
many topological changes occur in the interval (n1 a, n2 a) for
the deformations for each particle in Mγagg ,p by normalized
(n ,n )
by Êp 2 1 .

4. Numerical Computation and Results
Let us show the relevance between δp and γp by the
Monte-Carlo simulations following the algorithm mentioned
in Section 2. Using the algorithm in Section 2, we have ten
pictures of agglomerated particles for each γagg = 0, 0.3, 0.6
and 0.9, and for each p = 0.1, 0.2, 0.3 and 0.4 by letting
L = 2400 and ρ = 10 as in Figures 2 and 3.
On the thickening to compute the CADE, we use two
types thickening process,

In the persistent homology, the Betti number is handled

in general. Since the computational cost to the evaluation
type I:
 →   ,
of the Euler number is not so high and the Euler number

could be compared with the results in [7, 15], we consider
  
the behavior of the Euler numbers of Mt in this article.
type II:
 →   ,
More precisely though there is no guarantee that χ(Mt )
  
is equal to χ(M(i−1) ) for t ∈ [a(i − 1/2), a(i + 1/2)), we
handle χ(M(i−1) ); as mentioned above, in digital analysis,
we should basically neglect finer geometrical difference than such that we generate an octagon asymptotically and apthe pixel size resolution and we follow the principle. In the proximates the area of the disks; In other words, on the
complicated system, we believe that it is quite important
how many topology changes occur for the i-step, and the
difference of the Euler number can represent the behavior.
Table 1. The pattern of thickening
Further the agglomeration can be discriminated whether
the particles are connected or not. From (1), if L/ρ is suffisteps type radius
area n.of pixels
ciently large, even for γagg = 0 and small p, Ê(M0,p , na, 0)
1
II
0.5
0.785
1
does not vanish n > 0, in general, due to the randomness of
2
I
1.5
7.065
9
the configurations. Further the behavior Ê(M; n1 a, n2 a) of
3
I
2.5
19.625
21
n1 , n2 ∈ [0, ρ/a) is quite important since the agglomeration
4
I
3.5
38.465
37
suppresses the topology change in the interval as illustrated
5
II
4.5
63.585
69
in Figures 4 and 5.
6
I
5.5
94.985
97
Due to the randomness of the configurations and the ag7
I
6.5
132.665
129
glomeration, it is not so important whether the Euler num8
II
7.5
176.625
185
bers increase or decrease, but the topological change for the
9
I
8.5
226.865
229
deformation is quite important. We define the agglomera10
I
9.5
283.385
277
(n1 ,n2 )
in (3) more precisely
tion parameter δagg
(4)

(n2 ,n1 )
δagg
(M) =

α
(n2 ,n1 )
Êp(M)

(n ,n )

2
1
− Ê(M; n2 a, n1 a)),
(Eˆp(M)

thickening process, we use the deformation in digital process procedure for each point which is given in Table 1.

6

Further for each point, we consider the thickening:

Since our radius is 10 and our agglomeration algorithm is
characterized by the radius, the behavior of the distribution
0
0
in Figure 5 strongly depends on the regions ρ > 10 and
0
9
9
0
0
9
8
8
9
0
ρ ≤ 10. Figure 4 correspond to (10, 20] region and thus it
0
9
8
8
8
8
9
0
0
9
8
8
7
7
8
8
9
0
implies that our improved thickening algorithm works well
0
9
8
8
7
6
6
7
8
8
9
0
except the first thickening step of the γagg = 0.9 case.
0
9
8
7
6
5
5
6
7
8
9
0
0
9
8
7
6
5
5
6
7
8
9
0
Since the agglomeration in our algorithm means that the
0
9
8
7
6
5
5
6
7
8
9
0 ,
0
9
8
7
6
5
5
6
7
8
9
0
number of agglomerated particles is larger than the uniform
0
9
8
7
6
5
5
6
7
8
9
0
randomness γagg = 0. The variation of the Euler number
0
9
8
8
7
6
6
7
8
8
9
0
0
9
8
8
7
7
8
8
9
0
is related to the deformation in which disjoint clusters con0
9
8
8
8
8
9
0
0
9
8
8
9
0
nect due to the thickening. Agglomeration means that the
0
9
9
0
0
0
number of the disjoint clusters is less than that of uniform
randomness. The variation of the Euler number for the inWe computed the ten pictures for each p = 0.1, 0.2, 0.3, 0.4 creasing of the radius ρ > 10 in Figure 5 is suppressed for
and γagg = 0.0, 0.3, 0.6, 0.9 with ten random seeds. Figure 4 large γagg . Hence the dependence in Figure 4 is very natural
shows the CADE, Ê(Mp,γagg ; na, 0), for the n-th thickening except the first thickening step of the γagg = 0.9 case.
Further in the image processing procedure, we must pay
step for each p.
attention to the digitalized errors. We should recognize that
the first step contains some digitalized errors because the
1000
4000
behavior in Figure 4 is contradict with that in Figure 5; the
γ =0.1
γ =0.1
900
3500
γ =0.3
γ =0.3
800
γ =0.6
γ =0.6
behavior of the curves of γagg = 0.9 in Figure 5 are very
3000
γ
=0.9
γ
=0.9
700
2500
600
mild
over (10, 20) whereas the first steps of γagg = 0.9 in
500
2000
400
1500
Figure
4 rapidly increase.
300
1000
200
Hence
we are concerned with Ê(Mp,γagg ; ρ, a).
500
100
0
9
8
8
7
6
5
5
4
4
4
5
5
6
7
8
8
9
0

0
9
8
7
6
5
5
4
3
3
3
4
5
5
6
7
8
9
0

0
9
8
7
6
5
4
3
2
2
2
3
4
5
6
7
8
9
0

0
9
8
7
6
5
4
3
2
1
2
3
4
5
6
7
8
9
0

0
9
8
7
6
5
4
3
2
2
2
3
4
5
6
7
8
9
0

0
9
8
7
6
5
5
4
3
3
3
4
5
5
6
7
8
9
0

0
9
8
8
7
6
5
5
4
4
4
5
5
6
7
8
8
9
0

CADE

agg
agg
agg
agg

CADE

agg
agg
agg
agg

0
2

4
6
8
thickening step, n

10

0

2

(a)

4
6
8
thickening step, n

10

(b)

8000

γagg=0.1
γagg=0.3
γagg=0.6
γagg=0.9

3500
3000
CADE

2000
CADE

4000

γagg=0.1
γagg=0.3
γagg=0.6
γagg=0.9

2500

1500
1000

euler number

3000

2500

2000

2

4
6
8
thickening step, n

(c)

10

2000
0

-2000

-2000

2000
1500

0

5

10

15
20
radius

25

30

0

5

(a)

0
0

4000

0

500

0

γagg=0.0
γagg=0.3
γagg=0.6
γagg=0.9

6000

4000

1000
500

8000

γagg=0.0
γagg=0.3
γagg=0.6
γagg=0.9

6000

euler number

0

0

2

4
6
8
thickening step, n

(d)

8000

15
20
radius

25

30

25

30

(b)

euler number

8000

γagg=0.0
γagg=0.3
γagg=0.6
γagg=0.9

6000

Figure 4. The CADE Ê(Mp,γagg ; na, 0)
of the n-th thickening step for each
γagg ; Those of the volume fraction p =
0.1, 0.2, 0.3 and 0.4 are illustrated in (a),
(b), (c), and (d) respectively.

10

10

4000
2000

4000
2000

0

0

-2000

-2000
0

5

10

15
20
radius

(c)

γagg=0.0
γagg=0.3
γagg=0.6
γagg=0.9

6000
euler number

0

25

30

0

5

10

15
20
radius

(d)

On the other hand, though it is difficult to identify the
Figure 5. The Euler number vs the racenter points of the particles for given pictures, especially
dius ρ as the point pattern of Mγagg ,p,iS :
of the agglomerated case as shown in images (b) and (c) of
(a), (b), (c), and (d) illustrate the EuFigures 2 and 3, we know the data of the center points of the
ler numbers of the volume fraction p =
particles. Thus we can use the techniques of the statistical
0.1, 0.2, 0.3 and 0.4 respectively.
analysis for the spatial point patterns. Figure 5 displays the
Table 2 and Figure 6 show the dependence of the CADE,
global distribution of the Euler numbers of different radius
2
Ê(M
p,γagg ; ρ, a), on the agglomeration parameter γagg for
of a seed by using the software provided in [7, p.204] .
each p. They exhibit the negative correlations. The agglom2 http://www.maths.jyu.fi/ penttine/ppstatistics.
eration means that the number of agglomerated particles is

7

2500

p=0.1

p=0.2

2000
CADE

CADE

Table 2. CADE vs γagg :
900
800
700
600
500
400
300
200
100
0

1500
1000
500
0

0

0.2

0.4
0.6
γagg

0.8

1

0

0.2

(a)
3000

0.8

1

(b)
3500

p=0.3

p=0.4

3000

2500

2500

2000

CADE

CADE

0.4
0.6
γagg

1500
1000

2000
1500
1000

500
0.2

0.4
0.6
γagg

0.8

1

0

0.2

(c)

0.4
0.6
γagg

0.8

p=0.1

Êp ≡ Êp(ρ/a,1) :=

0.6

δagg

0.6
0.4
0.2

Ave
2705.3
1966.3
1066.6
582.1

Min
2043
1388
766
160
Min
2555
1798
905
464

(ρ/a,1)

0.4

Table 3. δagg vs γagg :

0
0.2

0.4
0.6
γagg

0.8

1

0

0.2

(a)

0.4
0.6
γagg

0.8

1

(b)

p=0.3

p=0.4

1
0.8

0.6

0.6

δagg

0.8

0.4
0.2

0.4
0.2

0

0
0.2

0.4
0.6
γagg

(c)

0.8

1

10
1 X
Ê(M0,p,iS ; ρ, a).
10 i =1

Thus we denote the agglomeration index δagg
by δagg as
in (3). Further in order that δagg corresponds to γagg , we
chose α = 1.2.
Figure 7 and Table 3 show the relation between δagg and
γagg ; both show that γagg is correlated to δagg and approximately recovers δagg up to the statistical fluctuation.

0.2

0

0

Min
2806
1910
1003
248

0.2
Max
2221
1675
886
364
0.4
Max
3072
2160
1188
702

S

p=0.2

1
0.8

1

Ave
2131.9
1501.2
821.9
244.0

Ê(Mp,0,iS ; ρ, a), of the uniform random configuration, i.e.,
γagg = 0 case,

(d)

0.8

0

Min
750
514
269
78

1

Figure 6. The CADE, Ê(Mp,γagg ; ρ, a),
vs γagg : (a), (b), (c), and (d) display
the states of the volume fraction p =
0.1, 0.2, 0.3 and 0.4 respectively.

δagg

Ave
2878.5
2035.2
1121.9
408.8

0.1
Max
879
657
340
204
0.3
Max
2954
2181
1262
507

0
0

δagg

Ave
805.7
582.9
317.6
141.9

500

0

1

p
γagg
0
0.3
0.6
0.9
p
γagg
0
0.3
0.6
0.9

0

0.2

0.4
0.6
γagg

0.8

(d)

Figure 7.
The agglomeration index
δagg vs γagg : (a), (b), (c), and (d) display the states of the volume fraction p
= 0.1, 0.2, 0.3 and 0.4 respectively.

1

p
γagg
0
0.3
0.6
0.9
p
γagg
0
0.3
0.6
0.9

Ave
0.000
0.332
0.727
0.989
Ave
0.000
0.352
0.732
1.030

0.1
Max
0.083
0.434
0.799
1.084
0.3
Max
0.030
0.404
0.782
1.097

Min
-0.109
0.221
0.694
0.896

Ave
0.000
0.355
0.737
1.063

Min
-0.031
0.291
0.674
0.989

Ave
0.000
0.328
0.727
0.942

0.2
Max
0.050
0.419
0.769
1.110
0.4
Max
0.067
0.402
0.799
0.994

Min
-0.050
0.257
0.701
0.995
Min
-0.163
0.242
0.673
0.889

In the statistical analysis of the spatial point patterns,
the Clark-Evans index is a well-established index which represents the agglomeration degree of a given point pattern,
though in general, it is very difficult to identify the center
larger than the uniform randomness γagg = 0 as mentioned points of the particles for a given picture, such as images
above. Since the agglomeration prevents the topological (b) and (c) of Figures 2 and 3; the problem is sometimes
changes on the thickening, we have the negative correlation ill-posed for the cases. Since we know the data of the cenin Figure 6.
ter points of the particles of every Mγagg ,p,iS , we illustrated
We, now, define the the average Êp of a “standard pat- the Clark-Evans index in Table 4 and Figure 8, which show
tern of volume fraction p” by the average of the CADE, that the Clark-Evans index represents our agglomeration

parameter γagg well. The correlation between the ClarkEvans index and δagg is displayed in Figure 9. It shows a
good negative-correlation for each volume fraction p.

1
clark-evans index

clark-evans index

1
0.9
0.8
0.7
0.6
0.5
p=0.1

0.4
0

0.9
0.8
0.7
0.6
0.5
p=0.2

0.4

0.2 0.4 0.6 0.8
γagg

1

0

(a)

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
-0.4

p=0.1
p=0.2
p=0.3
p=0.4

0

0.4 0.8
δagg

1.2

clark-evans index

1

0.9
0.8
0.7
0.6
0.5
p=0.3

0.4
0

0.9

Figure 9. The Clark-Evans index and δagg .

0.8
0.7
0.6
0.5
p=0.4

0.4

0.2 0.4 0.6 0.8
γagg

1

0

(c)

0.2 0.4 0.6 0.8
γagg

(d)

Figure 8.
The Clark-Evans index
vs γagg : (a), (b), (c), and (d) display
the states of the volume fraction p =
0.1, 0.2, 0.3 and 0.4 respectively.

Table 4. Clark Evans index vs γagg :
p
γagg
0
0.3
0.6
0.9
p
γagg
0
0.3
0.6
0.9

1

(b)

1
clark-evans index

0.2 0.4 0.6 0.8
γagg

clark-evans index

8

Ave
1.018
0.755
0.574
0.412
Ave
1.006
0.890
0.788
0.661

0.1
Max
1.030
0.776
0.599
0.419
0.3
Max
1.010
0.897
0.795
0.665

Min
0.999
0.733
0.561
0.402

Ave
1.011
0.837
0.703
0.556

Min
1.002
0.881
0.783
0.655

Ave
1.003
0.927
0.850
0.741

0.2
Max
1.022
0.850
0.709
0.569
0.4
Max
1.008
0.933
0.856
0.746

Min
0.995
0.828
0.693
0.547
Min
1.000
0.922
0.845
0.735

1

algorithm in [11], we constructed digital pictures of aggregated particles controlled by the agglomeration parameter
γagg ∈ (0, 1) as a Monte-Carlo simulation. By applying
the image processing procedure to the pictures, we showed
that δagg statistically reproduces γagg . Since we have the
data of the center points of the particles, we also computed
the well-established Clark-Evans index and showed that it
also represents γagg well. However though the methods in
the point process analysis including the Clark-Evans index
require the data of the configuration of the points, the determination of the center points of the particles for a given
picture is basically an ill-posed problem. Hence our method
has an advantage because we do not need to find the center
points in the computation of δagg . In other words, our purpose that we recover γagg for a given picture by means of
the digital image processing procedure is accomplished by
considering the deformation of the geometrical object. It
implies that we can measure the agglomeration in a given
picture of agglomerated particles.

In this article, we have investigated pictures whose volume fraction p is less than 0.5, because it is difficult to deal
with pictures with large volume fraction. It is expected to
find further natural index to discriminate the agglomeration with the large volume fraction, e.g., in terms of the
persistent homology [8].
Further in [5], T. Kaczynski, K. Mischaikow and M. Mrozek
5. Summary
studied the pattern analysis in a regular lattice using the cuIn this article, we introduced the novel geometrical index bical homology. They also investigated a topological propδagg , which is associated with the Euler number and is ob- erty of the time development of a complicated pattern govtained by an image processing procedure for a given digital erned by the Cahn-Hilliard equation by considering its time
picture of aggregated particles such that δagg represents the development of its Betti numbers [5]. It means that a topodegree of the agglomerations of the particles. Following the logical property of (geometrical or physical) deformation of

9

a complicated geometrical object is important in order to
describe the degree of its complication.

Acknowledgment
The authors are grateful to Professors Y. Fukumoto and
Y. Hiraoka for their critical and helpful comments, and
thank the referee for critical comments and for directing
their attention to the reference [5].
References
[1] R. Bott and L. W. Tu, Differential Forms in Algebraic Topology,
(GTM 82) Springer, Berlin, 1982.
[2] P. J. Clark and C. Evans, Distance to Nearest Neighbor as a
Measure of Spatial Relationships in Populations, Ecology, 35
(1954) 445-453.
[3] H. Edelsbrunner and J. Harer, Persistent homology: A survey, in
Surveys on Discrete and Computational Geometry. Twenty Years
Later, 257282 (J. E. Goodman, J. Pach, and R. Pollack, eds.),
Contemporary Mathematics 453, Amer. Math. Soc., Providence,
Rhode Island, 2008.
[4] L. Hui, R. C. Smith, X. Wang, J. K. Nelson and L. S. Schadler,
Quantification of Particulate Mixing in Nanocomposites, 2008
Annual Report Conference on Electrical Insulation Dielectric Phenomena, IEYASU, (2008) 317-320.
[5] T. Kaczynski, K. Mischaikow, M. Mrozek, Computational Homology (Applied Mathematical Sciences), Springer, New York,
2004.
[6] B. Kelvin, Electrical Units of Measurement, Popular Lectures
and Addresses Volume I, London: Macmillan and Co., 1889, pp.
73-74.
[7] J. Illian, A. Penttinen, H. Stoyan, D. Stoyan, Statistical Analysis
and Modelling of Spatial Point Patterns (Statistics in Practice),
Wiley, New York, 2008.
[8] R. MacPherson and B. Schweinhart, Measuring shape with topology, J. Math. Phys., 53 (2012) 073516(13 pages).
[9] K. R. Mecke and D. Stoyan, Morphological Characterization of
Point Patterns, Biometrical J., 47 (2005) 473-488.
[10] S. Matsutani, Y. Shimosako, and Y. Wang, Numerical Computations of Conductivity in Continuum Percolation for Overlapping
Spheroids, Int. J. Mod. Phys. C, 21 (2010) 709-729.
[11] S. Matsutani, Y. Shimosako, and Y. Wang, Numerical Computations of Conductivities over Agglomerated Continuum Percolation Models, Appl. Math. Modeling, 37 (2013) 4007-4022.
[12] R. E. Miles, Estimating aggregate and overall characteristics
from thick sections by transmission microscopy, J. of Microscopy,
107 (1976) 227-729.
[13] W. K. Pratt, Digital Image Processing, 2nd ed., Wiley, New
York, 1991.
[14] J. A. Sethian, Level Set Methods and Fast Marching Methods:
Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge Univ.
Press Cambridge, 1999.
[15] D. Stoyan, W. S. Kendall and J. Mecke, Stochastic Geometry
and its Applications, 2nd ed., Wiley, New York, 1995.
[16] A. Tscheschel and D. Stoyan, Statistical reconstruction of random point patterns, Comp. Stat. Data Anal., 51 (2006) 859-871.
[17] S. Weinberger, What is Persistent Homology?, Notices of the
AMS, 58 (2011) 36-39.

Shigeki Matsutani, Yoshiyuki Shimosako
Simulation & Analysis R&D Center,
Canon Inc., 3-30-2, Shimomaruko Ohta-ku,
Tokyo, Japan

