Barbara Kuijper
email: kuijperb@maasstadziekenhuis.nl

Jos J Th

Bas F Tans

Frans Van Der Kallen

Geert J Nollet

Lycklama

Marianne Nijeholt

De Visser

Tans

Bas F Van Der Kallen

Frans Nollet

Geert J Lycklama

Marianne De Visser

Root compression on MRI compared to clinical findings in patients with recent onset cervical radiculopathy

Keywords: Cervical radiculopathy, Magnetic resonance imaging (MRI), disc herniation, spondylotic foraminal stenosis

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Root compression without clinical substrate as a coincidental finding on MRI of the cervical spine is well known. [START_REF] Jensen | Prolapsed cervical intervertebral disc in male professional drivers in Denmark, 1981-1990. A longitudinal study of hospitalizations[END_REF][START_REF] Van Rijn | Observer variation in the evaluation of lumbar herniated discs and root compression: spiral CT compared with MRI[END_REF][START_REF] Van Rijn | Observer variation in MRI evaluation of patients suspected of lumbar disk herniation[END_REF][START_REF] Boden | Abnormal magnetic-resonance scans of the lumbar spine in asymptomatic subjects. A prospective investigation[END_REF] However, presence of a clinically evident cervical radiculopathy without root compression on MRI does also occur. [START_REF] Fortin | No clear winner in differing imaging modalities for cervical radiculopathy[END_REF][START_REF] Bartlett | A comparison of T2 and gadolinium enhanced MRI with CT myelography in cervical radiculopathy[END_REF][START_REF] Birchall | Evaluation of magnetic resonance myelography in the investigation of cervical spondylotic radiculopathy[END_REF] Knowledge of the occurrence of this phenomenon is important to take correct treatment decisions.

The most common causes of cervical root compression are narrowing of the foraminal space secondary to spondylarthrosis and herniated discs. [START_REF] Tanaka | The anatomic relation among the nerve roots, intervertebral foramina, and intervertebral discs of the cervical spine[END_REF][START_REF] Adams | Cervical spondylotic radiculopathy and myelopathy[END_REF][START_REF] Radhakrishnan | Epidemiology of cervical radiculopathy. A population-based study from Rochester, Minnesota, 1976 through 1990[END_REF] Large studies investigating the relationship between the radiological and clinical findings are lacking.

In this study we assess the relation between the clinically affected level of root compression and the level of root compression on magnetic resonance imaging (MRI) in a group of patients with recent onset cervical radiculopathy.

Methods

Subjects

In our prospective cohort study we included patients with a clinical diagnosis of unilateral cervical radiculopathy in whom efficacy of either a cervical collar or physiotherapy was compared with a wait and see policy. [START_REF] Kuijper | Cervical collar or physiotherapy versus wait and see policy for recent onset cervical radiculopathy: randomised trial[END_REF] The diagnosis cervical radiculopathy was made by a neurologist according to the following inclusion criteria: radiation of arm pain distal to the elbow, symptoms for less than one month, arm pain on a visual analogue scale of 40 mm or more, plus at least one of the following: (1) worsening of the arm pain by neck movements, (2) sensory symptoms in one or more adjacent dermatomes, (3) diminished deep tendon reflexes in the affected arm, or (4) muscle weakness in one or more adjacent myotomes. The clinically affected level was determined by the neurologist via structured history taking and physical examination data.

Written informed consent was obtained from all patients. The medical ethics committees of the participating hospitals approved the protocol.

MRI protocol

MRI was performed at 1.5 Tesla, and included sagittal proton density (PD) and T2 weighted turbo spin-echo imaging (TSE; TR/TE: 2900/23/182; echo train length (ETL): 5) and T1 weighted TSE (TR/TE 664/13; ETL: 3), using 3 mm thick slices and pixels of 1 mm 2 . Axial imaging consisted of a T2 weighted MEDIC sequence (TR/TE 1140/27) and T1 TSE (TR/TE 538/13, ETL: 3) using 3 mm thick slices and 1 mm 2 pixels. Axial slices were placed perpendicular to the vertebral bodies and were angulated in the same way. Axial imaging covered intervertebral spaces C4 to T1. Per patient four levels were investigated on both sides: C4-5 (root C5), C5-6 (root C6), C6-7 (root C7) and C7-Th1 (root C8).

MRI evaluation protocol

The MRI examinations were evaluated independently by two neuroradiologists who were blinded to clinical findings. For each cervical level, presence of herniated disc and bony foraminal stenosis by spondylarthrosis was assessed.

The probability of root compression, either caused by herniated disc or spondylotic foraminal stenosis, was scored on a five-point scale: 'definitely no root compression', 'possibly no root compression', 'indeterminate', 'possibly root compression' and 'definitely root compression'.

Presence of other abnormalities like spinal canal stenosis with or without cord compression, tumors, and other findings was recorded but not included in the analysis.

Analysis

Data on root compression were dichotomized as either root compression (possibly and definitely root compression) or no root compression (definitely and possibly no root compression and indeterminate).

For each patient we compared the level of root compression with that of the clinically affected root. As it is not always feasible to determine with certainty the level of root compression by history and physical examination, [START_REF] Kuijper | Degenerative cervical radiculopathy: diagnosis and conservative treatment. A review[END_REF] we separately recorded the presence of root compression one level higher or lower. We also examined the cause of compression: herniated disc, spondylotic foraminal stenosis or both. Data on the clinical signs of radiculopathy in these 78 patients are listed in table 1. Sensory abnormalities were found in 89.9% of the patients, diminished reflexes in 48.1% and muscle weakness in 29.1%. Two patients (2.6%) only experienced pain. The mean age was 47.4 year (range 24-70). There was no difference in age between patients with herniated discs and spondylotic foraminal stenosis.

Results

From

Discussion

To the best of our knowledge this is the first study evaluating the correlation between clinical signs and symptoms and the occurrence of root compression on MRI in a well defined population of patients with recent onset unilateral cervical radiculopathy, in a setting closely resembling clinical practice.

In almost three quarter of the 78 patients the clinical level corresponded with the level of root compression on MRI. Assessment of the level of root compression by neurological examination alone is reported to be difficult. [START_REF] Fisher | Electrophysiology of radiculopathies[END_REF][START_REF] Levine | Cervical Radiculopathy: Diagnosis and Nonoperative Management[END_REF][START_REF] Yoss | Significance of symptoms and signs in localization of involved root in cervical disk protrusion[END_REF][START_REF] Ozgur | Atypical presentation of C-7 radiculopathy[END_REF][START_REF] Slipman | Provocative cervical discography symptom mapping[END_REF][START_REF] Narayan | Treatment of degenerative cervical disc disease[END_REF] Inclusion of the root compressions one level above or below the clinically suspected level, led to an increase of the percentage of matching clinical and MRI levels by only 5 percent for both radiologists indicating that the localising value of the neurological examination is quite good.

Strikingly, we found false positive root compressions in 45 percent of the cases, defined as compressions on the contralateral, i.e., asymptomatic, side or at least two levels higher or lower than the clinically affected level. These asymptomatic root compressions were often present together with root compression at the clinically affected level or adjacent to this. In order to avoid unnecessary treatment of root compressions observed on MRI without clinical substrate, physicians should only interpret MRI results after a careful history and neurological examination.

Earlier studies reported a high percentage of cervical degenerative disease on MRI in asymptomatic patients, with a higher prevalence in older age. Teresi et al studied 100 MRI's of the cervical spine of patients who were investigated for laryngeal disease. Twenty percent of patients aged 45-54 years and 57 percent of patients older than 64 years had cervical disc protrusion without clinical symptoms. [START_REF] Teresi | Asymptomatic degenerative disk disease and spondylosis of the cervical spine: MR imaging[END_REF] Another study on asymptomatic subjects showed 'major' abnormalities (herniated disc, foraminal stenosis, disc space narrowing) in 28 percent of people over 40 years of age. [START_REF] Boden | Abnormal magnetic-resonance scans of the lumbar spine in asymptomatic subjects. A prospective investigation[END_REF] However, these studies only investigated asymptomatic patients, whereas our study included patients based on the presence of cervical radiculopathy at one level.

We also found a fair amount of false negative MRI's. MRI's were assessed as normal in 12.8 /15.4% of cases and in 9.0/10.3 % only asymptomatic roots were compressed. So, together we had 21.8/25.6 percent of false negative MRI's.

Patients were eligible for our study if they fulfilled strict inclusion criteria. We feel confident that they had indeed cervical radiculopathy because all patients except two had neurological deficits corresponding with radicular pain, and no other conditions emerged during a 6-months follow-up period. We had low percentages of surgery so we could not confirm our diagnosis surgically in most cases. [START_REF] Kuijper | Cervical collar or physiotherapy versus wait and see policy for recent onset cervical radiculopathy: randomised trial[END_REF] The cause of the root compression was most often herniated disc, particularly in those patients with unilateral MRI abnormalities only at the clinically affected level. This is noteworthy since it is often assumed that spondylotic foraminal stenosis is the most common cause of root compression in cervical radiculopathy. 8-10 12 Our study shows that foraminal stenosis was more frequently asymptomatic. The occurrence of herniated discs and spondylotic changes was not related to age.

We previously found a high interobserver agreement of 91% for MRI evaluation of root compression with a kappa score of 0.67 . This agreement was less for the cause of the compression, i.e., herniated disc (81%) and spondylotic foramen stenosis (82%) (unpublished data, manuscript accepted for publication in Clinical Radiology). In the present study the two neuroradiologists also disagreed more on assessment of herniated discs and foraminal stenoses than on the presence of root compression. Particularly the reported percentages of spondylotic foraminal stenosis differed (table 3). It is known that on MRI spondylotic foraminal stenosis is often more difficult to detect and CT-myelography techniques are probably more accurate. 5-7 20 Our neuroradiologists may have underreported spondylotic foraminal stenosis, although the 55 and 63 percent of patients in whom stenosis was found by radiologist 1 and 2 suggest otherwise (unpublished data, manuscript accepted for publication in Clinical Radiology).

Another limitation is that the neuroradiologists knew that all patients in the study had a clinical diagnosis of cervical radiculopathy. They may have been more aware of the possibility of root compression in general, resulting in higher percentages abnormal MRI's.

Because the radiologists had no information on the level and side of the radiculopathy the main results of our study on the relation between clinical and MRI findings seem reliable.

Conclusion and clinical implications

It is evident from this study that MRI findings in patients with cervical radiculopathy are only meaningful in a clinical context. False negative MRI results were encountered in almost one quarter, and false positive results in half of the patients. Therefore, cervical MRI is only useful when there is a clear picture of cervical radiculopathy.

  May 2005 to December 2006 82 patients with a clinical diagnosis of recent onset unilateral cervical radiculopathy underwent MRI of the cervical spine. Four scans were of poor quality because of movement artifacts by pain (n=2) and claustrophobia (n=2). 78 MRI's were fully examined.

Table 1

 1 Clinical signs and symptoms in 78 patients with suspected cervical radiculopathy

		78 patients of total
		group
		No (%)
	only pain	2 (2.6)
	muscle weakness	23 (29.1)
	diminished reflexes	38 (48.1)
	sensory abnormalities	70 (89.9)

In table 2 single and multiple root compression assessments are listed. The clinically affected root was judged as compressed in 73.1/69.2% (radiologist 1/2) of all cases. When we consider compression of the roots adjacent to the clinically affected level as symptomatic root compression, this percentage increased to 78.2/74.4%.

Single level, one sided root compression was found in 47.4% of patients. Most of these compressions were of the clinically affected root, according to table 2 in 35.9/34.6% (radiologist 1/2). Multiple root compression was seen in 39.7/37.2% of cases (radiologist 1/2), mostly a combination of clinically affected and unaffected roots. One or more clinically unaffected root compressions were reported in 44.9/46.2% of all patients. MRI's were assessed as normal in 12.8/15.4% of cases and in 9.0/10.3% only unaffected roots were compressed. We additionally recorded the clinical signs listed in table 1, for the group of patients on whom the 2 neuroradiologists agreed that the MRI showed no root compression at the clinical level or one level above or below. No difference was found between the clinical data of patients with and without root compression on MRI.

Table 2

 2 Single and multiple root compression assessed on MRI in 78 patients with suspected cervical

	radiculopathy
	total n=78	radiologist 1	radiologist 2
	Single root compression on MRI	n=37 (47.4)	n=37 (47.4)
	Compression of clinically affected root only	28 (35.9)	27 (34.6)
	Compression of the root one level higher or lower than	3 (3.8)	2 (2.6)
	that of the clinically affected root
	Compression of the root two levels higher or lower	1 (1.3)	1 (1.3)
	Compression of contralateral root	5 (6.4)	7 (9.0)
	Multiple root compression on MRI	n=31 (39.7)	n=29 (37.2)
	Compression of clinically affected root together with	2 (2.6)	1 (1.3)
	compression one level higher or lower
	Compression of clinically affected root and
	contralateral root compression	18 (23.1)	13 (16.7)
	Compression of clinical affected root together with
	one level higher or lower and contralateral root	9 (11.5)	13 (16.7)
	compression
	Compression of root one level higher and lower
	together with contralateral root compression	1 (1.3)	2 (2.6)
	Compression of multiple contralateral roots	1 (1.3)	0
	Normal MRI	10 (12.8)	12 (15.4)
	Herniated discs were the single cause of root compression in 41.1/34.6% and spondylotic
	foraminal stenosis in 35.9/21.7% of our patients (table 3). A combination of herniated disc
	and spondylotic foraminal stenosis was reported in 10.3/28.2%. Herniated discs were mainly

seen in the group with one compressed root, whereas spondylotic foraminal stenosis more often caused multiple root compressions, including compression of clinically unaffected roots (see table

3

).

Table 3

 3 Causes of root compression in 78 patient with suspected cervical radiculopathy who

	underwent MRI
	total n=78	radiologist 1	radiologist 2
	single root compression	n=37 (47.4)	n=37 (47.4)
	-by herniated disc	25 (32.1)	24 (30.8)
	-by spondylotic foraminal stenosis 10 (12.8)	4 (5.1)
	-by both	2 (2.6)	9 (11.5)
	multiple root compression	n=31 (39.7)	n=29 (37.2)
	-by herniated disc(s)	7 (9.0)	3 (3.8)
	-by spondylotic foraminal stenosis 18 (23.1)	13 (16.7)
	-by both	6 (7.7)	13 (16.7)
	Normal MRI	10 (12.8)	12 (15.4)


Kalle Åström

Luce Morin
email: luce.morin@imag.fr

Random Cross Ratios

The main result of this paper is the derivation of the probability density function of the cross ratio

when the numbers a, b, c and d are independent stochastic variables with identical rectangular distribution. The probability density function of the cross ratio is needed to study the effectiveness of invariant based recognition systems. A short background into invariant based recognition is given. The use of the result is illustrated by examples and a short comment is given on generalizing the result to other distributions on a, b, c and d. Although theoretical, this result is useful in understanding algorithms based on the cross ratio, e.g. in autonomous vehicle navigation, object recognition and reconstruction from images.

Introduction

The cross ratio of four real numbers is

k ¦ a § b § c § d© ¦ c a¦ d b¦ d a¦ c b (1)
This cross ratio is invariant under projective transformations. Geometric invariants like the cross ratio have become increasingly popular in Computer Vision, e.g. in reconstruction and recognition algorithms. One interesting application is to use the cross ratio in autonomous navigation of vehicles, cf. [START_REF] Åström | A correspondence problem in laser guided navigation[END_REF][START_REF] Åström | Automatic mapmaking[END_REF]. In this application vehicles are navigated by means of identical reflective beacons, which are placed on the walls of the factory. Since many of the beacons are collinear, the cross ratio of a collinear group of beacons is invariant with respect to the viewpoint of the vehicle. This enables fast recognition. Another use of the cross ratio is the recognition of planar features using the cross ratios of five coplanar points. Such objects are common in man-made environments. The use of invariants in recognition is briefly explained in section 2. In this short paper we address the problem of assessing the effectiveness of this approach. The main result is deriving the analytical expression for the probability density function and the cumulative distribution function of the cross ratio of four random numbers, in this case four independent numbers of identical rectangular distributions. Much to our surprise this probability function has a simple analytical expression. Some brief comments are also given on how this result can be extended to other input distributions.

Invariant based recognition

Object based recognition using viewpoint invariant features have become increasingly popular in computer vision, cf. [START_REF] Mundy | Geometric Invariance in Computer Vision[END_REF][START_REF] Mundy | Applications of Invariance in Computer Vision[END_REF]. The recognition problem can be stated as follows.

In an image, extract features ω and try to identify them with a number of known objects

Ω 1 § § Ω N .
One difficulty is that the extracted image feature ω, which is a projection of some 3D object features Ω, depends in a complicated way on viewpoint and camera calibration. It is possible to test whether a given object Ω i is compatible with image features ω. This is called matching. Matching is usually quite time consuming, but all known information about camera and viewing situation can be used. Similar to matching is the use of mutual invariants, cf. [START_REF] Heyden | Proximity measures in computer vision[END_REF]. These are functions f on a object/image pair whose values f ¦ Ω i § ω¨are zero or close to zero if the object Ω i is compatible with image features ω. Both matching and mutual invariants share a common disadvantage: A linear search through the list of objects

Ω 1 § § Ω N is needed.
Pure image invariants, on the other hand, is a way of getting around this time consuming linear search through all objects in the object data base. An image invariant is a function f on the image features only, such that f ¦ ω¨is the same for all images of one object. Using such an invariant feature, one can directly reduce the search of compatible objects to those Ω i , which have invariant value close to that of ω. Two practical questions arise.

1. What is meant by close? How can we ensure that we do not miss the correct object within some probability?

2. How much reduction is achieved?

After reducing the list of compatible objects these are examined more thoroughly using mutual invariants or matching. One of the simplest situations is the case of four collinear points. In this case an invariant, the cross ratio Eq. ( 1), was known already by Pappus 300 A.D. Even this simple case has by itself interesting applications, e.g. in autonomous vehicle navigation [START_REF] Åström | A correspondence problem in laser guided navigation[END_REF][START_REF] Åström | Automatic mapmaking[END_REF]. The analysis of the cross ratio is also needed to understand other invariants like the invariants of 5 coplanar points and the invariant of 6 non-coplanar points seen in two images. These are based on the cross ratio.

The question 1 above can be solved to reasonable degree of accuracy by standard techniques illustrated by the following example.

Example. Assume that four points ω © ¦ 10 § 123 § 223 § 310¨have been measured. Assuming that the errors have a small estimated standard deviation

σ a © σ b © σ c © σ d © 1 4, then the cross ratio is k¦ 10 § 123 § 223 § 310¨¡ 1
33 and the standard deviation σ k of the cross ratio can be approximated by

σ 2 k ¡ ¦ ∂k ∂a ¨2σ 2 a ¢ ¦ ∂k ∂b ¨2σ 2 b ¢ ¦ ∂k ∂c ¨2σ 2 c ¢ ¦ ∂k ∂d ¨2σ 2 d which gives σ k ¡ 0
014. So by reducing the matching to those objects with cross ratio within

the interval £ k 2 58 σ k § k ¢ 2 58 σ k¤ ¡ £ 1 29 § 1 36 ¤
, there is a low probability, approximately 0.01, that we miss the correct object.

The question 2 above is the question on how much reduction we can expect given such an interval. Are some cross ratios more common than others? A similar question is what happens if the image feature ω is not the projections of a known object but just four random image points. Are some cross ratios more probable than others? In the next section we try to answer this question by deriving the probability density function of four independent identically rectangularly distributed points. The reason for using rectangular distribution is based on the assumption that a random point in the image has uniform distribution over the whole screen.

Main Result

X © ¦ C A¨¦ D B¦ D A¦ C Bï s f X ¦ x© ¡ ¢ f 1 ¦ x¨¢ f 3 ¦ x¨if x £ 0 f 3 ¦ x¨¢ f 2 ¦ x¨if 0 £ x £ 1 f 2 ¦ x¨¢ f 1 ¦ x¨if 1 £ x (2)
where

f 1 ¦ x¨© 1 3 ¦ ¦ 2x 1¨ln¦ x x¤ 1 ¨ 2¨f 2 ¦ x¨© 1 3 ¦ ¦¥ x § 1¨ln ¥ x© § 2 ¥ 1¤ x¥ x¤ 1¨3 f3 ¦ x¨© 1 3 ¦ ¦¥ x¤ 2¨ln ¥ 1¤ x¤ 2x x 3 The corresponding cumulative distribution function is F X ¦ x© P¦ X £ x¨© ¡ ¢ F 1 ¦ x¨¢ F 3 ¦ x¨if x £ 0 1 3 if x © 0 1 2 ¢ F 2 ¦ x¨¢ F 3 ¦ x¨if 0 £ x £ 1 2 3 if x © 1 1 ¢ F 1 ¦ x¨¢ F 2 ¦ x¨if 1 £ x (3)
where

F 1 ¦ x¨© 1 3 ¦ x¦ 1 x¨ln¦ x¤ 1 x ¨ x ¢ 1 2 F2 ¦ x¨© 1 3 ¦ x¤ x ln ¥ x¤ 1¥ x¤ 1¨2 F3 ¦ x¨© 1 3 ¦ ¥ 1¤ x¨ln ¥ 1¤ x© § x x 2
The theorem is a consequence of the symmetry of the cross ratio with respect to permutations in A § B § C § D and the following lemma.

Lemma 3.1. If A, B, C and D are independent stochastic variables with identical rectangular distributions, then the conditional probability density function f X ordered ¦ x¨of the cross ratio, given that A

£ B £ C £ D, is f X ordered © 2¦ ¦ 2 x 1¨ln¦ x ¥ x¤ 1¨¨¤ 2¨if x 1 0 if x 1 ( 4
)
The corresponding cumulative distribution function is

F X ordered ¦ x¨© 2 ! 1 x¦ x 1¨ln ! x¤ 1 x " x " if x 1 0 if x 1 (5)
Remark. The functions

f 1 § f 2 § f 3 § F 1 § F 2 § F 3 of Theorem 3.
1 are obtained from F X ordered of Lemma 3.1 by simple transformations as described later.

Proof. (of the lemma)

Since the cross ratio is invariant under projective transformations there is no loss in generality to assume that A, B, C and D are rectangularily distributed between 0 and 1, i.e.

f A ¦ a© f B ¦ b¨© f C ¦ c© f D ¦ d© 0 if x 0 or x ¡ 1 1 if 0 £ x £ 1
We will first calculate the probability 1

F X ordered ¦ x© P¦ x £ X ¢ A £ B £ C £ D¨, using P¦ x £ X ¢ A £ B £ C £ D¨© P¦ x £ X § A £ B £ C £ DP ¦ A £ B £ C £ D¨
All of the 24 permutations are equally likely so P¦ A £ B £ C £ D¨© 1 24. Denote by P x the joint probability P¦ x £ X § A £ B £ C £ D¨, which can be calculated according to the definitions

P x © ¤£ ¥£ ¦£ ¥£ M x f ¥ A B C D¨¦ a § b § c § d¨da db dc dd
where

M x © ¦ a § b § c § d §¢ k¦ a § b § c § d¨ x § 0 £ a £ b £ c £ d £ 1
With this ordering the cross ratio is greater than one so P x © 1 24 for all x 1. In order to calculate the probability P x for other x we have to I. Determine the region M x .

II. Calculate the integral.

It turns out that the integral can be solved analytically and the solution has a remarkably simple expression.

I. Determination of

£ c £ d £ 1 the condition k 2 ¦ b § c § d¨ x is verified for b ¨¦ b min § c¨, where b min ¦ c § d § x© x 1 x c d c fulfills k 2 ¦ b min § c § d¨© x.
To recapitulate: We want to find those a, b, c and d with 0 £ a £ b £ c £ d £ 1, whose cross ratio k ¦ a § b § c § d¨is greater than x. According to the above for x 1 this set M x is

M x © ¡ ¢ ¦ a § b § c § d¨¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ 0 £ d £ 1 § 0 £ c £ d § b min ¦ c § d § x¨£ b £ c § 0 £ a £ a max ¦ b § c § d § x¨¢ £ ¤ For x 1 the region M x is the whole region ¦ a § b § c § d¨¢ 0 £ a £ b £ c £ d £ 1 .

II. Computation of the joint probability distribution. The probability P x is

P x © £ ¦£ ¥£ ¦£ M x f ¥ A B C D¨¦ a § b § c § d¨da db dc dd
Using rectangular distribution this becomes

P x © £ 1 0 £ d 0 £ c b min £ a max 0 1 da db dc dd (6)
Extensive calculations, by means of Maple, cf. Appendix A, give

P x © 1 24 ¦ 2x¦ x 1¨ln¦ x 1 x ¨¢ 2x 1¨
Finally F x ordered © 1 24 P x and f x ordered is the derivative of F x ordered with respect to

x. This gives Eqs. ( 4) and ( 5).

The functions (4) and ( 5) are the basis on which the distribution function F X and the probability density function f X is built. This part only depends on the symmetry of the cross ratio.

Proof. (of the theorem)

The cross ratio for an ordered four point configuration k¦ a § b § c § d¨is related in a simple way to the cross ratio of all permutations of the four points a § b § c § d . The 24 permutations can be divided into 6 groups of 4 each. In each such group the cross ratio of the permuted points is obtained from the cross ratio of the ordered points by one of the following functions. Assume that the stochastic variable W has probability density function f W ¦ w¨. A trans- formed stochastic variable Z, constructed as Z © g¦ W ¨, then has probability density function

g 1 : ¦ 1 § ∞¨ ¦ 1 § ∞¨g 1 ¦ x¨© x g 2 : ¦ 1 § ∞¨ ¦ 1 § ∞¨g 2 ¦ x¨© 1 ¢ 1 ¦ x 1g 3 : ¦ 1 § ∞¨ ¦ 0 § 1¨g 3 ¦ x¨© 1 x g 4 : ¦ 1 § ∞¨ ¦ 0 § 1¨g 4 ¦ x¨© 1 1 x g 5 : ¦ 1 § ∞¨ ¦ ∞ § 0¨g 5 ¦ x¨© 1 x g 6 : ¦ 1 § ∞¨ ¦ ∞ § 0¨g 6 ¦ x¨© 1 ¦ x 1¨
f Z ¦ z© f W ¦ g ¤ 1 ¦ z¨ §¢ d ¦ g ¤ 1 dz ¦ z¨¢ (7)
Using the transformation rule on f X ordered for g 1 § § g 6 gives

f X ¦ x© 1 6 Σ 6 i¢ 1 f X ordered ¦ g ¤ 1 i ¦ x¨ §¢ d ¦ g ¤ 1 i dx ¦ x¨¢ (8)
which can be simplified as Eq. ( 2). Finally the cumulative distribution function (3) is simply is the indefinite integral of the probability density function [START_REF] Åström | Automatic mapmaking[END_REF].

Example. Continuation of the first example.

Assume as in the previous example that we have measured a cross ratio k ¡ 1 33 and standard deviation σ k ¡ 0 014. If we have a large database of 1000 'random' objects, we can expect

1000 ¦ F X ordered ¦ 1 36¨ F X ordered ¦ 1 29¨¨¡ 45
of them to have a cross ratio in the interval £ 1 29 § 1 36 ¤ . Thus by using the cross ratio we directly reduce the number of objects that we have to match from 1000 to 45. These figures give an indication on the effectiveness of using the cross ratio in object recognition.

Numerical Illustration

We will now show what the probability distribution looks like and at the same time we will show results of a Monte Carlo Simulation. In the Monte Carlo simulation 300 000 samples of rectangularly distributed A, B, C and D was created using a pseudo random number generator. The cross ratios for these samples were computed. Both the distribution function and the probability function were estimated using this material. The result is shown together with a plot of the analytical expressions in Fig. 1.

Generalizations

In this section we will briefly discuss how the probability density function f X of X depends on the assumed input probability density function of the four points,

f A © f B © f C © f D © f input .
One immediate observation is that a projective transformation of the input distribution does not change f X since the cross ratio is invariant under projective transformations.

In this paper rectangularily distributed points have been used. This input distribution is relevant in several cases. In other cases, e.g. the cross ratio of four random directions in the plane, other distributions would be more appropriate.

The main results of this paper can easily be extended to other distributions. Typical examples are polynomial distributions with compact support. The calculations and the analytical expressions of the resulting probability function may, however, be quite complicated. Therefore, in Appendix A below, we present a Maple program that calculates the probability function of the cross ratio given a input distribution. As an example the input distribution

f input ¦ x© 2x if 0 £ x £ 1 0 if x 0 or x ¡ 1 gives f X ordered ¦ x© ! 8¦ 2x 1¨¦ 3 x 2 3 x 1¨ln¦ x x 1¨¢ 12¦ 4x 2 4x 1¨" 5 (9)
Although this is a completely different expression than Eq. ( 4), the general appearance of the two probability functions are quite similar.

Another input distribution, the Gaussian distribution, has been studied by Steve Maybank, cf. [START_REF] Maybank | Probabilistic analysis of the application of the cross ratio to model based vision[END_REF]. The appearance of the resulting probability function f X is also quite similar to that of Eq. ( 4).

Thus for most of these simple input distributions the appearance of the probability density function of the cross ratio has similar appearance, i.e. similar to Fig. 1. They are different but they share some important features, having roughly the same shape with logarithmic singularities at the 0 § 1 and in some sense a logarithmic singularity also at ∞. Thus for a rough assessment of the effectiveness of using the cross ratio, any of these distributions can be used.

Theorem 3 . 1 .

 31 If A, B, C and D are independent stochastic variables with identical rectangular distributions, then the probability density function f X ¦ x¨of the cross ratio

¢ 2 ∞

 2 M x . To obtain the boundaries of M x we first study the boundaries for a. Assume that b, c and d are fixed and ordered. The partial derivative with respect to a is ∂k ∂a© ¦ d b¦ d c¦ d a¨2 ¦ c b¨This is negative because of the ordering. The maximal cross ration is thus obtained fora © 0, k 2 ¦ b § c § d© k¦ 0 § b § c § d¨© c ¦ d bd ¦ c bä nd the minimal cross ratio is obtained for a © b, which gives k ¦ b § b § c § d¨© 1. Now the condition k ¦ a § b § c § d¨x can only be satisfied using some a ¨¦ 0 § b¨if k 2 ¦ b § c § d¨is greater than x and a ¨¦ 0 § a max ¨, where a max solves k ¦ a max § b § c § d¨© x. Some calculations givea max ¦ b § c § d § x© cb cdTo ensure that k 2 is greater than x we fix c and d and study k 2 as a function of b ¨¦ 0 § c¨.The partial derivative of k 2 with respect to b is ∂k 2 ∂b Since 0 b £ c £ d it follows that k 2 is strictly increasing. It attains the value 1 for b © 0 and k as b c. For any x 1 and all 0

Figure 1 :

 1 Figure 1: The probability density function and the cumulative probability distribution of the cross ratio of four points with independent identical rectangular distribution. The crosses illustrate the estimated probability function from the Monte Carlo Simulation and the graph illustrate the theoretical result.

Conclusions and Acknowledgements

The main contribution of this paper is the derivation of the analytical expression for the probability function of the cross ratio of four random collinear points. It was surprising that this was possible at all and that the result was so simple. We believe that this will help understanding the effectiveness of invariant based object recognition, used for example in autonomous vehicle navigation, and surveillance systems. It seems that the general appearance of the probability function does not depend on the distribution of the four points as long as they are of independent and identical distribution. Interesting extensions of this work would be to assess the effectiveness of other invariants used in vision.

The work is done within the ESPRIT-BRA project View-point Invariant Visual Acquisition.

Appendix A

This is the short maple script for calculating the distribution function F x ordered , for a given


Nicolas Andre

Towards the Embedding of On-Line Hand-Eye Calibration into Visual Servoing

This work is related to the visual servoing of a robot hand-mounted camera. The control law of this robotic system exhibits the necessity of determining the position/orientation of a camera reference frame with respect to the control reference frame. The hand-eye calibration consists in the determination of this transformation (i.e. a rotation and a translation).

We provide in this paper an on-line calibration method in two stages: rst, an initial estimation is computed with two self-calibration movements, then the estimated hand-eye transformation is updated using the controlled motions of the robot. As the latter generate low amplitude rotations, the classical formulations are not e cient any more. The homogeneous matrix equation AX=XB appearing in usual approaches is therefore reformulated into a linear system and used as a basis of a Kalman lter which is to be run in parallel with the visual servoing.

Preliminary simulation results were obtained and show the good behaviour of the method.

Introduction

To control a robot by visual servoing, one has to know or to estimate several parameters ruling the control laws. In the case of a visual system rigidly mounted onto a robot arm, some of these parameters are related to the robot itself (robot calibration), some to the visual system (camera calibration) and some to the spatial relation between the robot and the visual system (hand-eye calibration).

More speci cally, hand-eye calibration consists in determining the rigid transformation (i.e. a rotation and a translation) between a robot reference frame (usually, the end-e ector frame) and the visual system reference frame (usually, in monocular systems, the camera frame).

The current approaches 1, 2, 6, 7, 9, 10, 11, 12] operate within a standard o -line calibration scheme. However, we would like to free ourselves from the human intervention required by such approaches. The goal of this paper is therefore to propose an on-line calibration method, i.e. without a-priori knowledge or any calibration object (self-calibration), and during the accomplishment of the robotic task.

From the viewpoint of a self-calibration procedure, the current solutions su er from the drawback to need large motions with signi cant di erences between each other, i.e. the rotation axis associated to each displacement should be di erent from the rotation axes of the other movements. Moreover, the higher accuracy in the estimation of the hand-eye transformation is demanded, the more motions are required ! The method we propose consists of a self-calibration stage where only two prede ned movements are needed followed by an on-line updating stage using the controlled robot motions.

The previous approaches use either the axis/angle or the quaternion representations of a rotation to solve the well-known AX = XB homogeneous matrix equation, and thus, need to avoid robot and camera motions with low amplitude rotation. (An exception is 9] where a method is presented which avoids the AX = XB formulation, but requires a non-linear minimization performed with the Levenberg-Marquardt algorithm.) However, controlled robot motions will be of low amplitude and will only contain small rotations. In the axis/angle representation, this means that the robot motions will have a rotation angle close to 0 , where the rotation axis can not be de ned. It is therefore interesting to nd a new formulation of the homogeneous matrix equation which will lead to the building of a Kalman lter.

The remainder of this paper is organized as follows: section 2 is devoted to the problem formulation, section 3 contains a complete set of simulation results and section 4 is the conclusion.

Problem formulation

Visual servoing and hand-eye calibration

Recall the simplest visual servoing control law under the task-function approach with exponential conver- where Tc is the kinematic screw applied to the robot (expressed in the end-e ector frame), is a gain factor, ẽ is the task function (to be regulated to 0), L T is the interaction matrix (expressed in the end-e ector frame) and s (resp. s ) is the current (resp. desired) image primitive.

The interaction matrix L T is a function of the camera intrinsic parameters (namely, the scale factors u and v , and the angle between the image axes), a canonical interaction matrix expressed in the camera frame (L T c , see 4] for details) and the hand-eye transformation (X = (R X ; tX )):

L T = u u cotg 0 v sin L T c R X R X A s ( R X T tX ) 0 R X
where A s ( ) denotes the skew matrix.

The intrinsic parameters are allowed to be only roughly determined since they do not have a strong inuence on the stability of the system. However, handeye transformation has to be accurately estimated 3]. A way to obtain it is to consider a pair of robot and camera positions (see Figure 1).

As the robot performs the rigid motion B = (R B ; tB ), the camera undergoes a displacement de- termined by the rigid transformation A = (R A ; tA ). A and B are related by: AX = XB

(1) The latter equation is an homogeneous matrix equation and, thus, can be rewritten as one equation dealing with the rotational part and another with the

translational part: R X R A R X R B T = 0 (2) R X tB + tX R A tX = tA (3)
We assume for the rest of this paper that the robot motion is measurable and that the camera motion can be estimated.

Since controlled motions often contain small rotations, the purpose of the method proposed here is to deal with such ill-de ned cases. To achieve this goal, the homogeneous matrix equation (1) can be reformulated as a linear system of the form: Cx = ỹ where x is a 12-vector containing the 9 coe cients of the rotation matrix R X and the 3 coe cients of the translation tX .

New formulation

De ne ~i, ~j and k, the row vectors of the rotation R X , R X = 0 B @ ~iT ~jT kT 1 C A and x = ~iT ~jT kT tX T T .

With these notations, equations (2) and ( 3) are rewritten as:

I 9 M 0 9 3 N I 3 R A | {z } C x = 09 1 tA | {z } ỹ (4)
where I n is the (n n) identity matrix and the notation R A(p;q) means the coe cient located on the p th row and the q th column of matrix R A .

M = 0 @ R A(1;1) R B R A(1;2) R B R A(1;3) R B R A(2;1) R B R A(2;2) R B R A(2;3) R B R A(3;1) R B R A(3;2) R B R A(3;3) R B 1 A ; N = 0 B @
The previous expression is derived by expanding the matrix equations (2) and (3). This exhibits linear relations between the coe cients of vector x. Then, these relations just have to be collected into the linear system above.

This formulation does not ensure that the estimated row vectors ĩ, j and k will ful ll the orthogonality constraint on the rotation matrix R X . However, it has the advantage of o ering a linear expression of the problem, thus allowing the use of a Kalman lter. As for the orthogonality constraint, it will be discussed in section 2.5.

Kalman lter

A Kalman lter is built upon this model and the assumption that the rigid transformation X remains constant as:

xk+1 = xk + ṽ(k) ỹ(k) = C k xk + ñ(k) (5)
where the subscript k denotes the value after the k th robot/camera motion, ṽ represents the noise induced by modeling errors and ñ represents the noise in the measurements. x is de ned in (4) and will be referred to as state vector in the following. ỹ is called the measurement vector and is also de ned in [START_REF] Espiau | E ect of Camera Calibration Errors on Visual Servoing in Robotics[END_REF]. A solution of this lter is given by:

P(0) = G 0 K(k) = P(k)C(k) T C(k)P(k)C(k) T + G 2 1 P(k + 1) = P(k) + G 1 K(k) C(k)P(k)C(k) T + G 2 K(k) T x(k + 1jk) = I 12 K(k)C(k)]x(kjk 1) + K(k)ỹ(k)
where x(k + 1jk) is the estimate at time k + 1 of the state-vector x, knowing all the required information at time k. G 0 is an initial estimate for P (which is taken as I 12 ), G 1 is the covariance matrix associated to the noise ṽ in the state equation and G 2 is the covariance matrix associated to the noise ñ in the measurement equation.

An initial estimation is needed in this method. In order to cope with the desire of keeping the number of self-calibration movements low and still have a good precision for this initial estimation, the best thing to do is to use one of the classical methods with only 2 movements.

A model of noise

In the Kalman lter noises appear that we must model.

The following assumptions will be made: ñ and ṽ are uncorrelated, robot motion is free of noise (which is not too much of an assumption since the uncertainties in the robot motion are transferred into the unknown hand-eye transformation), and, nally, the camera motion estimation process is not required to ful ll the orthogonality constraint on the rotation matrix. The latter assumption allows one to use additive noise on each coe cient of the camera rotation matrix instead of additive noise on the rotation angles, thus simplifying the analysis.

As ñ and ṽ are uncorrelated, it is straightforward to set ṽ as a Gaussian noise with a diagonal (12 12) covariance matrix G 1 . Then, we can express the noise ñ in the measurement equation. To achieve this goal, rst insert noise in the camera motion: R A = RA + a tA = ta + ṽa where the notation ~ represents the estimated or measured value and where a and ṽa are uncorrelated Gaussian noises and have respective covariance matrices 1 (9 9) and 2 (3 3). Then, feed equations (2) and ( 3) with the latter two equations and obtain:

R X RA R X R B T = a R X R B T R X tB + tX RA tX = ta + ṽa + a tX
For the purpose of analysis, assume, temporarily, that R X and tX are known. Denote ṽa 0 = ṽa + a tX and ~ a the vector form of a R X R B T : 8(i; j) 2 f1; 2; 3g 2 ; ( ~ a ) 3(i 1)+j = ( a R X R B T ) ij Hence, ñ = ~ a T ; ṽa 0 T T . The covariance matrix associated to ñ according to the noise model is:

G 2 = V T 1 V V T 1 W W T 1 V 2 + W T 1 W where V = 0 @ R X R b T 0 0 0 R X R b T 0 0 0 R X R b T 1 A and W = 0 B @ tX T 0 0 0 tX T 0 0 0 tX T 1 C A:
As R X and tX are unknown and R B is time-varying, they are approximated by: R X I 3 , tX 0 and R B I 3 . Under this approximation, the covariance matrix G 2 is chosen as: G 2 = 1 0 0 2

From Kalman state vector to rigid transformation

After each iteration (i.e. robot/camera motion), the Kalman Filter above determines a new state vector xk . This state vector contains, as seen earlier, the current estimation of the 3 row vectors ( ~i, ~j and k) of the rotation matrix R X and the current estimation of the translation vector tX .

However, it is not granted that the estimated row vectors satisfy the necessary condition of orthonormality. As the raw estimates can not be used as an approximation of a rotation, an orthogonalization has to be made. In order to keep the optimal convergence of the Kalman Filter, this operation has to be applied outside the lter loop.

The orthonormalization step consists in nding a rotation matrix R which ts the best the row vectors ~i; ~j; k. One can use the Gram-Schmidt orthogonalization method but it is not numerically robust. One should therefore use either an SVD or a QR decomposition 8] or the minimization method presented in 5], which is the most robust (see Figure 2). The underlying idea of the latter method is to nd the rotation matrix R which satis es: This problem leads to a minimization problem in the quaternion space of the form: min q (q T Bq) [START_REF] Daniilidis | The dual quaternion approach to hand-eye calibration[END_REF]

Simulation results

Now that the orthonormality constraint is guaranted, the axis/angle or quaternion formalisms can be used to quantify the simulation results. In the remaining, we will consider: relative errors in translation, computed as k tX t X k=k tX k, errors in rotation expressed in the quaternion representation and computed as kq X qX k.

Sensitivity to measurement noise

In order to study the sensitivity of the method to noise in the measurements, 50 simulations were run with several noise levels. Noise in the rotation part and noise in the translation part have been studied separately.

Robot motions are de ned by choosing Roll-Pitch-Yaw angles according to Gaussian laws of standard deviation equal to 2 =100 and translation coe cients according to Gaussian laws of standard deviation equal to 1 length unit. Hand-eye transformations to be recovered are also chosen randomly with Roll-Pitch-Yaw angles between 0 and 2 and translation coe cients have standard deviation equal to 100 length unit.

Sensitivity to noise in the rotation part. Here, the noise in the translation measurements is kept constant (i.e. a measurement error in the camera translation ṽa is randomly chosen with 0.1 length unit stan- dard deviation, which is the value expected by the Kalman lter) and 6 cases were studied. They all correspond to a di erent choice of standard deviation while randomly choosing the measurement error in the camera rotation a . The chosen values for the standard deviations are respectively: 0, 0.01, 0.05, 0.1, 0.5, 1, while the Kalman lter is tuned for a 0.1 standard deviation (i.e. the correlation matrix is chosen diagonal with non-zero coe cents equal to 0.01). (Recall that this noise is added to a rotation matrix, which therefore has unit row vectors, i.e. its coe cients vary most of the time between -1 and 1.)

Figures 3 and4 show the results of these simulations.

Sensitivity to noise in the translation part. A Gaussian noise is added to each coe cient of the camera translation. Fifty simulations were run for 6 levels of noise (i.e. standard deviation equal to 0, 0.01, 0.05, 0.1, 0.5 and 1). These values are to be compared with the camera translations. However, it is not easy since they depend on the robot motions and the hand-eye transformation. Nevertheless, we can have an idea of the amplitude of the camera translations. Recalling (3) and approximating R A by the identity matrix in it leads to tA R X tB , thus considering that camera and robot translations have almost the same amplitudes. In the present case, they have standard deviation equal to 1.

Figures 5 and6 show the results of these simulations.

Sensitivity to initial estimate

This section deals with sensitivity to initial estimates.

In order to study it, the reference hand-eye transformation was disturbed by adding white noise to its Roll-Pitch-Yaw angles and to its translation coe cients.

The standard deviation of these disturbances are dened as a percentage of the standard deviation used in the determination of the hand-eye transformation.

Here again, fty simulations were run for each of the following error percentages: 0%, 5%, 10%, 50%, 75% and 100%. The results can be found in Figures 7 and8.

Realistic simulations

The previous section shows the good numerical behaviour of the proposed method, but it still remains to study more realistic cases. Aiming to this, one should choose a realistic initial estimation of the hand-eye transformation then use it with realistic robot motions. Computing an initial estimation. As discussed above, this method can not completely free itself from the classical approaches since it needs an initial estimation of the hand-eye transformation, but it requires a smaller number of calibration movements. Here, two methods are compared: the 'historical' one 11], based on the axis/angle representation of the rotation and 2] which is based on a dual quaternion representation of the rotation. The accurate non-linear method 7] is left apart since it uses a Levenberg-Marquardt minimization, thus taking too much time for an accuracy increase out of the present needs.

Tests have been run to determine which of the two methods should be used. They were done with various noise levels and, for each noise level, 500 simulations were realized. Each simulation consisted in choosing randomly a hand-eye transformation, computing a disturbed camera motion from it and two prede ned self-calibration moves, recovering the hand-eye transformation from robot and camera motions with both methods and, nally, comparing the results (see Figure 9). Realistic robot motions. In this experiment, a hand-eye transformation is chosen randomly. Then, realistic robot motions are chosen: according to the control law, the components of the kinematic screw decrease exponentially from their initial values; there- fore, an initial kinematic screw is chosen randomly as the starting point of a simulated convergent robot motion. As these robot motions are determined, the camera motions can be computed. Noise is then added to them and the obtained values are fed with the robot motions into the lter.

The whole process is preceded by 2 self-calibration movements which determine, according to 2] (since previous tests shows its higher robustness), an initial estimate of the hand-eye transformation.

The results of this experiment are shown in Figure 10.

Conclusion

This was the rst attempt to solve for the hand-eye calibration problem in the context of a hand-mounted camera. The linear formulation which has been proposed permits an on-line improvement of an initial estimation. The latter results from the use of one of the current approaches in a self-calibration scheme where the robot/camera set-up is only given two prede ned movements.

The simulation results are promising, showing the interest of such an approach. The model of noise should however be modi ed to cope with noise on the Roll-Pitch-Yaw angles rather than noise on each coefcient of the rotation matrix. Then, it will remain to experiment the method on a real robot, rst with prede ned robot motions then in a whole visual servoing application.

Figure 1 :

 1 Figure 1: Two positions of the same hand-eye system linked by a camera displacement A associated to a robot motion B via the hand-eye transformation X.

Figure 2 :

 2 Figure 2: Error in the recovering of a rotation matrix according to several methods: SVD decomposition (dashed line), QR decomposition (dash-dotted line) and minimization (solid line).

Figure 3 :

 3 Figure 3: Evolution of the relative error in the rotation (expressed as a quaternion) with various levels of noise in the rotation: 0, 0.01, 0.05, 0.1, 0.5 and 1. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dash-dotted line.

Figure 4 :

 4 Figure4: Evolution of the relative error in the translation vector with various levels of noise in the rotation: 0, 0.01, 0.05, 0.1, 0.5 and 1. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dash-dotted line.

Figure 5 :

 5 Figure5: Evolution of the relative error in the rotation (expressed as a quaternion) with various levels of noise in the translation: 0, 0.01, 0.05, 0.1, 0.5 and 1. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dash-dotted line.

Figure 6 :

 6 Figure6: Evolution of the relative error in the translation vector with various levels of noise in the translation: 0, 0.01, 0.05, 0.1, 0.5 and 1. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dashdotted line.

Figure 7 :

 7 Figure7: Evolution of the relative error in the rotation (expressed as a quaternion) with various error percentages in the initial estimates: 0%, 5%, 10%, 50%, 75% and 100%. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dash-dotted line.

Figure 8 :Figure 9 :

 89 Figure8: Evolution of the relative error in the translation vector with various error percentages in the initial estimates: 0%, 5%, 10%, 50%, 75% and 100%. Maximal and minimal errors are plotted with solid line, mean value errors with dashed line and RMS errors with dash-dotted line.

Figure 10 :

 10 Figure 10: Simulation with realistic robot motions: error in rotation (left) and in translation (right). The horizontal lines represent the initial errors after 2 calibration movements.


Roberto Cipolla

Edmond Boyer

3D Model Acquisition from Uncalibrated Images

In this paper we address the problem of recovering 3D models from uncalibrated images of architectural scenes.

We propose a simple, geometrically intuitive method which exploits strong rigidity constraints such as parallelism and orthogonality present in indoor and outdoor architectural scenes. We show how these simple constraints can be used to calibrate the cameras and to recover the 3 4 projection matrices for each viewpoint. The projection matrices are used to recover partial 3D models of the scene and these can be used to visualise new viewpoints.

Our approach does not need any a priori information about the cameras being used. A working system has been designed and implemented to allow the user to interactively build a 3D model from a pair of uncalibrated images from arbitrary viewpoints.

Introduction

Considerable e orts have been made to recover photorealistic models of the real world. The most common geometric approach is to attempt to recover 3D models from calibrated stereo images 16] or uncalibrated extended image sequences 21, 1, 17] by triangulation and exploiting epipolar 15] and trilinear constraints [START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF][START_REF] Seitz | Toward Image-Based Scene Representation Using View Morphing[END_REF]]. An alternative approach consists of visualisation from image-based representations of a 3D scene. This is the preferred method when trying to generate an intermediate viewpoint image given two nearby viewpoints since it does not need to make explicit a 3D model of the scene [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF][START_REF] Press | Numerical Recipes in C, The Art of Scienti c Computing[END_REF][START_REF] Faugeras | Representing Three-Dimensional Data as a Collection of Images and Fundamental Matrices for Image Synthesis[END_REF]12,[START_REF] Faugeras | Three-Dimensional Computer Vision: A Geometric Viewpoint[END_REF].

Facade 4] { one of the most successful systems for modelling and rendering architectural buildings from photographs { consists of a hybrid geometric and image-based approach. Unfortunately it involves considerable time and effort from the user in decomposing the scene into prismatic blocks, followed by the estimation of the pose of these primitives.

In this paper we propose a much simpler approach to construct a 3D model and generate new viewpoint images by exploiting strong constraints present in the scenes to be modelled. In the context of architectural environments, the constraints which can be used are parallelism and orthogonality. These constraints lead to very simple and geometrically intuitive methods to calibrate the cameras and to recover Euclidean models of the scene from only two images from arbitrary positions.

2 Outline of the algorithm A 3D model can be recovered from two or more uncalibrated images. Our algorithm consists of the four following stages (see gure 1).

1. We rst de ne a set of primitives { segments and cuboids { for which parallelism and orthogonality constraints are derived. These primitives are precisely localised in the image using the image gradient information. 2. The next step concerns the camera calibration: the intrinsic parameters of the camera are determined for each image. This is done by determining the vanishing points associated with parallel lines in the world. Three mutually orthogonal directions are exploited 2].

The motion (a rotation and a translation)

between the viewpoints is then computed. The motion combined with the knowledge of the camera parameters allows the recovery of the perspective projection matrices for each viewpoint. 4. The last step consists in using these projection matrices to nd more correspondences between the images and then to compute 3D textured triangles that represent a model of the scene. (1)

The projection matrix has 11 degrees of freedom and can be decomposed into the orientation and position of the camera relative to a the world co-ordinate system (a 3 3 rotation matrix R and a 3 1 translation vector T:

P = C R T (2)
and a 3 3 camera calibration matrix, C, corresponding to the following transformation:

C = 2 4 u s u 0 0 v v 0 0 0 1 3 5 ;
(

where u , v are scale factors; s is a skew parameter; and u 0 , v 0 are the pixel coordinates of the principal point (the intersection of the optical axis with the image plane 5]).

In general, the parameters of the projection matrix can be computed from the measured image positions of at least six known points in the world. This matrix can then be decomposed by QR decomposition into camera and rotation matrices 5].

In our approach the vanishing points corresponding to three mutually orthogonal directions can be used to determine:

1. the camera calibration matrix, C under the assumption of zero skew and known aspect ratio 2].

2. the rotation matrix R. 3. the direction of translation, T.

We show that the 8 degrees of freedom of the projection matrix for this special case can be determined from three vanishing points corresponding to the projections of 3 points at in nity and a reference point. The projection matrix can thus be recovered from the projection of at least one arbitrary cuboid. Applying the algorithm to two views allows the Euclidean reconstruction of all visible points up to an arbitrary scale.

Primitives

The rst step of the method consists in de ning image segments and cuboids by means of which the parallelism and orthogonality constraints can be derived. These constraints are then used to recover the intrinsic and extrinsic cameras parameters and hence the projection matrices for each view.

First, interest points are extracted in the reference views. This can be done by means of a Harris corner detector 8] and also interactively by hand since results from a corner detector may not be su cient (see gure 3). Such points are then used to de ne the image projections of primitives: segments and cuboids, which are present in the observed scene.

The de nition of a segment projection in an image is straightforward and requires two image points. The image of a cuboid is completely speci ed by six or more vertices. The interest points, extracted corners and user de ned corners which are used to de ne segments and cuboids, are not necessarily well localised in the image plane. An optimisation is therefore performed to improve their localisations. This optimisation takes into account image gradient information and modi es image point coordinates in order to position primitive edges along image edges. A vanishing point corresponds to the projection of the intersection of parallel lines at innity. A number of approaches have been proposed to localise precisely this intersection, from the simple calculation of a weighted mean of pairwise intersections 2] to more elaborate approaches involving noise assumption and nonlinear criteria 14] (see also 11] for line intersection estimation). However, all these approaches are based on criteria which take into account image distances. Though trying to minimise the distance from the vanishing point to the image lines is geometrically correct, it appears to be numerically unstable in the presence of noise. In contrast, our approach is based on a linear criterion which optimises a three dimensional direction, the dual of a vanishing point in 3D space.

Let l 1 ; ::; l n be a set of image lines corresponding to parallel space lines and let be the vanishing point de ned by these lines. Then if we suppose that l i and are expressed in homogeneous form in the image coordinate system (u; v; w), we have the following relation: l 1 ; ::; l n ] T = 0:

(4)

and the null vector of the 3 n matrix l 1 ; ::; l n ] gives a solution ^ for the vanishing point position. However, the minimised criterion does not necessarily have a physical meaning.

By choosing an appropriate normalisation for (l 1 ; ::; l n ), the solution ^ can minimise the sum of the Euclidean distances from ^ to l 1 ; ::; l n ].

As said before, this solution will not be completely satisfactory in the presence of noise due to poor conditioning of the matrix l 1 ; ::; l n ]. Now suppose that the intrinsic parameters of the camera are known. Then image lines and points can be expressed using normalised image coordinates, that is to say in the coordinate system (x; y) associated with the camera retinal plane (see gure 4). Let (L 1 ; ::; L n ) and V be the homogeneous representations of (l 1 ; ::; l n ), in the retinal plane. Since V belongs to all the image lines, we still have the relation: L 1 ; ::; L n ] T V = 0:

l l i j O L i V u v ν j V V P i P j L x y z
(5)

However, this relation has also an interpretation in the three dimensional coordinate system (x; y; z) associated with the camera (see gure 4). Indeed, L i is a vector orthogonal to the plane P i spanned by the camera projection centre O and the image line l i , and V is a three dimensional vector. This gives a physical meaning to the null vector V of L 1 ; ::; L n ] T which is the space direction closest to all planes (P 1 ; ::; P n ).

Finally, we compute the space direction V as the null vector of the matrix L 1 ; ::; L n ] T , where jL i j = 1. This is done using a singular value decomposition of L 1 ; ::; L n ] T 18]. Experiments show that estimating a space direction is more robust to noise perturbations than estimating an image point. We use this method in the calibration routine described in the next section.

In the case where the camera's intrinsic parameters are not known, an intermediate solution consists of using ( 4) with prenormalized image point coordinates 9] in order to improve the condition number of the matrix l 1 ; ::; l n ].

Camera Calibration

Having found the vanishing points we now show how they are used to derive constraints on the intrinsic parameters and the latter estimated in practice.

Constraints

A set of parallel line projections de ne a vanishing point in the image plane which is the projection of the intersection parallel lines at innity. Let V be the space direction of these parallel lines (a unit vector). The vanishing point depends only on the direction of the line in space and not its position. The relationship between vanishing point in homogeneous coordinates, and the corresponding space direction in the camera coordinate system is given by: = CV [START_REF] Faugeras | Three-Dimensional Computer Vision: A Geometric Viewpoint[END_REF] Now suppose that we have di erent vanishing points i ; j corresponding to distinct space directions V i ; V j separated by a known angle . Then:

V T i V j = cos and from (6): T i C ?T j T i C ?T j C ?1 j jC ?1 j j = cos [START_REF] Faugeras | Representing Three-Dimensional Data as a Collection of Images and Fundamental Matrices for Image Synthesis[END_REF] which provides a constraint on the matrix C ?T C ?1 and thus on the intrinsic parameters.

A particularly interesting case is when the space directions are orthogonal 23]. Then:

T i C ?T C ?1 j = 0 (8)
which gives a linear constraint on the parameters of C ?T C ?1 . It is easy to show both algebraically and geometrically that the following degenerate cases exist:

1. if both directions V i , V j are parallel to the image plane then their projections give information on the angle between image plane axes u and v only 1 .

2. if one of the directions V i , V j is orthogonal to the image plane, then their projections determin u 0 and v 0 but do not constrain the other parameters.

The faces of a cuboid de ne three orthogonal space directions, thus its projection leads to three constraints on the intrinsic parameters allowing the computation of three of them. As noticed before, degenerate cases exist. In particular, information on the parameters u and v can not be derived from a cuboid which has a face parallel to the image plance as shown in gure 5.

O Figure 5: The degenerate case when the cuboid's projection does not provide any constraint on the intrinsic parameters u and v . One of the cuboid directions is perpendicular to the image plane.

Recovering intrinsic parameters

Now suppose that the skew parameter s is zero, then, for two vanishing points i , j associated with two orthogonal space directions, (8) be-comes: and if (u i ; v i ; w i ) are the homogeneous image coordinates of i then:

(u i ? w i u 0 )(u j ? w j u 0 )

2 u + (v i ? w i v 0 )(v j ? w j v 0 ) 2 v + w i w j = 0: (10)
The expression above is linear in 1 2 u , 1 2 v and therefore leads to linear estimation of these parameters if (u 0 ; v 0 ) are known. Degenerate cases appears when space directions associated with vanishing points belong to the image plane.

A cuboid de nes three orthogonal directions, thus its projection allows, in general, the computation of three of the intrinsic parameters. In our current implementation, we determine the scale parameters from orthogonal directions using this method and under the assumption that the principal point, is known (assumed to be in the image centre). The vanishing points can, of course, be used to nd the principal point 2] but the following method was found to be more reliable.

Algorithm

The global procedure to evaluate the scale parameters is as follows:

1. Extract lines fl i g of image primitives (cuboid and segment projections).

2. Set (u 0 ; v 0 ) to the image centre coordinates.

Set ( u ; v ) initially to (1; 1).

3. Estimate the positions of vanishing points using intrinsic parameters and fl i g as explained in section 5.

4. Estimate scale parameters ( u ; v ) using [START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF].

5. Iterate steps 3 and 4 until the estimated values for the scale parameters converge.

Evaluation

We have tested this algorithm on synthetic data in order to evaluate the sensitivity of the estimated parameters to errors in the assumed position of the principal point (u 0 ; v 0 ). Figure 6 shows an example test image. Di erent images were generated with principal points at di erent distances from the assumed image centre (256; 256). The graph of gure 6 shows relative errors on the estimation of u and v . The maximum relative error is of 12:5% when the principal point is 85 pixels from the image centre. Large errors in the principal point can therefore be tolerated. Note also the di erence between relative errors on u and v . This is due to the fact that the overall error on these parameters is divided according to the relative proximity of vanishing points to one or the other image axes (i.e. according to the orientation of the cuboid's space directions in the camera coordi-nate frame).

Projection Matrices

In the previous sections we have explained how to estimate intrinsic parameters using the vanishing points. Once the intrinsic parameters are known, it is possible to estimate the relative position and orientation of the cameras, that is the rigid motion between views. This is used to compute consistent projection matrices for all views.

Two di erent approaches can be used. The conventional approach is to use many pairs of matched points to estimate a fundamental matrix and then to decompose this matrix into a rotation and translation once the intrinsic parameters are known [START_REF] Liebowitz | Metric Recti cation for Perspective Images of Planes[END_REF][START_REF] Harris | A combined corner and edge[END_REF]. Note that this method can not be used if the points lie on a cube { a critical surface 15].

An alternative method presented below is to use the vanishing points directly to estimate the rotation (since these are independent of camera translation) and then the translation from pairs of corresponding points using the epipolar constraint.

Rotation

Let 1 and 2 be two vanishing points in image 1 and 2 respectively, corresponding to the same space direction. Let V 1 = C ?1 1 1 =kC ?1 1 1 k and V 2 = C ?1 2 2 =kC ?1 2 2 k be the 3D unit vectors representing this space direction in the camera coordinate frames 1 and 2 respectively. Since vanishing points correspond to a space direction, they do not change their coordinates under viewer translation 2] but only by rotation of the viewpoint. From this observation, we derive:

V 2 = RV 1 and for n vanishing points: (V 1 2 ; ::; V n 2 ) = R(V 1 1 ; ::; V n 1 ): The rotation matrix can be estimated from at least three vanishing points in both images by SVD decomposition 11].

Translation

Let (x 1 ; x 2 ) be the projections of a space point X in the two images. By the epipolar constraints the vectors RC ?1 1 x 1 , C ?1 which gives a linear constraint on the translation parameters. Let i be the unit vector orthogonal to the epipolar plane associated with matched point (x i 1 ; x i 2 ): i = (RC ?1

1 x i 1 ) ^(C ?1 2 x i 2 )
k(RC ?1

1 x i 1 ) ^(C ?1 2 x i
2 )k and let ? be the matrix de ned by: ? = ( 1 ; ::; n ) T for n pairs of matched points. The eigenvector of the matrix ?? T corresponding to the smallest eigenvalue gives a solution for the direction of translation { the direction which is closest to all the epipolar plane in the least-squares sense. The translation is estimated up to a scalar and the ambiguity on its sign is resolved during the reconstruction stage.

3D reconstruction

From the previous step, we have computed motions between the di erent views, up to scale factors. These motions combined with the intrinsic parameters allow us to compute projection matrices for the di erent views involved. From these projection matrices, we can determine the epipolar geometry to help nd more point correspondences and then the 3D point positions. These points are then used in association with an image point triangulation to obtain 3D structure. This structure is rendered afterwards using a texture mapping procedure and the nal model is stored in standard VRML format. In this section, we summarise the di erent steps involve to produce such a model.

Reconstruction

Let P 1 and P 2 be the projection matrices for view 1 and 2 respectively. If we now take the coordinate system of the rst camera as the 3D space coordinate system then: P 1 = C 1 I 0 [START_REF] Hartley | Lines and points in three views and the trifocal tensor[END_REF] P 2 = C 2 R T (12) where R and T now represent the rigid motion between view 1 to view 2. There is an ambiguity in the sign of T. To solve this ambiguity, 3D points must be reconstructed and the sign of the translation is chosen to make sure that the depths are not negative. The estimation of the 3D coordinates of a point X is straightforward and can be obtained by the least-squares solution from at least two image points.

Since no metric information about the scene are available, reconstruction is performed up to a scale factor.

Triangulation and texture mapping

The projection matrices are used to reconstruct points individually. New primitives can be added to the original set of primitives, either automatically 1] or interactively. To compute a 3D structure, image points are then triangulated in one view. In order to preserve the information given by primitives segments, the triangulation performed is a constrained triangulation. More precisely, it is a constrained Delaunay triangulation 5] which preserve edges. The 2D triangulation leads to triangular 3D structure that represent the scene.

To render the reconstructed structure, we use a texture mapping algorithm in which a texture map is constructed for each triangular patch in the 3D model. Each element in the map is obtained by projecting the corresponding triangle point onto the image plane and by determining the intensity value by bilinear interpolation.

Experimental results

Experiments on synthetic and real data have been conducted. Figures 7 { 9 show some preliminary results for real images of the Downing College library in the University of Cambridge. These images were obtained by an Olympus digital camera. The calibration for these images was performed using two cuboids. Figure 8 shows these cuboids and the corresponding directions which determine the vanishing points. The geometry was computed using the vanishing points the 3D co-ordinates were then recovered. Figure 9 shows an example of the nal textured model. This model consists in fty textured triangles. Note that the reconstructed structure is visually consistent with the original one even if points are treated individually: reconstructed squares are squares and reconstructed right angles are right angles. This shows that the computed geometry is correct with respect to the imposed constraints. Furthermore, the realism of the 3D model is strongly dependent on the fact the constraints of parallelism and orthogonality are actually veri ed.

These results are extremely promising in scene modelling applications. [START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF]

Conclusions

The techniques presented have been successfully used to interactively build models of architectural scenes from pairs of uncalibrated photographs. The simple but powerful constraints of parallelism and orthogonality in architectural scenes can be used to recover very precise projection matrices with only a few point and line correspondences.

We plan to use these initial estimates of the projection matrices (and hence the epipolar geometry) to automatically match additional features and then to optimise the parameters of the camera motion and 3D structure by standard ray-bundle adjustment. Preliminary results are extremely promising.

3 .

 3 Figure 1: Outline of the algorithm.

Figure 2 :

 2 Figure 2: A cuboid and its projection which is de ned by ve image points and an image direction. The cuboid image edges de ne three vanishing points in the image plane.

Figure 3 :

 3 Figure 3: Example of primitive de nitions and localisations. Two cuboids and one segment have been interactively de ned using extracted corners.

Figure 4 :

 4 Figure 4: Vanishing point determination. The idea is to nd the viewing direction V of the vanishing point . This direction should belong to all the planes P i .

Figure 6 :

 6 Figure 6: The 512 512 test image and the graph of the relative error in the estimation of u and v as a function of the error in the principal point position.

Figure 7 :

 7 Figure 7: The original images from uncalibrated cameras.

Figure 8 :

 8 Figure 8: The two cuboids used for the calibration and the lines used to estimate the vanishing points.

Figure 9 :

 9 Figure 9: The 3D textured model.

Finding Vanishing PointsOnce an initial set of primitives has been dened, parallelism and orthogonality constraints are used through the vanishing points to constrain the intrinsic parameters.

is linked to the skew parameter by the relation s = ? cot = v 5]

x 2 and the direction of translation between the two views, T, are coplanar. Thus:(RC ?1 1 x 1 ) ^(C ?1 2 x 2 )] T = 0;

Acknowlegments

We acknowledge the nancial support of Ofce Workstation Limited (Edinburgh) and Panasonic, Japan.