Aline Baudry
email: aline.baudry@irisa.fr

Sylvain Gué Gan

Marie Babel

Taking caster wheel behavior into account in the kinematics of powered wheelchairs

Keywords: smart wheelchair, caster wheels, kinematics, identification, driving

The driving experience of an electric powered wheelchair (EPW) can be disturbed by unpleasant dynamic and kinematic effects of the caster wheels, particularly during maneuvers in narrow rooms and direction changes. In order to minimize this nasty behavior, we propose in this article a kinematic model of the wheelchair taking into account the effects of the caster wheels. The orientation of the caster wheels has been measured for different configurations: initial orientation, wheelchair velocity and user mass. The repeatability of the motions has been studied, and from these tests, their behavior has been modelled. This model has been then used to determine the wheelchair trajectory by using a more realistic kinematic. We have thus obtained a model allowing us to calculate more precisely the spaces that can be reach by powered wheelchairs, which are useful information to enhance control laws of their driving assistances.

INTRODUCTION

For the personal and cognitive development, it is important to explore the environment and to have social activities. However, this can be very difficult for people with reduced mobility. The use of a powered wheelchair can help them to be independent when they are unable to propel a manual wheelchair and when they need a third party for the most common tasks [START_REF] Cooper | Wheelchair selection and configuration[END_REF]. However, some of them cannot obtain the permission to drive an electric wheelchair because of visual and/or cognitive impairments. It has been shown that a large proportion of wheelchair users would benefit from using a smart wheelchair [START_REF] Simpson | How many people would benefit from a smart wheelchair?[END_REF]. For that reason, since the early eighties, research works have been conducted to design smart wheelchairs and to develop autonomous and collision avoidance assisted wheelchairs [START_REF] Simpson | Smart wheelchairs: A literature review[END_REF][START_REF] Ding | Electric powered wheelchair[END_REF].

Wheelchairs are usually modeled as a (2, 0) mobile robot [START_REF] Campion | Wheeled robots[END_REF]. The non-holonomic control problem induced from this model [START_REF] Morin | Motion control of wheeled mobile robots[END_REF] can be avoided with sensor-based approaches [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF], but the caster wheel kinematics is generally not taken into account. Yet the kinematic and dynamic effects of these wheels can be disturbing while driving. In particular, they can cause deviations from the desired trajectory during maneuvers and abrupt directions changes [START_REF] Ding | Analysis of driving backward in an electric-powered wheelchair[END_REF][START_REF] Guo | Influence of wheelchair front caster wheel on reverse directional stability[END_REF].

Hence, to assist the user to quickly and comfortably reach a location, it is important to take into account the caster wheel effect in the wheelchair dynamic behavior. These passively orientable wheels are designed to follow the wheelchair motion along non-straight paths. Their kinematic behavior has been studied in [START_REF] Chenier | An orientation estimator for the wheelchair's caster wheels[END_REF] for manual wheelchairs for only smooth motions, when they follow the wheelchair motion. A dynamic model estimator has been then implemented on an ergometer to simulate straight and curvilinear paths of a manual wheelchair [START_REF] Ché Nier | A new dynamic model of the manual wheelchair for straight and curvilinear propulsion[END_REF]. It has also been exploited in the development of assistances applied to manual wheelchairs [START_REF] Xi | Force sensorless power assist control for wheelchair on flat road using recursive least square with multiple forgetting[END_REF] but it has never been used on powered wheelchairs.

In addition, caster wheel dynamical effects induced by direction changes have not been modelled, whereas collisions mostly occurred in such situations [START_REF] Gaal | Wheelchair rider injuries : Causes and consequences for wheelchair design and selection[END_REF]. Hence, designing a maneuver assistive control law requires the definition of a related model. The dynamic behavior of the caster wheels depends on their frictions with the ground, and on the rolling resistance which oppose to their self-orientation motion [START_REF] Johnson | Dynamic modeling of an electric wheelchair[END_REF]. If the Pacejka's magic formula [START_REF] Pacejka | The magic formula tyre model[END_REF] provides a complete model of any wheel behavior, it requires a good knowledge about the nature of the ground. Nevertheless, to implement an assistive control-law taking into account the caster wheel behavior, it is not feasible to consider the multiple grounds encountered during the navigation because of the complexity to characterize their properties in real-time. In [START_REF] Lee | Orientation compensation for initially misaligned caster wheels[END_REF], a compensation scheme based on odometry measurement prevents the wheelchair to drift because of misaligned caster wheels. This assistance is mostly based on odometry and, for a low-cost and pragmatic assistance solution, we aim to develop a sensorless control law. For such a scheme, a behavior model is required.

This paper is organized as follows. First, section II presents the wheelchair kinematics and the caster wheels kinematics. The caster wheels' behavior during direction changes is defined in section III and applied in section IV to compute the wheelchair kinematics. Finally, section V discusses the prospects of a control law which could be used for maneuvering assistance of powered wheelchairs.

KINEMATIC MODELS OF THE WHEELCHAIR AND THE CASTER WHEELS

Powered wheelchairs are modeled as a (2,0) mobile non holonomic robot [START_REF] Campion | Wheeled robots[END_REF]. The kinematics resulting from the differential motorization of the wheelchair allows the driver to move along the longitudinal direction (𝑉) and to rotate (Ω) around z-axis. Thus two degrees of mobility make it possible to control the position and the orientation of the wheelchair.

Wheelchairs are equipped with caster wheels (from 2 up to 4, depending on the configuration) that clearly influence the trajectory. Indeed, contrary to fixed wheels (Fig. 1), passive caster wheels or off-centered orientable wheels mounted onto wheelchairs present an offset 𝑑 (Fig. 2) which allows their self-orientation. If these wheels present the advantage of following non-linear trajectories, during abrupt changes of directions (a backward motion after a forward motion for example), they are opposed to the movement of the wheelchair before self-orientation. Then, the desired kinematics of the wheelchair is not respected, and trajectory deviations as well as brutal accelerations are observed. This transition motion is followed by a smooth trajectory, where caster wheels are aligned with the trajectory induced by differential wheels. We aim then at modelling this transition state that corresponds to a sudden reorientation of a caster wheel.

Caster wheel velocity -steady state case

Let 𝐹 𝑊 (𝐹, 𝑋, 𝑌, 𝑍) be the frame attached to the wheelchair. Let 𝐹 𝐴 𝑖 (𝐴 𝑖 , 𝑥 i, 𝑦 i, 𝑧 i, ) and 𝐹 𝐶 𝑖 (𝐶 𝑖 , 𝑥 i, 𝑦 i, 𝑧 i, ) be respectively the frames located on the steering axis 𝐴 𝑖 of the caster wheel 𝑖 and the related ground contact point 𝐶 𝑖 (Fig. 3). Let 𝑉 and Ω be the linear and angular velocities of the wheelchair. The velocity of the caster wheels is obtained at the ground contact point C with respect to frame, and depends on the velocity of the wheelchair itself. The steady state longitudinal velocity 𝑉 ⃗ ∞ the caster wheel is given by

𝑉 ⃗ ∞ = 𝑉 ⃗ + 𝐴 𝑖 𝐹 ⃗⃗⃗⃗⃗⃗ × Ω ⃗⃗ (3)
By denoting (𝑃 𝑋 , 𝑃 𝑌 , 0) the position of the steering axis of the caster wheel, we obtain

𝑉 ∞ ⃗⃗⃗⃗ = ( 𝑉 -𝑃 𝑌 Ω 𝑃 𝑋 Ω 0 ). (4)
From this longitudinal velocity observed in steady state motion, we can calculate the angle 𝛼 ∞ = atan2(𝑃 𝑥 Ω, 𝑉 -𝑃 𝑦 Ω ) that the caster wheel should respect (Fig. 4).

Caster wheel velocity -transient case

During the transient state, the wheel velocity 𝑉 ⃗ 𝑖 is decomposed into a longitudinal velocity and a lateral velocity due to the angle (𝛼 ∞ -𝛼(𝑡)) (Fig. 4). In the frame 𝑅 𝐶 𝑖 , they are defined as

𝑉 ⃗ long / 𝑅 𝐶 = ‖𝑉 ⃗ ∞ ‖ × cos(𝛼 ∞ -𝛼(𝑡)) • 𝑥 (6) 𝑉 ⃗ lat / 𝑅 𝐶 = ‖𝑉 ⃗ ∞ ‖ × sin(𝛼 ∞ -𝛼(𝑡)) • 𝑦 . ( 7
)
In the frame of the wheelchair 𝐹 𝑊 , we then calculate the longitudinal velocity as

𝑉 ⃗ long / 𝑅 𝐹 = ( ‖𝑉 ⃗ ∞ ‖ × cos(𝛼 ∞ -𝛼(𝑡)) cos(𝛼(𝑡)) ‖𝑉 ⃗ ∞ ‖ × cos(𝛼 ∞ -𝛼(𝑡)) sin(𝛼(𝑡)) 0 ). ( 8
)

Wheelchair kinematics along with caster wheel effects

The kinematics of the wheelchair (𝑉(𝑡), Ω(𝑡)) can be determined from the combinations of the longitudinal velocity of each caster wheel, expressed in the wheelchair frame, so as

𝑉(𝑡) = 1 𝑛 ∑ 𝑉 ⃗ long 𝑖 / 𝑅 𝐹 𝑛 𝑖=1 (𝑡) • 𝑥 , (9)
Ω(𝑡) = 1 𝑛 ∑ 𝑉 ⃗ long 𝑖 / 𝑅 𝐹 𝑛 𝑖=1 (𝑡) × 𝐴 i 𝐹 ⃗⃗⃗⃗⃗⃗ cos (𝛽 - 𝜋 2 ) ( 10
)
where 𝛽 is the angle between 𝑉 ⃗ long 𝑖 / 𝑅 𝐹 (𝑡) 𝛽 (Fig.

CASTER WHEEL BEHAVIOR DURING DIRECTION CHANGES

From the previous section, it has to be noted that the lateral velocity of the wheel will create a rotational movement because of the offset 𝑑, thus inducing a reorientation of the wheel. A kinematic study is thus not enough at this level since we must take into account the efforts that apply at the wheel. However, as we do not want to take into account the contact properties with the ground because of the complexity of the related identification, we propose to realize a kinematic study of the behavior. We conducted experimentations to obtain a transfer function realistic enough to use it in a control law.

Data acquisition setup

The characteristics of the sensors used for the instrumentation of the wheelchair (Fig. 6) are:

• Caster wheel orientation: Two incremental encoders with 1000 pulses per revolution, for a 0.36 degrees resolution, have been mounted on rotation axis of the caster wheels. • Driving wheel angular velocity: Two incremental encoders with 12 pulses per revolution (1.23 degrees resolution) have been mounted on the motor axis. • Wheelchair trajectory: Laser tracker Radian from API with a resolution of 0:08 µm.

The data acquisition was realized by means of a Raspberry Pi 3 synchronized with the laser recording at the sampling rate f = 100 Hz. The control of the wheelchair was made thanks to a Tiva control card (Texas Instruments) connected to the CAN bus. During these experiments, we recorded the orientation of the caster wheels for 72 configurations of (𝛼 0 , 𝑉, 𝑀 𝑅 ) while operating a forward-backward movement. For each configuration, 5 data sets were recorded.

Repeatability

Before modeling the system, it is important to evaluate the repeatability of the observed behavior. This repeatability is verified if the correlation of the measurements for a given configuration is equal to 1. Figure 8 highlights in particular the non-repeatability of the behavior for the 0° initial orientation. It outlines that the caster wheels do not behave the same for the 5 successive experiments, which present different delays and rotation directions. The analysis of the trajectories leads to the same conclusion. As our measurements are repeatable for the 30°, 60° and 90° initial orientations, we can now model the caster wheel behavior.

A transfer function to model the behavior

We choose to model the caster wheels self-orientation with a transfer function to, afterwards, implement a dedicated control law. In a first step, we identify the 1 st , 2 nd and 3 rd order transfer functions corresponding to our measurements. The related transfer functions 𝐻 1 𝑠𝑡 (𝑠), 𝐻 2 𝑛𝑑 (𝑠) and 𝐻 3 𝑟𝑑 (𝑠) are defined as: where 𝐾 is a static gain, 𝑇 0 the time delay [sec], and 𝜏 1 , 𝜏 2 , 𝜏 3 the time constants [sec]. For the identification of the system, we use a trust-region-reflective algorithm on the function parameters (𝐾, 𝑇 0 , 𝜏 1 , 𝜏 2 , 𝜏 3 ) and use the root-mean-square between the data and the computed model as the error term. This identification gave us very good results (correlation > 93%) for the three different transfer functions (Table 1). As the 1st, 2nd and 3rd order transfers functions give then similar results and the addition of a time constant does not greatly increase the correlation with the data, we decided to model the caster wheel behavior with a simple 1 st order transfer function. We then applied this 1 st order transfer function to the previously defined kinematic model and we observe a good correlation of our model with the ground truth (Fig. 9).

DISCUSSION

In this study, we designed a kinematic model of the wheelchair that integrated the caster wheels effects during changes of direction. The behavior of the caster wheels has been modeled with a 1st order transfer function which appears as a good tradeoff between accuracy and complexity. This model will be used to enhance driving assistances of powered wheelchair. The aim will be to minimize the selfalignment time of the caster wheels as well as the jolt felt during maneuvers.

In this study, we used only two encoders to measure the orientation of the rear caster wheels, assuming that the behavior of the two front caster wheels was similar. In future works, we will use four encoders to investigate if there is an influence of the rear caster wheels on the front caster wheels. Another statement is that the behavior obtained with the 0° initial orientation of the caster wheels cannot be modeled because of the non-repeatability of the motion. Thus, relying on the space analysis [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF], the motion has to be forced to initially orientate the caster wheels and then use the knowledge that we have on their behavior.

CONCLUSION

The objective of this paper was to model the behavior of the caster wheels of powered wheelchairs while direction changes occurred, to enhance maneuver assistance control laws. Results have shown that the 1st order transfer function is sufficient enough to model the caster wheel behavior, leading to a simple parametric model. Thus, the real kinematics of the wheelchair with regard to the orientation of caster wheels could be computed, in order to perfectly control its motion.

In future works, we will focus on the design of a control law that assists the user during maneuvers and abrupt direction changes. In addition, we envisage to implement this model in a virtual reality powered wheelchair driving simulator [START_REF] Devigne | Design of an immersive simulator for assisted power wheelchair driving[END_REF], which takes into account the effects of the caster wheels.

Figure 1 .

 1 Figure 1. Fixed wheel

Figure 2 .

 2 Figure 2. Passively steerable caster wheel

Figure 3 .

 3 Figure 3. 𝐹 𝑊 , 𝐹 𝐴 𝑖 and 𝐹 𝐶 𝑖 frames

Figure 4 .Figure 5 .

 45 Figure 4. Kinematic relations in the front right caster wheel case

Figure 6 . 3 . 2 7 •

 6327 Figure 6. Instrumented wheelchair (Permobil M400) used for the experiment 3.2 Experimental protocol In the design of this experiment, we choose to give the following values to the study parameters: • Initial orientation of the caster wheels (α 0 ): (0°, 30°, 60°, 90°, 180°) as shown on Fig. 7 • Wheelchair velocity (V): (0.2, 0.4, 0.6) [m/sec] ; • User mass (M R relative mass): (0, 50) [Kg].During these experiments, we recorded the orientation of the caster wheels for 72 configurations of (𝛼 0 , 𝑉, 𝑀 𝑅 ) while operating a forward-backward movement. For each configuration, 5 data sets were recorded.

Figure 7 .

 7 Figure 7. The different initial caster orientations used during the experiments

Figure 8 .

 8 Figure 8. Rear caster wheel orientation evolutions and associated trajectories while going backward for 0.4 m/sec velocity and 0 Kg user mass for different initial caster orientation (𝛼 0 )

Figure 9 .

 9 Figure 9. Modeled and identified transfer functions together with ground truth data

Table 1 .

 1 Correlation between the transfer functions and the data

	Transfer	Correlation [%]
	function
	1 st order	93.680 ± 1.1698
	2 nd order	95.362 ± 1.1288
	3 rd order	95.957 ± 1.1523

Vol. 79, No. 4, December, 2018, pp. 168-172 Journal homepage: http://iieta.org/Journals/MMC/MMC_C

ACKNOWLEDGEMENTS

This work is carried out as part of the INTERREG VA FMA ADAPT project "Assistive Devices for empowering disAbled People through robotic Technologies http://adaptproject. com/index.php.". The Program is funded by the European Regional Development Fund (ERDF). Authors would also like to thank Eric Bazin, Franç ois Pasteau (INSA/IRISA), and Sylvain Rigaud.


L Devigne

F Pasteau

N Le Borgne

M Babel

T Carlson

P Gallien

Centre MPR Saint Hé lier Pole

Assisting power wheelchair driving on a sidewalk: A proof of concept

Keywords: power wheelchair, driving assistance, sensor-based control, powered mobility, negative obstacles

The use of a power wheelchair allows to maintain mobility by providing better access to daily activities and thus positive impact on the quality of life. However, driving a power wheelchair is a complex task, particularly within an environment consisting of negative obstacles (e.g. steps, sidewalk edges). In this context, falling accidents can occur while driving a power wheelchair on a sidewalk. Therefore, driving assistance is required to prevent from falling off a curb edge. In order to meet these expectations, we here propose a semi-autonomous shared control framework assisting the user while driving on a sidewalk. We present simulations as well as an experiment carried out with our system embedded on a standard wheelchair. In both cases, our method allows progressive velocity adaptation when approaching a curb edge resulting in the wheelchair avoiding the risk of falling. The obtained results thus provide a proof of concept of our method.

INTRODUCTION

Some people with motor disabilities may benefit from the use of a power wheelchair to increase their mobility and independence. However, the use of motorized mobility devices can lead to risks for others, oneself or the property [START_REF] Chen | Wheelchair-related accidents: relationship with wheelchair-using behavior in active community wheelchair users[END_REF].

In this context, many research teams are involved in the development of solutions to assist in the driving of electric wheelchairs. Indeed, driving assistance can be provided in the form of autonomous [START_REF] Burhanpurkar | Cheap or robust? The practical realization of selfdriving wheelchair technology[END_REF][START_REF] Baklouti | Autonomous wheelchair navigation with real time obstacle detection using 3D sensor[END_REF] or semi-autonomous [START_REF] Horn | Smart wheelchairs: Past and current trends[END_REF][START_REF] Ezeh | Probabilistic vs linear blending approaches to shared control for wheelchair driving[END_REF] assistance but are so far restricted to indoor use because they do not provide assistance with respect to negative obstacles.

However, negative obstacles constitute major physical barriers for safe power wheelchair driving. Indeed, a study in the form of regular users interviews, carers emphasized the difficulties encountered with curbsides [START_REF] Wang | Power mobility with collision avoidance for older adults: User, caregiver, and prescriber perspectives[END_REF]. In this context, driving assistance in outdoor environments would be expected to provide more safety while driving in an environment consisting of negative obstacles. Indeed, recent studies highlight the benefits of using motorized mobility devices in outdoor activities on user's self-esteem and sense of independence by providing easier access to outdoor activities thus maintaining daily-life activities and social interactions [START_REF] Löfqvist | Mobility and mobility-related participation outcomes of powered wheelchair and scooter interventions after 4months and 1-year use[END_REF][START_REF] Samuelsson | Powered wheelchairs and scooters for outdoor mobility: A pilot study on costs and benefits[END_REF].

In the robotics field, few methods have been proposed for assisting robots navigation in en environment constisting of negative obstacles [START_REF] Larson | Lidar based off-road negative obstacle detection and analysis[END_REF][START_REF] Murarka | Detecting obstacles and drop-offs using stereo and motion cues for safe local motion[END_REF] typically based on detection relying on stereo vision or expensive sensors such as laser rangefinders. Moreover, these solutions generally represent the presence of obstacles in an occupancy grid and are based on the definition of local environmental maps requiring odometry as well as high computational load.

In the field of power wheelchair navigation, only a few studies have focused on the detection of negative obstacles. Altough they use low-cost sensors and low-computational cost methods, the solutions presented in [START_REF] Simpson | The smart wheelchair component system[END_REF] and [START_REF] Ka | Intelligent single switch wheelchair navigation[END_REF] stop the power wheelchair when a drop-off (e.g. curbside) is detected. A detection of the edges of obstacles via stereo vision measurement is proposed in [START_REF] Coughlan | Terrain analysis for blind wheelchair users: Computer vision algorithms for finding curbs and other negative obstacles[END_REF] but this method has important limitations such as the absence of detection of certain sidewalk edges as well as false detections. In addition, the few proposed method only provide detection with no semi-autonomous or autonomous control law to assit the user while navigating with the wheelchair by another means than stopping when a dropoff is detected.

In this article, we present a low-cost semi-autonomous wheelchair control solution for assisting wheelchair navigation on a sidewalk. The proposed method is based on the same mathematical principle as the approach presented in [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF] which application has been until now restricted to indoor navigation assistance as the detection of negative obstacles was not ensured yet. This method has been tested and validated by wheelchair regular users within clinical trials whhich have been conducted within the Pôle Saint Hé lier rehabilitation center in Rennes [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF][START_REF] Babel | HandiViz project: clinical validation of a driving assistance for electrical wheelchair[END_REF]. In this paper, we propose to apply this principle to navigation on a sidewalk by considering detected drop-offs as virtual positive obstacles in the horizontal plane thus leading the wheelchair to progressively correct its trajectory to avoid the curbside the same way it avoids colliding with a wall. The control law is based on a reactive sensor-based servoing method which does not rely on the definition of a local map of the environment.

The paper is organized as follows. Section II presents the method for detecting negative obstacles. Section III presents the wheelchair velocity control method. The simulations of the proposed method are presented in Section IV. Finally, an experiment with a standard power wheelchair conducted in the UCL PAMELA facility is presented and commented on in Section V.

CURB DETECTION METHOD

In order to detect negative obstacles around the power wheelchair, we use distance sensors mounted on the wheelchair frame oriented towards the floor. Here, as the idea is to provide a proof of concept of the proposed method, the distance to the floor measured by each sensor is used to calculate the distance to the curbside in the horizontal plane. With such a detection, the wheelchair will avoid approaching a curbside the same way it already avoids to collide with a wall, thus bringing us back to the same configuration as in [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF].

Figure 1. Negative obstacle detection

As shown in Fig. 1, • let 𝑑 𝑖 the distance measured by the sensor 𝑠 𝑖 , • let ℎ 𝑠 𝑖 be the sensor 𝑠 𝑖 height from ground,

• let 𝜑 𝑖 be the orientation of the sensor with respect to the vertical, • let ℎ 𝑖 𝑛 be the height of the negative obstacle,

• let ℎ 𝑖 𝑝 be the height of the positive obstacle,

• let 𝐻 𝑖 𝑛 = ℎ 𝑠 𝑖 + ℎ 𝑖 𝑛 the sensor 𝑠 𝑖 height from negative obstacle level, • 𝐻 𝑖 𝑝 = ℎ 𝑠 𝑖ℎ 𝑖 𝑝 the sensor 𝑠 𝑖 height from positive obstacle level. The detection is such that we can distinguish 3 cases: 1) if we define hip as the minimum height of a positive obstacle we observe, then we detect a positive obstacle when , then, the distance to the obstacle is defined as 𝑥 𝑖 = 𝑑 𝑖 sin(𝜑 𝑖 ).

2) if we define 𝑖 𝑛 as the minimum height of a negative obstacle we observe, then we detect a negative obstacle when

𝑑 𝑖 > 𝑑 𝑖 𝑛 = ℎ 𝑠 +ℎ 𝑛 cos(𝜑 𝑖 ) = 𝐻 𝑖 𝑛 cos(𝜑 𝑖 )
, then, the distance to the obstacle is defined as 𝑥 𝑖 = ℎ 𝑠 𝑖 cos(𝜑 𝑖 ).

3) as long as the distance 𝑑 𝑖 measured by the sensor 𝑠 𝑖 is such that 𝑑 𝑖 𝑝 𝑑 𝑖 𝑑 𝑖 𝑛 , we can consider that there is no obstacle to detect. Then, the distance to the obstacle is defined as 𝑥 𝑖 = 𝑥 𝑚 where 𝑥 𝑚 is the maximum distance which a sensor can measure.

This detection method necessitates an extrinsic calibration of the sensors. Indeed, we use the sensor position information in our calculations.

SENSOR-BASED SERVING FRAMEWORK

In a previous paper, we presented a semi-autonomous approach providing adaptive assistance in the form of an intuitive obstacle avoidance [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF]. So far, we applied this method to indoor navigation only as we did not detect negative obstacles. Here, although the method relies on the same principles, we detail the calculation of the allowed and forbidden areas defined from sensor constraints in the wheelchair velocity domain. In particular, we emphasize the calculations for the case of corner sensors for which the geometric configuration allows us to estimate the orientation of the detected obstacle. This additional information allows the shared control scheme to achieve more comfortable trajectory corrections while approaching an obstacle.

The Modeling is given in previous work [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF]. We can define an interaction matrix 𝑳 𝒙 𝒊 for each sensor 𝒔 𝒊 such that 𝑥̇𝑖 = 𝐿 𝑥 𝑖 𝑣 𝑠 𝑖 .

(

) 1
The interaction matrix 𝐿 𝑥 𝑖 which is defined in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], is such that

𝐿 𝑥 𝑖 = [-1 tan(𝛼 𝑖 ) 0 0 0 𝑥 𝑖 tan(𝛼 𝑖 )] (2)
where 𝛼 𝑖 is the angle formed by the sensor 𝑠 𝑖 x-axis and the perpendicular to the detected obstacle (Fig. 2).

From the distance measured by the sensor 𝑠 𝑖 , we are not able to determine the orientation of the obstacle. Therefore, we assume that the obstacle detected by the sensor 𝑠 𝑖 is perpendicular to the sensor x-axis. Then, we approximate 𝛼 𝑖 0 for the computation of the interaction matrix 𝐿 𝑥 𝑖 . This leads to

𝐿 𝑥 𝑖 = [-1 0 0 0 0 0] (3)
To avoid collisions with obstacles, we constrain the distance 𝑥 𝑖 measured by the sensor 𝑠 𝑖 by a minimum value 𝜆𝑒 𝑠 𝑖 with 𝜆 > 0. We then get 𝐿 𝑥 𝑖 𝑠 𝑖 𝑇 𝑟 𝑢 -𝜆𝑒 𝑠 𝑖 .

(4)

Based on the modelling given in [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF], we find

𝐴 𝑖 𝑢 𝐵 𝑖 (5)
with 𝐴 𝑖 𝐿𝑥 𝑖 𝑠 𝑖 𝑇 𝑟 = [-𝑎 𝑖 -𝑏 𝑖 ] and 𝐵 𝑖 = -𝜆𝑒 𝑠 𝑖 . This defines a half-plane in the wheelchair velocity domain as shown on Fig. 2.

SIMULATIONS

In this section, we present the results of a simulation we performed as a first step to develop our curb-following solution.

Simulation conditions

As a part of our research activities, we developed a simulator for assisted power wheelchair driving This simulator has been designed with Unity3D and consists of a Virtual Environment in which we operate a virtual robotized wheelchair equipped with virtual distance sensors. The simulator is presented in more details in [START_REF] Devigne | Design of an immersive simulator for assisted power wheelchair driving[END_REF].

The virtual wheelchair (Fig. 3b) is here equipped with a total of 12 distance sensors placed at each corner of the virtual wheelchair (Fig. 3a). The virtual environment is displayed on a computer screen as it is only used to demonstrate the properties of the curb-following solution.

Simulation results

The virtual wheelchair starts from the position shown on Fig. 4 i.e. at the beginning of the sidewalk and parallel to the curb. The user's linear and angular inputs are constant and equal to the maximum speed of the wheelchair. This leads the wheelchair to go forward and right at maximum speed.

Figure 5. Simulation results

As shown on Fig. 5, the wheelchair starts from the beginning, no negative obstacle is detected thus the wheelchair is going towards the curb. Then, when the curb is detected, the linear velocity progressively decreases and the rotational velocity is modified in order to follow the curb. When the wheelchair drives parallel to the curb, the linear velocity increases as there is no danger to fall from the curb. The rotational velocity is maintained to zero as long as the wheelchair is detecting the curb on its right. Then, when the wheelchair arrives at the point where there is 90 ° turn to the right, it turns right and stops when it detects the drop-off at the end of the sidewalk in front of it.

EXPERIMENT

Experimental setup

We developed and tested our solution on a commercially available Typhoon II power wheelchair from Invacare. We had previously modified Invacare's DX controller to intercept the user input and adjust the control signals before sending them to the motors, allowing us to implement and test various shared control algorithms on the platform [START_REF] Ezeh | Probabilistic vs linear blending approaches to shared control for wheelchair driving[END_REF]. The algorithm itself runs on an embedded Beagle-Bone Black.

For the detection, we used Time-of-Flight (ToF) VL53L0X laser-ranging sensors from STMicroelectronics. Unlike conventional infrared sensors, these sensors are low-cost and provide accurate distance measurement whatever the target reflectance and can measure absolute distances up to 2 meters with a field of view of 25 degrees.

Our detection method necessitates an extrinsic calibration of the sensors as we use the sensor position information in our distance calculations. Here, we rigidly attached the detection modules to the wheelchair frame thus ensuring a fixed position of the sensors.

Experimental results

An example of results obtained from a trial on the sidewalk within the PAMELA facility is shown on Fig. 6. For each trial, we used the same initial conditions as during the simulations: the wheelchair starts parallel to the curb at the beginning of the sidewalk.

DISCUSSION

As shown on Fig. 6, the user first applies a rotational velocity to the wheelchair towards the curb. As long as no drop-off is detected by our system, the user rotational velocity input is directly applied to the wheelchair. Then, when a dropoff is detected, the rotational velocity is set to zero in order to avoid from falling down the curb. After that, the user maintains his input for a moment and then starts to move forward while keeping a constant rotational velocity input. The shared control rotational velocity output is then maintained to zero while the user linear velocity input is directly applied to the wheelchair as there is no risk of falling from the curb while driving forward. We can notice that at the end of the experiment, the linear velocity is abruptly set to zero while the rotational velocity is being directly applied to the wheelchair. This happened several time as the floor was irregular with gaps between the modules thus leading to false detections making the wheelchair to stop.

The experiment we performed within the PAMELA facility allowed us to assess the performances of our curb-following solution while embedded on a real wheelchair. Although we encountered false detections due to the floor irregularities, we were able to validate our detection method and to avoid from falling from the curb while driving. These results bring a proof of concept of our curb-following solution for assisted power wheelchair driving.

CONCLUSION

In this article, we have presented a generic shared control wheelchair navigation assistance solution to assist power wheelchair driving on a sidewalk. While the control law is based on the same principle as the power wheelchair collision avoidance assistance solution presented in [START_REF] Devigne | Low complex sensor-based shared control for power wheelchair navigation[END_REF], here the proposed innovative approach is based on a different detection method. Indeed, rather than detecting distances to obstacles, the proposed system measures distance with sensors oriented towards the floor. From these measurements, the distance from the wheelchair to the curb side is determined. These constraints are merged with user's input to provide a smooth and progressive shared control by means of a velocity adaptation applied to the wheelchair when approaching a curb edge. The behaviour of the system in simulation as well as embbeded on a standard power wheelchair gives a proof of concept of our method. Future works will address navigation assistance on slopes and ramps by calculating a dedicated interaction matrix for each sensor which would incorporate the local floor orientation. By doing so, we hope to expand the assistance possibilities of our method to complete outdoor navigation assistance. The method presented in this article is independent of the type os sensor used and has a low computationnaly cost, thus meeting smart power wheelchair solution design requirements.

Figure 2 .

 2 Figure 2. Half plane in the wheelchair velocity domain

Figure 3 .

 3 Figure 3. Wheelchair setup for the simulation As shown on Fig. 4, it consists of a L-shaped 15 cm height, 1.2 m wide and 5 m long sidewalk with a right turn at its end.

Figure 4 .

 4 Figure 4. Virtual environment for the simulation.B

Figure 6 .

 6 Figure 6. Experiment results: linear and angular user input(blue) and shared control output (red) velocities

Journal homepage: http://iieta.org/Journals/MMC/MMC_C


N Börlin

A Murtiyoso
email: arnadi.murtiyoso@insa-strasbourg.fr

P Grussenmeyer
email: pierre.grussenmeyer@insa-strasbourg.fr

F Menna
email: fmenna@fbk.eu

E Nocerino
email: nocerino@fbk.eu

MODULAR BUNDLE ADJUSTMENT FOR PHOTOGRAMMETRIC COMPUTATIONS

Keywords: Bundle adjustment, Camera model, Analytical Jacobians, Software, Photogrammetry, Tilt-shift lens

In this paper we investigate how the residuals in bundle adjustment can be split into a composition of simple functions. According to the chain rule, the Jacobian (linearisation) of the residual can be formed as a product of the Jacobians of the individual steps. When implemented, this enables a modularisation of the computation of the bundle adjustment residuals and Jacobians where each component has limited responsibility. This enables simple replacement of components to e.g. implement different projection or rotation models by exchanging a module. The technique has previously been used to implement bundle adjustment in the open-source package DBAT (Börlin and Grussenmeyer, 2013) based on the Photogrammetric and Computer Vision interpretations of Brown (1971) lens distortion model. In this paper, we applied the technique to investigate how affine distortions can be used to model the projection of a tilt-shift lens. Two extended distortion models were implemented to test the hypothesis that the ordering of the affine and lens distortion steps can be changed to reduce the size of the residuals of a tilt-shift lens calibration. Results on synthetic data confirm that the ordering of the affine and lens distortion steps matter and is detectable by DBAT. However, when applied to a real camera calibration data set of a tilt-shift lens, no difference between the extended models was seen. This suggests that the tested hypothesis is false and that other effects need to be modelled to better explain the projection. The relatively low implementation effort that was needed to generate the models suggest that the technique can be used to investigate other novel projection models in photogrammetry, including modelling changes in the 3D geometry to better understand the tilt-shift lens.

INTRODUCTION

The mathematical problem that is solved by the Bundle adjustment (BA) process includes a residual between two components. One residual component is the simulated projection of an object point according to the camera model at hand. The other component is computed from the corresponding image measurement 1 . As a typical presentation, consider equation (1) below (Luhmann et al., 2014, equation (4.94)) Ignoring ∆x and ∆y , equation (1) completely describes the projection of the object point (OP) (X, Y, Z) through a pinhole camera with principal point (x 0 , y 0 ) and camera constant z placed at (X0, Y0, Z0) and with an orientation given by the matrix R.

Another requirement of the bundle adjustment is the linearisation of equation ( 1) with respect to any parameter that is to be estimated (see e.g. Kraus (1993, Sec. 5.3.2), Wolf and Dewitt (2000, App. D-4), [START_REF] Mikhail | Introduction to Modern Photogrammetry[END_REF], App. C.3), or Luhmann et al. (2014, Sec. 4.4.2)). Some partial derivatives of equation ( 1) are given in equation ( 2) below (Luhmann et al., 2014, equation (4.96)). In equation (2), kX and kY are the respective numerators of equation ( 1), and N is the denominator.

What may not be immediately obvious from equation ( 1) is that the computation can be split into a sequence of basic operations. The goal of this paper is to show how that can be done and that, if each operation is treated as a module with a responsibility for computing both the result of the operation and its linearisation (Jacobians) with respect to any parameter, the computation of the analytical Jacobians can be greatly simplified.

x = x 0 + z r11(X -X0) + r21(Y -Y0) + r31(Z -Z0) r13(X -X0) + r23(Y -Y0) + r33(Z -Z0) + ∆x , y = y 0 + z r12(X -X0) + r22(Y -Y0) + r32(Z -Z0) r13(X -X0) + r23(Y -Y0) + r33(Z -Z0) + ∆y .

(1)

∂x ∂X = - z N 2 (r13kX -r11N ), ∂x ∂Y = - z N 2 (r23kX -r21N ), ∂x ∂Z = - z N 2 (r33kX -r31N ), ∂y ∂X = - z N 2 (r13kY -r12N ), ∂y ∂Y = - z N 2 (r23kY -r22N ), ∂y ∂Z = - z N 2 (r33kY -r32N ), (2)

SPLITTING THE RAYS

Basic functions

If we group the parameters as

p =   X Y Z   , p0 =   X0 Y0 Z0   , q =   kX kY N   , u = x y , u0 = x 0 y 0 , ∆u = ∆x ∆y , (3)
and introduce the basic functions

T3(p, p0) = p + p0, 3D translation (4a) L(R, v) = Rv, 3D rotation (4b)
H(q) = 1 q3 q1 q2 , Pinhole projection (4c) S(c, m) = cm, Scaling (4d) T2(n, u0) = n + u0, 2D translation, ( 4e
)
we may write equation (1) as the following composition:

u = g(p, p0, R, z , u0) ≡ T2(S(z , H(L(R T , T3(p, -p0)))), u0). (5)
In the definition of the functions in equation ( 4), generality has been favoured over typical usage to simplify future extensions. For instance, the function L(R, v) describes a 3D linear transformation whereas the typical usage is with a rotation matrix R.

Similarly, the translation T3(p, p0) is defined with a positive sign on p0 whereas the typical usage would be T3(p, -p0).

The algorithm corresponding to equation ( 5) is shown in Algorithm 1 with a graphical representation in Figure 1.

Linearisation

The linearisation of composed functions are governed by the chain rule (see Appendix A). For example, the Jacobian of the projection function g in equation ( 5) with respect to the object point Algorithm 1 Pinhole projection corresponding to equation (5).

1: procedure PINHOLENOJAC(p, p0, R, c, u0)

2:
a1 ← T3(p, -p0)

3: a2 ← L(R T , a1) 4:
a3 ← H(a2)

5:

a4 ← S(c, a3)

6:
a5 ← T2(a4, u0)

7:

return a5 8: end procedure coordinates p is the matrix product of five Jacobians

dg dp = dT2 dn n=S(z ,H(L(R T ,T 3 (p,-p 0 )))) • dS dm m=H(L(R T ,T 3 (p,-p 0 ))) • dH dq q=L(R T ,T 3 (p,-p 0 )) • dL dv v=T 3 (p,-p 0 ) • dT3 dp , (6)
where • indicates matrix multiplication and the subscripts indicate at which point the respective Jacobians are evaluated. A detailed inspection of equation ( 6) reveals that it is identical to equation (2).

To better highlight the chaining, a simplified notation that leaves the evaluation point implicit may be used:

dg dp = dT2 dS dS dH dH dL dL dT3 dT3 dp . ( 7
)

Modularisation

Algorithm 2 is an extension to Algorithm 1 to also compute the Jacobian according to Equation (6). In Algorithm 2, each function is assumed to compute both the function value proper and the Jacobians with respect to each parameter. From Algorithm 2, we see that each function will indeed have all the necessary information to compute the function value and its Jacobians at the proper Algorithm 2 Algorithm 1 extended to compute Jacobians. Each function is capable of computing the Jacobians with respect to each of its parameters. This example only shows the Jacobians necessary to compute dg dp according to equation (6). 1: procedure PINHOLEWITHJAC(p, p0, R, c, u0) Figure 2: The computational chain implemented in DBAT [START_REF] Börlin | Bundle adjustment with and without damping[END_REF]. In the photogrammetric formulation, the optical scaling in the camera results in image coordinates expressed in physical units, e.g. mm. The image coordinates are scaled from pixels to mm before the [START_REF] Brown | Close-range camera calibration[END_REF] polynomials are used to "correct" the measured image coordinates for lens distortion. The residual (thick circle) is computed as the difference between the projected ideal point and the corrected point. Figure 3: The DBAT implementation of the Computer Vision formulation of the [START_REF] Brown | Close-range camera calibration[END_REF] lens distortion model [START_REF] Tsai | A versatile camera calibration technique for high-accuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses[END_REF][START_REF] Heikkilä | A four-step camera calibration procedure with implicit image correction[END_REF][START_REF] Zhang | A flexible new technique for camera calibration[END_REF]. The Brown polynomials are used to "add" lens distortion to the projected ideal point in normalised units before the coordinates are scaled directly to pixels. The residual (thick circle) is computed between the measured point and the ideal projected point with added lens distortion. The same modules have been used as in Figure 2.

evaluation point. Thus, by adding the requirement that each function should be able to compute its own Jacobians, we may modularise the residual computations and use the functions as building blocks. Furthermore, as each module is self-contained, it is possible to validate the analytical Jacobian of each module independently.

Photogrammetric modules

In addition to the functions listed in equation ( 4), we add the computation of the rotation matrix R from the ω -φ -κ Euler angles (Förstner and Wrobel, 2004, eqs. (2.128)

-(2.130)) R(ω, φ, κ) = R3(κ)R2(φ)R1(ω), (8a)
the [START_REF] Brown | Close-range camera calibration[END_REF] lens distortion model

D(u, K, P ) = u + dr(u, K) + dt(u, P ), (8b)
and the 2D affine transformation

A2(u, b) = 1 + b1 b2 0 1 u. (8c)
In equation (8b), the vectors K and P contain the radial and tangential distortion coefficients, respectively. In equation (8c), the scalar b1 controls the aspect ratio and b2 controls the skew. For more details and derivation of the Jacobians, see Appendix B.

USAGE

The modular technique has previously been used in the opensource Damped Bundle Adjustment Toolbox (DBAT)2 package (Börlin andGrussenmeyer, 2013, 2016). In the 2016 paper, the technique was used to implement two bundle pipelines that used different adaptations of the [START_REF] Brown | Close-range camera calibration[END_REF] lens distortion model. In the Photogrammetric formulation, the Brown polynomials are used to "correct" for the effect of lens distortion on the measured image coordinates. In contrast, the formulation largely adopted by the Computer Vision community uses the same polynomials to "add" lens distortion to the ideal projection of object points [START_REF] Tsai | A versatile camera calibration technique for high-accuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses[END_REF][START_REF] Heikkilä | A four-step camera calibration procedure with implicit image correction[END_REF][START_REF] Zhang | A flexible new technique for camera calibration[END_REF]. The two pipelines are contrasted in figures 2 and 3.

In this paper, we have used the modular technique to investigate the effect of the relative ordering of an affine transformation and lens distortion on the calibration of a tilt-shift lens. Two extensions of the photogrammetric pipeline of Figure 2 was implemented in DBAT. The reference implementation (Model 2) has no affinity. In Model 3, the affinity was applied before lens distortion. In Model 4, the affinity was applied after. The corresponding functions are: where s is the pixel size and u the measured image coordinates. The models are illustrated in Figure 4. As an illustration of the (low) level of complexity needed, the difference in the computation of the Jacobian with respect to the principal point u0 was limited to a inserting a Jacobian at the proper place in the matrix product (in addition to a change in the evaluation points):

r2 = D(T2(S(s, u), -u0), K, P ), (9a
dr2 du0 = - dD du dT2 du0 , (10a)
dr3 du0 = - dD du dA2 du dT2 du0 , (10b)
dr4 du0 = - dA2 du dD du dT2 du0 . ( 10c
)

EXPERIMENTS AND RESULTS

Simulation experiment

The first experiment was constructed to investigate the effect, if any, the difference in affine-lens distortion ordering had. Two sets of synthetic data were generated, where simulated error-free image observations were generated by back-projection of known 3D object coordinates and with known exterior orientation parameters (EO), and camera calibration parameters (IO) of a strong self-calibration network. The network consisted of 24 cameras at varying roll angles. About 100 targets were simulated in a 3D configuration with a largest dimension of 1000 mm (Figure 5).

The following algorithms were used to simulate Model 3 and Model 4 (note that the order is reversed compared to Figure 4 as we are building image observations): Model 3 1. Collinearity equations (3D translation, 3D rotation, pinhole projection and optical scaling).

Introduce lens distortion (iterative, [mm]).

3. Convert from mm to pixels using a non-square pixel size corresponding to b1 = 0.01218.

Introduce the principal point [pixel].

Model 4 1. Collinearity equations (3D translation, 3D rotation, pinhole projection and optical scaling).

Apply an affine transformation of the image coordi-

nates corresponding to b1 = 0.01218.

Introduce lens distortion (iterative, [mm]).

4. Convert from mm to pixels using a square pixel size.

Introduce the principal point [pixel].

Both simulations used a skew (shear) parameter of b2 = 0. The algorithms were implemented in software developed in-house at FBK and not by DBAT.

Each synthetic data set was analysed by a self-calibration bundle adjustment using DBAT models 2, 3, and 4. The datum problem was solved by fixing the EO parameters of one camera and one coordinate of another. No control points were used and the prior weight for the image observations corresponded to a sigma of 0.1 pixels. The following parameters were estimated; the focal length, the principal point, the radial distortion parameters K1-K3, and the tangential distortion parameters P1-P2. For models 3 and 4, the affine parameters b1-b2 were also estimated. The quality of each analysis was evaluated in image space (σ0 and 2D image point RMS) and object space (3D RMSE between the true and estimated OP coordinates). The results are given in Table 1.

When the correct estimation model was used, the simulated b1 value was recovered to the number of available decimals and the internal and external residuals were effectively zero. When the wrong model was used, the residuals were significantly higher.

Camera calibration of a tilt-shift lens

In the second experiment, models 2, 3 and 4 were used to calibrate a camera with a tilt-shift lens. A tilt-shift lens allows the following movements of the lens [START_REF] Ray | Applied Photographic Optics[END_REF]):

• a tilt, i.e. a rotation of the optical axis around either the exit pupil or the center of the sensor plane,

• a shift, i.e. a translation of the optical axis, and

• a rotation, i.e. a rotation about the optical axis.

The data set from Nocerino et al. ( 2016) was used for the calibration. The data set was acquired by a Nikon D750 full-frame DSLR camera with a PC-E Micro NIKKOR 45mm f/2.8D ED tilt-shift lens (Figure 6) in two different configurations:

NORMAL The normal configuration where neither tilt nor shift was applied.

TILTED The lens was tilted in the vertical plane by 4 degrees.

The calibration target was a 3D photogrammetric calibration test object (Figure 7) with about 170 coded targets with a largest dimension of 900 mm. A high redundancy photogrammetric camera network was realised, consisting of up to 48 convergent images, with a diversity of camera roll angles to enhance the determinability of the IO parameters [START_REF] Fraser | Photogrammetric camera component calibration: A review of analytical techniques[END_REF].

The NORMAL data set was analysed and used as a reference for the processing of the TILTED data set. The TILTED data set was analysed by a self-calibration bundle adjustment in DBAT using models 2, 3 or 4. The datum and prior weights were the same as in the synthetic experiment. Furthermore, the parameters b1, b2, and P were individually included or excluded from the estimation (P1 and P2 were always estimated together). The quality of the estimation was evaluated as in the synthetic experiment with the results of the NORMAL data set used as the reference. The results are given in Table 2.

From Table 2, we may conclude that the difference between models 3 and 4 on real-world data is small. The internal and external residuals were small when the aspect parameter b1 and the decentering distortion parameters P1-P2 were included in the esti-mation. In those cases, the estimated b1 value was about 0.0019. This value is equal to the value in Nocerino et al. (2016, Table 3, column T-4S0R90) that was estimated by other software. In contrast, the inclusion of b2 had a negligible effect on the residuals.

DISCUSSION

In this paper we discuss how the residual computations used by bundle adjustment can be split into modules. If each module is responsible for computing Jacobians with respect to each parameter, in addition to the function value proper, Jacobians of complex expressions can be computed using the chain rule. Furthermore, the analytical Jacobians of each module can be validated independently of the other modules.

The modular technique has previously been used in the Damped Bundle Adjustment Toolbox (DBAT) to model the Photogrammetric and Computer Vision adaptations of the [START_REF] Brown | Close-range camera calibration[END_REF] lens distortion models [START_REF] Börlin | External verification of the bundle adjustment in photogrammetric software using the damped bundle adjustment toolbox[END_REF]. In this paper, the Photogrammetric pipeline was extended by an affine transformation. Two models with different placement of the affine transformation compared to lens distortion were implemented with minimal effort.

An experiment on synthetic image observations was performed.

The conclusion is that the relative placement of the affine transformation and lens distortion does matter and that DBAT was able to distinguish which model was used to generate the synthetic data.

The tilt-shift lens is a complex design whose projection model

is not yet rigorously supported by DBAT. In a previous paper, it was found that some of the deviation of the tilt-shift lens from the classic photogrammetric projection model (pin-hole with Brown lens distortion) could be absorbed by including the affine distortion parameter b1 in the estimation [START_REF] Nocerino | Experiments on calibrating tilt-shift lenses for close-range photogrammetry[END_REF]. In this paper, the calibration experiment was constructed to test the hypothesis that the order of the affine and lens distortion transformations could be modified to explain more of the tilt-shift distortion. The results suggest that the hypothesis was false -compared to other, possibly unmodelled effects, the ordering does not significantly contribute to the observed residuals. The results however support the conclusion in [START_REF] Nocerino | Experiments on calibrating tilt-shift lenses for close-range photogrammetry[END_REF] that the aspect parameter b1 and the decenter parameters P1-P2 are significant but the skew parameter b2 is not.

The primary goal of this paper was to show how the presented modularity of DBAT can be used to test novel projection models. The first experiment shows that the ordering does have an effect and is detectable if it is the dominant factor. The second experiment shows that the tested re-ordering has little effect, thus suggesting that the ordering is not the dominant factor. In both cases, the implementation effort to test the hypothesis was minimal. This suggests that DBAT can be used to test novel projection models in the future. Potential uses of the modularity include replacing a module to form another projection model. For example, the pin-hole module could be swapped by a fish-eye module, or rotation modules could be swapped between omega-phi-kappa and azimuth-tilt-swing.

Future work to better understand the tilt-shift lens include modelling changes in the 3D geometry of the lens with DBAT.

A. THE CHAIN RULE

According to the chain rule, if a univariate function h(x) is formed as the composition of two univariate, differentiable functions f (x) and g(x) as

h(x) = f (g(x)), (11)
the derivative h (x) of the composed function may be calculated as the product of the derivatives f (x) and g (x). With substitutions y = f (u), u = g(x), the derivative may be written as

dy dx = dy du du dx , (12)
where the re-appearance of the denominator du of one factor as the numerator of the next gives the chain rule its name. To make it clear where each derivative is computed, the computation of h (x) at x = c is usually written as

dy dx x=c = dy du u=g(c) du dx x=c , (13)
where the subscript is to be read "evaluated at".

The chain rule may be extended to multivariate functions. For instance, if

g(x(u, v), y(u, v), z(u, v))
is a function of (x, y, z) that themselves are functions of (u, v), the chain rule dictates that

is the product of the two Jacobians dg dt and dt dw .

B. JACOBIANS B.1 Preliminaries

This section mostly follows [START_REF] Magnus | Matrix Differential Calculus with Applications in Statistics and Econometrics[END_REF].

B.1.1 The vec(•) operator Given an m × n matrix

A = a1 • • • an =    a11 • • • a1n . . . . . . . . . am1 • • • amn    ∈ m×n , ( 17
)
the vectorisation operator vec(•) rearranges A into a column vector where the elements of A are in column-major order

vec(A) =    a1 . . . an    =      a11 a21 . . . amn      ∈ mn×1 . (18)
B.1.2 Jacobians of matrix functions The Jacobian of any scalar-or vector-valued function f with respect to a matrix argument A is implicitly assumed to be with respect to vec(A), i.e.

df (A) dA ≡ df (A) d vec(A) . (19)
Furthermore, the Jacobian of any matrix-valued function F (A) is implicitly assumed to be applied to vec(F (A)), i.e.

dF (A) dA ≡ d vec(f (A)) d vec(A) . (20)
B.1.3 The Kronecker product The Kronecker product C = A ⊗ B, where the matrices A ∈ m×n , B ∈ p×q , and C ∈ mp×nq , is defined as

C = A ⊗ B =    a11B • • • a1nB . . . . . . . . . an1B • • • annB    . (21)
B.1.4 The Jacobian of matrix expressions For matrices A ∈ k×l , B ∈ l×m , C ∈ m×n , the following identities hold:

vec(AB) = (Im ⊗ A) vec(B) (22a) = (B T ⊗ I k ) vec(A) (22b)
and

vec(ABC) = (C T ⊗ A) vec(B) (23a) = (In ⊗ AB) vec(C) (23b) = (C T B T ⊗ I k ) vec(A). (23c)
Equations ( 22) and ( 23) can be used to derive Jacobians of matrix products.

The Jacobians of T3 with respect to p and p0 are

dT3(p, p0) dp = dT3(p, p0) dp0 = I3. (25b)
If the translation is applied to m points stored as columns in P , we instead get

T3(P, p0) = P + p01 T m , (25c)
and

dT3(P, p0) dP = Im ⊗ I3 = I3m (25d) dT3(P, p0) dp0 = 1m ⊗ I3. (25e)
Similarly for a 2D point u and offset u0,

T2(u, u0) = u + u0, dT2(u, u0) du = dT2(u, u0) du0 = I2. (26)
B.2.3 3D linear transformation An arbitrary linear mation of a point p ∈ 3 can be formulated as

L(M, p) = M p, (27a)
with Jacobians given by ( 22)

dL(M, p) dM = p T ⊗ I3, dL(M, p) dp = M. (27b)
B.2.4 3D rotation matrix If the 3D rotation matrix is defined using the ω -φ -κ Euler angles (Förstner and Wrobel, 2004, eqs. (2.128)-(2.130)), and with the vector k = ω φ κ T , we have the following rotations

R1(ω) =   1 0 0 0 cos ω -sin ω 0 sin ω cos ω   , (28a)
R2(φ) =   cos φ 0 sin φ 0 1 0 -sin φ 0 cos φ   , (28b)
R3(κ) =   cos κ -sin κ 0 sin κ cos κ 0 0 0 1   , (28c)
R(k) = R(ω, φ, κ) = R3(κ)R2(φ)R1(ω), (28d)
with Jacobians given by Lucas (1963, eqs. (3)-( 9))

dR(k) dω = -R(ω, φ, κ)Px, (28e) dR(k
) dφ = -R3(κ)R2(φ)PyR1(ω), (28f) dR
(k) dκ = -PzR(ω, φ, κ), (28g)
dR(k) dk = dR(k) dω dR(k) dφ dR(k) dκ . (28h)
B.2.5 Pin-hole projection The pin-hole projection of a 3D point p to a 2D point is given by

H(p) = 1 p3 p1 p2 , (29a) dH(p
The vp function expands a scalar value x to a vector of its powers:

vp(x; n) =     x x 2 • • • x n     , dvp(x; n) dx =     1 2x 1 • • • nx n-1     , (32b)
The radial scaling function rs computes the inner product between a coefficient vector c and the power vector of the radial values squared: Given the helper functions above, the radial distortion part of Brown (1971, equation (20)) is reduced to dr(u, K) = urs(u, K),

Figure 1 :

 1 Figure 1: The computational chain corresponding to equation (5). The gray circles indicate bundle adjustment parameters. The white circles indicate component functions in equation (4). The arrows indicate how the parameters and results are propagated. The arrow labels indicate the name of the formal (input) parameter of the corresponding function.

  Jm) ← S(c, a3) Jm = dS dm m=a 3 6: (a5, Jn) ← T2(a4, u0) Jn = dT 2

Figure 4 :

 4 Figure 4: Two versions of the right-hand side of Figure2with different relative placement of the affine transformation with respect to lens distortion. In Model 3 (left), the affine distortion was applied before lens distortion. In Model 4 (right), the affine distortion was applied after lens distortion. Note that in both cases, the digital scaling uses square pixel sizes. Any non-unit aspect ratio is handled by the affine step.

Figure 5 :

 5 Figure 5: The synthetic network consisted of 24 cameras at varying roll angles. The simulated targets were in a 3D configuration with a largest dimension of 1000 mm.

Figure 6 :

 6 Figure 6: The Nikon D750 DSLR camera with the PC-E Micro NIKKOR 45 mm f/2.8D ED tilt-shift lens used in this paper.

Figure 7 :

 7 Figure 7: The 3D test object consisted of about 170 coded targets with a largest dimension of 900 mm.

B. 2

 2 Component functions B.2.1 The transpose The transpose of an m-by-n matrix A B(A) = A T (24a) is a permutation of the elements of A. The Jacobian is a permutation matrix known as the Commutation matrix Kmn dB dA = Kmn. (24b) B.2.2 Translation The translation of a point p ∈ 3 in by an offset c ∈ 3 is given by T3(p, p0) = p + p0.

  scaling The scaling of the 2D point u by a scalar k is given by S(k, u) = ku. where b are the affine parameters. The corresponding 2D affine transformation function A2 is defined as A2(u, b) = Am(b)u, To simplify the expressions below, we define a number of helper functions: The r function is the norm ("radius") squared of a vector: r(u) = u T u, dr(u) du = 2u T ,

  rs(u, c) = c T vp(r(u); |c|), (32c)drs(u, c) du = c T dvp(r(u); |c|) dr dr(u) du ,(32d)drs(u, c) dc = vp(r(u); |c|) T , (32e)where |c| is the number of elements of the vector c. Finally, the tangential scaling function ts computes the non-radial part of the tangential distortionts(u, p) = (u T uI2 + 2uu T )p,(32f)dts(u, p) du = 2pu T + 2p T uI2 + 2up T ,(32g)dts(u, p) dp = u T uI2 + 2uu T . (32h)

Table 1 :

 1 Assessment of the analysis of the synthetic data sets. The IO parameters columns indicate what IO parameters and how many (n) were included in the estimation. The point RMS is the average residual over all image observations. The 3D RMSE is the average error over all object points. The b1 column shows the estimated b1 value.

	Simulated DBAT	IO parameters		Internal assessment	External assessment
	dataset	model parameters	n	b1	σ0	Point RMS (pixels)	3D RMSE (µm)
	3	2	b1 b2 P	8	0	34.2847	4.56	620.60
		3	b1 b2 P	10 0.01218 00.0003	0.00	000.01
		4	b1 b2 P	10 0.01223 00.3773	0.05	008.26
	4	2	b1 b2 P	8	0	34.1330	4.54	616.59
		3	b1 b2 P	10 0.01212 00.3697	0.05	008.33
		4	b1 b2 P	10 0.01218 00.0003	0.00	000.01

Table 2 :

 2 Assessment of DBAT estimation models 2, 3, and 4 on the real-world data set. The IO parameters and assessment are as in Table1, except that the target coordinates computed from the NORMAL data set was used as the true data. The estimated b1 value was about 0.0019 for all green rows. The corresponding value of b2, when estimated, was about 10 -5 .

	DBAT	IO parameters			Internal assessment	External assessment
	model b1 b2 P	n	σ0 Point RMS (pixels)	RMSE (µm)
	2		P	8	5.9	0.79	82.3
				6	5.9	0.80	92.9
	3	b1 b2 P 10 1.1	0.15	11.9
		b1	P	9	1.1	0.15	11.9
		b1 b2		8	3.0	0.41	73.2
		b1		7	3.0	0.41	72.9
		b2 P	9	5.9	0.79	82.3
		b2		7	5.9	0.80	92.6
	4	b1 b2 P 10 1.1	0.15	11.9
		b1	P	9	1.1	0.15	11.9
		b1 b2		8	3.0	0.41	73.3
		b1		7	3.0	0.41	73.3
		b2 P	9	5.9	0.79	82.3
		b2		7	5.9	0.80	92.8

https://github.com/niclasborlin/dbat