- Email: [email protected]
Vesicle computers: Approximating a Voronoi diagram using Voronoi automata
Chaos, Solitons & Fractals 44 (2011) 480–489
Contents lists available at ScienceDirect
Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Vesicle computers: Approximating a Voronoi diagram using Voronoi automata Andrew Adamatzky a,⇑, Ben De Lacy Costello a,b, Julian Holley a,d, Jerzy Gorecki c, Larry Bull d a
Unconventional Computing Centre, University of the West of England, Bristol, United Kingdom Centre for Analytical, Material and Sensor Sciences, University of the West of England, Bristol, United Kingdom c Institute of Physical Chemistry PAN, Warsaw, Poland d Artificial Intelligence Group, University of the West of England, Bristol, United Kingdom b
a r t i c l e
i n f o
Article history: Received 15 April 2010 Accepted 27 January 2011
a b s t r a c t Irregular arrangements of vesicles filled with excitable and precipitating chemical systems are imitated by Voronoi automata – finite-state machines defined on a planar Voronoi diagram. Every Voronoi cell takes four states: resting, excited, refractory and precipitate. A resting cell excites if it has at least one neighbour in an excited state. The cell precipitates if the ratio of excited cells in its neighbourhood versus the number of neighbours exceeds a certain threshold. To approximate a Voronoi diagram on Voronoi automata we project a planar set onto the automaton lattice, thus cells corresponding to data-points are excited. Excitation waves propagate across the Voronoi automaton, interact with each other and form precipitate at the points of interaction. The configuration of the precipitate represents the edges of an approximated Voronoi diagram. We discover the relationship between the quality of the Voronoi diagram approximation and the precipitation threshold, and demonstrate the feasibility of our model in approximating Voronoi diagrams of arbitrary-shaped objects and in constructing a skeleton of a planar shape. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Manually designed regular or irregular arrangements of vesicles filled with excitable chemical media exhibit huge computational potential [6]. When vesicles are in close contact with each other (at least in terms of diffusion) joined via tiny pores, excitation waves can pass from one vesicle to its nearest neighbours. If level of excitation is controlled and the size of the vesicle is appropriate then excitation wave-fragments are able to maintain a constant shape as they travel within the BZ-vesicle. A wave-fragment passing from one vesicle to another is modulated by the restricted size of the connecting pore causing the wave fragment to contract. When two or more wave-fragments collide inside a vesicle they can annihilate (either as ⇑ Corresponding author. E-mail address: [email protected] (A. Adamatzky). 0960-0779/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2011.01.016
a direct result of the collision or via deflection into the vesicle wall at a position without a connecting pore), deviate, or multiply. When interpreting the presence or absence of wave-fragments in certain vesicles of a vesicle network as a values of a Boolean variable then we can implement all basic Boolean logical operations via the collisions between wave-fragments. In computer experiments we designed a binary adder in a hexagonal array of vesicles filled with excitable chemical media [6], built polymorphic logical gates (switching between XNOR and NOR) by using illumination to control outcomes of inter-fragment collisions [7] and geometry-modulated complex arithmetical circuits [25]. Our theoretical ideas and results are backed up by experimental results that exhibit that excitation transfer occurs between vesicles in a Belousov–Zhabotinsky mixture enclosed by a lipid membrane [23,32]. While in computer models constructing regular arrangements of uniform vesicles is effortless real-life
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
experiments prove much more difficult. Usually vesicles are different sizes, they do not form a hexagonal lattice as a rule, and they may be unstable. Coalescence transforms fine-grained networks of elementary vesicle-processors into a coarse-grained assembly of very large vesicle structures. What kind of computation can be undertaken utilising an irregular arrangement of non-uniform vesicles? We address the question by representing the vesicle assembles by automata networks and studying how the problem of planar subdivision, construction of a Voronoi diagram, can be solved using these vesicle-automata. We abstract vesicle assembles as planar Voronoi diagrams (Fig. 1) of planar sets, points of which are centres of the vesicles. Voronoi diagrams are routinely used as an approximation of arrangements of discs [22] and sphere packing [30,20,31]. The diagram is also used in structural analysis of liquids and gases [9], and protein structure [36], and to model dense gels [41] and inter-atomic bonds [24]. Voronoi diagrams are introduced in Section 2. We assume that every cell of a Voronoi diagram is a finite-state
481
automata that takes four states and updates its states depending on the states of its first and second order neighbours. We design a cell-state transition function which combines generalised, and highly-abstracted, properties of both excitable and precipitating chemical media: a local disturbance gives birth to quasi-circular waves of excitation whilst collisions between the waves lead to precipitation. The Voronoi automaton is defined in Section 3. The problem solved by Voronoi automata is the approximation of a Voronoi diagram. Approximated Voronoi diagrams are much more coarse-grained than the Voronoi diagram on which the excitable-precipitating automaton is built. To approximate a Voronoi diagram on a Voronoi automata we project a planar set onto the automaton lattice, thus cells corresponding to data-points are excited. Excitation waves propagate across the Voronoi automaton, interact with each other and form precipitate as a result of their interaction. Configuration of precipitate represents edges of approximated Voronoi diagram. In our model precipitation depends on a local density of excitation — the
Fig. 1. Examples of Voronoi representations in vesicle conglomerates: (a) rhodamin, vegetable oil, water and silicon oil (polydimethylsiloxane) mixtures, (b) blue food colouring and vegetable oil mixtures (c) rhodamin, water and silicon oil mixtures.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
precipitation threshold. For a low precipitation threshold the medium becomes cluttered with meaningless clusters of precipitate, for a high threshold few domains of precipitation are formed. Our quest for an optimal threshold of precipitation is described in Section 4. In Section 5 we show that the optimal threshold found can be used equally successfully to approximate Voronoi diagrams of arbitrary planar shapes and also in constructing skeletons of planar shapes. 2. Voronoi diagram Let P be a nonempty finite set of planar points. A planar Voronoi diagram [40] of the set P is a partition of the plane into such regions, that for any element of P, a region corresponding to a unique point p contains all those points of the plane which are closer to p than to any other node of P. A unique region vor(p) = {z 2 R2 : d(p, z) < d(p, m), m 2 R2, m – z} assigned to point p is called a Voronoi cell of the point p [37]. The boundary @vor(p) of the Voronoi cell of a point p is built of segments of bisectors separating pairs of geographically closest points of the given planar set P. A union of VD(P) = [p2P@vor(p) all boundaries of the Voronoi cells determines the planar Voronoi diagram [37]. A variety of Voronoi diagrams and algorithms of their construction can be found in [27,33]. Approximation of Voronoi diagrams with propagating patterns is based on time-to-distance transformation: to approximate a bisector separating planar points p and q we initiate growing patterns at p and q. The patterns travel the same distance from the sites of origination before they meet each other, The loci where the waves meet indicate sites of the computed bisector [1,2]. Precipitating reaction–diffusion chemical media have proved to be an ideal computing substrate for approximation of the planar Voronoi diagram [39,18,3]. A Voronoi diagram can be approximated in a two-reagent medium. One reagent a is dissolved in a suitable substrate (usually a gel), drops of another reagent b are applied to the sites corresponding to the planar points to be separated by bisectors. The reagent b diffuses through the substrate and reacts with reagent a. A coloured precipitate is produced in the reaction between a and b. When two or more waves meet, no precipitate is formed [18]. Thus uncoloured loci of the reaction–diffusion media represent the bisectors of the computed diagram. A range of chemical precipitating processors of this general type have been designed and working prototypes have been tested in laboratory conditions [39,16–19]. 3. Voronoi automata A Voronoi automaton is a tuple V ¼ hVðPÞ; Q ; N; u; f i, where B is a finite planar set, V(P) = {V(p) : p 2 P}, Q is finite set, N is a set of natural numbers and u : V(P) ? V(P)k is second-order neighbourhood, 0 < k < jPj, and f : Qk ? Q is a cell-state transition function. Excitable-precipitating Voronoi automata studied in the present paper are specified as follows. The cell state set has four elements, Q = {, +, , #}. Thus we assign three
excitation-related states – resting (), excited (+) and refractory () – to the cells, and one precipitate state #. Cells update their state in discrete time. A state of cell V(x) at time step t 2 N is denoted as V(x)t. All cells update their states in parallel using the same cell-state transition function. Let w(p), p 2 P, be a first-order neighbourhood of a Voronoi cell Vor(p) from V(P) : w(p) = {V(q) 2 V(P) : @V(p) \ @V(q)}, i.e. a set of Voronoi cells which have common edges with V(x). A second-order neighbourhood u(p) is a set of neighbours of first-order neighbours of T V(p) : u(p) = V(q)2w(V(p))w(V(q)). Examples of Voronoi cell
0.5 0.45 0 45 0.4 0.35 Ratioo of cells
482
0.3 0.25 0.2 0.15 0.1 0.05 0 1
3
5
7
9 11 13 15 17 19 21 23 25 27 29 31 33 35 Number of neighbours of Voronoi cell
Fig. 2. Structure and properties of neighbourhoods. (a) Example of Voronoi tessellation with second-order neighbourhood structures highlighted for four cells. In each case central cell p is filled with blue (black in grey-scale reproduction) colour, first-order neighbours of w(p) are filled with red (dark-grey) colour, and neighbours of w(p) are filled with green (light-grey) colour. For each blue (black) Voronoi cell a second-order neighbourhood is a set of red (dark-grey) and green (light grey) Voronoi cells. (b) Distribution of first order, or immediate, neighbours (dotted line), second order neighbours (dashed line) and size of second order neighbourhood (solid line) in Voronoi automata of 15 K cells. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
483
Fig. 3. Time lapsed contours of activation in Voronoi automaton with density / = 0.407 and activation threshold (a) g = 0.2 and (b) g = 0.3. Propagating wave front is converted to a contour every 10th step of simulation. Initially eastmost nodes are activated.
Fig. 4. Approximation of Voronoi diagram VðBÞ on Voronoi automaton V; g ¼ 0:4. Points of B are projected onto V at time step t = 0. Seven Voronoi cells of V are excited and generate quasi-circular excitation waves. Configuration of V at time step t = 1, when waves have just started to develop, is shown in (a). The waves propagate outwards from their initial stimulation sites and convert the Voronoi cells they are occupying into refractory states (bc). When two or more waves collide precipitation occurs. By the 44th step of the automaton development of the original excitation is extinct but domains of precipitate-cells represent edges of the approximated Voronoi diagram V. Resting cells are blank, excited cells red (dark grey), refractory cells grey, and precipitate cells black. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
484
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
neighbourhoods are shown in Fig. 2a. Distributions of numbers of immediate and second neighbours and sizes of second-order neighbourhood u are given in Fig. 2b. Predominant number, over 40%, of cells have six first-order (immediate) neighbours. Almost half of the cells have 13 or 14 second-order neighbours. For any cell size of second-order neighbourhood u is a sum of first-order neighbours and second-order neighbours. For Voronoi automata with high-density packing of 15 K cells, half of the cells have either 19, or 20 or 21 neighbours. If we consider only immediate neighbourhood then V can be, in principle, seen as a slightly distorted, hexagonal lattice. However, by adopting second-order neighbourhood we
are moving to a random structure with higher than hexagonal lattice coordinate number (a node in hexagonal lattice has 18-node second-order neighbourhood while dominating neighbourhood sizes in V are 19, 20 and 21). Transition from the excited to the refractory state is unconditional, i.e. takes place regardless of the states of the cell’s neighbours. A resting cell excites if it has at least one excited neighbour. In this particular model we take a refractory state and precipitate state as absorbing: once a cell takes either of these two states it does not update its state any longer. Neighbourhood sizes may differ between Voronoi cells therefore we use a state transition function, where a cell
Fig. 5. Quality of approximation depends on g. (a) Voronoi diagram VðBÞ of seven planar points set B constructed by Fortuna’s sweep line algorithm, (b)–(o) stationary configurations of precipitate in automaton V excited by B for different g = 0.2. . .0.525. In snapshots (b)–(o) black pixels symbolise Voronoi cells of V in precipitate state, cells in refractory states are blank.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
updates its state depending on a relative excitation in its neighbourhood. We assume that precipitation occurs in a resting cell when a ratio of excited neighbours to a total number of neighbours exceeds some threshold g 2 [0, 1]. Let r(V(x)t) be a number of excited cells in the cell V(x)’s second-order neighbourhood, rðVðxÞt Þ ¼ P t VðyÞ2uðVðxÞÞ jfVðyÞ : VðyÞ ¼ þgj then a cell updates its state by the following rule:
VðxÞtþ1
8 #; if VðxÞt ¼ and rt ðxÞ=mðxÞ > g > > > < þ; if VðxÞt ¼ and rt ðxÞ=mðxÞ > 1 ¼ > ; if VðxÞt ¼ þ > > : ; otherwise
Example of propagating wave-fronts are shown in Fig. 3. For g = 0.3 excitation wave-front is convex, shaping into planar by the end of propagation (Fig. 3b). In case of lower threshold, g = 1, wave near the edges of moves quicker then inside the diagram core (Fig. 3a). Thus wave front becomes concave. The part of wave front travelling along edges of diagram reaches the side opposite to the initial perturbation side quicker then the front propagating inside the diagram. Thus wave front becomes closed (Fig. 3a). The domain surrounded by a nested group of target waves, in the western part of the diagram in Fig. 3a, is the place where the excitation waves travelling from east collapse. Such patterns of excitation are somewhat universal and can be observed in various systems, even as calcium waves in coupled cells induced by internal noise [35]. In computational experiments we construct set P by filling a disc-container of radius 480 units with up to 15 K points. The points are packed at random but there is
485
at least 5 units distance between any two points. Voronoi diagram V – on which Voronoi automaton is constructed – is calculated by a classical sweep line algorithm [21].
4. Constructing Voronoi diagram on Voronoi automata Let B be a set of planar points on which a Voronoi diagram V on V; g ¼ 0:4, is approximated. We project B onto V and excite cells of V which are closer than 9 units to points of B. Excitation waves spread on V. The excitation waves collide and precipitation occurs near the sites of the waves’ collisions. Configuration of cells in precipitate state represents edges of VðBÞ. An example of excitationprecipitation dynamics in V which approximates VðBÞ, where B is a planar set of seven points is shown in Fig. 4. A scheme of the Voronoi diagram computed is shown in Fig. 8a. Even in a single example shown in Fig. 4 it is clear that V does not calculate VðBÞ precisely, it rather approximates it with discrete domains of precipitation, and also introduces some noise (sites of precipitation not coinciding with bisectors of VðVÞ). Only threshold of precipitation g can be varied in our model. How does the degree of approximation and signal to noise ratio depend on the value of g? To answer we undertook a series of computational experiments, illustrated – for some values of g – in Fig. 5. For a set B of seven points we constructed ‘classical’ planar Voronoi diagrams VðBÞ (Fig. 5a). Then for g = 0.2. . .0.525 we excited V with B, waited until V reaches its stationary configuration, where all cells are in either precipitate or refractory state, and compared the configurations of V and the image of VðBÞ. For each g we recorded
Fig. 6. Characterisation of approximation of VðBÞ (Fig. 5a) by automaton V: (a) q vs g and (b) m vs g. Markers (rhombs and circles) show actual values of m(g) and q(g), while solid lines represent a moving average trend line purely for guidance.
486
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
a fraction q of precipitate-cells in a final configuration of VðgÞ to a total number of precipitate-cells in the same configuration (Fig. 6a).
We found that while density m, calculated as a ratio of the number of precipitate-cells in a stationary configuration VðgÞ versus the number of precipitate-cells in
Fig. 7. Approximation of Voronoi diagram VðOÞ of planar shapes O in Voronoi automaton V. Voronoi cell of V in precipitate state are black, cells in refractory states are blank. Excited cells are red (grey). (a) representation of O, (b)–(f) development of V. VðOÞ is represented by black pixels in (f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
487
Fig. 8. Schematic representation of data objects and their Voronoi diagram or skeleton for precipitation-based approximation shown in (a) Fig. 4, (b) Fig. 7, and (c) Fig. 9.
configuration of Vð0:2Þ, exponentially decreases with increase of precipitation threshold g (Fig. 6b) degree of approximation q shows more intriguing behaviour (Fig. 6a). From g = 0.2 till g = 0.41 ‘signal-to-noise’ ratio q(g) grows almost exponentially. Trend of q(g) undergoes S-type transition between g = 0.41 and g = 0.453 with one peak q(0.44) and two cavities q(0.43) and q(0.452). As illustrated in Fig. 5 V with precipitation threshold in the zone around g = 0.4, onset of ‘strange’ behaviour of q(g), demonstrates best approximation of VðBÞ, as recognised by eye.
5. Arbitrary-shaped planar objects and contours The Voronoi automaton V copes well with data set O of planar finite shapes. To approximate diagram V of arbitrary planar shapes we project the objects of O (Fig. 7a) onto V in such a manner that at time t = 0 all cells of V covered (coinciding with pixels of) by shapes from O becomes excited (Fig. 7b). Waves of excitation travel across V and provoke precipitation during their collisions (Fig. 7c–e). When all cells take refractory or precipitate state the automaton’s configuration becomes stationary. The spatial distribution of precipitate represents a Voronoi diagram of the data shapes (see Fig. 7f and Fig. 8b). A precipitateenhancement of the original position of the data shape is a byproduct of V’s development. A skeleton of a planar contour is a set of centres of bitangent circles lying inside the contour [15]. Blum’s grass-fire algorithm for computing skeleton employs propagating patterns [12,13,15]: to compute a skeleton we set a contour on fire and let the fire spread, quenched points where the advancing fire-fronts collide represent the skeleton. Most known algorithms of skeletonisation are based on Blum’s approach: simulation of grass-fire [29], distance transform [38], analytical construction of medial axis and topological thinning [34,10]. To approximate a skeleton of a shape S (Fig. 9a) we project S onto V. Voronoi cells of V corresponding to black pixels of S are excited (Fig. 9b). Excitation waves propagate inside and outside of the domain S of initial excitation;
we can ignore behaviour of outward waves. Waves travelling inside the contour S trigger precipitation at the sites of the waves’ interaction (Fig. 9c–e). Distribution of precipitate (Fig. 9e) inside S approximates skeleton of S (Fig. 8c).
6. Summary We demonstrated that is possible to solve some problems of computational geometry on a discrete model of irregular arrangements of non-uniform vesicles filled with excitable chemical mixtures [6,7,25]. We have shown that a threshold of relative local excitation density can be used to parameterise excitation and precipitation dynamics, enhance results of the computation and reduce signalto-noise ratio. We proved that the same optimal threshold of precipitation works well not only for approximating Voronoi diagram of planar sets but also for arbitrary shapes and skeletonisation of planar contours. In terms of automata-network based computation we advanced our previous results on excitation in Delaunay automata [4] and b-skeletons [5,8]. The parameterisation developed could be used in designing experimental laboratory prototypes of finegrained compartmentalised excitable chemical processors, e.g. using micro-emulsion approach [26], and developing nano-scale massively parallel computers [11]. Namely, we have already designed computers models of lightsensitive BZ-vesicles packed into a hexagonal lattice. The vesicles implement a set of universal logical gates and also one-bit binary adder [6]. Our parameterisation approach is used to select levels of illumination, which in turn controls level of the medium’s excitability, and size of pores connecting BZ-vesicles, which determines propagation of wave-fragments. With regards to laboratory implementation in mind we must mention that our approach is not the only one by organically complement early developed techniques of automata-based simulation of reaction–diffusion chemical systems. For example, lattice gas automata is amongst highly used computer simulation tools [28,14]. The approach based on Voronoi automata, discussed in present
488
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
Fig. 9. Approximation of a skeleton of a star-contour S by V. Voronoi cell of V in precipitate state are black, cells in refractory states are blank. Excited cells are red (grey). (a) representation of a star-contour, (b)–(e) development of V. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
paper, is less accurate than that based on lattice gas automata however provides higher computational efficiency
with degree of accuracy high enough for rapid prototyping of massively-parallel chemical processors.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
Acknowledgement The work is part of the European project 248992 funded under 7th FWP (Seventh Framework Programme) FET Proactive 3: Bio-Chemistry-Based Information Technology CHEM-IT (ICT-2009.8.3). References [1] Adamatzky A. Reaction–diffusion algorithm for constructing discrete generalised Voronoi diagram. Neural Netw World 1994;9:6635–43. [2] Adamatzky A. Voronoi-like partition of lattice in cellular automata. Math Comput Model 1996;23:51–66. [3] Adamatzky A, De Lacy Costello B, Asai T. Reaction–diffusion computers. Elsevier; 2005. [4] Adamatzky A. On excitable Delaunay automata (Kybenertes) 2011, in press.. [5] Adamatzky A. On excitable beta-skeletons. J Comput Sci 2010;1:175–86. [6] Adamatzky A, Holley J, Bull L, De Lacy Costello B. On computing in finegrained compartmentalised Belousov–Zhabotinsky medium. Chaos Solitons Fract 2010. Available from: . [7] Adamatzky A, De Lacy Costello B, Bull L. On polymorphic logical gates in sub-excitable chemical medium. Int J Bifurcat Chaos; in press. . [8] Alonso-Sanz R, Adamatzky A. On beta-skeleton automata with memory. J Comput Sci. 2011;2(1):57–66. [9] Anikeenko AV, Alinchenko MG, Voloshin VP, Medvedev NN, Gavrilova ML, Jedlovszky P. Implementation of the Voronoi– Delaunay method for analysis of intermolecular voids. Lect Notes Comput Sci 2004;3045:217–26. [10] Attali D, Montanvert A. Computing and simplifying 2D and 3D continuous skeletons. Comput Vis Image Understand 1997;67:261–73. [11] Bandyopadhyay A, Pati R, Sahu S, Peper F, Fujita D. Massively parallel computing on an organic molecular layer. Nat Phys 2010;6. [12] Blum H. A transformation for extracting new descriptors of shape. In: Wathen-Dunn W, editor. Models for the perception of speech and visual form. MIT Press; 1967. [13] Blum H. Biological shape and visual science. J Theor Biol 1973;38:205–87. [14] Boon J-P, Eijnden EV, Hanon D. A lattice gas automaton approach to ‘‘turbulent diffusion’’. Chaos Solitons Fract 2000;11:187–92. [15] Calabi L, Hartnett WE. Shape recognition, prairie fires, convex deficiencies and skeletons. Am Math Mon 1968;75:335–42. [16] De Lacy Costello BPJ. Constructive chemical processors – Experimental evidence that shows this class of programmable pattern forming reactions exist at the edge of a highly nonlinear region. Int J Bifurcat Chaos 2003;13:1561–4. [17] De Lacy Costello BPJ, Adamatzky A. On multitasking in parallel chemical processors: experimental results. Int J Bifurcat Chaos 2003;13:521–33. [18] De Lacy Costello B, Hantz P, Ratcliffe N. Voronoi diagrams generated by regressing edges of precipitation fronts. J Chem Phys 2004;120:2413.
489
[19] De Lacy Costello BPJ, Jahan I, Adamatzky A, Ratcliffe NM. Chemical tessellations. Int J Bifurcat Chaos 2009;19:619–22. [20] Filatovs GJ. Delaunay-subgraph statistics of disc packings materials characterisation 1998;40(1):27–35. [21] Fortune S. A sweep line algorithm for Voronoi diagrams. In: Proceedings of the 2nd annual symposium on computational geometry; 1986. p. 313–22. [22] Gervois A, Annic C, Lemaitre J, Ammi m, Oger L, Troadec J-P. Arrangement of discs in 2D binary assemblies. Physica A 1995;218:403–18. [23] Gorecki J. Private communication; 2010. [24] Hobbs WL. Network topology in aperiodic networks. J Non-Crystal Solids 1995;192/193:79–91. [25] Holley J, Adamatzky A, Bull L, De Lacy Costello B, Jahan I. Modalities of Belousov-Zhabotinsky encapsulated vesicles. Nanocommunication Netw J; in press.. [26] Kaminaga A, Vanag VK, Epstein IR. A reaction–diffusion memory device. Angew Chem Int Ed 2006;45:3087–9. [27] Klein R. Concrete and abstract Voronoi diagrams. Berlin: SpringerVerlag; 1990. [28] Lawniczak A, Dab D, Kapral R, Boon J-P. Reactive lattice gas automata. Physica D 1991;47:132–58. [29] Leymarie F, Levine MD. Simulating the grass fire transform using an active contour model. IEEE Trans Pattern Anal Mach Intell 1992;14:56–75. [30] Lochmann K, Oger L, Stoyan D. Statistical analysis of random sphere packings with variable radius distribution. Solid State Sci 2006;8:1397–413. [31] Luchnikov VA, Gavrilova ML, Medvedev NN, Voloshin VP. The Voronoi–Delaunay approach for the free volume analysis of a packing of balls in a cylindrical container. Future Generat Comput Syst 2002;18:673–9. [32] NeuNeu: Artificial wet neuronal networks from compartmentalised excitable chemical media; 2010.
Contents lists available at ScienceDirect
Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Vesicle computers: Approximating a Voronoi diagram using Voronoi automata Andrew Adamatzky a,⇑, Ben De Lacy Costello a,b, Julian Holley a,d, Jerzy Gorecki c, Larry Bull d a
Unconventional Computing Centre, University of the West of England, Bristol, United Kingdom Centre for Analytical, Material and Sensor Sciences, University of the West of England, Bristol, United Kingdom c Institute of Physical Chemistry PAN, Warsaw, Poland d Artificial Intelligence Group, University of the West of England, Bristol, United Kingdom b
a r t i c l e
i n f o
Article history: Received 15 April 2010 Accepted 27 January 2011
a b s t r a c t Irregular arrangements of vesicles filled with excitable and precipitating chemical systems are imitated by Voronoi automata – finite-state machines defined on a planar Voronoi diagram. Every Voronoi cell takes four states: resting, excited, refractory and precipitate. A resting cell excites if it has at least one neighbour in an excited state. The cell precipitates if the ratio of excited cells in its neighbourhood versus the number of neighbours exceeds a certain threshold. To approximate a Voronoi diagram on Voronoi automata we project a planar set onto the automaton lattice, thus cells corresponding to data-points are excited. Excitation waves propagate across the Voronoi automaton, interact with each other and form precipitate at the points of interaction. The configuration of the precipitate represents the edges of an approximated Voronoi diagram. We discover the relationship between the quality of the Voronoi diagram approximation and the precipitation threshold, and demonstrate the feasibility of our model in approximating Voronoi diagrams of arbitrary-shaped objects and in constructing a skeleton of a planar shape. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Manually designed regular or irregular arrangements of vesicles filled with excitable chemical media exhibit huge computational potential [6]. When vesicles are in close contact with each other (at least in terms of diffusion) joined via tiny pores, excitation waves can pass from one vesicle to its nearest neighbours. If level of excitation is controlled and the size of the vesicle is appropriate then excitation wave-fragments are able to maintain a constant shape as they travel within the BZ-vesicle. A wave-fragment passing from one vesicle to another is modulated by the restricted size of the connecting pore causing the wave fragment to contract. When two or more wave-fragments collide inside a vesicle they can annihilate (either as ⇑ Corresponding author. E-mail address: [email protected] (A. Adamatzky). 0960-0779/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2011.01.016
a direct result of the collision or via deflection into the vesicle wall at a position without a connecting pore), deviate, or multiply. When interpreting the presence or absence of wave-fragments in certain vesicles of a vesicle network as a values of a Boolean variable then we can implement all basic Boolean logical operations via the collisions between wave-fragments. In computer experiments we designed a binary adder in a hexagonal array of vesicles filled with excitable chemical media [6], built polymorphic logical gates (switching between XNOR and NOR) by using illumination to control outcomes of inter-fragment collisions [7] and geometry-modulated complex arithmetical circuits [25]. Our theoretical ideas and results are backed up by experimental results that exhibit that excitation transfer occurs between vesicles in a Belousov–Zhabotinsky mixture enclosed by a lipid membrane [23,32]. While in computer models constructing regular arrangements of uniform vesicles is effortless real-life
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
experiments prove much more difficult. Usually vesicles are different sizes, they do not form a hexagonal lattice as a rule, and they may be unstable. Coalescence transforms fine-grained networks of elementary vesicle-processors into a coarse-grained assembly of very large vesicle structures. What kind of computation can be undertaken utilising an irregular arrangement of non-uniform vesicles? We address the question by representing the vesicle assembles by automata networks and studying how the problem of planar subdivision, construction of a Voronoi diagram, can be solved using these vesicle-automata. We abstract vesicle assembles as planar Voronoi diagrams (Fig. 1) of planar sets, points of which are centres of the vesicles. Voronoi diagrams are routinely used as an approximation of arrangements of discs [22] and sphere packing [30,20,31]. The diagram is also used in structural analysis of liquids and gases [9], and protein structure [36], and to model dense gels [41] and inter-atomic bonds [24]. Voronoi diagrams are introduced in Section 2. We assume that every cell of a Voronoi diagram is a finite-state
481
automata that takes four states and updates its states depending on the states of its first and second order neighbours. We design a cell-state transition function which combines generalised, and highly-abstracted, properties of both excitable and precipitating chemical media: a local disturbance gives birth to quasi-circular waves of excitation whilst collisions between the waves lead to precipitation. The Voronoi automaton is defined in Section 3. The problem solved by Voronoi automata is the approximation of a Voronoi diagram. Approximated Voronoi diagrams are much more coarse-grained than the Voronoi diagram on which the excitable-precipitating automaton is built. To approximate a Voronoi diagram on a Voronoi automata we project a planar set onto the automaton lattice, thus cells corresponding to data-points are excited. Excitation waves propagate across the Voronoi automaton, interact with each other and form precipitate as a result of their interaction. Configuration of precipitate represents edges of approximated Voronoi diagram. In our model precipitation depends on a local density of excitation — the
Fig. 1. Examples of Voronoi representations in vesicle conglomerates: (a) rhodamin, vegetable oil, water and silicon oil (polydimethylsiloxane) mixtures, (b) blue food colouring and vegetable oil mixtures (c) rhodamin, water and silicon oil mixtures.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
precipitation threshold. For a low precipitation threshold the medium becomes cluttered with meaningless clusters of precipitate, for a high threshold few domains of precipitation are formed. Our quest for an optimal threshold of precipitation is described in Section 4. In Section 5 we show that the optimal threshold found can be used equally successfully to approximate Voronoi diagrams of arbitrary planar shapes and also in constructing skeletons of planar shapes. 2. Voronoi diagram Let P be a nonempty finite set of planar points. A planar Voronoi diagram [40] of the set P is a partition of the plane into such regions, that for any element of P, a region corresponding to a unique point p contains all those points of the plane which are closer to p than to any other node of P. A unique region vor(p) = {z 2 R2 : d(p, z) < d(p, m), m 2 R2, m – z} assigned to point p is called a Voronoi cell of the point p [37]. The boundary @vor(p) of the Voronoi cell of a point p is built of segments of bisectors separating pairs of geographically closest points of the given planar set P. A union of VD(P) = [p2P@vor(p) all boundaries of the Voronoi cells determines the planar Voronoi diagram [37]. A variety of Voronoi diagrams and algorithms of their construction can be found in [27,33]. Approximation of Voronoi diagrams with propagating patterns is based on time-to-distance transformation: to approximate a bisector separating planar points p and q we initiate growing patterns at p and q. The patterns travel the same distance from the sites of origination before they meet each other, The loci where the waves meet indicate sites of the computed bisector [1,2]. Precipitating reaction–diffusion chemical media have proved to be an ideal computing substrate for approximation of the planar Voronoi diagram [39,18,3]. A Voronoi diagram can be approximated in a two-reagent medium. One reagent a is dissolved in a suitable substrate (usually a gel), drops of another reagent b are applied to the sites corresponding to the planar points to be separated by bisectors. The reagent b diffuses through the substrate and reacts with reagent a. A coloured precipitate is produced in the reaction between a and b. When two or more waves meet, no precipitate is formed [18]. Thus uncoloured loci of the reaction–diffusion media represent the bisectors of the computed diagram. A range of chemical precipitating processors of this general type have been designed and working prototypes have been tested in laboratory conditions [39,16–19]. 3. Voronoi automata A Voronoi automaton is a tuple V ¼ hVðPÞ; Q ; N; u; f i, where B is a finite planar set, V(P) = {V(p) : p 2 P}, Q is finite set, N is a set of natural numbers and u : V(P) ? V(P)k is second-order neighbourhood, 0 < k < jPj, and f : Qk ? Q is a cell-state transition function. Excitable-precipitating Voronoi automata studied in the present paper are specified as follows. The cell state set has four elements, Q = {, +, , #}. Thus we assign three
excitation-related states – resting (), excited (+) and refractory () – to the cells, and one precipitate state #. Cells update their state in discrete time. A state of cell V(x) at time step t 2 N is denoted as V(x)t. All cells update their states in parallel using the same cell-state transition function. Let w(p), p 2 P, be a first-order neighbourhood of a Voronoi cell Vor(p) from V(P) : w(p) = {V(q) 2 V(P) : @V(p) \ @V(q)}, i.e. a set of Voronoi cells which have common edges with V(x). A second-order neighbourhood u(p) is a set of neighbours of first-order neighbours of T V(p) : u(p) = V(q)2w(V(p))w(V(q)). Examples of Voronoi cell
0.5 0.45 0 45 0.4 0.35 Ratioo of cells
482
0.3 0.25 0.2 0.15 0.1 0.05 0 1
3
5
7
9 11 13 15 17 19 21 23 25 27 29 31 33 35 Number of neighbours of Voronoi cell
Fig. 2. Structure and properties of neighbourhoods. (a) Example of Voronoi tessellation with second-order neighbourhood structures highlighted for four cells. In each case central cell p is filled with blue (black in grey-scale reproduction) colour, first-order neighbours of w(p) are filled with red (dark-grey) colour, and neighbours of w(p) are filled with green (light-grey) colour. For each blue (black) Voronoi cell a second-order neighbourhood is a set of red (dark-grey) and green (light grey) Voronoi cells. (b) Distribution of first order, or immediate, neighbours (dotted line), second order neighbours (dashed line) and size of second order neighbourhood (solid line) in Voronoi automata of 15 K cells. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
483
Fig. 3. Time lapsed contours of activation in Voronoi automaton with density / = 0.407 and activation threshold (a) g = 0.2 and (b) g = 0.3. Propagating wave front is converted to a contour every 10th step of simulation. Initially eastmost nodes are activated.
Fig. 4. Approximation of Voronoi diagram VðBÞ on Voronoi automaton V; g ¼ 0:4. Points of B are projected onto V at time step t = 0. Seven Voronoi cells of V are excited and generate quasi-circular excitation waves. Configuration of V at time step t = 1, when waves have just started to develop, is shown in (a). The waves propagate outwards from their initial stimulation sites and convert the Voronoi cells they are occupying into refractory states (bc). When two or more waves collide precipitation occurs. By the 44th step of the automaton development of the original excitation is extinct but domains of precipitate-cells represent edges of the approximated Voronoi diagram V. Resting cells are blank, excited cells red (dark grey), refractory cells grey, and precipitate cells black. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
484
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
neighbourhoods are shown in Fig. 2a. Distributions of numbers of immediate and second neighbours and sizes of second-order neighbourhood u are given in Fig. 2b. Predominant number, over 40%, of cells have six first-order (immediate) neighbours. Almost half of the cells have 13 or 14 second-order neighbours. For any cell size of second-order neighbourhood u is a sum of first-order neighbours and second-order neighbours. For Voronoi automata with high-density packing of 15 K cells, half of the cells have either 19, or 20 or 21 neighbours. If we consider only immediate neighbourhood then V can be, in principle, seen as a slightly distorted, hexagonal lattice. However, by adopting second-order neighbourhood we
are moving to a random structure with higher than hexagonal lattice coordinate number (a node in hexagonal lattice has 18-node second-order neighbourhood while dominating neighbourhood sizes in V are 19, 20 and 21). Transition from the excited to the refractory state is unconditional, i.e. takes place regardless of the states of the cell’s neighbours. A resting cell excites if it has at least one excited neighbour. In this particular model we take a refractory state and precipitate state as absorbing: once a cell takes either of these two states it does not update its state any longer. Neighbourhood sizes may differ between Voronoi cells therefore we use a state transition function, where a cell
Fig. 5. Quality of approximation depends on g. (a) Voronoi diagram VðBÞ of seven planar points set B constructed by Fortuna’s sweep line algorithm, (b)–(o) stationary configurations of precipitate in automaton V excited by B for different g = 0.2. . .0.525. In snapshots (b)–(o) black pixels symbolise Voronoi cells of V in precipitate state, cells in refractory states are blank.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
updates its state depending on a relative excitation in its neighbourhood. We assume that precipitation occurs in a resting cell when a ratio of excited neighbours to a total number of neighbours exceeds some threshold g 2 [0, 1]. Let r(V(x)t) be a number of excited cells in the cell V(x)’s second-order neighbourhood, rðVðxÞt Þ ¼ P t VðyÞ2uðVðxÞÞ jfVðyÞ : VðyÞ ¼ þgj then a cell updates its state by the following rule:
VðxÞtþ1
8 #; if VðxÞt ¼ and rt ðxÞ=mðxÞ > g > > > < þ; if VðxÞt ¼ and rt ðxÞ=mðxÞ > 1 ¼ > ; if VðxÞt ¼ þ > > : ; otherwise
Example of propagating wave-fronts are shown in Fig. 3. For g = 0.3 excitation wave-front is convex, shaping into planar by the end of propagation (Fig. 3b). In case of lower threshold, g = 1, wave near the edges of moves quicker then inside the diagram core (Fig. 3a). Thus wave front becomes concave. The part of wave front travelling along edges of diagram reaches the side opposite to the initial perturbation side quicker then the front propagating inside the diagram. Thus wave front becomes closed (Fig. 3a). The domain surrounded by a nested group of target waves, in the western part of the diagram in Fig. 3a, is the place where the excitation waves travelling from east collapse. Such patterns of excitation are somewhat universal and can be observed in various systems, even as calcium waves in coupled cells induced by internal noise [35]. In computational experiments we construct set P by filling a disc-container of radius 480 units with up to 15 K points. The points are packed at random but there is
485
at least 5 units distance between any two points. Voronoi diagram V – on which Voronoi automaton is constructed – is calculated by a classical sweep line algorithm [21].
4. Constructing Voronoi diagram on Voronoi automata Let B be a set of planar points on which a Voronoi diagram V on V; g ¼ 0:4, is approximated. We project B onto V and excite cells of V which are closer than 9 units to points of B. Excitation waves spread on V. The excitation waves collide and precipitation occurs near the sites of the waves’ collisions. Configuration of cells in precipitate state represents edges of VðBÞ. An example of excitationprecipitation dynamics in V which approximates VðBÞ, where B is a planar set of seven points is shown in Fig. 4. A scheme of the Voronoi diagram computed is shown in Fig. 8a. Even in a single example shown in Fig. 4 it is clear that V does not calculate VðBÞ precisely, it rather approximates it with discrete domains of precipitation, and also introduces some noise (sites of precipitation not coinciding with bisectors of VðVÞ). Only threshold of precipitation g can be varied in our model. How does the degree of approximation and signal to noise ratio depend on the value of g? To answer we undertook a series of computational experiments, illustrated – for some values of g – in Fig. 5. For a set B of seven points we constructed ‘classical’ planar Voronoi diagrams VðBÞ (Fig. 5a). Then for g = 0.2. . .0.525 we excited V with B, waited until V reaches its stationary configuration, where all cells are in either precipitate or refractory state, and compared the configurations of V and the image of VðBÞ. For each g we recorded
Fig. 6. Characterisation of approximation of VðBÞ (Fig. 5a) by automaton V: (a) q vs g and (b) m vs g. Markers (rhombs and circles) show actual values of m(g) and q(g), while solid lines represent a moving average trend line purely for guidance.
486
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
a fraction q of precipitate-cells in a final configuration of VðgÞ to a total number of precipitate-cells in the same configuration (Fig. 6a).
We found that while density m, calculated as a ratio of the number of precipitate-cells in a stationary configuration VðgÞ versus the number of precipitate-cells in
Fig. 7. Approximation of Voronoi diagram VðOÞ of planar shapes O in Voronoi automaton V. Voronoi cell of V in precipitate state are black, cells in refractory states are blank. Excited cells are red (grey). (a) representation of O, (b)–(f) development of V. VðOÞ is represented by black pixels in (f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
487
Fig. 8. Schematic representation of data objects and their Voronoi diagram or skeleton for precipitation-based approximation shown in (a) Fig. 4, (b) Fig. 7, and (c) Fig. 9.
configuration of Vð0:2Þ, exponentially decreases with increase of precipitation threshold g (Fig. 6b) degree of approximation q shows more intriguing behaviour (Fig. 6a). From g = 0.2 till g = 0.41 ‘signal-to-noise’ ratio q(g) grows almost exponentially. Trend of q(g) undergoes S-type transition between g = 0.41 and g = 0.453 with one peak q(0.44) and two cavities q(0.43) and q(0.452). As illustrated in Fig. 5 V with precipitation threshold in the zone around g = 0.4, onset of ‘strange’ behaviour of q(g), demonstrates best approximation of VðBÞ, as recognised by eye.
5. Arbitrary-shaped planar objects and contours The Voronoi automaton V copes well with data set O of planar finite shapes. To approximate diagram V of arbitrary planar shapes we project the objects of O (Fig. 7a) onto V in such a manner that at time t = 0 all cells of V covered (coinciding with pixels of) by shapes from O becomes excited (Fig. 7b). Waves of excitation travel across V and provoke precipitation during their collisions (Fig. 7c–e). When all cells take refractory or precipitate state the automaton’s configuration becomes stationary. The spatial distribution of precipitate represents a Voronoi diagram of the data shapes (see Fig. 7f and Fig. 8b). A precipitateenhancement of the original position of the data shape is a byproduct of V’s development. A skeleton of a planar contour is a set of centres of bitangent circles lying inside the contour [15]. Blum’s grass-fire algorithm for computing skeleton employs propagating patterns [12,13,15]: to compute a skeleton we set a contour on fire and let the fire spread, quenched points where the advancing fire-fronts collide represent the skeleton. Most known algorithms of skeletonisation are based on Blum’s approach: simulation of grass-fire [29], distance transform [38], analytical construction of medial axis and topological thinning [34,10]. To approximate a skeleton of a shape S (Fig. 9a) we project S onto V. Voronoi cells of V corresponding to black pixels of S are excited (Fig. 9b). Excitation waves propagate inside and outside of the domain S of initial excitation;
we can ignore behaviour of outward waves. Waves travelling inside the contour S trigger precipitation at the sites of the waves’ interaction (Fig. 9c–e). Distribution of precipitate (Fig. 9e) inside S approximates skeleton of S (Fig. 8c).
6. Summary We demonstrated that is possible to solve some problems of computational geometry on a discrete model of irregular arrangements of non-uniform vesicles filled with excitable chemical mixtures [6,7,25]. We have shown that a threshold of relative local excitation density can be used to parameterise excitation and precipitation dynamics, enhance results of the computation and reduce signalto-noise ratio. We proved that the same optimal threshold of precipitation works well not only for approximating Voronoi diagram of planar sets but also for arbitrary shapes and skeletonisation of planar contours. In terms of automata-network based computation we advanced our previous results on excitation in Delaunay automata [4] and b-skeletons [5,8]. The parameterisation developed could be used in designing experimental laboratory prototypes of finegrained compartmentalised excitable chemical processors, e.g. using micro-emulsion approach [26], and developing nano-scale massively parallel computers [11]. Namely, we have already designed computers models of lightsensitive BZ-vesicles packed into a hexagonal lattice. The vesicles implement a set of universal logical gates and also one-bit binary adder [6]. Our parameterisation approach is used to select levels of illumination, which in turn controls level of the medium’s excitability, and size of pores connecting BZ-vesicles, which determines propagation of wave-fragments. With regards to laboratory implementation in mind we must mention that our approach is not the only one by organically complement early developed techniques of automata-based simulation of reaction–diffusion chemical systems. For example, lattice gas automata is amongst highly used computer simulation tools [28,14]. The approach based on Voronoi automata, discussed in present
488
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
Fig. 9. Approximation of a skeleton of a star-contour S by V. Voronoi cell of V in precipitate state are black, cells in refractory states are blank. Excited cells are red (grey). (a) representation of a star-contour, (b)–(e) development of V. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
paper, is less accurate than that based on lattice gas automata however provides higher computational efficiency
with degree of accuracy high enough for rapid prototyping of massively-parallel chemical processors.
A. Adamatzky et al. / Chaos, Solitons & Fractals 44 (2011) 480–489
Acknowledgement The work is part of the European project 248992 funded under 7th FWP (Seventh Framework Programme) FET Proactive 3: Bio-Chemistry-Based Information Technology CHEM-IT (ICT-2009.8.3). References [1] Adamatzky A. Reaction–diffusion algorithm for constructing discrete generalised Voronoi diagram. Neural Netw World 1994;9:6635–43. [2] Adamatzky A. Voronoi-like partition of lattice in cellular automata. Math Comput Model 1996;23:51–66. [3] Adamatzky A, De Lacy Costello B, Asai T. Reaction–diffusion computers. Elsevier; 2005. [4] Adamatzky A. On excitable Delaunay automata (Kybenertes) 2011, in press.
489
[19] De Lacy Costello BPJ, Jahan I, Adamatzky A, Ratcliffe NM. Chemical tessellations. Int J Bifurcat Chaos 2009;19:619–22. [20] Filatovs GJ. Delaunay-subgraph statistics of disc packings materials characterisation 1998;40(1):27–35. [21] Fortune S. A sweep line algorithm for Voronoi diagrams. In: Proceedings of the 2nd annual symposium on computational geometry; 1986. p. 313–22. [22] Gervois A, Annic C, Lemaitre J, Ammi m, Oger L, Troadec J-P. Arrangement of discs in 2D binary assemblies. Physica A 1995;218:403–18. [23] Gorecki J. Private communication; 2010. [24] Hobbs WL. Network topology in aperiodic networks. J Non-Crystal Solids 1995;192/193:79–91. [25] Holley J, Adamatzky A, Bull L, De Lacy Costello B, Jahan I. Modalities of Belousov-Zhabotinsky encapsulated vesicles. Nanocommunication Netw J; in press.