# Possibility of non-invasive blood pressure estimation by measurements of force and arteries diameter

Florent VEYE<fnref n="1" />, Sandrine MESTRE<fnref n="2" />, Antonia PEREZ-MARTIN<fnref n="3" />,<fnref n="4" />, and Jean TRIBOULET<fnref n="1" />,<fnref n="5" />

**Abstract**— Ultrasound examination is the first line procedure for the diagnosis and follow-up of cardiovascular diseases. Instrumenting an ultrasound probe with a force sensor may improve thenon-invasive measurement of arterial biomechanical parameters (diameter, pulsatility, intima-media thickness and flow-dependent dilation) by measuring and controlling the force exerted by the sonographer. We present here the results obtained with this approach coupled with image processing for the measurement of arterial diameter changes during the cardiac cycle and dynamic tests. The analysis of artery force feedback provided non-invasively an estimation of blood pressure waveform.

## I. INTRODUCTION

The non-invasive diagnosis and follow-up of cardio-
vascular diseases need reliable measurements of arterial
function. The evaluation of the arterial wall elastic modulus
can be based on distensibility measurement as the ratio of
arterial diameter and blood pressure changes during the
cardiac cycle [1-8]. This assessment can be improved by the
simultaneous measurement of the arterial wall thickness and
recording of the blood pressure curve at the same site [9,10].
Arterial diameter and wall thickness measurements require
high resolution ultrasound scanners and dedicated image
processing software. B-mode echographic images analysis is
generally performed with non-isotropic filtering to smooth
speckle noise, and level set methods for contour extraction
[11]. Nevertheless, these methods proved not robust enough
for this application. Therefore, we used simple but robust
operators to extract the arterial wall contour and estimate
vessel diameter and diameter changes.

The effect of the pressure exerted by the operator on the
ultrasound probe is generally not taken into account during
arterial measurements although it is liable to affect the
amplitude and the waveform of the diameter changes upon
which distensibility measurements are based [12]. The
validity of non-invasive arterial pressure measurement
(sphygmomanometry, oscillometry...) remains controverted
when applied to arteries others than the brachial artery, but is
widely used in this context.

Different, and often complex, approaches, have been
proposed to estimate the local arterial pressure from the
conventional measurement performed on the brachial artery
[13]. Most of these approaches are based on a priori
estimation of the modifications that undergo the pressure
wave during its propagation along arteries, while vascular
diseases, by their consequences on the biomechanics of
arterial wall, are likely to affect this propagation. Obtaining
the arterial blood pressure waveform at the exact site of
diameter measurement would drastically improve the
accuracy and the reliability of the arterial distensibility
measurement.

Our aim was to couple the measurement of the force
exerted on the ultrasound probe with image analysis for the
evaluation of arterial distensibility in order to obtain an
image of the local blood pressure waveform. In the first part
of this article, we present the corresponding hardware and
software developments. In the second part, we report the
clinical results we obtained; showing a correlation between
the arterial diameter changes and the force feedback curve.
To conclude, we discuss the potential applications of this
approach.

## II. MEASUREMENT SYSTEM

### A. Hardware and software developments

The acquisition system was composed mainly of a
computer equipped by a Matrox® analogic video acquisition
card and by a National Instruments® multifunction data
acquisition board. The video acquisition card was connected
to the analogic video output of the ultrasound scanner *via* a
coaxial cable. A force sensor (XFTC300, FPG Sensors, Les-
Clayes-Sous-Bois, France) was connected to a signal
amplifier/conditioner (ARD154, FPG Sensors, Les-Clayes-
Sous-Bois, France) whose output was connected to the data
acquisition board. The linear force sensor was fixed on a
sliding support receiving the linear probe of the ultrasound
scanner (Fig.1). Clinical results were obtained with a Logiq 9
duplex Doppler ultrasound system (GE Healthcare
Ultrasound, Milwaukee, USA) with a 10 MHz linear probe.

The workflow included video and force signal acquisition
and synchronized B-mode ultrasonographic image
processing for the measurement of arterial diameter and
diameter changes.

<footnote n="1">¹ LIRMM – UMR 5506 – UM2 – CNRS, 161 rue Ada 34095 Montpellier, France.</footnote>

<footnote n="2">² Internal and Vascular Medicine department – Montpellier University Hospital Center, 191 avenue du Doyen Gaston Giraud 34295 Montpellier, France.</footnote>

<footnote n="3">³ Vascular Medicine and Laboratory department – Nîmes University Hospital Center, Place du Pr R-Debré 30029 Nîmes, France.</footnote>

<footnote n="4">⁴ Dysfunctions of Vascular Interfaces Laboratory – EA2992 – UM1, Montpellier and Nîmes, France.</footnote>

<footnote n="5">⁵ Nîmes University, Place Georges Salan 30921 Nîmes, France.</footnote>
<figure><img src="image_58.png" /><figcaption>Fig. 1 Drawing showing the linear sensor implementation on the linear ultrasound probe.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>Fig. 2 Hardware architecture.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>Fig. 3 The successive stages of image processing.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Fig. 4 Image processing with edges detection and smoothing by polynomial approximation.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Fig. 5 Brachial artery diameter changes calculated on a sequence of images without and with polynomial approximation, in a normal subject.</figcaption></figure>

## B. Method

Data acquisition was controlled by a laboratory-written piece of software synchronising force and video data acquisition and processing. One force value was obtained and saved for each echographic image. Another piece of software processed the acquired video sequence to extract the arterial diameter as follows:

First, image calibration was performed to allow pixel to millimetres conversion. Thereafter, the operator had to select the region of interest (ROI) in which the following analyses were performed (Fig.3):

* Closing operation: erosion followed by expansion allowing to remove artefacts and to smooth out the vessel's contour using a 9x9 structuring element on neighbourhood.

* Adaptive threshold: The suitable threshold level was calculated on the grey scale along a line parallel to the ultrasound beam,

* A Canny's filter [14], with a 3x3 window width, was applied.

* Contour superimposition was performed on each image.

At the centre of the ROI, a 20 pixels width area was automatically drawn, and the results of the previous stages of image processing were further processed for smoothing by 2^nd^-order polynomial approximation in order to be a part of the best-fit parabola following the arterial wall. This process enabled to suppress residual artefacts (Fig.4).

This 20 pixels width allowed eliminating the effect of isolated acoustic artefacts, thus overcoming the pitfall of other techniques like echo-tracking while avoiding the excessive temporal smoothing that would happen with a larger area because of the propagation of the arterial pulse wave (e.g. on the brachial artery of a normal subject, the pulse wave propagates along the wall with a velocity about 7 m.s⁻¹).

The distance separating the two interpolated polynomials (i.e. proximal and far arterial wall) was then averaged to get the arterial diameter, and this calculation was repeated along the sequence of images (about 300 for a typical sequence). Fig.5 shows the variation of the arterial diameter on a full sequence with a mean value of 5.06 mm.

Frequency analysis of the resulting waveform yielded results a mode about of 1.2 Hz, i.e. the normal frequency of heart beats. Direct calculation and polynomial approximation gave quite similar results because most artefacts had been already corrected by the operation of mathematical morphology.
### III. CLINICAL TESTS

Clinical tests and measurements were performed in collaboration with vascular laboratory of the Nîmes University Hospital.:

* An Alpha 10 ProSound ultrasound scanner (Hitachi-Aloka, Tokyo, Japan) with a 10 MHz linear probe equipped with the above-mentioned force sensor.

* A Dinamap DPC320NR (GE Healthcare, Milwaukee, USA) electronic tensiometer.

* An articulated probe-holder mounted on a 3-axis adjustment device (MP-PH0001, Hitachi-Aloka, Tokyo, Japan).

* An ECG and breath-movement monitor with analog output.

When the sonographer handled the probe equipped with the force sensor without using the micrometric probe-holder for the examination of the brachial artery in normal volunteers, we were able to acquire and display simultaneously the force signal and the brachial artery diameter calculation results.

The results confirmed that arterial diameter measurements could be affected by the way the sonographer handled the probe, depending on the patient's size and posture, and that the force exerted on the probe could be quite different along the examination, and quite different between sonographers.

Therefore, the articulated probe-holder proved essential for accurate and reproducible measurements, allowing to maintain the ultrasound probe in the required position throughout the acquisition procedure, with a constant force. The probe-holder was used for all subsequent measurements.

#### A. Relationship between arterial diameter and force:

Fig. 6 shows the effects of increasing the force applied by on the ultrasound probe on the brachial artery diameter. Increasing the force F from 0 to 6.75 N resulted in a decrease in mean arterial diameter ($Φ$) from 4.6 to 3 mm in this subject.

Performing a polynomial approximation of both the force and the arterial diameter curves, and subtracting the result from initial data, allowed us to centre the waveforms (Fig. 7). We could then observe that with increasing of the force resulted in the diastole-to-systole diameter ($ΔΦ$) and force feedback ($ΔF$) changes during the cardiac cycle increased jointly by a factor 2, respectively from 0.15 to 0.3 mm and 0.175 to 0.350 N. These observations supported our hypothesis that the force applied on the ultrasound probe during the examination does influence the measurement of arterial diameter, and the force feedback represents an image of the local blood pressure waveform.

<figure><img src="image_63.png" /><figcaption>Fig. 6 Raw data showing a decrease in mean brachial artery diameter ($Φ$) when increasing the force (F) exerted on the ultrasound probe.</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Fig. 7 Centred data showing an increase in diastole-to-systole brachial artery diameter changes when increasing the force exerted on the ultrasound probe.</figcaption></figure>

#### B. Interest of exerting a constant force:

Using the probe-holder, we exerted a constant force on the ultrasound probe during a 20 s recording. This force was the minimal force allowing obtaining a stable and satisfying B-mode image of the brachial artery. The patient was first asked to breathe normally, then to maintain an apnoea. The first part of the graph showed 3 cycles of change in $ΔF$, whereas it remains stable during apnoea. The period of these cycles was roughly 4 s, and corresponded to the respiratory rate. The diastole-to-systole diameter ($ΔΦ$) and force feedback ($ΔF$) changes showed inverse variations during breathing (respectively 0.1mm and 0.072 N). This demonstrated that breathing affected diameter measurements probably by inducing slight arm movements.

<figure><img src="image_65.png" /><figcaption>Fig. 8 Raw data showing brachial artery diastole-to-systole diameter and force feedback changes during breathing as measured with a constant force exerted on the ultrasound probe.</figcaption></figure>
IV. CONCLUSION

The system we developed allowed simultaneous acquisition of ultrasound images and measurement of both the force exerted on the ultrasound probe by the sonographer and the force feedback resulting from arterial blood pressure pulsatility. We also developed image analysis software providing automatic and synchronized arterial diameter measurements. We were thus able to obtain measurements of the brachial artery diameter and diameter changes during the cardiac cycle, and showed the influence of the pressure exerted on the ultrasound probe by the operator on the arterial diameter and its cyclic variations. Our image analysis method has proved robust enough to correlate these force measurements with the arterial diameter and the amplitude of diameter changes during the cardiac cycle. We also demonstrated the ability of our system to capture the force feedbacks, providing an image of systemic blood pressure waveform. Moreover, we could show the influence of the breathing on arterial diameter changes.

Calibrating force feedback measurements would be
needed before this approach can be used to offer a non-
invasive method for providing the arterial blood pressure
waveform. This can be achieved by simultaneous systemic
blood pressure measurement by conventional method on the
brachial artery on the other side, but this approach would
have to be validated on a series of normal and pathologic
subjects.

REFERENCES

[1] J. M. Bland and D. G. Altman, "Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement," The Lancet, vol. 327, no. 8476, pp. 307-310, 1986.

[2] T. Gustavsson, R. Abu-Gharbieh, G. Hamarneh, and Q. Liang, "Implementation and Comparison of Four Different Boundary Detection Algorithms for Quantitative Ultrasonic Measurements of the Human Carotid Artery," in Cardiology 1997, vol. 24, no. December, pp. 69-72, 1997.

[3] R. H. Selzer, W. J. Mack, P. L. Lee, H. Kwong-Fu, and H. N. Hodis, "Improved common carotid elasticity and intimamedia thickness measurements from computer analysis of sequential ultrasound frames," Atherosclerosis, vol. 154, no. 1, pp. 185-193, 2001.

[4] F. Beux, S. Carmassi, M. V. Salvetti, L. Ghiadoni, Y. Huang, S. Taddei, and A. Salvetti, "Automatic evaluation of arterial diameter variation from vascular echographic images," Ultrasound in Medicine &amp; Biology, vol. 27, no. 12, pp. 1621-1629, 2001.

[5] B. Haluska, A. Sutherland, C. Case, R. Kennedy, and T. H. Marwick, "Automated edge-detection technique for measurement of brachial artery reactivity: a comparison of concordance with manual measurements," Ultrasound in Medicine & Biology, vol. 27, no. 9, pp. 1285-1289, 2001.

[6] M. Cinthio, T. Jansson, A. R. Ahlgren, K. Lindström, and H. W. Persson, "A method for arterial diameter change measurements using ultrasonic B-mode data.," *Ultrasound Med. Biol.*, vol. 36, no. 9, pp. 1504-12, Sep. 2010.

[7] T. Nilsson and S. Segstedt, "A robust and fast method for arterial lumen diameter and intima-media thickness measurements," IEEE Int. Ultrason. Symp. Proc., pp. 1678-1681, 2012.

[8] S. Graf, D. Craiem, and R. L. Armentano, "Non invasive assessment of carotid and femoral arterial pressure using B-mode ultrasound diameter waveforms.," *Conf. Proc. IEEE Eng. Med. Biol. Soc.*, vol. 2012, pp. 5610-3, Jan. 2012.

[9] W. Yousheng, S. Xiaodi, and C. Jianxin, "A Method of Analyzing the Strain of Arterial Wall," *2009 Int. Forum Inf. Technol. Appl.*, no. 1, pp. 217-220, May 2009.

[10] F. Sun, Z. Liu, Y. Li, and P. Babyn, "Improved T-snake model based edge detection of the coronary arterial walls in intravascular ultrasound images," ICBBE, no. 1, pp. 1-4, 2009.

[11] N. Santhiyakumari, P. Rajendran, M. Madheswaran, and S. Suresh, "Detection of the intima and media layer thickness of ultrasound common carotid artery image using efficient active contour segmentation technique.," Med. Biol. Eng. Comput., vol. 49, no. 11, pp. 1299-310, Nov. 2011.

[12] J. Triboulet, E. Nasr, and M. Dauzat, "Evaluation of the Influence of Probe Pressure on the B-Mode Ultrasound Measurement of Arterial Diameter," in EMBC'06: 28th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. New York (USA): IEEE, 2006, pp. 3831-3835.

[13] J. Su, J.-J. Wang, C.-C. Chen, and L.-C. Lin, "A Cuff-Based Approach for Assessing the Flow-Mediated Dilatation of Brachial Arteries," 5th Int. Conf. Bioinforma. Biomed. Eng., no. 2, pp. 1-4, May 2011.

[14] J. F. Canny, "A Computational Approach to Edge Detection." IEEE transactions on pattern analysis and machine intelligence, vol. 8, no. 6, pp. 679-98, Jun. 1986.
# Comparative study of two 3D reconstruction methods
for underwater archaeology

Arnaud Meline, Jean Triboulet and Bruno Jouvencel

**Abstract**— The underwater 3D reconstruction cartography has made great progress in the last decade. The work presented in this paper is about the analysis and 3D reconstruction of archeological objects. Using a calibrated single camera and an uncalibrated system, we propose to describe a method to perform the Euclidian 3D reconstruction of unknown objects. A comparison of two methods is presented and tested on synthetic and real underwater pictures. Filters are proposed to simulate underwater environment and inherent problems. Finally, robust and stable features have been extracted from underwater pictures and used to perform the 3D model.

## I. INTRODUCTION

In this paper, our aim is to analyze natural underwater scenes and especially the 3D cartography of submarine environments. To handle this problem, several techniques are widely used today. The originality of our project involves the fusion of two kinds of maps obtained with sensors of different resolutions. An Autonomous Underwater Vehicle analyzes the seabed with a lateral sonar to construct a first global map of the zone (fluvial archaeology or coastal oceanographic applications). This map is then decomposed into smaller cells representing a mosaic of the seabed. A second scanning is performed on particular cells using a second sensor with a higher resolution leading to a detailed 3D partial map. In our case, we aim at applying this method to submerged archaeological sites, where objects of interest (statues, plates, amphora) are detected in the global map, while 3D precise reconstruction would be necessary. This project includes two parts: the first part correspond to the 3D reconstruction of constrained submarine environment with a video camera. The Second part is the multimodal aspect of the problem. The work presented in this paper focuses on the first part of the project and deals with the 3D reconstruction of the submarine environment with a video camera.

Nowadays, even if vision can solve many problems, images analysis in case of submarine environments is quite complicated. The different methods used on the surface are not robust enough to changes produced by the underwater medium. So, the robustness of several methods decreases and their results become unstable.

When a scene or an object must be reconstructed in 3D, detection and matching points in the image are crucial parts for the accuracy of the model. If these steps contain errors, the 3D model will be of insufficient quality or completely wrong. Thus, we used a well known methodology

for 3D reconstruction and we have tested and compared the robustness of the algorithms which are sensitive to underwater disturbances. We propose to use surface images on which are applied filters simulating underwater characteristics. Therefore, we can recreate the underwater conditions and determine which setting defeats the reference methods used on the surface. From these results, we selected the best methods to achieve the 3D model of a submerged object.

This paper is organized as follows. Section II presents a survey on 3D reconstruction and problems caused by the underwater environment. In section III, a method to reconstruct underwater 3D scene is presented. Section IV describes the detector and the matching methods. In section V, the filters employed to simulate underwater characteristics are detailed. Section VI presents experimentations results: the comparison of detector methods used and the 3D reconstruction. Finally, conclusion and future work are presented in section VII.

## II. RELATED WORKS

In recent years, the community of computer vision widely studied the topic of 3D reconstruction from image sequences. Numerous methods exist, but they cannot be applied to all underwater images. However, they depend on the knowledge of the environment and the used system. For example, the work of Snavely [1], Bartoli [2] and Barazzetti [3] using unknown images to create 3D urban models. They have no information about the devices that took the pictures but their methods utilize a priori knowledge of the scene. By focusing on underwater reconstructions, we find the work of Brandou [4] who proposes to calibrate the camera in situ and to use a robot arm to move it. With this system, he knows the exact camera movement and simplifies the reconstruction. Espiau [5] works on the same images but proposes a method to reconstruct textured scenes from an uncalibrated camera. To achieve this problem, he proposes to use the scene rigidity to estimate the camera parameters. In a different way, Kim [6] proposes a SLAM approach to inspect the ship hull status. He determines the camera position in 3D using a calibrated monocular camera and a geometric model selection. Unlike Beall [7] or Negahdaripour [8] who uses a full calibrated stereo rig to pre-built the 3D trajectory or to inspect ship hull. Then, he simultaneously optimizes the camera trajectory and estimates a 3D dense reconstruction. In a same way, Hogue [9] combines a 3DOF inertial sensor and a calibrated stereo ring to estimate the trajectory and creates 3D dense map.

All these applications impose constraints (type of scene, knowledge of the camera movement, calibration parameters, shape of the object) which change the approach. Some

A. Meline, J.Triboulet and B. Jouvencel are with LIRMM, Univ. Montpellier 2 - CNRS, 161 rue Ada, 34392 Montpellier, France {meline, triboulet, jouvencel}@lirmm.fr

This work is supported thanks to ANR project CFLAM and UNIMES.
authors, like Hartley [10], Pratt [11], Szeliski [12] and Nicosevici [13] propose general methods to reconstruct a 3D scene. To obtain Euclidian structure of scenes, it is fundamental to know the intrinsic camera calibration parameters. Many methods exist to estimate these parameters on the surface but when the camera is in a submarine environment, it introduces optical problems. To avoid these problems, Pessel [14] proposed a self-calibration method. Other authors like Lavest and Rives [15] proposed to adjust the surface parameters to the underwater environment. We [16] proposed another approach that includes the environment change in the radial distortion model.

All the reconstruction methods need robust feature points to obtain the best 3D model. So the feature points detection and matching is a crucial step for 3D reconstruction. There are many detection points methods [12] such as Moravec corner detector, Kitchen and Rosenfeld detector, Forstner operator, Hessian detector, SURF, Ferns, Harris and Stephens detector, SIFT... All of these feature detectors have invariant characteristics in a spatial domain. But they have lacking quality results when images undergo important modifications like underwater images. Indeed, the presence of suspended particles, the energy absorbing of light ray and the variation of water refractive index make blur and add noise on images. All methods are not robust enough for such modifications. In next parts, we will compare a robust method (SIFT) with respect to a classical one (Harris and correlation).

### III. 3D RECONSTRUCTION METHOD

As the geometry of the scene and the camera motion are unknown, we used the method proposed by Hartley and Zisserman [10] to achieve 3D reconstruction. Fig. 1 resumes the different steps to obtain the 3D model. Each block will be detailed.

<figure><img src="image_66.png" /><figcaption>Fig. 1: 3D reconstruction diagram</figcaption></figure>

#### A. Camera calibration

The goal of our problem is to obtain an Euclidian reconstruction, so we have to calibrate our camera. To take into account the underwater parameters which affect the precision, we used the method presented in [16]. We proposed to include the underwater constraints in the camera distortion model. Thus, we used a pinhole model taking into account radial distortion terms. Then calibration is done using the Zhang's method [17]. The major advantage of this method is the utilization of a simple planar calibration target and the integration of underwater optical constraints.

#### B. Feature points detection and matching

As we saw in section II, the detection and matching of the feature points is crucial to obtain a representative 3D model of the observed scene. In this paper, Harris and SIFT method are used to obtain matching points and fundamental matrix. The motivation of those choices and detailed methods will be presented in section IV.

#### C. Euclidian Reconstruction

To perform the reconstruction, we use the essential matrix method proposed by [10]. The purpose of this method is to obtain the pose of the second camera with respect to the first one. We define the following projection matrix $P1 = K[I|0]$ for the first camera and $P2 = K[R|t]$ for the second. Then, we estimate the rotation matrix (R) and the translation vector (t). The essential matrix $E = K^tFK$ is calculated from fundamental (F) and calibration (K) matrices obtained previously. Then, calculating the Single Value Decomposition (SVD) of the essential matrix gives us the camera pose (R and t). Finally, the 3D reconstruction is done by a linear triangulation from the inliers found earlier and the projection matrices.

The 3D model is then treated to remove some aberrant points that are locally isolated. The mesh structure is achieved through a 3D Delaunay triangulation.

### IV. DETECTION AND MATCHING

This phase is divided in two parts. The first part shows how to extract feature points on images and the second part present how to realize the matching between all detected points and check if all matching points are correct using the epipolar geometry constraints. In this section, we present the two methods used to detect and match points, then we detail the verification and the removal of false matching (outliers).

#### A. Features points detection and matching

The feature points detection and matching is an essential step for 3D reconstruction. The quality of the 3D model is related to the precision of these elements.

There are many detection point methods but they do not have the same performances. Indeed, Harris can find points on objects, more specifically near the corner. This enables us to obtain the features points of the object's geometry.

SIFT is considered as one of the best performing detectors in many applications because of its robustness to scaling, rotations, translations and illumination changes. We decided to compare these two representative methods in our work.

1) Method 1: *Harris detector and correlation matching:* Harris and Stephens [18] proposed an operator to detect corners in image. They used a criterion based on corners and edges informations by calculating the image derivative along x and y. The values of this criterion is positive near a corner, negative near an edge and low in a region of constant intensity. The detected points are selected by choosing the highest values.

To match these points, we used a Normalized Cross Correlation (NCC) method. It consists in searching the corresponding points in the other image which have a maximum
of correlation value. To improve performance and reduce the computational time, we used a guiding approach for this matching. A method detecting the global motion between the two images is used to reduce the number of candidates to selected points. Then, only pairs of corresponding points in both directions are matched and one point can only have a single matched point in the other image. Thus we limit the selection of false matches.

**2) Method 2: Scale Invariant Feature Transform (SIFT):**
The SIFT method [19] extracts interest points after series of treatments (pyramid of images, Differences of Gaussian, finding extrema). Then each point becomes a Keypoint by assigning a descriptor vector with 128 dimensions. This invariant vector represents the gradient norm for eight different orientations in each zone of a 4x4 window around the point. This vector is then used to associate a point in the first image with one of the points of the second image. The matching is done with an Euclidian distance calculus.

### B. Suppress false matching with the epipolar geometry

Previous methods present some limits and can therefore provide false matches ("outliers"). To reinforce this step, we used the epipolar constraint to verify all the matching points and calculate the fundamental matrix.

The RANSAC algorithm [20] can classify the matched data into good and false matches (inliers and outliers) using the calculation of the fundamental matrix. The algorithm select randomly eight points and estimate the fundamental matrix from the "eight points method" presented by Hartley in [21]. Then, the distance between the projected point in the second image and the epipolar line is calculated for each point. If this distance exceeds a threshold, the point is rejected and the remaining items are grouped in a set. The process is repeated and the set with the maximum of elements is selected. Finally, the fundamental matrix is estimated with these points.

Once the fundamental matrix estimated and the inliers obtained, the movement between two successive images can be deduced (e.g. section III-C).

## V. UNDERWATER ENVIRONMENT PROBLEMS AND CORRESPONDING FILTERING

The work presented in [16] demonstrates that under water, Harris detector with correlation gives better results than SIFT. However, for the vision community, SIFT is one of the most efficient detector in many applications. Comparisons have been performed on surface images but few for underwater pictures [22], [23].

In order to find which underwater parameters have failed SIFT algorithm, we propose to take surface picture and add noise filters that represent each of these disturbances.

Water introduced variables effects on the signal propagation decreasing the images quality. The main causes of this degradation is due to the refractive index, the presence of particles and the energy absorption of light ray according to their wavelength (color loss).

Some authors like [24], [25] have created global model more or less detailed of these disturbances, but none has decomposed every phenomena. The use of these filters do not permit us to determine what decreases SIFT performances.

Presence of suspended particles (turbidity) causes refraction and reflection phenomena on light rays. To model this phenomenon, we proposed to use a speckle noise. Indeed, this effect is often studied by physicist who characterize the suspended particles from the speckle noise. They used an image of the speckle noise generated by the underwater environment to determine some particles characteristics likes velocity, density, etc... [26], [27]. Therefore, turbidity can be represented by a speckle noise.

In an other hand, the changing of the refractive index introduces a blur effect and decreases the dynamic components on the images. To recreate this effect many filters can fit, however we chose to use a Gaussian low pass filter [28].

Regarding the energy absorption of light ray, as depth increases, colors drop off one by one depending on their wavelength. Low frequency color like red disappears at first ~3m. Then orange color is lost ~5m. Blue color travels the longest in water due to its higher frequency. To model this phenomenon, an attenuation on the amplitude of colors channels is used.

# VI. EXPERIMENTATIONS AND RESULTS

## A. Underwater experimentations

A camera/camcorder was used for these experiments. The resolution used for pictures extracted from the video is 1280 × 720 pixels. A waterproof case fitted the device to take underwater pictures obtained at sea in shallow water (15m). We immersed various objects (a statue of a female bust, a statue of a fish, bowls, plates and a calibration grid) which we are able to compare with 3D scanned models. We also took pictures of the same objects in a pool to limit some noises and in surface to apply filters and simulate underwater pictures. The results presented in this document are principally made with images of the female bust. In this experiment, different pictures were taken in order to cover all the object.

### B. Detection and matching results

As we have seen previously, the detection and matching step is an essential point of the 3D reconstruction. Thus, the two methods presented in Section IV have been examined and compared in different situations.

#### 1) Sea environment:

Initially, we tested these methods on sea images in order to determine the best method for this kind of pictures. Table I summarizes the results obtained from four pairs of images. The first three pairs of images are selected to have a small displacement. The last image pair shows a greater movement. The detections thresholds are adjusted to get roughly the same number of feature points. For each pair, we compared: the number of feature points detected in both images, the number of inliers and the inliers percentage.

These experiments show that Harris method gives better results and more points (17% more) than the SIFT algorithm,
<table><caption>TABLE I: Comparison of inliers number obtained with the two methods in sea environment</caption><thead><tr><th rowspan="2"></th><th colspan="2">Harris</th><th colspan="2">SIFT</th></tr><tr><th>Detected points</th><th>Inliers</th><th>Detected points</th><th>Inliers</th></tr></thead><tbody><tr><td>Pair 1</td><td>902 I1<br>902 I2</td><td>305<br>34%</td><td>902 I1<br>1122 I2</td><td>162<br>18%</td></tr><tr><td>Pair 2</td><td>797 I1<br>797 I2</td><td>244<br>30%</td><td>797 I1<br>894 I2</td><td>146<br>18%</td></tr><tr><td>Pair 3</td><td>894 I1<br>894 I2</td><td>310<br>35%</td><td>894 I1<br>999 I2</td><td>148<br>17%</td></tr><tr><td>Pair 4</td><td>1148 I1<br>1148 I2</td><td>291<br>26%</td><td>1505 I1<br>1148 I2</td><td>186<br>16%</td></tr></tbody></table>

contrary to what is usually found for surface images. In case of larger movements, the results are decreasing because we reach the limits of the correlation method which is valid only for small displacement.

## 2) Surface images and filtering:

The SIFT method does not work well under water, so we proposed to find which characteristic of the underwater environment is problematic. We used surface images on which proposed filters (Section V) are added to simulate underwater conditions. Both detection methods are tested with only one filter in a first time, then all filters are added to create an synthetic underwater image.

Fig. 2 represents the colors histogram of a sea, a surface and a filtered image. The first filter used to test these methods simulates the water color absorption phenomenon. To determine the proportion of color attenuation, a colors histogram of sea and surface images was studied. The red plane has completely disappeared in the sea picture. However, green and blue planes have suffered a decrease of 15% and 13% with respect to the surface image. Others pictures were compared and the same results were obtained. These results are consistent with theoretical attenuations at this depth (~15m). Indeed, the red color disappears completely around 3m and 5m and the green and blue colors are less reduced. For this experimentation, images were taken in the same orientation and resolution.

The parameters are used to adjust the color attenuation filter on all pictures. Fig. 3 represents the sea image (a) and the obtained image (b) with the color filter (the color is roughly the same between these two pictures).

Afterwards, we identified the variance of the blur and the speckle in the sea images. These two noises are estimated together for two reasons. First, these noises are physically linked. Second, measurement and exact identification of the core noise in the images are not achievable accurately. Therefore, we proposed to calculate the local variance of images and determine simultaneously the parameters of the Gaussian and speckle filter. We analyzed the histogram and calculated the average of the local variance of the sea image. Then, we adjusted the filter parameters to obtain the same value of variance between the sea and the filtered images. We have found a variance of 0.02 for the speckle noise and a variance of 8 for a windows size of [9 × 9] for the Gaussian filter. Fig. 4 shows the results achieved for sea, surface and filtered image with the parameters determined above from

<figure><img src="image_67.png" /><figcaption>(a) Red histograms</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>(b) Green histograms</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(c) Blue histograms</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>Fig. 3: Sea picture (a) and surface picture filtered with color filter (b) the local variance histogram.</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>Fig. 4: Histogram of a surface, a filtered and a sea image.</figcaption></figure>

The parameters of both filters defined, we have applied them to the surface image and obtained an image almost identical to the sea image (Fig. 5).

After the filter estimation, their effects were tested on the two methods of detection/matching points. At first, we used two images whose displacement corresponds to the fourth pair of images used in table I. Table II summarizes the results achieved for the different methods according to the applied filters.

As these experiments show, SIFT gives better results than Harris in similar situation to the surface such as blur or color
<figure><img src="image_72.png" /><figcaption>Fig. 5: Sea picture (a) and surface picture filtered with speckle and blur (b)</figcaption></figure>

<table><caption>TABLE II: Comparison of the two methods with different filters applied to surface images</caption><thead><tr><th rowspan="2"></th><th colspan="2">Harris</th><th colspan="2">SIFT</th></tr><tr><th>Detected points</th><th>Inliers</th><th>Detected points</th><th>Inliers</th></tr></thead><tbody><tr><td rowspan="2">No one</td><td>617 I1</td><td>154</td><td>625 I1</td><td>212</td></tr><tr><td>617 I2</td><td>25%</td><td>617 I2</td><td>35%</td></tr><tr><td rowspan="2">Color</td><td>629 I1</td><td>143</td><td>692 I1</td><td>223</td></tr><tr><td>629 I2</td><td>22%</td><td>629 I2</td><td>36%</td></tr><tr><td rowspan="2">Blur</td><td>631 I1</td><td>140</td><td>631 I1</td><td>227</td></tr><tr><td>631 I2</td><td>22%</td><td>666 I2</td><td>36%</td></tr><tr><td rowspan="2">Speckle</td><td>581 I1</td><td>139</td><td>581 I1</td><td>36</td></tr><tr><td>581 I2</td><td>24%</td><td>593 I2</td><td>6%</td></tr><tr><td rowspan="2">Speckle & Blur</td><td>631 I1</td><td>152</td><td>631 I1</td><td>91</td></tr><tr><td>631 I2</td><td>24%</td><td>641 I2</td><td>15%</td></tr><tr><td rowspan="2">All</td><td>583 I1</td><td>153</td><td>608 I1</td><td>84</td></tr><tr><td>583 I2</td><td>26%</td><td>583 I2</td><td>15%</td></tr></tbody></table>

change (due to a non natural lightning). But when speckle noise is added, the quality of results decreases dramatically with only 6% of inliers. This deterioration is due to the descriptor part of SIFT algorithm. The descriptor is based on oriented gradients. The speckle noise adds a "granularity" on images, which significantly disrupts the oriented gradients.

Subsequently, we tested different combinations of filters and we have seen that whenever the speckle noise is present in the image, SIFT results correspond to those that were expected from the results of sea pictures. We realized that Harris does not suffer practically any deterioration whatever the filters used. Inliers rate remains around 25% in all cases.

To check if this trend was general, we made the same experiments as before on a group of twenty images with different objects and backgrounds. The minimum, maximum and average of the two methods inliers percentages are computed in Table III.

<table><caption>TABLE III: Statistics of results of the two methods with different filters</caption><thead><tr><th rowspan="2"></th><th colspan="3">Harris Inliers</th><th colspan="3">SIFT Inliers</th></tr><tr><th>Min.</th><th>Max.</th><th>Mean</th><th>Min.</th><th>Max.</th><th>Mean</th></tr></thead><tbody><tr><td>No one</td><td>22.3</td><td>27.6</td><td>24.8</td><td>34.4</td><td>37.3</td><td>35.5</td></tr><tr><td>Color</td><td>22.7</td><td>26.6</td><td>24.1</td><td>33.6</td><td>38.2</td><td>35.7</td></tr><tr><td>Blur</td><td>21.4</td><td>27</td><td>24.5</td><td>34.8</td><td>38.9</td><td>36.4</td></tr><tr><td>Speckle</td><td>21.3</td><td>26.2</td><td>23.9</td><td>5.6</td><td>7</td><td>6.2</td></tr><tr><td>Speckle & Blur</td><td>23.1</td><td>26.7</td><td>24.6</td><td>13.8</td><td>15.3</td><td>14.5</td></tr><tr><td>All</td><td>24.1</td><td>27</td><td>25.7</td><td>12.1</td><td>15</td><td>14.1</td></tr></tbody></table>

We conclude that the descriptor of SIFT method is not robust to speckle noise and thus is not suitable for underwater pictures. Contrary to the Harris/correlation method which is robust to any noises present in the water.

### 3) Pool environment:

In section VI-B.2, we have highlighted the weaknesses of SIFT descriptors to the speckle noise. In this section, we propose to test the two algorithms on pool pictures. Indeed, pool images are subject to the same noise model as sea images but turbidity has almost zero in pool. This medium characteristic allow us to test the two methods in water with a low speckle noise.

We tested both algorithms on twenty pictures of each type and made statistics from results. Table IV regroups statistics of inliers percentage of pictures took in various environments. The air, filtered and sea statistics were obtained with twenty images.

<table><caption>TABLE IV: Statistics of results of the two methods in different conditions (surface, filtered, sea and pool)</caption><thead><tr><th rowspan="2"></th><th colspan="3">Harris Inliers</th><th colspan="3">SIFT Inliers</th></tr><tr><th>Min.</th><th>Max.</th><th>Mean</th><th>Min.</th><th>Max.</th><th>Mean</th></tr></thead><tbody><tr><td>Air</td><td>22.3</td><td>27.6</td><td>24.8</td><td>34.4</td><td>37.3</td><td>35.5</td></tr><tr><td>Filtered</td><td>24.1</td><td>27</td><td>25.7</td><td>12.1</td><td>15</td><td>14.1</td></tr><tr><td>Sea</td><td>21.4</td><td>28.1</td><td>25.1</td><td>11.9</td><td>16.7</td><td>15.2</td></tr><tr><td>Pool</td><td>23.8</td><td>28.5</td><td>26.1</td><td>26.2</td><td>30.7</td><td>28.4</td></tr></tbody></table>

These results reinforce our idea that speckle noise disrupts the SIFT descriptor. Indeed, the SIFT inliers rate is higher in pool than in sea. This phenomenon is due to the low number of suspended particles in pool and therefore to a weak speckle noise. Results remain lower than surface because noise is still present in pool pictures.

### C. 3D reconstruction results

In this part, we applied the method described in Section III. The reconstruction is performed from two images selected from a video. The choice of the two consecutive images (Fig. 6) is made so that the displacement between images is neither too large nor too small to ensure sufficient movement and prevent loss of information on the object.

<figure><img src="image_73.png" /><figcaption>Fig. 6: Underwater pictures used for 3D reconstruction.</figcaption></figure>

As we demonstrated above, the method of Harris / correlation allows us to get more robust points in the image than SIFT. Therefore, the fundamental matrix estimation is realized from Harris points and the essential matrix is calculated. Finally, the transformation between two pictures is deduced. Fig. 7 represents the relative movement between two successive images deduced from matched points with Harris and SIFT methods.

However, SIFT points are not neglected. Even if the number of points is lower, those founded are mostly coherent and we can remove the last outliers through the epipolar constraint. These points are different from those found by
<figure><img src="image_74.png" /><figcaption>Fig. 7: Feature points and estimate movement from Harris (a) and SIFT (b) method on the first image.</figcaption></figure>

the Harris method, but they enable to increase the model resolution. On fig. 7, we can see a few numbers of outliers with the 2 methods. These points will create errors in the 3D models, but we can easily remove them because their representation in 3D is aberrant with respect to the structure. In general, they are isolated from the structure and thus easy to identify and remove.

To obtain the 3D model, a triangulation is performed using the inliers and the projection matrix found previously. Then, the mesh structure is achieved through a 3D Delaunay triangulation. Fig. 8 shows the 3D wireframe model.

<figure><img src="image_75.png" /><figcaption>Fig. 8: 3D wireframe model.</figcaption></figure>

## VII. CONCLUSIONS AND FUTURE WORK

The objective of this work is to obtain robust matching points to reconstruct 3D underwater archaeological objects with a single camcorder. To achieve this goal, we have to make sure that the feature points and matching are robust enough to underwater conditions. Three filters are determined and quantified to characterize the underwater environment (color, blur and speckle). We have drawn a parallel between two detecting/matching methods through various experiments which enable us to single out the most robust method. We concluded that the SIFT descriptor is not robust to speckle noise unlike Harris/correlation method. In the future, we plan to apply a despeckle filter on images and check if the performance of SIFT is increased.

Then we used an euclidian reconstruction method to create a 3D model from these robust points and camera parameters. The 3D model obtained is sparse because it is created from only two images, therefore we plan to enrich it with a 3D dense matching method. This will increase the accuracy of the 3D object model. Then, texture will be added on the model.

## REFERENCES

[1] N. Snavely, S. M. Seitz, and R. Szeliski, "Photo tourism: exploring photo collections in 3d," *ACM Trans. Graph.*, vol. 25, pp. 835–846, July 2006.

[2] A. Bartoli and P. Sturm, "Constrained structure and motion from multiple uncalibrated views of a piecewise planar scene," *Int. J. Comput. Vision (IJCV)*, vol. 52, pp. 45–64, April 2003.

[3] L. Barazzetti, "A trip to rome: physical replicas of historical objects created in a fully automated way from photos," in *the First int. conf. on Human-computer interaction, tourism and cultural heritage, HCITOCH'10*, (Berlin, Heidelberg), pp. 63–72, 2011.

[4] V. Brandou, A. Allais, M. Perrier, E. Malis, P. Rives, J. Sarrazin, and P. Sarradin, "3d reconstruction of natural underwater scenes using the stereovision system iris," in *OCEANS 2007 - Europe*, pp. 1–6, 2007.

[5] P. Espiau, F.-X. Rives, "Extracting robust features and 3d reconstruction in underwater images," in *OCEAN 2001*, vol. 4, 2001.

[6] A. Kim and R. Eustice, "Pose-graph visual slam with geometric model selection for autonomous underwater ship hull inspection," in *Proceedings of IROS'09*, (Piscataway, USA), pp. 1559–1565, 2009.

[7] C. Beall, B. Lawrence, V. Ila, and F. Dellaert, "3d reconstruction of underwater structures," in *Proceedings of IROS'10*, december 2010.

[8] S. Negahdaripour and P. Firoozfam, "An rov stereovision system for ship-hull inspection," *IEEE J. of Oceanic Engineering*, vol. 31(3), pp. 551–564, July 2006.

[9] A. Hogue, A. German, J. Zacher, and M. Jenkin, "Underwater 3d mapping: Experiences and lessons learned," *Computer and Robot Vision, Canadian Conference*, vol. 0, p. 24, 2006.

[10] R. Hartley and A. Zisserman, *Multiple View Geometry in Computer Vision*. Cambridge University Press, second ed., 2003.

[11] W. K. Pratt, *Digital Image Processing: PIKS Inside*. New York, NY, USA: John Wiley & Sons, Inc., 2001.

[12] R. Szeliski, *Computer Vision: Algorithms and Applications*. 2010.

[13] T. Nicosevici, N. Gracias, S. Negahdaripour, and R. Garcia, "Efficient three-dimensional scene modeling and mosaicing," *J. of Field Robot*, vol. 26, pp. 759–788, October 2009.

[14] N. Pessel, J. Opderbecke, and M.-J. Aldon, "Camera self-calibration in underwater environment," in *WSCG*, 2003.

[15] J.-M. Lavest, G. Rives, and J.-T. Lapresté, "Underwater camera calibration," in *6th European Conf. on Computer Vision (ECCV)*, 2000.

[16] A. Meline, J. Triboulet, and B. Jouvencel, "A camcorder for 3d underwater reconstruction of archeological objects," in *OCEANS'10*, (Seattle Etats-Unis), pp. 1–9, MTS/IEEE, 2010.

[17] Z. Zhang, "A flexible new technique for camera calibration," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, pp. 1330–1334, 2000.

[18] C. Harris and M. Stephens, "A combined corner and edge detection," in *The Fourth Alvey Vision Conference*, pp. 147–151, 1988.

[19] D. G. Lowe, "Distinctive image features from scale-invariant key-points," *Int. J. Comput. Vision*, vol. 60, no. 2, pp. 91–110, 2004.

[20] M. A. Fischler and R. C. Bolles, "Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography," *Commun. ACM*, vol. 24, pp. 381–395, 1981.

[21] R. Hartley, "In defence of the 8-points algorithm," *IEEE Trans. On Pattern Analysis and Machine Intelligence*, vol. 19(6), pp. 580–593, June 1997.

[22] K. Oliver, W. Hou, and S. Wang, "Image feature detection and matching in underwater conditions," in *Ocean Sensing and Monitoring II - SPIE 7678*, 2010.

[23] A. Sarafraz, S. Negahdaripour, and Y. Schechner, "Performance assessment in solving the correspondence problem in underwater stereo imagery," in *OCEANS'10*, (Seattle Etats-Unis), MTS/IEEE, 2010.

[24] W. Hou, "A simple underwater imaging model," *Optics letters*, vol. 34, no. 14, pp. 2688–2690, 2009.

[25] J. S. Jaffe, "Monte carlo modeling of underwater-image formation: validity of the linear and small-angle approximations," *Appl. Opt.*, vol. 34, pp. 5413–5421, Aug 1995.

[26] R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P. K. Dixon, and D. J. Durian, "Speckle-visibility spectroscopy: A tool to study time-varying dynamics," *Review of Scientific Instruments*, vol. 76, no. 9, pp. 93–110, 2005.

[27] R. Adrian and J. Westerweel, *Particle Image Velocimetry*. Cambridge University Press, 2010.

[28] R. Schettini and S. Corchs, "Underwater image processing: state of the art of restoration and image enhancement methods," *EURASIP J. Adv. Signal Process*, vol. 2010, pp. 14:1–14:7, January 2010.
# Texture Features Selection for Visual Servoing of the Beating Heart

Aurélien Noce, Jean Triboulet, Philippe Poignet and Etienne Dombre

LIRMM - UMR 5506 CNRS-UM2
161 rue Ada
34392 Montpellier Cedex 5 - France
{noce,triboule,poignet}@lirmm.fr

**Abstract**— Texture features are widely used for image classification and retrieval. They offer an efficient way to add prior knowledge in image processing, using an appropriate prior parametrization.

In this paper, we present the methods we used to select among numerous existing features the most adapted to deal with beating heart tracking, using an experimental database of images and data mining techniques. Then we introduce the way we plan to use this information to reinforce region tracking algorithms in order to track the motion of characteristic landmarks on the heart surface.

**Index Terms**— beating heart, vision, visual servoing, texture

## I. INTRODUCTION

Recent advances in robotics - like *Da Vinci* from *Intuitive Surgical* - have demonstrated the potential of robotized tools [15] for biomedical applications. Surgical tasks are performed with those tools, in many hospitals around the world, and this may be the beginning of a new era of surgical applications. But, as versatile as those robots can be, there is still a lot of research to be done to extend their capacities. A current trend in medical robotics is to concentrate on very specific operations on a robotics point of view, and provide the surgeons with new possibilities of surgical interventions.

### A. Introducing beating heart surgery

This trend has risen up new stakes: beating heart surgery is one of those new research fields where robots seem very promising. The main issue in heart interventions is that it is extremely difficult for the surgeon to perform precise tasks on the heart because of the speed of its movements. So in conventional procedures such as coronary bypasses, the heart is stopped and blood circulation is ensured by what medical specialists call “the pump”. This procedure has several drawbacks: there is a great deal of risks [10] in stopping the patient’s heart, who can suffer from neurological injuries. Furthermore, this method is not compatible with the use of extra-corporeal circulation, as a full sternotomy is necessary to connect the equipments.

To limit those drawbacks, a solution might be to provide the surgeon with a motion compensation device. This approach to robotic assisted surgery has gained a growing interest in the scientific community for the last past years [8], [20] and work is done to apply it to the beating heart.

The motion compensation of the heart's movement is a very hard task for many reasons. Firstly, the heart's movement has very high dynamics - especially when compared to other organ's motion, such as respiration - and can become erratic depending on what is done during the operation. So predicting and tracking its motion is an open problem, and algorithms [1], [18] to deal with this issue are still in development.

Another issue is to find one or more appropriate sensor(s) to feed the control loop of the robot. The medical technologies offer many possibilities, but some are more constraining than others: for example the use of CT-scans or X-Rays raise a lot of questions about the security of the patient and the medical team exposed to the radiation, the resolution of the system...

### B. Computer Vision

To ensure the providing of the necessary information to the control algorithm, another interesting approach is to use computer vision. This choice is well adapted to mini invasive surgery, which is considered by many as the future of heart surgery. Use of vision is motivated by major considerations. The first is the compacity of a vision-based solution compared to large medical equipements such as CT-scanners or X-ray systems. As endoscopic equipement is already present on MIS equipements, this kind of solutions may easily be embedded in existing robots.

Anyway, fast visual servoing is a difficult task, very demanding in the quality of sensors and algorithms [5]. Work of Nakamura [17] was a first attempt to vision-based motion compensation of the heart's motion, and since other attempts [9] have shown interesting results.

<figure><img src="image_76.png" /><figcaption>Fig. 1. Artificial landmarks</figcaption></figure>
But there are several limitations in those approaches, one of them being the use of visual markers - see Fig.1 - like laser spots to track motion in the image. Our goal is to avoid the use of markers on the heat surface which requires to extract more information from the images in order to perform robust tracking. The main motivation for not using markers on the heart surface is the practical complexity of the placement of those markers, amounting with the scientific interest in exploring the possibilities of computer vision. Considering the sequences we could obtain from medical sources, texture analysis seemed to be a good candidate to extract such information. But as many texture characterization features were available in literature, we needed to study which of them were the most appropriate for the task, and decided to check their efficiency on our data.

The paper will be presented as follows : section II will introduce the different texture features we studied, section III will present the method used for data analysis, and finally sections IV and V will conclude our study.

## II. INTRODUCING TEXTURE FEATURES

In order to integrate texture features in our application, we decided to perform a comparative study of some commonly used texture features found in literature. Since the work of Julesz [12], [13] on texture decomposition and Haralick's [11] on texture classification, many new characterization methods have been proposed for texture discrimination. Evaluations of some of those methods have been published [4], [14], [21], but not for our specific purpose. The following texture features were found in literature and evaluated in the context of heart images.

### A. First-order statistics

Considering the image as a stochastic field, the first approach to texture characterization is to evaluate statistical operators on the gray-level values. The operators can be evaluated directly on the gray-levels histogram: they don't take into account the inter-pixel dependences, but only the frequencies of the gray-levels in the image. In this study we computed most usual statistics as mean, variance, and also skewness and kurtosis:

$$ skewness = \sum_{i=0}^{N-1} \frac{(x_i - \mu)^3}{N\sigma^3} \quad (1) $$

$$ kurtosis = \sum_{i=0}^{N-1} \frac{(x_i - \mu)^4}{N\sigma^4} \quad (2) $$

### B. Gray Level Cooccurrence Matrix

The gray-level co-occurrence matrix - GCLM - is one of the most commonly used approach for texture characterization, and has been used successfully by many computer vision applications [2], [7]. Introduced by Haralick [11], it describes the relation between neighbor pixels in the image. For a given image, different co-occurrence matrices can be computed for a specific angle $\theta$ and distance $d$. Fig.2 illustrates how the GCLM is computed. For the study we

<figure><img src="image_77.png" /></figure>

Original image

<figure><img src="image_78.png" /></figure>

Computation for $\theta = 0$ and $d = 1$

<figure><img src="image_79.png" /></figure>

GLCM

<caption>Fig. 2. Computation of a GLCM with $d = 1$ and $\theta = 0$</caption>

$$ energy = \sum_{i,j} p(i,j)^2 \quad (3) $$

$$ correlation = \sum_{i,j} \frac{(i - \mu_i)(j - \mu_j)}{\sigma_i \sigma_j} p(i,j) \quad (4) $$

$$ variance = \sum_{i,j} (i - \mu)^2 p(i,j) \quad (5) $$

$$ homogeneity = \sum_{i,j} \frac{1}{1 + (i - j)^2} p(i,j) \quad (6) $$

$$ entropy = - \sum_{i,j} p(i,j) \ln(p(i,j)) \quad (7) $$

<caption>TABLE I
SOME MAJOR HARALICK TEXTURE FEATURES</caption>

used a normalized GLCM, which is less sensitive to global illumination.

Co-occurrence matrices are directional, so many GLCMs are computed along a single image, for different directions and distances, in order to obtain an almost-isotropic analysis. Strategies have been developed [3] to reduce the necessary computing. GLCMs are usually symmetrized, and the number of considered gray-levels is reduced. This last point is very important as the size of the co-occurrence matrix is $N_g \times N_g$, $N_g$ being the number of gray-levels. Considering less gray-levels can speed up the calculations, but may also reduce the efficiency of the texture features.

Once the GCLM is computed, we evaluate several texture features: Haralick's texture features presented on Tab.I<fnref n="1" />, but also classical features like cluster shade, cluster prominence, etc. All those features are evaluated for different resolutions, directions and distance.

### C. Gray Level Difference Matrix

This method, introduced by Weszka et al [21], has been widely use for analyzing medical images. The gray-level difference matrix is obtained by subtracting the original image with itself translated by a given displacement vector $d = (d_x, d_y)$:

$$ GLDM(i, j) = |I(i, j) - I(i + d_x, j + d_y)| \quad (8) $$

$$ p(k) = P(GLDM(i, j) = k) \quad (9) $$

Then statistics are evaluated from the frequencies $p(k)$ of gray-levels $k$ in the image. This approach provides very fast computation of inter-pixel statistics, and this can be even improved by reducing the number of gray-levels. Nevertheless, as with the co-occurrence matrix, it has to be computed in several directions for each image to

<footnote n="1">¹see [11] for details</footnote>
<figure><img src="image_80.png" /><figcaption>Fig. 3. Exemple computation of a GLDM</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>Fig. 4. run-length matrix computation</figcaption></figure>

give global information. Fig.3 illustrates the construction
process of the GLDM on a basic example.

Several features are evaluated on this matrix, among
which :

$$
contrast = \sum_{i=1}^{N_g} i^2 \cdot p_d(i) \qquad (10)
$$

$$
uniformity = \sum_{i=1}^{N_g} [p_d(i)]^2 \quad (11)
$$

$$
entropy = \sum_{i=1}^{N_g} p_d(i) \log(p_d(i)) \quad (12)
$$

$p_d(i)$ representing the probability of a pixel's luminance
being $i$ and $N_g$ the number of gray-levels considered.

D. Image filtering

Convolutions with several structuring elements were
performed to provide selective spectral information. 3 × 3
pixels four-neighborhood elements and also 5 × 5 pixels
elements morphological tools performs low pass filtering
on the image : This makes our analysis less influenced by
local discontinuities. Once the operations are performed we
evaluate the *mean* and *variance* of the modified gray-levels
matrices.

E. Run-Length Matrix

The run-length matrix stores the number of runs of gray-level *i*, of length *j*, in a specified direction. This approach gives higher-order statistics than co-occurrence matrix or gradients method. It involves a high amount of computation, but has been used with success in medial image computing [19], [22]. The Run-Length matrix computation is illustrated on Fig.4 : each color represents a run in the direction θ = 0° - the horizontal direction - and the Run-Length Matrix represents the count of runs for each possible gray-level.

We studied the efficiency of several classical features
evaluated on the run-length matrix. Tab.II explicits those
computations.

$$
\text{Short Run Emphasis} = \frac{1}{N_r} \sum_{i,j} \frac{r(i,j)}{j^2} \quad (13) \\
\qquad (N_r \text{ number of runs})
$$

$$
\text{Long Run Emphasis} = \frac{1}{N_r} \sum_{i,j} j^2 \cdot r(i,j) \quad (14)
$$

$$
\text{Gray-Level Non-Uniformity} = \frac{1}{N_r} \sum_{i} \left[ \sum_{j} r(i,j) \right]^2 \quad (15)
$$

$$
\text{Run Length Non-Uniformity} = \frac{1}{N_r} \sum_j \left[ \sum_i r(i,j) \right]^2 \quad (16)
$$

$$
\text{Run Percentage} = \frac{N_r}{N_p} \quad (17) \\
\text{(}N_p \text{ number of pixels)}
$$

<figure><img src="image_82.png" /><figcaption>Fig. 5. Example computation of Principal Axes</figcaption></figure>

F. Fourier transform

Spectral analysis has been performed on our data to
provide another approach to texture analysis. We used FFT
algorithm to compute the spectral components, but our
texture features are mainly based upon the power spectrum,
on which we evaluated energy, mean value, maximum
value, amplitude and variance.

III. COMPARISON OF TEXTURE FEATURES

A. *Principal Component Analysis*

Principal component analysis is a statistical method
for studying multidimensional data. It is widely used in
data mining because it eases the interpretation of high
dimensional data. PCA has been used in computer vision
to perform shape and object recognition [6], [16] and is
generally useful in pattern recognition. Our focus here is
to reduce the number of attributes in order to find which
of our texture features are best suited for processing on
heart images, PCA will then give us information about
the correlation between our variables and their ability to
discriminate our sample images.

Technically, PCA finds a rotation in the representation
space that maximizes the variance obtained along the first
axes. Fig.5 illustrates the extraction of Principal compo-
nents - represented by the two perpendicular lines - on two-
dimensional data. By doing so we reduce the correlation
between the projection axes and we can reduce the dimen-
sions of the problem. To find the principal components, the
covariance matrix is computed :
<figure><img src="image_83.png" /><figcaption>Fig. 6. Some of our test-images</figcaption></figure>

<figure><img src="image_84.png" /><figcaption>Fig. 7. Cumulative histogram of variances</figcaption></figure>

$$C_{i,j} = \frac{1}{N-1} \sum_{k=1}^{N} (X_{k,i} - \mu_{i})(X_{k,j} - \mu_{j}) \quad (18)$$

where X is the data matrix, n the number of variables, C the covariance matrix which has variances of the data on its diagonal and co-variances between data elsewhere. Then the theory tells us that the eigenvector of the covariance matrix - once reordered with decreasing eigenvalues - give us the new base of PCA.

Once the result of PCA is known, we can evaluate the quality of the representation among each axis using the coefficient :

$$q = \frac{\sum_{i=1}^{N} \lambda_{i}}{\sum_{i=1}^{N_{axes}} \lambda_{i}} \quad (19)$$

This coefficient helps in selecting how many Principal Components must be chosen to have a good representation of the dataset i.e. to choose the dimension of the decision space.

## B. Input Data

In order to generate the data for our analysis, we computed the feature vector from a set of 200 images. Those images come from experiments made at the experimental laboratory of the Medical University in Montpellier, and were selected to represent the regions of interest with different conditions of illumination. They have been classified into different sets before the study, so we can compare the result of our automatic classification procedure with the man-made classification. Our goal was to find a minimal set of features able to distinguish between the average surface of the heart and specific regions of interest that have a specific texture : the veins networks visible on the heart. Some of those manually-selected test-images are presented on Fig.6.

The feature vector is initially very large : we compute 117 features for each test-image. We perform Principal Component Analysis - Fig.7 and 9 - on the data. The following analysis is based on the result of this PCA.

## C. Features selection

The cumulative histogram of variances - Fig.7 - shows the percentage of total variance that is represented on the first Principal Components. It represents the Principal Components that hold more than 95% of total variance : here only 9 PCs are needed to represent nearly all of the

total variance of our dataset. As we can see on Fig.7, the first Principal Axis holds almost half of the total variance, and first three Principal Axes hold about 80% of the total variance : the major part of the discriminative potential of our features is concentrated in this 3-dimensional space.

1) First selection procedure: In order to select the best features to retain the overall variability, an approach was to study the contribution of the different features on the first three principal components. As the PCA performs a rotation of data X from original features basis to the principal components basis, we looked at the rotation matrix R - X<sub>PCA</sub> = R.X - and the coefficients of R express the contributions of the different features to the Principal Components. By looking at the R coefficients, we highlighted that a few of our feature had a very high influence on the Principal Components, which confirms the first impression given by the histogram of variance. We extracted a 28 features set which corresponds to the features necessary to evaluate more than 90% - in norm - of the first three PCs.

To complete the analysis we recomputed a PCA using those components as the new basis. The results of this PCA are given on Fig.8. Looking at the two-dimensional projections<fnref n="2" />, variability of data appears to have been preserved as there is no excessive grouping of plots and groups homogeneity remaining. On the other hand the cumulative variance histogram still shows a predominance the first Principal Component.

2) Second selection procedure: Another very informative representation is the correlation circle - Fig.9 - which shows the correlation<fnref n="3" /> of the features with the first two Principal Components, and consecutively shows the inter-correlation between the Principal Axes. The variables that appear nearby the circle are highly inter-correlated, while the one close to the origin are more linearly independent. For example the components in the very right of Fig.9 are highly inter-correlated, while the set of features closer to the center in the left of the figure is more independent.

As we try to reduce our features set and still preserve the information, we selected a set uncorrelated features in

<footnote n="2">²the different kind of plots on the projections represent the different preselected groups in our study</footnote>

<footnote n="3">³the indexes on the correlation circle are the indexes of the features in our study</footnote>
<figure><img src="image_85.png" /><figcaption>Fig. 8. PCA on our data with a reduced feature vector</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Fig. 9. the correlation circle</figcaption></figure>

order to have uncorrelated features. We choose 32 textures attributes, corresponding to the group visible on the upper left of Fig.9. The result of Principal Component Analysis that you can see on Fig.10 shows that we still obtain a good distribution of the features on the projections. A noticeable characteristic of this set of feature is the good repartition of the variance among the components, as you can see on the histogram of Fig.10 which proves that we have reduced the correlation between our variables.

### D. The predominant features

Considering the new set of features, we observe the predominance of *Gray-Level Co-occurrence Matrix* and *Run-Length matrix* approaches. The Run Length method is more represented in set 1, while the Gray Level Co-occurrence Matrix features are more present in the set with uncorrelated PCs.

Some of the features, like *cluster shade* and *cluster prominence* for the GLCM, or *Short Run Low Grey-Level Run Emphasis* and *Gray-Level Non-Uniformity* for the Run-Length method hold an important part of the overall entropy. Other features like *skewness* and *kurtosis* showed interesting possibilities, appearing in both sets of selected features, but first-order statistics stay behind higher order approaches.

On the other hand, Image Filtering did not show good results, and the FFT analysis showed mixed results : feature *mean of the power spectrum* is a part of the first set of

<figure><img src="image_87.png" /><figcaption>Fig. 10. PCA on our data with a reduced feature vector</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>Fig. 11. Entropy of the feature sets from the correlation circle method</figcaption></figure>

features and has a good contribution in the first Principal Component, but the other features appeared to be correlated with other variables.

## IV. EVALUATION OF THE NEW FEATURE VECTORS

### A. Entropy analysis

For a representation with $N_i$ images and $N_f$ features, the entropy is defined from the eigenvalues $\lambda_i$ of the $N_f \times N_f$ covariance matrix by the relations :

$$ 
\begin{gathered}
H = \bar{\lambda}_i \log \bar{\lambda}_i \qquad (20) \\
\text{with } \bar{\lambda}_i = \frac{\lambda_i}{\sum_{j=1}^{N_f} \lambda_j}
\end{gathered}
$$

An high $H$ value means that entropy is well balanced among the axis, so the best set will have the higher $H$ value. The Tab.III shows the comparison between the entropy of the feature sets introduced in section III-C : set 1 is the set obtained from the analyzes of the rotation operator and set 2 is the set obtained with the correlation circle . It also shows the maximum value for the correlation circle method : as seen on Fig.11, the entropy is maximum for the set of the 8 less intercorrelated features, called set 3 in the table.

As we can see on Tab.III, the first set does a better job in extracting entropy than the correlation-circle based methods do.
<table><caption>TABLE III
ENTROPY VALUES FOR THE SETS OF OUR STUDY</caption><thead><tr><td>Set</td><td>H value</td></tr></thead><tbody><tr><td>All features</td><td>3.0713e-017</td></tr><tr><td>set 1</td><td>0.6704</td></tr><tr><td>set 2</td><td>0.0011</td></tr><tr><td>set 3</td><td>0.2319</td></tr></tbody></table>

<table><caption>TABLE IV
AUTOMATIC CLASSIFICATION RESULTS</caption><thead><tr><td>Set</td><td>correctly classified</td></tr></thead><tbody><tr><td>All features</td><td>72%</td></tr><tr><td>set 1</td><td>88%</td></tr><tr><td>set 2</td><td>66%</td></tr></tbody></table>

As those vectors were designed to help distinction between regions of the heart, we tested our new feature vectors on a set of heart images to evaluate their ability to perform a correct classification of our images into two set: the images of the regions of interest and the ones of the rest of the heart surface. We performed geometric discrimination using the eulerian between our feature vectors. This approach is well adapted to integrate later in a region tracking algorithm.

Finally Tab.IV sums up the results of the classification process. The classification is not perfect but considering the complexity of the heart texture those results are encouraging. The first set of features still shows good results compared to the second one.

## V. CONCLUSION AND FUTURE WORK

In this paper, we presented a comparison between some commonly used texture features applied to beating heart images. The issue of defining the best texture vector for performing visual servoing is difficult, as the solution is not unique and the criteria of evaluation numerous. In this work we focused on the quality of representation more than the computational efficiency. This is a deliberate choice, as the security of the patient may depend on the quality of the servoing. On the other hand selecting the good attributes also contributes in decreasing computational need and improves overall performance.

The results we obtained show some big differences between the attributes. This reinforces our conviction that texture tools are not universal and must be used selectively. The results on PCA showed the predominance of co-occurrence and run-length approaches, while entropy analysis and classification tests showed the advantages of reduced sets of features. We can see that inter-pixel information is primordial in analysis heart textures, and only high order statistics manage to perform efficient recognition.

Future work includes the integration of our feature vector in the region tracking process, with joint use of correlation methods and our texture vector. Medical visual servoing is a very complex task, and only extensive testing will determine if our approach really leads to robust estimation.

## REFERENCES

[1] M. Cenk Cavusoglu, Jason Rotella, Wyatt S. Newman, Sangeun Choi, Jeff Ustin, and S. Shankar Sastry. Control algorithms for active relative motion cancelling for robotic assisted off-pump coronary artery bypass graft surgery. In *Proc. 12th Int. Conf. on Advanced Robotics*, Seattle, WA, USA, jul 2005.

[2] C.-W. CHAN, N. LAPORTE, and R. S. DEFRIES. Texture classi?cation of logged forests in tropical Africa using machine-learning algorithms. *NT. J. REMOTE SENSING*, 24(6):1401–1407, 2003.

[3] Clausi and Jernigan. A fast method to determine co-occurrence texture features. *IEEE Trans. on Geosciences and Remote Sensing*, 36(1):298–300, 1998.

[4] David A. Clausi. Comparison and fusion of co-occurrence, gabor and mrf texture features for classification of sar sea-ice imagery. *ATMOSPHERE-OCEAN*, 39(3):183–194, 2001.

[5] Loïc Cuvillon, Edouard Laroche, Jacques Gangloff, and Michel de Mathelin. Gpc versus $h_\infty$ control for fast visual servoing of a medical manipulator including flexibilities. In *ICRA*, 2005.

[6] Bruce A. Draper, Kyungim Baek, Marian Stewart Bartlett, and J. Ross Beveridgea. Recognizing faces with pca and ica. *Comp. Vision and Image Understanding*, 91:115–137, 2003.

[7] Aria E.H., Saradjian M.R., Amini J., and Lucas C. Generalized cooccurrence matrix to classify irs-1d images using neural network. In *Proc. XXth ISPRS Congress*, pages 117–122, Istanbul, Turkey, jul 2004.

[8] R. Ginhoux, J. Gangloff, M. de Mathelin, L. Soler, M. M. Arenas Sanchez, and Jacques Marescaux. Active filtering of physiological motion in robotized surgery using predictive control. *IEEE Trans. Robotics*, 2003.

[9] R. Ginhoux, J. A. Gangloff, M. F. de Matheli, L. Soler, Mara M. Arenas Sanchez, and J. Marescaux. Beating heart tracking in robotic surgery using 500 hz visual servoing. In *Medical Robotics, Navigation and Visualization*, Remagen, Germany, mar 2004.

[10] Mekhann GM, Goldsborough MA, and Borowicz LM. Predictors of stroke risk in coronary artery bypass patients. *Ann Thorac Surg*, 63:516–521, 1997.

[11] Haralick, R.M., Shanmugam, K., and Dinstein I. Textural features for image classification. *SMC*, 3(6):610–621, November 1973.

[12] B. Julez. Visual pattern discrimination. *IRE Trans. on Information Theory*, 8:84–92, 1962.

[13] B. Julez. Experiments in the visual perception of texture. *Scientific American*, 232(4):2–11, 1975.

[14] P. Kruizinga, N. Petkov, and S.E. Grigorescu. Comparaison of texture features based on gabor filters. In V. Roberto et al., editor, *Proc. 10th Int. Conf. Image Analysis and Processing*, pages 142–147, sep 1999.

[15] A. P. Kypson and Chitwood V. R. Jr. Robotic applications in cardiac surgery. *Int. Journal on Advanced Robotic Systems*, 1(2), 2004.

[16] Fernando De la Torre and Michael J. Blacky. Robust principal component analysis for computer vision. In *Int. Conf. on Computer Vision*, Vancouver, Canada, jul 2001.

[17] Y. Nakamura, K. Kishi, and H. Kawakami. Heartbeat synchronisation for robotic cardiac surgery. In *Proc. 2001 Int. Conf. Robotics and Automation*, Seoul, Korea, may 2001.

[18] T. Ortaiaer. Motion estimation in beating heart surgery. *IEEE Trans. BME*, 52(10), oct 2005.

[19] D. S. Raicu, J. D. Furst, D. Channin, D. H. Xu, and A. Kurani. A texture dictionary for human organs tissues’ classification. In *Proc. 8th World Multiconf. on Systemics, Cybernetics and Informatics*, Orlando, USA, jul 2004.

[20] A. Schweikard, G. Glosser, M. Bodduluri, M. Murphy, and J. R. Adler. Robotic motion compensation for respiratory motion during radiosurgery. *J. Comp. Aided Surgery*, 2000.

[21] Weszka, J. S., C. R. Dyer, and A Rosenfield. A comparative study of texture measures for terrain classification. *Trans. Systems, Man, and Cybernetics*, SMC-6(4):269–285, 1976.

[22] D.-H Xu, A.S. Kurani, J.D. Furst, and D.S. Raicu. Run-length encoding for volumetric texture. In *4th Int. Conf. on Visualization, Imaging and Image Processing*, 2004.