Geometry and Morphology of the Cosmic Web:
Analyzing Spatial Patterns in the Universe
Abstract
We review the analysis of the Cosmic Web by means of an extensive toolset based on the use of Delaunay and Voronoi tessellations. The Cosmic Web is the salient and pervasive foamlike pattern in which matter has organized itself on scales of a few up to more than a hundred Megaparsec. The weblike spatial arrangement of galaxies and mass into elongated filaments, sheetlike walls and dense compact clusters, the existence of large nearempty void regions and the hierarchical nature of this mass distribution are three major characteristics of the comsic matter distribution.
First, we describe the Delaunay Tessellation Field Estimator.Using the unique adaptive qualities of Voronoi and Delaunay tessellations, DTFE infers the density field from the (contiguous) Voronoi tessellation of a sampled galaxy or simulation particle distribution and uses the Delaunay tessellation as adaptive grid for defining continuous volumefilling fields of density and other measured quantities through linear interpolation. The resulting DTFE formalism is shown to recover the hierarchical nature and the anisotropic morphology of the cosmic matter distribution. The Multiscale Morphology Filter (MMF) uses the DTFE density field to extract the diverse morphological elements  filaments, sheets and clusters  on the basis of a ScaleSpace analysis which searches for these morphologies over a range of scales. Subsequently, we discuss the Watershed Voidfinder (WVF), which invokes the discrete watershed transform to identify voids in the cosmic matter distribution. The WVF is able to determine the location, size and shape of the voids. The watershed transform is also a key element in the SpineWeb analysis of the cosmic matter distribution. Finding its mathematical foundation in Morse theory, it allows the determination of the filamentary spine and connected walls in the cosmic matter density field through the identification of the singularities and corresponding separatrices. The first results of a direct implementation on the Delaunay tessellation itself are presented. Finally, we describe the concept of Alphashapes for assessing the topology of the cosmic matter distribution.
I Introduction: the Cosmic Web
The large scale distribution of matter revealed by galaxy surveys features a complex network of interconnected filamentary galaxy associations. This network, which has become known as the Cosmic Web [1], contains structures from a few megaparsecs^{1}^{1}1The main measure of length in astronomy is the parsec. Technically a parsec is the distance at which we would see the distance EarthSun at an angle of 1 arcsec. It is equal to 3.262 lightyears . Cosmological distances are substantially larger, so that a Megaparsec () is the regular unit of distance. Usually this goes along with , the cosmic expansion rate (Hubble parameter) in units of km/s/Mpc (). up to tens and even hundreds of Megaparsecs of size. Galaxies and mass exist in a wispy weblike spatial arrangement consisting of dense compact clusters, elongated filaments, and sheetlike walls, amidst large nearempty void regions, with similar patterns existing at earlier epochs, albeit over smaller scales. The hierarchical nature of this mass distribution, marked by substructure over a wide range of scales and densities, has been clearly demonstrated. Its appearance has been most dramatically illustrated by the recently produced maps of the nearby cosmos, the 2dFGRS, the SDSS and the 2MASS redshift surveys [2, 3, 4]^{2}^{2}2Because of the expansion of the Universe, any observed cosmic object will have its light shifted redward: its redshift . According to Hubble’s law, the redshift is directly proportional to the distance of the object, for : (with the velocity of light, and the Hubble constant). Because it is extremely cumbersome to measure distances directly, cosmologists resort to the expansion of the Universe and use as a distance measure. Because of the vast distances in the Universe, and the finite velocity of light, the redshift of an object may also be seen as a measure of the time at which it emitted the observed radiation. .
The vast Megaparsec cosmic web is one of the most striking examples of complex geometric patterns found in nature, and certainly the largest in terms of sheer size. Computer simulations suggest that the observed cellular patterns are a prominent and natural aspect of cosmic structure formation through gravitational instability [5], the standard paradigm for the emergence of structure in our Universe.
Ia Dynamical Evolution of the Cosmic Web
At least three “universal” characteristics of the resulting nonlinear cosmic matter distribution can be recognized in large Nbody simulations of cosmic structure formation such as the Millennium simulation [6]. One prominent aspect is the hierarchical clustering. The spatial cosmic matter distribution is marked by a large range of scales and densities. It is the product of an evolution in which the first objects to condense are small, with larger structures forming through the gradual merging of smaller structures. Another prominent feature is its anisotropic and weblike spatial geometry, the consequence of the gravitational tendency of overdensities to collapse in an anisotropic fashion. It finds its origin in the intrinsic flattening of the overdensities in the primordial density field [7], augmented by the anisotropy of the gravitational force field induced by the external matter distribution (i.e. by tidal forces). A third important characteristic is the dominant presence of Voids, large roundish underdense, often nearempty, regions. They form in and around density troughs in the primordial density field. They have an essential role in the organization of the cosmic matter distribution. Recently, their emergence and evolution has been explained within the context of hierarchical gravitational scenarios [8].
The Cosmic Web theory of Bond et al. [1] succeeded in synthesizing all relevant aspects into a coherent dynamical and evolutionary framework. Instrumental is the realization that the outline of the cosmic web may already be recognized in the primordial density field. The statistics of the primordial tidal field explains why the large scale universe looks predominantly filamentary and why in overdense regions sheetlike membranes are only marginal features. Of key importance is the realization that the rare high peaks which will eventually emerge as clusters are the dominant agents for generating the large scale tidal force field: it is the clusters which weave the cosmic tapestry of filaments [1, 9, 10]. They cement the structural relations between the components of the Cosmic Web and themselves form the junctions at which filaments tie up. This relates the strength and prominence of the filamentary bridges to the proximity, mass, shape and mutual orientation of the generating cluster peaks: the strongest bridges are those between the richest clusters that stand closely together and point into each other’s direction. Through the direct relation between the Cosmic Web and the primordial tidal field one may understand why the large scale universe looks predominantly filamentary and why in overdense regions sheetlike membranes are only marginal features [11].
The emerging picture is one of a primordially and hierarchically defined network whose weblike topology is imprinted over a wide spectrum of scales. Weblike patterns on ever larger scales get to dominate the density field as cosmic evolution proceeds, and as small scale structures merge into larger ones. Within the gradually emptying void regions, however, the topological outline of the early weblike patterns remains largely visible.
IB Tessellations and Web Analysis
Despite a large variety of attempts, as yet no generally accepted descriptive framework has emerged for the objective and quantitative analysis of the Cosmic Web. Despite the multitude of elaborate qualitative descriptions it has remained a major challenge to characterize its structure, geometry and topology. Many attempts to describe, let alone identify, the features and components of the Cosmic Web have been of a rather heuristic nature. The overwhelming complexity of both the individual structures as well as their connectivity, the lack of structural symmetries, its intrinsic multiscale nature and the wide range of densities that one finds in the cosmic matter distribution has prevented the use of simple and straightforward toolboxes.
In the observational reality galaxies are the main tracers of the cosmic web and it is mainly through the measurement of the redshift distribution of galaxies that we have been able to map its structure. Likewise, simulations of the evolving cosmic matter distribution are almost exclusively based upon Nbody particle computer calculations, involving a discrete representation of the features we seek to study. Both the galaxy distribution as well as the particles in an Nbody simulation are examples of spatial point processes in that they are discretely sampled and have an irregular spatial distribution.
For furthering our understanding of the Cosmic Web, and to investigate its structure and dynamics, it is of prime importance to have access to a set of proper and objective analysis tools. In this contribution we follow the finding that spatial tessellations – in particular Voronoi and Delaunay tessellations – generated by the discretely sampled cosmic density and velocity fields form an ideal basis for an elaborate toolset that succesfully deals with several challenging aspects of the analysis of the Cosmic Web:

Tessellations may be used to interpolate a sample of discrete and irregularly distributed values into a volumeweighted and volumecovering continuous field.

In case the point sample is a representative reflection of an underlying smooth and continuous density/intensity field, Voronoi tessellations can be used to estimate and reconstruct the density field.

Tessellations are highly sensitive to the morphology of the spatial point distribution. As a result they sensitively probe the key aspects  hierarchical matter distribution, anisotropic patterns, and voids – of the Cosmic Web without resorting to userdefined parameters or functions, and without affecting any of the other essential characteristics.
These issues are all addressed by the Delaunay Tessellation Field Estimator (DTFE). The DTFE technique, developed by Schaap & van de Weygaert [12], forms an elaboration of the velocity interpolation scheme introduced by Bernardeau & van de Weygaert [13] towards a multidimensional and fully adaptive scheme for the estimation and interpolation of density, velocity and other nonuniformly and discretely sampled quantities to yield a corresponding volumecovering continuous field.
The application of DTFE to the observed spatial distribution of galaxies, or that of particles in a computer Nbody simulation, produces a continuous density and velocity field which retains the crucial aspects of its hierarchical nature and anisotropic morphology. This may be directly appreciated from the three consecutive zoomins of the DTFE rendered cosmic web density field in fig. 6, and its comparison to the more conventional gridbased TSC interpolation scheme. Because of these characteristics, the DTFE density field forms the basis for a diverse set of tools for the analysis of various aspects of the cosmic matter distribution. In this contribution we will report on our work on the Multiscale Morphology Filter for detecting filaments and sheets, the Watershed Void Finder for detecting voids, the Cosmic Spine formalism for the analysis of the filamentary network of the Cosmic Web and, finally, the related Alphashape analysis of the topological structure of the cosmic matter distribution.
Ii DTFE: fundamentals
DTFE obtains optimal local estimates of the spatial density (see [14], sect. 8.5), while the tetrahedra of its dual Delaunay tessellation are used as multidimensional intervals for linear interpolation of the field values sampled or estimated at the location of the sample points (see [14], ch. 6). With respect to the aspect of interpolation, the Delaunay Tessellation Field Estimator is the linear version of nnneighbour interpolation techniques (see sect. IIB). Its core, however, is the additional aspect of being able to translate the generating spatial point distribution into a continuous density field.
The DTFE method has been first defined in the context of a description and analysis of cosmic flow fields which are sampled by a set of discretely and sparsely sampled galaxy peculiar velocities. Bernardeau & van de Weygaert [13] demonstrated the method’s superior performance with respect to conventional interpolation procedures. They also proved that the obtained field estimates involve those of the proper volumeweighted quantities, instead of the usually implicit massweighted quantities (see sect. IIE). This corrected a few fundamental biases in estimates of higher order velocity field moments.
The DTFE technique allows us to follow the same geometrical and structural adaptive properties of the higher order nnneighbour methods while allowing the analysis of truely large data sets and include the reconstruction of the most important aspect, that of the cosmological density field. For threedimensional samples with large number of points, akin to those found in large cosmological computer simulations, the more complex geometric operations involved in the pure nnneighbour interpolation still represent a computationally challenging task (see sect. IIB). To deal with the large point samples consisting of hundreds of thousands to several millions of points we chose to follow a related nnneigbhour based technique that restricts itself to pure linear interpolation.
Iia Voronoi and Delaunay Adaptivity
The primary ingredient of the DTFE procedure, and the related Natural Neighbour Interpolation methods, is the Delaunay tessellation of the particle distribution . The Delaunay and Voronoi cells adjust themselves to the spatial characteristics of the point distribution. This concerns its spatial resolution and its geometry.
DTFE, and its higherorder equivalent Natural Neigbhour Interpolation [15, 16, 17, 18, 19], exploit three properties of Voronoi and Delaunay tessellations (see eg. [20, 21]):

The tessellations are very sensitive to the local point density. DTFE uses this to define a local estimate of the density on the basis of the inverse of the volume of the tessellation cells.

The sensitivity of Voronoi and Delaunay tessellations to the local geometry of the point distribution. This allows DTFE to accurately trace anisotropic features, such as seen in the Cosmic Web.

The adaptive and minimum triangulation properties of Delaunay tessellations allow them to be used as adaptive spatial interpolation intervals for irregular point distributions. This has also been recognized in a
In sparsely sampled regions the Voronoi and Delaunay cells are large and the distance between natural neighbours, points that share a Voronoi simplex, is large. Not only the size, but also the shape of the Delaunay tetrahedra is fully determined by the spatial point distribution: if is anisotropic this will be reflected in the distribution. In fig. 2 we see that the Delaunay tessellation traces the density and geometry of the local point distribution at three consecutive spatial scales to a remarkable degree. The density and geometry of the local point distribution will therefore determine the resolution of the spatial interpolation and reconstruction procedures based on the Voronoi and Delaunay tessellation.
IiB Natural Neighbour Interpolation
Natural Neighbour Interpolation formalism is a generic higherorder multidimensional interpolation, smoothing and modelling procedure utilizing the concept of natural neighbours to obtain locally optimized measures of system characteristics. Its theoretical basis was developed and introduced by Sibson [16], while extensive treatments and elaborations of nninterpolation may be found in [17, 19]. Natural neighbour interpolation produces a conservative, artificefree, result by finding areaweighted weighted averages, at each interpolation point, of the functional values associated with that subset of data which are natural neighbors of each interpolation point. According to the nninterpolation scheme the interpolated value at a position is given by
(1) 
in which the summation is over the natural neighbours of the point , i.e. the sample points with whom the order2 Voronoi cells are not empty. Sibson interpolation is based upon the interpolation kernel to be equal to the normalized order2 Voronoi cell,
(2) 
in which is the volume of the potential Voronoi cell of point if it had been added to the point sample and the volume concerns the order2 Voronoi cell , the region of space for which the points and are the closest points. The interpolation kernels are always positive and sum to one. The resulting function is continuous everywhere within the convex hull of the data, and has a continuous slope everywhere except at the data themselves.
IiC DTFE Interpolation
The linear interpolation scheme of DTFE exploits the same spatially adaptive characteristics of the Delaunay tessellation generated by the point sample as that of regular natural neighbour schemes. For DTFE the interpolation kernel is that of regular linear interpolation within the Delaunay tetrahedron in which is located (see eq. 5),
(3) 
in which the sum is over the four sample points defining the Delaunay tetrahedron. Note that for both the nninterpolation as well as for the linear DTFE interpolation, the interpolation kernels are unity at sample point location and equal to zero at the location of the other sample points ,
(4) 
where is the location of sample point .
In practice, it is convenient to replace eqn. 3 with its equivalent expression in terms of the (linear) gradient inside the Delaunay simplex ,
(5) 
The value of can be easily and uniquely determined from the field values at the sample points constituting the vertices of a Delaunay simplex. Given the location of the four points forming the Delaunay tetrahedra’s vertices, , , and , and the value of the sampled field at each of these locations, , , and and defining the quantities
(6)  
as well as the gradient follows from the inversion
(7)  
Once the value of has been determined for each Delaunay tetrahedron in the tessellation, it is straightforward to determine the DTFE field value for any location by means of straightforward linear interpolation within the Delaunay tetrahedron in which is located (eqn. 5).
The one remaining complication is to locate the Delaunay tetrahedron in which a particular point is located. This is not as trivial as one might naively think. It not necessarily concerns a tetrahedron of which the nearest nucleus is a vertex. Fortunately, a very efficient method, the walking triangle algorithm [22, 23] has been developed. Details of the method may be found in [24, 20].
IiD DTFE Density Estimates
The DTFE procedure extends the concept of interpolation of field values sampled at the point sample to the estimate of the density from the spatial point distribution itself. This is only feasible if the spatial distribution of the discrete point sample forms a fair and unbiased reflection of the underlying density field.
It is commonly known that an optimal estimate for the spatial density at the location of a point in a discrete point sample is given by the inverse of the volume of the corresponding Voronoi cell (see [14], for references). Tessellationbased methods for estimating the density have been introduced by Brown [25] and Ord [26]. In astronomy, Ebeling & Wiedenmann [27] were the first to use tessellationbased density estimators for the specific purpose of devising source detection algorithms. This work has recently been applied to cluster detection algorithms by [28, 29, 30, 31, 32]. Along the same lines, Ascasibar & Binney [33] suggested that the use of a multidimensional binary tree might offer a computationally more efficient alternative. However, these studies have been restricted to raw estimates of the local sampling density at the position of the sampling points and have not yet included the more elaborate interpolation machinery of the DTFE and Natural Neighbour Interpolation methods.
The density field reconstruction of the DTFE procedure consists of two steps, the zerothorder estimate of the density values at the location of the points in and the subsequent linear interpolation of those zerothorder density estimates over the corresponding Delaunay grid throughout the sample volume. This yields the DTFE density field estimate .
An essential requirement for cosmological purposes of our interpolation scheme is that the estimated DTFE density field should guarantee mass conservation: the total mass corresponding to the density field should be equal to the mass represented by the sample points. Indeed, this is an absolutely crucial condition for many applications of a physical nature. Since the mass is given by the integral of the density field over space, this translates into the integral requirement
(8)  
with is the mass per sample point. It is straightforward to infer that if the zerothorder estimate of the density values would be the inverse of the regular Voronoi volume the condition of mass conservation would not be met. Instead, the DTFE procedure employs a slightly modified yet related zerothorder density estimate, the normalized inverse of the volume of the contiguous Voronoi cell of each point . For dimensional space this is
(9) 
The contiguous Voronoi cell of a point is the union of all Delaunay tetrahedra of which point forms one of the four vertices. It is straightforward to prove that with the contiguous Voronoi cell density estimator mass conservation is guaranteed (see [21], sect. 10).
IiE Volumeweighted and Massweighted fields
For relating measured quantities to theoretical predictions, one has to take account of the implicit filtering process involved in translating measurements (including that in computer simulations) to interpretable measures. When dealing with a probe of an underlying density field  e.g. the galaxy distribution in a cosmological context  there are two possibilities. Theoretically the cleanest measure is that of the volumeweighted quantity ,
(10) 
In practice, however, most techniques implicitly (and often unknowingly) yield the massweighted field averages.
(11) 
where is the adopted filter function defining the weight of a mass element in a way that is dependent on its position with respect to the position . It turns out that the use of the Voronoi and Delaunay tessellation guarantees the volumeweighted nature of the DTFE fields [13].
An additional crucial ingredient of any reconstruction procedure are its error characteristics. Within the limited context of the present contribution it is not feasible to include a detailed discussion of the noise and error characteristics of DTFE reconstructed density and velocity fields. For this we refer the reader elsewhere [20, 21].
Iii the DTFE procedure
The complete DTFE reconstruction procedure is shown in the the schematic diagram of fig. 5. It involves the following steps:

Point sample
Defining the spatial distribution of the point sample. 
Delaunay Tessellation
Construction of the Delaunay tessellation from the point sample. 
Field values point sample
Dependent on whether it concerns the densities at the sample points or a measured field value there are two options:
General (nondensity) field:
(Sampled) value of field at sample point. 
Density field


Field Gradient
Calculation of the field gradient estimate in each dimensional Delaunay simplex (: tetrahedron; : triangle) by solving the set of linear equations for the field values at the positions of the tetrahedron vertices,(12)

Interpolation.
The final basic step of the DTFE procedure is the field interpolation. The processing and postprocessing steps involve numerous interpolation calculations, for each of the involved locations . Given a location , the Delaunay tetrahedron in which it is embedded is determined. On the basis of the field gradient the field value is computed by (linear) interpolation,(13) In principle, higherorder interpolation procedures are also possible.

Processing.
Though basically of the same character for practical purposes we make a distinction between straightforward processing steps concerning the production of images and simple smoothing filtering operations on the one hand, and more complex postprocessing on the other hand. The latter are treated in the next item. Basic to the processing steps is the determination of field values following the interpolation procedure(s) outlined above.
Straightforward “first line” field operations are “Image reconstruction” and, subsequently, “Smoothing/Filtering”.
Image reconstruction.
For a set of image points, usually grid points, determine the image value: formally the average field value within the corresponding gridcell. 
Smoothing and Filtering:


Postprocessing.
The real potential of DTFE fields may be found in sophisticated applications, tuned towards uncovering characteristics of the reconstructed fields. An important aspect of this involves the analysis of structures in the density field. Some notable examples are:
Advanced filtering operations. Potentially interesting applications are those based on the use of wavelets [34].

Tracing the filamentary spine of the Cosmic Web, on the basis of a Morsetheory related analysis, the Cosmic Spine [39].

Halo detection in Nbody simulations [32].

The computation of 2D surface densities for the study of gravitational lensing [40].

In addition, DTFE enables the simultaneous and combined analysis of density
fields and other relevant physical fields. As it allows the simultaneous
determination of density and velocity fields, it can serve as
the basis for studies of the dynamics of structure formation in the cosmos.
Its ability to detect substructure as well as reproduce the morphology
of cosmic features and objects implies DTFE to be suited for assessing their
dynamics without having to invoke artificial filters.
Iv DTFE and the Cosmic Web
Figure 4 provides a nice impression of the versatility of the DTFE formalism. DTFE density and velocity fields may be depicted at any arbitrary resolution without involving any extra calculation: zoomins represent themselves a real magnification of the reconstructed fields. This is in stark contrast to conventional reconstructions in which the resolution is arbitrarily set by the users and whose properties are dependent on the adopted resolution.
It shows the density and velocity field, at two resolutions, in and around filaments and other components of the cosmic web in a GIF CDM simulation. DTFE not only it allows a study of the patterns in the nonlinear matter distribution but also a study of the related velocity flows. Because DTFE manages to follow both the density distribution and the corresponding velocity distribution into nonlinear features it opens up the window towards a study of the dynamics of the formation of the cosmic web and its corresponding elements [41]. Note that such an analysis of the dynamics is limited to regions and scales without multistream flows.
Iva DTFE characteristics
Within the cosmological context a major – and crucial – characteristic of a processed DTFE density and velocity field is that it is capable of delineating the three fundamental characteristics of the spatial structure of the Megaparsec cosmic matter distribution:

It outlines the full hierarchy of substructures present in the sampling point distribution, relating to the standard view of structure in the Universe having arisen through the gradual hierarchical buildup of matter concentrations (fig. 6).

DTFE also reproduces any anisotropic patterns in the density distribution without diluting their intrinsic geometrical properties. This is particularly important when analyzing the prominent filamentary and planar features marking the Cosmic Web (fig. 7).

A third important aspect of DTFE is that it outlines the presence and shape of voidlike regions. Because of the interpolation definition of the DTFE field reconstruction voids are rendered as regions of slowly varying and moderately low density values.
The success of DTFE in following these key characteristics has been substantiated by quantitative tests. In [42] we demonstrate DTFE’s ability to infer the fractal dimensions of a range of fractal point distributions: DTFE follows the density distribution throughout the available spatial and mass range. Additional tests in [43] have confirmed the ability of DTFE to reconstruct density fields with the correct shape morphology (inertia tensor), whether it related to distributions dominated by compact high density cluster peaks, elongated filaments or sheets. Its comparative performance with respect to other interpolation schemes may be assessed from fig. 7. Clearly, DTFE remains closer to the actual structure represented in the particle distribution.
In a thorough analysis of the interpolation and morphological performance of DTFE, Platen et al. [44] has demonstrated that DTFE produces better and more reliable results than Natural Neighbour interpolation and various Kriging schemes. Even while the latter do take into account longrange spatial correlations, the selfadaptive local character of DTFE appears not only to be better suited for uniformly sampled spatial datasets, but also for dealing with datasets with a radially varying selection function or datasets beset by systematic redshift distortions due to the peculiar velocities of galaxies^{3}^{3}3Redshift distortions: in addition to their cosmic expansion velocity, galaxies have velocities with respect to the expanding Universe. The Doppler redshift corresponding to these velocities add to the cosmic redshift to yield a redshift that is not precisely proportional to the distance of an object..
IvB Analysis of the Cosmic Web
Within its cosmological context, DTFE will meet its real potential in more sophisticated applications tuned towards uncovering morphological characteristics of the reconstructed spatial patterns. The true potential of DTFE and related adaptive random tessellation based techniques comes to the fore when we take DTFE reconstructions as the basis for a variety of “postprocessing” tools. A variety of recent techniques have recognized the high dynamic range and adaptivity of tessellations to the spatial and morphological resolution of the systems they seek to analyze.
We have been working on three different, yet mutually related, formalisms which use DTFE density fields as the basis for the identification of various key aspects of the Cosmic Web. The Multiscale Morphology Filter (MMF), introduced and defined by AragónCalvo et al. [36], attempts to detect weblike anisotropic features over a range of spatial scales and implicitly incorporates the hierarchical nature of the Megaparsec matter distribution. While incorporating the use of a more or less artificial scale filter, the Watershed Void Finder (WVF, [37]) and the Cosmic Spine technique [39] turn to the topological structure of the underlying density field to identify void regions and trace the pervasive weblike filamentary network. Both techniques are manifest translations of concepts from Morse Theory, which describes the topological structure of a density field on the basis of its singular points  minima, maxima and saddle points.
V the Multiscale Morphology Filter
The Multiscale Morphology Filter (MMF), introduced by AragónCalvo et al. [36], is used for the identification and characterization of different morphological elements of the large scale matter distribution in the Cosmic Web. An example of the morphological structure of the matter distribution of a cosmological Nbody simulation can be seen in fig. 8 [36]. The image shows the identified elongated filaments (dark grey), sheetlike walls (light grey) and clusters (black dots) in the weblike pattern of the simulation.
The Multiscale Morphology Filter (MMF) method has been developed on the basis of visualization and feature extraction techniques in computer vision and medical research [45]. The technology finds its origin in computer vision research and has been optimized within the context of feature detections in medical imaging. Frangi et al. [46] and Sato et al. [47] presented its operation for the specific situation of detecting the web of blood vessels in a medical image. This defines a notoriously complex pattern of elongated tenuous features whose branching make it closely resemble a fractal network.
The MMF dissects the cosmic web on the basis of the multiscale analysis of the Hessian of the density field. Its formulation in terms of its density field scalespace structure implicitly takes into account the hierarchical nature of the matter distribution. Scale Space analysis looks for structures of a mathematically specified type in a hierarchical, scale independent, manner. It is presumed that the specific structural characteristic is quantified by some appropriate parameter (e.g.: density, eccentricity, direction, curvature components). The data is filtered to produce a hierarchy of maps having different resolutions, and at each point, the dominant parameter value is selected from the hierarchy to construct the scale independent map. We refer to this scalefiltering processes as a Multiscale morphology filter.
Va Scale Space
Crucial for the ability of the method to identify anisotropic features such as filaments and walls is the use of a morphologically unbiased and optimized continuous density field retaining all features visible in a discrete galaxy or particle distribution. Unless one manages to work directly from the sampled particle or galaxy distribution, possibly via its Voronoi or Delaunay tessellation, it is therefore imperative to translate the discrete particle distribution into its DTFE density field, The morphological intentions of the MMF method render DTFE a key element for translating the particle or galaxy distribution into a representative continuous density field .
Following the determination of the DTFE density field , it is smoothed over a range of scales by means of a hierarchy of spherically symmetric Gaussian filters having different widths :
(14) 
where denotes a Gaussian filter of width :
(15) 
A pass of the smoothing filter attenuates structure on scales smaller than the filter width. The MMF scalespace analysis involves a discrete number of levels, . The basescale is taken to be equal to the pixel scale of the raw DTFE density map, while the smoothing radii of the other levels is chosen such that essential (sub)structures and features in the density field are resolved. We denote the level smoothed version of the DTFE reconstructed field by the symbol . The Scale Space itself is constructed by stacking these variously smoothed data sets, yielding the family of smoothed density maps : . In our implementation [36], we found it sufficient to use only levels, although this may be easily expanded.
VB Density Field Hessian & Eigenvalues
A data point can be viewed at any of the scales where scaled data has been generated. The crux of the concept is that the neighbourhood of a given point will look different at each scale. There are potentially many ways of making a comparison of the scale dependence of local environment. The MMF employs the Hessian Matrix of the local density distribution in each of the smoothed replicas of the original data,
where and is the Kronecker delta. In other words, the scale space representation of the Hessian matrix for each level is evaluated by means of a convolution with the second derivatives of the Gaussian filter, also known as the Marr (or, less appropriately, “Mexican Hat”) Wavelet.
At each point of a dataset in the Scale Space view of the data we can quantify the local “shape” of the density field in the neighbourhood of that point by calculating at each point the eigenvalues of the Hessian Matrix of the data values,
where the eigenvalues are arranged so that . The are coordinate independent descriptors of the behaviour of the density field in the locality of the point and can be combined to create a variety of morphological indicators: they determine the local morphological signal. The corresponding eigenvectors show the local orientation of the morphology characteristics.
Structure  ratios  constraints 

Blob  
Line  
Sheet 
VC Morphological Identification
Dependent on the eigenvalue signature and value ratios, we may classify at each scale whether we are dealing with a blob, filament or a sheet. Their eigenvalue signatures are listed in table I.
Following the sequence clusters filaments walls, we identify the regions and scales at which the local matter distribution follows the corresponding eigenvalue signature. For each of these morphological components  walls, filaments and blobs  the MMF defines a morphology response filter whose form is dictated by the particular morphological feature it seeks to extract and whose local value depends on the local shape and spatial coherence of the density field. Detailed expressions for the morphology response filter can be found in [36]. The morphology signal at each location is then defined to be the one with the maximum response across the full range of smoothing scales. By means of a density criterion the physically significant clumps amongst the blobs are found, while a percolation criterion is invoked to select the physically significant filaments.
Figure 8 shows the morphological segmentation of a 150 Mpc CDM simulation obtained with the MMF. The figure illustrates the large scale distribution of matter as an interconnected network of filaments (dark gray) defining the boundaries of walls (light gray), and the blobs (black) located at the intersections of the network. For clarity fig. 8 only show the largest structures. By restricting the number of structures included it is easier to identify the individual components of the Cosmic Web. Each morphological component is well differentiated and occupies, by construction, mutually exclusive regions. The ability of the MMF to identify the features of the matter distribution at different scales may be inferred from the variety of sizes of blobs.
Focussing on the filamentary component of the Cosmic Web, fig. 9 shows a slice of the simulation box in which MMF identified filaments have been compressed to delineate their spines (see AragónCalvo [35] for a description of the compression algorithm). Gray circles indicate the location of clusters with masses above M h. The figure presents the cosmic web as a network of interconnected filaments spread all over the simulation box. The clusters sit at the intersections or ”nodes” of the network [1]. The filamentary network permeates all regions of space, even the very underdense voids.
VD Assessment of MMF
Figures 8 and 9 testify of the potential of the MMF method for studying the structure of the Cosmic Web. Despite its virtues, we may identify various issues. One issue concerns the use of artificial, spherically symmetric filters for constructing Scale Space. In the light of the evident anisotropic nature of the cosmic matter distribution, this remains a questionable choice. Also, the choice of filter radii is rather artificial, and may be improved by finding filter radii more attuned to the particular spectrum of matter flucuations. The singling out of individual morphological components prevents a direct and natural connection of these elements in a pervasive structure. Finally, in the light of the Cosmic Web theory, with its stress on the key role of tidal field, one may question whether an analysis on the basis of the Hessian of the gravitational potential field would not make more sense.
Partly, these issues are alleviated within the context of the Cosmic Spine method (see sect. VII).
Vi the Watershed Void Finder
Voids are enormous regions with sizes in the range of Mpc that are practically devoid of any galaxy and usually roundish in shape. Forming an essential ingredient of the Cosmic Web [1], they are surrounded by elongated filaments, sheetlike walls and dense compact clusters. Voids are distinctive and striking features of the cosmic web (see fig. 10), yet identifying them and tracing their outline within the complex geometry of the galaxy and mass distribution in galaxy surveys and simulations has proven to be a nontrivial issue.
There have been extensive searches for voids in galaxy catalogues and in numerical simulations. The fact that voids are almost empty of galaxies means that the sampling density plays a key role in determining what is or is not a void. There is not an unequivocal definition of what a void is and as a result there is considerable disagreement on the precise outline of such a region (see e.g. [48]). Moreover, voidfinders are often predicated on building void structures out of cubic cells [49] or out of spheres (e.g. [50, 51]). Such methods attempt to synthesize voids from the intersection of cubic or spherical elements and do so with varying degrees of success. Because of the vague and different definitions, and the range of different interests in voids, there is a plethora of void identification procedures. The “voidfinder” algorithm of [52] has been at the basis of most voidfinding methods. However, this succesfull approach will not be able to analyze complex spatial configurations in which voids may have arbitrary shapes and contain a range and variety of substructures.
The watershedbased WVF algorithm of Platen et al. [37] aims to avoid issues of both sampling density and shape. This new and objective voidfinding formalism has been specifically designed to dissect in a selfconsistent manner the multiscale character of the void network and the weblike features marking its boundaries. The WVF is defined with respect to the DTFE density field, assuring optimal sensitivity to the morphology of spatial structures and an unbiased probe over the full range of substructure in the mass distribution.
Via Watershed Transform
The Watershed Void Finder (WVF) is based on the Watershed Transform [53, 54, 55], a commonly used method in Image Analysis. It is a concept defined within the context of mathematical morphology, and was first introduced by Beucher & Lantuejoul [53]. It is widely used for segmenting images into distinct patches and features.
The basic idea behind the WST stems from geophysics, where the word watershed found its origin in the analogy of the procedure with that of a landscape being flooded by a rising level of water. The WST is used to delineate the boundaries of separate domains, i.e. basins into which yields of e.g. rainfall will collect. Its name may be appreciated from following the analogy of a surface in the shape of a landscape. After the surface is pierced at the location of each of the minima, and as the waterlevel rises, a growing fraction of the landscape will be flooded by the water in the expanding basins. Ultimately basins will meet at the ridges corresponding to saddlepoints in the density field.
In many practical situations, including that of the cosmological density fields, we are dealing with a landscape of a smooth Cdifferentiable density field. The mathematical framework for the topological analysis of such surfaces is Morse theory. In the case of a proper Morse function, all critical points are nondegenerate with distinct values. Based on the location and nature of the singular points – minima, maxima or saddle points – and the corresponding integral lines or slope lines in the density field (see e.g. fig. 11),
(18) 
one may infer a variety of spatial segmentations (see e.g. [56, 57, 58, 59, 60, 61]). This may involve regions connected to the maxima, regions connected to the minima, or regions connecting a maximum to a particular minimum via the integral lines of the field. The latter is the MorseSmale complex of the manifold (also see sect. VII on the Cosmic Spine). The slope lines that define the boundaries of adjacent valleys are called watersheds, while watershed lines are the set of slope lines emanating from saddle points and connecting to a local maximum or minimum.
The segmentation of the image is one in regions of influence, defined by a minimum and those points whose topographic distance to any of the minima in the landscape is minimal,
(19) 
In this definition, the integral denotes the field pathlength along all paths in the set of all possible paths, . This concept of distance is related to the geodesics of the surface : the path of steepest descent, which specifies the track a droplet of water would follow as it flows down a mountain surface. An illustrative example of a landscape and its segmentation into watershed basins surrounded by watershed lines can be seen in fig. 10.
ViB Cosmological Watersheds
Extrapolating the application of the watershed transform to other areas of interest, its implementation also represents a practical instrument for the segmentation of surfaces and volumes defined by the topological structure of cosmological density fields. When studying the topological and morphological structure of the cosmic matter distribution in the Cosmic Web, it is convenient to draw the analogy with a landscape (see fig. 10). Valleys represent the large underdense voids that define the cells of the Cosmic Web. Their boundaries are sheets and ridges, defining the network of walls, filaments and clusters that defines the Cosmic Web.
The WVF focusses specifically on the valleys, which the Watershed Transform identifies with the cosmic voids in the density field. The watershed algorithm holds several advantages with respect to other voidfinders:

Within an ideal smooth density field (i.e. without noise) it will identify voids in a parameter free way. No predefined values have to be introduced. In less ideal, and realistic, circumstances a few parameters have to be set for filtering out discreteness noise. Their values are guided by the properties of the data.

The watershed works directly on the topology of the field and does not reply on a predefined geometry/shape. By implication the identified voids may have any shape.

The watershed naturally places the divide lines on the crests of a field. The void boundary will be detected even when its boundary is distorted.

The transform naturally produces closed contours. As long as minima are well chosen the watershed transform will not be sensitive to local protrusions between two adjacent voids.
ViC Discrete Watershed Transform
There are many different implementations of the watershed transform and the development of such algorithms is still an active field in computer image analysis. An important characteristic of the nearly all modern scientific images is their discrete nature. This concerns their spatial discreteness, sampled at discrete intervals, and their discrete intensity levels. While their discreteness allows the use of highly efficient algorithms, it does also involve complications due to the intrinsic difficulty in identifying saddle points from a discrete local neighborhood and of accurately extracting slope lines. On the other hand, by removing faint features the discretization also helps to remove some artefacts of the watershed transform without the need of pre or postprocessing.
Several methods have been devised in an attempt to alleviate the limitations imposed by discrete images for the extraction of critical points. Among them, the discrete watershed transform algorithm [62] represents a simple and elegant way of identifying the watershed separatrices and it can be shown to converge to the continuous case. The watershed transform emulates the flooding of valleys. As these are flooded the newly formed “lakes” increase their area until they touch each other. The points where two or more lakes converge are marked and the algorithm continues until all the pixels in the image have been flooded. At the end of the flooding process the image will be segmented into individual regions or catchment basins sharing a local minima. The resulting watershed transform is defined by the set of all points marked as the dividing boundaries between two or more basins.
ViD the WVF algorithm
We have implemented a version of the discrete watershed algorithm of Beucher [62], and incorporated it in the WVF void identification procedure. The WVF consists of the following steps:

DTFE: Given a point distribution (Nbody, redshift survey), the Delaunay Tessellation Field Estimator is used to define a continuous density field throughout the sample volume.

Grid Sampling: For practical processing purposes the DTFE field is sampled on a grid. The optimal grid size has to assure the resolution of all morphological structures while minimizing the number of needed gridcells.

RankOrdered Filtering: The DTFE density field is adaptively smoothed by means of Natural Neighbour Maxmin and Median filtering. This involves the computation of the median, minimum or maximum of densities within the contiguous Voronoi cell, the region defined by a point and its natural neighbours .

Contour Levels: The image is transformed into a discrete set of density levels. The levels are defined by a uniform partitioning of the cumulative density distribution.

Pixel Noise: With a morphological opening and closing of 2 pixel radius we further reduce pixel by pixel fluctuations.

Field Minima: The minima in the smoothed density field are identified as the pixels (grid cells) which are exclusively surrounded by neighbouring gridcells with a higher density value.

Flooding: The flooding procedure starts at the location of the minima. At successively increasing flood levels the surrounding region with a density lower than the corresponding density threshold is added to the basin of a particular minimum.

Segmentation: Once a pixel is reached by two distinct basins it is identified as belonging to their segmentation boundary. By continuing this procedure up to the maximum density level the whole region is segmented into distinct void patches.

Hierarchy Correction: A correction is necessary to deal with effects related to the intrinsic hierarchical nature of the void distribution. The correction involves the removal of segmentation boundaries whose density is lower than some density threshold.
ViE WVF by example: Voids in SDSS
A telling impression of the practical implementation of the watershed voidfinding method may be obtained from fig. 12. The watershed procedure has been applied to the SDSS DR7 galaxy survey.
The Sloan Digital Sky Survey is the largest campaign for mapping the spatial distribution of galaxies in the local universe. It has imaged 25% of the sky: with the latest DR7 data release it completed the survey of objects over 8423 sq. deg. area on the sky. Around 230 million objects were imaged in 5 photometric bands, out to apparent magnitude (r band, nm)^{4}^{4}4In astronomy, the brightness of objects is quantified in terms of magnitude. Originally defined by Hipparchus, and in accordance with the logarithmic sensitivity of the eye, it concerns the logarithm of the brightness of the objects, . Note that a lower magnitude implies a brighter object. The human eye can see stars brighter than . The star Vega defined the zeropoint, , while Sirius, the brightest star on the sky (except for the Sun) has . A galaxy with a magnitude is really very faint !. Of 928,567 galaxies, those brighter than , the redshift (distance) was measured.
Fig. 12 shows the DTFE density field inferred from the discrete galaxy distribution in a wide slice through the DR7 sample. It concerns a region which goes out to a redshift , ie. out to a distance of Mpc. The galaxies are superimposed as blue dots. Note the lower spatial resolution of the survey at larger distances: as their distance increases only the brightest galaxies can be seen, resulting in a reduced spatial sampling.
The blue circles indicate the position and size (but not shape, and outline) of the voids in the SDSS. WVF identified voids in the SDSS. The WVF voids are compared with those found by the voidfinder of Hoyle & Vogeley [50], indicated as green circles. The location of the larger voids match quite well, although there are differences concerning the size of the identified voids. A more substantial mismatch exists for smaller voids.
ViF Assessment of WVF
With respect to the other void finders the watershed algorithm has several advantages. Because it identifies a void segment on the basis of the crests in a density field surrounding a density minimum it is able to trace the void boundary even though it has a distorted and twisted shape. Also, because the contours around well chosen minima are by definition closed the transform is not sensitive to local protrusions between two adjacent voids. The main advantage of the WVF is that for an ideally smoothed density field it is able to find voids in an entirely parameter free fashion.
Vii the Cosmic Spine & Morse Theory
The SpineWeb method [39] addresses the topological structure of the landscape, defined by its singularities and their mutual connections, that of minima, maxima and saddle points. It invokes the local adjacency properties of the boundaries between the watershed basins to trace the critical points in the density field and the separatrices defined by them. The separatrices are classified into walls and the spine, the network of filaments and nodes in the matter distribution. The SpineWeb formalism is intimately related to Morse theory, in which it finds a solid mathematical foundation.
In the discussion of the WVF (sec. VI), we have seen that the watershed transform is a key instrument for the segmentation of a density field, and as such is also ideally suited for tracing the boundaries between the identified segments. With this in mind we have invoked the watershed transform in a more extensive framework for studying both identity and topology of the cosmic web and its various constituents. The result is the SpineWeb method, a complete framework for the identification of voids, walls and filaments. Based on the watershed segmentation of the cosmic density field, the method invokes the local properties of the regions adjacent to the critical points, which define its separatrices. An inspection of the flow lines in a 2D density field, and the corresponding separatrices, forms a telling illustration of the idea (fig. 11).
Viia SpineWeb and the Cosmic Web
The SpineWeb technique is based on the following key aspects of the topological structure of the Cosmic Web:

The Cosmic Web is an interconnected system of walls, filaments and clusters. High resolution Nbody simulations indicate that all the elements of the Cosmic Web are interconnected. As mass resolution increases, one finds structures that form the connection between otherwise apparently isolated objects.

Walls are the defining boundaries between adjacent voids. The analogy between the watershed transform defining the boundary between underdense basins and walls defining the boundaries between voids is one of the major justifications of the method presented here.

Filaments run between clusters and across the intersection of walls. Membranes permeate the space between adjacent filaments.
An important aspect of the SpineWeb formalism is the practical role of the watershed transform in computing the Morse complex of a density field. In sect. VI we have shortly described how the flow field of a landscape can be used towards defining a variety of spatial segmentations (see e.g. [56, 57, 58, 59, 60, 63, 61, 64]). These are based on the connection between the singularities  maxima, minima and sadddle points  in the landscape by means of the slope lines. Within this context, the watershed lines are the set of slope lines emanating from saddle points and connecting to a local maximum or minimum, with the ones connecting to maxima defining the watershed boundaries around the valleys.
The Morse theoretical basis of the SpineWeb formalism relates it to the skeleton analysis of the Cosmic Web [65, 66, 67], which is also based on Morse theory [68, 69]. It formed the basis for the development of an elegant and mathematically rigorous tool for filament identification, which in the meantime has also been extended towards tracing a range of morphological features [70]. Its present implementation refers to only one specific scale. The SpineWeb procedure is being generalized to a fully multiscale formalsim.
ViiB SpineWeb procedure
The SpineWeb procedure is based on an implementation of the discrete watershed transform code [62, 71]. The code starts by identifying and labeling all the local minima present in the density field as the pixels which have the lowest value among its 26 neighbours. Subsequently, for each voxel the neighbour with lowest density is determined. By pursuing this, the procedure iteratively traces the maximum gradient paths until they converge to their local minima. This topographical distance algorithm yields a fast segmentation of space into locally connected underdense regions.
In the second step we identify the voxels embedded in the watershed transform. For this, a local immersion algorithm is used. All the voxels located at the boundaries between two or more regions are identified, Subsequently, this subset of voxels is sorted in density, followed by the application of the standard immersion algorithm. The advantage of this twostep process is that there is no need for sorting the complete density field. Instead, only the voxels in the boundaries around the watershed segments are addressed.
In the final step, each boundary voxel is assigned a morphological label as void, wall or filament on the basis of their local neighbourhood, ie. the 26 adjacent voxels to each boundary element. The assigned morphology follows a set of simple but effective criteria based on the number of adjacent voids/basins to a given boundary voxel (see fig. 13),
(20) 
It is important to emphasize that these criteria for identifying filaments and walls from the watershed transform relies only on the topology of the density field.
ViiC Filaments, Sheets and the Spine
The identification of walls with locations marked by two adjacent voids is straightforward: they are sheets that collect matter from their surrounding voids. The identification of filaments depends strongly on the resolved cellular structure of the Cosmic Web. Filaments that are related to voids whose outline is visible in the density map are succesfully traced. However, when the density field representation is too coarse, or when the substructure of voids is insufficiently resolved, the procedure tends to produce a variety of less pronounced and smaller filamentary features that have not been identified as spinal filaments. In the hierarchical multiscale extension of the SpineWeb method this is directly addressed.
Figure 14 shows the full 3D network of filaments (red) and walls (blue) in a cosmological Nbody simulation. It concerns a CDM universe inside a box of 200 Mpc, restricted to the dark matter particles. Initial conditions were generated on a grid with , , and . Lowerresolution versions of , and particles were generated from the same initial conditions. The depicted spine structure focuses on the largest filaments and walls in the particle distribution, following from the analysis of the lowresolution dataset. This lower resolution corresponds to a cutoff scale of roughly hMpc.
For visualization purposes, the walls and filaments are separately smoothed with a Gaussian kernel of voxels. Figure 14 clearly shows the threedimensional nature of the filamentwall network and the close relation between both: filaments define an interconnected web and walls “fill” the spaces in between filaments forming a watertight network of voids. Walls can be seen as extended membranes with somewhat irregular surfaces on small scales. On larger scales walls can be considered semiplanar structures, although it is possible to find examples of walls with more convoluted shapes reflecting equally complex local environments. Filaments (here an individual filament is defined as the filament between two intersections points) also present small scale features but overall they can be considered semilinear, even though one can find very twisted structures. Short filaments appear to be more straight than long ones.
ViiD the Hierarchical Spine
An important aspect of our SpineWeb method is that it is an intrinsically scalefree method, starting from a scalefree reconstruction of the density field. The DTFE density field guarantees an optimal and unbiased representation of the hierarchical nature and anisotropic morphology of cosmic structure. Having guaranteed the capability of invoking a full scalefree Scalespace representation of cosmic structure, the watershed procedure not only traces the outline of filaments and sheets, but may also be extended towards doing so over a range of scales in order to address their hierarchical structure.
At present, the spine formalism has indeed been extended to a genuine multiscale formalism. This allows the identification of the spine over a range of scales. The spine segments on a range of scales are subsequently related to each other, allowing the construction of hierarchical filament, wall and void trees (see upcoming publication [72]).
The hierarchical spine procedure has demonstrated that seemingly unrelated substructures at one scale, do actually show up in the spine of higher resolution realizations of the density field. This is a direct reflection of the presence of the cosmic web’s topological outline in the primordial density field, and assures the fact that filaments always lie at the intersection of walls and form the edges of empty cells.
ViiE Tessellation Spines
While in the present implementation we determine the spine via the DTFE density field sampled on a regular grid, a highly promising innovation will be to circumvent the density field reconstruction and gridbased spineweb identification procedure entirely by working directly on the Delaunay tessellation defined by the point sample itself.
From the DTFE discussion in sect. II we have learnt that the Delaunay/Voronoi tessellation contains all necessary information on the local density and geometry in the sample. Instead of first determining the density first, we could define local criteria for the Delaunay tetrahedra for tracing the singularities and the separatrices in the corresponding field. It will render the Spineweb procedure independent of an artificial gridscale, grid geometry and save the consumption of CPU time for computing the DTFE field. The extraction of the complete hierarchy of weblike structures could be accomplished in one calculation, that of the tessellation.
Various visualization and image processing studies have been seeking to define Morse theory and watershed segmentation algorithms that work directly via the Delaunay tessellation of the discrete dataset [60, 59, 73, 74, 75, 64].
The high potential of this strategy is illustrated by fig. 15. It involves the a first test implementation of the SpineWeb procedure directly on the Delaunay grid [76]. Instead of resorting to the 26 neighbours on a regular grid, the morphological evaluations (eqn.20) for each point in the sample concerns its natural neighbours, the points with whom it shares a Delaunay tetrahedron. The two panels of fig. 15 shows the points in the filaments and sheets of the Cosmic Web, along with their Delaunay (edge) connections. The intricate weblike features are clearly represented in unprecedented detail, both the rather tenuous thin membranes and the dense and abundantly branching filaments.
Viii Topology & Alphashapes
A final tessellationbased technique for the analysis of aspects of the Cosmic Web relates to its topological aspects. Alpha shape is a description of the (intuitive) notion of the shape of a discrete point set. Alphashapes of a discrete point distribution are subsets of a Delaunay triangulation and were introduced by Edelsbrunner and collaborators [77, 78, 79, 80]. Alpha shapes are generalizations of the convex hull of a point set and are concrete geometric objects which are uniquely defined for a particular point set. Reflecting the topological structure of a point distribution, it is one of the most essential concepts in the field of Computational Topology [81, 82, 83].
If we have a point set and its corresponding Delaunay triangulation, we may identify all Delaunay simplices – tetrahedra, triangles, edges, vertices – of the triangulation. For a given nonnegative value of , the alpha complex of a point sets consists of all simplices in the Delaunay triangulation which have an empty circumsphere with squared radius less than or equal to , . Here “empty” means that the open sphere does not include any points of . For an extreme value the alpha complex merely consists of the vertices of the point set. The set also defines a maximum value , such that for the alpha shape is the convex hull of the point set.
The alpha shape is the union of all simplices of the alpha complex. Note that it implies that although the alpha shape is defined for all there are only a finited number of different alpha shapes for any one point set. The alpha shape is a polytope in a fairly general sense, it can be concave and even disconnected. Its components can be threedimensional patches of tetrahedra, twodimensional ones of triangles, onedimensional strings of edges and even single points. The set of all real numbers leads to a family of shapes capturing the intuitive notion of the overall versus fine shape of a point set. Starting from the convex hull of a point set and gradually decreasing the shape of the point set gradually shrinks and starts to develop cavities. These cavities may join to form tunnels and voids. For sufficiently small the alpha shape is empty.
Following this description, one may find that alpha shapes are intimately related to the topology of a point set. As a result they form a direct and unique way of characterizing the topology of a point distribution. A complete quantitative description of the topology is that in terms of Betti numbers and these may indeed be directly inferred from the alpha shape. The first Betti number specifies the number of independent components of an object. In this context may be interpreted as the number of independent tunnels, and as the number of independent enclose voids. The Betti number effectively counts the number of independent dimensional holes in the simplicial complex.
Applications of alpha shapes have as yet focussed on biological systems. Their use in characterizing the topology and structure of macromolecules. The work by Liang and collaborators [84, 85, 86, 87] uses alpha shapes and betti numbers to assess the voids and pockets in an effort to classify complex protein structures, a highly challenging task given the 10,00030,000 protein families involving 1,0004,000 complicated folds. Given the interest in the topology of the cosmic mass distribution [88, 89, 90], it is evident that alpha shapes also provide a highly interesting tool for studying the topology of the galaxy distribution and Nbody simulations of cosmic structure formation. Directly connected to the topology of the point distribution itself it would discard the need of userdefined filter kernels.
In a recent study Vegter et al. [82] computed the alpha shapes for a set of GIF simulations of cosmic structure formation (see fig. 16). On the basis of a calibration of the inferred Minkowski functionals and Betti numbers from a range of Voronoi clustering models their study intends to refine the knowledge of the topological structure of the Cosmic Web.
Acknowledgement
We are very grateful to Willem Schaap for his contributions to substantial parts of the work presented in this review. We also thank Gert Vegter for very useful discussions on various aspects of the presented work, and Danny Pan for permission to use fig. 12. Finally, we wish to express our gratitude to the organizers of the ISVD09 conference, foremost F. Anton, for the invitation to present this review.
References
 [1] J. Bond, L. Kofman, and D. Pogosyan, “How filaments are woven into the cosmic web,” Nature, vol. 380, pp. 603–606, 1996.
 [2] M. Colless and 2dF consortium, “The 2df galaxy redshift survey: Final data release,” pp. 1–32, 2003, astroph/0306581.
 [3] M. Tegmark et al., “The threedimensional power spectrum of galaxies from the sloan digital sky survey,” Astrophys. J., vol. 606, pp. 702–740, 2004.
 [4] J. Huchra et al., “The 2mass redshift survey and low galactic latitude largescale structure,” in Nearby LargeScale Structures and the Zone of Avoidance, ser. ASP Conf. Ser. Vol., A. Fairall and P. Woudt, Eds., vol. 239. Astron. Soc. Pacific, San Francisco, 2005, pp. 135–146.
 [5] P. Peebles, The Large Scale Structure of the Universe. Princeton Univ. Press, 1980.
 [6] V. Springel et al., “Simulations of the formation, evolution and clustering of galaxies and quasars,” Nature, vol. 435, pp. 629–636, 2005.
 [7] J. Bardeen, J. Bond, N. Kaiser, and A. Szalay, “The statistics of peaks of gaussian random fields,” Astrophys. J., vol. 304, pp. 15–61, 1986.
 [8] R. Sheth and R. van de Weygaert, “Much ado about nothing: a hierarchy of voids,” Mon. Not. R. Astron. Soc., vol. 350, pp. 517–538, 2004.
 [9] R. van de Weygaert and E. Bertschinger, “Peak and gravity constraints in gaussian primordial density fields: An application of the hoffmanribak method,” Mon. Not. R. Astron. Soc., vol. 281, pp. 84–118, 1996.
 [10] R. van de Weygaert and J. R. Bond, “Clusters and the theory of the cosmic web,” in A PanChromatic View of Clusters of Galaxies and the LargeScale Structure, ser. Lecture Notes in Physics, M. Plionis, O. LopezCruz, and D. Hughes, Eds. SpringerVerlag, 2008, no. 740, pp. 335–407.
 [11] D. Y. Pogosyan, J. Bond, L. Kofman, and J. Wadsley, “Cosmic web: Origin and observables,” in Wide Field Surveys in Cosmology (14th IAP meeting), S. Colombi and Y. Mellier, Eds. Editions Frontieres, 1998, pp. 61–67.
 [12] W. Schaap and R. van de Weygaert, “Continuous fields and discrete samples: reconstruction through delaunay tessellations,” Astron. Astrophys., vol. 363, pp. L29–L32, 2000.
 [13] F. Bernardeau and R. van de Weygaert, “A new method for accurate estimation of velocity field statistics,” MNRAS, vol. 279, pp. 693–711, 1996.
 [14] A. Okabe, B. Boots, K. Sugihara, and S. N. Chiu, Spatial Tessellations, Concepts and Applications of Voronoi Diagrams. Chichester: John Wiley, 2000.
 [15] R. Sibson, “A vector identity for the dirichlet tessellation,” Math. Proc. Cambridge Phil. Soc., vol. 87, pp. 151–155, 1980.
 [16] ——, “A brief description of natural neighbor interpolation,” in Interpreting Multivariate Data, V. Barnett, Ed. Wiley, Chichester, 1981, pp. 21–36.
 [17] D. Watson, Contouring: A Guide to the Analysis and Display of Spatial Data. Pergamom Press, Oxford, 1992.
 [18] J. Braun and M. Sambridge, “A numerical method for solving partial differential equations on highly irregular evolving grids,” Nature, vol. 376, pp. 655–660, 1995.
 [19] N. Sukumar, The Natural Element Method in Solid Mechanics, 1998, phD thesis, Northwestern University.
 [20] W. E. Schaap, the Delaunay Tessellation Field Estimator, 2007, phD thesis, Groningen University.
 [21] R. van de Weygaert and W. E. Schaap, “The cosmic web: Geometric analysis,” in Data Analysis in Cosmology, ser. Lecture Notes in Physics, V. Martínez, E. Saar, E. MartínezGonzalez, and M. PonsBordería, Eds. SpringerVerlag, 2008, no. 665, pp. 289–419.
 [22] C. Lawson, “Software for surface interpolation,” in Mathematical Software III, J. Rice, Ed. Academic Press, New York, 1977, vol. 3.
 [23] S. Sloan, “A fast algorithm for constructing delaunay triangulations in the plane,” Adv. Eng. Software, vol. 9, pp. 34–55, 1987.
 [24] M. Sambridge, J. Braun, and H. McQueen, “Geophysical parameterization and interpolation of irregular data using natural neighbours,” Geophys. J. Int., vol. 122, pp. 837–857, 1995.
 [25] G. Brown, “Point density in stems per acre,” New Zealand Forestry Service Research Notes, vol. 38, pp. 1–11, 1965.
 [26] J. Ord, “How many trees in a forest,” Mathematical Scientist, vol. 3, pp. 23–33, 1978.
 [27] H. Ebeling and G. Wiedenmann, “Detecting structure in two dimensions combining voronoi tessellation and percolation,” Phys. Rev. E, vol. 47, pp. 704–710, 1993.
 [28] M. Ramella, W. Boschin, D. Fadda, and M. Nonino, “Finding galaxy clusters using voronoi tessellations,” Astron. Astrophys., vol. 368, pp. 776–786, 2001.
 [29] R. S. J. Kim et al., “Detecting clusters of galaxies in the sloan digital sky survey. i. monte carlo comparison of cluster detection algorithms,” Astron. J., vol. 123, pp. 20–36, 2002.
 [30] C. Marinoni, M. Davis, J. Newman, and A.L.Coil, “Threedimensional identification and reconstruction of galaxy systems with fluxlimited redshift surveys,” Astrophys. J., vol. 580, pp. 122–143, 2002.
 [31] P. Lopes, R. de Carvalho, R. Gal, S. Djorgovski, S. Odewahn, A. Mahabal, and R. Brunner, “The northern sky optical cluster survey. iv. an intermediateredshift galaxy cluster catalog and the comparison of two detection algorithms,” Astron. J., vol. 128, pp. 1017–1045, 2004.
 [32] M. C. Neyrinck, N. Y. Gnedin, and A. J. S. Hamilton, “Voboz: an almostparameterfree halofinding algorithm,” Mon. Not. R. Astron. Soc., vol. 356, pp. 1222–1232, 2005.
 [33] Y. Ascasibar and J. Binney, “Numerical estimation of densities,” MNRAS, vol. 356, pp. 872–882, 2005.
 [34] V. M. V., J.L. Starck, E. Saar, D. Donoho, S. R. S.C., P. de la Cruz, and S. Paredes, “Morphology of the galaxy distribution from wavelet denoising,” Astrophys. J., vol. 634, pp. 744–755, 2005.
 [35] M. AragónCalvo, Morphology and Dynamics of the Cosmic Web, 2007, phD thesis, Groningen University.
 [36] M. AragónCalvo, B. Jones, R. van de Weygaert, and J. van der Hulst, “The multiscale morphology filter: identifying and extracting spatial patterns in the galaxy distribution,” Astron. Astrophys., vol. 474, pp. 315–338, 2007.
 [37] E. Platen, R. van de Weygaert, and B. Jones, “A cosmic watershed: the wvf void detection technique,” Mon. Not. R. Astron. Soc., vol. 380, pp. 551–570, 2007.
 [38] M. C. Neyrinck, “Zobov: a parameterfree voidfinding algorithm,” Mon. Not. R. Astron. Soc., vol. 386, pp. 2101–2109, 2008.
 [39] M. AragónCalvo, E. Platen, R. van de Weygaert, and A. Szalay, “the spine of the cosmic web,” Mon. Not. R. Astron. Soc., p. 16 pp., 2009, subm.
 [40] M. Bradač, P. Schneider, M. Lombardi, M. Steinmetz, L. Koopmans, and J. Navarro, “The signature of substructure on gravitational lensing in the cdm cosmological model,” Astron. Astrophys., vol. 423, pp. 797–809, 2004.
 [41] E. RomanoDíaz and R. van de Weygaert, “Delaunay tessellation field estimator analysis of the pscz local universe: density field and cosmic flow,” Mon. Not. R. Astron. Soc., vol. 382, pp. 2–28, 2007.
 [42] W. E. Schaap and R. van de Weygaert, “Dtfe tracing the cosmic web: Hierarchical substructure,” Astron & Astrophys., p. 20 pp., 2009, subm.
 [43] ——, “Dtfe tracing the cosmic web: Complex geometric patterns,” Astron & Astrophys., p. 15 pp., 2009, subm.
 [44] E. Platen, M.A. AragónCalvo, R. van de Weygaert, B. Jones, and G. Vegter, “the sdss density field reconstruction, the methods,” Mon. Not. R. Astron. Soc., p. 23 pp., 2009, to be subm.
 [45] L. Florack, B. ter Haar Romeny, J. Koenderink, and M. Viergever, “Scale and the differential structure of images,” Image and Vision Computing, vol. 10, pp. 376–388, 1992.
 [46] A. Frangi, W. Niessen, K. Vincken, and M. Viergever, “Multiscale vessel enhancement filtering,” in Medical Image Computing and ComputerAssisted Interventation  MICCAI’98, ser. Lecture Notes in Computer Science. SpringerVerlag, 1998, no. 1496, pp. 130–137.
 [47] Y. Sato, D. Nakajima, H. Atsumi, T. Koller, G. Gerig, S. Yoshida, and R. Kikinis, “3d multiscale line filter for segmentation and visualization of curvilinear structures in medical images,” IEEE Medical Image Analysis, vol. 2, pp. 143–168, 1998.
 [48] S. Shandarin, H. Feldman, K. Heitmann, and S. Habib, “Shapes and sizes of voids in the lambda cold dark matter universe: excursion set approach,” Mon. Not. R. Astron. Soc., vol. 376, pp. 1629–1640, 2006.
 [49] G. Kauffmann and A. Fairall, “Voids in the distribution of galaxies  an assessment of their significance and derivation of a void spectrum,” Mon. Not. R. Astron. Soc., vol. 248, pp. 313–324, 1991.
 [50] F. Hoyle and M. Vogeley, “Voids in the point source catalogue survey and the updated zwicky catalog,” Astrophys. J., vol. 566, pp. 641–651, 2002.
 [51] S. Patiri, J. BetancortRijo, F. Prada, A. Klypin, and S. Gottloeber, “Statistics of voids in the twodegree field galaxy redshift survey,” Mon. Not. R. Astron. Soc., vol. 369, pp. 335–348, 2006.
 [52] H. ElAd and T. Piran, “Voids in the largescale structure,” Astrophys. J., vol. 491, pp. 421–435, 1997.
 [53] S. Beucher and C. Lantuejoul, “Use of watersheds in contour detection,” in International Workshop on Image Processing: Realtime Edge and Motion Detection/Estimation. CCETT/IRISA, 1979, pp. 17–21.
 [54] F. Meyer and S. Beucher, “Morphological segmentations,” J. Visual Comm. Image Rep., vol. 1, pp. 21–46, 1990.
 [55] S. Beucher and F. Meyer, “The morphological approach to segmentation: the watershed transformation,” in Mathematical Morphology in Image Processing, Marcel Dekker, E. Dougherty, Ed. New York, 1993, ch. 12, pp. 433–481.
 [56] A. Cayley, “On contour and slope lines,” The London, Edinburgh and Dublin Philosophical Magazine and Journal of Science, vol. 18, pp. 264–268, 1859.
 [57] J. Maxwell, “On hills and dales,” The London, Edinburgh and Dublin Philosophical Magazine and Journal of Science, vol. 40, pp. 421–425, 1870.
 [58] H. Edelsbrunner, J. Harer, and A. Zomorodian, “Hierarchical morsesmale complexes for piecewise linear 2manifolds,” Discrete and Computational Geometry, vol. 30, pp. 87–107, 2003.
 [59] H. Edelsbrunner, J. Harer, V. Natarajan, and V. Pascucci, “Morsesmale complexes for piecewise linear 3manifolds,” in Proc. 19th ann. symp. Computational geometry, V. Barnett, Ed., 2003, pp. 361–370.
 [60] E. Danovaro, L. D. Floriani, P. Magillo, M. M. Mohammed, and E. P. E., “Morphologydriven simplification and multiresolution modelling of terrains,” in Geographic Information Systems: Proc. 11th ACM intern. symp. Adv. geographic information systems, 2003, pp. 63–70.
 [61] A. Gyulassy, V. Natarajan, V. Pascucci, P. Bremer, and B. Hamann, “Topologybased simplification for feature extraction from 3d scalar fields,” in Proc. IEEE Conf. Visualization, 2005, pp. 535–543.
 [62] S. Beucher, “Watersheds of functions and picture segmentation,” in ICASSP82, 1982, pp. 1928–1931.
 [63] A. Gyulassy, P.T. Bremer, V. Pascucci, and B. Hamann, “Hierarchical morsesmale complex in 3d,” 2004, manuscript.
 [64] A. Gyulassy, V. Natarajan, V. Pascucci, and B. Hamann, “Efficient computation of morsesmale complexes for threedimensional scalar functions,” IEEE Transactions on Visualization and Computer Graphics, vol. 13, no. 6, pp. 1440–1447, 2007.
 [65] D. Novikov, S. Colombi, and O. Doré, “Skeleton as a probe of the cosmic web: the twodimensional case,” Mon. Not. R. Astron. Soc., vol. 366, pp. 1201–1216, 2006.
 [66] T. Sousbie, C. Pichon, S. Colombi, D. Novikov, and D. Pogosyan, “the 3d skeleton: tracing the filamentary structure of the universe,” Mon. Not. R. Astron. Soc., vol. 383, pp. 1655–1670, 2008.
 [67] T. Sousbie, C. Pichon, H. Courtois, S. Colombi, and D. Novikov, “The threedimensional skeleton of the sdss,” Astrophys. J., vol. 672, pp. 1–4, 2008.
 [68] S. Colombi, D. Pogosyan, and T. Souradeep, “Tree structure of a percolating universe,” Phys. Rev. Lett., vol. 85, pp. 5515–5519, 2000.
 [69] D. Pogosyan, C. Pichon, C. Gay, S. Prunet, J. Cardoso, T. Sousbie, and S. Colombi, “The local theory of the cosmic skeleton,” Mon. Not. R. Astron. Soc., pp. 673–706, 2009.
 [70] T. Sousbie, S. Colombi, and C. Pichon, “The fully connected ndimensional skeleton: probing the evolution of the cosmic web,” Mon. Not. R. Astron. Soc., vol. 393, pp. 457–477, 2009.
 [71] J. B. T. M. Roerdink and A. Meijster, “The watershed transform: Definitions, algorithms and parallelization strategies,” Fundamenta Informaticae, vol. 41, pp. 187–228, 2000.
 [72] M. AragónCalvo, E. Platen, R. van de Weygaert, and A. Szalay, “the hierarchical structure of the cosmic spine,” Mon. Not. R. Astron. Soc., 2009, in prep.
 [73] P.T. Bremer, B. Hamann, H. Edelsbrunner, and V. Pascucci, “A topological hierarchy for functions on triangulated surfaces,” IEEE Transactions on Visualization and Computer Graphics, vol. 10, pp. 385–396, 2004.
 [74] P. Magillo, E. Danovaro, L. D. Floriani, L. Papaleo, and M. Vitali, “Extracting terrain morphology: A new algorithm and a comparative evaluation,” in GRAPP (GM/R), 2007, pp. 13–20.
 [75] E. Danovaro, L. D. Floriani, M. Vitali, and P. Magillo, “Multiscale dual morse complexes for representing terrain morphology,” in GIS 2007: Proc. 15th Intern. Symposium Adv. Geographic Inform. Systems ACM. ACM, 2007, pp. 1–8.
 [76] M. AragónCalvo, E. Platen, R. van de Weygaert, and A. Szalay, “Tessellating to the cosmic spine,” Mon. Not. R. Astron. Soc., 2009, in prep.
 [77] H. Edelsbrunner, D. Kirkpatrick, and R. Seidel, “On the shape of a set of points in the plane,” IEEE Trans. Inform. Theory, vol. 29, pp. 551–559, 1983.
 [78] E. Muecke, Shapes and Implementations in threedimensional geometry, 1993, phD thesis, Univ. Illinois UrbanaChampaign.
 [79] H. Edelsbrunner and E. Muecke, “Threedimensional alpha shapes,” ACM Trans. Graphics, vol. 13, pp. 43–72, 1994.
 [80] H. Edelsbrunner, D. Letscher, and A. Zomorodian, “Topological persistance and simplification,” Discrete and Computational Geometry, vol. 28, pp. 511–533, 2002.
 [81] T. Dey, H. Edelsbrunner, and S. Guha, “Computational topology,” in Advances in Discrete and Computational Geometry, B. Chazelle, J. Goodman, and R. Pollack, Eds. American Mathematical Society, 1999, pp. 109–143.
 [82] G. Vegter, R. van de Weygaert, E. Platen, N. Kruithof, and B. Eldering, “Alpha shapes and the topology of cosmic large scale structure,” Mon. Not. R. Astron. Soc., 2009, in prep.
 [83] A. Zomorodian, Topology for Computing, ser. Cambr. Mon. Appl. Comp. Math. Cambr. Univ. Press, 2005.
 [84] H. Edelsbrunner, M. Facello, and J. Liang, “On the definition and the construction of pockets in macromolecules,” Discrete Appl. Math., vol. 88, pp. 83–102, 1998.
 [85] J. Liang, H. Edelsbrunner, P. Fu, P. Sudhakar, and S. Subramaniam, “Analytical shape computation of macromolecules: I. molecular area and volume through alpha shape,” Proteins: Structure, Function, and Genetics, vol. 33, pp. 1–17, 1998.
 [86] ——, “Analytical shape computation of macromolecules: Ii. inaccessible cavities in proteins,” Proteins: Structure, Function, and Genetics, vol. 33, pp. 18–29, 1998.
 [87] J. Liang, C. Woodward, and H. Edelsbrunner, “Anatomy of protein pockets and cavities: measurement of binding site geometry and implications for ligand design,” Protein Science, vol. 7, pp. 1884–1897, 1998.
 [88] J. G. III, M. Dickinson, and A. Melott, “The spongelike topology of largescale structure in the universe,” Astrophys. J., vol. 306, pp. 341–357, 1986.
 [89] K. Mecke, T. Buchert, and H. Wagner, “Robust morphological measures for largescale structure in the universe,” Astron. Astrophys., vol. 288, pp. 697–704, 1994.
 [90] J. Schmalzing, T. Buchert, A. Melott, V. Sahni, B. Sathyaprakash, and S. Shandarin, “Disentangling the cosmic web. i. morphology of isodensity contours,” Astrophys. J., vol. 526, pp. 568–578, 1999.