1.

Introduction

Two-photon laser scanning microscopy¹ has become a principal technique for high-resolution fluorescence imaging in various biological tissues because it provides intrinsic optical sectioning and exceptional depth penetration (for reviews see Refs. 2, 3, and 4). Imaging depths in the cortex of more than 500 μm are now routinely achieved and image acquisition from 1 mm inside a mouse neocortex has been demonstrated recently.⁵ Combined with techniques for labeling individual neurons or sparse populations of neurons, e.g., dye loading via intracellular pipettes^{6 7} or the expression of fluorescent proteins,^{8 9 10} two-photon microscopy can resolve neurons with high resolution in vivo, i.e., within the intact brain of living animals.^{6 11} Thus, 3-D fluorescence images of neurons can be obtained including their entire dendritic morphology within their native environment.

These advances in imaging technology are prerequisites for the automatic reconstruction of neuronal morphologies. An automatic reconstruction would allow fast, high-throughput determination of characteristic anatomical features, for instance the dendritic branching pattern of different neuronal cell types. This is in contrast to standard manual reconstruction techniques, which are time-consuming and highly dependent on the experience of the human anatomist.¹² They also suffer from scaling problems due to shrinkage in fixed tissue. Automatic reconstruction would furthermore help to establish large databases of neuronal morphologies for biophysical modeling of cellular and neural network signal processing.

A number of different approaches for automatically reconstructing neuronal morphologies have been reported previously.^{13 14 15} Recently, Maravall et al.;¹⁶ automatically reconstructed a large number of layer 2/3 pyramidal neurons for an analysis of experience-dependent plasticity of their dendritic branching pattern. However, results in general strongly depend on the quality of the image data and in most cases require preprocessing. One of the major obstacles for developing an automatic reconstruction algorithm is the noise inherent in low-level fluorescence images. For example, using two-photon microscopy for in vivo imaging, both excitation light and fluorescence light are increasingly scattered with imaging depth, causing a reduction in signal-to-noise ratio and making it difficult to fully resolve thin, weakly fluorescent neural processes [Fig. 1(a)]. As a result, simple thresholding procedures for image segmentation might erroneously insert gaps into dendritic branches, preventing the reconstruction of a fully connected dendritic tree. Therefore, preprocessing of the raw fluorescence data in order to increase the signal-to-noise ratio while preserving dendritic structure is an essential prerequisite for automatic segmentation and subsequent morphological reconstruction.

Figure 1

Comparison of the anisotropic diffusion filtering with Gaussian blur. (a) One slice of a data stack from a two-photon scan of a pyramidal cell in layer 2/3 of rat somatosensory cortex in vivo. This slice is at a depth of 280 μm inside the brain. The loading pipette is visible at the top of the image next to the cell body. (b) Filtering with anisotropic diffusion. The filter was used with two time steps, a time step size of τ=2.0 and a scanning range of δ=10 μm. (c) Data after Gaussian blur with $\sigma =2.8\approx 2\sqrt{(2)}$ (measured in voxel).

One way to preprocess the raw data is filtering. A wide range of filters exist in image processing. The most basic filters calculate an average brightness value in a region around a central voxel. Other, more sophisticated filters use spectral analysis to extract signals within a defined bandwidth, such as low- or high-pass filters. Both methods show a close connection to the theory of partial differential equations.¹⁷ Actually the well-known Gaussian blur is an excellent low-pass filter.¹⁸ But none of these methods are sensitive to the local structure of the processed data.

One of the first approaches to include information about the data into the filter was made by the direction-pyramidal decomposition method.¹⁹ Alternatively, wavelet shrinkage has been used for preprocessing neuronal image data.^{20 21} One of the advantages of the wavelet approach is the use of multiple scales. A shortcoming is that the filtering threshold in the wavelet space cannot be determined automatically. One of the reasons is that the proper threshold depends on the local noise statistic within the image. Because this threshold varies throughout the data set an optimal value is difficult to find and thus no optimal denoising strategy can be given in general. A more general way to take the data structure into account is the use of diffusion filters which have a long tradition in image processing.¹⁷ They have mostly been used as convolution filters like Gaussian blurring. Nonlinear diffusion filters were first used by Perona and Malik in 1987.²² Since then many specific filters have been used to address a wide range of problems. For instance, an anisotropic diffusion filter controlled by local properties of the data was used by Lenzen²³ to reconstruct DNA structures.

In our case the filter has to prepare the data for segmentation. The filtering process aims to (1) separate noise and signal, (2) close apparent gaps in the structure, and (3) preserve dendritic diameters. For this purpose we designed an anisotropic diffusion filter which is now implemented in our software toolbox, NEURA (NEUron Reconstruction Algorithm). Here, we present this diffusion filter designed specifically for three-dimensional data of nerve cells. The filter is tuned to delete disturbances (remove noise) and to bridge open structures while preserving dendritic diameters. In this respect, it performs much better than Gaussian smoothing (see Fig. 1 for a first impression). The primary goal of our algorithm is to facilitate automatic reconstruction of neuronal morphology, for example in order to import them into the NEURON simulation environment²⁴ for in silico experiments.

2.

Material and Methods

2.3.

Structure Detection

Fluorescence images often are too noisy to use a gradient criterion to control the diffusion direction. More information is needed to reliably detect the object structure. In the case of filtering cellular processes, e.g., neuronal dendrites, we would like to use strong diffusion parallel to the main axis of the dendrite but not perpendicular to it. Thus we have to find a way to detect the axis of the dendrite locally. Motivated by Lenzen and Rumpf²³ we decided to use the moments of inertia.

2.3.1.

Moment of inertia

To determine the local structure of the data in a three-dimensional image stack we consider the gray value function as a density function of a real body. Then we can calculate the (local) moments of the virtual body, by choosing an integration volume B_δ(x₀) around the voxel of interest. The parameter δ represents the size of the integration volume and is referred to as scanning range.²⁸ The local moments are defined by:

• • mass:

  Eq. (6)

  $$M^0\left({\mathit{x}}_0\right)={\displaystyle {\int }_{{\mathit{B}}_{\delta }\left({\mathit{x}}_0\right)}}u\left(x\right)dx,$$

• • center of mass:

  Eq. (7)

  $$M^1\left({\mathit{x}}_0\right)=\frac{1}{M^0\left({\mathit{x}}_0\right)}{\displaystyle {\int }_{{\mathit{B}}_{\delta }\left({\mathit{x}}_0\right)}}u\left(x\right)xdx,$$

• • moment of inertia:

  Eq. (8)

  $$M^2\left({\mathit{x}}_0\right)={\displaystyle {\int }_{{\mathit{B}}_{\text{δ}}\left({\mathit{x}}_0\right)}}u\left(x\right)[x-M^1({\mathit{x}}_0\left)\right]\otimes [x-M^1({\mathit{x}}_0\left)\right]dx,$$

with ⊗ being the outer product.

The eigenvectors of the moment of inertia are the main axes of inertia. The eigenvalues contain information about the spatial structure. The size of the integration volume, the scanning range, is a critical parameter (see below).

2.3.2.

Example

The eigenvectors and eigenvalues of the moment of inertia of a hexahedron (see Fig. 2) can easily be calculated. The tensor calculated at the origin is:

Eq. (9)

$$M^2\left(x\right)=\frac{1}{12}{\displaystyle \left(\begin{array}{ccc}{a}^{3}bc& 0& 0\\ 0& a{b}^{3}c& 0\\ 0& 0& ab{c}^3\end{array}\right)}\text{.}$$

The corresponding eigenvectors and the eigenvalues are given by:

Eq. (10)

$${\mathit{V}}_1\left(0\right)={\displaystyle \left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)},{\mathit{V}}_2\left(0\right)={\displaystyle \left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)},{\mathit{V}}_3\left(0\right)={\displaystyle \left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)},$$

Eq. (11)

$${\alpha }_1=\frac{1}{12}{a}^{3}bc,{\alpha }_2=\frac{1}{12}a{b}^{3}c,{\alpha }_3=\frac{1}{12}ab{c}^3$$

Figure 2

Direction of the eigenvectors of the moment of inertia of a hexahedron. The structure is assumed to be solid with a constant density function. a,b,c: the three sides, V₁,V₂,V₃: the main axes of inertia.

2.3.3.

Geometry classification

Following Lenzen²³ we define the following variables to quantify the size of the eigenvalues:

Eq. (12)

$${\mathit{c}}_1=\frac{{\alpha }_1-{\alpha }_2}{\Sigma {\mathit{\alpha }}_{\mathit{i}}},{\mathit{c}}_2=\frac{2({\alpha }_2-{\alpha }_3)}{\Sigma {\mathit{\alpha }}_{\mathit{i}}}\text{and}{\mathit{c}}_3=\frac{3{\alpha }_3}{\Sigma {\mathit{\alpha }}_{\mathit{i}}},$$

the α_i are sorted by α₁>α₂>α₃. 1. Remark. The c_i in Eq. (12) are normalized in the following way: ∑c_i=1 and 0⩽c₁,c₂,c₃⩽1. Consequently, the c_i can be visualized by a state triangle; see Fig. 3. A high value of c₁ means that the local structure resembles a cylinder, a high value of c₂ a plane, and a high value of c₃ an isotropic structure. The parameters c_i can be further used for geometry classification as shown in the following example.

Figure 3

Visualization of the normalized values c_i using a state triangle.

2.3.4.

Example

We calculate the eigenvectors, eigenvalues, and the c_i for the structure from Fig. 4 in the voxel next to the gap. The integration volume is plotted around the voxel.

Figure 4

A simple linear structure to analyze the behavior of the eigenvalues and c_i respectively to the scanning range δ. The structure is again assumed to be solid with a constant density value. The structure is three voxel high (y direction) and three voxel deep (z direction).

Table 1 shows how the values of c₁, c₂, and c₃ change according to the integration size δ. The dominating eigenvector (DEV) changes when the integration size grows and finally they converge as can be seen from Table 1. For small integration volumes the algorithm tends to detect an isotropic structure (δ={3,5,7}) while for integration volumes large enough to reach over the gap (δ⩾9) a cylindrical structure is identified. Therefore it is important to know the scale of the structure to be detected. Consequently the algorithm can be tuned to detect small structures or large structures. For example, for enhancing dendritic branches we usually choose an integration size of 10 μm. For enhancement of smaller structures such as dendritic spines a smaller scanning range of about 3 μm has to be used.

Table 1

Dependence of the ci on the integration size δ in pixels.
scanning range	δ=3	δ=5	δ=7	δ=9	δ=11
c1	0.00	0.00	0.23	0.75	0.84
c2	0.53	0.00	0.00	0.00	0.00
c3	0.47	1.00	0.77	0.25	0.16
DEV	y,z	x,y,z	x	x	x

2.3.5.

Construction of D(u)

We are now in a position to construct the permeability tensor D(u). We define D(u) such that diffusion occurs only parallel to the axis of the tube, but not perpendicular to it. The dominating eigenvector V₁ of the moment of inertia gives us the local main direction of the structure. So we separate the diffusion direction ∇u into two, (∇u)_p and (∇u)_t, with (∇u)_p⋅V₁=0.

Technically speaking this means to transform the vector ∇u into the eigenspace of the diffusion tensor and multiply the first component with unity, the others with a nearly zero constant ε (which is necessary for well-posedness). So we finally set D from Eq. (5) to:

Eq. (13)

$$D\left(u\right)=B{\displaystyle \left(\begin{array}{ccc}1& 0& 0\\ 0& ϵ& 0\\ 0& 0& ϵ\end{array}\right)}{B}^T$$

with B=(V₁,V₂,V₃) and ε>0 [with V_i the eigenvectors; compare with Eq. (10)]. This means that the filter always causes diffusion into the direction of the largest eigenvalue of the local tensor of inertia, even if the local mass is small or nearly isotropic. We found this to be the optimal strategy when preparing neuronal image data for subsequent segmentation (see discussion in Sec. 4). However, for other applications of the filter it may be advantageous to switch off local diffusion if the data within the scanning range is close to isotropic, or to apply planar diffusion if c₂ is large.

3.

Results

3.1.

Testing the Filter on Artificial Data

We designed a simple data stack containing a Y-like structure to test how the filter works under controlled conditions. In Fig. 5(a) one slice of this 65³ voxel data stack is illustrated. The Y has three gaps, each of them being three voxels wide.

Figure 5

Test data stack. A simple test data stack with a Y-like structure was created to analyze the behavior of the filter. The structure is assumed to be solid with a density value of 1. The diameter of the Y-like structure is three voxels. The background is defined with density value zero. The structure has three gaps, each of which is three voxels wide. (a) One slice of the initial data. The red lines mark the area where the data for the signal profiles shown in (d), (e) were taken from. (b) Nonlinear anisotropic diffusion. (c) Linear isotropic diffusion. The normalized signal profiles demonstrate the effect of the anisotropic filter and the Gaussian smoothing. (d) Signal perpendicular to the branch. Anisotropic filtering conserves the signal profile. However, Gaussian blurring distorts the signal structure. Sharp edges are lost. (e) Signal along the branch. Both filters fill the gap in nearly the same way.

We applied the anisotropic diffusion filter to this data set, with a scanning range of 10 voxels [Fig. 5(b)]. As expected the filter preserved the diameters of the Y branches but filled in the gaps. For comparison we also applied a simple isotropic diffusion to the same data set [Fig. 5(c)]. Although the linear filter closed the gaps similar to the anisotropic filter, the diameters of the branches were severely broadened. Because of the enlargement of the structure and the conservation of the mean gray value the peak signal level of the isotropically filtered image is reduced compared to the anisotropic case.

To further quantify the differences between anisotropic and isotropic filtering, we measured the spatial intensity profiles parallel and perpendicular to one branch [Figs. 5(d) and 5(e)]. Broadening of this structure was analyzed by calculating the full width at half maximum (FWHM) for the raw data set, and for the anisotropically and isotropically filtered data set, respectively (using linear interpolation). The anisotropic filter maintained the half width of the branch (FWHM _raw =6, FWHM _anisotropic =5.84, measured in pixels). In contrast the FWHM was increased with linear filtering (FWHM _gauss =6.91.) . Because of the step-like nature of the artificial data the anisotropic diffusion even caused a slight reduction of the FWHM. The spatial profiles parallel to the branch demonstrate how both filters fill the gap in nearly the same way [Fig. 5(e)]. This is expected because along detected linear structures anisotropic diffusion and Gaussian blurring convolve the data nearly equally.

3.2.

Testing the Filter on Two-photon Images

We next applied the anisotropic diffusion algorithm to 3-D fluorescence image stacks of pyramidal neurons in rat neocortex obtained using two-photon microscopy.⁶ Neurons were filled with a fluorescent dye via the whole-cell patch pipette. The data typically consisted of several hundreds of slices taken at 2 μm focal increments, each consisting of 256×256 pixels (8-bit depth). A maximum intensity side projection of a raw data set is shown in Fig. 6(a). In addition, two example slices from different focal planes are shown [Figs. 6(b) and 6(c)] to illustrate how noisy the raw images are. In particular, thin basal dendrites close to the soma are difficult to resolve. Following anisotropic diffusion the signal-to-noise ratio was improved and the thin dendrites can easily be separated from background [Figs. 6(d)–6(f)]. To gain an impression of how the filter works we again compared the filter result to the Gaussian blur [Figs. 6(g)–6(i); σ=2.8 voxel] . The difference between isotropic and anisotropic filtering is particularly evident at the basal dendrites. After Gaussian filtering it is nearly impossible to see thin dendrites in slices deep inside the tissue, whereas the anisotropic diffusion further enhances the signal-to-noise ratio of these fine structures.

Figure 6

A layer 2/3 pyramidal neuron was imaged using two-photon microscopy. The data stack consists of 216 slices each 2 μm thick. Each slices has 256² pixels. One pixel is 1 μm in x direction and 1 μm in y direction. (a) Sagittal view of a maximum intensity projection of the raw data. (b), (c) Two slices of the raw data. The loading pipette is visible on the left side of the images. (d), (e), (f) Same data after anisotropic diffusion filtering. The broken red lines in (d) indicate the depth where the two slices shown are taken from. In (g), (h), (i) a Gaussian blurring is applied to the raw data.

For automatic reconstruction of the dendritic tree it is very important how dendritic diameters are affected by the filter. To quantify the diameters we again measured the FWHM of seven dendrites of different sizes. Some of them were parallel to the coordinate axis, some were oblique (see Fig. 7). The effects of these filters are shown in Table 3. The anisotropic filter has little effect on dendritic diameter (on average they shrink by a factor of 0.92±0.45) whereas Gaussian blurring almost doubles the diameter on average (expansion by a factor of 1.78±0.78) . In Fig. 9(a) the (normalized) signal profile of dendrite number 5 is shown.

Figure 7

Close-up views of dendritic branches. The conservation of the diameter of the dendrites was one of the major goals for the design of the filter. In this figure the dendrites which were used to measure the change in the diameter are shown. See Table 3 for results of the measurement of the full width at half maximum.

Table 3

Changes in diameter (FWHM) after anisotropic diffusion filtering or Gaussian blurring compared to raw size (in μm): f/raw is the ratio of the FWHM after anisotropic diffusion to the raw FWHM, and G/raw is the ratio of the FWHM after Gaussian blur to the raw FWHM. Seven different dendrites have been evaluated.
dendrite	1	2	3	4	5	6	7	mean	s.d.
raw	1.45	0.65	2.15	1.25	1.3	2.05	2.3
filtered	1.35	1.25	1.55	1.0	0.95	1.35	1.5
Gauss	2.1	2.25	2.95	2.3	2.35	2.8	2.65
f/raw	0.93	1.92	0.72	0.80	0.73	0.66	0.65	0.92	0.45
G/raw	1.45	3.46	1.37	1.84	1.81	1.37	1.15	1.78	0.78

We next examined whether the filter bridges gaps in two-photon microscopic data. Therefore we first enlarged one slice of the data set and measured the signal profile along one dendrite. In Fig. 8(a) the same raw data as in Fig. 6 are shown on an expanded scale. Figure 8(b) shows the data after anisotropic filtering. The box indicates the areas which are presented in Figs. 8(c) and 8(d). In Fig. 8(c) there are discontinuities in the dendritic branch as bright pixels are separated by dark pixels. After anisotropic filtering these discontinuities are closed [see Figs. 8(d) and 9(b)]. The filter effect is clearer when looking at a projection. Figures 8(e) and 8(f) are maximum intensity projections of 10 slices. Dendrites, which were difficult to follow in the raw data set, are easy to find in the filtered data, for instance, the dendrite in the upper left corner.

Figure 8

Zooming into one slice illustrates further how the filter enhances the quality of the data. One slice of the raw data from Fig. 6 is presented in (a). The red arrows mark the line where the signal profiles along the dendrite [see Fig. 9(b)] were measured. The blue box indicates the area from which panel (c) and panel (d) respectively were taken from. (b) Data after nonlinear anisotropic diffusion filtering. (c) and (d) Initial data (filtered data) at higher magnification. (e) and (f) Maximum intensity z projections of ten slices of the raw and the filtered data respectively are illustrated.

Figure 9

Conservation of dendritic diameters and closure of gaps in the two-photon fluorescence data stack. The plot in (a) shows the normalized signal profile perpendicular to dendrite 5 in Fig. 7. Again the anisotropic diffusion filtering conserves the signal profile whereas the Gaussian blurring causes a widening of the signal. In (b) the closing of a gap is demonstrated by showing the (normalized) signal along the dendrite marked in Fig. 8(a) by red arrows. There is nearly no difference between anisotropic diffusion and blurring.

3.3.

Application to Binary Data

One of the purposes of the filter is to preprocess the data for segmentation of objects. Even with this preprocessing some gaps may remain after a segmentation process. In this case it is convenient to apply the anisotropic diffusion filter a second time on the segmented (i.e., binary) data. An example application of the filter on binary data is shown in Fig. 10. The filter bridges the remaining gaps while preserving the diameter of the dendrites.

Figure 10

Nonlinear anisotropic filter applied to binary data. During the reconstruction process filtering on binary data is needed. These figures illustrate how the filter can be used on binary data, too. The initial binary data was obtained from the same raw data as used in Fig. 6. (a) Initial binary data. The data was filtered by nonlinear anisotropic diffusion with τ=2.0, two time steps, and δ=10 μm. The result is shown in (b). The same region as in Fig. 6 was chosen for analyzing the data at higher magnification (c), (d). The gaps in (c) are closed in (d). The diameter is nearly unchanged.

To illustrate how the filter facilitates automatic reconstruction of the dendritic morphology, we compared reconstructions of the same layer 2/3 pyramidal neuron without and with prior filtering (Fig. 11). Filtering included a first pass of anisotropic diffusion filtering on the raw data, subsequent segmentation, and a second pass of the filter on the segmented (binary) data. The unfiltered and filtered data were then both segmented, reconstructed using the same algorithm and parameter settings, and visualized using the NEURON simulation environment.²⁴ In both reconstructions the main apical dendrites are well detected. However, the fine basal dendrites and the thin side branches of the neuron are missing in the case of no pre-processing [Fig. 11(a)] whereas these structures are present following application of the anisotropic diffusion filter [Fig. 11(b)]. This is particularly evident for the fine basal dendrites deep within the neocortex. A detailed discussion of results obtained with different reconstruction procedures will be presented in a forthcoming paper.

Figure 11

Automatic reconstruction of neuronal morphology from two-photon fluorescence data (same layer 2/3 pyramidal neuron as in Fig. 6). The morphology of soma and dendrites is presented using the simulation program NEURON.²⁴ (a) Result of segmentation and reconstruction without prior anisotropic diffusion filtering. (b) Result of filtering, segmentation, and a second pass of the filter over the segmented (binary) data, followed by the same segmentation and reconstruction algorithms with the same parameter settings as in (a).

4.

Discussion

We have presented a method for filtering two-photon microscopy data of neuronal morphologies. The method is based on anisotropic diffusion in three dimensions. The key idea of the algorithm is the use of the local moments of inertia to reliably detect the dimensionality (isomorph, planar, or linear) and the orientation of the morphological structures, particularly of dendrites. The filter is controlled by adapting the diffusion tensor in a way to produce smoothing along but not perpendicular to the structure. We have demonstrated that anisotropic filtering preserves dendritic diameters in a set of artificial test data and in two-photon microscopic images of neocortical neurons. In addition it bridges apparent gaps in images of dendrites, which result from poor signal-to-noise ratio, by smoothing along the dendritic axis similar to a Gaussian blurring filter. It should be noted that dendrite diameters can only be resolved accurately if they are larger than the optical resolution of the microscope system. This is a fundamental limitation, which is not overcome by anisotropic diffusion filtering. For thin dendrites the diameters obtained therefore represent an overestimate of the true diameters because they are convolved with the point spread function of the microscope. However, in any case anisotropic diffusion filtering does not introduce an additional systematic bias in dendritic diameters.

The anisotropic diffusion filter was implemented using a semi-implicit scheme for time discretization and a finite-volume method for spatial discretization. The semi-implicit time scheme is advantageous since each time step can be used with a larger time step size without losing stability.²⁹ In contrast, explicit schemes for time disretization are limited in the time step size by the Courant-Friedrich-Levy condition.²⁹ Smaller time steps have to be used so that the number of time steps required to achieve the same filtering increases. As a result the time to solution is typically longer for an explicit scheme compared to an implicit scheme. Due to the complexity of the filtering algorithm, anisotropic diffusion filtering of large data sets is time consuming. But in some parts of the algorithm parallelization is possible, for instance for the calculation of the tensor of inertia or the assembly of the stiffness matrix [see Eq. (17)]. In the future, such parallelization will substantially reduce the computation time.

Anisotropic diffusion filtering offers a convenient way to enhance the quality of fluorescence image stacks without using information about the microscope system used for imaging. The signal-to-noise ratio is enhanced such that the differentiation between foreground and background (segmentation) can easily be done by a simple local thresholding algorithm. The smooth signal profile along structures and the large gradient between structure and background allow even the detection of small dendrites several hundred micrometers deep inside the cortex. Anisotropic diffusion filtering therefore represents an excellent starting point for automatic reconstruction of neurons. In addition, filtering and segmentation can be iterated since the filter does not depend on the depth of the data. It works both on grayscale and on binary data. The advantage of our method compared, for instance, to the method of wavelet shrinkage is that no a-priori information about the given data set is needed. Thus, for a specific structure a general filter algorithm can be used, offering the possibility of automated filtering of large numbers of data sets. Nevertheless it will be interesting to compare the present approach to the wavelet shrinkage method used in Ref. 21. Mra´zek et al.;³⁴ showed for the one-dimensional case that the Perona-Malik diffusivity or the Weickert diffusivity have a corresponding wavelet shrinkage. This means that certain diffusion processes can be described by wavelet shrinkage. In the future it may be possible to extend this connection between anisotropic diffusion filtering and wavelet transformations to the three-dimensional case.

For those who are interested in trying the anisotropic diffusion filter by themselves we offer a sample program in Ref. 35.

We are currently developing a new software toolbox, NEURA, designed for automatic reconstruction of neuronal morphologies, which includes the anisotropic diffusion for data preprocessing. Preliminary results indicate that the combination of the anisotropic filter with a sophisticated segmentation algorithm can deliver a good restoration of the neuron shape. Even fine dendrites can be tracked using graph tracking algorithms, yielding a reconstruction of the branching pattern of the neuron. Here the next challenge will be to apply the reconstruction algorithm to the finer axonal arborisation of cortical neurons. The overall algorithm for neuron reconstruction implemented in NEURA will be discussed in a forthcoming paper.