PLEASE SCROLL DOWN FOR ARTICLE
On: 4 January 2011
Access details: Access Details: [subscription number 907850282]
Publisher Taylor & Francis
Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House,
37-41 Mortimer Street, London W1T 3JH, UK
Computer Methods in Biomechanics and Biomedical Engineering
Publication details, including instructions for authors and subscription information: http://www.informaworld.com/smpp/title~content=t713455284
An approach to automatic blood vessel image registration of
microcirculation for blood flow analysis on nude mice
Wen-Chen Lina; Chih-Chieh Wua; Geoffrey Zhangb; Tung-Hsin Wuc; Yang-Hsien Lind; Tzung-Chi
Huangd; Ren-Shyan Liue; Kang-Ping Lina
a Department of Electrical Engineering, Chung Yuan Christian University, Chungli, Taiwan b
Department of Radiation Oncology, Moffitt Cancer Center, Tampa, FL, USA c Department of
Biomedical Imaging and Radiological Science, National Yang-Ming University, Taipei, Taiwan d
Department of Biomedical Imaging and Radiological Science, China Medical University, Taichung, Taiwan e Department of Nuclear Medicine, Taipei Veterans General Hospital, Taipei, Taiwan
First published on: 15 November 2010
To cite this Article Lin, Wen-Chen , Wu, Chih-Chieh , Zhang, Geoffrey , Wu, Tung-Hsin , Lin, Yang-Hsien , Huang, Tzung-Chi , Liu, Ren-Shyan and Lin, Kang-Ping(2010) 'An approach to automatic blood vessel image registration of microcirculation for blood flow analysis on nude mice', Computer Methods in Biomechanics and Biomedical Engineering,, First published on: 15 November 2010 (iFirst)
To link to this Article: DOI: 10.1080/10255842.2010.497489 URL: http://dx.doi.org/10.1080/10255842.2010.497489
Full terms and conditions of use: http://www.informaworld.com/terms-and-conditions-of-access.pdf
This article may be used for research, teaching and private study purposes. Any substantial or systematic reproduction, re-distribution, re-selling, loan or sub-licensing, systematic supply or distribution in any form to anyone is expressly forbidden.
The publisher does not give any warranty express or implied or make any representation that the contents will be complete or accurate or up to date. The accuracy of any instructions, formulae and drug doses should be independently verified with primary sources. The publisher shall not be liable for any loss, actions, claims, proceedings, demand or costs or damages whatsoever or howsoever caused arising directly or indirectly in connection with or arising out of the use of this material.
An approach to automatic blood vessel image registration of microcirculation for blood flow
analysis on nude mice
Wen-Chen Lina, Chih-Chieh Wua, Geoffrey Zhangb, Tung-Hsin Wuc, Yang-Hsien Lind, Tzung-Chi Huangd*, Ren-Shyan Liueand Kang-Ping Lina
a
Department of Electrical Engineering, Chung Yuan Christian University, Chungli, Taiwan;bDepartment of Radiation Oncology, Moffitt Cancer Center, Tampa, FL, USA;cDepartment of Biomedical Imaging and Radiological Science, National Yang-Ming University, Taipei, Taiwan;dDepartment of Biomedical Imaging and Radiological Science, China Medical University, Taichung, Taiwan;eDepartment of
Nuclear Medicine, Taipei Veterans General Hospital, Taipei, Taiwan (Received 23 March 2010; final version received 27 May 2010)
Image registration is often a required and a time-consuming step in blood flow analysis of large microscopic video sequences in vivo. In order to obtain stable images for blood flow analysis, frame-to-frame image matching as a preprocessing step is a solution to the problem of movement during image acquisition. In this paper, microscopic system analysis without fluorescent labelling is performed to provide precise and continuous quantitative data of blood flow rate in individual microvessels of nude mice. The performance properties of several matching metrics are evaluated through simulated image registrations. An automatic image registration programme based on Powell’s optimisation search method with low calculation redundancy was implemented. The matching method by variance of ratio is computationally efficient and improves the registration robustness and accuracy in practical application of microcirculation registration. The presented registration method shows acceptable results in close requisition to analyse red blood cell velocities, confirming the scientific potential of the system in blood flow analysis.
Keywords: microcirculation; blood flow; image registration
1. Introduction
Microcirculation is crucial in maintaining tissues and organs, due to delivering oxygen and nutritious materials. Nowadays, it is well known that microcirculation has a lot to do with every kind of critical disease (Mchedlishvili 1998). In order to identify its characteristics and its relationship to each disease, it is essential to measure the blood flow accurately. Furthermore, microcirculatory study is important for understanding the amount of the blood supplied in the local tissue, such as in the tumour of a laboratory mouse. The examination of tumour microvascular mechanisms, such as quantification of tumour vascular response by monitoring the blood flow information at each local capillary, is an effective way to identify specific tumours, and has been studied extensively in mouse models (Jain and Ward-Hartley 1984; Hori et al. 1991). Several studies on tumour microcirculation by intravital video microscopy (IVVM) have been presented to allow dynamic observation of red blood cell (RBC) activity in the microcirculation of intact organs and tissues in live animals (Sugii et al. 2002; Iga et al. 2006).
Large serial sets of microscopic images are usually required for analysing microcirculatory videos in vivo in which various imaging modalities provide different and unique information. There are several kinds of methods of
RBC velocity estimation using space – time diagrams. They have been comprehensively investigated both experimen-tally and theoretically (Jeong et al. 2006; Dobbe et al. 2008). However, movement of the subject caused by regular heartbeats can result in unstable images that hamper the vessel recognition and velocity measurements. In order to obtain stable images for blood flow analysis, frame-to-frame image matching as a preprocessing step is an effective solution to the motion problem. Furthermore, image registration is a prerequisite and time-consuming step. A typical method for solving the registration problem is usually a mathematical optimisation process, using a cost (or similarity) function to quantify the quality of the alignment of the two images for given geometric transformations including translation, affine and projection. In practice, only the translation component of motion was determined neglecting any rotation effects for a continuous series of microscopic vessel images. It usually relies on the alignment in the horizontal and vertical directions for which the use of an optimisation method is crucial for obtaining automatic, rapid, robust and precise registrations. Presently, most attention has been focused on some aspects of the problem, such as defining measure of match or objective functions and certain search methods (Maes et al. 1999; Matsopoulos et al. 1999, 2003; Jenkinson and Smith 2001).
ISSN 1025-5842 print/ISSN 1476-8259 online
q2010 Taylor & Francis
DOI: 10.1080/10255842.2010.497489 http://www.informaworld.com
*Corresponding author. Email: tzungchi.huang@mail.cmu.edu.tw iFirst article, 2010, 1–12
The metric evaluates the similarity between the two images. Several image similarity measures have been proposed. They can be classified depending on the theoretic foundations:
. Correlation measures. The intensity values of each image are analysed and the alignment is achieved when a certain correlation measure is minimised (or maximised). Usually, a priori information is used in these fast metrics.
. Intensity occurrence. These measures depend on the probability of each intensity value and are based on information theory.
. Geometrical features. A segmentation process detects some features and then they are aligned. These methods obtain volume accuracy and they do not have high computational cost. Nevertheless, there is a great dependence on the initial segmentation results.
The optimiser searches the maximum (or minimum) value as the metric varies the spatial transformation. A lot of numerical methods have been developed in order to obtain the global extreme of a non-analytical function. The following methods are the most used in medical image registration: Powell’s method, gradient descent, one-plus-one evolutionary and the simplex method. The choice of method will depend on the implementation criteria and the measure features (smoothness, robustness, etc.).
In this paper, an applicable automatic registration scheme is proposed and applied in the evaluation of existent metric for the capillary blood flow video sequence in tumour regions of nude mice with severe combined immune deficiency (SCID). We allow employing a strategy with highly accurate and rapid registration, with its objective function based on the feature-based or intensity-based matching technique. In addition, to examine the blood cell flow velocity in each local capillary, the microscopic system without fluorescent labelling provides precise and continuous quantitative data of blood flow rate in individual small vessels.
2. Materials
2.1 Animal origin and tumour
Female nude mice with SCID, 7 – 12 weeks old and with an average weight of 16 – 17 g, were used for tumour blood flow measurements and vital microscopic observations under anaesthesia (Figure 1). After tumour cells were implanted and grown in the back leg muscle, the nude mice were imaged for capillary blood flow in the tumours using the M320 system (JMC Corporation, Kyoto, Japan).
2.2 Microscopy and image acquisition
Many different blood flow measurement techniques have been proposed (Sugii et al. 2002; Iga et al. 2006; Jeong
et al. 2006). The electromagnetic blood flowmeter and ultrasonic Doppler flowmeter are useful for larger blood vessel measurement, but not suitable for microcirculation. In general, IVVM involves the use of a fluorescence microscope in a living animal for real-time observing, monitoring and recording, and can be applied in the quantitative analysis of specific variables and events. However, some detrimental effects as given may arise with fluorescence microscopy on living tissues, below. (1) The fluorophore itself may interfere with the signalling
pathway or alter cellular function in some way. (2) The excitation light itself may damage the living
tissue, which may affect the behaviour of the sample or even cause its death.
(3) The effect arises as a result of the combination of fluorophores, and excitation light is clearly not an option with live cell imaging.
Our microscopic system without fluorescent labelling provides precise and continuous quantitative data of blood flow rate in individual small vessels. The wavelength and illumination are key parameters of the light source in a microcirculation imaging system. The penetration capa-bility of light depends on the filter wavelength. Generally, a longer wavelength will result in deeper penetration under skin but with worse resolution (see Equation (1)). Therefore, the wavelength of the light source should be chosen properly to obtain correct and clear images which can be used to analyse blood flow or morphologic features accurately. Absorption of blood under the skin is another important consideration that will influence the contrast or clarity of the image. According to the absorption spectrum of haemoglobin, the 532-nm wavelength is strongly absorbed by oxygenated haemoglobin, giving it a superior probing ability for RBCs in capillaries. Theoretically, a green light source (500 – 550 nm) is more suitable for
(a) (b)
Figure 1. (a) Nude mouse with SCID. Tumour cells implantation is within the red circles, (b) the microcirculation in the tumour region observed after an easy surgical treatment of surface.
achieving good image quality: rAiry¼ 1:22
l
2 £ NA; ð1Þ
wherelis the wavelength of light and NA is the numerical aperture of the objective lens (Wu et al. 2007).
The technique of window scanning is used in this application. Commonly used techniques such as horizontal, vertical or diagonal pixels scan are much more efficient for sampling spatial information on the surface of the tumour target. A comparison of the two scanning techniques, horizontal and vertical, is shown in Figure 2. Both scanning techniques move on the constructed geometry signature with a distance of 1 mm per line. The geometry of the scanning is able to gather the amount of data needed to achieve retrieval.
2.3 Capillaries microcirculation on nude mice
Microcirculation includes arterioles, capillaries and venules; the capillaries lie or connect between the arterioles and venules. Capillaries form extensive branching net-works, which are the smallest and most numerous blood vessels in the body, dramatically increasing the surface area available for the rapid exchange of molecules. In microcirculation, true capillaries are thin vessels which allow one RBC to pass. Pre-capillaries are vessels lacking complete coats, intermediate just to the arterial side of a capillary, and are about 3 – 5 times the diameter of a capillary. They not much different from a capillary in other aspects. Moreover, pre-capillaries and capillaries branch off from metarterioles and terminal arterioles. Within medical anatomical considerations of microcirculation, capillary microcirculation in the skin of mice can be observed and displayed using the microscopic system without fluorescent labelling, as shown in Figure 3.
3. Registration methodology
In this study, image registration algorithms of intensity-based and feature-intensity-based matching techniques are
evaluated. The registration process is divided into three major steps, as illustrated in Figure 4. The first step is image preprocessing. In this step, a clear image is chosen as a reference by quality assessment, and the preprocessing property and down-sampling are used to produce images from fine to coarse resolution for the floating images and reference image. The second step is a coarse initial alignment using a fast optimisation search which is based on many values for similar images in video sequences. In order to achieve a highly accurate registered image, the final step is to fine-tune the alignment between the original reference and floating images according to the initial values.
3.1 Image preprocessing of video sequence 3.1.1 Quality assessment of reference image
The quality of the reference image is essential for the accuracy of image registration. Selecting good reference images from a sequence of frames can efficiently depress the false rejection rate and false acceptance rate of the recognition system. However, the evaluation function of image quality should have several characteristics:
(a) (b) (c)
Figure 2. Microscopic system for lighting and a sampling scheme. (a) M320 (JMC Corporation, Kyoto, Japan) using a closed-circuit video system consisting of a CCD video camera, an LCD monitor and a video recorder. The video recorder acquires real-time blood flow video of microcirculation at £ 380 magnification, with a special resolution of 1.42 mm and a image sampling rate of 30 frames per second, and the static image pixel matrix size is 720 £ 480. Two scanning directions on the tumour target are shown in this illustration, (b) horizontal and (c) vertical directions.
(a) (b) (c)
(d) (e) (f)
Figure 3. The microcirculation in capillaries on nude mice is observed as shown in this illustration. Various patterns of micro-vessels are shown, such as (a) breathing micro-vessels, (b) capillary, (c) pre-capillary and (d) – (f) complex capillaries.
(a) agonic, (b) single peak, (c) high sensitivity, (d) high SNR and (e) low computation. The common evaluation methods for focusing are as follows: (Winkler 1999; Wang et al. 2004; Wei et al. 2005):
(1) High frequency component. Perform the fast Fourier transform on an image and extract the high-frequency component as the focus evaluation, which has a large computation cost and is not suitable for real-time measurement.
(2) Smoothness. The smoothness between the adjacent pixels is calculated. The smoothed pixels are summed as the focus evaluation function. The maximum value is the focus criterion. This method has low sensitivity.
(3) Threshold integration method. Choose a threshold according to the grey distribution. Sum the grey values above the threshold. The maximum value is the focus criterion. This method has low precision. (4) Grey difference method. This method makes use of the absolute value of adjacent pixels as the focus criterion function. To each pixel of an image, the step summation on the x-axis and y-axis is calculated in
Equation (2), DðiÞ ¼X ðx;yÞ {j fiðx; yÞ 2 fiðx; y 2 1Þj þ j fiðx; yÞ 2 fiðx 2 1; yÞj}; where i ¼ 0; 1; . . . ; m and m [ Z: ð2Þ
When the image is focused, i ¼ 0; 1; . . . ; m, D(i) is the maximum value. From observation, motion blur is mainly caused by factor-oriented abnormity which affects the image quality. Since motion blur occurs when the pixels overlap along the moving direction, it leads to the grey level difference between neighbouring points getting smaller. So the grey level difference can be used to evaluate the motion blur. We can then choose the clearest image as reference by searching this maximum of D(i). An example is shown in Figure 5.
3.1.2 Down-sampling and preprocessing property Usually, this is a common approach for registration. The original images are down-sampled to speed up the registration convergence. For instance, in this study, the image pixel matrix is up to 720 £ 480 pixels. With a scale factor of two, an image is down-sampled to about 360 £ 240 pixels prior to the processing stage. Sometimes, the Gaussian filter is employed as the preprocessing property to smooth the down-sampled images, and it can improve the registration performance and avoid local extremes using low pass smoothing.
3.2 Measurement of metric
Serial image sequences of the same modality are acquired, and they have the same resolution in space and grey scale. The highest-quality image is selected as the reference (denoted with R) whereas the others are referred to as floating (denoted as F). The floating images are aligned with the reference image in space by translating the floating image. In general, the alignment is solved either manually or by computer-based (semi-)automatic registration techniques. Manual registration is shown to be error-prone, time consuming and probably biased. Computer-based methods use either intensity-based registration or feature-based registration. Intensity-based registration with different similarity metrics has been generally adopted; these include correlation measures and intensity occurrence. Feature-based matching techniques do not use the grey values to describe matching entities, but use image features derived by a feature extraction algorithm. The form of the description as well as the type of features used for matching depends on the task to be solved. Evaluation of registration strategies for large serial sets of microscopic images is strongly related to
Image quality assessment (Microcirculatory video sequence) Measurement of matching Measurement of matching Geometric transformation Geometric transformation
Optimisation local search by Powell’s method Slight search by direct method Registration results (a) (b) (c) Reference image Uniform random generator Image pre-processing Image pre-processing Down-sampling Down-sampling
Figure 4. Typical block diagram of automatic registration. (a) Assessment of reference image and preprocessing applied to video sequences. (b) Coarse and fast alignment by an optimisation search. (c) Precise alignment by fine searching with the search range depending on the specific ratio of down-sampling.
the work of measurement of metric. A brief introduction of matching metric is summarised as follows:
3.2.1 Correlation measures
. Sum of absolute differences (SAD). It uses the absolute summed grey-level differences between R and F images to minimise the distance given by the following equation:
SAD ¼X x;y
jrðx; yÞ 2 f ðx; yÞj; ð3Þ
where rðx; yÞ and f ðx; yÞ are the grey levels of reference and floating images, respectively, and the sum is taken over the size of the window.
. Variance of ratio (VR). If rðx; yÞ is the grey value of a pixel at (x,y) in the reference image and f ðx; yÞ is the value of the corresponding pixel in the floating image, let ratioðx; yÞ ¼ rðx; yÞ=f ðx; yÞ. The regis-tration technique is based on the assumption that ratioðx; yÞ is maximally uniform across pixels when the two images are accurately registered. Ifdratiois the standard deviation of ratioðx; yÞ over all templates within the reference and mratio is the mean value of ratioðx; yÞ over the template within the floating image, the algorithm uses dratio=mratio as a measure of how well the two image sets are
registered. A detailed description of this objective function is provided in the literature (Woods et al. 1992).
. Normalised cross-correlation (NCC). Using image correlation to measure the similarly of matching between the reference frame and the floating frame, the cross-correlation function equation is described as follows:
where fmðx; yÞ and gmðx; yÞ indicate the intensity at coordinates R and F in the window in successive frames, f and g show the average intensity in windows, (k,l) is the moving vector of the window between frames and N is the number of pixels in the window. By scanning on the second image for the highest point of the ratio, the displacement volume can be determined. In this method, measurement based on image intensity was made on a pixel basis and results could only be obtained under scattering conditions. Continuous results can be obtained by formulating correlation values between pixels with the maximum ratio and adjacent pixels using the quadratic equation. The coordinate that has the maximum value is regarded as the one with maximum correlation. A detailed description of this objective function is provided in the literature (Avants et al. 2008). 0 50 100 150 200 250 300 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
Index of video sequence i
Sun of difference D(i)
Evaluation of image quality
(a) (b) (c)
Figure 5. (a) Image quality in video sequence is evaluated, and (b) then the best clear image with the highest value is chosen as the reference image (c) images with low values are defocused or motion blurred ones.
CCORðk; lÞ ¼ ð1=NÞPN21 m¼0ðfmðx; yÞ 2 fÞðgmðx þ k; y þ lÞ 2 gÞ ð1=NÞPN21m¼0ðfmðx; yÞ 2 fÞ2ð1=NÞP N21 m¼0ðgmðx þ k; y þ lÞ 2 gÞ2 1=2; ð4Þ
3.2.2 Intensity occurrence
. Mutual information (MI). MI (Huang et al. 2003) is an information theoretic measure and was proposed for use in image registration by two independent groups. The basic concept behind the use of this measure is to find a transformation. The MI of two images R and F is then defined by summing all grey value pairs at corresponding positions, regarded as random vari-ables X and Y and their intensity values at a certain coordinate in the images as the joint outcome of a random experiment. The MI between X and Y, denoted as IðX; YÞ, is a measure for the statistical dependency between both variables. The MI between two random variables X and Y is defined as
IðX; YÞ ¼ HðXÞ 2 HðXjYÞ ¼ HðYÞ 2 HðYjXÞ ¼X
x[X X y[Y
pðx; yÞ log pðx; yÞ
pðxÞpðyÞ: ð5Þ MI represents the amount of information that one random variable, the output of the channel, gives (or contains) about a second random variable, the input of the channel, how much the knowledge of X decreases the uncertainty of Y. Therefore, I(X,Y) is a measure of the shared information between X and Y.
3.2.3 Geometrical features matching
Feature-based matching is a powerful method for medical image registration, such as retinal images. A wide range of operators for image enhancement and detection has been described in the literature (Matsopoulos et al. 1999). In this similar case of vascular images, after preprocessing of feature extraction two binary images of the vascular vessels are obtained, BR and BF. The coarse alignment is based on localised correlation of binary templates in binary images. The calculation of the correlation is performed using an XOR operation between BRand BFdata in order to find all pixels set to different values. Here, XOR is tested by developing an objective function defined as follows:
M ¼ arg maxXORðBR; TðBFÞÞ
H £ W ; ð6Þ
where H £ W is the size of the overlap. The search pattern is moved over the binary image by the search method and the correlation can be mathematically formulated as the maximisation. The number of bits agreeing at each location is taken as a final measure of the result.
3.3 Optimisation search
Typically search algorithms may use either a direct method or a mathematical method. The use of the former is
restricted almost completely to applications relying on very sparse information. The mathematical search methods that exist are suitable for global or local optimisation (Matsopoulos et al. 2003). Some global optimisation methods exist (such as genetic algorithms and simulated annealing). Many of the methods require a very large number of iterations to satisfy statistical convergence criteria. It is, therefore, the speed of the global optimisation methods that is the limiting factor for their successful application to this problem (Jenkinson and Smith 2001). Although some existing global optimisation techniques may be viable, the majority of existing registration methods opt for local optimisation methods, such as Powell’s direction set method (PDSM), gradient descent, downhill simplex methods and so on, in order to increase speed. PDSM was found to often yield the best optimal solution among these methods (Xu and Dony 2004). It only requires evaluations of the cost function itself and not of its gradient. Despite its strengths, its performance is found to be strongly dependent on the initial solution, and the global solution is not guaranteed. Generally, this drawback can be improved using a combined simplex (Plattard et al. 2000), a method hybrid genetic scheme (Hsiao et al. 2001) and a multi-resolution algorithm (Jenkinson and Smith 2001) to solve the desired global minimum, which is hidden among many local minima.
Powell’s direction set algorithm is selected to perform the optimal search in this study. It is quadratically convergent, which makes it one of the fastest methods for solving non-linear minimisation problems. A key idea of Powell’s method is to select a starting point Pðp1; . . . ; pnÞ in parameter space with a uniform random generator, which avoids getting stuck in local minima as much as possible, and minimise the cost function along some chosen direction vector n using a 1-D minimiser such as bracketing or Brent’s parabolic method. The most appreciated scheme for choosing successive directions is the conjugate gradient method, which we now briefly describe. When the cost function is minimised along the direction n, the gradient in the obtained minimum is necessarily perpendicular to that direction. If it were not, that would mean that the projection of the gradient to that direction is non-null, which would in turn mean that the point was not really a minimum. The cost function may be expanded in a Taylor series around that minimum:
f ðpÞ ¼ f ðPÞ þX i ›f ðPÞ ›pi piþ 1 2 X i X j ›2f ðPÞ ›pi›pj pipj þ . . . ¼ f ðPÞ 2 bp þ1 2p †Hp; ð7Þ
where b is the negative gradient of f and H is the Hessian matrix of second partial derivatives of f at P. The gradient
is then simply expressed as
7f ðPÞ ¼ Hp 2 b: ð8Þ From Equation (8), we may immediately deduce the change in the gradient along the chosen direction:
dð7f ðPÞÞ ¼dðHp 2 bÞ ¼ HðdpÞ: ð9Þ To adopt the best possible direction from the found minimum P, the method must seek an orthogonal direction to the former direction vector, which obviously has to point along the direction of the gradient. If the former direction is denoted with n and the new direction is denoted with m, then
ndð7f ðPÞÞ ¼ nHm ¼ 0: ð10Þ When Equation (10) holds for vectors n and m, they are said to be conjugate. For as long as the minimisation is done only in conjugate directions, a single minimisation along a given direction is necessary, which implies quadratic convergence. Note that there is no need to compute the gradients at any point, only the orthogonal condition is used. This enables PDSM to preserve its derivative-less nature.
Our implementation closely follows the algorithm described in Lu et al. (1995). For Powell’s method, we consider a two-dimensional minimisation function f ðx; yÞ. In principle, we can always minimise the function in the directions of x and y successively. However, this could take many steps to reach the global minimum point. Obviously, what we need is a better set of directions. Powell’s method uses the averaged direction as a new direction to search for the minimum in each iteration. After each iteration, in which all directions in the set are optimised in turn, the overall distance moved in parameter space in that iteration is taken as a new direction. Owing to differences in image resolution and parameters used in object function, different parameters are considered and optimised over influence optimisation performance and registration robustness. As optimisation in image registration is to maximise similarity, similarity metric values, as functions of transformation parameters, comprise the objective function, henceforth denoted as f(x,y). Powell’s optimis-ation problem is formulated as minimisoptimis-ation problem and, thus, without the loss of generality, it is understood that for image registration, the goal is to minimise 2f ðx; yÞ.
4. Validation methodology
Image preparations of simulation video sequences are performed using a static image of capillaries and a nude mouse subject. Three different condition studies are simulated: video sequences are acquired at many different levels of motion blur, noisy translation and reflection
illumination effect. The purpose of this set of experiments is to evaluate the accuracy of the present alignment algorithms, and to assess its quantitative effects in actuality and speed in operation.
4.1 Blood flow change at vascular tree
In order to extract the precise area of vessels, the first step of processing is to increase virtual visibility. A mathematical morphology operation is performed. The virtual image is filtered using a greyscale opening with a flat structuring element (SE) of diameter larger than the maximum width of the retinal vessels. The opened image is then subtracted from the original image according to the following equation:
Ienchanced¼ Ioriginal2 ðIoriginal†SÞ
¼ Ioriginal2 ðIoriginalQðSÞ%SÞ; ð11Þ
where S is the flat SE, Q is greyscale erosion operation and % is greyscale dilation operation. Greyscale erosion operation is the process of placing the SE under the grey level of image I. In the resultant enhanced image (Figure 6(b)), it is noted that the morphological operation enhances vessel visibility in accordance with a diameter greater than the largest width of the virtual vessels.
In the second step to yield description of vessel-like patterns, a simplistic matched filter via Gaussian smooth-ing and vascular detection is convoluted based on the theory described in the literature (Matsopoulos et al. 1999). The result of vascular tree segmentation using Gaussian matching filter is shown in Figure 6(c). Moreover, blood flow change with time in vessels is different between similar frames. It is assumed that one RBC is 8 – 11 mm in diameter and it should fill a circle 8 – 11 pixels in diameter
(a)
(d) (e)
(b) (c)
Figure 6. Segmentation of vascular tree original greyscale image (a), after morphological enhancement (b) vascular tree segmentation using Gaussian matching filter(c), the simulation result of random blood flow at vascular tree (d). (e) An enlargement of the rectangular region in (d), with the original image on left side and simulated blood flow on right.
and can be randomly located anywhere in the vessel. The random cells are placed by Gaussian distribution function using uniform random numbers, in which it is possible to use standard deviation for dominating volumes of various density vessels. The intensity of one RBC ranges randomly from 0 to 10 levels of greyscale, which closely approximates the actual images as shown in Figure 6(d) and (e).
4.2 Displacement for simulation and measurement of alignment errors
Within this work, the simulated displacement of virtual blood flow images was first to verify our registration method. The mechanical simulations gave random positions from 0 to 50 pixels as a displacement of time. From these positions at each time step, a 2D magnified image sequence was simulated using microscopic image formation equations that consider motion blur, noisy translation and reflectance illumination effects. After simulating with the eccentric condition in actuality, the distance between registered image and reference image was determined. A least squares approach was applied to construct a single measurement of Euclidean distance using the root mean square error (RMSE) function. The distance was excluded from the analysis. The mean error, standard deviation of error and maximum error were calculated as accuracy estimation overall. The correlation coefficient (CC) is also calculated to serve as the robustness estimation for the trials registration. All displacements for simulation are implemented from 50 randomly generated transformations as described in this section.
4.3 Motion blur
Many types of motion blur can be distinguished, all of which are caused by relative motion between the camera and the object. This can be in the form of a translation, a rotation, a sudden change of scale or some combination of these. In this study, only the important type, translation, was included. When the object translates at a constant horizontal velocity during the exposure interval, the discrete equivalent point-spread function makes use of the blurring distance L, which is the number of additional points in the image resulting from a single point the in original scene. The blur image caused by regular heartbeats can be simulated at horizontal and vertical velocities given by hði; j; k; lÞ ¼ hði 2 kÞ ¼ 1 Lþ1; o # i 2 k # L; j ¼ l 0 otherwise ( ð12Þ
4.4 White noise in image
Real images are often degraded by some random noise. Noise can occur during capture or transmission. White noise is frequently applied as the approximation of worst degradation. A special case of white noise is Gaussian noise. Gaussian noise in image is fabricated by the density function, which is a very good approximation to realistic noise. Test-moving images were generated by adding different levels of Gaussian noise to the original images, and transforming the noise-corrupted images according to random transformations. The adopted Gaussian noise has zero mean with standard deviationl.
4.5 Surface reflectance and illumination effect When light is incident on a boundary interface between two different media with a microscopically smooth and flat material surface, such as float glass, the incident and reflected light rays make the same angle with the normal to the reflecting surface, producing specular reflection. If the material surface is ‘rough’, not microscopically smooth, diffuse reflections will occur. The reflection of light by a surface can be studied using optics. Two approaches were applied in the study of reflection.
Considering reflectance illuminated by a parallel light source, we assume that the shape of the threads can be approximated by a Gaussian function as a distribution function of the radiance around the direction of ideal reflection, LðxÞ ¼ k exp 2x 2 2s2 : ð13Þ
The standard deviation is used as a parameter to describe the degree of the surface, because the specular reflection and diffuse reflection are defined between low and high degree, respectively. The parameter k which denotes the intensity of illumination amplitude is created using uniform random numbers from 1 to 3 levels.
Figure 7 demonstrates the influences of various conditions as described above in the simulations.
5. Results and discussion
There are a lot of aspects to be evaluated with the present method of matching, including the robustness, the accuracy and the speed using simulated and real data. To investigate the characteristic of the present alignment algorithms, validation experiments were performed using simulation images. All of the present methods are run in Cþ þ on a personal computer (Intel Core CPU, 2.33 GHZ, RAM 4G). The quantitative analysis results are summar-ised in Tables 1 – 4. The variation of registration errors is small for different controlled values, indicating that all
Table 1. Estimations by motion blur effect at vertical velocity (n ¼ 50). RMSE and CC values are listed.
Matching metric SAD VR NCC MI GFM
L ¼ 10 r 1.000 1.000 1.000 1.000 1.000 RMSE (CC) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) L ¼ 20 r 1.00 1.00 1.00 1.00 0.999 RMSE (CC) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.02 ^ 0.14 (1.00) L ¼ 30 r 0.999 0.998 1.000 1.000 0.971 RMSE (CC) 0.04 ^ 0.20 (1.00) 0.48 ^ 0.50 (1.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 2.82 ^ 1.04 (5.00) Total r 0.999 0.999 1.000 1.000 0.999 RMSE (CC) 0.01 ^ 0.07 (1.00) 0.16 ^ 0.17 (1.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.95 ^ 0.40 (5.00) (a) (b) (c) (d)
Figure 7. Examples of simulation effects. (a) Displacement. Top panel is original image and bottom panel is original image with displacement. (b) Motion blur effect. Top panel is motion blur effect with L ¼ 20 horizontal velocity and bottom panel is motion blur effect with L ¼ 20 vertical velocity. Original images with white noise as shown in (c), top panel is SNR ¼ 29, bottom panel is SNR ¼ 21. Illumination effect as shown in (e), k ¼ 1, s ¼ 150 in top panel, k ¼ 2, s ¼ 150 in bottom panel.
Table 2. Estimations by motion blur effect at horizontal velocity (n ¼ 50). RMSE and CC values are listed.
Matching metric SAD VR NCC MI GFM
L ¼ 10 r 1.000 1.000 1.000 1.000 1.000 RMSE (CC) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) L ¼ 20 r 1.000 1.000 1.000 1.000 0.998 RMSE (CC) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 1.48 ^ 0.90 (3.16) L ¼ 30 r 0.991 0.998 1.000 1.000 0.960 RMSE (CC) 1.39 ^ 0.20 (1.41) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 3.16 ^ 0.54 (4.12) Total r 0.997 0.999 1.000 1.000 0.986 RMSE (CC) 0.46 ^ 0.07 (1.41) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 1.55 ^ 0.48 (4.12)
registration methods are accurate. For simulated images, the MI-based matching method is used in image alignment because of its insensitiveness to illumination changes and allows for fully automated alignment. Although segmen-tation, feature extraction or any other type of preproces-sing are not required for the calculation of similarity measures, MI has the shortcoming of complex compu-tation. Feature-based matching techniques do not use the grey values to describe matching entities, but image features are derived by a feature extraction algorithm. The approach assumes that point correspondence is available between the images. The binary images of vessel structures are obtained using a Gaussian filter for the registration process which maximises the computed
similarity measure. Feature-based matching computation is faster than MI, but it requires sophisticated image processing for feature extraction and its matching depends on the robustness of feature detection. As a result, the image matching precision is not as good as that of the intensity-based matching with motion blur effect. Furthermore, the matching stability is greater with lower complexity in operation and shorter matching time between the reference image and the target image for the methods such as the SAD, VR and NCC. This study finds that NCC is the most effective metric with average errors less than 1 pixel, and the CC greater than 1 for most cases, with three simulation effects: motion blur, noise and illumination effect. Table 5 lists the performance measures Table 3. Estimations by noise effect (n ¼ 50). RMSE and CC values are listed.
Matching metric SAD VR NCC MI GFM
SNR ¼ 62 r 0.999 1.000 1.000 1.000 1.000 RMSE (CC) 0.02 ^ 0.14 (1.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) SNR ¼ 29 r 0.999 0.999 1.000 0.998 0.999 RMSE (CC) 0.02 ^ 0.14 (1.00) 0.04 ^ 0.20 (1.00) 0.00 ^ 0.00 (0.00) 0.38 ^ 0.51 (1.41) 0.04 ^ 0.20 (1.00) SNR ¼ 21 r 0.999 0.999 0.999 0.998 0.999 RMSE (CC) 0.28 ^ 0.49 (1.41) 0.25 ^ 0.45 (1.41) 0.22 ^ 0.44 (1.41) 0.87 ^ 0.65 (3.16) 0.06 ^ 0.24 (1.00) Total r 0.999 0.999 0.999 0.998 0.999 RMSE (CC) 0.11 ^ 0.26 (1.41) 0.10 ^ 0.22 (1.41) 0.10 ^ 0.22 (1.41) 0.42 ^ 0.39 (3.16) 0.03 ^ 0.15 (1.00)
Table 4. Estimations by illumination effect on surface reflectance (n ¼ 50). RMSE and CC values are listed.
Matching metric SAD VR NCC MI GFM
k ¼ 1, s ¼ 100 r 0.999 1.000 1.000 1.000 1.000 RMSE (CC) 0.11 ^ 0.39 (1.41) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) k ¼ 1, s ¼ 150 r 0.998 1.000 1.000 1.000 1.000 RMSE (CC) 0.14 ^ 0.43 (1.41) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) k ¼ 2, s ¼ 100 r 0.666 0.998 1.000 1.000 1.000 RMSE (CC) 23.87 ^ 41.95 (123.37) 0.30 ^ 0.46 (1.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) 0.00 ^ 0.00 (0.00) k ¼ 2, s ¼ 150 r 0.415 0.528 1.000 0.097 1.000 RMSE (CC) 49.99 ^ 52.21 (126.81) 18.75 ^ 41.55 (126.09) 0.00 ^ 0.00 (0.00) 22.77 ^ 41.96 (110.92) 0.00 ^ 0.00 (0.00) Total r 0.769 0.882 1.000 0.774 1.000 RMSE (CC) 18.53 ^ 23.74 (126.81) 4.76 ^ 10.50 (126.09) 0.00 ^ 0.00 (0.00) 5.69 ^ 10.49 (110.92) 0.00 ^ 0.00 (0.00)
of time consumed described previously. For comparison of time consumed between our proposed measurement methods, we find that Powell’s search algorithm was able to effectively reduce the calculation time. Especially for the more complex calculation method MI, the time reduction ratio 71.94% is much higher than that of other proposed methods. Based on the results of simulation, we think that the VR method is simple, fast and robust to low levels of all the effects. Consequently, the NCC method is suggested to replace the VR method in practical application of microcirculation registration. After this simulation, the study used a set of images which were acquired using the microscopic system. All images were good in quality and did not require a pretreatment. Compared with the NCC method in Table 6, very good performance was reached through the use of the VR method with the same matching accuracy of the NCC method. Nevertheless, VR outperformed SAD with the illumination effect. Therefore, on the basis of the comparison, a general conclusion is reasonable.
6. Conclusion
This paper describes a microscopic system without fluorescent labelling to provide continuous quantitative blood flow data in tumour small vessels. Thus far, few studies in blood flow imaging have performed experimen-tation in mesenteric microvessels. We also propose a new applicable framework for automatic and rapid image registration to register unstable images of blood flow for large digital video sequences in the microscopic system analysis. For image sequences stabilisation, in which low resolution images, by down-sampling method, are used first as the pre-registration, the intensity-based alignment is adopted as a fine registration by linear translation between serial consecutive frames. The performance properties of matching metric are evaluated through image registration using simulated images. In this paper, an automatic medical image registration method based on Powell’s optimisation search method is implemented. This method is translation invariant, and has low calculation redundancy. By a systematic survey and rigorous evaluation of different methods, we conclude that the matching method of VR is computationally efficient and improves the registration robustness and accuracy in practical application of microcirculation registration. The presented registration method showed acceptable results in close requisition to analyse RBC velocities, confirming the scientific potential of the system and indicating that this automatic registration method can also be useful for many biotechnological analyses and evaluations.
Acknowledgements
This research was supported by a grant provided by the National Research Program for Genomic Medicine, Taiwan (Molecular and Genetic Imaging Core, NSC98-3112-B-010-012) and by China Medical University research grant (CMU98-N2-21).
References
Avants BB, Epstein CL, Grossman M, Gee JC. 2008. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegen-erative brain. Med Image Anal. 12:26 – 41.
Dobbe JGG, Streekstra GJ, Atasever B, van Zijderveld R, Ince C. 2008. Measurement of functional microcirculatory geometry and velocity distributions using automated image analysis. Med Biol Eng Comput. 46:659 – 670.
Hori K, Suzuki M, Tanda S, Saito S. 1991. Characterization of heterogeneous distribution of tumor blood flow in the rat. Cancer Sci. 82:109 – 117.
Table 6. Estimations of registration result for practical methods. Microscopic image (n ¼ 300) from real video cases. RMSE and CC values are listed.
Matching metric SAD VR
Case 1 r 0.999 0.999 RMSE (CC) 0.03 ^ 0.18(1.00) 0.04 ^ 0.18(1.00) Case 2 r 0.997 0.998 RMSE (CC) 0.06 ^ 0.26(1.41) 0.13 ^ 0.34(1.41) Case 3 r 0.991 0.996 RMSE (CC) 1.40 ^ 3.96(34.66) 0.31 ^ 0.51(2.24) Case 4 r 0.981 0.995 RMSE (CC) 9.52 ^ 7.74(26.08) 0.35 ^ 0.50(2.24) Case 5 r 0.979 0.993 RMSE (CC) 52.04 ^ 25.18(114.11) 0.45 ^ 0.58(3.00) Total r 0.989 0.996 RMSE (CC) 12.61 ^ 7.46(114.11) 0.25 ^ 0.42(3.00) Table 5. Time consumed (n ¼ 150 s).
Matching metric SAD VR NCC MI GFM
Direct search 103.26 149.85 246.42 4611.57 381.63
Powell’s search 64.59 83.59 134.68 339.92 179.48
Hsiao C-T, Chahine G, Gumerov N. 2001. Application of a hybrid genetic/Powell algorithm and a boundary element method to electrical impedance tomography. J Comput Phys. 173(2):433 – 454.
Huang H-L, Lin H-D, Chung B-T, Lin K-P. 2003. Structural and functional image registration using maximum mutual information. J Med Biol Eng. 23:111 – 118.
Iga AM, Sarkar S, Sales KM, Winslet MC, Seifalian AM. 2006. Quantitating therapeutic disruption of tumor blood flow with intravital video microscopy. Cancer Res. 66(24), 15 December.
Jain RK, Ward-Hartley K. 1984. Tumor blood flow-character-ization, modifications, and role in hyperthermia. IEEE Trans Sonics Ultrasonics. 31(5):504 – 525.
Jenkinson M, Smith S. 2001. A global optimisation method for robust affine registration of brain images. Med Image Anal. 5:143 – 156.
Jeong JH, Sugii Y, Minamiyama M, Okamoto K. 2006. Measurement of RBC deformation and velocity in capillaries in vivo. Microvasc Res. 71:212 – 217.
Lu L, Brown BH, Barber DC, Leathard AD. 1995. A fast parametric modelling algorithm with the Powell method. Physiol Meas. 16(3):A39 – 47.
Maes F, Vandermeulen D, Suetens P. 1999. Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information. Med Image Anal. 3(4):373 – 386.
Matsopoulos GK, Mouravliansky NA, Delibasis KK, Nikita KS. 1999. Automatic retinal image registration scheme using global optimization techniques. IEEE Trans Inf Technol Biomed. 3(1):47 – 60.
Matsopoulos GK, Delibasis KK, Mouravliansky NA, Asvestas PA, Nikita KS, Kouloulias VE, Uzunoglu NK. 2003. CT-MRI automatic surface-based registration schemes combin-ing global and local optimization techniques. Technol Health Care. 11:219 – 232.
Mchedlishvili G. 1998. Disturbed blood flow structuring as critical factor of hemorheological disorders in microcircula-tion. Clin Hemorheol Microcirc. 19(4):315 – 325.
Plattard D, Soret M, Troccaz J, Vassal P, Giraud JY, Champleboux G, Artignan X, Bolla M. 2000. Patient set-up using portal images: 2D/2D image registration using mutual information. Comput Aided Surg. 5:246 – 262.
Sugii Y, Nishio S, Okamoto K. 2002. In vivo PIV measurement of red blood cell velocity field in microvessels considering mesentery motion. Physiol Meas. 23(2):403 – 416.
Wang Z, Lu L, Bovik AC. 2004. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 13(4):600 – 612.
Wei Z, Tan T, Sun Z, Cui J. 2005. Robust and fast assessment of iris image quality. Adv Biomet :464 – 471.
Winkler S. 1999. Issues in vision modeling for perceptual video quality assessment. Signal Processing. 78:231 – 252. Woods RP, Cherry SR, Mazziotta JC. 1992. Rapid automated
algorithm for aligning and reslicing PET images. J Comput Assist Tomogr. 16(4):620 – 633.
Wu C-C, Lin K-P, Chung B-T. 2007. A study of image quality analysis for a cutaneous capillary imaging system. J Med Biol Eng. 27(3):110 – 115.
Xu X, Dony RD. 2004. Differential evolution with Powell’s direction set method in medical image registration. IEEE Int Symp Biomed Imaging. 1:732 – 735.
