Structural Analysis of the SDSS Cosmic Web I.Nonlinear Density Field Reconstructions

We investigate the ability of three reconstruction techniques to analyze and investigate weblike features and geometries in a discrete distribution of objects. The three methods are the linear Delaunay Tessellation Field Estimator (DTFE), its higher order equivalent Natural Neighbour Field Estimator (NNFE) and a version of Kriging interpolation adapted to the specific circumstances encountered in galaxy redshift surveys, the Natural Lognormal Kriging technique. DTFE and NNFE are based on the local geometry defined by the Voronoi and Delaunay tessellations of the galaxy distribution. The three reconstruction methods are analysed and compared using mock magnitude-limited and volume-limited SDSS redshift surveys, obtained on the basis of the Millennium simulation. We investigate error trends, biases and the topological structure of the resulting fields, concentrating on the void population identified by the Watershed Void Finder. Environmental effects are addressed by evaluating the density fields on a range of Gaussian filter scales. Comparison with the void population in the original simulation yields the fraction of false void mergers and false void splits. In most tests DTFE, NNFE and Kriging have largely similar density and topology error behaviour. Cosmetically, higher order NNFE and Kriging methods produce more visually appealing reconstructions. Quantitatively, however, DTFE performs better, even while computationally far less demanding. A successful recovery of the void population on small scales appears to be difficult, while the void recovery rate improves significantly on scales>3 h-1Mpc. A study of small scale voids and the void galaxy population should therefore be restricted to the local Universe, out to at most 100 h-1Mpc.


INTRODUCTION
Over the past thirty years a clear paradigm has emerged as large redshift surveys opened the window onto the distribution of matter in our Local Universe: galaxies, intergalactic gas and dark matter exist in a wispy weblike spatial ar-⋆ E-mail: weygaert@astro.rug.nl rangement consisting of dense compact clusters, elongated filaments, and sheetlike walls, amidst large near-empty void regions, with similar patterns existing at higher redshift, albeit over smaller scales. The Cosmic Web is the fundamental spatial organization of matter on scales of a few up to a hundred Megaparsec, scales at which the Universe still resides in a state of moderate dynamical evolution Peebles (1980); Zel'Dovich (1970); Bond et al. (1996). Its appearance has been most dramatically illustrated by the recently produced maps of the nearby cosmos, the 2dF galaxy redshift survey (2dFGRS), the Sloan Digital Sky Survey (SDSS) and the 2MASS redshift surveys Colless et al. (2003); Tegmark et al. (2004); Huchra et al. (2005).
According to the standard lore of structure formation, structures emerged from small perturbations in the primordial field of Gaussian density and velocity perturbations. Under the force of gravity these fluctuations grow and cluster to become the present day observed structures. At large scales the density field has been evolving (quasi)-linearly and still retains much information on cosmological parameters and structure formation. The linear density field provides an abundance of probes from which cosmological information can be extracted. Prominent probes are the clustering of galaxies, which has been used to infer the underlying primordial power spectrum of density fluctuations (Peebles 1980;Percival et al. 2001;Tegmark et al. 2004), the temperature and polarisation anisotropies in the Cosmic Microwave Background (e.g. Smoot et al. 1992;Spergel et al. 2003) and the shearing of galaxy images by the gravitationally lensed photon paths through the inhomogeneous matter distribution (Mellier 1999;Massey et al. 2007;Hoekstra & Jain 2008).
While these (quasi-)linear cosmological probes have yielded an impressive amount of cosmological information, the exploitation of the pronounced nonlinear patterns of the cosmic web towards probing cosmological parameters and cosmic structure formation has been less fortuitous. Even though the morphology, shape and other statistical characteristics of the quasi-linear cosmological density field forms a direct reflection of the structure assembly process in the Universe, on the corresponding small scales nonlinear growth has significantly altered and erased some of the essential cosmological information. The absence of an objective and quantitative procedure for identifying and isolating clusters, filaments and voids in the cosmic matter distribution has been a major obstacle in investigating the structure and dynamics of the Cosmic Web. The overwhelming complexity of the individual structures and their connectivity, the huge range of densities and the intrinsic multi-scale nature prevent the use of simple tools that may be sufficient in less demanding problems. However, various interesting new approaches and methods have been forwarded in the past few years, often based on ideas stemming from image processing, mathematical morphology, and medical imaging (Aragón-Calvo et al. 2007, 2010Sousbie 2011;Way 2011;Genovese et al. 2010, e.g.).
This study is the first in a series in which we systematically investigate the web-related structures found in the Sloan Digital Sky Survey (SDSS). It involves a systematic program in which we explore the cosmography of the Local Universe, assess the statistical characteristics of the density field, identify and categorize the voids and filaments within the SDSS galaxy redshift sample, and study the biasing of the galaxy population with respect to the mass distribution and the dependence of galaxy properties on the large scale environment.
The intention of this paper is the reconstruction of the underlying (nonlinear) density field by translating the spatially irregularly distributed and discrete galaxy sample in the SDSS galaxy redshift survey into a representative den-sity field. The density field reconstructions described in this study will form the basis of a series of studies in which we address the statistical and topological properties of the cosmic density field in our Local Universe and in which we analyze the cosmography, void population and spinal structure of the local Cosmic Web in the Sloan Digital Sky Survey. The cosmographic description of the nearby Universe includes the identification and cataloguing of the filaments, voids and clusters in the SDSS galaxy distribution. The key motivation for the density field reconstruction techniques should therefore be that it allows us to probe the complex quasi-linear (and nonlinear) structures that we find in galaxy redshift surveys. This involves the ability to reproduce the distinct anisotropic -filamentary and wall-like -features that make up the Cosmic Web, as well as the ability to trace the hierarchical substructure of the cosmic matter distribution and the ability to identify the void population that forms one of its most salient features.
We are specifically interested in the performance of the fast and efficient linear tessellation-based DTFE density field reconstruction technique (Schaap & van de Weygaert 2000;Schaap 2007;van de Weygaert & Schaap 2008). In addition, we investigate two additional higher-order techniques which are potentially suited for representing the quasi-linear density field of the cosmic web, the local natural neighbour interpolation (NNFE) and a nonlocal Kriging interpolation technique, the Natural Lognormal Kriging formalism. By means of an extensive comparison we evaluate which aspects of the density field are best reproduced by either DTFE, NNFE or Natural Lognormal Kriging, and which of these methods is best suited to function for further structural analysis.
In order to be able to assess the reliability of the results of our structural tools, it is crucial to understand the details and errors in reconstructed density fields. We will therefore present a detailed comparison study between three different reconstructions methods, and assess their density and topological errors over a range of scales. This will range form the small nonlinear scales, at 1 h −1 Mpc, to a scale of 10 h −1 Mpc, which represents the transition from quasilinear to the linear regime.
We focus on the translation of the galaxy positions in the SDSS DR6 galaxy redshift survey to a representative density field within the survey volume. Currently, SDSS encloses the largest and deepest contiguous region of the nearby Universe mapped by a galaxy redshift survey. This makes the SDSS an ideal data sample for a full threedimensional density field reconstruction. A first impression of the resulting DTFE density field map of a region in the SDSS DR6 survey is shown in figure 1. Our techniques will be generally applicable to any uniform galaxy redshift survey. Even though survey limits and scales are different from that of the SDSS, it will be straightforward to carry over the results to other redshift surveys.

Inferring Cosmic Density Fields
The richest source of data for investigating the intricate weblike cosmic matter distribution is the distribution of galaxies in galaxy redshift surveys. We shall make the assumption that the galaxy distribution is a representative tracer of the underlying mass density field and of the underlying struc- Figure 1. A visualisation of the (DTFE) SDSS density field (1 h −1 Mpc), the contour levels are divided roughly into the overdense and the underdense regime. Both the galaxies (blue dots) and the density represent a slice of 12 h −1 Mpc thickness. Some of the most prominent features have been named, like the Boötes SuperVoid (Kirshner et al. 1981) and the large supervoid (BS SuperVoid) identified by Bahcall & Soneira (1982). Also the largest overdense structure, the (Coma) Great Wall, and the location of the Hercules supercluster within the Wall are indicated.
ture. Different samples selected from such catalogues may of course trace different structures.
We set out to explore high fidelity methods for reconstruction of the density field from the discrete galaxy distribution. We require accurate local density estimates that can be used to reliably compute structural and topological indicators. Computational efficiency is an important factor since our reconstruction technique has to be able to deal with galaxy samples of a million of galaxies as well as with N-body simulations comprising orders of magnitude more particles.
There are two immediate and important implications and complications with respect to our intention of including the (quasi)-nonlinear components in the density field. The first one is that we have to be aware of the noise components that tend to be included in the reconstructions as well (Schaap 2007;van de Weygaert & Schaap 2008). The second one is that nonlinear data are typically not well behaved, marked by strong gradients in the density field. It means that we have to take special care of those locations marked by non-linearities in the data.

Spatial Point Processes and Continuous Density Fields
We pursue the reconstruction of a density field from a spatial point process consisting of irregularly distributed points. We assume that the local intensity of the points is a fair tracer of the density and that these values are samples of a continuous underlying density field. In doing so we appreciate that differently defined samples may trace different structures. The reconstruction problem we address here is a data processing problem. The interpretation of the results is an astronomical process.
When the spatial point process is defined by the galaxy distribution, the situation will be more complex. The cosmic web is marked by a distinct luminosity, colour and morphology segregation. A strong and systematic trend of galaxy morphology with density has been established a few decades ago by Dressler (1980). Early type galaxies preferentially in rich groups and clusters and late types galaxies residing mainly in filaments and walls (e.g. Giovanelli et. al 1986). Other strong systematic clustering trends with luminosity and colour have also been established, such as in the complete SDSS DR7 galaxy sample (Zehavi et. al. 2010). A fair sample of galaxies should in principle take this intrinsic segregation into account. We are pursuing this in a forthcoming publication. In this study, mainly intent on establishing our reconstruction technology, we will consider the total galaxy population. The reconstructed density maps are therefore unlikely to represent a fair reflection of the underlying dark matter matter network, but will nonetheless convey the overall pattern of the large scale structure.
The reconstruction consists of two fundamental steps. The first is the estimation of the local galaxy density. The second is the interpolation of these density values to obtain a continuous spatial density field.
Following straightforward grid based interpolation methods and more sophisticated adaptive filter techniques, there has been a substantial investment into developing more advanced techniques. Examples of recently introduced techniques to follow the multiscale nature of the Cosmic Web is the Multiscale Morphology Filter route involves the wavelet reconstruction by Martínez et al. (2005).
Here we concentrate in particular on the density estimation and the subsequent field interpolation, in anticipation of the post-processing and feature detection studies of the SDSS DR6/DR7 survey, the subject of the additional papers in this series. To deal with the shot-noise and necessary selection effects in the resulting density fields, the analysis will involve filtering, a necessary post-processing step.

Density Estimation
The first step in the reconstruction procedure is density estimation. We show how this fits into our general scheme in Fig. 2. A similar, but slightly different outline can be found in Kitaura & Enßlin (2008).
In the literature we may find a large variety of density estimators (e.g. Silverman 1986;Sain 2002;Katkovnik 2002;Miller 2003). Usually they involve filters of some kind. Dependent on the filter kernel, density estimates may have a local character or include the information from a wider range of points in a point distribution. Filter kernels may have a rigid scale and shape, or they may be adapting themselves to the local point density. An example of a global estimator is a Gaussian kernel, which takes along the information -be it weighted -from distant points. A well-known example of rigid local estimators are the CIC and TSC grid interpolation formalisms (e.g. Hockney & Eastwood 1981). A considerably more flexible and adaptive filter kernel, frequently used in current N-body studies, is that of the spline-based interpolation recipes used in Smooth Particle Hydrodynamics codes (e.g. Monaghan 1992) whose scale is determined by the distance to the K th nearest neighbour.
Recently, various local and adaptive density estimators have been shown to provide highly favourable results for complex cosmological matter distributions, marked by a large dynamics range of scales and density values and intricate geometric patterns. One of these exploits the adaptivity of the kD tree (Ascasibar & Binney 2005;Sharma & Steinmetz 2006). We also note the use of the Epanechnikov kernel estimator in the identification of superclusters Einasto et al. (2007).
For the three reconstruction techniques addressed in this study, we use the local DTFE density estimate introduced by Schaap & van de Weygaert (2000) (see appendix B). It sets the density value at a sample point proportional to the inverse volume of the contiguous Voronoi cell around a sample point, i.e.. the sum of the Delaunay tetrahedra to which the sample point belongs (Schaap & van de Weygaert 2000;Schaap 2007;van de Weygaert & Schaap 2008). Tessellation-based methods are based on the realization that the optimal estimate for the spatial density at the location of a point xi in a discrete point sample P is given by the inverse of the volume of the corresponding Voronoi cell Okabe et al. (2000). They have been introduced by Brown (1965) and Ord (1978) and were first used in astronomy, for the specific purpose of devising source detection algorithms, by Ebeling & Wiedenmann (1993).
In an extensive study, Schaap (2007) demonstrated that the DTFE estimates are substantially better than those of the rigid grid based CIC and TSC techniques. Also, it outperforms the adaptive SPH spline performance, in particular in areas with substantial density gradients.

Interpolation
Interpolation on randomly scattered points aims to approximate a continuous function constrained by the available data points. A wide variety of approaches have been put forward, each with their own advantages and disadvantages. These methods can be roughly divided into two categories, global and local methods. Local methods have the benefit that they are fast and able to deal with large data-sets. Global methods tend to produce smoother interpolated functions but are computationally more expensive.
An overview and discussion of spatial interpolation techniques can be found in Press (2007) (chapter 3.4) and in Watson (1992). We refer to Lombardi (2002) for a detailed review of the statistical properties of a large number of techniques, and to Franke (1982) and Amidror (2002) for a detailed comparison between various methods.
RBF and Kriging interpolation methods use predefined kernels to interpolate the field. The kernel in RBF methods is a basis function that spans the space of all the interpolating functions. Kriging interpolation uses spatial correlations between sample points, with its kernel being equal to the corresponding covariance function. The method, introduced by Matheron (1963), is the best linear unbiased estimator of a density value given a set of measured sample points at irregularly spaced points. Given the commonly accepted fact that the primordial density perturbation field is Gaussian, we may therefore see the Kriging interpolator as the natural choice for reconstructing the cosmological density fields which emerged out of these primordial Gaussian circumstances.
There is a distant relationship of Kriging interpolation to Wiener filtering techniques (Wiener 1949;Rybicki & Press 1992). However, Wiener filtering is based on a different philosophy than Kriging, in that it includes a model for the noise and is evaluated in Fourier space. Also, classical Wiener filtering is predicated on an underlying Gaussian distribution. As a result, it has the serious disadvantage of suppressing or substantially diluting nonlinear structures of interest. More advanced recent developments and applications of Wiener filters to the reconstruction of the density distribution have largely remedied its capacity for reconstructing the density distribution (Kitaura & Enßlin 2008;Kitaura et al. 2009Kitaura et al. , 2010. A rather different approach is advocated in the triangulation-based interpolation techniques. Both the linear Delaunay Tessellation Field Estimator (DTFE Schaap & van de Weygaert 2000;Schaap 2007;van de Weygaert & Schaap 2008) and the higher-order Natural Neighbour Field Estimator (Sibson 1981;Watson 1992;Braun & Sambridge 1995) use the neighbourhood relationships defined by the Voronoi and Delaunay tessellation of the point sample to establish a fully adaptive, irregular and local interpolation grid. DTFE uses the Delaunay triangulation to reconstruct in a self-adaptive, mass conservative and parameter free way the underlying spatial (density) distribution. In combination with the DTFE tessellation-based sample points density estimates (see previous sect. 1.2.1), the DTFE interpolation leads to a volume-covering density field which has been shown to recover the hierarchical as well as the anisotropic morphology of the Cosmic Web (Schaap 2007;van de Weygaert & Schaap 2008).

DTFE, NNFE and Kriging
The DTFE Delaunay Tessellation Field Estimator method will be compared to two other techniques having a similar potential for a proper reconstruction of the cosmic web. The NNFE Natural Neighbour Interpolation technique shares the local nature of DTFE, but involves higher order interpolations. The third formalism is Natural Lognormal Kriging is a version of Kriging interpolation, and thus involves a nonlocal higher-order interpolation methodology.
The crucial step of importance in our investigation is the interpolation step (see diagram Fig. 2). Using the same initial sample point density estimates, the first step of the reconstruction procedure, we can assess the relative merits of DTFE, NNFE and Kriging interpolators.

Outline of this study
We start by describing the data samples in section 2, the SDSS survey sample and a mock galaxy sample which mimics the SDSS, obtained from the Millennium simulation. The local DTFE density estimate at the location of each of the sample galaxies is the subject of section 1.2.1. In the subsequent section we present and describe each of the three interpolation techniques investigated in this study, DTFE in subsection 4.1, NNFE in subsection 4.2 and Lognormal Kriging in subsection 4.3. The comparison throughout the paper is based on one specific sample, the density field reconstruction of the SDSS mock galaxy sample. Following a discussion in sect. 5 of the qualitative appearance of the density maps by each of the three methods, we turn to an intensive quantitative error and quality analysis of the density field in sect. 6. In sect. 7 this is followed by an investigation of the topological structure of the weblike galaxy distribution in the survey, mainly based on the void population as traced by the Watershed Void Finder. Section 8 presents the density field reconstruction of the SDSS data sample and in sect. 9.we summarize this study.
Following this first paper in a series will be a statistical study of the density field, focusing in particular on the one-point probability density function. Subsequently, we will present and discuss the cosmography of the reconstructed local Universe. Later studies will analyze the void population in the SDSS density field, and concentrate on the properties of galaxies as function of the large scale environment as characterized by our technique.

THE DATA
Our study is based on two major datasets. The principal one is the genuine SDSS DR6 galaxy redshift survey. For the purpose of understanding the errors and artefacts in the density field reconstruction we use a set of galaxy mock catalogues that model this SDSS dataset. The mock samples are obtained from the Millennium simulation (Springel et al. 2005).

The SDSS galaxy sample
For our analysis we use the main galaxy sample in the North Galactic Cap from the 6th data release of the Sloan Digital Sky Survey (SDSS) (Strauss et al. 2002;York et al. 2000;Stoughton et al. 2002;Adelman-McCarthy et al. 2008).
The SDSS DR6 data release consists of various contiguous regions. We restrict ourselves to the largest contiguous region, the northern strip of the North Galactic Cap (NGC) region. The sample was retrieved from the SDSS "casjobs server" (www.sdss.org) using the SQL query interface. Relevant properties were downloaded and most of the post-processing was done on a workstation. We did not attempt to assign galaxies with missing redshift due to fiber collisions. This gives a lower density estimated at the position of the missing redshift. However, given the sampling scheme used to resample the density field to a regular grid, the probability of sampling one of the affected areas is very low. The problem occurs only in the high density regions of the catalogue.
The spectroscopic SDSS galaxy sample is almost complete between a Petrosian magnitude limit of mr = 14.5 and mr = 17.77. We assume that the completeness of the sample does not vary significantly, even though there are in fact some angular variations (Blanton et al. 2003).
For our study, we select the SDSS galaxies that are located within a comoving box of 600 h −1 Mpc. In terms of the survey coordinates (X,Y,Z) (see app. A for definition), the observer is located at (X,Y,Z)=(300.,0.300.) h −1 Mpc and the centre of the northern strip is rotated to lie parallel to the Y-axis starting at (X, Z) = (300, 300) h −1 Mpc. In Fig. 3 the galaxies are plotted in the XY-projection (top) and the XZ-projection (bottom).
Within the (600 h −1 Mpc) 3 volume there are a total of 311474 galaxies with a magnitude less than the magnitude limit mr = 17.77.

The SDSS mock samples
The Millennium simulation (Springel et al. 2005) was used to construct galaxy mock catalogues which emulate the SDSS galaxy redshift sample. They are needed to get estimates of the errors induced by the magnitude selection, redshift distortions as a result of peculiar velocities and errors resulting from the survey mask.
We  The same (X,Y,Z) coordinate system and box size of 600 h −1 Mpc have been used for the mock galaxy samples. The mock galaxy catalogues are constructed as follows: • Periodically tile the Millennium cube to obtain enough cosmic volume.
• Calculate the redshift of the model galaxies wrt. the observer • Compute the apparent magnitude of each model galaxy from its absolute magnitude and redshift.
• Select the model galaxies brighter than mr = 17.77.
• Add the peculiar velocity to the Hubble redshift to obtain the total redshift • Apply the observational mask of the DR6-NGC sample to decide whether the model galaxy is included in the mock catalogue.

Redshift Space Distortions
In this study we assess redshift space density maps as well as partially corrected "real" space density maps.
Redshift space surveys like the SDSS are beset by distortions in the estimated distance. The result of large coherent cosmic flows, infall velocities onto clusters and highly nonlinear "thermal" velocities within clusters, these redshift distortions can have a dramatic effect on the estimated distances and reshape the large scale matter distribution. However for our purposes here, which is testing and comparing the methods, the presence or absence of fingers of God is immaterial.
Hence, in this paper, which deals with the Mock SDSS catalogues, we do not correct for the redshift distortions induced by large scale coherent cosmic flows. In particular the DTFE density reconstruction is locally adaptive and so the reconstruction is not corrupted by the presence of elongated radial features; this can be seen in figure 4.
The redshift distortions will, however, be addressed in a following paper in which we analyse the real SDSS catalogue.

Magnitude-vs. Volume-limited Samples
For the analysis of the SDSS sample, we extract two different samples from the full SDSS galaxy sample. These are a volume limited (i.e. absolute magnitude limited) sample and an apparent magnitude limited sample. Each sample is used for different aspects of our analysis.

Volume Limited sample
A volume limited galaxy sample is defined in order to assure a uniform galaxy coverage over the survey volume. A volume limited sample consists of a subset of galaxies which are homogeneously sampled throughout the sample volume. It has the advantage that each sample galaxy is an equal weight tracer of the underlying density field, and the resulting field will be statistically uniform. Our volume-limited sample has a distance limit of 300 h −1 Mpc and includes all galaxies brighter than Mr < −20.45, roughly representing the galaxies brighter than L * .
While the uniformity of the volume-limited samples assures a straightforward error assessment for any analysis, it has the disadvantage of losing the high spatial resolution represented by fainter galaxies nearby, therefore it does not necessarily have the smallest error.

Magnitude Limited sample
The magnitude limited sample contains all 311474 SDSS galaxies brighter than mr = 17.77. While a magnitude limited sample takes along all sampled information, one needs to correct for the inhomogeneous selection process.
A characteristic of magnitude-limited surveys is the change of intrinsic spatial resolution as we proceed out to larger distances. Potentially, this could be a serious issue when galaxies would be biased in a very complicated fashion (e.g. higher order biasing). This would render it very difficult to infer the density field at large distances, where only the most luminous objects would remain visible. By default, we therefore assume that all galaxies -independent of their luminosity-are a fair tracer of the density field.
Following this assumption, we correct for the dilution as a function of survey depth by weighing each sample galaxy by the reciprocal w(z) of the radial selection function ψ(z) at the distance of the galaxy. For the SDSS, the selection function ψ(z) as a function of redshift z is well fitted by the expression forwarded by Efstathiou & Moody (2001) where zr is the characteristic redshift of the distribution and β specifies the steepness of the curve. The corresponding number density N (z) of galaxies at redshift z is where A is a normalization constant and the z 2 term represents the increase of volume as function of z. The resulting galaxy weights w(z) are When incorporating w(z) the weights into the density field reconstruction, the normalization of the resulting density field may be achieved by modelling the details of the selection function and calculating the appropriate normalisation constant. We chose to follow the alternative of calculating the average of the reconstructed density field and subtracting it from the reconstruction. This is a simple and straightforward procedure, and perfectly valid as long as the volume of the sample is large enough to enable the estimate of the true average.

Mock Samples
From the mock catalogues we select three different samples.
The full mock catalogue comprises all Millennium (semianalytical) model galaxies within the DR6-NGC mask. The magnitude limited mock sample includes all model galaxies resulting from the procedure described above. The third sample is a volume limited set with an absolute magnitude Mr < −20.45.
To appreciate the differences between the magnitudelimited and the volume-limited galaxy sample, and also the full mock catalogue within the SDSS volume, the three mock catalogues are shown in Fig. 4.
The contours superimposed on the corresponding galaxy distributions are the resulting NNFE density field contours (see sect. 4.2).

LOCAL DTFE DENSITY ESTIMATE
Throughout this study we use the local DTFE density estimate, following the definition by Schaap & van de Weygaert (2000). In appendix B one may find more details of the DTFE procedure which we followed (for an extensive review see van de Weygaert & Schaap 2008). Implicitly, and for simplicity, we assign to each sample galaxy the same mass mi, i.e. the density value is predicated on the number density of galaxies.
The sample point DTFE density value is inversely proportional to the volume of of the local neighbourhood as defined by the Voronoi tessellation of the spatial galaxy distribution. (Schaap & van de Weygaert 2000) argued that the inverse of the volume of the contiguous Voronoi cell is the proper density estimate, assuring mass conservation for the subsequent linear interpolation step. The contiguous Voronoi cell, sometimes dubbed umbrella in the computational geometry literature, is the region defined by all Delaunay tetrahedra of which a given sample point is a vertex and which it shares with its natural neighbours. A twodimensional illustration of a contiguous Voronoi cell of a point is shown as the surrounding hatched region in the lefthand frame of Fig. 5.
The density value at each sample point is determined following the construction of the Delaunay triangulation (see e.g. Delaunay 1934; Aurenhammer 1991) 1 . Within the triangulation we identify for each sample point i all Ni neighbouring tetrahedra Tj (Fig. 5), which together constitute the contiguous Voronoi cell Wi ∪j Tj. Summation of the individual tetrahedral volumes V (Tj) yields the volume of the contiguous Voronoi cell, For the three-dimensional SDSS sample volume, the resulting DTFE estimate of the density fi at sample point i is (see equation B2), where the weight w(zi) is the sample selection weight at the galaxies' redshift zi (equation 3). Note that the factor four takes account of the fact that in three dimensions each sample point belongs to four tetrahedra. In practice, the density of all particles is calculated by looping in sequence over all Delaunay tetrahedra.

Shot noise errors
A local density estimator such as the contiguous Voronoi cell has the advantage of being very sensitive to the signal. However, this also implies them to have a high sensitivity to shot noise present in the data.
For appreciating the influence of shot noise in the DTFE density estimates, we turn to the probability distribution of the estimated intensity λ for a 3D Poisson point process with intensity λ = 1. Schaap (2007) found that it can be very well approximated by The first observation is that the DTFE density estimator is unbiased, as the mean of p( λ) is equal to one. An estimate of the error involved is that of the variance, σ 2 λ = 1/5, and is equal to ∼ 57%. Evidently, the distribution is non-Gaussian with a long tail extending towards high density values (cf. Fig. 7).

Centroidal Voronoi tessellations
The high density tail of the DTFE density estimate, and the implied shot noise level, can be suppressed or regularised by using the centroidal Voronoi tessellation (CVT) (see Lloyd 1982;Browne 2007). For a CVT the generating point distribution is such that the generating points are the mass centres of the resulting Voronoi cells.
The calculation of a CVT is usually done by means of an iterative procedure known as Lloyd iteration. Starting with an originally random point distribution, the centre of mass of the corresponding Voronoi cells is computed. Subsequently, the points are displaced to these centres. After a sequence of iteration steps, the resulting point distribution tends to converge to a proper CVT constellation. Effectively, the points have been repelling each other.
An impression of the CVT iteration procedure can be obtained from Fig. 6. Involving an initial point distribution, the resulting intensity distribution p( λ) for the tessellations obtained after zero, two, four, six, eight and 10 Lloyd iterations is shown in the righthand frame. Clearly, a CVT involves a much more regular distribution: after four iterations p( λ) has turned into a narrow and near symmetric distribution whose high-end tail is almost absent (Fig. 7). Potentially, a CVT might therefore help to suppress the density estimate error and its asymmetric distribution.

DTFE noise: a case study
A visual impression of the shot noise involved in the density estimate is provided by Fig. 8. It concerns a random sample of 4500 points distributed according to an anisotropic Gaussian distribution. This configuration is more representative for what may be expected in the real galaxy distribution.
Densities were estimated according to the pure DTFE procedure (see equation 6) and following a Lloyd CVT procedure of 5 iterations. The top row of Fig. 8 shows the density contours of the original peak and the raw DTFE and DTFE/CVT field reconstructions. The raw DTFE reconstruction (centre) is highly irregular compared to the original contours (left), while the CVT contours are much more regular (right). The linear profile along the y-axis (at x=0.5, bottom left panel) emphasizes the visual impression: the DTFE profile (orange) is marked by salient peaks which reflect the high density tail of p( λ), while the DTFE/CVT   . DTFE and DTFE/CVT reconstruction of an anisotropic Gaussian peak. Sampling the peak by 4500 random points, the top row shows the contour levels of the reconstructions: original (left), DTFE (centre) and DTFE/CVT with 5 Lloyd iterations. The resulting density maps after Gaussian smoothing, with R f = 0.05, is shown in the central row. The bottom row shows linear density profiles through the resulting mass distribution. Left: linear profiles along the y − axis at x = 0.5, for the original (black), DTFE (orange) and DTFE/CVT (blue). Central and Right: following the same colour schame, linear profiles through the filtered density field reconstruction, along the y-axis (at x = 0.5) and the x-axis (at y = 0.5).
profile (blue) appears to adhere considerably better to the original profile.
To appreciate the average trend in the density reconstruction, we filter the original, DTFE and DTFE/CVT fields with a Gaussian filter (R f = 0.05). In addition to the suppressed shotnoise, the resulting density level plots, shown in the central row of Fig. 8, reflect the way in which the reconstructions affect the shape of the Gaussian peak.
While the DTFE reconstruction retains the shape of the original, entirely in accordance with the findings by Schaap (2007), the DTFE/CVT appears not to do so. This impression is confirmed by the two linear profiles, along the y-axis, at x = 0.5 and along the x-axis, at y = 0.5, shown in the bottom central and righthand panels. While the DTFE reconstruction is able to follow accurately the shape of the original, the DTFE/CVT reconstruction displays considerable deviations.
We may therefore conclude that while CVT appears to suppress the shotnoise effects on small scales, it is not able to follow the morphology of the mass distribution on larger scales. Since we will be filtering our fields in a similar way, this result leads us to pursue our study on the basis of the pure DTFE density estimate procedure without the CVT regularization.

INTERPOLATION METHODS
The key aspect in our investigation is the interpolation step of each of the three reconstruction methods. Here we describe the interpolation steps of the investigated methods, the DTFE, NNFE and Kriging interpolators.

Interpolation Method I: Delaunay Tessellation Field Estimator (DTFE)
The most straightforward way to interpolate and/or reconstruct a density field is by linear interpolation between neighbouring data points.  & Schaap 2008) is the multidimensional equivalent of simple piecewise one-dimensional linear interpolation from an irregularly distributed set of points. DTFE generalizes the concept of natural interpolation interval to any dimension D by adopting the Delaunay tetrahedra of a multidimensional point set as such. It uses the adaptive and minimum triangulation properties of Delaunay tessellations to use them as adaptive spatial interpolation intervals for irregular point distributions (Bernardeau & van de Weygaert 1996).
Once the Delaunay tessellation has been constructed, and the densities at each sample point determined (see sect. 3), we determine the density gradient ∇f j within each Delaunay tetrahedron Tj from the density values (f (r0), f (r1), f (r2), f (r3)) at its four vertices at location r0, r1, r2, r3.
Using the density gradients in the Delaunay tetrahedra, the DTFE density value at any point r can be calculated by determining in which tetrahedron it is located and subsequently computing its density estimate f ( r) from the linear equation, To obtain an image of the density field, one calculates these density estimates at each of the voxel locations of the image grid. For a more detailed outline of the DTFE method we refer to section B. An impression of a DTFE interpolated field in a cosmological context is presented in Figure 9 (top righthand panel). It concerns a density field reconstruction from a dataset extracted from a Millennium mock sample (top lefthand panel). The galaxy selection follows the distant observer approximation, i.e. following parallel lines of sight, and assumes the magnitude limit mr = 17.77 of the SDSS redshift survey.
The resulting DTFE density field is the level map in the top righthand panel. DTFE recovers the fine small-scale structures and at the same time adapts itself to the larger scale structures at greater distances. It also reveals the linear interpolation artifacts, the triangular shaped low-intensity wings. These are especially noticeable when the data points are sparse. We must note that these wings are not significant in mass, but arise when one takes a lower dimensional (1 or 2) section through the data.

Interpolation Method II: Natural Neighbour Field Estimator (NNFE)
The DTFE is a piecewise linear interpolation (C 0 ) method. In a sense it is a linear version of a larger class of tessellation based interpolation methods. Of these, Natural Neighbour interpolation (Sibson 1981;Watson 1992;Braun & Sambridge 1995) is the most well known higher order tessellation based method (for more details see van de Weygaert & Schaap 2008).
The Natural Neighbour Interpolation formalism is a generic higher-order 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 (1981), while extensive treatments and elaborations of nn-interpolation may be found in Watson (1992); Sukumar (1998). As has been demonstrated by telling examples in geophysics (Braun & Sambridge 1995) and solid mechanics (Sukumar et al. 1998;Sukumar 1998) NN methods hold tremendous potential for gridindependent analysis and computations.
According to the Sibson natural neighbour interpolation, the interpolated value f ( r) at a position r is given by in which the summation is over the natural neighbours i of the point r amongst the data points (see Fig. 5, righthand frame). Note that a slight movement of the interpolation point will evoke a different set of natural neighbours. The Sibson natural neighbour interpolation uses areabased (or volume in 3D) interpolation weights λnn,i( r). These are determined from the volumes of the order-2 Voronoi cells V2( r, rj ). To understand the concept, imagine we virtually insert the location r in the spatial sample point distribution. Around its location a new Voronoi cell V( r) is delimited (see Fig. 5, where the cell is traced by the blue edges). The virtual cell V( r) overlaps with the original Voronoi cells Vj of its natural neighbours j. The order-2 Voronoi cells V2( r, rj ) are the regions of overlap, and define the region of space for whom r and rj are the closest "nuclei".
According to Sibson interpolation the interpolation kernel λnn,j ( r) is equal to the normalized area of the order-2 Voronoi cell, in which A( r) = j A2( r, rj ) is the area of the virtual Voronoi cell of point r and A2( r, ri) the area of the order-2 Voronoi cell V2( r, ri).
Evidently, the closer one moves a point r to a sample point rj , the more Voronoi cell V( r) will overlap with the original Voronoi cell Vj , and the larger the volume A2( r, ri) of the order-2 Voronoi cell V2( r, ri) becomes, thus increasing the weight λnn,i( r). When r finally coincides with one of the natural neighbours, the order-2 Voronoi cell will be identical to the old Voronoi cell at this point. The interpolated field value f ( r) will then be equal to the field value at that point.
Notice that the interpolation weights λnn,i are always positive and sum to one, This property is called partition of unity. The resulting function f (r) is continuous everywhere within the convex hull of the data, and has a continuous slope everywhere except at the data themselves. At the position of the vertices the derivative of the interpolant is discontinuous. In one dimension DTFE and NNFE are exactly the same. When the data-points are given on a regular grid, the NNFE reduces to the more familiar bi-linear (2d) or trilinear (3d) interpolation schemes. Our NNFE implementation is that of Eldering et al. (2006), a three-dimensional adaption of the two-dimensional version available in the CGAL library.
The NNFE density field reconstruction for the same Millennium mock sample as described in section 4.1 is shown in the bottom lefthand panel of Fig. 9. Some of the peaks in the regions with sparse sampling appear somewhat anisotropic. This is a consequence of the discontinuous derivative at the sample point. The overall resulting NNFE density field is well-behave and smooth, without the artifacts that beset the DTFE reconstruction.
A drawback of the DTFE and the NNFE methods is that neither take into account the existing spatial correlations that characterize the cosmological density field. These are explicitly taken into account by the Kriging interpolation technique.

Interpolation Method III: Kriging Interpolation
By basing itself on the covariance function of the density field, Kriging naturally includes the global spatial correlations of the field. The method was named by Matheron (1963) after D. G. Krige, who started the development of the method (see Cressie 1990, for a historical overview). The interpolator has the property that it is a best linear unbiased estimator (Cressie 1988(Cressie , 1993. Most applications of Kriging stem from the field of geostatistics, where Kriging found its origin. The applications concern measurements at irregularly scattered points which have to be translated into, for example, gold, ore or oil field reconstructions or into altitude maps. There is a distant relationship of Kriging interpolation to Wiener filtering techniques (Wiener 1949;Rybicki & Press 1992;Zaroubi et al. 1995;Zaroubi 2002;Erdogdu et al. 2004Erdogdu et al. , 2006Kitaura & Enßlin 2008). However, Wiener filtering is based on a different philosophy than Kriging, in that it includes a model for the noise and is evaluated in Fourier space. The retrieved field therefore corresponds to an optimally filtered field over a range of unknown scales. The filter scale is dictated by the locally estimated noise, with more noise corresponding to a larger amount of smoothing. An additional disadvantage for recovering nonlinear weblike features of classical Wiener filtering is that it is predicated on an underlying Gaussian distribution in the construction of the least squares estimator for filtering the data. While advantageous for the purpose of ascertaining the exclusive presence of significant features, it therefore has the serious disadvantage of suppressing or substantially diluting nonlinear structures of interest and may have difficulties in reconstructing the intricacies of the nonlinear structures.
More advanced recent developments and applications of Wiener filters to the reconstruction of the density distribution have revived its potential for reconstructing the density distribution. For an extensive and in-depth overview of these developments we refer to Kitaura & Enßlin (2008) (also see Kitaura et al. 2009Kitaura et al. , 2010.

The Kriging formalism
Interpolation can be viewed as estimating the field value f at location r by means of a weighted linear combination of nearby known data points f (ri); The main idea of Kriging, as originally formulated, is to calculate the values of weights λi that minimise the error with respect to the data according to the mean square variation, where E is the expectation over the specified quantity. We show in Jones et al. (2011) that this criterion can be replaced with the weaker requirement that the data and the errors be orthogonal in a statistical sense. It follows that the statistical distribution of the field f (r) need not be Gaussian distributed in order to achieve optimal reconstruction of the density field via Kriging. The Kriging equations for the weights λj are N j=1 C(ri, rj) λj = c(ri, r), where the matrix elements of C are given by and the vector elements c(ri, r) by From this linear system of N equations, it is straightforward to determine the N unknown weights λi.
While the matrix has to be inverted only once, the weights λi have to be specifically computed for each interpolation site r. After the weights have been determined, one can directly obtain the interpolated field values f from equation (11).

The Kriging Variogram
Usually, the covariance function E(f (ri)f ( rj)) depends only on distance, d = |ri − rj |. In geostatistics, this spatial dependence of random field is usually characterised by means of a variogram γ(r1, r2). The variogram is the mean square variation of the field values as function of distance, which for a stationary random field reduces to The variogram is related to the covariance function c(h), For practical purposes it is preferable to use a functional form for the variogram. There is variety of such variogram models (see Cressie 1993, for a detailed description). We use which represents a good fit for the variogram measured from a a Millennium SDSS mock galaxy survey (see sect. 2.2). To estimate the variogram for this galaxy distribution, we used the estimator (Cressie 1993), We base our estimate on a large number Np of randomly chosen locations within the sample volume.

Lognormal density fields
A lognormally distributed density field f has distribution where S is the S = log(1 + σ) and σ 2 = f 2 is the variance of the density field, (see Coles & Jones 1991).
What might simply be referred to as the "Lognormal Kriging" procedure uses the logarithm of the density field value to transform the density field data, Since the density ρ =ρ(1 + f ) is always positive, the application of the lognormal approach is valid everywhere and guarantees a positive definite reconstruction of the density. The final interpolation values are obtained by taking the inverse transformation, ie. the exponential of the interpolated data values, λi log (f (ri)) .
We will use the logarithmic value of the DTFE-interpolated field in what follows. For such a field, figure 10 shows the variogram for a DTFE interpolated field (red) and for fields Gaussian smoothed on scales R f = 1, 3, 6, 10 h −1 Mpc. The fitted variogram model parameters (equation 18) are listed in table 1.
In the second paper of this series on SDSS density field reconstructions, we will demonstrate that the galaxy distribution is indeed very well modelled by a lognormal distribution at scales in excess of 3 h −1 Mpc.

Localized Kriging
The value of N to be used in equation (11) has so far not been defined. We chose the local neighbourhood to be the tetrahedral natural neighbourhood. In 3 dimensions this is the union of all natural neighbours of the four vertices of the Delaunay tetrahedron in which a point r is located. This choice exploits the self-adaptivity to density and local shape of the Delaunay triangulation, and does not suffer the adverse effects mentioned above for the options of distance or number of neighbour selection (also cf. the discussions in Schaap 2007;van de Weygaert & Schaap 2008). Our experiments indeed confirm that the tetrahedral natural neighbourhood choice is superior to that of the 2 options listed above.
We found that in 3 dimensions the tetrahedral natural neighbourhood on average contains approximately 57 particles. One may extend the neighbourhood by adding a third or even more layers around it. An additional third layer would involve an average of 284 neighbours: however, the overall quality of the field reconstruction is not significantly better than with two layers despite a substantial increase in computational effort. See Jones et al. (2011) for more details on this.
Thus in our key equation (22) we use the value of N that is the number of first and second layer vertices surrounding each point.

QUALITATIVE DENSITY COMPARISON
For an assessment of the performance of the three reconstruction techniques, we turn to the mock catalogues modelling the SDSS DR6 galaxy survey sample. On the basis of the knowledge of the underlying density field of the mock samples, we will be able to infer absolute statements about the quality of the reconstructions.
In this section will first address the visual appearance of the density field reconstructions, which will allow a qualitative and global judgement on their ability to reproduce the true density field. A quantitative error and correlation analysis will follow in the subsequent section 6, while the topological properties of the reconstructions will be discussed in sect. 7.

Maps of the density field
Density contour maps for the three reconstructions are shown in Fig. 11 and Fig. 12. The maps concern thin slices through the density field reconstructions (and thus have zero intrinsic width).
The first set of three maps concern the reconstruction on the basis of the full Millennium mock galaxy sample within the boundaries of the SDSS DR6 North Galactic Cap region. The second set concerns the same region, but for the volume-limited mock sample.
To distinguish underdense and overdense regions, the underdense regions are indicated by means of a colour contour map while the overdense regions are represented by black contours. The contour levels in all maps are the same.
The DTFE, NNFE and Kriging maps all display the same global structure. The overall appearance of both the DTFE, NNFE and Kriging density maps is strikingly similar. In the higher density regions the differences are minor, and mainly concern the coherence and anisotropy of the reconstructed filamentary (and sheetlike) features. Qualitatively speaking, the lower density regions do reveal more differences between the three methods.
As expected, the volume-limited maps fail in reproducing the small-scale structure, while they trace the overall weblike outline seen in the full sample maps.

DTFE map
The DTFE map is remarkably accurate in outlining the tenuous weblike filamentary features, in particular in the case of the full sample reconstruction. Amongst the three reconstructions, the DTFE one looks more crispy than the NNFE and Kriging maps. It is slightly more capable in tracing the thin filamentary and sheetlike features, while one might have a slight worry with respect to the correct reproduction by NNFE and Kriging of the shape of filaments and walls. Also, we find that the DTFE maps are clearly marked by higher density contrasts, both within the overdense regions as well as with respect to the underdense regions.
The downside of the detailed structural reconstruction of DTFE is the more erratic nature of the DTFE contours, marked by sharp artifacts. These artifacts are particularly prominent in the field reconstruction of the volume-limited sample. They are manifestations of the linear interpolation method: when two neighbouring grid cells are located in different triangles, the field would appear to be discontinuous. The latter is mostly an impression, as sampling at finer scales would show it is just continuous. These artifacts occur in situations where the point sample density is considerably sparser than the size of the gridcells. This is also the reason why these artifacts are more prominent in the DTFE reconstruction of the volume-limited sample than in that of the full sample.

NNFE and Kriging
The reconstructed density maps of the two higher order schemes, NNFE and Kriging, have a considerably smoother appearance than the DTFE maps. The larger number of neighbours involved in the NNFE and Kriging reconstructions translate into the slightly more roundish contours of these structures. This also reveals itself in the absence of artifacts such as seen in the DTFE maps.
Part of the differences between the smoother higher order maps and the DTFE maps is an expression of the number of points, and field values, involved in the interpolation step (cf. eg. equation 22). DTFE uses 4 points, the vertices of a Delaunay tetrahedron. NNFE involves on average 17 natural neighbours, while the Natural Kriging scheme invokes 57 neighbours.
The smoother nature of the NNFE and Kriging maps is therefore an expression of the somewhat lower information content of these filtered maps. As a result, they also have less noisy low density regions than those seen in the DTFE maps.

Anisotropic Structure and Features
One of the crucial benefits of DTFE is that it is able to identify anisotropic features, like walls and filaments, and successfully reproduce their shape and morphology (see Schaap 2007;van de Weygaert & Schaap 2008).
From the density maps we see that NNFE and Kriging find the same filamentary structures. Overall, the impression is that DTFE and NNFE produce maps in which the cosmic web is more coherent than in the Kriging maps: the Kriging map mass concentrations have a slight tendency to break up more easily into clumps. This is true for both the full sample maps as well as the volume-limited sample maps. One exception, in the volume-limited map, seems to be the filamentary extension running from (X, Y ) = (170, 160) h −1 Mpc to (X, Y ) = (300, 225) h −1 Mpc. One reason for the somewhat more fragmentary nature of the Kriging maps is its use of a radially symmetric covariance function, while DTFE and NNFE are based on kernels that adapt to the local shape.

Underdensities & Voids
When turning to the underdense regions, we find that in the full density map both DTFE and Kriging delineate them at high contrast levels. The voids in the NNFE map the voids have a lower contrast. This is partially the result of the larger NNFE neighbourhood radius in low density areas. In this respect, the Kriging correlations assure a better performance. DTFE remains sensitive on behalf of its highly local character.
A comparison between the void population in the full sample reconstruction and that in the volume-limited sample reconstruction reveals a considerable contrast between the results for the different reconstruction methods. None of the volume-limited reconstructions contain the small voids   visible in the full sample maps. This is a reflection of the absence of such depressions in the diluted point sample. Of the three reconstructions, the contrast between the two DTFE maps is less distinct than that between the two NNFE and two Kriging maps. DTFE at least manages to trace the large voids at (X, Y ) = (180, 140) h −1 Mpc, at (X, Y ) = (350, 200) h −1 Mpc and at (X, Y ) = (460, 120) h −1 Mpc. Kriging and NNFE hardly manage to find the latter in the volume-limited maps, while the huge void complex near (X, Y ) = (350, 200) h −1 Mpc is a largely uniform moderate underdensity in the Kriging map. Only the DTFE map re-veals its true nature, a region marked by several deep voids embedded in a larger moderate undensity.

Density Profiles
To appreciate the small-scale details in the density field reconstructions, figure 13 displays a linear profile, along a radial distance of 500 h −1 Mpc, through the density field reconstructions. All profiles are taken along the Y-axis, at (X, Z) = (300, 300) h −1 Mpc. Also cf. Fig. 11 and Fig. 12.
The panels, from top to bottom, concern the DTFE, NNFE and Kriging reconstructions. The black lines are linear profiles through the full sample reconstructions, the superimposed red profiles concern the volume-limited mock sample reconstructions. We also show the linear profiles through these density field reconstructions, Gaussian filtered on a scale of R f = 10 h −1 Mpc. The gray lines are the linear profiles through the filtered full sample density field, the orange lines those through the filtered volume-limited sample density field.
The comparison between the linear profiles through the NNFE and Kriging reconstructions on the one hand, and the DTFE reconstruction on the other hand, leads to the following observation: • Underdense regions of NNFE and Kriging are less noisy. A nice example of this is the underdense region at X = 200 h −1 Mpc.
• Maxima tend to be wider in the higher order schemes.
• At larger radial distance, the topology of the unfiltered reconstructions is much smoother.
In the cosmological context, the density is approximately lognormally distributed (see Coles Kitaura et al. 2010), which impels us to work the log of the inferred density field. Our adaption of the Kriging method to tesselations of lognormally distributed discrete point data is presented in detail in a separate paper (Jones et al. 2011), where we also discuss the relationship of this procedure with the Constrained Random Field formalism developed by Bertschinger (1987) (also see Hoffman & Ribak 1991;Rybicki & Press 1992;Sheth 1995;van de Weygaert & Bertschinger 1996). We adopt the nomenclature used in the geostatistical literature.

QUANTITATIVE DENSITY FIELD ANALYSIS
For the quantitative comparison and error evaluation of the three density field reconstruction methods, we are particularly interested in the ability to recover the underlying density field from a magnitude-limited or volume-limited survey.
In this section we compare the DTFE, NNFE and Natural Lognormal Kriging reconstructions on the basis of magnitude-limited or volume-limited mock samples with that of the corresponding density fields determined from the full galaxy samples. The latter are galaxy redshift space samples that we would hypothetically obtain if we were able to observe all galaxies in the Millennium mock sample. The magnitude-and volume-limited samples are obtained from this Millennium mock sample by imposing the observational specifications of the SDSS survey (see sect. 2.1). An impression of the differences in structural resolution is provided by the visual comparison of the NNFE full, magnitude-limited and volume-limited sample reconstructions at the beginning of this study, in Fig. 4.
We will restrict the error analysis to the redshift space density maps, and not assess the errors introduced by peculiar velocity distortions. These are investigated in detail in the follow-up paper.
Most of our analysis focuses on the quality of the density field reconstructions obtained from magnitude-limited samples, unless specifically stated otherwise.

Magnitude-limited survey reconstructions: Correlation Diagrams
The first comparison between the magnitude-limited sample density field reconstruction and the full sample density field concerns a purely local point-to-point comparison. This test involves an inspection of the correlation diagrams of the density field value δ f ull (r) of the full sample density field at location versus the corresponding density value δ mock (r). Since this is a strictly local comparison, and does not distinguish between systematic nonlocal offsets or random field fluctuations, we try to get an impression of environmental effects by simply studying a sequence of filtered density fields. Of each density field we study four Gaussian filtered versions, at filter radii of R f = 0.0, 1.0, 6.0 and 10.0 h −1 Mpc. These scales represent the transition from the non-linearity (1 h −1 Mpc) to quasi-linear and linear scales at 10.0 h −1 Mpc.
If the survey-based reconstruction were perfect, the correlation diagram should be a perfect one-to-one line. The correlation diagrams show the level of scatter, and reveal whether the reconstruction errors are dependent on density and as well expose the presence of systematic offsets. Instead of a pure scatter diagram, we depict the density of the pairs (δ f ull , δ mock ) by means of contour levels. The contour levels in Fig. 14 depict the number density of pairs in logarithmic bins of size ∆ log(1+δ f ull )×∆ log(1+δ f ull ) = (0.02×0.02). The levels run from one pixel per bin (black) to the maximum number density (white), in logarithmic steps of ≈ 2.5. In each frame the black diagonal line indicates the exact one-to-one relation between full sample density field and the density field following from the magnitude limited sample.

DTFE correlation diagrams
For all filter radii, we find that over two orders of mag-nitude the correlation diagrams centre around a strict linear one-to-one relation. For the smaller filter radii (R f = 0.0 h −1 Mpc and R f = 1.0 h −1 Mpc) the relation appears to be slightly offset. This is a reflection of the necessarily small volumes probed by the galaxy survey at distances close to the observer, introducing cosmic variance effects through the locally estimated density normalization (see section 6.3). Evidently, the scatter around the linear 1-1 relation decreases as the filter radius increases. For both R f = 0.0 h −1 Mpc and R f = 1.0 h −1 Mpc, we find a more substantial scatter in the low density regions compared to the higher density areas. The scatter over the full density range turns into a more and more uniform level as we go to larger filter radii. Indeed, for R f > 3.0 h −1 Mpc and 0.1 < δ < 1 the agreement is very good.
A related additional trend is that from a slight upturn at low density values for R f = 1.0 h −1 Mpc towards a slight downturn at larger filter radii. In other words, at larger filter radii the reconstructions seem to have a systematic bias towards more underdense values.
However, we need to be careful in drawing general conclusions with respect to the low and high density extremes. In particular, for the large filter radii, these tend to be dominated by only a few rare objects. Because the offsets are relatively minor, we do not consider it a serious problem.

DTFE, NNFE and Kriging correlation diagrams
To compare the performance of the three methods, Fig. 15 shows the correlation diagrams for the density field reconstructions at a filter scale of R f = 3.0 h −1 Mpc. Each of the reconstructions shows an almost perfect one-to-one relation: all three methods yield unbiased reconstructions.
We cannot detect any large differences between the methods. This suggests that the remaining deviations of the reconstructed density fields are due to the initial DTFE density estimate. Again, for R f > 3.0 h −1 Mpc and 0.1 < δ < 1 the agreement is very good.
Alternative density estimators might lead to further improvements.

Magnitude-limited survey reconstructions: Intrinsic Smoothing Scale
A characteristic of the magnitude-limited survey is the change of intrinsic spatial resolution as we proceed out to larger distances and the galaxy sample includes only the most luminous objects. While one may correct the density estimates for the accompanying dilution (see sec. 3), the loss of small scale resolution and the accompanying geometric resolution is impossible to correct. This affects any study of filaments and walls on the basis of the reconstructed density field maps. The correction for this dilution effect is complicated by another effect, the intrinsic density-dependent resolution of the sampled galaxy surveys. Higher density regions are sampled by more galaxies, and are therefore more resolved than the low density void regions.
We may evaluate the intrinsic resolution of the density maps in terms of the intrinsic smoothing scale Rint. Effectively, it is related to the characteristic galaxy separation at Figure 16. Intrinsic Smoothing Scale and Density. Black contours: the correlation diagram between the R f = 3.0 h −1 Mpc filtered DTFE density field of the full galaxy sample (abscissa) and the unsmoothed DTFE density field reconstruction on the basis of the magnitude-limited survey (ordinate). Superimposed on the correlation diagram between the unfiltered full sample density field (abscissa) and the unsmoothed magnitude-limited mock sample density field (ordinate). The colours indicate the number of voxels that occupy the corresponding position in the correlation diagram, with white indicating the density pairs with highest occurrence and subsequent contour levels at logarithmic spacings of ∼ 2.5. a given redshift and position. It is rather straightforward to incorporate the effect of the increasing dilution as a function of redshift z, on the basis of the radial selection function ψ(z) (see equation 1). However, the considerable density contrasts in the inhomogeneous matter distribution render if far from trivial to find a correct expression for Rint(z). The mean galaxy separation is biased to high density regions and seriously underestimates the correct value Rint(z): the intrinsic smoothing scale is smaller than average in higher density areas, and larger in low density regions.

Intrinsic Smoothing Scale: Correlation Diagram
To appreciate the effect and the dependence of intrinsic smoothing length on density, we evaluate the correlation between the filtered (DTFE) density field of the full galaxy sample and the DTFE density field reconstruction on the basis of the magnitude-limited sample.
In particular, we are interested in the question in how far the filtered density field reconstructions are representative. We may assume they are as long as the corresponding filter scale R f > Rint(z). By evaluating the filtered density fields we may obtain an idea of the intrinsic smoothing scale on the basis of the density field error analysis: a sudden rise of error would indicate that Rint(z) > R f . By means of the black contours in Fig. 16 we show the correlation diagram of the full sample density field filtered at R f = 3.0 h −1 Mpc versus the unfiltered density field on the basis of the magnitude-limited galaxy sample. It is superimposed on the (colour) correlation diagram between the unfiltered full sample density field and the unfiltered density field reconstruction from the magnitude-limited mock sample. While the latter has substantial scatter over the full density range, the correlation between the filtered field and the magnitude-limited field appears to be much tighter. It confirms the impression that reconstructed fields resemble density fields that are filtered on a particular scale. The reconstructed density maps in Fig. 11 and Fig. 12 form a telling illustration of this observation.
The most outstanding feature of the filtered field correlation diagram is the curved nature of the contours, quite different from a regular linear relation. In the higher density regions we find that the densities in the magnitude-limited sample reconstruction are systematically biased to higher values. In the low to moderate density areas the trend is reverse: the densities of the raw survey density field tend to be somewhat lower than in the filtered field. This is a direct illustration of the density dependent nature of the intrinsic smoothing length Rint(z). Because the effective smoothing length is small in high density regions the low density regions get relatively over-smoothed and so the density is relatively over-estimated. More formally, it is an expression of the greater information content in the higher density regions of galaxy redshift surveys.

Filtered density field reconstructions
The filtering of a raw reconstructed density field will smooth out the high density values. The corresponding mass is redistributed to lower density regions. We argue that this can only yield correct density distributions if the nonlinear features in the mass distribution were reproduced at the correct position and with the correct amplitude in the raw sample density field reconstruction. Information from high density features and nonlinear objects is crucial for obtaining the correct large scale density field.

Radial Error Analysis
To quantify the reconstruction errors, we compare the local density values f (r) mock of the magnitude-limited survey reconstruction with that of the density f (r) of the full sample. We evaluate the absolute error ǫ1(r), as well as the relative error ǫ2(r),

Error Profiles
In figure 17 we plot the absolute errors ǫ(r) along the same radial direction as in Fig. 13. The corresponding DTFE density profile is plotted in the top panel, while the subsequent three panels depict the ǫ1 profiles for the DTFE, NNFE and Natural Lognormal Kriging density fields. The green lines are the error profiles of the unfiltered Figure 17. Linear Error profiles. The absolute error ǫ 1 (equation 23) as a function of distance, along the same axis as the density profiles in Fig. 13. Upper panel: the density profile through the DTFE density field (see Fig. 13 for description). The subsequent frames are the error profiles for the three reconstruction methods. The green line represents the error of the unsmoothed density field, the purple line that of the 10.0 h −1 Mpc filtered field. 2nd frame from top: DTFE density field; 3rd frame from top: NNFE density field; bottom: Natural Lognormal Kriging density field.
density field (black line in the upper panel), the purple lines are the error profiles of the 10 h −1 Mpc smoothed field (gray line in the upper panel). We find that the errors of all three reconstruction techniques are fairly similar, with a slightly increasing trend with distance for the unfiltered density field. At a distance of approximately 100 h −1 Mpc it is of the order of unity. Beyond this distance the errors are characterised by wide peaks that are mainly due to the magnitude-limited surveys undersampling of the density field at large distances. The errors in the 10 h −1 Mpc filtered field remain small over the whole reach of the linear profile, averaging an error level of around 10 percent. At intermediate distances, 200 h −1 Mpc < R < 400 h −1 Mpc, the DTFE errors are somewhat lower than those of the other methods. For all three methods, the largest errors are found at close distances. This is attributed to the survey mask, which is relatively thin in our immediate cosmic environment, ie. at distances R < 100 h −1 Mpc. Note that the filtered field errors at large distances appear to converge to the error profile of the Figure 18. Radially averaged error profiles. Top: the absolute error profile ǫ 1 (r) (solid line) and relative error ǫ 2 (r) (dashed line) for the magnitude-limited survey DTFE density field reconstruction. Going from dark to light blue, the profiles represent the error profiles through a sequence of filtered density fields: R f = 1.0 h −1 Mpc (black), R f = 3.0 h −1 Mpc, R f = 6.0 h −1 Mpc and R f = 10.0 h −1 Mpc (light blue). Bottom: The radially averaged relative error profile ǫ 2 (r), for the DTFE density field reconstruction (solid lines), the NNFE density field reconstruction (dashed lines) and the Natural Lognormal Kriging reconstruction (dot-dashed lines). The colours correspond to the same Gaussian filtered fields as specified for the top panel.
unfiltered field reconstruction. It is a reflection of the fact that the intrinsic smoothing length becomes comparable to the filter radius.

Radially averaged error profiles
By averaging the errors of the magnitude-limited survey density field reconstructions over radial shells we obtain an idea of error trends as function of distance. To minimize edge effects we exclude locations within 15 voxels from the edge of the survey volume.
The radially averaged error profiles for the DTFE density field reconstruction are shown in the top panel in Fig. 18. It concerns the absolute error ǫ1 (solid line) and relative error ǫ2 (dashed line) for four different filtered DTFE density fields. These are Gaussian filtered fields at filter scales of R f = 1.0, 3.0, 6.0 and 10.0 h −1 Mpc.
At nearby distances R < 100 h −1 Mpc the average error trend is that of decreasing absolute and relative errors. This is a reflection of the small and unrepresentative survey volumes involved. The fact that this effect appears to be most prominent for the largest filter scale, R f = 10.0 h −1 Mpc, is indicative. At this scale the reconstruction involves only a few independent wavevectors that constitute the resulting density field. Errors are enhanced by the way in which the weight factors w(z) (equation 3) are determined. We do this by fitting the selection function on the basis of the data (see sect. 2.4). At close distances this fit is heavily influenced by the presence or absence of large superstructures, which easily evokes systematic offsets in the local density estimate. It would probably be better to deal with such systematics on the basis of volume limited samples at close distances R < 150 h −1 Mpc.

DTFE, NNFE and Kriging error profiles
The bottom panel of Fig. 18  For all three methods the errors at the smallest filter scale, R f = 1.0 h −1 Mpc, are substantial, hardly ever better than 50%. The errors decrease rapidly with filter scale, such that at scales R f > 6 h −1 Mpc the density field can be reconstructed with reasonable accuracy. Beyond a distance of 100 h −1 Mpc the relative errors throughout the whole survey volume do not exceed the 10% level.
Overall, we find strikingly small differences between the three methods. None performs distinctly better than any of the others. On a more detailed level, we may observe two differences. Firstly, we see that the Kriging errors at small scales are relatively high. Secondly, DTFE appears to perform somewhat better at nearly all scales, except the smallest scale R f = 1.0 h −1 Mpc. In terms of density errors, the relatively simple and direct DTFE method seems to work best.

Volume limited survey density field reconstructions
A volume limited survey circumvents the resolution complications encountered in magnitude limited surveys. Instead of having to correct for the steadily decreasing resolution at larger distances, volume limited surveys have the advantage of a uniform sample resolution. This renders the geometric and topological analysis of structures considerably more straightforward. Moreover, in addition to removing the problem of varying spatial resolution, the major advantae of using volume-limited samples is that one is using the same type and luminosity of galaxies at different distance. The major disadvantage of volume-limited galaxy samples is their relatively low resolution, resulting from the necessity to uniformly sample galaxies throughout the sample volume. In essence, it involves a trade-off between spatial resolution and a sufficiently large and representative sample volume.
To test the performance of density field reconstructions on the basis of a volume-limited galaxy redshift survey, we take a volume limited sample from the DeLucia mock sample based on the Millennium simulation De Lucia & Blaizot (2007), restricted to a volume in between redshifts z = 0.02− 0.1 and containing galaxies brighter than Mr = −20.45.
The resulting galaxy sample can be seen in the bottom panel of Fig. 4 (along with the corresponding NNFE density field contours). With respect to the magnitude-limited sample (central panel), let alone the full sample, it clearly lacks Figure 19. Density field Correlation Diagrams.Left: Correlation diagram between full sample DTFE density field reconstruction (abscissa) and the DTFE density field reconstruction on the basis of the magnitude-limited sample (ordinate). Right: Correlation diagram between full sample DTFE density field reconstruction (abscissa) and the DTFE density field reconstruction on the basis of the volumelimited sample (ordinate). All fields are Gaussian filtered on a scale R f = 3.0 h −1 Mpc. Figure 20. Radially averaged error profiles. Top: radially averaged profile of absolute error ǫ 1 (r). Bottom: radially average profile of relative error ǫ 2 (r). Solid lines, both panels: the error of the DTFE density field reconstruction on the basis of the magnitudelimited sample. Dashed lines, both panels: the error of the DTFE density field reconstruction on the basis of the volume-limited sample (within the limits of the volume-limited sample volume). Each panel shows the error profiles for three different filter scales: R f = 1.0 (orange), 3.0 (red) and 10.0 h −1 Mpc (black). spatial resolution. This may also be appreciated from the three corresponding density field reconstructions in Fig. 12. The loss of small scale details is obvious: close inspection and comparison with the full sample reveals the loss of fine filamentary features separating small voids. They either merge into larger surrounding overdensities or get lost in enlarged void regions. Figure 19 compares the correlation diagrams for the DTFE density field reconstructions on the basis of magnitude (left) and volume limited samples (right). Overall, they occur rather similar, although there are some important details in which they differ. One particular one concerns the more extended maximum of the correlation diagram in the case of the magnitude limited sample density field. In other words, the errors of the magnitude limited sample density field are usually smaller than those for the volume limited sample density field. Figure 20 compares the profiles of the radially averaged absolute errors ǫ1(r) (equation 23, and of the radially averaged relative errors ǫ2(r) (equation 24) for the DTFE density field reconstructions on the basis of the magnitude-limited (solid lines) and on the basis of the volume-limited surveys (dashed lines). The error profiles are assessed for three Gaussian filter radii, R f = 1.0, 3.0 and 10.0 h −1 Mpc.

Radially averaged error profiles
The first observation from Fig. 20 is that the errors of the volume-limited sample are uniformly distributed throughout the survey volume, as might be expected for a statistically uniform sample. We also find that the errors of the magnitude-limited sample density field reconstructions are systematically lower than those for the volume limited sample. This remains so up to the edge of the volume limited sample, at R ≈ 300 h −1 Mpc, where the sampling density of the volume limited and magnitude limited sample are equal.

DTFE, NNFE and Kriging error profiles
When comparing the error profiles for the three different density field reconstructions we hardly find any significant differences. This is confirmed by Fig. 21, which present the radially averaged relative error profiles ǫ2(r) for the DTFE, NNFE and Kriging reconstructions at three different filter radii.
The uniformity of the errors appears to suggest that a major source for the observed errors has to be found in the density estimate itself, rather than in the interpolation technique. If there are any differences, we may argue that these concern a slightly better performance of DTFE.

TOPOLOGICAL ANALYSIS
To probe the global pattern of the mass distribution we turn to the topological structure of the SDSS survey. This is largely dependent on the higher order correlations in the density field. The error analysis described in the previous subsections would not necessarily be able to detect key differences in the large scale topology. In this subsection we will seek to evaluate the quality of the topological renderings of the density field, where we focus on the structure defined by the voids in the cosmic desnsit field.
There are various approaches to studying the topological structure of the large scale mass distribution. One option for characterizing the global topology of the cosmic matter distribution is in terms of four Minkowski functionals (Mecke et al. 1994;Schmalzing et al. 1999). These are solidly based on the theory of spatial statistics and also have the great advantage of being known analytically in the case of Gaussian random fields. In particular, the genus of the density field has received substantial attention as a strongly discriminating factor between intrinsically different spatial patterns (Gott et al. 1986(Gott et al. , 1989Park et al. 1992;Hoyle et al. 2002;Gott et al. 2008;Zhang et al. 2009). An attempt to extend the scope of the Minkowski functionals towards locally defined topological measures of the density field has been developed in the SURFGEN project defined by Sahni and Shandarin and their coworkers (Sahni et al. 1998;Shandarin et al. 2004). The main problem for these formalisms remains the user-defined, and thus potentially biased, nature of the continuous density field inferred from the sample of discrete objects.
Here we specifically address the topology of the SDSS galaxy distribution on the basis of the void population. To this end, we segment the galaxy distribution into void patches by means of the Watershed Void Finder (Platen et al. 2007). We test the watershed segmentation of the DTFE density field obtained from the magnitude-limited mock sample by comparing it to the watershed segmentation of the full galaxy sample density field. The topological errors are quantified according to the watershed segmentation of both fields.

Watershed Void Finder (WVF)
The Watershed Void Finder (WVF) is an implementation of the Watershed Transform for segmentation of images of the galaxy and matter distribution into distinct regions and objects and the subsequent identification of voids (Platen et al. 2007).
The basic idea behind the watershed transform finds its origin in geophysics. It delineates the boundaries of the separate domains, the basins, into which yields of e.g. rainfall will collect. The analogy with the cosmological context is straightforward: voids are to be identified with the basins, while the filaments and walls of the cosmic web are the ridges separating the voids from each other.
With respect to the other void finders the watershed algorithm has several advantages. Because it is 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.
The WVF consists of eight steps, which are extensively outlined in Platen et al. (2007). For the success of the WVF it is of utmost importance that the density field retains its morphological character. To this end, the two essential first steps relate directly to DTFE, which guarantees the correct representation of the hierarchical nature, the weblike morphology dominated by filaments and walls, and the presence of voids (van de Weygaert & Schaap 2008). Because in and around low-density void regions the raw density field is characterized by a considerable level of noise, a second essential step suppresses the noise by an adaptive smoothing algorithm which in a consecutive sequence of repetitive steps determines the median of densities within the contiguous Voronoi cell surrounding a point. The determination of the median density of the natural neighbours turns out to define a stable and asymptotically converging smooth density field fit for a proper watershed segmentation. The sub-  sequent central step of the WVF formalism consists of the application of the discrete watershed transform on this adaptively filtered density field.
A related tessellation-based method for void identification, ZOBOV Neyrinck (2008), does yield similar results as WVF (see Colberg et al. 2008). It demonstrates the successful application of tessellation-based techniques to identify structures within the cosmic matter distribution. In addition to the WVF and ZOBOV there is an array of void identification procedures (see e.g. Kauffmann & Fairall 1991;El-Ad & Piran 1997;Hoyle et al. 2002;Arbabi-Bidgoli & Müller 2002;Plionis & Basilakos 2002;Patiri et al. 2006;Colberg et al. 2005;Shandarin et al. 2006;Colberg et al. 2008). The "voidfinder" algorithm of El-Ad & Piran (1997) has been at the basis of most voidfinding methods. However, this succesful approach is not able to analyze complex spatial configurations in which voids may have arbitrary shapes and contain a range and variety of substructures, which lies at the heart of our analysis. Figure 22 shows the watershed segmentation generated by the mock catalogue. It is marked by the watershed boundaries superimposed on the contour maps of the density field of the full galaxy sample. These boundaries are visible as the white or thick black cell edges (dependent on the width of the boundary). The softer contour levels indicate the overdense regions in the R f = 2.0 h −1 Mpc filtered full sample density field.

WVF void population maps
For an impression of the ability to infer the correct watershed segmentation -and void population -from the magnitude-limited survey, we compare it with the segmentation for the full galaxy sample. In Fig. 23 the watershed boundaries of the latter are marked by red edges, while the ones for the magnitude-limited sample are marked by the black lines.
At distances up to R ∼ 200 h −1 Mpc there is overall a reasonable agreement between the two void segmentations. Beyond that distance, we find that the larger segments of the magnitude-limited sample -a consequence of the decreasing structural resolution of the survey -encompass several smaller void segments from the full sample. Beyond that distance we also find some regions with strong differences between the two segmentations. The segmentation within the large voids in the two zoom-ins illustrate this. In the magnitude-limited survey segmentation, the weaker bridges between the smaller voids visible in the full sample have vanished.
It is a clear illustration of the fact that at large distances the only topological information retained in the magnitudelimited density field is the skeleton defined by the strongest and most overdense filaments and walls, locations traced by the brightest galaxies. The more tenuous filigree of smaller filaments within the low density regions is lost.
7.3 Topological Error Definition: Figure 24. Definition of the false split measure. The four gray circles represent the patches from the reconstructed segmentation, the black circles would refer to a patch in the original segmentation. The gray shaded areas A j belong to the areas of the three reconstructed patches that have an intersection with the black circle. The dashed areas represent intersections between the original and the reconstructed segmentation A j∩O .

False-splits and False-mergers.
We evaluate the performance of the magnitude-limited survey watershed segmentation by comparing its void patches with those of the full galaxy sample.
As described extensively in Platen et al. (2007), the errors can be classified, to first order, by false splits and false mergers. A false split is the situation in which a void segment from the reference field splits into two or more voidpatches. The reverse situation is that of a false merger, where two void segments in the reference field merge into one voidpatch in the segmentation of the magnitude-limited survey.

False Split measure
We define a measure for the false split errors. Starting from the void segmentation in the full sample reconstruction, we find for each void segment the overlapping void patches in the magnitude limited segmentation. The relative fraction of overlap of these survey void patches with the original void segment is used as a measure of significance of the overlap. If there are two or more significant overlaps, we classify the configuration as a false split. Figure 24 illustrates the above. The large black circle represents a voidpatch in the original segmentation, while the four gray circles are void segments in the survey segmentation. The gray shaded areas belong to the areas of the three reconstructed patches j that intersect the black circle, while the resulting dashed areas represents the intersection of the original with each of these void patches. If the surface area of circle j is equal to Aj , and that of the overlap with the black original segment Aj∩O, then the false split measure decides on whether this is a significant overlap or not. We deem this to be so if f F S j 0.6. In Fig. 24, void1 and void2 would correspond to significant overlaps. Void3 would be excluded as such and would only correspond to a slight shift of the boundary. Because of the two significantly overlapping voids, this configuration would be classified as a false split.

False Merger measure
An almost equivalent measure can be defined for a false merger. Since by symmetry, a false merger can be considered a false split in the full (original) sample segmentation, we may simply reverse the definition of the false split measure. If l is a voidpatch in the original field, with area/volume A l , and its relative overlap with a large voidpatch in the magnitude-limited survey is we deem it a significant overlap when f F M l > 0.6. When there are at least two original void segments that overlap significantly with one in the magnitude-limited survey, we may consider it a false merger.

Correctly identified voids
Having defined the false split and false merger errors, we may specify the meaning of a correctly identified patch. A correct void patch in the magnitude-limited survey is a void segment which overlaps for at least 60% with a void segment in the original field, as well as the other way around. These two conditions prevent a void from being either a false split or a false merger. Figure 23 shows the spatial distribution of the topological errors. The image contains the watershed segmentation for the full galaxy sample, indicated by the black solid lines. These are superimposed on the watershed segmentation for the magnitude-limited survey, indicated by the red solid lines.

Spatial Distribution Topological Errors
The false mergers are indicated by orange patches. The dark red patches represent the false splits, while the correctly reproduced patches are represented as white cells. It is directly apparent that false mergers are far more abundant than false splits.
A visual comparison between the two segmentations (also cf. Fig. 22) reveals the disappearance of void boundaries seen in the original full sample matter distribution. They disappear within the large void regions found in the magnitude limited survey. In other words, these minima are absorbed by one large encompassing void. It is a result of the tenuous and usually underdense nature of the walls in these regions, so that only a few galaxies have to fall out of the survey to evoke a merging of voids.
Nonetheless, the coherence of the large persistent watershed lines remains strong throughout the volume. The defining skeleton of the cosmic mass distribution remains largely intact in a magnitude-limited survey.

Topological Error Characteristics
The quality of the void representation of the survey can be inferred from the percentage of correctly identified voidpatches, as well as from the percentage of topological errors, ie. the total number of false splits and false mergers. Here we assess the correct identifications and topological errors as a function of distance R. We accomplish this by counting the number correctly identified void patches in radial shells, along with the number of error patches in the same shells.
Note that the percentage of correct and of incorrect void identification do not always add up to 100 percent, since topological errors may be far more complex than just a false split and/or false merger. This happens with multiple additions or disappearances of void-walls.  In the case of the 3 h −1 Mpc field we see a gradual and continuous decrease of correct void identifications as we move outward, while there is a corresponding increase of erroneous identifications. The 10 h −1 Mpc filtered shows an increase of correct identifications at a very close range. This is a consequence of the rapidly rising ability to outline the corresponding larger underdensities as we expand towards a larger volume starting from the small nearby cosmic volume.

Filter Scale & Void finding
The DTFE reconstruction has been studied in somewhat more detail in Fig. 26. It plots the percentage of correctly and incorrectly identified void segments at a range of filter scales, running from For the larger filter scales, we find that the topology is reliably reconstructed within a range of R ≈ 100 − 300 h −1 Mpc, while for the smaller filter scales this seems to at a much closer range from R ≈ 0 − 200 h −1 Mpc. It is also interesting to see that on small scales, over the whole radial range, the number of correct identifications remains rather low and hardly exceeds 60%. On filter scales of 6 h −1 Mpc and 10 h −1 Mpc the performance is certainly better and consistently acceptable over a 200 h −1 Mpc range. Only beyond R ∼ 300 − 400 h −1 Mpc the number of correct identifications drops below the 50% mark.

Topology Range
In general, we find that the number of correctly identified voids decreases at larger distances, corresponding with a related increase of topological errors with distance. Dependent on the filter scale, we may define a distance Rmax out to which DTFE/WVF manages to identify a reasonable amount of voids from a magnitude-limited survey. From Fig. 26 we find: Within this range, the fraction of erroneous void identifications remains well below unity: Rmax is roughly estimated from Fig. 26 from the scale R at which the solid curves (correct identification) and dashed curves (erroneous identification) cross. Also, overall we find a slightly better performance by the DTFE reconstruction. Even though the higher order NNFE and Kriging techniques produce smoother void regions, this does not result in a higher number of correctly identified voidpatches.

SDSS-DR6 DENSITY RECONSTRUCTION
Following the extensive analysis discussed in previous sections, we arrive at the application of the assessed technologies to the real world galaxy distribution in the 6th data release of SDSS. These resulting density maps form the starting point of the extensive statistical study -starting with a study of the one-point density distribution function and galaxy bias -presented in the subsequent papers discussing the density field and its cosmography.
The DTFE, NNFE and Kriging density field reconstructions based on the (magnitude-limited) galaxy distribution in the 6th data release of SDSS are shown in Fig. 28. The density field has been reconstructed on a 512 3 grid, representing a spatial resolution of 1.17 h −1 Mpc.
For reference, in Fig. 27 we show the spatial distribution of the SDSS DR7 galaxies on which the reconstruction is based. They are superimposed as black dots on the DTFE density field.
Besides some minor differences, all three density maps show the same prominent features. The higher order NNFE and Kriging maps look smoother than the DTFE map. It is easy to recognize the almost one-to-one correspondence between the DTFE and NNFE map, with the DTFE map having the appearance of more noisy version. The Kriging map not only contains the same structures and features, it also has a more coherent appearance in which we can more easily recognize the global weblike morphology of the density field. Both the DTFE and NNFE maps have a more fragmented appearance.

SUMMARY AND DISCUSSION
This study is the first in a series in which we analyze the structure and topology of the Cosmic Web as traced by the Sloan Digital Sky Survey. In this study we investigate the ability of three reconstruction techniques to analyze and investigate weblike features and geometries in a discrete distribution of objects. The three methods are the Delaunay Tessellation Field Estimator (DTFE), Natural Neighbour Field Estimator (NNFE) and a local "natural" implementation of Kriging, the Natural Lognormal Kriging technique. DTFE and NNFE are based on the local geometry defined by the Voronoi and Delaunay tessellations of the galaxy distribution. The Kriging formalism is adapted and optimized for the the approximate lognormal density distribution encountered in the mildly nonlinear cosmic web and is based on the logarithm of the measured density values. Also, we have chosen to restrict the evaluations to a localized neighbourhood, the tetrahedral natural neighbourhood based on the Delaunay triangulation of the point sample.
The three reconstruction methods are analysed and compared using mock magnitude-limited and volumelimited SDSS redshift surveys, obtained on the basis of the Millennium simulation. We investigate error trends, biases and the topological structure of the resulting fields. The reconstructed density fields, mainly from the magnitudelimited mock survey samples but also from volume-limited ones, are compared with the density field of the total simulation galaxy sample. The differences between the various field reconstructions are investigated on the basis of an error analysis, mostly involving point-to-point comparisons. Environmental effects are addressed by evaluating the density fields on a range of Gaussian filter scales. With respect to the topology of the survey fields, we concentrate on the void population identified by the Watershed Void Finder (Platen et al. 2007). The void population in the full density field and in the magnitude-limited survey density fields are compared. The number of false mergers -in which original voids emerge as a part of a larger void in the survey field -and of false splits -in which an original void splits up in one or more voids in the survey field -forms the basis of the topological quality evaluation.
By investigating the quality of the resulting density field estimates, over a range of scales and in different environments, as well as the more global topological structure of the weblike network, we wish to identify and understand the qualities of these techniques for the different purposes and post-processing steps which are the subject of the following papers in this series. The following observations were made: • In most tests, DTFE, NNFE and Kriging have largely similar density and topology error behaviour.
• Cosmetically, higher order NNFE and Kriging methods produce more visually appealing reconstructions.
• Quantitatively, DTFE performs (marginally) better. Part of this at first sight surprising finding is a consequence of the higher sensitivity of the higher order NNFE and Kriging interpolation to intrinsic errors in the galaxy sample. An additional factor is the smaller natural neighbourhood of DTFE and NNFE with respect to the 3 to 4 times larger neighbourhood of Natural Lognormal Kriging, which restricts density errors to smaller volumes.
• With respect to the topological properties of the reconstructed density fields, it has become clear that a successful recovery of the void population on small scales is rather difficult. On these scales, the removal of only a couple of void galaxies leads to the spurious merging of observed voids.
• The void recovery rate improves significantly at filter scales > 3 h −1 Mpc. The immediate repercussion is that a proper analysis of small scale voids, and void galaxies, has to be necessarily restricted to the local Universe out to at most 100 h −1 Mpc. As the environmental influences on the galaxy formation process seem to be mostly determined on these scales (Park et al. 2007), our study within this project, subject of a forthcoming paper in this series, will restrict itself mainly to this volume.
A variety of technical improvements of our DTFE, NNFE and Kriging implementations may lead to a better performance. One immediate option is to invoke an image grid which has a more natural character for the galaxy survey context of our study than the cubic grid we have used in the evaluations described in this paper. Because of the flexibility of their definition, the Kriging formalism will be form a promising context for further adaptations and improvements. Of immediate importance is the implementation of non-local techniques to optimize the matrix inversion calculations while retaining the influence of large scale correlations and predictor-corrector methods to deal with oscillatory instabilities. While we have not yet investigated its performance, results of studies in other fields emphasize the importance of evaluating the performance of Radial Basis Function techniques as a possible alternative.
The DTFE, NNFE and Kriging density field reconstructions form the basis of a series of studies in which we analyze the cosmography, void population and spinal structure of the local Cosmic Web in the Sloan Digital Sky Survey. Given the ease and efficiency of calculation, and its good quantitative behaviour, for most of these studies we use the DTFE formalism. However, the Natural Lognormal Kriging results looks very promising and appears to produces a well-behave coherent weblike density map of the SDSS survey. With the large advantage of controlling the error behaviour and properties of the reconstructed map, along with the large potential for extensions and optimizations of the method, the Natural Lognormal Kriging maps will play a dominant role in our study of the Cosmic Web.
The density field reconstructions within the SDSS DR6, and subsequently SDSS DR7, volumes will be subjected to a statistical study in the second paper of this series. We will focus in particular on the 1-point probability function of the SDSS density field. With the help of corresponding magnitude-limited mock catalogues we will infer in how far the lognormal probability function, as well as related higher order versions, form a proper description of its statistical character. Also, the mock catalogues will allow us to assess the bias of the galaxy population in high density areas and, in particular, the low density void region. A cosmographic description of the reconstructed Local Universe, along with the identification of filamentary supercluster complexes, voids and supervoid complexes is the subject of the third paper in this series. The study of the void population in the SDSS density field, concerning a detailed assessment of the void size and shape properties, forms an important rationale behind the development of the density field reconstructions described in this study. This has become an interesting area of research, in particular as recent studies have emphasized the potential of extracting cosmological information from cosmic voids, in particular that concerning the dark energy equation of state (Lee & Park 2009;Lavaux & Wandelt 2009;Biswas et al. 2010). In recent years, there has also been a strong interest in large scale environmental influences on galaxy properties and on the galaxy formation process. The tools described in the present study, will allow an assessment on the basis of a properly defined density field on quasi-linear scales.
• Field values point sample The density values at the sampled points are determined from the corresponding Voronoi tessellations. The estimate of the density at each sample point is the normalized inverse of the volume of its contiguous Voronoi cell Wi of each point i. The contiguous Voronoi cell of a point i is the union of all Delaunay tetrahedra of which point i forms one of the four vertices (see Fig. 5 for an illustration). We recognize two applicable situations: + Uniform sampling process: the point sample is an unbiased sample of the underlying density field. Typical example is that of N -body simulation particles. For D-dimensional space the density estimate is, .
with mi the "mass" of sample point i. This situation concerns the "full" mock galaxy samples and the volumelimited galaxy samples.
+ Systematic non-uniform sampling process: sampling density according to specified selection process quantified by an a priori known selection function ψ(x), varying as function of sky position (α, δ) as well as depth/redshift. For D-dimensional space the density estimate is, f (xi) = (1 + D) mi ψ(xi) V (Wi) .
This situation is relevant for the magnitude-or fluxlimited SDSS and mock galaxy samples.
• Field Gradient Calculation of the field gradient estimate ∇f |m in each D-dimensional Delaunay simplex m (D = 3: tetrahedron; D = 2: triangle) by solving the set of linear equations for the field values at the positions of the (D + 1) tetrahedron vertices, The final basic step of the DTFE procedure is the field interpolation. The processing and post-processing steps involve numerous interpolation calculations, for each of the involved locations r. Given a location r, the Delaunay tetrahedron m in which it is embedded is determined. On the basis of the field gradient ∇f |m the field value is computed by (linear) interpolation, f (r) = f (ri) + ∇f m · (r − ri) . (B4) • Processing. 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 post-processing on the other hand. 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. For that purpose in this study we use the Monte Carlo approach: approximate the integral by taking the average over a number of (interpolated) field values probed at randomly distributed locations within the gridcell around an image point. The final estimate is obtained by averaging over the interpolated field values within a gridcell.
For image reconstruction we need to assure ourselves that we obtain a sensible density estimate within a voxel element. In a spatially irregular sample, the Voronoi cell defines the natural Nyquist interval. To avoid aliasing, the number of interpolation points therefore needs to oversample the voxels of the image grid. One may also take the alternative and exact option of piecewise integrating the density of the Delaunay tetrahedra that are (partially) overlapping with the image voxel. For DTFE this can be accomplished exact and relatively fast, although the computational and geometric aspects are far from trivial.

+ Smoothing and Filtering:
Linear filtering of the field f : convolution of the field f with a filter function W f (r, y), usually user-specified, fs(r) = f (r ′ ) W f (r ′ , y) dr ′ (B5) • Post-processing. 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. This can be finding voids, identifying cosmic structures or advanced statistical analysis of the density field.