2.1.

Sample Preparation

All cells were filled in Wistar rats (P28 to P31; Charles River Laboratory). Experiments were carried out in accordance with the animal welfare guidelines of the Max Planck Society.

Briefly, animals were anesthetized with urethane or isoflurane. Body temperature was maintained at $37°\mathrm{C}$ by a servo-controlled heating blanket. A metal post for positioning the head was attached to the skull overlying the cerebellum by dental acrylic. A craniotomy was made overlying the left barrel cortex (each $<0.5{\mathrm{mm}}^2$ ) or the ventral posteriomedial nucleus of the thalamus $(<2.0{\mathrm{mm}}^2)$ . Thalamic craniotomies were centered at $3.5\mathrm{mm}$ posterior to bregma and $3.0\mathrm{mm}$ lateral of the midline. The dura was removed for some of the experiments. The craniotomy was kept covered with artificial cerebral spinal fluid [aCSF; in millimoles: 135 NaCl, 5.4 KCl, 1.8 $\mathrm{Ca}{\mathrm{Cl}}_2$ , 1.0 $\mathrm{Mg}{\mathrm{Cl}}_2$ , and 50 $N$ -(2-hydroxyethyl)piperazine- $N^{\prime }$ -2-ethanesulfonic acid (HEPEs); pH 7.2].

Cells were filled with biocytin either extracellularly using juxtasomal recording and electroporation²² or via whole-cell recording.²³ In both cases, we used patch pipettes, pulled from unfilamented borosilicate glass on a Sutter P-97 puller (Sutter Instruments). The outside diameter of the shank entering the brain varied from $25\text{to}75\mu \mathrm{m}$ , and the tip opening had an inside diameter less than $1\mu \mathrm{m}$ .

Whole-cell recording pipettes were tip-filled with (in millimoles) 135 K-gluconate, 10 HEPEs, 10 phosphocreatin-Na, 4 KCl, 4 adenosine triphosphate–Mg, 0.3 guanosine triphosphate, and 0.2% biocytin (pH 7.2, osmolarity 291). Pipettes were inserted perpendicular to the pia under high pressure $\left(200\text{to}300\mathrm{mbar}\right)$ , and cells were searched for blindly under low pressure $\left(20\text{to}30\mathrm{mbar}\right)$ .²³ Recordings were made in bridge mode for $30\text{to}60\min$ during which biocytin passively diffused into the cells. Seal resistance was $>1\mathrm{G}\Omega$ , access resistance was $1\text{to}100\mathrm{M}\Omega$ , and spike height and overall $V_m$ were stable throughout the recording. No holding current was used.

Juxtasomal fillings were made with pipettes filled with aCSF containing 1% biocytin. No pressure was applied to the pipette at any time. Once a cell was isolated extracellularly, current pulses ( $1\text{to}5\mathrm{nA}$ , $200\mathrm{ms}$ on, $200\mathrm{ms}$ off) were applied continuously for several minutes.²²

At least $1\mathrm{h}$ was allowed to pass after cell filling and prior to tissue fixation to ensure sufficient diffusion of the biocytin throughout the axonal arbor. The animal was deeply anesthetized and perfused transcardially with phosphate buffer followed by 4% paraformaldehyde. Cortex and thalamus were cut tangentially and coronally, respectively, in $100\text{-}\mu \mathrm{m}$ vibratome sections. Biocytin in these sections was stained with the chromogen $3,3^{\prime }$ -diaminobenzidine tetrahydrochloride using the avidin-biotin-peroxidase method.¹⁶ Sections were sometimes counterstained for cytochrome oxidase.²⁴ Processed sections were then mounted on slides and coverslipped with Mowiol (Hoechst, Austria). For comparison with our reconstruction approach, four well-filled cells from four different rats were manually reconstructed using a Neurolucida system with a $100\times$ oil immersion objective [Olympus $100\times$ Plan; 1.25 numeric aperture (NA).]

2.2.

Image Acquisition

A standard TLB microscope (Olympus BX-51, Olympus, Japan) equipped with a motorized $x\text{-}y\text{-}z$ stage (Maerzhaeuser Wetzlar, Germany) was used for image acquisition [1 in Fig. 1a]. A $546\text{-}\mathrm{nm}\pm 5\text{-}\mathrm{nm}$ bandpass illumination filter (CHROMA AF-analysentechnik, Germany), which is attached to the diaphragm of the lighthouse, provides quasi-monochromatic illumination of the specimen [2 in Fig. 1a]. This bandpass filter minimizes the chromatic aberrations of the imaging system and simplifies the theoretical description of the optical pathway from a polychromatic to a monochromatic one.

The light is transmitted by a high NA (1.4 NA) oil immersion condenser (Olympus, Japan) [3 in Fig. 1a], ensuring parallel illumination of the specimen under “Koehler”-conditions.²⁵ The specimen is imaged by a $100\times$ high NA oil immersion objective (Olympus $100\times$ UPLSAPO; 1.4 NA) in combination with a $0.5\times$ television (TV)-mount (Olympus U-TV0.5XC-3) or a $40\times$ objective (Olympus $40\times$ UPLFLN; 1.3 NA) in combination with a nonmagnifying TV mount [4 in Fig. 1a]. The immersion oil has a refractive index of $n_{\mathit{\text{oil}}}=1.516$ similar to glass.

The stage is navigated in three spatial directions by OASIS-4i-controller hardware and software (Objective Imaging Ltd., Cambridge, United Kingdom). It allows the acquisition of large mosaic images¹⁹ at different focal planes. Mosaic in this context refers to a two-dimensional (2D) image of overlapping tiles (i.e., adjacent fields of view) that are aligned automatically and then stitched during the image acquisition, resulting in a large composite image [Fig. 2a]. The user-defined scan area is automatically divided into a series of overlapping fields of view that we call tiles. For each tile location, a stack of images is acquired using optical sectioning²⁰ with a typical separation of $0.5\mu \mathrm{m}$ between the focal planes. For each focal plane, the corresponding tiles are then stitched together, resulting in a 3D stack of 2D mosaic images for each focal position [Fig. 2b]. This process is executed by the Surveyor Software (Objective Imaging Ltd.).

The images are recorded by a Q ICAM Fast 1394 camera (QImaging) [5 in Fig. 1a] equipped with a CCD chip, which in combination with the $100\times$ objective and the 0.5 TV mount yields an $x∕y$ sampling of $92\mathrm{nm}\text{per}\text{pixel}$ ( $116\mathrm{nm}$ per pixel if the $40\times$ objective and a nonmagnifying TV mount are used). Because the illumination is limited to one wavelength, eight-bit gray value images, rather than red-green-blue color images, are acquired. Other cameras, such as the Retiga 2000R Fast 1394 Mono Cooled (QImaging), were tested as well and showed no significant improvement in resolution or signal-to-noise ratio.

To guarantee a similar dynamic range of the gray values for mosaic image stacks from different brain slices and animals, the exposure time of the CCD camera is set semiautomatically by the Surveyor Software. Therefore, the mean gray value in a typical field of view within the scan area is calculated and is set to be 190. Typical here refers to a field of view without any stained neuron somata or blood vessels. This highly important feature of the image acquisition is the basis for an optimal deconvolution and, therefore, a robust neuron tracing. However, significant exposure gradients within the scan area are a considerable constraint to our method. The consequences and limitations for neuron tracing will be discussed later.

In summary, image acquisition results in a high-resolution 3D mosaic image. The typical image stack for neuron tracing that we use is $1\mathrm{mm}\times 1\mathrm{mm}\times 100\mu \mathrm{m}$ . However, the mosaic pattern can be adjusted to smaller or larger areas. Using the $100\times$ objective in combination with the 0.5 TV mount or the $40\times$ objective with a nonmagnifying TV mount, the sampling is $92\times 92\times 500{\mathrm{nm}}^3$ or $116\times 116\times 500{\mathrm{nm}}^3$ per voxel. Hence, the data volume for such a stack is approximately $15\mathrm{GB}$ .

2.3.

Hardware and Software Requirements for Automated Image Processing

In general, the data size is between 10 and $30\mathrm{GB}$ per section. Therefore, a fast image processing pipeline based on multi processor computing that can handle such large data sets was developed. The algorithms are written in $C++$ .²⁶ The raster image file input/output, iteration through a raster image and the dilation as well as the closing filter, use the ITK Image Processing Library.²⁷ The algorithms were parallelized by applying the OPENMP standard.²⁸ They are executed on AMD dual-core $64\text{-}\text{bit}$ Opteron servers, equipped with either 4 central processing units (CPUs) and $32\mathrm{GB}$ memory (DELTA Computer Products GmbH, Reinbek, Germany) or 8 CPUs and $64\mathrm{GB}$ memory (fms-computer.com, Netphen, Germany).

Once the mosaic image stack has been saved to disk, a software daemon (a custom written PERL script) detects new data on the hard drive and starts the processing pipeline. Additional image stacks are added to a queue and processed sequentially. The script initiates the subsequent image processing steps: inversion of the gray values, subdivision into bricks, deconvolution, down sampling, and finally the 3D neuron tracing [Figs. 1b, 1c, 5]. Therefore, the automatic tracing routine can run $24\mathrm{h}$ a day, because the queue of scanned image stacks is sequentially processed. The next step needing user interaction is the manual postprocessing of the final 3D reconstructions.

2.4.

Image Restoration

Diffraction as well as interference phenomena decrease the image quality and resolution of image stacks acquired through a TLB microscope. The resolution along the optical axis ( $z$ axis) is usually insufficient for an automatic tracing of faint axonal branches surrounded by more prominent neuronal structures.¹¹ Image restoration (deconvolution) can partly compensate for this limited resolution. Deconvolution is the most important step for the performance of the automatic tracing procedure.²⁹

The image formation is described by a convolution of the 3D intensity distribution within the object with a point spread function (PSF).³⁰ The PSF (Fig. 3 , see Color Plate 1) describes the optical properties of the imaging system. The convolution results in the 3D intensity distribution within the image. The intensity distribution within the object is then restored by deconvolving the recorded image stack using the PSF. Deconvolution is essentially an inversion of the image formation process. Thus deconvolution yields the best results provided the PSF can be determined exactly.

Fig. 3

Illustration of PSF modeling. (a) The $x\text{-}y$ view of the cone of light from a dendrite. Its intensity distribution is approximately symmetric with respect to the focal plane and the optical axis. (b) Modeled PSF obtained by treating the TLB microscope as if it would be a fluorescent microscope that suffers from no primary aberrations. These assumptions have been verified (Ref. 21). (High resolution image available online only. 3 10.1117/1.2815693.3 3 )

2.4.2.

Subdivision into bricks and linear deconvolution

The deconvolution is carried out by the HUYGENS software.³⁴ With the assumptions described above, the PSF is calculated in a wide-field fluorescence mode (Fig. 3), which is suggested by the HUYGENS software. An important parameter is the refractive index $\left(n\right)$ of the specimen. Aberration mea-surements²¹ yield a uniform value of $n_{\mathit{\text{specimen}}}=1.44$ for the specimen that is a $100\text{-}\mu \mathrm{m}$ -thick brain slice embedded in Mowiol $n_{\mathit{\text{mowiol}}}=1.49$ . Because the simplifications of the image formation process are justified, the imaging system is well described by the modeled PSF (Fig. 3).

This PSF is applied to the recorded image stack by a linear Tikhonov-Miller restoration filter.³⁵ It is derived from a least squares approach. This approach is based on minimizing the squared difference of the acquired image and a blurred estimate of the original object. The resultant image stack is improved significantly in terms of resolution and signal-to-noise ratio (Fig. 4 ). This deconvolved imaged stack is then of sufficient quality to apply an automatic tracing of neuronal structures (Figs. 5, 6, 7 ).

Fig. 4

Illustration of image restoration (deconvolution). (a) The $x\text{-}z$ view of various cones of light from dendrites and axons from one brain slice of $100\text{-}\mu \mathrm{m}$ thickness. The intensity is measured along the $z$ axis through a dendrite (white vertical line). (b) The $x\text{-}z$ view of the image stack after the application of a linear deconvolution filter. The longitudinal intensity is measured along the same line. (c) Intensity line plots before and after deconvolution. The signal-to-noise ratio is improved significantly, and the dendrite’s depth is decreased by a factor of more than 2.5. (one voxel is $0.5\mu \mathrm{m}$ wide.) (See Color Plate 1.) (High resolution image available online only. 4 10.1117/1.2815693.4 4)

Fig. 6

Image processing. (a) Maximum intensity $z$ projection of inverted image stack cropped from the mosaic image shown in Fig. 2b (box). Clearly visible is a fragmented axon with intensely filled boutons (swellings with likely synaptic contact sites) and a weak stain between them. (b) Maximum intensity $z$ projection of deconvolved, inverted image stack. The image quality is improved in terms of signal-to-noise ratio and resolution, especially along the $z$ direction. (c) Maximum $z$ projection after intermediate voxels (gray) being tested with the first part of the local threshold property function. (d) Maximum $z$ projection after intermediate voxels (gray) being tested with the second part of the local threshold property function. (e) Maximum $z$ projection after hit or miss transformation. Rectangular frame masks of increasing size (1-to-3-voxel radius) delete small and isolated artifacts. (f) Maximum $z$ projection after dilation and closing. A spherical structuring element tends to smooth the foreground objects and fills small gaps in the contour. (g) Maximum $z$ projection illustrating the end-line locations (white) that are local distance maxima from an inner seed compartment. These end-line locations will not be erased and, therefore, guarantee connectivity of objects during thinning. (h) Maximum $z$ projection illustrating the result of iterative layer-by-layer thinning. The end-line locations (white) remain unchanged. (i) Maximum $z$ projection illustrating the validation of the thinned graph. The shortest distance from one end-line location to all others will be calculated (gray value coded). (j) Maximum $z$ projection illustrating the validation of the thinned graph. Compartments and intercompartmental connections that are not part of a shortest connection (gray) will be erased. (k) Maximum $z$ projection illustrating the pruning of the thinned and validated graph. Because end-line locations were local intensity maxima, morphological swellings (boutons, spines) result in short branches that will be pruned if shorter than $3\mu \mathrm{m}$ . (l) Maximum $z$ projection after midline extraction. The prior fragmented axon is transformed to the thin less fragmented line representation. The scale bar in the lower left corner applies to the panels [(a),(a1) to ((l),(l1)]. Figures (a1) to (11) show maximum $x\text{-}y$ projections of the same image stack. Figures (a2) to (12) show enlargements to illustrate the fragmented boutonlike axon filling. The scale bar in the lower right corner applies to the panels (a2) to (l2). (High resolution image available online only. 6 10.1117/1.2815693.6 6 )

Fig. 7

Image processing. Illustration of the image processing pipeline for axonal branches. For detailed figure legends, see Fig. 6. (High resolution image available online only. 7 10.1117/1.2815693.7 7 )

Using the Tikhonov-Miller filter, the HUYGENS software is not capable of processing images larger than $2\mathrm{GB}$ at once. Therefore, the mosaic image stack is divided into smaller stacks of appropriate size (e.g., $3100\times 3100\times 120$ voxels), overlapping by 100 voxels, that can be processed by HUYGENS. These substacks, hereafter called bricks, are regrouped during the neuron tracing, resulting in a deconvolved and segmented mosaic image stack.

2.5.

Neuron Tracing

2.5.1.

Raster image–based segmentation

A major feature of image acquisition is the semiautomated setting of the exposure time of the CCD camera to control the image histogram as described previously. This feature is critical to the downstream neuron tracing algorithms. Their parameters were adjusted and systematically tested to derive the best possible tracings from such images. This issue of neuron tracing will be discussed later.

The down-sampled 3D image bricks are now subjected individually to segmentation algorithms (Figs. 5, 6, 7). The purpose is to separate voxels that represent neuronal structures from voxels that are part of the background.

Local threshold filtering

A target image with the dimensions of a deconvolved image brick is created and initialized with gray value zero. The voxels of the deconvolved bricks are then separated into three groups.

One group, below a lower threshold is set to background (gray value 0 in the deconvolved brick). A second group of voxels belonging to potential neuronal structures, with values above an upper threshold, is assigned as foreground (gray value 255 in the target image). The foreground consists of disjointed voxel regions that will be referred to as foreground objects. Voxels with values between the two thresholds form the third group of intermediate value voxels. The intermediate voxels are tested for local features defined by a local property function to decide whether they belong to foreground objects or the background.

The lower threshold is determined by calculating the intensity distribution of the deconvolved image brick. The deconvolution produces a thin unimodal histogram with the background clustered near 0, structures in the high range, and some remaining intermediate gray values. The lower threshold is taken as the histogram’s mean value plus 1.5 standard deviations (SDs). This gray value is usually between 1 and 15. It was derived after systematic testing and essentially deletes the background noise from unstained tissue.

The upper threshold assigns the prominent and most intense structures to the foreground. The value of the upper threshold is determined by calculating the histogram of the maximum $z$ projection of the deconvolved image brick. The upper threshold is then taken as the mean value plus 3.0 SDs. Systematic testing yielded that this is the best value to detect dendrites, well-filled axons, and the prominent parts of fragmented filled axons [Figs. 5c, 6c, 7c]. This upper threshold usually has a gray value between 30 and 60. The remaining voxels between the background (1 to 15) and the foreground (30 to 60) margins are referred to as the intermediate voxels.

The local property function comprises two steps. First, each intermediate voxel is set as the center of an $11\text{-}\times \text{-}11$ -voxel mask in the $x\text{-}y$ plane and the mean intensity of the voxels, regardless of group (background voxels were already set to 0), within this 2D neighborhood is calculated. If the centered intermediate voxel has a gray value that is larger than this mean intensity plus an epsilon value of 5 gray values, it is set to an intermediate gray level of 125 [Figs. 5c, 6c, 7c] in the target image. Otherwise it is set to background (gray value 0 in the target image). The optimal value of $ϵ$ was derived by systematic testing. The group of intermediate voxels consists usually of isolated artifacts or dim bridges between bright structures. These dim bridges are often found between axonal boutons (swellings of the axon that are likely sites of synaptic contact [Figs. 6(a2) to 6(e2)] and should therefore be part of the foreground. The resultant target image comprises three gray values: 0, 125, and 255 for the background, intermediate and foreground voxels, respectively.

The second part of the local property function inspects the intermediate voxels of the target image for connectivity to a foreground object. Systematic testing suggested that an intermediate voxel is set to the foreground if 10% of the voxels from a $17\times 17$ mask in the $x\text{-}y$ plane centered on this intermediate voxel are part of a foreground object, otherwise it is set to background [Figs. 5d, 6d, 7d].

Hit or miss transformation

The application of the above local threshold filter results in a binary image. This is then subjected to a hit or miss transformation³⁶ with rectangular frame masks of increasing size as structuring elements. The transformation is applied to every image plane. Isolated foreground objects that are completely surrounded by a frame are converted to background. Beginning with a radius of one voxel and increasing the frame size subsequently to three voxels, small, isolated artifacts that were introduced by the linear deconvolution are removed [Figs. 5e, 6e, 7e].

Dilation and closing

Next, dilation and closing filters³⁶ are applied. The dilation filter bridges gaps between the axonal boutons that have not been closed by the local threshold filter. Finally, a closing filter is applied [Figs. 5f, 6f, 7f]. Its geometrical interpretation is that a “sphere” rolls along the outside boundary of a foreground object within the image. Therefore, it tends to smooth sections of contours and fuses narrow breaks as well as long thin gulfs, eliminates small holes, and fills small gaps in the contour.³⁶ The 3D structuring element (sphere) for the closing and dilation has a radius of three voxels. Larger radii could result in a fusion of objects that should not be connected, and smaller radii would have no significant effect on the foreground objects. The radius of three, therefore, proved to be the best compromise after systematic testing.

Object labeling

The brickwise segmentation yields binary image stacks. Regrouping of these segmented bricks leads to a binary 3D image stack with the dimensions of a down-sampled mosaic image stack. This large 3D binary stack is subjected to an object labeling algorithm. Once a foreground voxel is detected, a region growing algorithm³⁶ fuses all connected foreground voxels to a subregion and assigns a unique integer label. The binary image is thus transformed to an image of $N$ individually labeled subregions, each representing a foreground object. Labeled objects are sorted according to their number of voxels. The largest foreground object is labeled as number 1, and the smallest object is assigned number $N$ . The background voxels are labeled as 0.

Raster to vector image conversion

Because the neuron occupies only a small fraction of the scanned image volume, the 3D raster image representation requires inappropriately large data storage. A more sophisticated vector data folder is created that stores only foreground objects but keeps the 3D information of the image. Figure 8 illustrates schematically the architecture of the new data folder. The voxels of each foreground object are replaced by vectors, hereafter called compartments. Each compartment stores the 3D coordinates of the prior voxel and could in the future store additional morphological information, such as local properties (e.g., radii, surface distance).

Fig. 8

Structure of the vector image. (a) Double linked list of objects. Each object represents an individual island of foreground subregions in the segmented raster image. The background is not stored. (b) Items of the object list are double linked lists of compartments. Each compartment represents one voxel of a foreground subregion (object) of the prior image. (c1) A compartment is realized as a vector. The 3D coordinates, as well as additional information, can be stored. (c2). The 3D topology of the image is maintained by intercompartmental smart pointers referencing the neighboring compartments and storing their relative positions (e.g., top left). These connections are used for navigating the vector image. (High resolution image available online only. 8 10.1117/1.2815693.8 8 )

Furthermore, each compartment is linked to its neighboring compartments by intercompartmental smart pointers, which allow access of neighboring compartments (e.g., top-left, bottom, left). Therefore, the implicit neighborhood representation by voxel coordinates of the 3D raster image is mapped to an explicit neighborhood construct in the compartment representation. The neighborhood pointers between the compartments preserve the 3D topology of any object. Furthermore, fast navigation through the object is possible using these intercompartmental smart pointers.

The compartments representing one object are grouped in a “double linked” list, which guarantees fast access to each compartment. Therefore, the 3D raster image of $N$ labeled foreground objects is converted to $N$ compartment lists. This transformation replaces the 3D raster image of labeled foreground objects by a double linked list, where each list item represents one of these disjointed foreground objects by a list of topology preserving compartments (vector image representation).

2.5.2.

Vector image–based midline extraction

Representing objects by their main structure or skeleton (approximate midline) is a commonly used technique in image processing. A fast and reliable way to calculate the skeleton of an object is thinning. Generally, thinning is a layer-by-layer (boundary voxel/compartment) erosion until only a unit-width skeleton is left.^{37, 38, 39} Therefore, three different classes of compartment topologies are defined: 6-adjacent $\left(N6\right)$ 18-adjacent $\left(N18\right)$ , and 26-adjacent $\left(N26\right)$ .³⁸ Two neighboring compartments are $N6$ if their Euclidean distance is equal to 1, $N18$ if their Euclidean distance is between (or equal to) 1 and the square root of 2, and $N26$ if their Euclidean distance is between (or equal to) 1 and the square root of 3³⁹ (Fig. 9 , see Color Plate 1).

Fig. 9

Neighborhood in 3D images, nonsimple and simple compartments. Panels (a) to (c) show a visualization of different neighborhood topologies in 3D images. Compartments are visualized as boxes centered on their position in image space. The central compartment is always visualized in red. All other compartments are presented in cyan. The boxes have an edge length of 1. In panels (a) to (c), the boxes are slightly moved away from the center box for better visualization. (a) The six adjacent neighborhood $\left(N6\right)$ around a central compartment. Compartments are adjacent if their boxes share one face. (b) The 18 adjacent neighborhood $\left(N18\right)$ . Compartments are adjacent if their boxes share at least one edge. (c) The 26 adjacent neighborhood $\left(N26\right)$ . Compartments are adjacent if their boxes share at least one point. Panels (d) to (f) show different possible configurations of compartments adjacent to each other. (d) This panel shows a central compartment, which is essential to preserve connectivity. If the red compartment is removed, the two cyan compartments are no longer connected. (e) As in (d), the red compartment is essential to preserve connectivity. Furthermore, this example illustrates the definition of an intersection compartment. (f) This panel shows a simple center compartment. If one removes the red compartment, the three remaining compartments are still connected. This is an example of a template, which is used to test during the thinning algorithm whether a compartment is removable, or not. (High resolution image available online only. 9 10.1117/1.2815693.9 9) Color Plate 1

A thinning algorithm has to obey the following demands that are derived from a 2D definition by Seul :⁴⁰ (1) Connected objects must thin to connected line structures. (2) The thinned lines should be minimally 26-connected. Two compartments are 26-connected if they are connected by a chain of $N26$ adjacent compartments. (3) Approximate end-line locations should be maintained. (4) The result of thinning should approximate the midlines of the structures. (5) Extraneous short branches introduced by thinning should be minimized.

To address these demands, we use the template matching algorithm described by Jonker.⁴¹ It is one of the most efficient thinning algorithms and can be briefly summarized as follows. First, a set of end-line locations is determined to represent the topology of the object. The end-line locations will not be removed. All other compartments, starting from the object boundary, are tested to check whether or not their removal affects the 26-connectivity of the object (Fig. 9). If compartment removal has no affect on the connectivity, this compartment is termed simple.

To determine whether a compartment is simple or not, its $N26$ neighborhood is compared with all possible $N26$ neighborhood templates that preserve local connectivity (Fig. 9). If the compartment neighborhood matches one of these templates, the removal of this compartment will not affect the local connectivity and the global connectivity will be preserved as well.⁴¹ Hence, this compartment is assigned to be simple and is removed. Whether or not a compartment is simple may change after deleting compartments. Therefore, the algorithm is iterative and finishes when no more compartments are removed (i.e., the thinning algorithm is an idempotent algorithm. The implementation of this thinning approach is described in detail below.

Detection of end-line locations

The end-line locations are determined by calculating the compartment, which is most distant to the object boundary (compartments with less than 26 neighbors). This is done by calculating the Euclidean distance map⁴² of the object. The compartment with the highest distance number is selected and called the seed point. In the case of several compartments sharing the highest distance value, one of them is selected randomly. Now the Euclidean graph distance³⁹ from each compartment in the object to the seed compartment is computed. This is the shortest connection along the graph between two compartments. Compartments with a local maximal graph distance are assigned as line endings [Figs. 5g, 6g, 7g].

Thinning

One disadvantage of thinning is the possibility of inward erosion, which could potentially create holes in the object. To prevent these artifacts, an iterative approach is used.³⁸

Using the 6-adjacent $\left(N6\right)$ and 18-adjacent $\left(N18\right)$ topologies defined above, boundary layers are assigned that are peeled by an iterative algorithm:

• 1. Collect the set of $N6$ boundary compartments (compartments where at least one of the $N6$ neighbors is missing).

• 2. Delete all simple compartments (according to the $N26$ topology) in this set, starting with those who have the lowest number of neighbors (i.e., intercompartmental pointers).

• 3. Collect the set of $N18$ boundary compartments (compartments where at least one of the $N18$ neighbors is missing).

• 4. Delete all simple compartments (according to the $N26$ topology) in this set, starting with those who have the lowest number of neighbors.

• 5. Repeat steps 1 to 4 until no more simple compartments are found [Figs. 5h, 6h, 7h].

  A detailed description of such an algorithm can be found in He ³⁸ and Jonker.⁴¹

Graph validation and pruning

Imaging and segmentation can lead to artifacts such as loops and clusters in the thinned midline representation of the objects, which have to be removed by postprocessing. Loops are removed by selecting one end-line location and calculating the shortest path from this end-line location to all other end-line locations in the object. Compartments or intercompartmental pointers that are not used by any path within this validation step [Figs. 5i, 5j, 6i, 6j, 7i, 7j] are removed from the object.

Because the approximate end-line locations were local distance maxima from an inner seed point, swellings within the neuronal branches will result in a short subbranch that we regard as artificial [Figs. 5k, 6k, 7k]. Therefore, these short branches are removed within this pruning step from the object if the distance from the according end-line location to the first intersection (Fig. 9) compartment is shorter than $3\mu \mathrm{m}$ (derived by systematic testing). Nonartificial axonal subbranches below this length threshold will hence be pruned as well. However, for most scientific problems, this error is negligible.

2.6.

Editing and Splicing

The image processing, so far, transformed the deconvolved raster image into a vector image representation of thinned foreground objects. Therefore, each object is represented as a 3D graph with end-line locations (endings). Furthermore, compartments that are connected to more than two neighbors will be set to be an intersection (Fig. 9). This set of endings, intersections, and normal compartments is converted to an ASCII (American Standard Code for Information Interchange) file that can be either visualized in AMIRA ⁴³ or NEUROLUCIDA. Up to this point, all processing steps are automated. Final editing and splicing of the reconstructions is done manually.

2.6.2.

Splicing

The results of this semiautomated reconstruction pipeline are 3D tracings of neuronal branches from adjacent, $100\text{-}\mu \mathrm{m}$ -thick sections. To obtain a complete reconstruction of a cell, the reconstructions from different sections have to be aligned and spliced (Fig. 10 ). Fortunately, large radial blood vessels running perpendicular to the pial surface retain their positions within each tangentially cut brain slice, rendering them ideal as position reference points. Therefore, the blood vessel pattern is extracted during the automatic tracing procedure. Blood vessels are detected via a region-growing algorithm in a maximum intensity $z$ projection of the original image stack, prior to the deconvolution. Voxels with values close to 0 (maximum transparent regions) are used as seed points for a region-growing algorithm. The region-growing algorithm seeks dark and almost circular shapes (shape number algorithm)³⁶ and labels these structures as blood vessels. According to the vessel pattern, the reconstructions are aligned coarsely. More accurate alignment is then done by matching the upper branch endings from a lower slice with the lower branch endings from an upper slice. Once all branches are aligned with their counterparts from the adjacent section, the nearest ending points are spliced (Fig. 10). Starting at the slice including the soma, progressively more sections will be connected resulting in a complete 3D cell reconstruction [Fig. 11a ]. If no soma is present (e.g., the cortical part of the thalamocortical axonal arbor), the splicing starts at the deepest (most distant from the pia surface) section [Fig. 11b].

Fig. 10

Splicing of adjacent brain slices. (a) Three-dimensional visualization of automatically reconstructed axonal branches located in two adjacent brain slices. A coarse alignment is made using the blood vessel pattern (circular profiles). Detailed alignment is made by connecting the branch endings. (b) and (c) Enlargements of sample branches from the two sections in (a). The fitting branch endings are spliced and ultimately result in a fully connected neuronal arbor. (High resolution image available online only. 10 10.1117/1.2815693.10 10 )

Fig. 11

Semiautomatic reconstruction of a pyramidal neuron in layers 2 and 3 and a thalamo-cortical axonal arbor. (a) Final version of an automatically reconstructed pyramidal neuron (in layers 2 and 3) in the rat barrel cortex. Twelve tangentially cut brain sections were scanned. After automatic image restoration and tracing, the user manually postprocessed the reconstruction by erasing artifacts and splicing axonal fragments within a section. After that, the 12 slices were aligned and the branch endings were spliced, resulting in the 3D graph representation of the neuron. (b) Final version of a cortical axon from a thalamocortical neuron. Twenty-eight tangentially cut brain sections were scanned and then automatically traced. Axonal branches were manually edited, aligned, and spliced to result in a fully connected 3D axonal arbor. The axon was reconstructed within $5\text{days}$ , comprising $22\mathrm{h}$ of human interaction. (High resolution image available online only. 11 10.1117/1.2815693.11 11 )

An international patent application is pending.⁴⁴

2.7.

Fractal Box Counting for Quantification of Relative Shape Deviation

For comparing the shape of automatically reconstructed branches with manual Camera Lucida–based reconstructions, a relative measure was used.⁴⁵ The aim was to obtain the percentage of relative deviation in shape at various resolutions. The shape comparison was applied to 2D $x\text{-}y$ projections of the traced branches. This reduction to projections is necessary because the embedding medium (Mowiol) may expand or contract over time depending on the ambient room temperature. After several hours at room temperature, the tissue tends to shrink in the $z$ direction (along the optical axis). After cooling the slide, this shrinkage is reversed (Table 1 ). Thus a comparison of the 3D morphology could be confounded by shrinkage and expansion of the tissue.

First, a bounding box is put around the $x\text{-}y$ projection of an individual branch [Fig. 12a ]. Second, this box is divided into eight subboxes, and the numbers of boxes that are intersected by the branch are counted [Fig. 12b]. In subsequent steps, each box is divided into eight boxes⁴⁵ until a box size between 200 and $300\mathrm{nm}$ is reached. This is done for both the manually and automatically traced branches.

Fig. 12

Fractal box counting method 1. (a) A bounding box around an automatically reconstructed axonal branch. The same box was put around the manual branch. (b) Each box is subsequently divided into eight subboxes and the number of structure intersected boxes is counted. In the illustrated case, six of eight boxes are intersected by the axon. (High resolution image available online only. 12 10.1117/1.2815693.12 12 )

If two reconstructions would match perfectly in shape, the number of boxes intersected by the branch would be the same within each dividing step. If one branch differs in shape from its counterpart, the number of intersected boxes will be different from a certain box size onward. We checked the difference in the number of intersected boxes at a box size (resolution) of $5\mu \mathrm{m}$ , $1\mu \mathrm{m}$ , $500\mathrm{nm}$ , and $300\mathrm{nm}$ (Fig. 13 ), which is nearing the physical resolution limit of the imaging system (Rayleigh Criterion).²⁵ The number of intersected boxes from the automatically reconstructed branch at the highest resolution $\left(300\mathrm{nm}\right)$ is set to be 100% of the branch. Differences in the number of intersected boxes will be given as percentages with respect to this value. This is illustrated in Fig. 14 , where the differences in structure-intersected boxes are shown for an increasingly cropped branch [Fig. 14(a1) to 14(a1d)].

Fig. 13

Results from the box counting method. The overlay of an automatically reconstructed branch (blue) and its manual counterpart (red) are analyzed at resolutions of $5\mu \mathrm{m}$ , $1\mu \mathrm{m}$ , $500\mathrm{nm}$ , and $300\mathrm{nm}$ . Identical boxes will be assigned a green color. When more than two braches match, then more green boxes are observed. At low resolution, the branches overlap almost entirely; whereas, at high resolutions, additional boxes within the automatic version increase. This depicts the quite sophisticated fractal box counting method. (High resolution image available online only. 13 10.1117/1.2815693.13 13 )

Fig. 14

Fractal box counting method 2. (a) An axonal branch reconstructed by the semiautomated system is shown (1). Progressively parts of the branch were erased, indicated in red [(1a) to (1d)]. The manual version of the branch is shown in (2). (b) The decreasing number of boxes becomes significant at a box size of approximately $10\mu \mathrm{m}$ . (c) Differences between the cropped branches and the original one (a1) in log-log plots. The difference in the number of boxes increases progressively. The number of boxes at the lowest resolution from the original branch is taken to be 100% of the structure. Percentages are plotted as horizontal lines. This gives an estimated measure for the relative deviation in shape between the original and the cropped branches at various resolutions. (High resolution image available online only. 14 10.1117/1.2815693.14 14) Color Plate 2

3.1.

Image Processing

An important constraint for the image processing is based on the exposure time setting of the CCD camera before the image acquisition process. The mean gray value in a typical field of view within the scan area is calculated and set to be 190. Therefore, the dynamic gray value range for neuronal structures is maximized and comparable within image stacks acquired from different brain slices. Furthermore, we avoid clipping of gray values,³⁴ which prevents artifacts from being introduced by deconvolution.

After acquisition of the mosaic image, the gray values are inverted and the large 3D mosaic image is subdivided into 3D bricks that can be processed by the HUYGENS software in combination with a linear Tikhonov-Miller filter. Because the deconvolution is performed on individual bricks, the gray value ranges for neuronal structure can be slightly shifted between the bricks. This is because HUYGENS estimates the background for each individual brick. The brickwise local threshold filtering compensates for this effect. Hence, the regrouped segmented mosaic image shows a uniform (not brickwise) distribution of foreground objects.

Connected subregions of foreground voxels within the segmented raster image are converted to lists of compartments with explicit intercompartmental neighborhood connections (vector image). The background is usually ${10}^3$ to ${10}^4$ times larger than the sum of the foreground voxels. Therefore, the storage size and the processing time decrease significantly after the data conversion from raster to vector image. Instead of processing gigabytes of data, the data set is reduced to a few hundred megabytes.

Template-based thinning reduces each object to its skeleton (approximate midline) [Figs. 5h, 6h, 7h]. The resulting graph is further validated and corrected until all thinning artifacts are erased. This correction usually creates smooth midline representations [Figs. 5l, 6l, 7l], neglecting morphological swellings such as dendritic spines or axonal boutons.

The tracing of poorly stained axons can result in a fragmented reconstruction of these structures (Figs. 5, 6, 7). Further small artifacts from unstained tissue, mainly astrocytes, are traced as well. Splicing of fragmented structures and erasing of artifacts during the manual postprocessing result in a 3D reconstruction of all neuronal branches within the brain section (Fig. 1).

Fig. 5

Image processing. (a) Maximum intensity $z$ projection of inverted image stack cropped from the mosaic image shown in Fig. 2b (box). Clearly visible is a fragmented axon with intensely filled boutons (swellings with likely synaptic contact sites) and weak stain between them. (b) Maximum intensity $z$ projection of deconvolved, inverted image stack. The image quality is improved in terms of signal-to-noise ratio and resolution, especially along the $z$ direction. (c) Maximum $z$ projection after intermediate voxels (gray) being tested with the first part of the local threshold property function. (d) Maximum $z$ projection after intermediate voxels (gray) being tested with the second part of the local threshold property function. (e) Maximum $z$ projection after hit or miss transformation. Rectangular frame masks of increasing size (1-to-3-voxel radius) delete small and isolated artifacts. (f) Maximum $z$ projection after dilation and closing. A spherical structuring element tends to smooth the foreground objects and fills small gaps in the contour. (g) Maximum $z$ projection illustrating the end-line locations (white) that are local distance maxima from an inner seed compartment. These end-line locations will not be erased and, therefore, guarantee connectivity of objects during thinning. (h) Maximum $z$ projection illustrating the result of iterative layer-by-layer thinning. The end-line locations (white) remain unchanged. (i) Maximum $z$ projection illustrating the validation of the thinned graph. The shortest distance from one end-line location to all others will be calculated (gray value coded). (j) Maximum $z$ projection illustrating the validation of the thinned graph. Compartments and intercompartmental connections that are not part of a shortest connection (gray) will be erased. (k) Maximum $z$ projection illustrating the pruning of the thinned and validated graph. Because end-line locations were local intensity maxima, morphological swellings (boutons, spines) result in short branches that will be pruned if shorter than $3\mu \mathrm{m}$ . (l) Maximum $z$ projection after midline extraction. The prior fragmented axon is transformed to the thin less fragmented line representation. (m) Maximum $z$ projection of final edited and spliced graph. The remaining background objects were manually erased and the axonal fragments were spliced. The scale bar in the lower left corner applies to all panels. (High resolution image available online only. 5 10.1117/1.2815693.5 5 )

3.2.

Comparison to State-of-the-Art Technique

The state-of-the-art technique for the 3D reconstruction of neuron morphology, especially of the largely extending axonal arbors, is the computer-aided interactive Camera Lucida technique (e.g., Neurolucida).¹³ Therefore, we compared the performance of the semiautomated reconstruction method to the manual reconstruction approach, as implemented with the Neurolucida system.

First, we checked whether manually traced branches were detected by the automated system and vice versa. We quantified whether branches were missed by either of the reconstruction approaches. In a second comparison step, the relative deviation in shape between manually and automatically traced neuronal structures was determined. Finally, the time frame for reconstructing complete cells with their entire axonal arbor was compared for the conventional method systems. The comparison was done for 54 slices containing 4 different filled cells (from 4 different rats). They were randomly chosen from a pool of already manually reconstructed cells. Each manual reconstruction was done by a different individual. Users manually traced a pyramidal neuron in layers 2 and 3, two star pyramids from layer 4, and one ventroposteromedial (VPM) thalamocortical axon.

3.2.1.

Number of identical branches

Figure 15 shows a comparison between an automated reconstruction of a section containing parts of a thalamocortical axon and its manually reconstructed counterpart. Example slices for the 4 neurons can be found in Figure 16 .

Fig. 15

Comparison of automated and manual tracing of axonal arbors. (a) Three-dimensional visualization of a tangentially cut brain slice containing a filled thalamocortical cell. The reconstruction is done by the semiautomated system. All branches are part of the axonal arbor. (b) Three-dimensional visualization of the manual reconstruction of the thalamocortical cell. The twofold enlargement shows that the two reconstructions are almost identical. Statistics from branch counting of 22 brain sections from this cell showed that the automated system traced 95% of the manually reconstructed branches. The manual tracer found 96% of its counterparts obtained by automated tracing. (High resolution image available online only. 15 10.1117/1.2815693.15 15)

Fig. 16

Comparison of brain slices for four different rats. The automatically traced brain slices can be seen in the left column, the manual ones in the right column. All scale bars are $100\text{-}\mu \mathrm{m}$ wide and apply to the left and right panels. (a) and (b) The $x\text{-}y$ projections of a slice from a thalamocortical cell. The comparison yielded that the performance of the manual and the automatic tracing were 96% identical. (c) and (d) The $x\text{-}y$ projections of a slice from a layer 4 cell. The comparison yielded that the performance of the manual system was worse than the automatic one. (e) and (f) The $x\text{-}y$ projections of a slice from a layer 4 cell. The comparison yielded that the performance of the manual system was worse than the automatic one. (g) and (h) The $x\text{-}y$ projections of a slice from a layer cell in layers 2 and 3. The comparison yielded that the performance of the manual system was much worse than the automatic one, due to the inexperience of the human tracer and some artifacts from unintentionally stained cells. (High resolution image available online only. 16 10.1117/1.2815693.16 16)

Each automatically reconstructed cell was inspected for missing braches that have been traced by the manual system and vice versa. The amount of manually reconstructed branches that were missed by the automatic system was similar for all four traced neurons. For the thalamocortical VPM cell, 5.4% were missing (Fig. 17 ), 3.5% and 3.4% for the layer 4 cells, and no manually reconstructed branch was missing for the cell in layers 2 and 3. In the case of branches that were missing in the manual tracing but found by the automated system, the result was less homogeneous. Approximately 4% of automatically reconstructed branches were missing in the manual version of the thalamocortical VPM cell, 29% and 47% were missed for the layer 4 cells, and almost 70% of the automatically traced branches were not found manually for the neuron in layers 2 and 3.

Fig. 17

Comparison of automated and manually traced neuron reconstruction. (a) and (b) Tangential view of 22 brain slices from a thalamocortical axonal arbor reconstructed automatically (a) and manually (b). (c) and (d) Coronal view of the 22 brain slices from the thalamocortical axon tree. Clearly visible are the cutting planes between the individual brain sections that need to be spliced. The manual tracing was done by an experienced person (T. K.). Comparison shows that the automatic system (c) performs as well as a well-trained human using a Camera Lucida system (d). (High resolution image available online only. 17 10.1117/1.2815693.17 17)

For the 4 cells, the automated method missed 29 of 682 manually reconstructed branches in 54 sections from 4 different rats. In other words, the new reconstruction technique found 96% of the manually reconstructed branches. The manual tracings were missing 268 of 923 automatically reconstructed branches from the same 54 slices. This means only 71% of the neuronal structures traced by the automated method could be found in the manual reconstructions. All of the missing branches, independent of whether they were missed by the manual or automatic system, were part of the axonal arbor. The dendritic branches were traced equally reliably by the two methods.

3.2.2.

Fractal box counting

The fractal box counting method was applied to 167 randomly picked branches from the 54 slices that were reconstructed by both methods. The frequency of the same percentage in the difference of the number of intersected boxes between an automatically traced branch and its manual counterpart is shown as a histogram (Fig. 18 , see Color Plate 1). For each of the four box sizes, a histogram of intersected boxes was made and fitted with a Gaussian curve. The mean values of these Gaussian fits yielded the average deviation in shape at various resolutions. The relative deviation in shape is almost 0% at the lowest resolution, meaning a box size of $5\mu \mathrm{m}$ . With increasing resolution (decreasing box size), the deviation increases to a value of around 2% at a box size of $1\mu \mathrm{m}$ and further to more than 3% at a resolution of $500\mathrm{nm}$ . At the highest measured resolution of $300\mathrm{nm}$ , the relative shape difference was around 5%. However, it can be seen that some branches deviate a lot more in shape than the mean values suggest. These large deviations were mainly observed for branches from dendritic trees.

Fig. 18

Relative shape deviation. A bounding box is put around the $x\text{-}y$ projection of an individual branch. This box is iteratively divided into eight subboxes and the numbers of boxes that are intersected by the branch are counted (Fig. 12). This is done until a box size between 200 and $300\mathrm{nm}$ is reached. The relative difference in the number of boxes that are intersected by automated and manually reconstructed branches is shown for various box sizes (resolution) (Fig. 13). The $y$ axis refers to the frequency of how often the difference between the boxes was of a certain percentage. For box sizes of $1\mu \mathrm{m}$ , most branches differ by approximately 2%. At the highest measured resolution of $300\mathrm{nm}$ , the relative deviation is approximately 5%. (High resolution image available online only. 18 10.1117/1.2815693.18 18)

4.1.

Recovery

The results showed that approximately 96% of the manually reconstructed axon branches were traced by the automated system. For each of the 4 cells, the performance of the new reconstruction technique is of a similar quality in terms of tracing branches that were found in their manual counterparts. The 4% of structure details that were not recovered fall into three groups.

The first group ( $<1\%$ in the sample of four cells) consists of faint axonal branches. An indicator for faintness is that these branches are not local intensity maxima. Therefore, the deconvolution algorithm tends to weaken these branches instead of amplifying their intensity.

The second group (less than 1%) comprises small axonal endings that are not distinguishable from reconstructed artifacts and, therefore, are erased by the user during manual editing and splicing. To minimize the loss of small axonal endings, we superimpose the maximum projection of the original image stack with the 3D reconstruction.

The third group comprises branches that are not located within the scanning area. By default, an area of $1\times 1{\mathrm{mm}}^2$ is scanned. For some sections, this area is not sufficient as the axons spread over a larger area. For these particular brain sections, the missing parts of the slice can be scanned. The branches missed initially are simply added to the rest of reconstructed branches, using NEUROLUCIDA for editing.

Summarizing these three groups, one can state that the 0 to 5% of missing automatically reconstructed branches consist either of small endings, faint branches, or not-scanned branches. By additionally scanning adjacent areas and inspecting the projection image, the number of missed branches can be minimized. Therefore, the performance of the semiautomatic reconstruction pipeline is limited by the first group, which implies the prerequisite that neuronal structures are local intensity maxima.

4.2.

Comparison with Manual Tracing

Only 71% of the automatically traced branches were found by the manual tracers. The number of missing structures varied widely between 4% and almost 70%. The missed branches can be grouped into two cases.

The first group consists of branches that were filled and stained but do not belong to the cell of interest. In some cases, surrounding neurons take up the tracer molecule, and some of their dendrites or their main axon will be traced.

The second group comprises neuronal structures that were missed by the human tracer. The user usually starts tracing at the soma. If an axonal branch ending is missed in that section, the human tracer tends to ignore these areas in the subsequent sections and, therefore, sometimes misses large parts of the axonal branching patterns. In some cases, this error becomes obvious during splicing and the human tracer has to go back to the tissue and to identify the connections missed previously. The size of the group of missing branches closely correlates to the experience of the user. Only training, double checking, and anatomical knowledge will minimize the loss of branches and improve the tracing quality.

Comparing the two groups, an explanation for the large deviation from 4 to 70% in missed branches can be given. The cortical part of the axonal arbor from the thalamocortical cell is far away from the biocytin injection site (soma) in the thalamus.¹³ Therefore, no artificial branches from other cells were found in cortex. Furthermore, this cell was traced manually by an experienced person in the laboratory (T. K.). Both facts result in the small number (4%) of missed branches.

In contrast, the manual reconstruction of one cell in layers 2 and 3 was missing almost 70% of the automatically reconstructed axonal branches. This probably has two reasons. First, layers 2 and 3 are closer to the site of cell filling. Therefore, some artificial branches from unintentionally filled cells were traced automatically and left out intentionally by the manual tracer. The second reason is that this cell was the first reconstruction of axonal morphology this individual ever attempted.

In conclusion, the manual reconstruction of neurons and their axonal arbors using a computer-aided Camera Lucida system (e.g., Neurolucida) depends strongly on the abilities of the individual tracer, and reconstructions will be somewhat subjective.⁴⁶ In contrast, the semiautomated tracing system delivers reliable and objective reconstructions even in the hands of new users (S. E.). S. E. was able to reconstruct the axonal arbor presented in Fig. 11b with no significant differences compared to a reconstruction of this neuron done by an experienced tracer (T. K.) (data not presented).

The automated system fails if the staining of the axon is insufficient, if the background is too high or in the presence of a strong gradient in thickness or contrast within the tissue. The latter will result in different optimal camera exposure times depending on the location of the field of view within the scanning pattern. Because the exposure is set to one value for the entire pattern, this may result in a failure of the subsequent deconvolution and neuron tracing. However, we are currently trying to compensate for this problem by applying of an exposure mask that compensates for inhomogeneous tissue.

Another limitation is that, due to the small diameters and tortuous paths of axons combined with physical resolution limits of the system, two close lying axons may occasionally be mistakenly connected by the automatic tracing algorithm. Because of the stereotyped nature of these errors, we may be able in the future to develop an algorithm to correct them. Nevertheless, because action potentials are propagated by saltatory conduction and the axon is surprisingly relatively electronically compact,⁴⁷ modelers and experimentalists are often unconcerned with the precise locations of all branch points. Rather, the focus is usually on the length and spatial organization of the axonal arbor. Thus, despite this limitation, our current system already measures the two main parameters of interest with high reliability and accuracy.

Indeed, the box counting method yielded small deviations in shape between manually and automatically reconstructed branches. This difference mainly arises from the larger number of points used by the automatic system to represent the graph. At low resolution, this difference is almost 0. At high resolutions $\left(300\mathrm{nm}\right)$ , close to the physical resolution limit, the relative deviation was approximately 5%. Therefore, depending on the scientific aim, these slight deviations can be irrelevant.

Nevertheless, Fig. 18 shows clearly that some branches, in particular at the highest resolution, deviate substantially from the 5% mean value given by the Gaussian fit. This group consists of dendritic branches. Because the dendrites are not as smooth as axonal branches, their automated reconstruction will also not be perfectly smooth. Reconstructed spine heads that were not removed during the graph validation and/or pruning, result in quite large relative shape deviations (Fig. 19 ). However, because the spines are part of the dendrites, the relative difference in shape can usually be neglected.

Fig. 19

Relative shape deviation for dendrites. (a) Maximum intensity projection of original image stack showing some dendritic branches. (b) Automatic tracing of dendritic branches. Spine fragments and thinning artifacts yield a rugged image. (c) Manual tracing of the same braches shown in (b). (See Color Plate 1.) (High resolution image available online only. 19 10.1117/1.2815693.19 19)

Table 1

Mowiol shrinkage. Qualitative measurement of the shrinkage effect of tissue embedded in Mowiol. After 3h at room temperature, the specimen showed shrinkage of 2to4μm along the optical axis. After a few hours at 4°C this process is reverted. By keeping the room air-conditioned and using the narrow illumination bandpass filter that protects the specimen from any infrared light, the specimen shrinkage is slowed down. However, once the specimen is at thermal equilibrium with the surroundings, the shrinkage is extremely slow that it can be regarded static over the scanning period.

4.3.

Reconstruction Time

An important issue is the reconstruction time needed to trace a neuron (e.g., the pyramidal neuron in layers 2 and 3 takes 10 to 30 brain slices). Because the automatic tracing is running $24\mathrm{h}$ a day and the human interaction time decreases from approximately $60\text{to}90\text{to}10\text{to}20\mathrm{h}$ , the reconstruction time decreases from weeks to days.

4.4.

Comparison to Other Methods

Reconstructing the branching pattern of neurons is a challenging task. In the literature, several methods for different imaging and processing approaches have been described. ^{1, 2, 3, 4, 5, 6, 7, 8, 9, 11} Most of the methods available today are based on confocal or two-photon imaging. This is due to the good signal-to-noise ratio and the high resolution along the optical axis for both of these imaging techniques. Rodriguez⁹ showed that laser scanning microscopy in combination with sophisticated software algorithms is even capable of resolving the neuronal morphology of dendrites at spine level. Their method is capable of reconstructing the dendritic arbor and visualizing the surface of the spines. By scanning several adjacent image stacks, the field of view is increased, but still their method lacks the possibility of reconstructing the axonal arbor. This is due to the limited field of view of laser scanning microscopes, the bleaching of the tissue,⁴⁸ and the lack of suitable mosaic systems for this imaging technique.

Another method that has been developed to analyze axonal arbors is the method described by Broser ¹⁵ Here lentiviral gene transfer is used to transfect a population of pyramidal neurons. The axonal arbors of these cells are subsequently stained. Two-dimensional mosaic microscopy in combination with extended focus imaging (EFI) is used to resolve the axons. Automatic 2D traces then quantify the axonal length in different areas of the tissue.

The appropriate reconstruction method strongly depends on the scientific task. To analyze spine shapes and distributions along a dendrite, the electron or confocal microscope methods described by Denk and Horstmann⁴⁹ and Rodriguez ⁹ might be the most suitable. For determination of statistical connectivity patterns, the 2D EFI method proposed by Broser ¹⁵ could be used. However, if the precise 3D extension of the axonal arbor is important, one can either use the time-consuming Camera Lucida technique or one can use the alternative method presented here.

4.5.

Future Improvements

We found that the human influence on the reconstruction during the splicing of neuronal branches from individual sections may result in ambiguous tracings. To further minimize the subjective issue of pattern recognition, we currently are developing a semiautomated interactive 3D editing and splicing tool (AMIRA extension: SPATIAL GRAPH)⁴³ in cooperation with the Konrad-Zuse Institute, which will replace the 2D editing environment of NEUROLUCIDA we used so far. The interface file format is the .hoc format (NEURON).^{50, 51} After editing and splicing in AMIRA, such a .hoc file can then be integrated into morphological databases or used for simulations in NEURON or NEURODUNE.⁵² Furthermore, an automated extraction of approximate dendritic radii will be added, because this is essential information for the simulation of neuron properties.

4.6.

Conclusions

We conclude that the semiautomated method of 3D tracing of neuronal morphology is a reliable, objective, and less time-consuming alternative to the current manual systems. Furthermore, the quality of reconstruction depends less on the abilities of the human tracer. Our method should facilitate the creation of large databases of neuron morphologies. This will yield interesting possibilities for simulation of information processing in small neuronal circuits. Furthermore, axonal outgrowth under various conditions of development and experience can now be more easily quantified in three-dimensions and statistically analyzed.