Moment algorithms for blood vessel detection in infrared images of laser-heated skin

Abstract Block-detection and three-dimensional moment algorithms are applied to determine the presence and orientation of cylinders in a three-dimensional image. The proposed method detects blood vessels in a three-dimensional tomographic image constructed by infrared radiometry. Based on the moment principle, the method can be employed to determine the centroid, diameter, and orientation of arbitrarily shaped blood vessels for application in the laser treatment of hypervascular skin lesions.


INTRODUCTION
Successful laser treatment of hypervascular lesions is based on proper selection of the irradiation parameters [l-3]. Kimel et al. [4], and Nelson et al. [5] have shown laser pulse duration should approximately equal the thermal relaxation time, z, of the targeted blood vessels.
Three-dimensional reconstruction of laser heated discrete subsurface chromophores in human skin recorded by infrared imaging radiometry is being investigated to determine the initial space dependent temperature distribution in skin to pulsed laser exposure [6, 71. By imaging the emitted radiation onto an infrared focal plane array, useful information regarding the laser irradiated skin can be derived [S]. The value of z, is directly proportional to the squared diameter (d) of the target blood vessel and inversely proportional to thermal dzflusiuity of skin, z, = d2/16X, where X is thermal diffusivity and governs the rate of diffusive heat propagation in the material (1.1 x 10-l cm'/s) [9].
We model a blood vessel as a cylinder and use the block-detection and three-dimensional moment algorithms to estimate centroid, diameter, and orientation. This information is important for proper selection of laser pulse duration. Results presented here combined with access to lasers with user-specified pulse durations (0.25-15 ms) may provide a methodology for improved PWS (port wine stains) treatment.
In Section 2, moments for blood vessel detection and their properties are described. Algorithms for block detection and three-dimensional moment calculation are given in Section 3. Experimental results and conclusions are given in Sections 4 and 5, respectively.

MOMENTS
To detect a pattern in an image, an invariant feature is selected independent of size, position, and orientation [lO-161. The mathematical foundation of an invariant feature is based on calculus of algebraic invariants [17]. Hu [l l]

Three-dimensional moments
Given an object with a three-dimensional density function, f(x,y,z), the (p + q + r)th order moments are defined in terms of a Riemann integral as where r, is the normal distance to axis i. i, i = x,y,z, and p,q,r = 0,1,2,... The integration extends over domain off.
For an object with a finite extent (i.e. bounded domain), the integration extends over the object volume. The second order moments (e.g. (p = 2, q = 0, r = 0) moment of inertia) about x, y, and z axes are where rt = y2 + z2, rf = x2 + z2, and r! = x' + y'. Definition 1. An object centroid is defined by coordinates X, 7, and 2 given by2 = m,oo/m~, j = m,,,/m,, and Z = moo,/m~, where mOoO is the object volume. Central object moments with density functionf (x,y,z) are defined as ix +r sss +* Ppqr = (x -X)Q -y)4(z -z)lf(x,y,z)dxdydz. --r -r _I A uniqueness theorem states that iff(x,y,z) is piecewise continuous and bounded in a finite region in R3 space, then moments of all order exist implying that the application of a small set of low-order moments may be used to distinguish different patterns.
We consider binary images so that f(x,y,z) is either 0 or 1. In a three-dimensional image, x, y, and z correspond to row (first dimension), column (second dimension), and frame (third dimension) of the voxel, respectively. One representation of a digital image is a collection of voxels with associated intensity values. For industrial applications, region segmentation or edge detection are used to transform images into a binary representation, where each voxel is either 0 (white) or 1 (black). Equation (1) is written as a summation, and the (p + q + r)th order moment of a three-dimensional image with n rows, m columns, and I frames is (2) Here xyk, y,jk, and z# are the black pixel coordinates, and rn.n.1 is the total number of voxels.
Similarly, central moments of a three-dimensional binary image are

CYLINDER DETECTION ALGORITHMS
In this section, algorithms for detecting the presence, diameter, and orientation of a cylindrical object in a three-dimensional image are presented.

Block-detection algorithms
An algorithm for detecting connected black pixels in an n x m binary image is described. Connectivity among pixels can be defined in terms of their adjacency. Figure 1 illustrates a pixel with eight neighbors.
Two black pixels (x,,J+) and (xz,yz) are an g-neighbor if max(lx, -%I, IYI -Yzl) 5 1  In this paper we use g-neighbors. Since a three-dimensional image has blocks in different frames of the image, a block-detection algorithm identifies three-dimensional blocks. Since mloo = mOIO = 0, the centroid is at (O,O,c/Z). If a given block image contains n x m x 1 voxels, the (p + q + r)th order moment are calculated (Equation

Algorithms for calculating cylinder diameter and orientation
In this section, algorithms to calculate cylinder orientation and diameter in a three-dimensional image are presented. An important cylinder property is described in Lemma 1.
Lemma 1. In a cylinder with x2/a' + y2/u2 = 1, 0 < z I c and a < c, the second moment is minimum along the major axis and maximum along the minor axis; these moments are given by   With the assumption that a < c, then moo2 < mzoo.
For a cylinder m,,a = ml,,, = mo,, = 0. Definition 2. The angles between minor, intermediate, and major cylinder axes and positive direction of the x, y, and z axes, CC, /3, and y, are defined as the cylinder orientation with respect to a given coordinate system (Fig. 2).

R=[;:]=[E;]. r3 cos y
The cylinder orientation in terms of different order moments can be calculated. The direction R can be found by determining the eigenuector corresponding to the minimum eigenvalue which represents the central moment along the major axis of the cylinder.  Fig. 7. Blocks detected by the block-detection algorithm in frame 3.
where I is a diagonal unit matrix. By assuming that a coordinate system originating at the centroid, where the z-axis is coincident with the major axis of the cylinder, the roots of the characteristic equation are In this configuration, moments and central moments are equal, and roots of the characteristic equation are the second moments as given in Lemma 1 (Equation (4)).
for blood vessel detection in infrared images 365 Solving for a and c a= J 3 V where 2a is the diameter and c is the length of the cylinder. The eigenvector R corresponding to the minimum eigenvalue (A,) gives the orientation of the cylindrical object. To determine the

EXPERIMENTAL RESULTS ON BLOOD VESSEL DETECTION
Infrared imaging radiometry uses a high speed infrared focal plane array (IR-FPA) camera to measure temperature changes induced in human skin exposed to pulsed laser radiation. In practice, a pulsed laser is used to produce transient heating of the object under study (Fig. 3). Heat generated due to light absorption by subsurface blood vessels in skin diffuses to the surface and results in increased infrared emission levels which is measured by a fast IR-FPA. If a pulsed laser source is used to irradiate the skin, an immediate increase in infrared emission will occur due to optical absorption by hemoglobin contained within the blood vessels. The infrared signals collected by each detector element are digitized by a 5 MHz, 12-bit A/D converter giving a frame rate of 217 Hz. An infrared imaging radiometry record of skin in response to pulsed laser exposure is composed of sequence of 150 gray-level frames (128 x 128) which are stored and accessed using a computer.

Processing program
Commercial software (Application Visualization System (AVS), Waltham, MA) is used to process input images. Figure 4 illustrates frame 1 of the input image. Dark pixels represent the blood vessels.
Each gray-level image undergoes a binary transformation using a threshold set at 85% of the average pixel value for each frame. The threshold was calculated using the background level. Figure 5 shows frame 1 after the gray-scale to binary transformation where numeral one represents a black voxel. 4.1.1. Block-detection algorithm. Simulation programs for detecting blocks and calculating moments are written in C programming language. Figs 6 -8 represent blocks detected by the block-detection algorithm in frames 1, 3 and 5, respectively. 4.1.2. Blood vessel centroid. A time sequence of infrared emission frames may be viewed as a three-dimensional image, and moments can be applied to calculate the blood vessel centroid detected by the block detection algorithm. Table 1 lists the centroid coordinates of the first 15 blocks calculated using Equation (3). 4.1.3. Diameter and orientation of blood vessels. In the analysis of input image, roots of the characteristic equation (Equation (5)) were determined and used to calculate the diameter of blood vessels [Eq.(6)]. Table 2 illustrates the diameter (in pixels) of the first 15 detected blood vessels. Orientation for the same vessels is presented in Table 3.

CONCLUSION
Moments of a three-dimensional tomographic image constructed by infrared radiometry contain information which can be used to calculate the centroid, diameter, and orientation of individual blood vessels. Blood vessels, modeled as cylinders, are identified by a block-detection algorithm. Three-dimensional moment algorithms are used to compute diameter and orientation of blood vessels. This information may be important for proper selection of laser pulse duration for improved PWS treatment. Although the algorithms used images derived by infrared radiometry, this procedure can be implemented on images generated by other modalities (e.g. MRI). Results demonstrate the utility of the moment concept in medical image processing.