Abstract
Data sets imaged with modern electron microscopes can range from tens of terabytes to about one petabyte. Two new tools, Ssecrett and NeuroTrace, support interactive exploration and analysis of large-scale optical and electron-microscopy images to help scientists reconstruct complex neural circuits of the mammalian nervous system.
An adult human brain has more than 100 billion neurons and 250 trillion neural connections. Connectomics is a research field to develop methods—from data acquisition to analysis—for reconstructing a nervous system’s detailed, complex neural connections.1 Further complicating this study, neuronal cells can be extremely small and closely packed: dendritic spines with diameters of 50 nanometers and synaptic clefts of 20 nm. High-resolution imaging plays a central role in identifying and resolving nanometer-scale neuronal cells. Optical and electron microscopes are common neuroscience tools for acquiring high-resolution images. A typical optical microscope can attain resolutions of 200 nm per pixel, whereas a modern electron microscope (EM) can attain resolutions up to 5 nm per pixel.
High-resolution imaging has opened the door to reconstructing detailed neural connections, but existing methods can’t effectively manage and process these large data sets. For example, with a 5 nm resolution in the x-y plane and a 30 nm slice thickness, the EM scan of a 1 mm3 brain tissue sample is about one petabyte of raw data. Automated scanning devices, such as Harvard University’s Atlum (Automatic Tape-Collecting Lathe Ultramicrotome), can produce raw image data at a rate of about one terabyte per day. So, one immediate problem is how to efficiently store and retrieve such huge data sets. Data storage should provide a reliable mechanism to both store raw data streams and enable data access in arbitrary locations with minimal latency.
Another challenge is processing and manipulating these large data sets to extract scientifically meaningful information in a reasonable processing time. The processing includes image filtering, segmentation, and visualization. Many techniques that work well on relatively low-resolution data modalities, such as computerized tomography (CT) and magnetic resonance imaging (MRI), just won’t work on large-scale EM data sets. For example, existing vascular-segmentation techniques work reasonably well for axon segmentation on optical-microscopy images but might not apply directly to feature-rich, high-resolution EM images. Similarly, many deformable image-registration methods developed for CT and MRI might not work well for aligning large stacks of optical- or electron-microscopy data.
Many well-known image-processing tools and libraries can’t handle extremely large data sets at interactive rates, and most existing tools and algorithms require fitting the entire input data into the main memory for processing. In addition, the common segmentation practice in optical and electron microscopy is a time-consuming, laborious manual process that becomes a major workflow bottleneck as data sizes increase.
Automated Tools for Interactive Workflows
We’ve developed two software tools to cope with ever-increasing data sizes and to give neuroscience researchers more flexible, interactive workflows.2
The first, Ssecrett (Serial Section Reconstruction and Tracing Tool), is a volume-exploration and manual axon-tracing tool for extremely large-scale neuroscience data sets based on 2D slices. It has a client-server architecture that places the data set on the server side and lets the client request an arbitrary 2D cross-section view of it through a socket connection. This simple solution follows the standard multiresolution approach of dynamically loading data on demand. Because users can view only a small amount of data at any one time, the Ssecrett server loads data only as it’s needed.
The second tool, NeuroTrace, focuses on interactive segmentation and high-quality 3D visualization of high-resolution EM data sets. It combines local 2D set segmentation with 3D tracking to extract the 3D tubular structure of neural processes without the burden of processing the entire input data set. Our visualization method uses an on-demand technique for local image filtering and volume rendering to achieve interactive performance. An out-of-core multiresolution octree data structure supports large-scale volume rendering, and a batch-processing framework enables offline segmentation of multiple neural processes. We also use the massively-parallel-processing capability of modern GPUs to accelerate many time-consuming and computationally demanding tasks.
Ssecrett and NeuroTrace fit nicely into the unified, semiautomatic neural-process segmentation and visualization workflow. Ssecrett provides an intuitive, interactive method for quickly inspecting large input data without requiring high computing power on the client side. It lets users find the region of interest with minimal time and effort. Once the user has selected the region of interest, NeuroTrace can extract the detailed 3D geometry of a neural process and visualize it using high-quality 3D rendering techniques.
For related work, see the sidebar on large-volume data visualization and processing.
Ssecrett
As data increases to petabyte sizes, copying or moving it in its entirety becomes difficult or even impossible. To facilitate large-volume data sharing, our prototype client-server architecture supports distributed connectome databases and peer-to-peer data sharing. Connectome volumes reside on the computer that gives the most convenient access and best network performance. Clients can share links and view data without ever needing to know where it actually resides.
Preprocess Dicing
One challenge in handling large data sets is formatting the input volume to make any part or resolution quickly and easily accessible. Our microscope systems initially store the volumes as a series of large TIFF images. Extracting an arbitrary subregion from a TIFF file is difficult, because neighboring pixels along the z-axis are in different files. Furthermore, to produce a high-level zoomed-out view, applications must load almost all the data.
Dicing is a preprocessing stage we developed to transform the input volume into small 3D blocks. The blocks are subsampled recursively and combined to create a hierarchy of resolutions. This structure guarantees that the server can load a subset of blocks for any view, whether high or low magnification. The availability of low-resolution blocks also supports progressive updates that the server can access and load quickly. Thus, if a client application has to wait before loading a high-resolution file, Ssecrett can first render a low-resolution image to keep the application responsive.
Still, our current dicing implementation remains time-consuming and I/O bound. For example, dicing an 8-Gbyte raw file takes about four hours. Most of the time is spent reading the images, compressing blocks, and writing files. We can reduce these times by running the algorithm in parallel on a distributed file system, but dicing a large volume would still require a significant amount of time. To view new stacks interactively, we extended the prototype so that it could view multiple volumes simultaneously. A subset of the volume (a stack of images from just-scanned slices) can be relatively small, so the dicer can process it quickly. Once diced, the subset can be visualized as part of the larger, already-processed volume. This feature also allows the system to position a new stack and warp it independently, which in turn supports integration of data alignment and stitching on the boundaries with the Ssecrett system.
The Ssecrett Server
We wrote the Ssecrett server in C++, making extensive use of the Visualization Toolkit (VTK). We created a vocabulary of fewer than 40 messages for communication between the client and server. Examples include messages such as LOAD_FILE (ask the server to load the indicated file), REQUEST (a specific image request), and IMAGE (a response from the server to the client indicating that the requested image follows). The server can run on Unix/Linux, Mac, or Windows systems.
The Ssecrett system supports many simultaneous users by starting a dedicated server for each user. This is the most stable approach but not the most memory efficient. If multiple users are viewing the same volume, each server will load its own copy of the data. Fortunately, the architecture doesn’t preclude more memory-efficient options. A single server can be easily extended to connect to multiple clients. Another option takes advantage of modern multicore processors by having multiple servers share a single cache that uses system shared memory. Either approach could help collaboration between multiple remote users who are viewing and discussing the same volume.
The server can run locally (automatically started up when the client is started up) or remotely, as we discussed earlier. At Harvard, we typically have a server launcher running on several different platforms (with access to different diced data sets), waiting for a client to connect and request a server on that machine. The server launcher acknowledges the client’s request (confirms receipt of the request), starts a server process, confirms the server is running, and passes the connection information to the client. The client can then connect directly to the server to request a data set, and the data-set visualization or analysis commences. Meanwhile, the server launcher returns to its listening state, waiting for another client to connect.
The Ssecrett client can display arbitrary off-axis resampled images. To minimize the data sent via the socket, the server resamples the volume and sends the client just compressed 2D images. This approach avoids sending larger 3D blocks. Additionally, the server doesn’t read or send more information than the client will display for a given request. If the client requests an 800 × 600-pixel image of the entire top slice of a volume, which might be 20k × 20k pixels, the server will read only low-resolution blocks from disk, delivering higher resolution only when the client can use the information.
The prototype can also compress the images with JPEG to make them smaller for transmission. The computational cost of compressing images is significant, so compression is skipped when the client is local to the server and network bandwidth isn’t an issue. With the client in Albany, New York, and the server in Boston, we can get about 5 frames per second in an 800 × 600 viewing window with only 165 Kbytes-per-second bandwidth. The frame rate varies considerably with the viewing window’s size—decreasing as the window size increases. By comparison, viewing of off-axis images results in negligible degradation in frame rate.
Ssecrett Client
We created a simple C++ prototype client viewer, built using FLTK (Fast Light Toolkit, a cross-platform C++ GUI toolkit) and VTK. We envision a variety of clients connecting to the server to request blocks of the volume. Ssecrett manages communication between the client and server via a vocabulary that lets any independent segmentation application easily access the central server. We regularly run the client on both Mac and Windows systems.
Our prototype client provides a simple manual tracing capability, whereby the user drops “bread crumbs” in the volume via the mouse or keyboard. These bread crumbs are connected to give a track or skeleton tree representing axons or other structures. Tracks can have multiple branches, merge with other tracks, and include notes associated with the track in general or specific locations along it. Furthermore, the notes are easily searchable, so one user might label data-set regions as “weird” or “check this out” and save the information to a database managing information about the data set. Or the user could just email another user the tracks file, a small compact format describing just the tracks and notes. The other user could then pull up previous researchers’ work, search for “weird,” and quickly iterate through the areas in question. Furthermore, a note acts much like a bookmark, saving the camera orientation when the note is created so that returning to a note gives the user the option of viewing it from the same perspective.
The client provides several view options, including a 2D-slice view (see Figure 1a), a 2D-slice view that includes a small 3D view of the slice within the entire volume (see Figure 1b), and a view that gives three orthogonal slice views in addition to a 3D view indicating the slice locations within the volume—similar to the MRI and CT viewing-software formats (see Figure 1c). In the near future, we plan to implement the client using Qt instead of FLTK, which will give the client a more polished feel as well as greater flexibility and functionality in the interface. This will also enable direct integration of the Ssecrett client into the NeuroTrace GUI. Like FLTK, Qt is a cross-platform C++ GUI toolkit.
Figure 1.
Three viewing options in the Ssecrett client: (a) a 2D slice view, (b) a 2D slice view with a 3D view of the slice within the volume, and (c) three orthogonal slice views with a 3D view. Users can interact with the data in multiple formats, which faciliates understanding a slice’s location within a full volume
Ssecrett Prototype Test Results
We’ve tested the Ssecrett prototype on several confocal and EM images. Figure 2a shows a full view of an EM slice, and Figure 2b shows a zoomed-in view of the full-view area outlined in red**). The Ssecrett system can progressively transmit images of coarse to fine resolution per request and, on our local network, zoom in and out through image stacks at interactive rates. We’ve also traced multiple myelinated axons in EM images using Ssecrett’s manual tracing (see Figure 3). Each color segment represents the cross section of the tube created according to the bread crumbs dropped by the user, denoting the segmented axon’s center axis. Figure 1c shows a similar axon-tracing result on a stack of confocal microscope images.
Figure 2.
The Ssecrett system’s zoom capability: (a) a full view of a slice and (b) a zoomed-in view of the region marked as a red rectangle in (a). Ssecrett displays the full-slice image quickly by loading large, low-resolution blocks. Similarly, it displays the zoomed-in view quickly because it loads only small blocks (containing the full-resolution data) for the area of interest.
Figure 3.
Manually traced axons using the Ssecrett system. The system can display the traced axons in 2D or 3D and can associate notes with each trace at specific locations. Notes also act as bookmarks to recall the view at the time the note is assigned.
To test the Ssecrett prototype’s scalability, we simulated a petabyte volume by duplicating an 8-Gbyte volume 50 times along each axis. Although the files weren’t actually duplicated on disc, the prototype read them as if they were separate files. As we expected, the system performed as well on this test as it did on smaller volumes. It’s possible that the operating system cached files, making the read process more efficient, but using off-axis sampling reduced this effect’s significance.
The prototype uncovered two areas for improving system performance. First is reducing the time to load and generate the highest-resolution image for a new slice. When no blocks are in cache, it can take up to 10 seconds to load all the blocks necessary to generate the full-resolution slice. We’re considering several ways to deal with this issue. We could use smaller leaf blocks, reducing the volume of data we need to load to generate an image. We could also implement a more sophisticated caching algorithm that prefetches blocks and anticipates future requests. Finally, we could run a distributed database in parallel to exploit parallel disk I/O.
The second performance issue is the slow rendering for very large viewing windows. Viewing full slices (as shown in Figure 2) with an 800 × 600 window takes a little over 5 frames per second, but it drops to closer to 2 fps for a 1,000 × 1,000 window (for a 500 × 500 window, it’s 8 fps). We can improve performance increasing the JPEG compression ratio, but this introduces obvious artifacts. An alternative solution is to use a video compression algorithm such as MPEG4 to get better compression from correlation between successive frames.
The biggest challenge we’ve encountered so far is efficiently accessing the original volume to reformat into the multiresolution diced files. The original data was arranged in a series of 20k × 20k TIFF images. We processed slabs of 50 images at a time. This approach is reasonably efficient but breaks down when images are too large. When the size of the 50 images overwhelms the memory available on the dicing computer, we switch to dicing 50 × 50 rows of pixels. This approach addresses memory limitations but requires opening and closing each image file many times, so it takes too long to process large volumes. The delay is especially pronounced when the raw input files are on a remote file system. The combination of opening files multiple times and a one-pixel overlap in the cubes increases the dicing time so that it’s much longer than a single read of the data. Processing multiple subvolumes helps, but we’re working toward an image-based dicing strategy that reads and processes images one at a time.
NeuroTrace
The NeuroTrace tool for 3D volume visualization includes user-controllable, semiautomatic segmentation. Our approach maintains system scalability to large data sets and high-performance parallel computing architectures. We described an earlier version of NeuroTrace elsewhere.2 Here, we focus on its recently added out-of-core volume-rendering and batch-processing capabilities.
The NeuroTrace Framework
We implemented NeuroTrace in C++, OpenGL, Nvidia CUDA (Compute Unified Device Architecture), and Qt. We use CMake as a cross-platform make tool for easily porting the existing application to different operating systems. We’ve now built NeuroTrace for Windows and Linux.
The NeuroTrace framework consists of four main conceptual modules: the core, segmentation module, visualization module, and GUI.
The core
This module manages data handling and storage, including out-of-core data loading, octree generation, and cache management. We use a separate loader thread for opening and loading data sets. When the thread has finished loading a new data block, the ray caster adds it interactively to the visualization. This streaming approach lets us speed up the volume renderer’s initial startup time by loading and displaying low-resolution blocks first, while the loader thread in the background still reads in the higher-resolution blocks. By encapsulating all data-access functions into our data-management layer, we can switch from local data storage to a network-based file system without any changes in the segmentation and visualization modules.
The segmentation module
This module accesses the original data via the core module. The segmentation results are propagated back and passed on to the interactive visualization system. The NeuroTrace framework treats this module as a black box.
The visualization module
This module handles 3D visualization, including direct EM data volume rendering, on-demand denoising and edge detection, and implicit surface ray casting of segmented processes. It accesses the data via the core module and provides functions such as 2D-slice and 3D rendering at interactive frame rates, transfer functions, and clipping planes. For binary segmentation information, this module also offers two-level volume rendering for segmentation masks.
The GUI
We used the Qt framework to develop the user interface, which is completely decoupled from the core, segmentation, and visualization modules. This lets users switch from one GUI framework to another if necessary. The flexible layout uses dockable views that users can adjust individually, depending on their primary requirements. A transfer function editor supports both 1D and 2D transfer functions, and adjustable 2D-slice and 3D views let users activate or hide views, according to their preferences.
Additionally, we keep two tables of all changeable GUI visualization parameters, one on the user-interface side and one at the renderer side. This lets us synchronize the two tables at an arbitrary time in a thread-safe way. Using a separate thread for the user interface keeps the application responsive, even during time-consuming processing tasks. We plan to extend NeuroTrace to a client-server application, further decoupling the user interaction and processing.
Neural-Process Segmentation
NeuroTrace’s segmentation method combines 2D neuron cross-section segmentation and 3D center-point tracking to reconstruct a 3D tubular surface of neural processes. This approach uses only a small fixed-size 2D subregion of the 3D input data, so it scales well to arbitrarily large data sets. We can treat the subregion retrieval process as a black box, integrated into the NeuroTrace core as a high-level API, without knowing the actual implementation details. For example, we can retrieve a 2D region through the network if the data is on a server.
For 2D neuron cross-section segmentation, we use the modified-active-ribbon model originally proposed by Amelio Vázquez-Reina and his colleagues.3 The model consists of two deformable interfaces (that is, level sets) that push or pull against each other according to various internal and external forces until they converge to the cell membrane’s inner and outer wall, respectively, as shown by ϕ1 and ϕ2 in Figure 4. This makes the method more robust in noisy and feature-rich EM images than other level-set approaches. Some of the forces we use are
the data-dependent force to move the interfaces toward the membrane boundary,
the ribbon-consistency force to keep the distance constant between two level-set interfaces,
the curvature force to maintain the smoothness of the interfaces, and
the image-correspondence force to robustly initialize the interface locations on subsequent slices.
Figure 5 shows four stages in the deformation of the active ribbon.
Figure 4.
The active-ribbon model, consisting of two deformable interfaces (ϕ1 and ϕ2) interacting with each other. Because the two interfaces shrink-wrap the cell membrane between them, the proposed model is more robust to noise and incomplete structures.
Figure 5.
Neural cell membrane segmentation using the active ribbon. From left to right: a user-drawn initial active ribbon (green) shrinks to the correct membrane boundary. Inner and outer level-set curves are moving according to the force fields defined on the image domain and the curves until they converge to the cell membrane.
To segment arbitrarily oriented neural processes in 3D space, we implemented a 3D tracking algorithm that follows the neural process’s centerline. The main assumptions behind our 3D tracking method are that neural processes are topologically equivalent to 3D tubes and that the center of each tube cross section lies on the 3D tube’s centerline. So, after we find the membrane of the cell cross section using active-ribbon segmentation, we can trace the neural process’s centerline in 3D space by connecting the center point of each cell cross section.
We propose a two-step method of estimation and correction. In the estimation step, the tangent direction Vt at the last center point is computed using a one-sided finite-difference method. We also keep the previous tracking direction Vp. The new tracking direction is then the weighted average between those two vectors: V = ωVp + (1 − ω) Vt (Figure 6a). The weight ω controls the amount of history used to determine the current tracking direction, providing smoother transitions between slices. Once we have a new tracking direction, we can use simple extrapolation to estimate a temporary new center position of the next slice. We then project the local frame of the previous slice onto the new plane defined by the temporary new center and the tracking direction. A new 2D region of interest is sampled from the input data for the segmentation process.
Figure 6.
A 3D segmentation. To create this segmentation result, we combine 3D centerline tracking with the 2D active-ribbon segmentation: (a) weighted extrapolation for neuron cross-section, center-point tracking and (b) a traced axon using NeuroTrace. The green rings are the segmented active ribbons in 2D slices; the red line is the traced axon’s centerline.
In the correction step, we replace the temporary center by the correct center of the segmented neural membrane. Figure 6b shows a real example of 3D centerline tracking and axon segmentation using our 3D centerline tracking method.
To speed the time-consuming active-ribbon-segmentation process, we implemented the GPU multiphase-level-set solver using Nvidia CUDA. The solver employs an iterative, adaptive-level-set partial differential equation (PDE) solver using a block-based update scheme. The implementation minimizes the communication between the CPU and GPU and reuses the data on the GPU by safely updating each active block multiple times without updating the active block list. To maintain the correct ribbon topology, we need to recompute the signed distance field for each level set when the active block list is updated. To compute image correspondence force, we use a GPU-accelerated deformable-image-registration method to find a 2D vector field that maps one image to another. More detail on this method is available elsewhere.2
Batch and Parallel Processing
Because a data set includes a huge number of neural-processes and our segmentation method can run automatically without user intervention, we implemented batch processing for segmenting multiple neural processes. This way, users can draw multiple initial level sets, one per neural process, and let them run sequentially. Once the batch processing ends, the user can examine the result and make any necessary edits to rerun the segmentation and improve the result.
To implement batch processing, we treat each neural-process segmentation as an independent process that keeps track of the current status variables and dynamically allocates memory to store the segmentation result. Our segmentation method stores each segmentation process on disk as a cache and subsequently loads only the current segmentation process into main memory. Therefore, users can launch many segmentation processes concurrently as long as each segmentation result can fit in main memory. The batch processing switches between segmentation processes and performing the segmentations sequentially.
Another benefit of storing each segmentation process is the ability to save the entire workspace and restore it later. The enormous data size and numerous neural connections can require days or even weeks to complete a segmentation of neural processes in a reasonably large input data set. To be able to continue the segmentation work for a longer period of time, NeuroTrace provides a way to save the snapshot of the current segmentation on disk. The approach saves all cached segmentation processes and some additional information, such as a pointer to the input data, the labeled volume, and the current segmentation process index.
To minimize the disk and memory usage, we convert two floating-point images of the level-set results into a simple but compact representation (see Figure 7a). The main idea is to create a binary image that each pixel value represents either inside or outside the active ribbon. Then we store only the list of indices for the inside pixels and later restore a full binary image by marking the corresponding pixels as inside. From the recovered binary image, we can restore floating-point level-set images by assuming that the level set passes through the center of the zero-crossing edge (see Figure 7b). This approach can cause a slight shift from the original zero-level set, but it doesn’t make a big difference in a realistic setup because we’re interested only in pixel-level segmentation results.
Figure 7.
A compact representation of the active ribbon. (a) Two floating-point images (left) are converted into a binary image (middle); only the pixels in two level sets—that is, on the active ribbon—are stored as a list of indices (right). (b) A zero-level set defined by a signed distance image (left) is converted into a binary image (middle) and then restored to a floating-point distance image, which results in a slight shift from the original zero-level set location (right).
To run multiple segmentation processes in parallel, we can distribute processes across different GPUs. The GPUs handle the assigned processes in parallel. Each GPU accesses the main data server to retrieve the 2D subregion image for running the level-set segmentation, but it doesn’t need to store the entire input data set. So, our parallel segmentation process can easily scale to various hardware configurations and data sizes, such as multiple GPUs in a shared-memory system for small and medium data sets or a GPU cluster system using a fast network for extremely large data sets. The implementation of such a client-server heterogeneous cluster system is a subject for future research.
Volume Visualization
NeuroTrace’s volume visualization system is based on CUDA, Nvidia’s C-like programming environment for GPU programming. We implemented single-pass ray casting,4 in which we use OpenGL for rasterization and CUDA kernels for the actual ray casting. CUDA is a versatile platform for scientific computing because it enables full utilization of the huge parallel programming power of today’s graphics hardware. This lets us render the volumetric data set, display the ongoing segmentation, and perform on-the-fly filtering on the volume’s visible parts, while still achieving interactive frame rates. For large data visualization, we integrated support for octree rendering into our ray caster.
NeuroTrace visualizes high-resolution EM data, which is extremely dense, heavily textured, and low in signal-to-noise ratio. To visualize the complex structure of interconnected nerve cells, we added support for on-demand filtering and edge detection of a volume’s visible parts. Additionally, we can display the neural-process-segmentation results by combining direct volume rendering with implicit surface ray casting in a single rendering pass. To reduce the memory footprint, we compress the segmentation results.
On-demand filtering and edge detection
We designed our approach to on-demand filtering and edge detection to explore a volume quickly and highlight regions of interest (such as myelinated axons). Figures 8a, b, and c show some results in an example EM data set. For noise removal, we implemented Gaussian, mean, nonlinear-median, and bilateral and anisotropic diffusion filters in 2D and 3D with variable neighborhood sizes. For EM data, we found bilateral and diffusion filtering to be the most useful.
Figure 8.
Visualization of electron microscopy (EM) volumes to depict neural processes. Our approach combines on-the-fly filtering (smoothing, edge detection) with volume ray casting: (a) on-demand edge detection in an EM data set, (b) opacity modulation based on the local-histogram-based edge metric, (c) denoising and edge enhancement, and (d) 3D visualization of segmented neural processes embedded in the volume. Together, these techniques provide a 3D context for all processing stages, combining 3D volume data with the implicit surfaces resulting from segmentation.
To modulate the EM data’s opacity during volume rendering, we use a local histogram-based edge metric.2 To find edges, we calculate the difference of two histograms in a voxel’s neighborhood several times along different directions.5 A large difference between histograms is a good indicator of a large brightness change and, therefore, an edge.
To cope with large data sets, we based our algorithm on a hierarchical blocked-volume representation (an octree) and a dynamic caching system in GPU and CPU memory. We start by creating a list of all visible blocks for the current viewpoint. In the next step, a CUDA kernel performs the filtering on all visible blocks and saves the results directly in a GPU memory cache. During a final ray-casting step, the values resulting from filtering and edge detection are combined with the original data values to highlight important regions and enhance the volume boundaries.
On the other hand, we use the inverted edge values to suppress less important tissue and more homogeneous regions. To store the on-the-fly computed blocks from the noise-removal or edge-detection pass, we allocate two caches directly on the GPU. After the algorithm determines all visible blocks for the current viewpoint, our cache management system automatically decides which blocks actually need to be computed. Blocks that are already in the cache don’t need recomputation. Blocks that have already been computed but aren’t visible in the current viewpoint are kept in the cache for later reuse. If cache memory gets low, unused blocks are flushed from the cache in least-recently-used (LRU) order.
Visualization of segmented neural processes
We visualize the original EM volume data together with the segmented neural processes in a combined volume rendering. Our system displays the segmented axons or dendrites as semitransparent isosurfaces, which the CUDA-based ray caster renders in the same ray-casting pass as the original data’s direct volume rendering. Keeping our large data set sizes in mind, we store the segmentation results in a compact format that represents each axon as a simple list of elliptical cross sections.
During standard front-to-back ray casting, we compute implicit surfaces on the fly from the set of ellipses. At each sample along the ray, we test for potential intersections with isosurfaces. If an intersection is detected, the isosurface’s color and opacity are composited with the previously accumulated volume-rendered part, and volume rendering continues behind the isosurface. If the user interactively assists the neural-process segmentation, the axons are continuously updated in the 3D view. This lets the user stop and adjust the segmentation process whenever necessary. Figure 8d displays several segmented axons in a 3D volume rendering.
Out-of-core data management
NeuroTrace’s out-of-core data management system comprises a three-level memory hierarchy (see Figure 9).
Figure 9.
A three-layer cache hierarchy. The GPU cache holds at least the current set of blocks required for ray casting; the CPU cache contains a larger superset of the blocks currently on the GPU; and the octree cache stores the entire hierarchy in local or network storage. This approach provides two different sets of working blocks in GPU and CPU memory. The CPU cache enables fast paging of blocks into the smaller working sets of the GPU cache at PCI Express transfer rates, while itself holding a larger working set to help hide the much-higher disk- or network-access latency.
On the first level, the GPU cache manager maintains the GPU cache, which consists of a large 3D cache texture that’s always resident in GPU memory and contains all the volume blocks that are currently needed for ray casting. The cache is managed using an LRU scheme, and the ray caster can freely mix and match volume-octree blocks of different resolutions.
The CPU manager maintains the second-level CPU cache, which consists of a large pool of blocks that constitute a superset of the blocks resident in the GPU cache. Blocks in the CPU cache are flushed only when the cache runs out of space; otherwise, the blocks are kept as long as possible to enable fast re-downloading of blocks that have been flushed from the GPU cache. Again, an LRU scheme governs cache management. The CPU cache manager can also prefetch blocks from disk to make them promptly available to the GPU cache manager when requested.
The octree manager maintains the third-level out-of-core octree cache, which consists of all octree nodes comprising the entire volume and its multiresolution hierarchy. All octree nodes at all resolution levels are stored on disk in a huge binary file. The octree manager creates the initial octree cache and loads specific blocks on demand, as requested by the CPU cache manager. Loading of requested blocks proceeds in the background using a separate thread that lets the ray caster proceed while more data are being loaded. Octree nodes are currently loaded from local disk storage, but the memory hierarchy we use is conceptually independent of whether blocks are loaded from disk or over the network. We’re extending the octree manager to request and load blocks on demand over a network connection. All other system components of the system will be completely unaware of the load location of the octree nodes.
NeuroTrace User Study and Discussion
To evaluate NeuroTrace’s usability and accuracy, we conducted an informal user study.2 Both novice and expert users liked its ease of use and accurate automatic segmentation. Batch processing has since added greater flexibility to the workflow. Figure 8d shows an example axon segmentation result.
We achieve interactive frame rates for our on-demand filtering and edge-detection implementation. However, we would like to extend our on-demand filtering framework. We’re working on a general multiresolution histogram framework, which would let us interactively compare the histograms of different bricks in different resolutions to further enhance our opacity-modulated EM data rendering.
The large z-slice distance of high-resolution EM data causes it to exhibit high anisotropy, which poses a challenge to high-quality data visualization. Using multiresolution approaches, however, we can alleviate this problem by first downsampling along only the volume’s x and y dimensions. After a few downsampling steps, the volume becomes isotropic, and we can proceed by downsampling in all three dimensions.
Figure 10 shows a top-down view of a volume consisting of 43 Gbytes of raw data. The entire octree cache for this volume consists of a single file of roughly 75 Gbytes, which includes all resolution levels and duplicated neighbor voxels between adjacent blocks to avoid interpolation artifacts at block boundaries. In contrast to viewing only 2D cross sections, viewing the entire volume reveals 3D structures more easily. Figure 10 also clearly shows the aligned layering of individual scanned EM slices that constitute the volume. So, only an interior region of the whole scan defines a full volume, whereas the borders show varying regions of the originally unaligned slices after alignment. This is an unavoidable side effect of the acquisition process. We can easily cut out the fully defined volume using axis-aligned clipping planes. However, we usually enable the entire unclipped volume, which lets users view all acquired image data and visually inspect alignment accuracy.
Figure 10.
An interactive overview using semitransparent volume rendering of a large hippocampus data set comprising about 43 Gbytes of raw data (resolution 14,176 × 10,592 × 308 pixels). The process that properly aligns the original scanned 2D slices places them in correct relation to each other in 3D, which shows up clearly at the volume’s borders, where varying regions are missing in different slices.
Volume-rendering performance using the data set and transfer function depicted in Figure 10 is between 10 and 30 fps on an Nvidia Geforce GTX 285, for a viewport resolution of 800 × 600 pixels. Table 1 gives performance results for different zoomed-in, top-down views, which are rendered using the octree levels that yield a roughly one-to-one mapping of rendered voxels to displayed pixels. The table also lists what percentage of the displayed octree level and the entire tree, respectively, are actually in the view frustum.
Table 1.
Rendering performance of top-down views for the 43-Gbyte hippocampus volume.
| Level | Frames per second | No. of 323 blocks | Visible percentage of octree level (%) |
Visible percentage of entire octree (%) |
|---|---|---|---|---|
| 0 | 10 | 443 × 331 × 10 | −0.5 | −0.4 |
| 1 | 19 | 222 × 166 × 10 | 0.4–1.1 | 0.1–0.2 |
| 2 | 26 | 111 × 83 × 10 | 2.0–2.3 | 0.1 |
| 3 | 30 | 56 × 42 × 5 | 11–55 | 0.1–0.3 |
| 4 | 30 | 28 × 21 × 3 | 76–100 | 0.1 |
Using a one-to-one mapping of octree voxels to output-image pixels effectively bounds the memory needed for any given view. For example, in a top-down view of level 0 of the octree, less than 0.5 percent of the nodes in this level are actually in the view frustum. The visible percentage is very low for all zoomed-in views and increases only when the view is zoomed out so much that most or all of the volume becomes visible. The volume is subdivided into blocks of 323 voxels. For the 14,176 × 10,592 × 308 hippocampus volume, the octree comprises 10 levels, from a single 323 block for the entire volume at the tree root (level 9) to 443 × 331 × 10 blocks for the highest resolution (level 0).
Although we introduce Ssecrett and NeuroTrace as independent software tools in this article, we envision them merging into a single program in the near future. Our first step in that direction is to investigate a common data API based on NeuroTrace’s current out-of-core data management system. This API will work on different systems and arbitrary data scales seamlessly. We’re also working on high-quality 3D volume visualization on a GPU cluster system using a distributed ray-casting method. Our long-term goal is to develop a real-time, large-scale neuroscience data acquisition and processing system.
Related Work in Large-Volume Data Visualization and Processing
The most widely used neural-image-processing tool is Reconstruct.1 It has a simple 2D-slice viewer with basic paintbrush editing functions to manually draw boundaries of neuronal cells in both optical- and electron-microscopy images. Ju Lu and his colleagues extended Reconstruct with a region-growing method for semiautomatic segmentation.2 Some other open source and commercial software packages, such as the NeuronJ plug-in for ImageJ (http://rsbweb.nih.gov/ij), Bitplane’s Imaris (www.bitplane.com), and Amira (www.amiravis.com), provide automatic segmentation and tracking. However, none of them can handle arbitrarily large data sizes with complex 3D neural structures or data retrieval over the network, as Ssecrett does (see the main article).
Unlike optical microscope images, electron microscope (EM) images have little research addressing their processing. Elizabeth Jurrus and her colleagues3 and Jakob Macke and his colleagues4 proposed automated slice-to-slice segmentation algorithms using deformable models on EM images. However, they didn’t discuss their methods’ scalability to large data sets, and their methods apply only to certain EM image types, such as serial block-face scans.
Traditional research in scalable volume visualization and processing has focused mainly on multiresolution approaches and parallel volume rendering on clusters. Eric LaMar and his colleagues were among the first to introduce a hierarchical (octree) bricking scheme for GPU-based volume rendering.5 Recently, Enrico Gobbetti and his colleagues introduced single-pass octree ray casting on the GPU.6 At each frame, they encode a compact octree indexing structure into a small texture and perform stackless octree ray casting. Christoph Müller and his colleagues introduced a GPU-based cluster system for volume rendering of large data.7 The system employs a sort-last rendering scheme that first subdivides the volume into a kd-tree and then distributes it to the cluster nodes, which perform GPU-based ray casting on the individual bricks. However, most previous approaches deal with data in the gigabyte range, which is several orders of magnitude lower than our data requirements.
References
J.C. Fiala, “Reconstruct: A Free Editor for Serial Section Microscopy,” J. Microscopy, vol. 218, no. 1, 2005, pp. 52–61.
J. Lu, J.C. Fiala, and J. W. Lichtman, “Semi-Automated Reconstruction of Neural Processes from Large Numbers of Fluorescence Images,” PLoS ONE, vol. 4, no. 5, 2009, pp. e5655.
E. Jurrus et al., “Axon Tracking in Serial Block-Face Scanning Electron Microscopy,” Medical Image Analysis, vol. 13, no. 1, 2009, pp. 180–188.
J.H. Macke et al., “Contour-Propagation Algorithms for Semi-automated Reconstruction of Neural Processes,” J. Neuroscience Methods, vol. 167, no. 2, 2008, pp. 349–357.
E. LaMar, B. Hamann, and K. Joy, “Multiresolution Techniques for Interactive Texture-Based Volume Visualization,” Proc. IEEE Visualization 99, ACM Press, 1999, pp. 355–362.
E. Gobbetti, F. Marton, and J. Guitan, “A Single-Pass GPU Ray Casting Framework for Interactive Out-of-Core Rendering of Massive Volumetric Datasets,” Visual Computer, vol. 24, nos. 7–9, 2008, pp. 797–806.
C. Müller, M. Strengert, and T. Ertl, “Optimized Volume Raycasting for Graphics-Hardware-Based Cluster Systems,” Eurographics Symp. Parallel Graphics and Visualization (EGPGV 06), Eurographics Assoc., 2006, pp. 59–66.
Acknowledgments
This work was supported in part by the US National Science Foundation under grant PHY-0835713, the Austrian Research Promotion Agency FFG, Vienna Science and Technology Fund WWTF, the Harvard Initiative in Innovative Computing, the US National Institutes of Health under grant P41-RR12553-10 and U54-EB005149, Microsoft Research, and Nvidia.
Biographies
Won-Ki Jeong is a research scientist at Harvard University’s Initiative in Innovative Computing and its School of Engineering and Applied Science. His research interests include image processing, scientific visualization, and general-purpose GPU computing in biomedical image analysis. Jeong has a PhD in computer science from the University of Utah. He received an Nvidia Fellowship in 2007. He’s a member of the ACM. Contact him at wkjeong@seas.harvard.edu.
Johanna Beyer is a postdoctoral fellow at King Abdullah University of Science and Technology (KAUST) and a research scientist at the VRVis Research Center. Her research interests include large-data visualization, parallel visualization, and GPU-based volume rendering for neuroscience and neurobiology. Beyer received her PhD in computer science from the Vienna University of Technology. Contact her at johanna.beyer@kaust.edu.sa.
Markus Hadwiger is an assistant professor of computer science at King Abdullah University of Science and Technology (KAUST) and the principal investigator of the WWTF Scale-VS research project on the scalability and confluence of scientific visualization and interactive segmentation at the VRVis Research Center. His research interests are petascale visual computing, scientific visualization, volume rendering, and general GPU techniques. Hadwiger has a PhD in computer science from the Vienna University of Technology. He’s a coauthor of Real-Time Volume Graphics (A.K. Peters, 2006). Contact him at markus.hadwiger@kaust.edu.sa.
Rusty Blue is a research scientist at Kitware. His research interests include real-time reconstruction algorithms for 3D electrical-impedance imaging, parallel-programming techniques, and high-performance-computing platforms. Blue has a PhD in biomedical engineering from Rensselaer Polytechnic Institute. Contact him at rusty.blue@kitware.com.
Charles Law is vice president of Kitware. His research interest is in parallel visualization of large data. Law has a PhD in neuroscience from Brown University. Contact him at charles.law@kitware.com.
Amelio Vázquez-Reina is a PhD candidate in computer science at Tufts University and a researcher at Harvard University’s Initiative in Innovative Computing and Center for Brain Science. His research interests are statistical machine learning, computer vision, and biomedical imaging. Vázquez-Reina has an MSc in electrical and computer engineering from the University of Alcala, Spain. He’s a student member of IEEE. Contact him at amelio@seas.harvard.edu.
R. Clay Reid is a professor of neurobiology at Harvard Medical School. His laboratory studies information processing in the mammalian visual cortex with a combination of physiological and anatomical approaches. His current research employs laser-scanning microscopy to image neural activity and electron microscopy for large-scale (multiterabyte) 3D reconstructions of neural circuits. Reid has an MD from Cornell University Medical College and a PhD in Biophysics from the Rockefeller University. Contact him at clay_reid@hms.harvard.edu.
Jeff Lichtman is a professor of molecular and cellular biology and a member of the Center for Brain Science at Harvard University. He’s a developmental neurobiologist whose research interests focus on ways in which information about the world become instantiated in the neural circuits of developing animals. Lichtman has a PhD in neurobiology from the Washington University School of Medicine, St. Louis. Contact him at jeff@mcb.harvard.edu.
Hanspeter Pfister is the Gordon McKay Professor of the Practice in Harvard University’s School of Engineering and Applied Sciences. His research interests are at the intersection of visualization, computer graphics, and computer vision. Pfister has a PhD in computer science from the State University of New York at Stony Brook. He’s a senior member of the IEEE Computer Society and member of ACM, ACM Siggraph, and Eurographics. Contact him at pfister@seas.harvard.edu.
Footnotes
Selected CS articles and columns are also available for free at http://ComputingNow.computer.org.
Contributor Information
Won-Ki Jeong, Harvard University.
Johanna Beyer, King Abdullah University of Science and Technology.
Markus Hadwiger, King Abdullah University of Science and Technology.
Rusty Blue, Kitware.
Charles Law, Kitware.
Amelio Vázquez-Reina, Harvard University.
R. Clay Reid, Harvard University.
Jeff Lichtman, Harvard University.
Hanspeter Pfister, Harvard University.
References
- 1.Sporns O, Tononi G, Kötter R. The Human Connectome: A Structural Description of the Human Brain. PLoS Computational Biology. 2005;vol. 1(no. 4) doi: 10.1371/journal.pcbi.0010042. doi:10.1371/journal.pcbi.0010042. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 2.Jeong W-K, et al. Scalable and Interactive Segmentation and Visualization of Neural Processes in EM Datasets. IEEE Trans. Visualization and Computer Graphics. 2009;vol. 15(no. 6):1505–1514. doi: 10.1109/TVCG.2009.178. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 3.Vázquez-Reina A, Miller E, Pfister H. Multiphase Geometric Couplings for the Segmentation of Neural Processes. Proc. IEEE Conf. Computer Vision and Pattern Recognition; IEEE CS Press; 2009. pp. 2020–2027. (CVPR 09) [Google Scholar]
- 4.Krüger J, Westermann R. Proc. IEEE Visualization 03. IEEE CS Press; 2003. Acceleration Techniques for GPU-Based Volume Rendering; pp. 287–292. [Google Scholar]
- 5.Martin D, Fowlkes C, Malik J. Learning to Detect Natural Image Boundaries Using Local Brightness, Color, and Texture Cues. IEEE Trans. Pattern Analysis and Machine Intelligence. 2004;vol. 26(no. 1):530–549. doi: 10.1109/TPAMI.2004.1273918. [DOI] [PubMed] [Google Scholar]