168	L. Domanski et al.

branches, their length, or the number of branching points. Branching structures typically possess a hierarchical organization, and it is important to capture and deliver information at each level. This higher order processing is described in Sect. 8.2.5.

There is renewed excitement nowadays around computational platforms such as Graphics Programming Units (GPUs) to speed up algorithms. This is particularly the case for high-throughput applications and for 3D and 4D datasets, where execution time still represents a major bottleneck. Section 8.3 describes the implementation of our linear feature detection algorithm on the GPU.

The technology presented in this contribution has been implemented in a WindowsTM package called HCA-Vision. For the user’s benefit, Sect. 8.4 outlines the software architecture and the main capabilities of this user-friendly tool (www.hca-vision.com).

Section 8.5 describes several representative applications in biological imaging. The first example demonstrates that even subtle phenotypic changes in the arbors of neurons in culture can be characterized if a sufficient number of neurons are analyzed. In this type of application, automation is central to reaching statistical significance. The importance of maintaining linear feature continuity is also highlighted by this example.

We have also used HCA-Vision with success to characterize the phenotype of astrocytes in response to drug-induced stress. In this application, fluorescently stained bundles of intermediate filaments were traced rather than neurites. Intermediate filaments are typical of the many polymers present within the cells, and we expect this application to grow in importance in the future.

In our final example, we describe how we used our sensitive linear feature detector to separate closely adjoining bacteria under minimal contrast, thus complementing the capability of more mainstream edge detectors.

Section 8.6 concludes this chapter by giving insight into recent developments, as well as by speculating about future directions.

8.2 Methods

There are many different algorithms suitable for linear feature detection. Good reviews are provided in [3,4]. Sun and Vallotton {Sun, 2009 #2}developed a fast linear feature detection method using multiple directional nonmaximum suppression (MDNMS). This approach is particularly intuitive and suitable for a nonspecialist audience. We will summarize the main steps of this algorithm in this section and describe its implementations on the GPU in Sect. 8.3.

8.2.1 Linear Feature Detection by MDNMS

Linear features are thin objects across which the image presents an intensity maximum in the direction of the largest variance, gradient, or surface curvature

8 High-Throughput Detection of Linear Features: Selected Applications...

169

Fig. 8.1 Illustration of linear windows at four different directions. “x” indicates the centre of the linear windows. The other pixels in the window are shown as small circles. The length of the linear window is 7 pixels in this illustration

(i.e., perpendicular to the linear feature). This direction may be obtained by the use of computationally expensive Hessian-based detectors or by using matched or steerable filters. Rather than searching for the local direction of linear features, we use 1D nonmaximum suppression (NMS) [5] in multiple directions to identify candidate pixels on a linear feature.

NMS is the process of removing all pixels whose intensity is not the largest within a certain local neighborhood, that is searching for local maxima. The shape of this local neighborhood is usually a square or rectangular window. For our purpose, we choose this local neighborhood as an orientated linear window. A number of approaches are available to obtain the directional local maximum within a onedimensional window. For example, a local maximum can be obtained by using basic morphological operators to check whether a pixel value in the input image is equal to that in the dilated image [6].

Figure 8.1 shows four examples of linear windows at angles equal to 0, 45, 90, and 135◦. Additional directions, such as 22.5, 67.5, 112.5, and 157.5◦ can also be used. NMS is performed successively for each direction at every pixel. The result is independent of the order of the scanning process. The longer the linear window, the lower the number of candidate pixels detected in an image. The ideal scenario occurs when the linear window crosses the linear feature at right angles as this will produce maximum contrast.

The outputs of the directional NMS are binary images, as opposed to grayscale images produced by most directional filter methods [3]. This eliminates the need for an arbitrary thresholding step. Linear features are detected by the union of MDNMS responses, that is we combine the NMS responses at each orientation.

The algorithms for linear feature detection using NMS can be easily extended for 3D images by using 3D linear windows. We used either 3 or 9 directions in 3D, although additional directions may be necessary, depending on the data.

170	L. Domanski et al.

Fig. 8.2 Cross section through a linear feature

8.2.2 Check Intensities Within 1D Window

Bona fide linear features are characterized by an approximately symmetric intensity profile across the feature, as opposed to edges which show a step-like profile. This translates into approximately equal intensity values around the central pixel on both sides of the linear window. Figure 8.2 illustrates the parameters that characterise the shape of a profile. Imax is the value of a local directional maximum at the centre pixel of the linear window. Iaverage1 and Iaverage2 are the average intensity values over the two sides of the local maximum within the local window; and Idiff1 and Idiff2 are the differences between the maximum value and these two average values. For a genuine linear feature, both Idiff1 and Idiff2 should be large. The parameter Idiff1 and Idiff2 can be used to control the sensitivity of the algorithm. Lowering the values of Idiff1 and Idiff2 increases the number of linear features detected.

8.2.3 Finding Features Next to Each Other

Some images present linear features in close proximity to one another. The NMS process will detect only one feature in each linear window. Figure 8.3 shows an example where four local maxima lie within the linear window. We can extend the NMS process so that multiple local maxima are detected for any particular linear window. Once a local maximum is found at the center of a linear window, a second local maximum search is performed within the original window using a smaller window size. The use of a smaller window size allows additional local maxima to be found within the original window. The newly found local maxima are ordered based on their maximum intensity as shown in Fig. 8.3. To prevent detecting false peaks due to noise, we also check that the intensity does not deviate significantly from Imax and that it conforms to the symmetry condition defined in Sect. 8.2.2. For example, in Fig. 8.3, Imax 3 and Imax 4 are two local maxima, but they are not strong enough to be retained.

8 High-Throughput Detection of Linear Features: Selected Applications...

171

Fig. 8.3 Multiple local maxima within a linear window

Fig. 8.4 Endpoint (red) and its neighborhood (green) for gap linking. The shortest path is shown in blue

8.2.4 Gap Linking for Linear Features

The union of the multiple responses of NMS at different directions generates many small objects, which may not belong to genuine linear features. We fit an ellipse to each object (set of connected pixels) and remove them if the major axis is too small. Alternatively, a simple pixel count can be used.

The continuity of the linear features may be broken due to noise or weakness of linear features. This results in small gaps in the combined responses of NMS. We restore continuity by joining linear feature endpoints to neighboring linear features through a shortest path. Here, a shortest path should be understood as a connection, which displays a high average intensity along its length. A standard technique called dynamic programming limits the combinatorial explosion that would otherwise be associated with the exploration of all candidate paths [7]. The computational overhead also remains small because the operation is only performed on small gaps – typically below 20 pixels in length. Figure 8.4 illustrates the process of gap linking from an endpoint to the skeleton in its neighborhood. Note that the domain shown in green in Fig. 8.4 is first transformed into polar coordinates prior to calculating the shortest path.