A fast method for the detection of vascular structure in images, based on the continuous wavelet transform with the Morlet wavelet having a low central frequency
Abstract
A manual measurement of blood vessels diameter is a conventional component of routine visual assessment of microcirculation, say, during optical capillaroscopy. However, many modern optical methods for blood flow measurements demand the reliable procedure for a fully automated detection of vessels and estimation of their diameter that is a challenging task. Specifically, if one measure the velocity of red blood cells by means of laser speckle imaging, then visual measurements become impossible, while the velocitybased estimation has their own limitations. One of promising approaches is based on fast switching of illumination type, but it drastically reduces the observation time, and hence, the achievable quality of images. In the present work we address this problem proposing an alternative method for the processing of noisy images of vascular structure, which extracts the mask denoting locations of vessels, based on the application of the continuous wavelet transform with the Morlet wavelet having small central frequencies. Such a method combines a reasonable accuracy with the possibility of fast direct implementation to images. Discussing the latter, we describe in details a new MATLAB program code realization for the CWT with the Morlet wavelet, which does not use loops completely replaced with elementbyelement operations that drastically reduces the computation time.
a]Eugene B. Postnikov b]Maria O. Tsoy b]Maxim A. Kurochkin b]Dmitry E. Postnov a]Kursk State University, Department of Theoretical Physics, Radishcheva st. 33, Kursk 305000, Russia b]Department of Optics and Biophotonics, Saratov State National Research University, Astrakhanskaya st. 83, Saratov 410012, Russia \authorinfoCorresponding author: E.B. Postnikov. Email: postnicov@gmail.com
1 Introduction
Being the conventional component of some routine methods of medical screening, say, during capillaroscopy [1], a measurement of blood vessels diameter becomes the key technique for everyone who needs to quantify the functioning of vascular networks. While laserspeckle imaging (LSI) generally shows only the relative changes of red blood cells (RBC) velocity [2], some other methods, like particle image visualisation (PIV) can deliver its absolute values. The computational scheme of both these methods implies the selection of regions of interest (ROI), to avoid the useless and computationally expensive processing of “empty” image segments. Here the segmentation of the micrivascular network together with the fast and reliable determination of margins of vessels in images is highly demanded. While there are the few published attempts to estimate stationary vessel diameter from laser speckle data [3], the problem becomes even more complicated if one needs to monitor the microcirculatory network dynamics, say in response to different stimuli [4]. One of the promising approaches is based on the fast switching of illumination type, but it drastically reduces the observation time, and hence, the quality of images. For LSI applications, recent attempt is known to avoid such problems [5].
The wavelet transform, which is a powerful tool for the signal and image processing being especially adjusted to the detection of multiscale patterns [6], attracted an attention in biomedical image processing even at first stages of its development [7]. Now the application of wavelet methods on parity with other image filtering approaches is widespread and included into standard mathematical software as MATLAB [8].
Being initially aimed for image compressing, these waveletbased methods are focused mainly on the discrete wavelet transform, which inherits classical filtering algorithms. In addition, there usually operate in the domain of real functions that gives multiscale expansions, which require an additional processing to reveal localized structures, see the reviews in books cited above. Some known attempts to utilize the discrete complex wavelet transform, which follow the pioneering work [9], meet certain difficulties operating with signals without the binary hierarchical subdivision. At the same time, the complex continuous wavelet transform is more flexible in the sense of arbitrary continual spatial shifts and scales but conventionally used mainly for the detection of periodic pattern but not to identify the individual spots and curves.
In our work, we explore how the spatialdomain continuous wavelet transform can be used on the playground of the pattern recognition typical for the processing of images of microcirculatory network. Specifically, we apply it for the automated determination of a “recommendedtoprocess” part of both the surrogate spatial pattern and for the real vascular network image obtained from the chorioallantoic membrane (CAM) of chicken embryo. The obtained results show the applicability of suggested approach and suggest the way for further its development.
2 Method
2.1 Premises of the approach
The integral representing the continuous wavelet transform (CWT),
(1) 
has a direct mathematical interpretation as a correlation between the analysed function and the wavelet (the asterisk in (1) denotes a complex conjugation) calculated for the given shift and scale . Thus, its maxima correspond to shapes of the function’s features, which locally resemble the wavelet’s shape in a best way.
Among a variety of wavelet functions, the standard Morlet wavelet defined in the amplitude norm as
(2) 
is one of the most popular for the local spectral analysis, i.e. for revealing local periodicities. Such an application is based on the shape of this wavelet (2), which is the harmonic function modulated by the Gaussian. It exhibits a number of oscillations for the conventionally used central frequencies , see Fig. 1(A). These large values are based on the socalled admissibility condition , which is satisfied for (2) only approximately: .
On the other hand, it has been shown recently [10, 11] that the admissibility condition is not a necessary property, and its violation does not prevent an exact invertibility of the transform. Therefore, we can feel free to consider smaller values of the central frequency, which result in low and nonoscillating shapes of the wavelet (2).
I should be pointed out that the complex conjugated expression (2) can be expanded in the Tailor series and truncated for small up to the first terms as
(3) 
One can see that the real part of (3) being substituted into Eq. (1) provides a Gaussian (diffusion) smoothing of the processed function , i.e. its denoising; and the imaginary part evaluates differentiation of the smoothed signal, i.e. edge detection. See Fig. 1(C) for a typical example of the shape of these sliding windows.
At the same time, the basic property of the CWT with the Morlet wavelet applied to the pure harmonic function, is the localisation of the signal’s frequency as a maximum of the Gaussian representing CWT’s modulus with respect to the scale:
This property is fulfilled for arbitrary . Note also that the series of wavelet coefficients corresponding to exactly coincides with a magnitude of the analysed signal as well as a phase factor of the CWT. Thus, the determination of automatically provides values of the desired features too.
Therefore, the CWT with small central frequencies simultaneously denoises localized maxima of the signal bounded by some edges and allows for determining their scales simply detecting CWT modulus maxima.
2.2 Algorithm of image processing
The proposed algorithm is formulated as follows:

determine number of pixels in the imageâs row, which give the scaling of frequencies of the Fourier transform (here we use the positive frequencies only, whence rounds towards plus infinity);

to fix a small central frequency, e.g. and an appropriate matrices of the proper scales and Gaussian factors ;

for each row

i) to find the Fourier coefficients of the deviations of pixels intensities from its average value in a row via the Fast Fourier Transform (FFT); compute their elementwise Gaussian filtering , and the inverse Fourier transform of the results doubled to take into account only half of frequencies used; this results in a matrix of the CWT coefficients , .

ii) to find a set of maxima for each array with a fixed and to reconstruct the envelopes of vasa’s cross section pictures by the summation and taking of real part of the result;

the obtained image should be binarised with respect zero (or small threshold above zero) level; the positive values will be mapped into 1’s, which provide the desired mask identifying vessels locations.
For better accuracy, the procedure can be repeated along columns of digital images (in this case vasa situated more horizontally will be detected; they may be missed during the rowbyrow processing) and both processed images combined to obtain a more accurate mask.
2.3 Speeding up CWT computation with MATLAB
It should be pointed out that the algorithm described above includes a number of stringbystring and rowbyrow loops being applied to 2D images processing. This fact naturally induces its practical realization using software utilizing internal matrixbased methods, which provides loopless fast computations, fro example MATLAB can be considered in such role.
We propose new realization of the CWT with the Morlet wavelet, which replaces for
loopbased algorithms utilized in the function included into the conventional packages such as Wavelet Toolbox by MathWorks and WaveLab by the Stanford University.
Then the program code (and comments marked by % sign, which describe all variables) reads as follows:
% Input parameters: % u  input data (a row array) % N  length of u % omega0  the central frequency's value % a  scale values % Na  length of a % Output parameter: % w  the resulting complexvalues matrix of the wavelet trasform of the size N*Na L=ceil(N/2);%1 omega_=2*pi*(0:L)/(N1);%2 AA=repmat(a',1,L+1);%3 OMEGA2=repmat(omega_,Na,1);%4 OMEGAs2=OMEGA2.*AA;%5 WINDOW2=exp((OMEGAs2omega0).^2/2);%6 F=fft(u);%7 FF2=repmat(2*F(1:L+1),Na,1);%8 CNV2=WINDOW2.*FF2;%9 w=ifft(CNV2,N,2);%10
The line %2
form a row array of frequencies, for which fft
function gives the Fourier coefficients; it is converted into 2D matrix, which contains L+1
rows and Na
columns in the line %4
by repetition of columns containing frequency values. The same procedure generates 2D matrix of scales of the same shape (the line %3
).
As a result, the set of scaled frequencies and windows (the lines %5
and %6
)
are formed by the elementbyelement multiplication (marked as .*) instead of any timeconsuming loop. This part of code is the same for all rows of an image, therefore, can be calculated one time and placed in the beginning of the full program.
The second part (the lines %7
–%10
)
directly calculates the CWT using the predefined common parameters. These lines of code contain the Fourier transform via FFT algorithm (%7
),
forming 2D matrix corresponding to its positive frequencies with the number of columns equal to the number of used scales (%8
),
sliding Gaussian filtering of the spectrum (%10
)
realized as an elementbyelement multiplication, and, finally, the inverse Fourier transform. Note that we kept two times reduced size of all matrices up to this moment. However, the special feature of MATLAB’s ifft
function provides an opportunity not only to invert the the whole matrix by columns but also make them padded with zeros up to the required number N
.
3 Results and their discussion
For the testing purposes we first created the surrogate image being composed of (i) a random pattern, (ii) stripes of different width with reduced brightness, and (iii) an intensity gradient applied from top to bottom. Both the surrogate image and the processing results are illustrated in Fig. 2. Obviously, the simple thresholdbased method to extract the stripes location and shape, which is shown in panel (b), fails due to the presence of gradient of overall brightness, which leads to the complete loss of a bottom part of the image. The waveletbased technique (c) shows much better result in this respect. Note, we intentionally did not provided any smoothing or interpolation to reduce the hatched structure which is visible of the image in panel (c) in order to show raw result of the transform. Such smoothing can be done during the subsequent composition with results of vertical scan.
In order to apply the method to the real microcirculatory network, we use the image obtained from CAM of 11day old chicken embryo. The preparation procedure is quite simple and can be found everywhere. The opened egg shell was fixed vertically and illuminated with green light. The optical image was recorded in gray scale. The arbitrary chosen single frame is shown in Fig. 3(left panel). It was processed with the algorithm described above with the central frequency , see the corresponding wavelet shape in Fig. 1(B). Although this value does not correspond to a single bellshaped curve like Fig. 1(C), the practically one full oscillation presented there provides an effective way to detect a localised bright regions obtained by the cross section of the image. A series of numerical experiments with different have demonstrated an effectiveness of such a choice for the present ranges of vessels diameters and spatial noise scales. It can be explained if one note that the chosen shape is quite close to the classical Laplacian filter [8] and its derivative. Thus, there is actually a kind of combination of Laplacian denoising and the zerocrossing edge detection algorithm realized simultaneously since the wavelet is a complex function.
Fig. 3(right panel) shows the obtained mask (binarized image). One can see the different degrees of visualization depending on vessel size. This is a sign of the promising method features, not yet fully implemented. Also note, that not all the vessels are mapped in proportionally reasonable sizes. Definitely, all this is dependent on the vessel orientation and, thus, the corresponding method development is in order.
4 Conclusion and outlooks
In conclusion, we have shown a new promising way for the fast and computationally effective processing of image data obtained from the microscopy of microcirculatory network of chicken embrio. Esentially the same approach can be used in the case of processing not just optical images but some calculated distributions of quantities relevant to the measured microcirculation dymanics. We believe that the capability of fast and reliable processing will be demanded in many applications. Specificaly, we intend the further use of the described technique for the purposes of pattern separation during the estivation of dynamic changes in vessel barrier functions. While in specific vascular networks like brain or retina, vessels show high protections from penetration of different chemicals, in available biomodels such leakage can be considerable, and, thus, more easy to be quantified, but at the same time need more precize determination of vessels locations and boundaries. The proposed technique is expected to show good results for this task.
Acknowledgements
This work was funded by Russian Science Foundation, project 161510252.
References
 Maricq, H. R. and Carwile LeRoy, E., “Patterns of finger capillary abnormalities in connective tissue disease by âwidefieldâ microscopy,” Arthritis and Rheumatism 16(5), 619–628 (1973).
 Tuchin, V. V., ed., [Handbook of Coherentdomain Optical Methods: Biomedical Diagnosis, Environmental Monitoring, and Materials Science ], Springer (2013).
 Liu, Q., Li, Y., Lu, H., and Tong, S., “Realtime high resolution laser speckle imaging of cerebral vascular changes in a rodent photothrombosis model,” Biomed. Opt. Express 5, 1483â1493 (2014).
 Neganova, A., Postnov, D., BringsJacobsen, J., and Sosnovtseva, O., “Laser speckle analysis of retinal vascular dynamics,” Biomed. Opt. Express 7, 1375–1384 (2016).
 Postnov, D. D., Tuchin, V. V., and Sosnovtseva, O., “Estimation of vessel diameter and blood flow dynamics from laser speckle images,” Biomedical Optics Express 7, 2759–2768 (2016).
 Mallat, S., [A wavelet tour of signal processing ], Academic press (1999).
 Aldroubi, A. and Unser, M., eds., [Wavelets in medicine and biology ], CRC press (1996).
 Semmlow, J. L. and Griffel, B., [Biosignal and Medical Image Processing ], CRC Press (2014).
 Kingsbury, N., “Image processing with complex wavelets,” Philos. Trans. Royal Soc. A: Math., Phys. Eng. Sci. 357(1760), 2543–2560 (1999).
 Lebedeva, E. A. and Postnikov, E. B., “On alternative wavelet reconstruction formula: a case study of approximate wavelets,” Royal Soc. Open Sci. 1, 140124 (2014).
 Postnikov, E. B., Lebedeva, E. A., and Lavrova, A. I., “Computational implementation of the inverse continuous wavelet transform without a requirement of the admissibility condition,” Appl. Math. Comput. 282, 128–136 (2016).