Franck Plestan
email: franck.plestan@irccyn.ec-nantes.fr.

Yuri Shtessel

Vincent Bregeault

Y Shtessel ++

V Brégeault

A Poznyak

New methodologies for adaptive sliding mode control

HAL is

Introduction

The sliding mode control is a very popular strategy for control of nonlinear uncertain systems, with a very large frame of applications fields [START_REF] Slotine | Applied nonlinear control[END_REF][START_REF] Utkin | Sliding mode in control in electromechanical systems[END_REF]. Due to the use of discontinuous function and high control gain, its main features are the robustness of closed-loop system and the finite-time convergence. However, its design requires the knowledge of uncertainties bound, which could be, by a practical point-of-view, a hard task; it often follows that this bound is over-estimated, which yields to excessive gain. Then, the main drawback of the sliding mode control, the well-known chattering phenomenon (for its analysis, see [START_REF] Boiko | Analysis of chattering in continuous sliding-mode controllers[END_REF][START_REF] Boiko | Analysis of chattering in systems with second order sliding modes[END_REF]), is important and could damage actuators and systems. A first way to reduce the chattering is the use of a boundary layer: in this case, many approaches have proposed adequate controller gains tuning [START_REF] Slotine | Tracking control of nonlinear system using sliding surfaces, with application to robot manipulators[END_REF]. A second way to decrease the chattering phenomenon is the use of higher order sliding mode controller [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF][START_REF] Bartolini | On multi-input chattering-free secondorder sliding mode control[END_REF][START_REF] Levant | Universal SISO sliding-mode controllers with finite-time convergence[END_REF][START_REF] Levant | Principles of 2-sliding mode design[END_REF][START_REF] Laghrouche | Higher order sliding mode control based on integral sliding surface[END_REF][START_REF] Plestan | A new algorithm for high order sliding mode control[END_REF]. However, in these both control approaches, knowledge of uncertainties bound is required. As the objective is the not-requirement of the uncertainties bound, an other way consists in using adaptive sliding mode, the goal being to ensure a dynamical adaptation of the control gain in order to be as small as possible whereas sufficient to counteract the uncertainties/perturbations. As recalled previously, this problem is an exciting challenge for applications given that, in many cases, gains are also over-estimated, which gives larger control magnitude and larger chattering. In order to adapt the gain, many controllers based on fuzzy tools [START_REF] Munoz | An adaptative sliding mode controller for discrete nonlinear systems[END_REF][START_REF] Tao | Adaptive fuzzy sliding mode controller for linear systems with mismatched time-varying uncertainties[END_REF] have been published ; however, these papers do not guarantee the tracking performances. In [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF], gain dynamics directly depends on the tracking error (sliding variable): the control gain is increasing since sliding mode is not established. Once it is the case, gain dynamics equals 0. The main drawback of this approach is the gain over-estimation with respect to uncertainties bound (see for example in [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF]-Section III). Furthermore, this approach is not directly applicable, but requires modifications for its application to real systems: thus, the sign function is replaced by a saturation function whom the boundary layer width affects accuracy and robustness. Furthermore, no boundary layer width tuning methodology is provided. A method proposed in [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] in order to limit the switching gain must be mentioned. The idea is based on use of equivalent control: once sliding mode occurs, disturbance magnitude is evaluable and allows an adequate tuning of control gain. However, this approach requires the knowledge of uncertainties/perturbations bounds and the use of low-pass filter, which introduces signal magnitude attenuation, delay, and transient behavior when disturbances are acting. In [START_REF] Hall | Sliding mode disturbance observer-based control for a reusable launch vehicle[END_REF], a gain-adaptation algorithm is proposed by using sliding mode disturbance observer. The main drawback is that the knowledge of uncertainties bounds is required to design observer-based controller. The objective of the current paper is to propose new methodologies for control of a class of uncertain nonlinear systems for which uncertainties are bounded but this bound (which is finite) is not known. The first controller is based on an gain adaptation law derived from coupling of [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF][START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] methods, whereas the second one is using an original gain adaptation law. In the both control algorithms, the gain is dynamically tuned in order to ensure the establishment of a sliding mode; once this sliding mode is achieved, the gain is adjusted in order to get a "sufficient" value in order to counteract the perturbations and uncertainties. Section 2 states the problem, displays some recalls on sliding mode and presents a review of [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF][START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] by addressing their deficiencies. Section 3 displays the both proposed adaptive sliding mode control algorithms. In order to prove the feasibility of the approaches, both simulations (control of tutorial system and position control of electropneumatic actuator model) are presented in Section 5.

Preliminaries

Problem statement

Consider the nonlinear uncertain system

ẋ = f (x) + g(x) • u (1)
with x ∈ X ⊂ IR n the state vector and u ∈ IR the control input. Functions f (x) and g(x) are smooth uncertain functions and are bounded for x ∈ X ; furthermore, f (x) contains unmeasured perturbations term and g(x) = 0 for x ∈ X (i.e. system (1) is controllable for all x ∈ X ). The control objective consists in forcing the continuous function σ(x, t), named sliding variable, to 0. Supposing that σ admits a relative degree equal to 1 with respect to u, one gets

σ = ∂σ ∂x ẋ + ∂σ ∂t = ∂σ ∂t + ∂σ ∂x • f (x) Ψ(x, t) + ∂σ ∂x • g(x)
Γ(x, t)

•u = Ψ(x, t) + Γ(x, t) • u (2)
Functions Ψ(x, t) and Γ(x, t) are supposed to be such that

|Ψ| ≤ Ψ M , 0 < Γ m ≤ Γ ≤ Γ M (3)
for x ∈ X . It is assumed that Ψ M , Γ m and Γ M exist but are not known. The objective of the paper is to propose a new sliding mode controller u(σ, t) with the same features as classical SMC (robustness, finite time convergence) but without any information on uncertainties and perturbations (appearing in f (x)); this latter is only known to be bounded. Furthermore, this objective allows to ensure a global stability of closed-loop system whereas the classical way (with knowledge of uncertainties bounds) only ensures its semi-global stability. In the sequel, definitions of ideal and real sliding mode are recalled.

Definition 1 [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] Consider the nonlinear system (1), and let the system be closed by some possibly-dynamical discontinuous feedback. Variable σ is continuous function, and the set

S = {x ∈ X | σ(x, t) = 0},
called "sliding surface", is non-empty and is locally an integral set in the Filippov sense [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF],

i.e. it consists of Filippov's trajectories of the discontinuous dynamical system. The motion on S is called "sliding mode" with respect to the sliding variable σ.

In real applications, an "ideal" sliding mode as defined in Definition 1 can not be established. Then, it is necessary to introduce the concept of "real" sliding mode.

Definition 2 [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] Given the sliding variable σ(x, t), the "real sliding surface" associated to (1) is defined as (with δ > 0)

S * = {x ∈ X | |σ| < δ}. (4)
Definition 3 [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] Consider the non-empty real sliding surface S * given by (4), and assume that it is locally an integral set in the Filippov sense. The corresponding behavior of system (1) on ( 4) is called "real sliding mode" with respect to the sliding variable σ(x, t).

Discussion of adaptive sliding mode control solutions and motivations

Two adaptive sliding mode controllers have been given in [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF] and [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] and are going now to be discussed. These two approaches follow two different ways : the first does not require the knowledge of uncertainties/perturbations bound and consists in increasing the control gain since a sliding mode is established whereas the second one requires the knowledge of the bound and is using the equivalent control concept in order to evaluate and to minimize the control gain. Then, a same objective can be claimed, which is the dynamical adaptation of the gain in order to counteract perturbations and uncertainties.

Theorem 1 ( [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF]) Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by

u = -K(t) • sign(σ(x, t))
with the adaptation law K(t) given by

K = K • |σ(x, t)| (5)
with K > 0 and K(0) > 0, then there exists a finite time t F ≥ 0 so that a sliding mode is established in system for all t ≥ t F , i.e. σ(x, t) = 0 for t ≥ t F .

Discussion. The main feature of this approach is that it does not require a priori the knowledge of control gain. However, from K-dynamics, it yields that when σ = 0, K = 0.

In this case, the gain K is clearly over-estimated with respect to uncertainties, which induces larger chattering (see simulations results in [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF]-Section III). Furthermore, this methodology is applicable only for "ideal" sliding mode, the objective σ = 0 being reachable. In case of real sliding mode, this latter objective is not reachable which induces that K gain is always increasing. In [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF], authors propose to modify K-dynamics by introducing a boundary layer neighboring the sliding surface σ = 0. It means that accuracy has to be sacrificed in order to apply the previous controller and that the control gain is still over-estimated. Also, with discrete-time measurements, adaptation law (5) inevitably implies the tending of K to infinity, if the sampling interval is separated from zero.

Theorem 2 ( [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF]) Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by

u = -K(t) • sign(σ(x, t))
with the adaptation law K(t) given by

K(t) = K • |η| + χ
with K > 0, χ > 0 and η the average of sign(σ) obtained through a low pass-filter

τ • η + η = sign(σ(x, t)) with τ > 0, if K ≥ | Ψ Γ |, (6)
then there exists a finite time t F > 0 so that the sliding mode is established for all t ≥ t F .

Discussion. The main feature of this controller is the adjustment of the control gain by using the equivalent control concept. It implies that chattering is decreasing. However, Kadaptation law requires the knowledge of uncertainties/perturbations bounds. Furthermore, the use of low-pass filter introduces in the closed-loop system dynamics and transient phenomena in case of perturbations. Note also that methodologies for the tuning of parameters τ and χ have not been detailed except their positivity and the fact that the time constant of the low pass filter τ must be sufficiently small. However, their tuning could engender very different behaviors of closed-loop system.

Motivations. The main objective of the current paper consists in providing new adaptive sliding mode controllers such that

• uncertainties/perturbation bounds exist but are not known.

• the gain adaptation law does not overestimate uncertainties/perturbations magnitude and then, the obtain control magnitude is reasonable.

A first solution to the problem under interest in this paper consists in combining the two previous gain adaptation laws: the first one would be used since a sliding mode is established, whereas the second allows to decrease the control gain once sliding mode is established. A second solution is an original adaptive real sliding mode control law that allows not to overestimate the gain value without using concept of equivalent control.

Main results

In this section, two adaptive sliding mode control algorithms that establish the sliding modes in SISO systems with bounded perturbations which bounds are unknown, are presented. The first algorithm combines the adaptive schemes formulated in Theorems 1-2 -Section 2, while mitigating the deficiencies of the combined gain-adaptation schemes. The second adaptive sliding mode control algorithm is completely original. It does not require the estimation of the perturbations via equivalent control as in [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] (Theorem 2) and does not overestimate the control gain as in [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF] (Theorem 1). Furthermore, the second adaptive-gain sliding mode control algorithm requires smaller amount of tuning parameters than the first algorithm, and is developed in the real sliding mode context.

First adaptive sliding mode control law

Consider the following controller

u = -K • sign(σ(x, t)) (7)
with the gain K(t) defined by (ǫ being a positive parameter)

• If |σ| = 0, K(t) is the solution of K = K1 • |σ(x, t)| (8)
with K1 > 0 and K(0) > 0.

• If σ = 0, K(t) reads as

K(t) = K2 • |η| + K3 τ η + η = sign(σ(x, t)) (9)
with K2 = K(t * ), K3 > 0 and τ > 0. t * is the largest time value such that, by denoting t * -the time just before t * , σ(x(t * -), t * -) = 0 and σ(x(t * ), t * ) = 0.

Discussion. By supposing that |σ(x(0), 0)| = 0, the adaptive sliding mode control law ( 7)-( 9) works as follows

• The gain K(t) is increasing due to the adaptation law (8) up to a value large enough to counteract the bounded uncertainty with unknown bounds in (1) until the sliding mode starts. Denote the time instant when the sliding mode starts for the first time as t 1 .

• As sliding mode has started, i.e. σ(x(t), t) = 0, from t = t 1 , K(t) follows the gain adaptation law [START_REF] Hall | Sliding mode disturbance observer-based control for a reusable launch vehicle[END_REF]. Then, gain K(t) is adapted through [START_REF] Hall | Sliding mode disturbance observer-based control for a reusable launch vehicle[END_REF] with K2 = K(t 1 ) . Note that this strategy will allow to decrease the gain and then to adjust it with respect to the current uncertainties/perturbations.

• However, if the varying uncertainty/perturbation exceeds the value K2 = K(t 1 ), then the sliding mode will be destroyed and σ(x(t), t) becomes not equal to zero. Next, the gain adaptation will happen in accordance with [START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF]. The gain K(t) will be increasing until the sliding mode occurs again at the reaching time instant t 2 .

• As sliding mode has occured and σ(x(t), t) = 0 from t = t 2 , K(t) now follows the gain adaptation law ( 9) with K2 = K(t 2 ). And so on ...1

Lemma 1 Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by ( 7),( 8),( 9), the gain K(t) has an upper-bound, i.e. there exists a positive constant K * so that

K(t) ≤ K * , ∀t > 0.
Proof of Lemma 1 is in Appendix section.

Theorem 3 Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by ( 7),( 8),( 9), there exists a finite time t F > 0 so that a sliding mode is established for all t ≥ t F , i.e. σ(x, t) = 0 for t ≥ t F .

Proof. The proof is composed by two steps. The first step concerns σ = 0 which yields K-dynamics described by [START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF], whereas the second one concerns σ = 0 with K-dynamics described by [START_REF] Hall | Sliding mode disturbance observer-based control for a reusable launch vehicle[END_REF].

• Suppose that σ(0) = 0: in this case, K(t)-dynamics read as [START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF]. Consider the following Lyapunov candidate function, with K * a positive constant

V = 1 2 σ 2 + 1 2γ (K -K * ) 2 (10)
One has

V = σ • Ψ -σ • Γ • K • sign(σ) + 1 γ (K -K * ) • K1 • |σ| ≤ Ψ M • |σ| -Γ m • K • |σ| + 1 γ (K -K * ) • K1 • |σ| = Ψ M • |σ| -Γ m • K • |σ| + Γ m • K * • |σ| -Γ m • K * • |σ| + 1 γ (K -K * ) • K1 • |σ| = (Ψ M -Γ m • K * ) • |σ| + (K -K * ) • -Γ m • |σ| + K γ • |σ|
From Lemma 1, there always exists K * > 0 such that K(t) -K * < 0 for all t > 0. It yields

V = -(Γ m • K * -Ψ M ) β σ •|σ| --Γ m • |σ| + K γ • |σ| β K •|K -K * | (11)
There always exists K * and γ such that K * > Ψ M Γm and γ < K Γm , which yields β σ > 0 and β K > 0. Then, one gets

V = -β σ • √ 2 |σ| √ 2 -β K • 2γ |K -K * | √ 2γ ≤ -min{β σ √ 2, β K 2γ} |σ| √ 2 + |K -K * | √ 2γ ≤ -β • V 1/2 (12)
with β = √ 2 min{β σ , β K √ γ}. Therefore, finite time convergence to a domain σ = 0 is guaranteed from any initial condition |σ(0)| > 0 , and the reaching time t r can be easily estimated as

t r ≤ 2V (0) 1/2 β .
• Suppose now that σ(0) = 0: Theorem 1 of [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] claims that, if K is large enough with respect to uncertainties/perturbations effects, then sliding mode control [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF] with the gain adaptation algorithm (9) allows keeping trajectories of system (1) on the sliding surface σ = 0. Given that trajectories of system (1) reach σ = 0, it means that K is large enough as required by Theorem 1 of [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF]. Then, the sliding mode is established in system (1) for all t ≥ t r . Theorem 3 is proven.

Implementation issues. The proposed adaptive-gain sliding mode control algorithm in ( 7)-( 8)-( 9) is not ready for the practical implementations. In fact, it is not possible to reach the objective σ = 0 due to sampled computation, noisy measurements or other non-idealities.

That is why it is important considering the implementation of the previous controller in a real sliding mode context. The following modifications of the gain-adaptation algorithm ( 7)-( 8)-( 9) is proposed2

• If |σ(x, t)| > ǫ > 0, K(t) is the solution of K = K1 • |σ(x, t)| (13)
with K1 > 0 and K(0) > 0.

• If |σ(x, t)| ≤ ǫ, K(t) reads as

K(t) = K2 • |η| + K3 τ η + η = sign(σ(x, t)) (14)
with K2 = K(t * ), K3 > 0 and τ > 0. t * is the largest time value such that, by denoting t * -the time just before t * , |σ(x(t * -), t * -)| > ǫ and |σ(x(t * ), t * )| ≤ ǫ.

Corollary 1 Given the nonlinear uncertain system (1) with the sliding variable σ(x, t)) dynamics (2) controlled by ( 7),( 13),( 14), there exists a finite time t F > 0 so that a real sliding mode is established for all t ≥ t F , i.e. |σ(x, t)| < ǫ for t ≥ t F .

Proof. The proof is composed by two steps. The first step concerns |σ| > ǫ which yields K-dynamics described by [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF], whereas the second one concerns |σ| ≤ ǫ with K-dynamics described by [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]. Then, the proof follows the same way than the proof of Theorem 3 based on Lyapunov analysis when |σ| > ǫ, and on Theorem 1 of [START_REF] Lee | Chattering suppression methods in sliding mode control systems[END_REF] when |σ| ≤ ǫ.

Second adaptive sliding mode control law

The first adaptive sliding mode control law uses concept of equivalent control which introduces low-pass filter dynamics with τ parameter that is not easy to tune [START_REF] Utkin | Sliding mode in control in electromechanical systems[END_REF]. The controller displayed in this section does not estimate the boundary of perturbation and uncertainties. But, there is an eminent price to do that: the new strategy guarantees a real sliding mode only. Consider the following controller

u = -K • sign(σ(x, t)) (15)
with the gain K(t) defined such that

K = K • |σ(x, t)| • sign(|σ(x, t)| -ǫ) if K > µ µ if K ≤ µ (16)
with K(0) > 0, K > 0, ǫ > 0 and µ > 0 very small. The parameter µ is introduced in order to get only positive values for K. In the sequel, for discussion and proof, and without loss of generality but for a sake of clarity, one supposes that K(t) > µ for all t > 0.

Discussion. Once sliding mode with respect to σ(x, t) established, the proposed gain adaptation law (16) allows the gain K declining (while |σ(x, t)| < ǫ). In the other words, the gain K will be kept at the smallest level that allows a given accuracy of σ-stabilization. Of course, as described in the sequel, this adaptation law allows to get an adequate gain with respect to uncertainties/perturbations magnitude.

Lemma 2 Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by ( 15),( 16), the gain K(t) has an upper-bound, i.e. there exists a positive constant K * so that

K(t) ≤ K * , ∀t > 0.
Proof of Lemma 2 is in Appendix section.

Theorem 4 Given the nonlinear uncertain system (1) with the sliding variable σ(x, t) dynamics (2) controlled by ( 15),( 16), there exists a finite time t F > 0 so that a real sliding mode is established for all t ≥ t F , i.e. |σ(x, t)| < δ for t ≥ t F , with

δ = ǫ 2 + Ψ 2 M KΓ m . ( 17
)
Proof. The proof is based on Lyapunov's approach and shows that, when |σ| > ǫ, then control strategy ensures that |σ| < ǫ in a finite time. Furthermore, it is proven that as soon as σ reaches the domain | σ | ≤ ǫ, it stays in the domain | σ | ≤ δ defined by ( 17) for all consecutive time. Therefore, the proof shows that the real sliding mode is established in finite time in the domain | σ | ≤ δ.

Consider the following Lyapunov candidate function

V = 1 2 σ 2 + 1 2γ (K -K * ) 2 (18)
One has

V = σ • Ψ -σ • Γ • K • sign(σ) + 1 γ (K -K * ) • K • |σ| • sign (|σ| -ǫ) ≤ Ψ M • |σ| -Γ m • K • |σ| + 1 γ (K -K * ) • K • |σ| • sign (|σ| -ǫ) = Ψ M • |σ| -Γ m • K • |σ| + Γ m • K * • |σ| -Γ m • K * • |σ| + 1 γ (K -K * ) • K • |σ| • sign (|σ| -ǫ) = (Ψ M -Γ m • K * ) • |σ| + (K -K * ) • -Γ m • |σ| + K γ • |σ| • sign (|σ| -ǫ) Introduce parameter β K > 0 as V = (Ψ M -Γ m • K * ) • |σ| + (K -K * ) • -Γ m • |σ| + K γ • |σ| • sign (|σ| -ǫ) +β K • |K -K * | -β K • |K -K * | (19)
From Lemma 2, there always exists K * > 0 such that K(t) -K * < 0 for all t > 0. It yields

V = -(-Ψ M + Γ m • K * ) β σ > 0 •|σ| -β K • |K -K * | --Γ m • |σ| + K γ • |σ| • sign (|σ| -ǫ) -β K • |K -K * | ξ (20)
Then, one gets

V = -β σ • |σ| -β K • |K -K * | -ξ = -β σ • √ 2 |σ| √ 2 -β K • 2γ |K -K * | √ 2γ -ξ ≤ -min{β σ √ 2, β K 2γ} |σ| √ 2 + |K -K * | √ 2γ -ξ ≤ -β • V 1/2 -ξ (21)
with β = √ 2 min{β σ , β K √ γ}.

• Case 1. Suppose that |σ| > ǫ. ξ is positive if

-Γ m • |σ| + K γ • |σ| -β K > 0 ⇒ γ < K • ǫ Γ m • ǫ + β K (22)
From ( 21), one gets

V ≤ -β • V 1/2 -ξ ≤ -β • V 1/2 (23)
It is always possible to choose γ such that the previous inequality fulfills. Therefore, finite time convergence to a domain |σ| ≤ ǫ is guaranteed from any initial condition |σ(0)| > ǫ.

• Case 2. Suppose now that |σ| < ǫ. Function ξ in (20) can be negative. It means that V would be sign indefinite, and it is not possible to conclude on the closed-loop system stability. Therefore, |σ| can increase over ǫ. As soon as |σ| becomes greater than ǫ, V ≤ -β•V 1/2 and V starts decreasing. Apparently, decrease of V can be achieved via increase of K allowing |σ| to increase before it starts decreasing down to |σ| ≤ ǫ. Without loss of generality, let estimate the overshoot when σ 0 = σ(0) = ǫ + and K 0 = K(0) = K(0) > 0: considering the "worst" case (with respect to uncertainties/perturbations), one has

σ = Ψ M -K • Γ m K = K • |σ| (24)
Then, it yields

σ(t) = σ 0 cos( KΓ m t) + Ψ M -K 0 • Γ m KΓ m • sin( KΓ m t) K(t) = σ 0 K Γ m sin( KΓ m t) + K 0 - Ψ M Γ m cos( KΓ m t) + Ψ M Γ m (25)
Then, one gets

σ(t) = σ 2 0 + (Ψ M -K 0 • Γ m ) 2 KΓ m sin KΓ m t + Θ σ K(t) = σ 2 0 K Γ m + K 0 - Ψ M Γ m 2 sin KΓ m t + Θ K + Ψ M Γ m (26)
It appears from (26) that, when σ 0 = ǫ + → ǫ, the maximum value δ of σ M reads as

δ = ǫ 2 + Ψ 2 M KΓ m (27)
In conclusion, σ converges to the domain |σ| ≤ ǫ in a finite time, but could be sustained in the bigger domain |σ| ≤ δ. Therefore, the real sliding mode exists in the domain |σ| ≤ δ.

On ǫ-tuning

The objective of this section consists in providing a methodology for the tuning of parameter ǫ for the both algorithms in case of sampled controllers. In fact, the choice of parameter ǫ has to be made by an adequate way because a "bad" tuning could provide either unstability and control gain increasing to infinity, or bad accuracy for closed-loop system as described in the sequel.

• If parameter ǫ is too small, and due to large gain K and sampling period T e , system trajectories are such that |σ| never stays lower than ǫ. From K-dynamics (13) or ( 16), it yields that gain K is increasing, which induces larger oscillation, and so on ...

• If parameter ǫ is too large, system trajectories are such that, in spite of large gain K and sampling period T e , |σ| is evolving around ǫ, it follows that controller accuracy is not as good as possible.

Note that ǫ should rather be too large than too small because, in the first case, closed-loop system trajectories are globally uniformly ultimately bounded even if accuracy is bad.

Parameter ǫ has to be tuned such that, as long as K(t) is greater than |Ψ/Γ|, |σ| < ǫ. In this case, given that gain K(t) is sufficient to counteract the perturbations, there is no reason that σ increases over ǫ. Suppose that, at t = t 1 , |σ(t 1 )| ≤ ǫ and

K(t 1 ) ≥ Ψ(t 1 ) Γ(t 1 ) . ( 28
)
It is clear that control gain K is sufficient to ensure

|σ(t 1 + T e )| ≤ ǫ.
From Euler's formula, one has

σ(t 1 + T e ) ∼ σ(t 1 ) + [Ψ(t 1 ) -Γ(t 1 ) • K(t 1 ) • sign(σ(t 1 )] • T e (29)
From (28), one has

|Ψ| ≤ |Γ| • K. It yields |Ψ(t 1 ) -Γ(t 1 ) • K(t 1 ) • sign(σ(t 1 ))| ≤ 2Γ M • K(t 1 ).
Furthermore, as gain K(t) fulfills (28), σ-sign is the opposite of σ sign. It means that, if 0 ≤ σ(t 1 ) ≤ ǫ (resp. -ǫ ≤ σ ≤ 0), function σ will decrease (resp. increase). Then, in order to guarantee that |σ(t 1 + T e )| will not exceed ǫ, the "worst" case is σ(t 1 ) = 0. From (29), it yields

|σ(t 1 + T e )| ≤ 2 Γ M • K(t 1 )T e (30)
Given that it must be ensured |σ(t 1 + T e )| ≤ ǫ, and as the best accuracy and the closed-loop system are the both objectives, ǫ has to be chosen as a time-varying function equal to the smallest admissible value

ǫ(t) = 2 Γ M • K(t 1 )T e . (31)
Discussion. Equation (31) means that the knowledge of bounds of Γ is required for the choice of ǫ in case of sampling controller. However, note that this latter equation gives the best tuning. As said previously, ǫ can be tuned to a "too" large value which guarantees the stability but with a worst accuracy. Furthermore, the equation (31) gives at least a methodology for ǫ-computation: no such information is given in [START_REF] Huang | Adaptive sliding-mode control for nonlinear systems with uncertain parameters[END_REF] whereas a such parameter is required for the boundary layer in case of practical implementation of the controller. However, in case of many practical applications, the knowledge of Γ bounds is not really required. In fact, when applied, the controller design procedure is very often the following:

uncertain functions read as Ψ = Ψ N om + ∆Ψ and Γ = Γ N om + ∆Γ where Ψ N om and Γ N om are the nominal known functions and ∆Ψ and ∆Γ their unknown parts. Furthermore, in most of the cases, one has |Ψ N om | ≥ |∆Ψ| and |Γ N om | ≥ |∆Γ|. Then, by supposing that Ψ N om is invertible, the control law

3 u = -Ψ N om + v Γ N om is applied to system (2) which gives σ = ∆Ψ + Ψ N om Γ N om ∆Γ + 1 - ∆Γ Γ N om • v (32)
As

|Γ N om | ≥ |∆Γ|, one gets |1 - ∆Γ Γ N om | ≤ 2. From (31), it yields ǫ(t) = 4K(t)T e .
5 Simulations

Tutorial

Consider the following uncertain system been made by supposing that σ(0) = 10 and T e = 0.0001 sec. From (31), parameter ǫ reads as ǫ(t) = 2K(t) • T e .

σ = Ψ(t) + u (33)
Control algorithm 1 in eqs. ( 7)-( 13)- [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]. Control input u and gain adaptation law K(t) respectively read as in eqs. ( 7)-( 13)- [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]. Gain is initialized at K(0) = 10, its dynamics being tuned with K1 = 1000, K3 = 1 and τ = 0.1 sec.

Control algorithm 2 in eqs. ( 15)- [START_REF] Levant | Principles of 2-sliding mode design[END_REF]. Control input u and gain adaptation law K(t) respectively read as in eqs. ( 15)- [START_REF] Levant | Principles of 2-sliding mode design[END_REF]. Gain is initialized at K(0) = 10, its dynamics being tuned with K = 1000 and µ = 0.1. 7)-( 13)- [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]. Top-Left. Sliding variable σ versus time (sec). Top-Right. Zoom on sliding variable σ (solid line), ǫ(t) and -ǫ(t) (dotted lines) versus time (sec). Bottom-Left. Control input u versus time (sec). Bottom-Right. Gain K(t) (solid line) and perturbation Ψ(t) (dotted line) versus time (sec).

Figures 2 and3 display simulations results of both previous controllers applied to system (33). It appears that both control laws yield to very similar results. Based on simulation plots, one can conclude that both adaptation algorithmes provide the controller gain K(t) to follow closely the perturbation Ψ(t) which profile and boudary are not known a priori (see bottom-right plots in Figures 2 and3). Therefore, the control gain K(t) is not over-estimated and control chattering is minimal. Also, the sliding variable σ is robustly constrainted to |σ| < ǫ(t) with exception of sparce spikes.

Electropneumatic actuator

The electropneumatic system under interest is a double acting actuator (Figure 4) composed by two chambers, denoted P (as positive) and N (as negative). The air mass flow rates entering the two chambers are modulated by two three-way servodistributors controlled with two electrical inputs of opposite signs (u and -u). The pneumatic jack horizontally moves a load carriage of mass M 15)-( 16). Top-Left. Sliding variable σ versus time (sec). Top-Right. Zoom on sliding variable σ (solid line), ǫ(t) and -ǫ(t) (dotted lines) versus time (sec). Bottom-Left. Control input u versus time (sec). Bottom-Right. Gain K(t) (solid line) and perturbation Ψ(t) (dotted line) versus time (sec).

Model

Following standard assumptions on the pneumatic part of the electropneumatic system [START_REF] Shearer | Study of pneumatic processes in the continuous control of motion with compressed air[END_REF][START_REF] Cloy | Discharge characteristics of servo valve orifices[END_REF][START_REF] Brun | Control of an electropneumatic actuator, comparison between some linear and nonlinear control laws[END_REF][START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF][START_REF] Laghrouche | Higher order sliding mode control based on optimal approach of an electropneumatic actuator[END_REF], one gets a nonlinear dynamic model for the whole system ṗP = krT V P (y) [q mN (u, p p ) -

S P rT p P v] ṗN = krT V N (y) [q mP (-u, p N ) + S N rT p N v] v = 1 M [S P p P -S N p N -bv -F f -F ext ] ẏ = v (34)
with y the load carriage position, v its velocity and p P and p N the pressures of P and N chambers. V X (X = P or N ) is the volume in the chamber X, q mX (u X , p X ) the mass flow rate provided from the servodistributor X, k the polytropic constant, r the perfect gaz constant, T the supply temperature and S X the piston area in the chamber X. The term F f represents all the dry friction forces which act on the moving part in presence of viscous friction (b.v) and an external force only due to atmospheric pressure (F ext ).

The model of mass flow rate delivered by each servodistributor can be reduced to a static function described by two relationships q mP (u, p P ) and q mN (-u, p N ). The two first equations of (34) concern the pneumatic part of the system and are obtained from the state equation of perfect gases, the mass conservation law and the polytropic law under the assumptions given above. The two last equations describe the mechanical part and are derived from the fundamental mechanical equation applied to the moving part. In order to get an affine nonlinear state model, the mass flow rate static characteristic issued from measurements [START_REF] Sesmat | Static characteristics of a three way servovalve[END_REF] reads as [START_REF] Belgharbi | ÒAnalytical model ofthe owstageofapneumatic servo-distributor for simulationandnonlinear controlÓ[END_REF] q mP (u, p P ) = ϕ P (p P ) + ψ P (p P , sgn(u))

• u q mN (-u, p N ) = ϕ N (p N ) -ψ N (p N , sgn(-u)) • u (35)
with functions ϕ P and ψ P (resp. ψ N and ϕ N ) 5 th -order polynomial functions with respect to p P (resp. p N ). Denoting x = [p P p N y v] T , system (34) can be written as nonlinear system (1). Let define X as the physical domain It yields that, for x ∈ X , system dynamics are bounded under a bounded control input u. Two kinds of uncertainties are taken into account: uncertainties due to the identification of physical parameters, and perturbations. Viscous and dry friction coefficients have been identified, but the determination of their variations around their nominal values is a hard task. For example,the dry friction coefficient is difficult to identify because depending on the track surface quality (thus the piston position), the seal wear, the working conditions (temperature, pressure, quality of air) .... The mass flow rate delivered by each servodistributor has been approximated by polynomial functions (35); the uncertainties on ϕ(•) and ψ(•) have been evaluated to ±15% ±5% respectively. Finally, during the load moving, the total mass in displacement can evolve from 17 kg until 47 kg.

Control design and simulations

The aim of the control law is to get a good accuracy in term of position tracking for the desired trajectory defined in Figure 5. Following the previous section, consider the sliding variable

σ = λ 2 • (y -y d (t)) + 2λ • ( ẏ -ẏd (t)) + (ÿ -ÿd (t)) (36)

Conclusion

This paper proposes two new methodologies for adaptive sliding mode controller design. The both algorithms allow establishing of the sliding mode via the sliding mode control laws with gain adaptation without a priori knowing uncertainties/perturbations bounds while the adaptive gains values are both not over-estimated. The first algorithm is based on evaluation of uncertainties/perturbation by using equivalent control concept that requires employment of low-pass filter. The second adaptive control law does not estimate the boundary of pertur-bations/uncertainties and yields establishing a real sliding mode. The efficacy of these new strategies has been confirmed on a tutorial example as well as by controlling the electropneumatic actuator. Future works will concentrate on extending the developed methodologies to both MIMO uncertain nonlinear systems and to systems with higher order adaptive-gain sliding mode control, as well as on experimental validating the advanced algorithms.

Note that this gain has a finite value by the absolute continuity property of K(t). From t = t 1 , given K-dynamics, gain K is large enough to make the sliding variable σ decreasing. Then, it yields that, in a finite time t 2 , σ = 0 and K(t 2 ) admits a bounded value. It yields that there always exists a positive constant K * such that K * > K(t), for all t ≥ 0.

Proof of Lemma 2. Suppose that |σ(x, t)| > ǫ. From K-dynamics, and given that functions Ψ and Γ are supposed bounded, it follows that K is increasing and there exists a time t 1 (see Figure 8) such that

K(t 1 ) = | Ψ(t 1 ) Γ(t 1 )
|.

From t = t 1 , given K-dynamics, gain K is large enough to make the sliding variable σ decreasing. Then, it yields that, in a finite time t 2 (see Figure 8), |σ| < ǫ. It yields that gain K is decreasing from t 2 , gain K being at a maximum value at t = t 2 . From K-dynamics, it yields that there exists a time instant t 3 > t 2 (see Figure 8) such that

K(t 3 ) = | Ψ(t 3 ) Γ(t 3 ) |.
From t = t 3 , gain K is not large enough to counteract perturbations and uncertainties as it is decreasing. It yields that there exists a time instant t 4 > t 3 such that |σ(t 4 )| > ǫ. The process then restarts from the beginning. By the assumptions (3), the gains K(t i ) remain bounded uniformly on t i . In fact,

K(t i ) = | Ψ(t i ) Γ(t i ) | ≤ Ψ M Γ m := K * *
and, hence, there always exists a finite constant K * such that

K * ≤ K * *
which proves the desired result.

Function Ψ(t) is an uncertain bounded function described by Figure 1 Figure 1 :

 11 Figure 1: Unknown term Ψ versus time (sec).

Figure 2 :

 2 Figure 2: Control algorithm 1(7)-(13)-[START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]. Top-Left. Sliding variable σ versus time (sec). Top-Right. Zoom on sliding variable σ (solid line), ǫ(t) and -ǫ(t) (dotted lines) versus time (sec). Bottom-Left. Control input u versus time (sec). Bottom-Right. Gain K(t) (solid line) and perturbation Ψ(t) (dotted line) versus time (sec).

Figure 3 :

 3 Figure 3: Control algorithm 2 (15)-[START_REF] Levant | Principles of 2-sliding mode design[END_REF]. Top-Left. Sliding variable σ versus time (sec). Top-Right. Zoom on sliding variable σ (solid line), ǫ(t) and -ǫ(t) (dotted lines) versus time (sec). Bottom-Left. Control input u versus time (sec). Bottom-Right. Gain K(t) (solid line) and perturbation Ψ(t) (dotted line) versus time (sec).

Figure 4 :

 4 Figure 4: Electropneumatic system

X

  = x | 1 bar ≤ p P ≤ 7 bar, 1 bar ≤ p N ≤ 7 bar, -200 mm ≤ y ≤ 200 mm, |v| ≤ 1 m.sec -1 .

Figure 6 :

 6 Figure 6: CONTROL ALGORITHM 1. Top. Current (solid line) and desired (dashed line) position trajectories (m) versus time (sec). Center. Control input u (V ) versus time (sec). Bottom. Control gain K versus time (sec).

Figure 7 :

 7 Figure 7: CONTROL ALGORITHM 2. Top. Current (solid line) and desired (dashed line) position trajectories (m) versus time (sec). Center. Control input u (V ) versus time (sec). Bottom. Control gain K versus time (sec).

It yields that t * defined just after (9) equals first-of-all t1, then t2. Furthermore, in implementation case, t * is on-line determined by comparing, at each time t, σ(x(t -), t -) and σ(x(t), t) = 0.

The parameter ǫ can be time-varying, as in the next section. Its tuning is detailed in Section 4.

In[START_REF] Castro-Linares | High order sliding mode observer-based control[END_REF], it has been shown that, for such class of uncertainties, this control law allows to decrease the influence of perturbations and uncertainties, and then to reduce the control gain.

In[START_REF] Brun | Control of an electropneumatic actuator, comparison between some linear and nonlinear control laws[END_REF], stability of zero-dynamics has been numerically established.

It can be shown that ΓNom is invertible in the physical domain of the application[START_REF] Laghrouche | Higher order sliding mode control based on optimal approach of an electropneumatic actuator[END_REF][START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF].

This value and the following represent parameters variations with respect to their nominal values.

This work has been made while Franck Plestan was on sabbatical at UAH. This stay has been

with λ > 0. As system (34)-(35) admits a relative degree equal to 1 for σ with respect to u, one gets

As described in [START_REF] Laghrouche | Higher order sliding mode control based on optimal approach of an electropneumatic actuator[END_REF][START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF], functions Ψ(•) and Γ(•) can be written as the sum of a "nominal" part (Ψ N om and Γ N om ) and an "uncertain" one (∆Ψ and ∆Γ), i.e.

From [START_REF] Laghrouche | Higher order sliding mode control based on optimal approach of an electropneumatic actuator[END_REF][START_REF] Girin | Robust control of an electropneumatic actuator: application to an aeronautical benchmark[END_REF], functions Ψ and Γ are bounded for all x ∈ X ; furthermore, for all x ∈ X , Γ N om > 0. Then, Ψ and Γ fulfill (3). It yields that the control law reads as

with control input v reading as previous control algorithms. Parameter ǫ reads as ǫ(t) = 4K(t) • T e . Parameters have been tuned as

Furthermore, for Algorithm 1, one has

For Algorithm 2, one has µ = 0.1. Following simulations have been made by making variations of load mass (+20%) 6 and uncertainties on ϕ and ψ (-20 %). Figures 67display, for the both control laws, the actuator position with respect to desired trajectory (Top), the control input u (Center), and the gain K (Bottom). It appears that the both strategies yield to quite similar results. Note that, with the proposed gain tuning, Algorithm 1 induces lower magnitudes of control and gain in this time interval; in fact, Algorithm 2 gain is increasing to a larger value and takes time to decrease to similar values of Algorithm 1. However, to authors' experience, implementation of Algorithm 2 is clearly more easy.

Appendix

Proof of Lemma 1. Suppose that σ(x, t) = 0. From K-dynamics (8), and given that functions Ψ and Γ are supposed bounded, it follows that K is increasing and there exists a time t 1 such that


Frédéric Mora

Lilian Aveneau

Livia Oana

Djamchid Apostu

Ghazanfarpour

O Apostu

D Ghazanfarpour

Lazy visibility evaluation for exact soft shadows

come

Lazy visibility evaluation for exact soft shadows

Introduction

Shadows and especially soft shadows play an important role in visual perception as they unveil spatial relationships between objects. This makes soft shadows essential to realistic rendering. However, they are very expensive to compute since they require solving the visibility of an area light source from any visible point in a scene. The accuracy of any solution to this problem determines the soft shadows quality and, as a consequence, the realism of the scene. Dedicated algorithms were developed to reduce the computational complexity inherent to soft shadows. Interactive solutions exist, but they make crude assumptions leading to approximate or plausible soft shadows. They are forced to sacrifice the correctness of the shadows in order to reach interactive frame rates. Geometry-based approaches can lead to fine results, using silhouette edges [START_REF] Lehtinen | An improved physically-based soft shadow volume algorithm[END_REF][START_REF] Forest | Accurate shadows by depth complexity sampling[END_REF] for example. While the accuracy of the visibility computation is crucial for soft shadows, the wide majority of algorithms are based on sampled light sources. As a consequence, whatever their accuracy, they are still sensitive to noise according to the sampling strategy and density.

In this paper, we present a new high quality soft shadow algorithm based on exact visibility computed from the area light source. Our method relies on an equivalence relation allowing to group lines according to the geometry they intersect [START_REF] Pellegrini | Handbook of Discrete and Computanional Geometry -second edition[END_REF]. We build on this theoretical framework to provide a novel and practical approach to lazily compute the exact visibility from a polygon. This visibility data allows to solve efficiently and exactly the visibility of an area light source for each point to be shaded, taking advantage of the visibility coherence between neighbour points. Since the direct illumination is analytically evaluated using exact visibility, the shadows are optimal and alias-free. While exact approaches suffer from implementation and robustness issues, our method is easy to implement, scalable and numerically robust. In addition, our results show that it outperforms an optimized ray tracer.

The paper is organized as follows: The first section outlines the works related to this paper. The second section presents the geometrical and mathematical framework we use. Moreover, we demonstrate an 1 essential theorem to our algorithm. The third section details our exact from-polygon visibility algorithm which is the basis to solve the visibility of an area light source. At last, in the section 4, we test our algorithm's behaviour in different configurations. Before concluding, we discuss these results, we analyse some limitations and we consider several perspectives.

Previous works

Many contributions on soft shadows can be found in the literature. A survey has been presented by Hasenfratz et al. [START_REF] Lapierre | A survey of real-time soft shadows algorithms[END_REF]. Recently, Eisemann et al. [START_REF] Eisemann | Casting shadows in real time[END_REF] provide another interesting lecture. In this section, we only focus on widely used or recent physically-based algorithms able to compute very high quality or exact soft shadows. In addition, since the core of our method relies on an exact from-polygon visibility algorithm, we also focus on this topic.

Soft shadows in ray and beam tracing

Ray tracing easily supports soft shadows using shadow rays and stochastic sampling of the light sources. As a counterpart, the results are sensitive to noise. This requires increasing the number of samples in order to avoid visual artifacts. But this drastically slows down the rendering. Packet [START_REF] Benthin | Interactive rendering with coherent ray tracing[END_REF] and frustrum [RSH05, WIK * 06] traversal techniques improve ray-tracing efficiency. They take advantage of the spatial coherence and have efficient implementations, using SIMD instructions. Moreover, they can be designed for a specific graphic architecture [START_REF] Benthin | Efficient ray traced soft shadows using multi-frusta tracing[END_REF]. Despite these improvements, soft shadows remain a challenging task. If the sample distribution can be improved using importance or interleaved sampling, the results remain sensitive to noise.

Beam tracing [START_REF] Heckbert | Beam tracing polygonal objects[END_REF] methods are suited for soft shadows computation. A single beam is cast to the light to solve exactly its visibility. As a consequence, results are not sensitive to noise. Beam tracing is usually considered as lacking robustness, which can be improved by adaptive methods [START_REF] Ghazanfarpour | A beam tracing with precise antialiasing for polyhedral scenes[END_REF]. Overbeck et al. [START_REF] Overbeck R | A real-time beam tracer with application to exact soft shadows[END_REF] propose a robust beam-triangle intersection to handle more complex models with soft shadows.

Packets, frustrum or beam techniques take advantage of the spatial coherence from each point to be shaded. Our algorithm uses the spatial coherence from the light source. This is an important difference since the shading of the points is not an independent operation anymore.

Silhouette based soft shadows

Möller and Assarson [START_REF] Akenine-Möller T | Approximate soft shadows on arbitrary surfaces using penumbra wedges[END_REF][START_REF] Assarsson | A geometry-based soft shadow volume algorithm using graphics hardware[END_REF] propose the soft shadow volume algorithm on graphics hardware. It extends the shadow volume algorithm [START_REF] Crow | Shadow algorithms for computer graphics[END_REF] using penumbra wedges. A wedge is the place where the light source view can change because of a silhouette edge. Given a point inside a wedge, the relevant silhouette edge is projected onto the light to determine the light samples visibility. However, overlapping silhouette edges lead to overestimating the shadow and silhouette edges are underestimated since they are computed only from the center of the light source. As a consequence, soft shadows are approximated.

The soft shadow volume algorithm is extended to offline ray-tracing by Laine et al. [LAA * 05]. All the silhouette edges are taken into account and the algorithm deals properly with the overlapping ones. Special rules determine the light samples affected by each silhouette edges. Then a single ray is shot to recover the samples visibility. Lehtinen et al. [START_REF] Lehtinen | An improved physically-based soft shadow volume algorithm[END_REF] remark that this algorithm can be over-conservative and less efficient than a classical ray-tracer. They change the data structure used to store the silhouette edges, improving both the algorithm efficiency and its memory consumption.

Forest et al. [START_REF] Forest | Accurate shadows by depth complexity sampling[END_REF] build on the works of [START_REF] Assarsson | A geometry-based soft shadow volume algorithm using graphics hardware[END_REF] and [START_REF] Lehtinen | An improved physically-based soft shadow volume algorithm[END_REF]. They remove the shooting ray step and provide a dedicated implementation on graphics hardware allowing several frames per second on moderate scenes using an average of 16 light samples per light sources.

Research on silhouette based algorithms has lead to very efficient solutions. However, they require triangles connectivity with consistent orientation. Since these methods are object-based, visibility computations are accurate. But they are still subject to aliasing because of the light source sampling. The sampling strategy has a significant impact as shown in [START_REF] Forest | Accurate shadows by depth complexity sampling[END_REF].

Other methods

Laine et al. [START_REF]Hierarchical penumbra casting[END_REF] change the pre-processing order of ray tracing. Instead of finding a triangle hit by a shadow ray, they search all the shadow rays intersecting a given triangle. Notice that the same idea was used to provide interactive GPU algorithms based on triangles [START_REF] Eisemann | Visibility sampling on gpu and applications[END_REF][START_REF] Sintorn | Sample-based visibility for soft shadows using alias-free shadow maps[END_REF] or extended to silhouette edges [JHH * 09]. However, approximations are made to reach such frame rates. In this work, we also focus on the sets of lines intersecting triangles. However, we use an analytical representation of those sets.

Exact from-polygon visibility

Exact solutions to visibility problems are important to better understand its underlying nature, and design new algorithms. However, exact from-polygon (and thus from-region) visibility is very complex, because it is a four-dimensional problem [START_REF] Drettakis | The 3d visibility complex[END_REF]. There are few solutions in the literature. They all rely on the same theoretical framework described by Pellegrini [START_REF] Pellegrini | Handbook of Discrete and Computanional Geometry -second edition[END_REF] in the Plücker space, a 5D space of lines. His analysis provides a worst case complexity of O(n 4 log n) in space and time, where n is the number of triangles. This underlines the complexity inherent to the visibility in space. To remain practicable, a very special care must be given to the algorithm design. As an example, a first approach is proposed by Mount and Pu [START_REF] Mount | Binary space partitions in plücker space[END_REF]. But their experiments do not exceed 14 random triangles, because, as stated by the authors, they reach the theoretical complexity upper bound.

The two first practicable algorithms are provided by Nirenstein et al. [START_REF] Nirenstein | Exact from-region visibility culling[END_REF] and Bittner [START_REF] Bittner | Hierarchical Techniques for Visibility Computations[END_REF] (see [START_REF] Aveneau | Coherent and exact polygon-to-polygon visibility[END_REF] for details about the differences between the two solutions), and further improved by Haumont et al. [START_REF] Haumont | A low dimensional framework for exact polygon-to-polygon occlusion queries[END_REF] and Mora et al. [START_REF] Aveneau | Fast and exact direct illumination[END_REF] respectively. It is worth underlying a major difference between the two approaches initiated by Nirenstein and Bittner respectively :

• Nirenstein and Haumont design their algorithms to prove that two polygons are visible. Thus the algorithm output is a boolean value. Whenever possible, the computation of the whole visibility set is avoided since the process can be terminated as soon as a visibility is found.

• Mount, Bittner and Mora algorithms aim to compute and store all the visibility information. In this case, the output is a partition of the Plücker space encoded in a Binary Space Partitioning tree.

Despite this difference, they all rely on Computing Solid Geometry in the Plücker space which is computationally very expensive. As a consequence, they are only used as a pre-process step. In addition 5D CSG operations are very complex to implement and prone to numerical instability. Thus, the reliability and the scalability of these solutions is restricted.

The soft shadow algorithm presented in this paper relies on a novel approach to capture the visibility from a polygon. It also leads to the construction of a BSP tree in the Plücker space. But in contrast to the previous works, the tree is built lazily and at run time, not as a pre-process step. Visibility is computed on-demand when and where it is required, according to the image resolution. It does not rely on any expensive and complicated 5D CSG operation. This makes the algorithm very easy to implement, efficient and computationally robust, whereas previous methods are complex and suffer from numerical instability.

Geometrical basis

In this section we introduce the geometrical knowledge underlying the work presented in this paper. In particular, we demonstrate a result on the orientation of the lines stabbing two polygons. This result is important for our algorithm.

Plücker's coordinates

The Plücker space is a five dimensional projective space, denoted P 5 , well known in computer graphics as an efficient solution for dealing with real 3D lines [START_REF] Shoemake | Plücker coordinates tutorial[END_REF]. We only recall here the properties used in this paper. Let's consider an oriented 3D line l going through two distinct points P and Q of coordinates (p x , p y , p z ) and (q x , q y , q z ) respectively. The line l maps to the Plücker point denoted π l ∈ P 5 . The six The relative orientation of two lines is given by the sign of the side operator. Under: Since side(l, r) = h l (π r ), applying the side operator comes down to testing the position of the Plücker point of one line against the dual hyperplane of the other line coordinates of π l , denoted (l 0 , l 1 , l 2 , l 3 , l 4 , l 5 ) are defined by: l 0 = q xp x l 3 = q z p yq y p z l 1 = q yp y l 4 = q x p zq z p x l 2 = q zp z l 5 = q y p xq x p y Any Plücker point is in bijection with its dual hyperplane:

h l (x) = l 3 x 0 + l 4 x 1 + l 5 x 2 + l 0 x 3 + l 1 x 4 + l 2 x 5 = 0
with x ∈ P 5 with coordinates (x 0 , x 1 , x 2 , x 3 , x 4 , x 5 ). Thus, any line l in the 3D space can be mapped to the Plücker space as a point π l or its dual hyperplane h l . We now recall the definition of the so-called side operator: side(l, r) = l 3 r 0 + l 4 r 1 + l 5 r 2 + l 0 r 3 + l 1 r 4 + l 2 r 5 where l and r are two 3D lines, (l 0 , l 1 , l 2 , l 3 , l 4 , l 5 ) and (r 0 , r 1 , r 2 , r 3 , r 4 , r 5 ) are respectively the coordinates of their Plücker point, π l and π r . The sign of the side operator dertermines the relative orientation of the two lines. In particular, two lines are incident (or parallel) if their side operator equals zero. We can notice that: side(l, r) = h l (π r ). This leads to the geometrical interpretation of the side operator as depicted in Figure 1.

On the lines intersecting the same triangles

A direct application of the side operator is an easy and robust line-triangle intersection test. A line intersects a triangle if its orientation is consistent in respect to the lines spanning the triangle edges. This is depicted in Figure 2. Beyond this intersection test, we can notice that the 3 hyperplanes related to the triangle edges

a πa l 0 h l 0 l 1 h l 1 l 2 h l 2 R 3 P 5
Figure 2: Left: a triangle in 3D space. l 0 , l 1 , l 2 are the lines spanning the triangle's edges. a is a line stabbing the triangle. Right: h l 0 , h l 1 and h l 2 are the Plücker hyperplanes mapped from l 0 , l 1 and l 2 . π a is the Plücker point mapped from the stabbing line a. π a has a consistent orientation in respect to l 0 , l 1 and l 2 . From a geometrical point of view, π a lies at the intersection of the halfpsaces induced by h l 0 , h l 1 and h l 2 . Thus, these 3 hyperplanes provide an analytical representation of all the lines stabbing the triangle.

divide the Plücker space into cells and one of them contains all the lines (i.e. their Plücker points) stabbing the triangle, while the other cells contain the lines missing the triangle.

Pellegrini [START_REF] Pellegrini | Handbook of Discrete and Computanional Geometry -second edition[END_REF] develops a more general approach. He uses the Plücker space as a framework to provide theoretical bounds on various problems involving lines. Let S be a set of several triangles (or convex polygons) and L S be the lines spanning the triangles edges in S. Each line in L S can be mapped to Plücker space as a hyperplane. This builds an arrangement in Plücker space: A decomposition of the space into cells by a set of hyperplanes. All the points in a same cell satisfy the following property: They all have the same sign in respect to its bounding hyperplanes. This is illustrated in Figure 3 with 2 triangles. Thus, in Plücker space, all the lines (i.e. their Plücker points) belonging to the same cell intersect the same subset of triangles in S. Notice that this subset can be empty if the lines miss all the triangles in S.

The decomposition of Plücker space into cells allows to group lines together according to the subset of triangles they intersect. This defines an equivalence relation on lines. As a consequence, each cell corresponds to an equivalence class (sometimes called isotopy class [START_REF] Pellegrini | Ray-shooting and isotopy classes of lines in 3-dimensional space[END_REF] or orientation class [CEG * 96]). To sum up, the Pellegrini approach allows an exact and analytical representation of all the sets of lines generated by a set of triangles S, using the set of lines L S .

Orientation of lines intersecting two polygons

We focus on the set of lines intersecting two convex polygons. We demonstrate the following theorem:

Theorem 1 Let A and B be two convex polygons with n and m vertices respectively. We define vvset

= {v i j , i ∈ [1, n], j ∈ [1, m]}

the set of the lines v i j defined by one vertex of A and one vertex of B. Let l be any line and q a line intersecting A and B:

side(l, x) ≥ 0, ∀x ∈ vvset ⇒ side(l, q) ≥ 0 side(l, x) ≤ 0, ∀x ∈ vvset ⇒ side(l, q) ≤ 0
In other words, if all the vertex-to-vertex lines of two polygons have a positive (resp. negative) orientation in respect to any line l, all the stabbing lines of the two polygons will have a positive (resp. negative) orientation in respect to l. To demonstrate this theorem, we use a vector notation for Plücker's coordinates: Given two points P and Q, the Plücker coordinates of (PQ) is the vector couple ( u, v) such as:

(PQ)    u = -→ PQ = (l 0 , l 1 , l 2 ) v = -→ OP × -→ OQ = (l 3 , l 4 , l 5 )
Where O is the origin in 3D space. For simplicity, we will omit the symbol , simply denoting the vector couple (u, v). Using this notation, we now formulate the side operator as the sum of two dot products. For example, given two lines l and l ′ with Plücker coordinates (u, v) and (u ′ , v ′ ) respectively, then : side(l, l ′ ) = u.v ′ + v.u ′ . Proof: Let A and B be two convex polygons with n and m vertices in the 3D space. Using barycentric coordinates, any point a on A (resp. b on B) can be defined by a combination of its vertices (a 1 , ...,

a n ) (resp. b 1 , ..., b m ): a = n ∑ i=1 α i a i with n ∑ i=1 α i = 1, with ∀i α i ≥ 0 b = n ∑ i=1 β i b i with n ∑ i=1 β i = 1, with ∀i β i ≥ 0
a and b define any line q stabbing A and B, and the Plücker point of q is (ab, oa × ob). Thus the relative orientation of a line q and any line l (u, v) is:

side(q, l) = u.(oa × ob) + v.ab
Since the orientation of two lines is relative, we can assume without loss of generality that l goes through the origin of the 3D space (otherwise it is always possible to translate l). Thus v (the last 3 Plücker coordinates of l) becomes null and:

side(q, l) = u.(oa × ob) = u.(∑ n i α i oa i × ∑ m j β j ob j )
Since the cross product is distributive over addition:

side(q, l) = u.(∑ m j ∑ n i (α i oa i × β j ob j )) 6
Next, by compatibility with scalar multiplication:

side(q, l) = u.(∑ m j β j ∑ n i α i (oa i × ob j ))
At last, since the dot product is distributive over addition:

side(q, l) = ∑ m j β j ∑ n i α i u.(oa i × ob j ) = ∑ m j β j ∑ n i α i (side(a i b j , l))
Since α i ≥ 0 and β j ≥ 0, the sign of side(q, l) depends only on the sign of side(a i b j , l). In particular, if they are all positive (resp negative) then side(q, l) is positive (resp negative).

Using this theorem, we can determine if the set of lines stabbing two polygons has a consistent orientation in respect to a given line. The visibilty algorithm presented in the next section uses this essential result.

It can be proved that the Plücker points of the lines in vvset are the vertices of the smallest convex polyhedron in the Plücker space containing the stabbing lines of A and B. A demonstration can be found in [START_REF] Charneau | Plücker space and polygon to polygon visibility computation with Geometric Algebra[END_REF][START_REF] Aveneau | Guide to geometric algebra in practice[END_REF], but it requires advanced knowledge in Geometric Algebras. For the understanding and the correctness of this work, the property demonstrated in this paper is sufficient.

Algorithm

Since equivalence classes are continuous sets of lines that hit/miss the same subset of triangles, they represent coherent paths through the scene independently of any viewpoint. As a consequence, two lines belonging to the same equivalence class are spatially coherent. We use this property to build a Plücker space partition representing the visibility of an area light source.

Overview

We consider an area light source L, a triangle T , and define their occluders O as the triangles intersecting their convex hull. To represent the light source visibility, we focus on the sets of lines intersecting L and which either miss all the occluders in O, or hit at least one occluder in O. Thus, we define:

• A visible class: Any equivalence class representing a set of lines that do not intersect any occluders.

• An invisible class: Any equivalence class representing a set of lines that intersect at least one occluder.

• An undefined class: An equivalence class that is not yet found as visible or invisible.

Our algorithm builds a BSP tree in Plücker space, providing a hierarchical representation of the equivalence classes generated by the occluders. Each leaf represents one of these three classes. The algorithm is lazy: The BSP tree is grown on-demand depending on when and where visibility information is needed. The construction only relies on two operations: Inserting and merging an occluder into the tree. The former allows to find the leaves affected by an occluder in the tree, while the latter grows the tree by replacing a leaf with an occluder's equivalence classes.

Merging an occluder:

Let o be an occluder, we note bsp(o) the BSP representation of the equivalence classes generated by o (see Figure 4). The tree grows each time an occluder is merged: The visible/invisible classes generated by the occluder are added to the tree. Merging o into the tree consists in replacing a leaf by the root of bsp(o). As a consequence, each inner node contains a hyperplane corresponding to an occluder's edge.

Inserting an occluder:

Occluders are inserted into the tree and located into the leaves (and thus the classes) they may affect. Inserting an occluder o comes down to testing the relative orientation of a hyperplane (from an inner node) and the lines occluded by o. Since the occlusion created by an occluder o is the set of lines intersecting both L and o, we use Theorem 1 to determine the orientation of the occluded lines with respect to the hyperplane. Figure 5 details the three cases that occur. The BSP tree construction is driven by the visibility queries in order to compute only the required equivalence classes. The following section presents this process and gives the details of our visibility algorithm. Next, we describe how it is used in our soft shadow framework.

l 0 l 1 l 2 h 0 h 1 h 2 o

Visibility algorithm

We The first one with a positive orientation with respect to h, the second one with a negative orientation. Figure 6 gives an illustration. Notice that the orientation of p is coherent with the orientation of h. As a consequence, if L is split by p, we can compute the two relevant polygons L ∩ p + and L ∩ p -so that they represent respectively the positive and negative subset of lines with respect to h. The visibility algorithm (Algorithm 1) finds the equivalence classes related to the visibility of L from xyz. It subdivides L into several convex fragments so that each one of them belongs to a single equivalence class. It starts with all the occluders associated to a single undefined leaf. For each inner node, we compute the plane defined by xyz and the line stored in the node (line 18). Then, L is tested against this plane to determine the orientation of the lines stabbing L and xyz (line 19). If L lies in the positive (resp. negative) half-space of the plane, then all the lines have a positive (resp. negative) orientation, and the algorithm continues in the left (resp. right) subtree (line 20 or 22). Otherwise, L is split against the plane and the algorithm continues recursively in both subtrees with the relevant parts of L (line 25). When a fragment reaches a leaf, two alternatives occur:

• The leaf has no occluders, therefore it is either a visible or an invisible class. In the former case, the fragment is a convex part of L which is visible from xyz, thus it is returned (line 29). Otherwise, the fragment is invisible from xyz and it is discarded. • The leaf has some occluders, the class is undefined and we cannot answer the query without further developing the tree.

In the latter case, the algorithm chooses a random occluder (line 31) ro among the occluders associated with the current leaf. This occluder is merged into the tree by replacing the leaf by bsp(ro), the BSP representation of the four classes generated by ro (line 32). Next, the remaining occluders are inserted into bsp(ro) which is now a new part of the tree. The procedure at line 1 achieve an occluder insertion as describe in Figure 5. When all the occluders have been inserted in bsp(o), the algorithm continues from the root of bsp(o) until it finds the equivalence class for the current light fragment.

Key points

This section underlines significant points of the algorithm.

Random selection of the occluders:

The algorithm's efficiency is related to the balance of the tree. To develop the BSP tree, the algorithm chooses an occluder randomly. Obviously, some choices may lead to a more balanced tree than others. However this is not predictable. In fact, we have tested different heuristics, trying to make a "good choice". Unfortunately, all of them achieved poor improvements compared to the extra computational cost. Moreover, their behaviour can be very different according to the nature of the scene. We have reached the conclusion that a random choice gives better results, and more important, it has a consistent behaviour whatever the rendered scene. This is similar to the choice of the pivot in the well known quicksort algorithm: Although a random pivot is not the optimal choice, it leads to the best performance in practice whatever the entry is.

Conservative insertion:

If all the lines occluded by an occluder o do not have the same orientation in respect to a hyperplane, o is inserted in both subtrees of the relevant node. As a consequence the insertion is conservative. This step could be computed exactly as in [START_REF] Haumont | A low dimensional framework for exact polygon-to-polygon occlusion queries[END_REF][START_REF] Aveneau | Fast and exact direct illumination[END_REF]. But, as explained in section 2, this is expensive and prone to numerical errors. Our algorithm avoids these problems to remain simple and robust.

Tree growth :

Notice that not all the queries will develop the BSP tree. We can expect some queries to develop new equivalence classes, at least at the begining since the BSP tree is empty. And we can expect queries to take advantage of the previous computations thanks to the visibility coherence. This is a key point of the algorithm's efficiency.

Soft shadows framework

To illustrate the efficiency and the reliability of our exact visibility algorithm, we plug it into a ray-tracing rendering software for computing high quality soft shadows. Algorithm 2 describes the process.

Algorithm 2

The following pseudocode illustrates how our visibility algorithm is plugged in a ray tracer software to analytically compute soft shadows end while 14: end for • Using the primary rays, all the image points are grouped together with respect to the triangle they belong to. This builds a list of visible triangles (line 1).

• Multiple lights are handled successively and independently (line 2)

• For each visible triangle, an empty BSP tree is created and associated with its set of occluders (line 8).

• Next, for each image point, the algorithm 1 is used to compute the visible parts of the light (line 10).

In our framework, we consider area light sources with a uniform emission, therefore we compute direct illumination analytically (line 11) by integrating over the visible parts of the light [START_REF] Nishita | Continuous tone representation of three-dimensional objects taking account of shadows and interreflection[END_REF].

This framework is designed to allow an efficient implementation. Since one BSP tree is built per light and per visible triangle, and since visible triangles are successively handled, BSP trees are developed successively and independently. This improves the memory coherence and avoids switching many times between BSP trees. This also limits the memory consumption since each BSP tree is deleted as soon as all its related image points are shaded. Moreover, the implementation can be easily multithreaded: A thread gets a visible triangle from the list, shades its image points and starts over until the list becomes empty. In this case, the triangle list access (line 6) has to be protected.

Implementation details

In Section 4, occluders are defined as the geometry intersecting the convex hull of a triangle and a light source. In practice, the occluders selection relies on a shaft culling approach as described in [START_REF] Haines | A shaft culling tool[END_REF]. Let's consider the bounding box of a visible triangle and the bounding box of an area light source. The shaft is defined as the convex hull of the two bounding boxes plus the triangle's plane support and the light's plane support. Any triangle intersecting the shaft is considered as a possible occluder. While this definition can lead to a conservative occluder set, it can be computed efficiently. At last, compared to Algorithm 1, a stack is managed to avoid recursive system calls. Except this detail, the implementation is straightforward, it uses single floating point precision and does not use SIMD instructions.

Results

All tests were run on a 2.67 GHz Intel Core i7 920 processor with 4GB of memory. For comparison purpose, all pictures were rendered at 1280x720 pixels with one primary ray per pixel. Four sets of results are presented, three sets testing the global performance of our soft shadow framework, and a last one giving an insight on Algorithm 1 behaviour.

Comparisons on time and quality

The first set of results compares our method to ray-traced soft shadows both at comparable time and comparable quality. The ray tracer implementation is similar to [START_REF] Benthin | Interactive rendering with coherent ray tracing[END_REF]: It relies on an optimized SAH kd-tree, uses SIMD instructions to trace four rays simultaneously, and supports multithreading to render several parts of the picture in parallel. Ray-traced soft shadows are computed using groups of 4 shadow rays, and an uncorrelated stratified sampling of the area light sources. Since both our method and the ray-tracer support multithreading, all computations are run using 4 threads. We use four scenes to test our method in different configurations. Despite its moderate geometrical complexity (26 673 triangles), the T-Rex scene is challenging for our approach because it presents difficult and complex shadows due to a long rectangular light source. This means that the light source visibility is complex and this is precisely what is computed by our algorithm. The modified Sponza Atrium with the statue of Neptune (115 737 triangles) and Conference (282 873 triangles) are significant and detailed models with different kinds of shadow complexities. At last, Soda Hall (2 147K triangles) is used to test the scalability of our approach on a massive model with heterogeneous geometry.

Figure 7 details these results, and presents the time spent in soft shadows computation for each method. At comparable time, ray traced soft shadows are always noisy. At comparable quality (i.e. the noise is not noticeable anymore in the stochastic shadows), our method is always faster. In addition, this is very noticeable on a complex case such as the T-rex scene, because the stochastic approach requires a very high number of shadow rays to remove almost all the visible noise. In contrast, the exact visibility algorithm produces high quality results in a few seconds.

Figure 8 presents the computation measurements (memory and time consumption) for our exact visibility algorithm. Since we compute a BSP tree in Plücker space for each visible triangle from the camera, the memory consumption varies during the process according to the building of the BSP tree. As a consequence, our results report the maximum memory load reached by our soft shadow framework. In any case, the memory footprint is low. The lazy evaluation of the visibility driven by the visibility queries allows to focus on the equivalence classes related to the shadows and avoid the computation of useless visibility data. In addition, our framework is designed to avoid building too many BSP trees simultaneously (at most one per thread) which helps to keep the memory low. The total time can be subdived in three steps: The occluder selection, the BSP tree initialization (including the computation of the hyperplanes and the vvset for each occluder), and the visibility queries used for shading the image points. We can notice that the computation time is clearly dominated by the visibility queries i.e. calls to the BSP_query function (see Algorithm 1) which is the core of our method.

Aside the comparison with the stochastic approach, these results demonstrate the robustness of our analytical visibility algorithm in many different configurations. As expected, the sensitivity of our algorithm to the visual complexity of the light source is confirmed. For instance, the T-Rex scene required roughly the same amout of time than the Sponza and Neptune scene or than the Conference scene.

About visibility coherence

A crucial property of our approach is its ability to take advantage of the visibility coherence between image points. To test this ability, the same pictures were rendered again using a modified version of our visibility algorithm: Between each visibility query (i.e. between each BSP_query call), the related BSP tree is reset to its root node associated with its former occluder set (this additional operation is excluded from the timings). In such a case, each query is "the first one" and we prevent all queries from taking advantage of the previous ones. Figure 8 presents (two last columns of the table) the computation times compared to the timings obtained with the true version: The loss of efficiency is considerable. This demonstrates the capacity of the visibility algorithm to benefit from the visibility coherence, which is a key point to its efficiency.

Increasing the area of the light

The soft shadows complexity also depends on the size of the area light source. This is intrinsic to the soft shadows problem and it will inevitably affect our algorithm since it relies on the visibility coherence of the light source. If its area increases, the visibility coherence may decrease and lead to a loss of efficiency. As a consequence the second part of our tests investigates our framework's behavior when the area light source increases. These tests are run on the Conference model, selected for its significant number of triangles as well as the wide range of shadows cast in the scene. The tests start using a small area light source whose size is progressively increased until it becomes 100 times larger. For each size of light, the time and maximum memory used by our method is measured. Figure 9 sums up the results. As expected, it shows a loss of efficiency, both in time and memory, as the size of light source grows (left and right graph of Figure 9). By extent, there is inevitably a critical light size where the time and the memory consumption would become a problem, in particular the memory since it is a limited resource. However, our tests show that, even with the largest area light source, we are far from such a point. In addition, our approach remains competitive with a stochastic solution: At comparable quality, 512 samples are required for the largest light and the computation takes 91 seconds against 39 seconds using our framework (2.33 times slower).

Increasing the light number

Generally, a scene has several light sources. As detailed by Algorithm 2, our implementation supports multiple lights which are handled successively. Thus, it is interesting to test our framework behaviour in such a case. The Conference model is used again and rendered with 2 to 36 area light sources. All the lights have the same size and cast roughly the same "amount" of shadow. Figure 10 graph shows that the time consumption is linear in respect to the number of area light sources. This is the expected behaviour since our framework evaluates the contribution of each light source one after another. In addition this makes the memory load independent from the number of lights, as shown on the right graph. However we can notice a significant growth of the memory load when the number of lights is more than 12. Indeed, the 12 first lights are above the conference table while the 24 other lights are mainly above the chairs, casting more complex shadows. The memory growth is independent from the number of lights in the scene, but it is coherent with the increase in the visibility complexity.

The previous set of results shows that oversized area light sources can become a limitation for our approach, mainly because of the memory consumption. The present tests give a solution: Any huge light source, even with a critical size, could always be treated as the union of several smaller lights. And the result is always guaranteed noise free, which is specialy interesting with huge light sources since they require a high number of shadow rays using a stochastic method.

Focus on the BSP_query behaviour

Previous results demonstrate the global efficiency and robustness of this work, but they do not highlight the visibility algorithm (Algorithm 1) behaviour, i.e. how the visibility data and the computational cost evolve in respect to the visibility queries. Figure 11 gives such an insight. It focuses on the construction of a single BSP tree from the Conference model. This tree was selected because it is representative of the BSP_query behaviour. It has 2 896 occluders and a significant number of queries (17 878) are performed. In particular, most of the image points are located in the umbra or penumbra, which represent the most complex cases for our visibility algorithm.

At first, the tree grows quickly because there is no visibility data and the algorithm has to develop it to answer the visibility queries. The timings show the extra computational cost required for this construction. As a second step, the tree growth slows down drastically because the previously computed equivalence classes can be re-used and only need to be completed from time to time. As a consequence the computational cost falls down. This is the global behaviour of the visibility algorithm. In addition, notice that the image points are shaded in the scanline order. This allows handling consecutive points which are likely to share the same visibility data. This is noticeable locally: An "expensive" query is always followed by "cheaper" queries, taking advantage of the previous computational effort.

Discussions and future works

In this section we discuss our algorithm and point out some limitations and some issues we would like to address in the future. The abscissa corresponds to the number of queries. The left ordinate is the time in milliseconds and is related to the continuous curve, while the right one is the memory consumption in MB and is related to the dashed curve. Above the graphic, a half-tone illustration of the light visibility for each query. Black means invisible while white is fully visible.

The lazy construction of the visibility data is an important feature of our algorithm. This requires to query the visibility from some points to drive the computations. On one hand, this allows our approach to fit the image resolution and to focus the computational effort where and when it is needed. On the other hand, it may be inappropriate to solve some visibility problems. For example it is not suitable to proove that the light source is invisible from a polygon or to compute exact Potentially Visible Sets. Indeed, it is not possible to query all the points on a polygon or on the PVS boundaries. However, using an adaptative sampling of the surfaces for example, our approach can still be used to query the visibility from each samples and produce an accurate solution to these problems.

In this paper, a BSP tree is built for each triangle with image points to shade and makes the most of the visibility coherence between those points. This is an easy solution to group image points. However, if the geometric resolution is very high in respect to the image resolution, it can lead to very few image points per triangle and thus, a loss of efficiency, because the visibility data may be dropped before being re-used. To overcome this problem, image points should be grouped according to their neighborhood relationship instead of the triangle they belong to. This will avoid the problem by giving control over the group size and thus, over the number of queries per BSP tree.

As a future work, we think our implementation could be improved using a parallel computing architecture such as CUDA or OpenCL to benefit from a higher number of threads on the graphic hardware. This will require fixing some technical issues, in particular about memory management. Indeed, dynamic allocations from the graphic device are not natively supported, which is not convenient for the lazy construction of data structures.

In this work, the visibility data is dropped as soon as it is not needed anymore. Instead, we are thinking it could be saved to be used again for other renderings. In addition, we plan to study how the visibility data could be simplified and compressed. This may further improve the memory footprint and the visibility queries efficiency.

At last, we also plan to introduce depth information in the visibility data, which means grouping lines according to the first triangle they intersect. This would allow extending our approach to other problems and applications such as ambient occlusion.

Conclusion

In this paper, we have presented a new approach to compute exact from-polygon visibility. The core of this method is an original algorithm based on equivalence classes of lines. It encodes lazily the visibility from a polygon, and allows taking advantage of the visibility coherence between successive visibility queries. In addition, the algorithm implementation is easy, and uses only basic geometric operations. As a consequence, it is fast and robust in contrast to previous works on exact from-polygon visibility.

Using this core algorithm, we have proposed an alias-free solution to analytically compute high quality soft shadows. We have tested our approach in various configurations, including different models, different sizes and different numbers of area light sources. These results demonstrated the reliability of our method.

To conclude, we hope this work illustrates that exact visibility, especially exact from-polygon visibility, is not necessarily incompatible with robustness and performance. Moreover, since visibility is a recurrent question in computer graphics, soft shadows is not the only problem that could benefit from the core of our method. As discussed in the previous section, we are interested in several issues and we intend to continue this work.

Figure 1 :

 1 Figure1: Above: The relative orientation of two lines is given by the sign of the side operator. Under: Since side(l, r) = h l (π r ), applying the side operator comes down to testing the position of the Plücker point of one line against the dual hyperplane of the other line

Figure 3 :

 3 Figure 3: Left: 2 triangles and 3 lines a, b, c in various configurations. Right: The arrangement of hyperplanes (illustrated by 6 2D lines) mapped from the 6 triangles edges. They divide the Plücker space into cells. Filled cells are set of lines intersecting at least one triangle. π a , π b and π c are the Plücker points mapped from a, b and c. They are located in the cells they belong to, according to the triangle(s) they stab. For example, π b has a consistent orientation in respect to the 6 hyperplanes since b intersects the two triangles. Its relevant cell holds all the lines intersecting the two triangles.

Figure 4 :Figure 5 :

 45 Figure4: Let l 0 , l 1 , l 2 be the lines defined by the edges of an occluder o, and h 0 , h 1 , h 2 their dual hyperplanes in Plücker space. We can build a BSP tree whose leaves are the four equivalence classes generated by the triangle: Three visible classes (left leaves) and one invisible class (right leaf). We note bsp(o) such a representation.

Figure 6 :

 6 Figure 6: The point xyz and the line h define a plane that sets apart the lines having a positive or negative orientation in respect to h

1

  : build visible_triangle, the triangle list visible from the camera 2: foreach light L in the scene do 3: triangle_list ← visible_triangle 4: // the following loop parallelization is straightforward 5: while triangle_list is not empty do 6: remove a triangle T from the triangle_list 7: select the occluders O of T and L 8: initialize a BSP tree root node n associated with O 9: foreach image point xyz on T do 10: visible_parts ← BSP_query( n, xyz, L)

Figure 7 :

 7 Figure 7: From top to bottom: T-Rex, Sponza with Neptune, Conference and Soda Hall. The middle column shows the results computed with our algorithm. The left column presents the same pictures computed at comparable time and the right one at comparable quality, both using shadow rays.

Figure 8 :

 8 Figure8: The time and memory consumption using our algorithm. The Max Size column (MS) is the maximum memory load reached during the process. The Shaft Culling column (SC) gives the time percentage spent to select the occluders. The Init column (I) gives the time percentage spent to initialize each BSP tree. The Queries columns (Q) gives the time percentage spent to query the BSP trees for shading each image point. The Total column (T) gives the time in seconds for the whole process. MV Time column (MV) presents the result obtained using a Modified Version of our framework, where all queries are prevented from taking advantage of the others. The last column gives the falling off factor (FO) between our framework (Total) and its modified version (MV Time).

  presents the results. The left

Figure 9 :

 9 Figure 9: Increasing the light size. Left: the time consumption (top) and the maximum memory load (down) reached during the computations. Right: 3 pictures from our tests with the smallest, the medium and the largest area light source (100 times larger than the smallest one).

Figure 10 :

 10 Figure 10: Increasing the light number. Left: the time consumption (top) and the maximum memory (down) load reached during the computations. Left: 3 picture from our tests with 2, 16 and 36 area light sources.

Figure 11 :

 11 Figure11: Lazy construction of a representative BSP tree. The abscissa corresponds to the number of queries. The left ordinate is the time in milliseconds and is related to the continuous curve, while the right one is the memory consumption in MB and is related to the dashed curve. Above the graphic, a half-tone illustration of the light visibility for each query. Black means invisible while white is fully visible.

  Visibility algorithm: Given a point xyz on T , the algorithm answers the query "which parts of L are visible from xyz ? L is used to drive the BSP tree construction and may be split into several fragments representing homogeneous sets of lines from xyz. The fragments reaching visible classes are the visible parts of L from xyz

	3:	if n is not an invisible class then
	4:	n.occluders ← n.occluders ∪ o
	5:	end if
	6: else if orientation( o.vvset, n.hyperplane) > 0 then
	7:	insert( n.le f t, o )
	8: else if orientation( o.vvset, n.hyperplane ) < 0 then
	9:	insert( n.right, o )
	10: else
	11:	insert( n.le f t, o )
	12:	insert( n.right, o )
	13: end if
	14:
	15: BSP_query(Node n, Polygon L, Point xyz ) return Polygons
	16: loop
	17:	while n is not a leaf do
	18:	Plane p ← makePlane( xyz, n.hyperplane )
	19:	if position( p, L ) > 0 then
	20:	n ← n.le f t
	21:	else if position( p, L ) < 0 then
	22:	n ← n.right
	23:	else
	24:	// split the light and work recursively
	25:
	26:	end if
	27:	end while
	28:	if n.occluders is empty then
	29:
	30:	else
	31:	ro ← random occluder in n.occluders
	32:
	34:	insertOccluder( n, o )
	35:	end for
	36:	end if
	37: end loop

Algorithm 1 1: insertOccluder(Node n, Occluder o) 2: if n is a leaf then return BSP_query(n.le f t, L ∩ p + ) ∪ BSP_query(n.right, L ∩ p -) return (n.class is visible) ? L : ∅ n ← root of bsp(ro) 33: for o in n.occluders, o = ro do

Acknowledgments: Sponza Atrium by Marko Dabrovic. Neptune model is provided courtesy of Laurent Saboret (IMATI, INRIA) by the AIM@SHAPE Shape Repository.