Automation of Hessian-Based Tubularity Measure Response Function in 3D Biomedical Images

Dzyubak, Oleksandr P.; Ritman, Erik L.

doi:https://doi.org/10.1155/2011/920401

International Journal of Biomedical Imaging

On this page

Abstract Introduction Background Methods Results Discussion Acknowledgments References Copyright Related Articles

Special Issue

Modern Mathematics in Biomedical Imaging

View this Special Issue

Research Article | Open Access

Volume 2011 | Article ID 920401 | https://doi.org/10.1155/2011/920401

Automation of Hessian-Based Tubularity Measure Response Function in 3D Biomedical Images

Oleksandr P. Dzyubak¹and Erik L. Ritman¹

Academic Editor: Yangbo Ye

Received04 Aug 2010

Revised12 Oct 2010

Accepted10 Dec 2010

Published22 Feb 2011

Abstract

The blood vessels and nerve trees consist of tubular objects interconnected into a complex tree- or web-like structure that has a range of structural scale 5 m diameter capillaries to 3 cm aorta. This large-scale range presents two major problems; one is just making the measurements, and the other is the exponential increase of component numbers with decreasing scale. With the remarkable increase in the volume imaged by, and resolution of, modern day 3D imagers, it is almost impossible to make manual tracking of the complex multiscale parameters from those large image data sets. In addition, the manual tracking is quite subjective and unreliable. We propose a solution for automation of an adaptive nonsupervised system for tracking tubular objects based on multiscale framework and use of Hessian-based object shape detector incorporating National Library of Medicine Insight Segmentation and Registration Toolkit (ITK) image processing libraries.

1. Introduction

Humans (indeed all biological multicellular organisms) are made of multiscale hierarchy of structures ranging from subcellular structures (10⁻⁷ m) to cells (10⁻⁵ m), basic functional units (the smallest aggregation of diverse cells that behaves like the parent organ 10^-4 m), organs (10⁻²–10⁻¹ m), and bodies (10⁰ m). In addition to the local scale variation, biological structures are also characterized by shape. For example, the blood vessels are tubular objects interconnected into a complex network and have a range of structural scale (5 μm diameter capillaries to 3 cm aorta). This large-scale range presents two major problems; one is just making the measurements, and the other is the exponential increase of component numbers with decreasing scale. Quantitative analysis of such systems, for example, blood vessel trees and networks of neurological dendrites and axons, seem to be best measured from 3D images of those structures. With the remarkable increase in the volume imaged by, and resolution of, modern day imagers, the practical problem now is the extraction of this multiscale data from those large, detailed image data sets. Thus, there is a need for an automatic tool which mines data across both space and scale to capture local information about the objects which describes the feature and now becomes associated with its appropriate position and scale.

Without prior information for a scale description of the image content, an image has to be studied at all scales. The basis for design of an automatic tool for such description could be derived from human perception model [1]. The human eye comprises a large number of individual receptors (over 150 million rods and cones). Such “imager” has no prior information about input, and therefore, it is designed to extract the information by applying sampling apertures at a wide range of sizes simultaneously [1, 2]. Since the information from a single individual sensor is almost meaningless, the sampling should be done not by the individual rods and cones in a human eye or detectors in imagers but by the sensor neighborhoods. Such sensor neighborhoods implement fundamental “multiscale perception” at different scales simultaneously but at this point no memory or analysis is involved yet. The very first level of analysis starts when grouping and storing the information from the local neighborhoods into meaningful sets (so-called “blurring signal”) and establishing interconnections between neighborhoods. Since a priori size of the features in the object (signal) is unknown, the “blurring” is used by both humans and machines for adaptive multiscale representation of meaningful signal features of a different size. Within the next step when we extract an object, certain properties have been already attributed to it, and the rest of the information is now considered as nonobject and often called “noise” or “background”. As both noise and object are always parts of the perception process, there is no way to separate the noise from the object if models of the object and of the noise are absent. To distinguish noise from object, we need to model the image content. That could be achieved by using certain mathematical operators, “feature detectors”, interacting with the data followed by analysis of the results of “feature extraction”. Since objects in images comprise features with varying size, it is quite logical to perform feature extraction at the scale which matches the local feature size, “its native scale”. A recipe for automatic feature extraction in multiscale framework can be given as follows: (i) build multiscale representation by smoothing the image at each scale; (ii) choose an appropriate “feature detector” to compute local structural properties (e.g., gradients, curvatures, flows); (iii) compute “local extrema” of a “feature detector response” function; (iv) find the strongest local response of the structural properties which is considered as feature identification.

2. Background and Principle

To build the multiscale representation of the image, the proper aperture (windowing) function as an operator should be chosen. Formalism for the scale-space representation was introduced by Witkin [3] and further developed by Koenderink [4]. The idea of this approach is to generate a one parameter family of smoothed images obtained by convolving the original image with a Gaussian kernel of size In case of Gaussian kernels, the kernel size is called variance, which is related to standard deviation as . The parameter in this family represents the scale at which the finer image structures are still “perceptible” whereas the spatial structures with size smaller than will be smoothed out as shown in Figure 1. As pointed out by Koenderink [4] and Hummel [5], this family of smoothed images may also be derived as the solutions of the diffusion equation If is constant, it reduces to the isotropic diffusion equation , and the linear spatial scale-space representation can be generated using Gaussian (continuous) or Binomial (discrete) kernels [2–4, 6–10]. Data sampled in the temporal domain (e.g., the movie frames are samples taken at regular intervals) can also be scaled in similar to the spatial domain fashion. To treat a multiscale context over the temporal domain, it was suggested to use the Poisson kernels [11–13]. Perona and Malik [14] extended the scale-space concept to nonlinear scale-spaces based on nonlinear diffusion formulation with nonconstant diffusion coefficient ; (see (2)). Comprehensive overview of nonlinear scale-spaces based on parabolic partial differential equations could be found in [15, 16]. For such nonlinear cases, the Bessel scale-space can be built [17].

(a)

(b)

(c)

(d)

The scale spaces could be generated by using various kernel functions. Recently, wavelets and their applications to signal and image processing have attracted attention of the scientists in many fields. A very good collection of papers on the wavelet theory and its applications can be found in the book by Heil and Walnut [18]. The relationship between the wavelets and the scale spaces was demonstrated by Mallat and Hwang [19]. In [20–22], methods for generating scale-space representations based on wavelets were suggested, and the results of application of wavelets were very promising. As an extension of the Gaussian and wavelet approaches, Wang and Lee [23] proposed the scale-space representation derived from B-splines.

Even though different kernel functions have been proposed to generate the scale-space representations, the Gaussian kernels remain the best candidates so far [1–10]. Such “uniqueness” of the Gaussian kernels was shown by Babaud et al. [24] and is based on a priori scale-space constraint formulated by Witkin [3], Koenderink [4], and Yuille and Poggio [6]: “No new feature points (no “spurious detail”) should be created with increasing scale”. Florack et al. [25, 26] extended constraints and formulated the mathematical requirements known as “axioms for an uncommitted visual front-end” [27] which should be satisfied for the systems without any a priori knowledge about inputs: (i) linearity: no “feedback” from the system; (ii) spatial shift invariance: no preferred spatial or temporal location; (iii) isotropy: no preferred spatial or temporal orientation; (iv) scale invariance: no preferred size or scale of the aperture.

Another motivation to use the Gaussian scale spaces has some support from neurophysiological and psychological experiments which has shown that receptive field profiles in the mammalian retina and visual cortex can be well modeled by sums of Gaussian [28] and Gaussian derivative components [29–31].

All these properties make the linear Gaussian scale spaces the best choice for development of the automatic unsupervised systems for multiscale signal and image analysis, when there is no in advanced information available concerning preferable scales. The recipe which allows adaptively choosing the proper local scale parameter at every geometrical location was suggested by Lindeberg [32, 33]. This recipe comprises a two-step procedure: (i) convolving the original image with a Gaussian kernel of size ; (see (1)) and (ii) analyzing a response function from derivatives or some (possibly nonlinear) combination of derivatives of that convolution [32] where and are orders of derivatives.

The strongest response of such function (with respect to σ) over scales then indicates the proper Gaussian scale probe (Gaussian observation kernel) with width corresponding to object feature size has been found [34]. Due to the commutative properties of the convolution and taking derivative operations, the order of operations in the procedure above could be changed to convolving the original image with the operators constructed from the derivatives of Gaussian (so called “Gaussian derivatives”). The distinguishing characteristic of such operators is a combination of the opposing properties, localization, and optimal response to noise [35–37]. To demonstrate this principle, we modeled the original image with an object with Gaussian intensity profile of kernel size and convolved this image with the first- (Figure 2) and second- (Figure 3) order Gaussian derivatives with kernel sizes varying over the range (in the figures only are shown). Then, we measured intensity values of response functions (see (3)), across the object and extracted the strongest responses for all scales of the Gaussian probe. The strongest responses to both first and second Gaussian derivatives with kernel probe varied over the range and convolved with the original image with feature size are shown in Figure 4. These plots demonstrate the property of the response functions such that it has the strongest response when Gaussian probe size approaches the object size.

(a)

(b)

(c)

(d)

(a)

(b)

(c)

(d)

(a)

(b)

Unfortunately, the amplitude of Gaussian derivative operators tends to decrease with increasing scale due to the fact that with increasing scale, the response is increasingly smoothed. This gives more preference to smaller scales. To compensate such increase and, thus, improve accuracy of the automatic scale selection, Lindeberg [10, 32, 34] suggested using the so-called -parameterized normalized derivatives This method of scale selection allows feature detectors to find such points in the image that the -normalized operator response has an extremum with respect to both position and scale.

In the case of tubular-like structures in an image, ridge detection with automatic scale selection can be done using a second derivative of Gaussian kernel function [38, 39]. Depending on values of the -parameter, the detected object features can be quite different. Analyzing the influence of the γ-parameter on feature detection with automatic scale selection, Lorenz et al. [40] chose to be 1.5 which worked well for variety of intensity line profiles (e.g., Gaussian, bar-, triangle-like). However, for elongated structures with bar-like intensity profiles (such intensity profiles can be found in high-quality imagers with narrow point spread functions [41, 42] or imagers that use deconvolution preprocessing algorithms [43, 44]), the ridge detector creates false responses at small scales (basically these are “edge responses at small scales”) as depicted in Figure 5. The line intensity profiles demonstrating the problem are shown in Figure 6. Since in automatic approaches there is no the preferred scale (all scales should be treated equally) in an image, a number of solutions have been proposed to avoid or suppress these false responses in scale space. Koller et al. [45] suggested applying a nonlinear operator that combines the response of two edge-detectors on both sides of a hypothetical ridge. Lorenz et al. [40] used an edge-indicator to suppress the response to edges. Lindeberg [34] used a hybrid approach taking the useful properties from both the scale-space height ridge and the second derivative scale-space ridge. While studying the problem of the influence of the -parameter on feature detection with automatic scale selection, Majer [38, 39] derived the g-normalization parameter value from the statistical approach based on a white noise sampling model. In these studies, he concluded that the ridges generated by a second-order Gaussian derivative operator do not suffer from the false responses to edges if the value is used.

(a)

(b)

(c)

(d)

(a)

(b)

Early vision perception occurs at all scale simultaneously and can be modeled by generating the image scale-space representation that introduces an additional variable, spatial scale size [1–4]. At this point, such early vision system is fully ignorant of the geometry. As soon as the local scales are established, the early analysis of an image starts with analysis of intensity variations and directions by means of spatial derivatives to reveal local image structure. For example, for elongated objects, the derivative value along the object is close to zero whereas the derivatives across the object are large negative values and their ratio is close to unity. Following the scale-space ideas developed so far, a complete hierarchical set of scaled differential operators has to be used [25]. The Gaussian derivative operators described earlier in this paper constitute the natural differential operations on a given scale. Thus, a set of Gaussian derivatives and their combinations could be used for very complex object models and analysis. Some well-known combinations of the high-order derivatives have special names like Hessian, Laplacian, and so forth and are used to build special functions for identification of certain shape patterns in images. Such mathematical functions model human and machine perception and are ultimately used in unsupervised object tracking systems [1, 2, 35, 36]. For this purpose, the Hessian-based multiscale object enhancement filters were developed [4, 25, 46–48]. In these anisotropic filters, the pixel intensity transformation is locally governed by the “objectness measure” functions which are built using the combination of local Hessian eigenvalues and calculated in multiscale framework [45, 49–54]. The responses of the filters based on these functions are computed at different scales and are expected to have a maximum value at a scale corresponding to the width of the object [32, 54]. Such selectivity to the object shape along with capability to adaptively choose the optimal scale allows these filters to extract the looked-for objects at their native local scales. As noted by others [32, 55], this is especially important for tubular-like objects axis tracking applications.

The Hessian matrix (or simply Hessian) is the square matrix composed of second-order partial derivatives of some scalar-valued multivariable function , and it describes the local curvatures of this function . Assuming continuity of the second-order derivatives, the mixed derivatives do not depend on the order of differentiation (e.g., , etc.). The Hessian is then a symmetric matrix which for a 3D image can be written a matrix (see (5)) In (5), , , , and so on. Due to symmetric property of the Hessian, for 3D images only, six out of nine values have to be calculated. Let the eigenvalues of the Hessian be , , and with their corresponding eigenvectors , , and . If the eigenvalues ordered as , then the eigenvector gives the direction of the maximum of the second derivative.

Following the scale-space ideas described earlier, the partial second derivatives of the image in the Hessian have to be replaced by the γ-parameterized normalized Gaussian derivatives (see (4)), convolved with the image which results in that now the eigenvalues , , and become adjusted to the local size of the tubular object in an image.

Various algorithms for multiscale tubular object tracking and enhancement were developed depending on the way Hessian eigenvalues are combined in the objectness measure function. For instance, Sato et al. [49–51] suggested the “objectness” measure function which used only two out of three values of Hessian eigenvalues In (6), and (for a bright line on dark background), , and . That filter showed good performance in object enhancement in noisy environment. However, it did not have a parameter to control noise (background) suppression. Frangi et al. [52] extended the “objectness” measure function so that it includes the combination of all three Hessian eigenvalues and a factor with a parameter which controls noise suppression (see more details later in this paper). That function reflects shape and scale of the objects, has a single maximum on the center of the vessels segments, and has bell-shape close to Gaussian with width proportional to object size. This approach showed very good performance and has “de facto” become a basis for building even more sophisticated hybrid filters. Manniesing et al. [53] used this response function to develop an effective denoising filter, where the image intensity transformation is based on anisotropic “diffusion” governed by the “objectness” measure function in multiscale framework. Such an approach, which incorporates the shape pattern analysis along with multiscale data representation, would give us an extremely powerful tool to model artificial system learning. For neuron network reconstruction from 3D confocal microscope images, the tubularity measure function was used to design a statistical learning system for training a classifier and generating the probability that a given structure belongs to the tubular-like object [56, 57]. In pulmonology, the algorithms based on tube detectors were effectively used for airway and lung vascular tree reconstruction from 3D CT images [58, 59]. If properly normalized, the multiscale tube detectors could be used to build various cost and propagation functions required in the level set and fast marching segmentation algorithms [60, 61].

To take the full advantage of the power of the multiscale shape detector filter in object tracking algorithms applied for a large variety of medical applications, in this work we focus on the process of automation of this filter.

The filter itself has many control parameters which can be separated into several groups: Brightness measure (objects are bright relative to background); objectness measure (shape description), scale description (range plus scale step function), and background noise suppression parameter (Frobenius norm scale factor) [52]. Parameters in all groups, except the last one, describe general properties of the object itself so they do not depend on the imaging system characteristics. Since in vascular studies the object brightness, tubularity, and range of diameters are known beforehand, those parameters can be chosen in advance and then fixed. Hence, the only parameter which prevents the algorithm to be fully automated is control of the noise suppression. This parameter depends on the acquisition system and imaging conditions; therefore, it has to be experimentally found for each image set. Such a procedure is very compute expensive.

We present our development of the multiscale Hessian-based tubular object-tracking filter with automatic selection of the parameter used for suppression of background noise. That finalizes the automation of the filter. In our approach, the information required for the parameter calculation is acquired from the image being processed thus it automatically takes into account all the individual properties of the particular image such as voxel size and noise level. This allows for increased automation as well as parallel processing—thereby greatly decreasing processing time.

3. Methods

3.1. Images

For our studies, we used both gray-scale images numerically derived and acquired by scanners. The modeled images were programmed so to model environment with certain features. Tubular objects with various widths were placed amid different background: Gaussian random noise, nontubular objects, background with noise, and polynomial varying intensity. For simulations of images degraded by noise, we used the C++ classes contributed to the ITK Insight-Journal by Lehmann [62]. These noise simulation classes are implemented with multithread support and are based on the Mersenne Twister uniform pseudorandom number generator which has a period (2¹⁹⁹³⁷–1), 632-dimensional equidistribution, and up to 32-bit precision [63]. Thus, this generator could be considered as a “true random” which results in that generated noise does not produce any “neighborhood artifacts” or periodic patterns. The noise was generated with various Standard Deviations SD = 25.0, 50.0, 100.0, and Mean M = 0.0. We also used biomedical images of the heart acquired by the micro-CT scanner [64], Cerebellar Climbing Fibers [65–67], and Hippocampal CA3 Interneuron [67]. Specimen H61 (coronary artery branch within a human heart wall) was a methyl methacrylate cast prepared as described previously [68]. A cast of that coronary arterial tree was scanned with an isotropic voxel size of 0.018 mm and voxel CT image volume.

3.2. Micro-CT Scanner

The custom-made micro-CT scanners generate images up to isotropic voxels down to 4 μm on a side [64].

3.3. Server for Image Processing

To be able to process large images using the developed algorithms, we built a specialized server with four 64 bit AMD Opteron 8350 Quad Core 2.0 GHz CPUs and 128 GB memory. The server is located in a server room and it is accessible in multiuser mode through the local network using remote clients.

3.4. Algorithms

For our software development, we used the library of C++ classes from the National Library of Medicine Insight Segmentation and Registration Toolkit (ITK) [69–71]. The library was compiled with multithread support based on the POSIX thread (Pthreads) model [72] using 64 bit C++ compiler GCC 4.3.2 [73–75] and installed on 64 bit Debian Linux 2.6.26 [76, 77].

3.4.1. Automation of Multiscale Shape Detector Response Function

The developed multiscale shape detector filter is based on the objectness measure function suggested by Frangi et al. [52] and the C++ classes contributed to the ITK Insight-Journal [78–80]. After thoroughly conducted studies and tests, we found that the C++ classes contributed by Antiga [80] satisfy our purposes the best; therefore, our further developments are based on those classes.

Let be the eigenvalue of the Hessian matrix at voxel x ordered such that (we drop the dependency to x). In the case of the ideal bright tubular structure, the voxels should satisfy the following relation for eigenvalues ; ; and for bright objects both and must be negative.

Frangi et al. [52] proposed to use the eigenvalues to define Vesselness measure as below: where , , and are the parameters that control the sensitivity of the filter to the measures , , and . The measures have the following meaning. is used to distinguish between plate-like and line-like patterns. is used to derive a blob-like pattern. The measures and are gray intensity level invariant and capture only the topological information of the objects in the image. The choice of the control parameters and defines the object pattern to be studied. For example, if tubular-like objects are chosen, the parameters are fixed to , [54].

The parameter (see (7)), controls background noise suppression in the Hessian-based object enhancement filter. The Frobenius Hessian matrix norm is chosen as a measure to distinguish background noisy pixels. Since the parameter strongly depends on the individual properties of particular image such as voxel size and noise level, for every new study, the optimal parameter value should be experimentally found again by trial. The very wide search range of the optimal parameter value in concert with a highly time-expensive calculation to derive the Hessian taken in multiscale framework makes this algorithm very labor-inefficient especially for large 3D biomedical images with high resolution. For example, one trial run to process the 16 bit gray-scale micro-CT image on our server took about an hour; thus, the interactive search for the parameter might take hours for the user who operates the program interactively.

The method for automating the selection of the parameter for suppression of background noise uses a scalar function (nondirectional) of the image voxels the Laplacian of the image. The Laplacian is a well known operator in image processing which is easy to calculate [81–83]. Ultimately, the Laplacian calculates the trace of Hessian matrix or, equivalently, the sum of its eigenvalues () making it invariant with respect to a change of tensor basis. This characterization can be used to steer the control parameter responsible for noise suppression. The schematic description of complete algorithm is below.(1)For each voxel in the image, calculate the Laplacian.(2)In the calculated Laplacian array, find the maximum value of Laplacian, that is, .(3)Take one tenth of that maximum value of Laplacian, that is, .(4)Assign the calculated value to parameter .

In this approach, the information required for parameter calculation is acquired solely from the image being processed; thus, it automatically takes into account all the individual properties of the particular image such as a voxel size and noise level.

3.4.2. Objective Measures for Perceptual Quality Evaluation of Images

There are two measures commonly used for objective evaluation of the perceptual quality of images: Mean Square Error (MSE) and Peak Signal to Noise Ratio (PSNR) [84–90]. These quality measures are defined as follows. Consider two images being compared and , where is the number of points (pixels) in the data sets and and are intensity levels in the images. Let x be an “ideal image” and be a “degraded image”. The MSE measure is then defined as and PSNR is, respectively, where is the dynamic range of allowable pixel intensities. As follows from the definitions above, the lower the values of MSE, the lower the error and the higher the PSNR, the better quality of processed image.

The MSE and PSNR measures were developed using C++ classes from the ITK library [69] and contribution [62].

4. Results

4.1. Manual Mode

In manual mode, the control parameters have to be provided by the operator before running the code. If the result is not acceptable, the operator has to change parameters and rerun the code again. The MSE and PSNR measures were used for objective quality evaluation of the processed image. As an “ideal image”, there was used a modeled image comprised three tube-like objects with Gaussian intensity profiles of different width. Then, the ideal image was degraded by the random Gaussian noise with Standard Deviation SD = 100.0 and Mean M = 0.0 and processed in manual mode by the algorithm. The background suppression parameter value was sampled over a wide range 10–500 to surely cover prospective optimal control parameter. The results are depicted in Figure 7. It can be seen that both measures indicate the optimal value is about which corresponds to MSE value 0.00214 and PSNR value 26.69. If run in automatic mode, the algorithm finds even more precise the optimal value which corresponds to MSE value 0.00203 and PSNR value 26.93. These results demonstrate capability of the developed algorithm in automatically finding the optimal control parameter for background suppression.

4.2. Fully Automatic Mode

To test efficiency of the algorithm in fully automatic mode, we used the modeled images described above. The control parameters were chosen and fixed: brightness “on”, objectness measure “tubular”, scale range “1–30”, scale steps “20”, step function “logarithmic”, noise suppression mode “automatic”. Scale range was chosen wide enough to cover all possible diameters. The step function was made “logarithmic” so as to emphasize finer scales.

4.2.1. Background with Gaussian Random Noise and Nontubular Objects

First, we processed the images with curved tubular objects with various widths which were placed amid nontubular objects. The images were degraded by Gaussian random noise with SD = 25.0, 50.0, and M = 0.0. As can be seen from Figure 8, the algorithm automatically finds the optimal parameters and successfully tracks tubular objects.

(a)

(b)

(c)

(d)

4.2.2. Objects in Images with Nonlinear Background

We also tested performance of the algorithm for neurological confocal microscopy image processing affected by tiling, shading, Gaussian noise, nonlinear background, and so forth (see, e.g., Figure 9). The modeled images were programmed so to present three parallel tubular objects aligned in the Z-direction with different diameters and distance between objects along with gray-scale distribution described by Gaussian intensity profiles. The polynomial background was added with and without Gaussian random noise (Figure 10). In Figure 11, there are depicted intensity profiles across the tubular objects with Gaussian intensity profiles in the processed images with various polynomial noisy backgrounds added SD = 0.0, 25.0, 50.0, and Mean M = 0.0. As could be seen, the algorithm effectively suppresses the background and successfully extracts the tubular images.

(a)

(b)

(a)

(b)

(c)

(d)

(a)

(b)

(c)

4.3. Speed Efficiency of Fully Automatic Mode

To evaluate time efficiency of our automatic method, we processed the (16 bits per voxel) micro-CT image with 0.018 cubic mm size in multithreaded mode using a Linux-based server with the 16 processors and 128 GB memory as described above. First, we ran the automatic algorithm to find the parameter . Then, we ran nonautomatic (the found parameter was manually plugged into a program) version. The amount of CPU time spent for a single run in non-automatic mode is 58 minutes whereas for automatic mode is 64 minutes, which is just about 10% longer. If the server is able to allocate enough memory for the run, the time spent by the algorithm in fully automated mode is comparable with time for the nonautomatic mode.

4.4. Biomedical Images

The result of applying our automatic algorithm to the micro-CT image of the coronary arteries in a heart in concert with the nonprocessed image is shown in Figure 12 using Maximum Intensity Projection (MIP) images. Before preprocessing (a) there are the “blobs” of white material due to a contrast medium accumulation in a cardiac chamber. Those “blobs” are selectively removed by the algorithm as can be seen in right panel (postprocessing).

(a)

(b)

In Figure 13, there are the not processed (a) and processed (b) MIP images of cerebellar climbing fibers. As expected, the tubular measure filter effectively suppresses the fixed pattern background and noise and “delineates” the tubular objects.

(a)

(b)

We also explored the efficiency of the algorithm as an initial filter in the central line extraction pipe line. The original image Sample H61 was processed with the developed filter and then segmented using the region growing connected threshold algorithm [69–71]. The seed for segmentation was allocated at the beginning of the main tree root thus the side trees were excluded. Afterwards, a center line was extracted using the 3D thinning approach based on the C++ classes submitted to the insight journal by Homann [91]. Images along with the extracted center lines are depicted in Figure 14. As could be seen from the figure, the algorithm effectively suppresses background and delineates the tree. The centerlines are extracted correctly as well.

(a)

(b)

5. Discussion and Future Work

We have presented a method for automation of adaptive nonsupervised system for tracking tubular objects that is based on analysis of local structures performed in multiscale framework. The designed filter has demonstrated a great potential for complete automation and showed very good performance in both background noise suppression and tubular object tracking.

The developed approach can be used in the reconstruction pipeline right after image deconvolution operation. Even though the convolution operator will reconstruct the object features at finer scales, those features will appear in increased noise environment which in return might require additional postprocessing for noise suppression yet to preserve extracted features.

Another application is the object feature extraction pipeline. This filter can be used as a preprocessing filter for vessel enhancement and background noise suppression right before segmentation or immediately in the segmentation algorithms itself, for instance, in the family of segmentation algorithms which require distributed seeds [70, 71]. By using the thresholded output of the vessel enhancement filter as a seeder, it can increase the speed efficiency of the segmentation process considerably.

Since the response function is built using exponents, with proper normalization this function can be considered as a probability function with values distributed over the interval “0.0-1.0”. In this case, after processing, the output image holds voxels with values of probability of the event that “a voxel belongs to the object with tubular shape”. These probabilities can be used in many ways. The most traditional way is to rescale it back to a gray-scale image. Although such images do not keep a proper intensity calibration, they still can be used for morphometric analysis. If calibration is of concern, the probabilities could be converted to a mask for sampling the original micro-CT image from which the calibration could be recovered.

Since the filter generates the response function with only one maximum across scale space at a scale that is proportional to the diameter of the tubular object and that maximum is located at the center of the object, the probability image is more suitable to construct various cost functions. The images with cost functions can further be used as the “feature image” in various image processing pipelines, for instance, such as in flux-driven centerline extraction algorithms [92, 93], level-set and fast-marching segmentation algorithms [61, 69–71], and so forth. As the multiscale vessel enhancement filter is very robust against noise, it can be superior over traditional approaches like, for example, in a filter pipeline “segmented image, distance map, and cost function”, since it can directly generate the cost function avoiding steps for producing segmentation and distance map. In addition, the multiscale space feature can be used to build the cost function in multiscale representation and use it for multiscale vessel tracking as suggested in [54].

Acknowledgments

This work was supported in part by NIH Grant no. EB 000305. The authors wish to thank Ms. Diane R. Eaker for her help with the micro-CT reconstructions and Point Spread Function deconvolution and Ms. Delories C. Darling for formatting this manuscript. The authors are exceedingly grateful to Professor Giorgio Ascoli and Professor German Barrionuevo for the permission to use their images in our work.

References

D. Marr, Vision: A Computational Investigation into the Human Representation and Processing of Visual Information, The MIT Press, Cambridge, MA, USA, 2010.
B. M. ter Haar Romeny, Front-End Vision and Multi-Scale Image Analysis: Multi-scale Computer Vision Theory and Applications, Written in Mathematica, Springer, Berlin, Germany, 1st edition, 2003.
A. P. Witkin, “Scale-space filtering,” in Proceedings of the 8th International Joint Conference on Artificial Intelligence, pp. 1019–1022, Karlsruhe, Germany, 1983.
View at: Google Scholar
J. J. Koenderink, “The structure of images,” Biological Cybernetics, vol. 50, no. 5, pp. 363–370, 1984.
View at: Google Scholar
R. A. Hummel, “Representations based on zero-crossings in scale-space,” in Readings in Computer Vision: Issues, Problems, Principles, and Paradigms, M. A. Fischler and O. Firschein, Eds., pp. 753–758, Morgan Kaufmann, Los Altos, Calif, USA, 1987.
View at: Google Scholar
A. L. Yuille and T. A. Poggio, “Scaling theorems for zero-crossings,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, no. 1, pp. 15–25, 1986.
View at: Google Scholar
J. J. Koenderink and A. J. van Doorn, “Generic neighborhood operators,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 6, pp. 597–605, 1992.
View at: Google Scholar
M. J. Florack, The syntactical structure of scalar images, Ph.D. thesis, Department of Medical Physics, Utrecht University, Utrecht, Netherlands, 1993.
T. Lindeberg, “Scale-space for discrete signals,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 3, pp. 234–254, 1990.
View at: Publisher Site | Google Scholar
T. Lindeberg, Scale-Space Theory in Computer Vision, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1994.
J. J. Koenderink, “Scale-time,” Biological Cybernetics, vol. 58, no. 3, pp. 159–162, 1988.
View at: Publisher Site | Google Scholar
T. Lindeberg and D. Fagerström, “Scale-space with casual time direction,” in Proceedings of the 4th European Conference on Computer Vision (ECCV '96), B. Buxton and R. Cipolla, Eds., vol. 1, pp. 229–240, Cambridge, UK, 1996.
View at: Google Scholar
T. Lindeberg, “Linear spatio-temporal scale-space,” in Proceedings of the First International Conference on Scale-Space Theory in Computer Vision (Scale-Space '97), B. M. ter Haar Romeny, L. Florack, J. Koenderink, and M. Viergever, Eds., pp. 113–127, Springer, Berlin, Germany, 1st edition, 1997.
View at: Google Scholar
P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990.
View at: Publisher Site | Google Scholar
J. Weickert, “A review of nonlinear diffusion filtering,” in Proceedings of the First International Conference on Scale-Space Theory in Computer Vision (Scale-Space '97), B. Buxton and R. Cipolla, Eds., pp. 3–28, Springer, 1997.
View at: Google Scholar
L. M. J. Florack, A. H. Salden, B. M. ter Haar Romeny, J. J. Koenderink, and M. A. Viergever, “Nonlinear scale-space,” Image and Vision Computing, vol. 13, no. 4, pp. 279–294, 1995.
View at: Google Scholar
B. Burgeth, S. Didas, and J. Weickert, “The Bessel Scale-Space,” in Proceedings of the 1st International Workshop Deep Structure, Singularities, and Computer Vision (DSSCV '05), O. F. Olsen, L. Florack, and A. Kuijper, Eds., vol. 3753 of Lecture Notes in Computer Science, pp. 84–95, Springer, Maastricht, The Netherlands, June 2005.
View at: Google Scholar
C. Heil and D. F. Walnut, Fundamental Papers in Wavelet Theory, Princeton University Press, 2006.
S. Mallat and W. L. Hwang, “Singularity detection and processing with wavelets,” IEEE Transactions on Information Theory, vol. 38, no. 2, pp. 617–643, 1992.
View at: Publisher Site | Google Scholar
D. Doering and A. Schuck, “A novel method for generating scale space kernels based on wavelet theory,” Revista de Informática Teórica e Aplicada, vol. 15, no. 2, pp. 121–137, 2008.
View at: Google Scholar
T. T. Vu and M. Tokunaga, “Wavelet and scale-space theory in segmentation of airborne laser scanner data,” in Proceedings of the 22nd Asian Conference on Remote Sensing (ACRS '01), vol. 1, pp. 176–180, Singapore, November 2001.
View at: Google Scholar
Y. Ferdman, C. Sagiv, and N. Sochen, “Full affine wavelets are scale-space with twist,” in Proceedings of the 1st International Conference on Scale-Space and Variational Methods in Computer Vision (SSVM '07), F. Sgallari, A. Murli, and N. Paragios, Eds., pp. 1–12, Springer, Ischia, Italy, June 2007.
View at: Google Scholar
YU. P. Wang and S. L. Lee, “Scale-space derived from B-splines,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 20, no. 10, pp. 1040–1055, 1998.
View at: Google Scholar
J. Babaud, A. P. Witkin, M. Baudin, and R. O. Duda, “Uniqueness of the Gaussian kernel for scale-space filtering,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, no. 1, pp. 26–33, 1986.
View at: Google Scholar
L. M. J. Florack, B. M. ter Haar Romeny, J. J. Koenderink, and M. A. Viergever, “Scale and the differential structure of images,” Image and Vision Computing, vol. 10, no. 6, pp. 376–388, 1992.
View at: Google Scholar
L. M. J. Florack, B. M. ter Haar Romeny, J. J. Koenderink, and M. A. Viergever, “Linear scale-space,” Journal of Mathematical Imaging and Vision, vol. 4, no. 4, pp. 325–351, 1994.
View at: Publisher Site | Google Scholar
R. Duits, L. Florack, J. De Graaf, and B. Ter Haar Romeny, “On the axioms of scale space theory,” Journal of Mathematical Imaging and Vision, vol. 20, no. 3, pp. 267–298, 2004.
View at: Publisher Site | Google Scholar
T. Wennekers and K. Suder, “Fitting of spatio-temporal receptive fields by sums of Gaussian components,” Neurocomputing, vol. 58–60, pp. 929–934, 2004.
View at: Publisher Site | Google Scholar
R. A. Young, “The Gaussian derivative model for spatial vision: I. Retinal mechanisms,” Spatial Vision, vol. 2, no. 4, pp. 273–293, 1987.
View at: Google Scholar
R. A. Young, R. M. Lesperance, and W. W. Meyer, “The Gaussian derivative model for spatial-temporal vision: I. Cortical model,” Spatial Vision, vol. 14, no. 3-4, pp. 261–319, 2001.
View at: Publisher Site | Google Scholar
R. A. Young and R. M. Lesperance, “The Gaussian derivative model for spatial-temporal vision: II. Cortical data,” Spatial Vision, vol. 14, no. 3-4, pp. 321–389, 2001.
View at: Publisher Site | Google Scholar
T. Lindeberg, “Feature detection with automatic scale selection,” International Journal of Computer Vision, vol. 30, no. 2, pp. 79–116, 1998.
View at: Google Scholar
T. Lindeberg, “Scale-space theory: a basic tool for analyzing structures at different scales,” Journal of Applied Statistics, vol. 21, no. 2, pp. 225–270, 1994.
View at: Google Scholar
T. Lindeberg, “Edge detection and ridge detection with automatic scale selection,” International Journal of Computer Vision, vol. 30, no. 2, pp. 117–154, 1998.
View at: Google Scholar
D. Marr, “Early processing of visual information,” Philosophical transactions of the Royal Society of London B, vol. 275, no. 942, pp. 483–519, 1976.
View at: Google Scholar
D. Marr and E. Hildreth, “Theory of edge detection,” Proceedings of the Royal Society of London, vol. 207, no. 1167, pp. 187–217, 1980.
View at: Google Scholar
J. Canny, “A computational approach to edge detection,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, no. 6, pp. 679–698, 1986.
View at: Google Scholar
P. Majer, “The influence of the gamma-parameter on feature detection with automatic scale Selection,” in Proceedings of the 3rd International Conference on Scale-Space and Morphology in Computer Vision, M. Kerckhove, Ed., pp. 245–254, Springer, 2001.
View at: Google Scholar
P. Majer, “On the influence of scale selection on feature detection for the case of linelike structures,” International Journal of Computer Vision, vol. 60, no. 3, pp. 191–202, 2004.
View at: Publisher Site | Google Scholar
C. Lorenz, I.-C. Carlsen, T.-M. Buzug, C. Fassnacht, and J. Weese, “A multi-scale line filter with automatic scale selection based on the Hessian matrix for medical image segmentation,” in Proceedings of the First International Conference on Scale-Space Theory in Computer Vision (Scale-Space '97), B. M. ter Haar Romeny, L. Florack, J. Koenderink, and M. Viergever, Eds., pp. 152–163, Springer, 1997.
View at: Google Scholar
E. L. Ritman, “Micro-computed tomography—current status and developments,” Annual Review of Biomedical Engineering, vol. 6, pp. 185–208, 2004.
View at: Publisher Site | Google Scholar
A. Faridani and E. L. Ritman, “High-resolution computed tomography from efficient sampling,” Inverse Problems, vol. 16, no. 3, pp. 635–650, 2000.
View at: Google Scholar
S. T. Witt, C. H. Riedel, M. Goessl, M. S. Chmelik, and E. L. Ritman, “Point spread function deconvolution in 3D micro-CT angiography for multiscale vascular tree separation,” in Physics of Medical Imaging, vol. 5030 of Proceedings of SPIE, pp. 720–727, 2003.
View at: Google Scholar
D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Magazine, vol. 13, no. 3, pp. 43–64, 1996.
View at: Google Scholar
T. M. Koller, G. Gerig, G. Szekely, and D. Dettwiler, “Multiscale detection of curvilinear structures in 2-D and 3-D image data,” in Proceedings of 15th International Conference on Computer Vision (ICCV '95), E. Grimson, S. Shafer, A. Blake, and K. Sugihara, Eds., pp. 864–869, Springer, 1995.
View at: Google Scholar
B. M. ter Haar Romeny, L. Florack, J. Koenderink, and M. Viergever, Eds., Scale-Space Theory in Computer Vision, Springer, Berlin, Germany, 1st edition, 1997.
R. Kimmel, N. Suchen, and J. Weickert, Eds., Scale Space and PDE Methods in Computer Vision, Springer, Berlin, Germany, 1st edition, 2005.
S. R. Aylward and E. Bullitt, “Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction,” IEEE Transactions on Medical Imaging, vol. 21, no. 2, pp. 61–75, 2002.
View at: Publisher Site | Google Scholar
Y. Sato, S. Nakajima, H. Atsumi et al., “3D multi-scale line filter for segmentation and visualization of curvilinear structures in medical images,” in Proceedings of the 1st Joint Conference Computer Vision, Virtual Reality and Robotics in Medicine and Medical Robotics and Computer-Assisted Surgery Grenoble (CVRMed-MRCAS '97), J. Troccaz, E. Grimson, and R. Mosges, Eds., pp. 213–222, Springer, 1997.
View at: Google Scholar
Y. Sato, S. Nakajima, N. Shiraga et al., “Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images,” Medical Image Analysis, vol. 2, no. 2, pp. 143–168, 1998.
View at: Google Scholar
Y. Sato, C. F. Westin, A. Bhalerao et al., “Tissue classification based on 3D local intensity structures for volume rendering,” IEEE Transactions on Visualization and Computer Graphics, vol. 6, no. 2, pp. 160–180, 2000.
View at: Google Scholar
A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” in Proceedings of the 1st International Conference Medical Image Computing and Computer-Assisted Intervention (MICCAI '98), W. M. Wells, A. Colchester, and S. L. Delp, Eds., pp. 130–137, Springer, Cambridge, Mass, USA, October 1998.
View at: Google Scholar
R. Manniesing, M. A. Viergever, and W. J. Niessen, “Vessel enhancing diffusion. A scale space representation of vessel structures,” Medical Image Analysis, vol. 10, no. 6, pp. 815–825, 2006.
View at: Publisher Site | Google Scholar
O. Wink, W. J. Niessen, and M. A. Viergever, “Multiscale Vessel Tracking,” IEEE Transactions on Medical Imaging, vol. 23, no. 1, pp. 130–133, 2004.
View at: Publisher Site | Google Scholar
S. D. Olabarriaga, M. Breeuwer, and W. J. Niessen, “Evaluation of Hessian-based filters to enhance the axis of coronary arteries in CT images,” International Congress Series, vol. 1256, pp. 1191–1196, 2003.
View at: Google Scholar
A. Santamaria-Pang, T. S. Bildea, C. Colbert, P. Saggau, and I. A. Kakadiaris, “Towards segmentation of irregular tubular structures in 3D confocal microscope images,” in Proceedings of the Workshop in Microscopic Image Analysis and Applications in Biology (MIAAB '06), pp. 78–85, Copenhagen, Denmark, October 2006.
View at: Google Scholar
G. Gonzalez, F. Aguet, F. Fleuret, M. Unser, and P. Fua, “Steerable features for statistical 3D dendrite detection,” in Proceedings of 12th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI '09), vol. 5762, pp. 625–632, Springer, London, UK, September 2009.
View at: Google Scholar
C. Bauer, T. Pock, H. Bischof, and R. Beichel, “Airway tree reconstruction based on tube detection,” in Proceedings of the 2nd International Workshop on Pulmonary Image Analysis, M. Brown, M. de Bruijne, B. van Ginneken, and A. Kiraly, Eds., pp. 203–214, CreateSpace, London, UK, September 2009.
View at: Google Scholar
M. Sonka, H. Shikata, G. McLennan, and E. A. Hoffman, “Segmentation of pulmonary vascular trees from thoracic 3D CT images,” International Journal of Biomedical Imaging, vol. 2009, Article ID 636240, 11 pages, 2009.
View at: Publisher Site | Google Scholar
S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988.
View at: Google Scholar
J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, New York, NY, USA, 2nd edition, 1999.
G. Lehmann, “Noise simulation,” 2010, http://hdl.handle.net/10380/3158.
View at: Google Scholar
M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Transactions on Modeling and Computer Simulation, vol. 8, no. 1, pp. 3–30, 1998.
View at: Publisher Site | Google Scholar
S. M. Jorgensen, O. Demirkaya, and E. L. Ritman, “Three-dimensional imaging of vasculature and parenchyma in intact rodent organs with X-ray micro-CT,” American Journal of Physiology, vol. 44, no. 3, pp. H1103–H1114, 1998.
View at: Google Scholar
I. Sugihara, H. S. Wu, and Y. Shinoda, “Morphology of single olivocerebellar axons labeled with biotinylated dextran amine in the rat,” Journal of Comparative Neurology, vol. 414, no. 2, pp. 131–148, 1999.
View at: Google Scholar
K. M. Brown, D. E. Donohue, G. D'Alessandro, and G. A. Ascoli, “A cross-platform freeware tool for digital reconstruction of neuronal arborizations from image stacks,” Neuroinformatics, vol. 3, no. 4, pp. 343–359, 2005.
View at: Publisher Site | Google Scholar
“The digital reconstruction of axonal and dendritic morphology (DIADEM 2010) competition,” 2010, http://diademchallenge.org/.
View at: Google Scholar
M. Zamir, “Tree structure and branching characteristics of the right coronary artery in a right-dominant human heart,” Canadian Journal of Cardiology, vol. 12, no. 6, pp. 593–599, 1996.
View at: Google Scholar
“Insight Segmentation and Registration Toolkit (ITK),” http://www.itk.org/.
View at: Google Scholar
L. Ibanez and W. Schroeder, The ITK Software Guide 2.4, Kitware, Clifton Park, NY, USA, 2005.
T. S. Yoo, Insight into Images: Principles and Practice for Segmentation, Registration, and Image Analysis, AK Peters, Wellesley, Mass, USA, 1st edition, 2004.
B. Nichols, D. Buttlar, and J. B. Farrell, Pthreads Programming: A POSIX Standard for Better Multiprocessing, O’Reilly Media, Sebastopol, Calif, USA, 1st edition, 1996.
B. J. Gough and R. M. Stallman, An Introduction to GCC, Network Theory, Bristol, UK, 1st edition, 2004.
W. von Hagen, The Definitive Guide to GCC, Apress, Berkely, Calif, USA, 2nd edition, 2006.
R. M. Stallman and GCC Developer Community, Using the Gnu Compiler Collection: A Gnu Manual for GCC Version 4.3.3, GNU Press, Boston, Mass, USA, 2009.
“Debian Linux Operating System,” www.debian.org.
View at: Google Scholar
M. Krafft, The Debian System: Concepts and Techniques, Pap/Dvdr ed. No Starch Press, San Francisco, Calif, USA, 2005.
A. Enquobahrie, L. Ibanez, A. Aylward, and E. Bullitt, “Vessel Enhancing Diffusion Filter,” http://hdl.handle.net/1926/558.
View at: Google Scholar
R. Manniesing, M. Viergever, and W. Niessen, “Vessel Enhancing Diffusion,” 2009, http://hdl.handle.net/1926/1515.
View at: Google Scholar
L. Antiga, “Generalizing vesselness with respect to dimensionality and shape,” 2007, http://hdl.handle.net/1926/576.
View at: Google Scholar
R. C. Gonzalez and R. E. Woods, Digital Image Processing, Prentice Hall, Upper Saddle River, NJ, USA, 3rd edition, 2007.
M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision, Brooks/Cole, Pacific Grove, Calif, USA, 2nd edition, 2007.
J. R. Parker, Algorithms for Image Processing and Computer Vision, John Wiley & Sons, Hoboken, NJ, USA, 1996.
Z. Wang and A. C. Bovik, Modern Image Quality Assessment, Synthesis Lectures on Image, Video & Multimedia Processing, Morgan & Claypool, San Rafael, Calif, USA, 1st edition, 2006.
I. Avcibaş, B. Sankur, and K. Sayood, “Statistical evaluation of image quality measures,” Journal of Electronic Imaging, vol. 11, no. 2, pp. 206–223, 2002.
View at: Publisher Site | Google Scholar
Q. Huynh-Thu and M. Ghanbari, “Scope of validity of PSNR in image/video quality assessment,” Electronics Letters, vol. 44, no. 13, pp. 800–801, 2008.
View at: Publisher Site | Google Scholar
Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
View at: Publisher Site | Google Scholar
Z. Wang, A. C. Bovik, and L. Lu, “Why is image quality assessment so difficult?” in Proceedings of the IEEE International Conference on Acoustic, Speech, and Signal Processing (ICASSP '02), vol. 4, pp. IV/3313–IV/3316, Orlando, Fla, USA, May 2002.
View at: Google Scholar
A. M. Eskicioglu and P. S. Fisher, “Image quality measures and their performance,” IEEE Transactions on Communications, vol. 43, no. 12, pp. 2959–2965, 1995.
View at: Publisher Site | Google Scholar
Z. Wang and A. C. Bovik, “Mean squared error: lot it or leave it? A new look at signal fidelity measures,” IEEE Signal Processing Magazine, vol. 26, no. 1, pp. 98–117, 2009.
View at: Publisher Site | Google Scholar
H. Homann, “Implementation of a 3D thinning algorithm,” 2007, http://hdl.handle.net/1926/1292.
View at: Google Scholar
S. Bouix, K. Siddiqi, and A. Tannenbaum, “Flux driven automatic centerline extraction,” Medical Image Analysis, vol. 9, no. 3, pp. 209–221, 2005.
View at: Publisher Site | Google Scholar
X. Mellado, I. Larrabide, M. Hernandez, and A. Frangi, “Flux driven medial curve extraction,” 2007, http://hdl.handle.net/1926/560.
View at: Google Scholar

Copyright

Copyright © 2011 Oleksandr P. Dzyubak and Erik L. Ritman. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1956

Downloads

1312

Citations