Automated Detection and Measurement of Isolated Retinal Arterioles by a Combination of Edge Enhancement and Cost Analysis

Pressure myography studies have played a crucial role in our understanding of vascular physiology and pathophysiology. Such studies depend upon the reliable measurement of changes in the diameter of isolated vessel segments over time. Although several software packages are available to carry out such measurements on small arteries and veins, no such software exists to study smaller vessels (<50 µm in diameter). We provide here a new, freely available open-source algorithm, MyoTracker, to measure and track changes in the diameter of small isolated retinal arterioles. The program has been developed as an ImageJ plug-in and uses a combination of cost analysis and edge enhancement to detect the vessel walls. In tests performed on a dataset of 102 images, automatic measurements were found to be comparable to those of manual ones. The program was also able to track both fast and slow constrictions and dilations during intraluminal pressure changes and following application of several drugs. Variability in automated measurements during analysis of videos and processing times were also investigated and are reported. MyoTracker is a new software to assist during pressure myography experiments on small isolated retinal arterioles. It provides fast and accurate measurements with low levels of noise and works with both individual images and videos. Although the program was developed to work with small arterioles, it is also capable of tracking the walls of other types of microvessels, including venules and capillaries. It also works well with larger arteries, and therefore may provide an alternative to other packages developed for larger vessels when its features are considered advantageous.


Introduction
Pressure myography is widely used to monitor the responses of small arteries and veins to physical and chemical stimuli. This technique involves mounting a small isolated vessel segment between two glass cannulae and pressurising the vessel to an appropriate transmural pressure. By transferring the myograph to the stage of a microscope fitted with a video camera, changes in vasomotor activity can then be continuously imaged throughout the experiment. Pressure myography has been instrumental in our current understanding of the myogenic response (i.e. the intrinsic property of blood vessels to respond dynamically to changes in intraluminal pressure) and in characterising the pharmacological actions of drugs and other vasoactive compounds on the vascular system [1], [2]. Importantly, it has also provided crucial insights into the pathogenesis of vascular dysfunction in a range of different disease states, including, for example, hypertension, diabetes, obesity and stroke [3], [4], [5], [6], [7].
A key technical aspect of any myography-based study is the ability to accurately measure changes in vessel diameter over time. Traditional methods of manually measuring inner and outer vessel diameters remain popular [8], [9], [10], [11], but are timeconsuming and prone to user bias [12]. Several automatic or semiautomatic algorithms have been developed to assist with these measurements. The SoftEdge Myocyte Cell Length Acquisition Module [13], [14], [15], [16] detects and tracks vessel walls within a couple of user-defined windows after a threshold value and a 'crossing condition' have been set. Another algorithm, developed by Kim et al. [17], involves the user defining a line perpendicular to the vessel and extending beyond the outer walls for intensity analysis at a high-contrast region of the vessel before carrying out Otsu's thresholding. A more sophisticated algorithm, VesselTrack [18], [19], [20], detects the abluminal edges of the vessel walls in two small user-defined windows and tracks vessel movements in these regions using iterative regression. The adluminal edges are subsequently estimated from the detected outer values. Other packages that have also been used to study isolated vessels are MyoView (DMT) [21], [22] and the program Mary (Nihil, Lund, Sweden) [23]. In general, all of these programs rely on the use of thresholding methods which are well-suited to the analysis of arteries and veins where there is normally a high level of contrast between the vessel wall and other areas of the image.
Arterioles are small blood vessels that act as the major vascular resistance elements controlling blood flow from arteries to capillaries, and the properties of these vessels are of particular interest in improving our understanding of how local tissue perfusion and capillary pressure are regulated [24]. Although small vessels (,50 mm in diameter) are unsuitable for use with commercially available pressure myography systems, we have recently developed methods to carry out myography studies on isolated arterioles from the rat retina [25]. Like small arterioles from many other vascular beds, these vessels are surrounded by just a single layer of smooth muscle cells and consequently the vessel walls appear non-uniform, with no unambiguous change in contrast relative to the external background or intraluminal space [26]. This renders detection algorithms based solely on thresholding methods inadequate. We present here a free and open source ImageJ plug-in, MyoTracker, for the automatic tracking and measurement of the diameter of small isolated retinal arterioles in image files. The method used is based on a combination of edge enhancement and cost analysis. Thresholding is only used in the algorithm as a check on the detected borders. Our software can run on individual images and on stacks of thousands of images and is also amenable to the analysis of other vessel types, including larger arterial and arteriolar vessels, venules and capillaries.

Ethics Statement
All animal tissue was obtained by schedule 1 methods in accordance with the Animals (Scientific Procedures) Act 1986 and with the agreement of the Queen's University of Belfast Animal Welfare and Ethical Review Body for which a specific project licence is not required.

Vessel Preparation
Male Sprague-Dawley rats (8-12 weeks of age; 200-250 g) were euthanized using CO 2 in accordance with UK legal and local institutional requirements. Retinas were placed in low Ca 2+ Hanks' solution containing (mmol/L): 140, NaCl; 6, KCl; 5, Dglucose; 0.1, CaCl 2 ; 1.3, MgCl 2 ; 10, HEPES (pH 7.4 with NaOH) and mechanically triturated using a Pasteur pipette. The tissue was pipetted into a recording bath mounted on an inverted microscope and isolated retinal arterioles, venules and capillaries (,5-40 mm diameter) were identified as previously described [25], [27]. In one set of experiments, larger bovine arteries were used. Bovine eyes were obtained from a local abattoir and transported back to the laboratory in low Ca 2+ Hanks' solution at 1uC. The retinas were removed and the arteries mechanically isolated using the methods described above. All experiments were carried out at 37uC.
Following introduction of the pipette, the vessel was superfused with Hanks' solution containing 2 mmol/L Ca 2+ (added to the above solution as CaCl 2 ) for 10-15 min, allowing the pipette to seal to the inner vessel wall. Intraluminal pressure was regulated by changing the height of a fluid reservoir connected to the inflow cannula and monitored using a pressure transducer. Individual vessels were viewed under a 20x, NA 0.4 objective and images (saved as BMP images of 128061024 pixels; 8-bit; 1.2 MB) captured at a rate of 140 images per minute using a MCN-B013-U USB camera. Acquisition was carried out using custom software implemented in Delphi. To convert pixels to microns, 102 pixels were found to be equivalent to 27 microns.
In separate experiments, vasoconstrictor responses of isolated, non-pressurised, retinal arterioles to 10 mM caffeine or 10 nM endothelin-1 (Et-1) were imaged using the same videomicroscopy set-up. Drugs were delivered via a gravity-fed multi-channel perfusion manifold connected to a single outlet needle (350 mm in diameter) that was positioned adjacent to the vessel of interest. Figure 1 illustrates the different steps used by MyoTracker to measure the diameter of an isolated retinal arteriole ( Fig. 1 A). The software automatically detects the walls of the vessel, draws two lines along their centres, and measures the overall diameter as the mean of the vertical distances between these lines, at approximately one pixel intervals. A cost function was developed to assist with the drawing of the lines. This function balances out two main constraining factors: (a) the distance between the pixel being drawn and initially estimated points at the start and end of the vessel, and (b) the intensity changes (from darker to lighter values) along the vertical line where the pixel is located. More specifically, the cost function can be written as:

Explanation of the Algorithm
where ds is defined as the distance of a pixel to a starting point and de as the distance to an end point (i.e. distances to the lowest intensity points at the start and end of the vessel over five columns), and i is the intensity of the pixel normalised by the average intensity over a number of pixels above and below it vertically. The function F is defined as follows: The advantages of using this cost function are (1) its efficiency of processing, (2) it achieves global optimisation using all the pixels on the vessel, and (3) it is sufficiently flexible, in that it only depends on the resulting terminal points and the intensity of the pixels. Given the high variability in intensity values observed along the walls of these small vessels, however, using the cost function alone usually gave 'noisy' lines (i.e. with many jumps between pixels at different intensities close to each other). This can be seen in the top right panels in Figure 1 (Fig. 1 A i , blue lines). To improve on this, several processing steps are carried out: Step 1. Two filters are used to enhance the adluminal and abluminal edges of the vessel walls ( Fig. 1 A ii ). The first filter, with kernel [3, 1, -3], is used to enhance the adluminal edge of the top wall and the abluminal edge of the bottom wall ( Fig. 1 B for an explanation of the methodology followed to choose these kernel values). These filters provide a new vessel image with four enhanced edges, two on each wall, which facilitate the subsequent use of the cost function.
Step 2. Start and end points are automatically estimated by finding the darkest regions at both ends of the vessel (averaging over five columns; Fig. 1 B). Using these start and end points as initial parameters, the cost function is then used to draw lines along the four enhanced edges obtained in step 1, giving two lines per wall ( Fig. 1 C, yellow lines).
Step 3. Lines are drawn along the centre of each wall using the four lines detected in step 2. To do this, a number of points per wall (dependent on the length of the vessel under analysis) are selected along the length of the vessel and the averages between the two lines corresponding to each wall at those points calculated. Joining these averaged values provides two new lines, corresponding to the detected vessel walls (Fig. 1 D).
Step 4. Validation of the two lines obtained in step 3 is carried out with the assistance of a binary image. The binary image is obtained using an initial Gaussian filter (with sigma value of 2), before applying an automated threshold determined using Otsu's method [29], dilating the resulting image twice with a 363 kernel and filling holes, both in the normal and the inverted images ( Fig. 1  E). As can be seen in Figure 1, the image obtained by applying thresholding methods is not good enough to unambiguously determine the centre of the vessel walls with certainty, as several other structures also appear as part of the mask (white regions of the binary image). This image, however, is good enough to check whether the already obtained lines fall within the mask, as an additional validation.
Step 5. After the lines are validated (i.e. they were contained in the mask), they are further adjusted using a modified version of the binary image from step 4. This modified image is obtained by further dilating the mask and by applying a Gaussian filter to the image ( Fig. 1 F). Each pixel in the two estimated lines is then shifted up or down until the lightest intensity values present in the new smoothed binary image are reached, bringing them closer to the centre of the wall in those cases where they might have been drawn very close to the edge.
Step 6. The final lines are used to measure the diameter of the vessel (Fig. 1 G). This final measurement is obtained by averaging the vertical distances between the two lines over all pixels from left to right.
To estimate the computational complexity of the algorithm, the image to be processed is assumed to have a size M * N, and a kernel with size (2R+1) * (2R+1). The computational complexity can then be obtained for each step of the algorithm as follows: (1) Step 1: O(MNR 2 ). (2) Step 2: O(MN). (3) Step 3: O(S) (S -sample points along the vessel). (4) Step 4: O(MN(2*9+P 2 )+L 4 ) (P is the kernel size of the Gaussian filter, and L is the number of gray levels). (5) Step 5: O(MNP 2 ). (6) Step 6: O(Q) (Q is the number of pixels on one line). Finally, the total computational complexity of the proposed algorithm can be expressed as: The algorithm implements automatic restarting and fitting of parameters when needed (for steps 2 to 4). Briefly, when the lines drawn in step 3 do not pass the validation test (step 4), several parameters are adjusted (up to a maximum value) to provide a better fit of the data and the detection is restarted (back to step 2). Figure 2 gives a lower level illustration of the fitting process. The following parameters are fitted by the program: Window. Divides the image into increasingly narrower horizontal subsections of width equal to that of the original image (ranging from 2 to 6).
Level. Divides the image into increasingly narrower vertical subsections of height equal to that of the original image (ranging from 1 to 5).
Smoother. Determines the number of points to be used in the drawing of the central lines at step 3. The initial number of points depends on the length of the vessel (a default distance between points of 25 pixels is initially set, with the maximum number of points limited to 40).
Pixel tolerance. Determines the maximum number of pixels outside the binary mask that can be allowed during validation at step 4 (given as a percentage of the overall length of the vessel; ranging from 1 to 50%).
The goal of these parameters is to increase the probability of finding subsections within the image where good lines can be detected. If the values for all the parameters are exhausted and no good lines are found, the image fails and no final measurement is provided.

Presentation of the Results
MyoTracker runs independently on each image/slice. After the last slice is analysed, the program optionally displays the detected vessel walls as an overlay on the image for verification if desired, and provides diameter values in a results table and/or plot. A   manual is also available, which describes the outputs and adjustable parameters used by the algorithm (File S1). However, fine-tuning of parameters is normally not necessary. All results presented in this paper were obtained using the default parameter values, unless otherwise stated.

Testing the Program
An ImageJ macro was developed for manual vessel analysis to assist validation of the MyoTracker algorithm. The macro allowed for tests to be carried out over a period of time so that tiredness did not have a negative impact on the accuracy of the measurements. The diameters of retinal arterioles from 102 test images were manually measured by 3 different researchers. For this, each user was required to draw two lines (using the segmented lines tool in ImageJ) per image along the centre of both top and bottom walls and the macro determined the distance between them. These measurements were subsequently compared to the automatic measurements obtained by MyoTracker (ImageJ version 1.48a).

Processing Times
Another macro was developed to assist with the analysis of the performance of the algorithm, i.e. the processing time required by the algorithm to analyse both individual images and videos containing different numbers of slices. The macro was run on all 102 images from the dataset and on several videos, and averaged times (over 10 different runs for each image and video) were obtained on a PC with an IntelH Core TM i7-2600 and 16 GB of RAM.

Statistics
To assess the agreement between automatic and manual measurements, the standard deviation of the measurement error (calculated as the difference between automatic and manual values) was determined as previously described [30]. Bland-Altman plots were also constructed as an additional method of assessment [31], [32], with the manual measurements used as the reference or 'gold standard' method [33]. The plots were performed using MedCalc for Windows, version 12.7.5.0 (Med-Calc Software, Ostend, Belgium). Summary data is expressed as means 6 SD. In all cases, the accepted significance level was set at 0.05. Figure 3 shows a sample of 8 images taken from a full dataset of 102 retinal arteriole images (Fig. 3 A). This subset illustrates the different widths, lengths, positions, focus, contrasts and intensities considered during testing. Each of these factors contributed to the overall complexity of the analysis. In Figure 3 there are also two binary images obtained from panels #8 and #45 (Fig. 3 B, top and bottom images, respectively). These give an indication of the structures that are detected when carrying out threshold analysis alone. The lines detected by MyoTracker are superimposed on the binary images (Fig. 3

B, top and bottom, yellow lines).
Manual and automatic measurements were carried out in the full dataset containing 102 images. The results are presented in Table 1. The mean values for manual and automatic measurements were 26.868.7 and 26.868.8 mm, respectively. The mean of the measurement error for the MyoTracker software versus the manual measurements was 0.02 mm, whereas its standard deviation was 0.52 mm. This value was lower than any of the standard deviations obtained by comparing each of the 3 manual measurements with the mean of the other 2 in each case (0.9, 0.88 and 0.76 mm).
Bland-Altman plots were also used to assess the agreement between measurements. Figure 4 shows these plots for both raw differences (Fig. 4 A) and also differences relative to the width of the vessels (Fig. 4 B). As can be seen, in almost all images errors were within 61.5 mm, or 65% of the width of the vessels. The only image with an error higher than 5% was image #8 (error = ,6%; Fig. 4 B). This value, however, was close to the best relative error seen when comparing all the manual measurements with each other for this image (5.33, 10.01 and 15.33%). Both plots show that the mean of the errors (solid lines in Fig. 4 A and B) was within the confidence intervals of an error of 0, indicating that the automatic measurements were not significantly

Automatic Measurements on Videos
To investigate the variability observed in the automatic measurements during video recordings, MyoTracker was used to measure the diameter of retinal arterioles in recordings containing around 100 images each. Only vessels under steady-state conditions and with no visually observable constrictions or dilations were used for this analysis. The values obtained in this analysis had standard deviations ranging from 0.03 to 0.15 mm (with a mean of 0.08 mm, n = 5). In all cases the variability remained below 1% of the width of the vessels (ranging from 0.24 to 0.66%, with a mean of 0.35%).

Analysis of Vessel Movements
Tests were carried out on videos showing changes in the diameter of arterioles, both constricting and dilating to different external stimuli. The top panel of Figure 5 shows the time course of initial dilation of an isolated and cannulated retinal arteriole during pressurisation from 0 to 40 mmHg (black arrow) and the subsequent development of myogenic response over a period of 7 minutes (Fig. 5 A). The vessel dilated from 28.2 to 30.9 mm (steady-state values) at a rate of over 6 mm per second before constricting by about 2 mm (Fig. 5 A). The time course of the constriction of a different arteriole to application of the vasoconstrictor peptide, endothelin-1 (10 nM), is shown in the middle panel (Fig. 5 B, black line). In this case, the diameter of the vessel was reduced from 41.7 to 34.1 mm (steady-state values) at an approximate rate of 0.57 mm per second. Finally, the time course of two fast constrictions and subsequent slow dilations to caffeine (10 mM) is shown in the panel at the bottom (Fig. 5 C, black lines). In the first application of caffeine, the diameter of the vessel was reduced from 24.6 to 15.3 mm at a rate of 2.17 mm per second. In all cases the program was able to track the changes.

Measuring other Kinds of Vessels
Apart from arterioles, other types of retinal vessels were also tested using the MyoTracker software. We initially tested isolated retinal venules and capillaries where the adluminal and abluminal edges of the vessel walls were less obvious due to the absence of smooth muscle cells. Figure 6 shows examples of isolated retinal venules (Fig. 6 A) and capillaries (Fig. 6 B) of different sizes. The detected lines are shown superimposed on top of each of the vessels (Fig. 6 A and B, yellow lines). As can be seen, the program is capable of detecting the walls of these vessels, generally ignoring other elements present in and around the vessels, such as blood cells. We also tested the program with larger bovine arteries (Fig. 6  C). As can be seen, Myotracker was also able to track the walls of these vessels (Fig. 6 C, yellow lines).

Processing Times
The processing times required by MyoTracker to run on each of the 102 images contained in the dataset were recorded and appear in Table 2. Each of the times represents the average of 10 independent runs of the program for each image. The average processing time was 19 ms, ranging from 7.7 to 100 ms. The variability observed between the times for the 10 runs on the same image was always below 6 ms (with mean standard deviation of 1.7 ms). In this dataset, there was no correlation between mean times and the size of the images (R = 0.35). It is expected, however, that processing time will increase with complex images, where the algorithm might need to iterate over some of the fitting parameters to find appropriate values in the analysis (see Methods). Analysis carried out on 'image stacks' (equivalent to video files) with increasing number of images showed that, in general, processing time grows linearly with the number of images (Table 3).

Discussion
Previous studies concerning myogenic mechanisms in small isolated arterioles (,50 mm in diameter) have, at least in part, been hampered by the lack of appropriate software for accurately automating the analysis of vessel diameter. We have presented   Table 1). An advantage of our software compared to other programs is that MyoTracker works in a fully automated way, with no need for the user to take any manual measurements, define any initial windows or mark the walls of the vessel for analysis prior to running the program. Moreover, this program measures and tracks the diameter of the vessel over the full vessel segment present in the image, as opposed to carrying out the measurements on a small window of interest. Tracking of vessel diameters was tested with videos of images taken during experiments and the program was able to track both fast and slow changes along the full segments of the vessels and over the full duration of the experiments (Fig. 5). Given these advantages, this software may also be of interest for use with larger arteries and veins and we have shown that MyoTracker works well with such vessels (Fig. 6). There is also increasing evidence that abluminal pericytes on small post-capillary venules and capillaries are capable of constricting to modulate vessel diameter, and thus, may contribute to the control of local tissue blood flow and capillary pressures [34], [35], [36], [37]. We have shown that MyoTracker is also able to detect the borders of these vessels and therefore may be useful in studying control mechanisms throughout the microvascular system (Fig. 6).
Despite its advantages, MyoTracker does also have some limitations. Firstly, it requires the vessel to be orientated horizontally within the field of view and that both ends of the vessel reach the ends of the image. Although images can be cropped and rotated to accomplish this, in our experience, vessels can occasionally change their angle during experimentation (e.g. during pressurisation), adding a small inclination-dependent error to the measurements in some, but not all slices. As explained in the user manual (File S1), additional parameters have been incorporated into the software to help circumvent this problem. Briefly, instead of calculating the distance between the vessel walls vertically, in videos where this problem exists the shortest distance between the detected lines can be calculated on a pixel-by-pixel basis using a distance transform. Another limitation is that MyoTracker yields diameter data from the centre of the vessel walls rather than the adluminal or abluminal edges. This design choice was motivated by the need to find an alternative detection method not dependent on thresholding algorithms. Although cost analysis and edge enhancement provided such an alternative, it was necessary to combine the adluminal and abluminal detected lines (i.e. use them to correct each other) to add consistency and reliability to the detection algorithm. However, given that the main interest in many experiments is the tracking of vessel diameters with the goal of identifying relative changes (i.e. constriction or dilations) over time, detecting the middle of the vessel walls can serve this purpose.
To conclude, we present here a program to facilitate the measurement and tracking of the diameter of small isolated retinal arterioles during myography experiments. Our program has been tested with images taken from experiments on isolated retinal arterioles, venules and capillaries of different sizes. However, it is also likely to work on vessels from other microvasculature tissues, as they all present similar detection problems. MyoTracker may also be used to study larger vessels when its features are considered advantageous when compared to other software packages.

Algorithm Availability
The plug-in implementation is included as supporting information, together with a user manual (File S1). The datasets of images used for this study are also provided (Dataset S1 and S2). Figure S1 Measurement errors for different parameters in the algorithm. The error is shown as the SD of the difference between automatic and manual measurements in each case. The automatic values were obtained by the algorithm with either, only cost function and no further correction, or cost function and 6 different corrections: no kernel (using the original image with fitting parameters as outlined in Fig. 2 [4,1,24]). The measurements were carried out in 101 of the 102 images in the dataset (image #8 was left out of the analysis because the detection failed with no kernel and with k1). The plot shows a significant improvement (decrease) in the values of the measurement errors with increasing correction (P,0.001 between only cost function and the last 3 kernels, k2, k3 and k4), up to a point where no further significant improvement was detected (P,0.01 between k2 and k3; P.0.05 between k3 and k4). Thus, k3 was selected as the default kernel for the algorithm in this study. Significance was estimated using One-way Repeated Number of slices included in the video under analysis. Each image was a copy of image #1 (of size 2796147 pixels) from the dataset in Table 2. b Averaged (mean 6 SD) processing times (in ms) required to automatically analyse each video from 10 independent runs. doi:10.1371/journal.pone.0091791.t003

Supporting Information
Measures ANOVA with Newman-Keuls Multiple Comparison post-hoc test. (TIF) File S1 MyoTracker's code and manual. The code folder contains both the source code (JAVA File) and the Executable Jar File for the plug-in. This allows the execution of the software within the FIJI (Image-J) environment, as explained in the manual. Copyright and GNU General Public License files are also included in this folder. The documentation folder includes the manual with instructions regarding the installation and use of the software. (ZIP) Dataset S1 First set of 51 images of retinal arterioles used for the analysis. (ZIP) Dataset S2 Second set of 51 images of retinal arterioles used for the analysis. (ZIP)