Antoine Courtay

Mickaël Le Gentil

Olivier Berder

Arnaud Carer

Pascal Scalart

Olivier Sentieys

Zyggie: A Wireless Body Area Network platform for indoor positioning and motion tracking

Nowadays, there is a high demand for human and/or objects monitoring/localizing in the context of applications like Building Information Modeling (BIM), automated drone missions, contextual visits of museum or sports monitoring for instance. While for outdoor positioning accurate and robust solutions (i.e. GPS) exist for many years, indoor positioning is still very challenging. There is also a need of gesture/motion tracking systems that could replace video solutions. We propose in this paper a hardware/software platform named Zyggie that combines both Ultra Wide Band (UWB) technology and Received Signal Strength Indicator (RSSI) for low power accurate indoor positioning and Inertial Measurement Unit (IMU) utilization for motion tracking. Very few industrial/academic existing solutions can simultaneously perform indoor positioning and motion tracking and none of them can do both under low power, low cost and compacity constraints addressed by our platform. As Zyggie has the capability to estimate distances w.r.t other platforms in the environment and quaternions (which represents the attitude/orientation) users can test/enhance state of the art algorithms for positioning and motion tracking applications.

I. INTRODUCTION

Knowing the posture/location of a person in its environment becomes an important issue for applications providing services based on these criteria. Many solutions have been proposed going from satellite based techniques (i.e. GPS), video tracking or radio-electrical waves combining. Some of them use optical signals such as infrared, laser or video tracking approaches, some use radio signal information (Time Of Flight (TOF) or RSSI) and others, data provided by IMU. To choose the adequate component for the platform design and obtain the best trade-off for considered applications, a comparison in terms of accuracy, deployment facility, size, power consumption and price has to be done. As both indoor positioning and motion tracking need low cost, portable and low power systems which can be deployed in minutes, radio and IMU are the best candidates. The main contribution of this paper is a platform named Zyggie that combines hardware components such as UWB chip for TOF distances measurement, 2.4 GHz transceiver for RSSI sensing and an IMU for motion estimations. Several Zyggie platforms are used to form a Wireless Sensor Network (WSN) which medium access has to be scheduled. To do so a custom Time Division Multiple Access Medium Access Control (TDMA MAC) has been developed to get the best trade-off between consumption and performance. Another contribution is the data processing chain which is composed of some filtering and calibration phases to format data for positioning and motion tracking algorithms. Two GUIs have been developed for the this purpose. First a Matlab one for testing/validating algorithms, then an Android tablet one for displaying on a 3D avatar the results in real-time. The rest of this paper is organized as follows. Section II summarizes state of the art platforms that propose either positioning and/or motion tracking capabilities. Section III introduces the hardware of our new developed platform named Zyggie. Data exchanges provided by Zyggie are then explained in Section IV using a custom TDMA MAC protocol. In Section V output data from Zyggie are used to evaluate positioning and motion tracking applications (section V-A & V-B). Finally, conclusions are drawn and future works are commented in Section VI. • Price of the system. Vicon system [START_REF]Vicon system[END_REF] is known as a reference in motion tracking applications due to its millimeter accuracy. Several infrared cameras are used to track retroreflective markers. This opticalpassive solution allows a long period of use but leads to significant computation resources and power needs. Kinect [START_REF] Zhang | Microsoft kinect sensor and its effect[END_REF] is a low cost general public markerless solution that embeds infrared and RGB cameras. Depth matrix coupled with image analysis enable the system to track in 3D space the attitude of a person. As Kinect has a unique camera viewpoint, it can lead in some cases to an information loss (i.e. arm hidden by the user). Unlike video tracking solutions Xsens is only based on IMU. High quality hardware and complex processing algorithms needed for motion tracking applications lead to a good accuracy but high power consumption. ShotTracker [3] is a basket system monitoring which allows the positioning of players and the ball in addition to scoring capability. This solution, based on UWB technology, is only dedicated to the previously mentioned application due to the fact that the system is fixed to the building. Bluetooth Low Energy (BLE), ZigBee or other Radio Frequency (RF) technologies [START_REF] Dortz | Wifi fingerprint indoor positioning system using probability distribution comparison[END_REF]- [START_REF] Aulery | Radio signature based posture recognition using wbsn[END_REF] can use received signal power to get distance information between network's nodes. Obtaining this RSSI is not consuming but is inaccurate (more or less 2.5 m) due to its environment perturbations' sensitivity. This makes it not suitable for some of the applications. The global view of Table I highlights the lack of a mobile, low cost, low power system with reasonable accuracy for both targeted applications. This lack is filled with Zyggie platform, detailed in Section III.

II. RELATED WORK

III. THE Zyggie PLATFORM

A. UWB Technology

UWB systems allow to emit low power and ultra-short duration pulses and provide applications with following advantages:

• High precision TOF ranging at centimeter level,

• Low electromagnetic pollution and reduced interference to other systems, • Small size and low power consumption. TOF computation is explained as following: a message is sent from the emitter to the receiver and after some processing time the receiver replies back. The TOF from the emitter to the receiver is the difference between the time that the emitter received the response from the receiver and the emission time divided by 2, by taking care of the response time of the receiver. Thus, we can compute the distance d depending on TOF with the following formula

d = T OF.c = T re -T em -T r 2 .c , (1)
where T em is the emission time, T re the response reception time, T r the response delay of the receiver, and c the celerity which is assumed to be the same as the light in air. Due to some drifts on the emitter and receiver frequencies/clocks, an error is present on the measurements. In order to minimize this error this scheme can be repeated twice.

There are only few vendors of UWB chipsets among which we can cite BeSpoon, DecaWave [START_REF]Decawave uwb chipset[END_REF], Blinksight and Alereon.

As DecaWave chips are the most open for experimentation (frequency band, datarate, transmitting power programming), they were chosen for our prototype (see Section III-B). The Zyggie node is a Wireless Sensor Node with UWB radio capabilities that was developed in our lab. Its architecture is represented on Fig. 1

B. Zyggie Node architecture

and composed of:

• A Silicon Labs EFM32 ARM CortexM4 micro-controller which is the Master Control Unit (MCU) of the node and also controls the UWB radio chip, • An ATMEGA2564RFR2 which drives the 2.4 GHz radio for data exchanges and network synchronization, • A DecaWave DW1000 UWB radio chip [START_REF]Decawave uwb chipset[END_REF] for TOF measurements, • An InvenSense IMU for acceleration, angular speed and magnetic field captures. Our node is ultra compact (37 mm by 63 mm) and only weights 32 g (including battery) which makes it a real good candidate for all applications where weight and size are critical. Some experiments with UWB radio enabled at full speed were done leading to a battery life time of around 5 hours which is enough for sportsman or in building applications.

IV. Zyggie MAC PROTOCOL

The network is composed of identical hardware Zyggie nodes. One of them is the central node (coordinator) connected to a computer/tablet, the others are on-body nodes and, if necessary, off-body fixed nodes for positioning. Coordinator aggregates data from all Zyggie nodes (Z i with i ∈ [1, N ] and N the number of nodes) to process and plot them on GUI Fig. 2. Custom TDMA MAC presented in Section V-B. A custom TDMA MAC protocol allocates time slots to each node to allow it sending its sensed data. All nodes have to sense movement data like IMU information or system absolute orientation, RF messages energy called RSSI and distance estimation with TOF measurements. Our MAC is based on the IEEE 802.15.4 communication protocol with a beacon-enable network. The coordinator broadcasts periodically a beacon frame to synchronize nodes and share network information. As shown on Fig. 2, a full mesh communication is started to get RSSI information. In its RSSI time slot, each node sends a short broadcast frame (orange TX slots) that others will listen (red RX slots) to get a complete RSSI matrix. As the IMU measurement step starts (purple time slot), the MCU gets data (red IMU Processing slots) and optionally performs an Extended Kalman Filtering (EKF) to produce a quaternion. Each node in its time slot (Data frame time slot) sends these data via the 2.4 GHz radio transceiver to the coordinator. During this star communication phase, nodes also process TOF ranging measurements. First, the UWB chip has to be awaken (green Awaking slot), then it can wait to be used (yellow Idle slots) and finally it does its ranging (orange & Red TX/RX slots) before returning to sleep. In the example of Fig. 2 the node 2 is selected to do ranging toward others, therefore a distances' vector can be sent to the coordinator at its next Data frame. A sequence of selected ranging nodes has to be programmed to get a full distances matrix after N beacons. Previously defined phases are scheduled to get the best trade-off between consumption and performance, that is the reason why transceivers spend a minimum time into receiving mode. On the other hand, an adaptive data-rate is implemented depending on the number of nodes, to deal with the processing and communication time.

V. PROOF OF CONCEPT

As previously said, our platform can output data that can be used by users to test algorithms in the context of positioning and motion tracking. To be user friendly, we have also developed a GUI modeling an avatar in 3D space to visualize the results. These two aspects are now going to be explained in more details.

A. Algorithms performance evaluation

Algorithms take part of the processing chain as shown in Fig. 3. On the chain at the top of the figure, UWB distances allow to deal with positioning application. On the chain at the bottom of the figure, computed information coming from IMU (i.e. quaternions) and UWB distances allow to deal with motion tracking application. For both applications, it exists state of the art algorithms (coloured boxes at the end of the chains) but users can also develop and/or optimize these algorithms. Blocks involved in these chains are now going to be presented. As we are able to obtain a measured distance between two nodes the next step is to compute several distances between a mobile node and some fixed nodes named anchors. The process chain starts with some corrections depending on antenna delays and hardware design. Then as TOF measurements can be noisy, some filtering is done to remove aberrant values. Finally, to obtain the coordinates of the mobile a localization algorithm is used (blue Positioning box in the chain) among we can cite centroid, trilateration, multilateration, Gauss-Newton [START_REF] Blumenthal | Weighted centroid localization in zigbee-based sensor networks[END_REF]- [START_REF] Nocedal | Numerical Optimization: 10.3 Algorithms for Nonlinear Least-Squares Problems[END_REF] and so on. Users can also test their algorithms with data provided by Zyggie (i.e. UWB distances and/or RSSI). For instance, we experimented 2D positioning of a walking man using Zyggie with a custom low computation noise-aware algorithm based on trilateration. For this experimentation, data were acquired in a room of nearly 243 m 2 (26.4 m length by 9.2 m width). The ideal trajectory was marked on the floor and one step by second (around 3 km/h) was done following this trajectory. In order to test our algorithm a six anchors configuration is used. Fig. 4 represents the computed trajectory by our algorithm in red and the ideal one in green. 2) Motion tracking: First, we implemented the proprietary InvenSense EKF on the 32bits MCU in the Zyggie platform to send less information via RF. This data fusion algorithm converts raw IMU data (i.e. acceleration, rotation and magnitude) into quaternions to represent absolute system orientation on 3 axes. Then, user needs to do 2 specific postures for Initial and North calibrations (the first two boxes in the bottom chain of Fig. 3). Classical "T" and "I" postures were chosen to measure arm movements and calibrate user's posture. When arms are raised to form a "T", all quaternions are saved to achieve the computation of

q initjoint = q rawjoint .q -1 0joint , (2)
where joint represents the joint name associated to a node, q initjoint the initialized output orientation, q rawjoint the measured orientation and q 0joint the initial quaternion equals to q -1 rawjoint only during the first calibration posture. After the second posture done, we select one arm node (i.e. left forearm) to resolve the System (3) and find the q north quaternion.

q T = q -1 north .q T .q north q I = q -1 north .q init LF .q north , (3)
where q T is the quaternion to form the "T" pose, q I the wanted quaternion to form the "I" pose and q init LF the measured and initialized quaternion for the "I" pose. This system is in the form of the Wahba's problem therefore a solution can be found with Singular Value Decomposition (SVD) as described in [START_REF] Bates | Gcvpack routines for generalized cross validation[END_REF]. Solutions can be found with the other arm nodes to get a single quaternion, therefore reducing the measuring noise due to both the sensor and user approximate posture. Finally, depending on the targeted application, user's movement or posture can be processed easily with data (i.e. quater-nions and/or RSSI and/or UWB distances) provided by Zyggie (orange Motion Tracking box in the chain).

B. Graphical User Interface

These algorithms were implemented into two GUIs for different needs. First of all, Matlab was used to make acquisitions and work on raw data. It allows us to display in real-time raw data like RSSI and UWB distances table and accelerometer, gyroscope and magnetometer graphics. Other features are available for posture and movement recognition with Principal Component Analysis (PCA) classification algorithm [START_REF] Aulery | Low-complexity energy proportional posture/gesture recognition based on wbsn[END_REF]. Targeted applications are in medicine for functional rehabilitation of a patient or in Human-Machine interaction with gesture. Then, an Android tablet application was developed for demonstrations to show the system ability in a mobile context. Features like IEEE 802.15.4 channel scanning at 2.4 GHz and 3D avatar tracking were added. Positioning algorithms were also included to update avatar position in a 3D environment as displayed on Fig. 5.

VI. CONCLUSION

We have proposed in this paper a new low power, low cost hardware/software platform named Zyggie for indoor positioning and motion tracking applications. Zyggie combines UWB technology for accurate indoor positioning and IMU for motion tracking; RSSI metrics can also be inputted for both applications. Our prototype has been validated in real use cases (indoor human tracking). An Android application has also been developed for real-time visualization on a tablet. Finally Zyggie offers the opportunity for users to test their own developed algorithms in the Matlab environment.

Fig. 1 .

 1 Fig. 1. Zyggie WSN Hardware illustration.

Fig. 3 .

 3 Fig. 3. Zyggie data processing chain

Fig. 4 .

 4 Fig. 4. Results for indoor walking context

Fig. 5 .

 5 Fig. 5. Android application: blue dots are the anchors for positioning and the tennis player is also equipped with nodes for motion tracking.

Table I

 I lists most relevant solutions for either indoor positioning and/or motion tracking applications. These solutions are compared in terms of:

	• Deployment which represents the system setup facility
	and the flexibility toward the application,
	• Size of wearable devices (if needed) and/or the associated
	system infrastructure,

• Indoor Positioning and/or Motion Tracking capabilities, • Accuracy of measurements, • Power consumption,

TABLE I INDOOR

 I POSITIONING AND MOTION TRACKING SOLUTIONS COMPARISON

	Technology	Indoor Positioning	Motion Tracking	Accuracy	Deployment	Size	Power	Price
	Video w/ marker [1]	Yes	Yes	+++	---	---	---	---
	Video w/o marker [2]	Yes	Yes	-	+	+++	-	+++
	IMU [10]	No	Yes	++	-	+	--	-
	UWB [3]	Yes	No	++	--	+	+	-
	RSSI [4]-[8]	Yes	No	---	-	++	+++	++
	UWB / IMU / RSSI	Yes	Yes	+	++	++	++	++

ACKNOWLEDGMENTS

This work has been conducted in the context of the Cominlabs Labex BoWI project [17] funded by the "Investissement d'avenir" program.


Maria Letizia Raffa

Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Raffaella Rizzoni
email: raffaella.rizzoni@unife.it

Maria Letizia Raffa
email: maria-letizia.raffa@u-pec.fr

Derivation of a model of imperfect interface with finite strains and damage by asymptotic techniques. Application to masonry structures

Keywords: Bonding, Asymptotic analysis, Finite strains, Damage, Imperfect interface, brick/mortar interface

come

Introduction

In the last decades, the safeguard of ancient masonry monumental structures has animated many researchers to develop new constitutive models and innovative tools of structural analysis able to simulate the complex response of such a structure subject to static and dynamic loads.

The complexity of mechanical behavior of masonry is due to its composite nature. Masonry is an assembly of blocks of artificial or natural origin joined by dry or mortar joints. The presence of the mortar joints influences the mechanical response of the structure and provides elements of weakness. Moreover, it is well-established that brick/mortar interfaces constitute privileged zones for microcracks nucleation and propagation.

To effectively model the masonry material, two main approaches are usually adopted: a macromodelling approach which considers an equivalent fictitious material where the presence of the mortar is smeared out (e.g., [START_REF] Anthoine | Derivation of the in-plane elastic characteristics of masonry through homogenization theory[END_REF][START_REF] Lourenço | Continuum model for masonry: Parameter estimation and validation[END_REF][START_REF] Milani | Homogenised limit analysis of masonry walls, part I: failure surfaces[END_REF][START_REF] Medeiros | Numerical modelling of non-confined and confined masonry walls[END_REF][START_REF] Rots | Numerical simulation of cracking in structural masonry[END_REF][START_REF] Sacco | A nonlinear homogenization procedure for periodic masonry[END_REF]) and a micromodelling approach which analyzes the structure as an assembly of blocks connected by interfaces (e.g., [13,[START_REF] Lourenço | Computational strategies for masonry structures[END_REF]26,[START_REF] Lotfi | Interface model applied to fracture of masonry structures[END_REF][START_REF] Zucchini | A micro-mechanical model for the homogenisation of masonry[END_REF]). The latter approach turns out to provide a more realistic description of the masonry, because it allows to adopt suitable constitutive models for each components constituting the masonry material more advantageous to the study case [START_REF] Fileccia Scimeni | The interphase model applied to the analysis of masonry structures[END_REF][START_REF] Giambanco | The interphase finite element[END_REF].

Because of the increasing interest in describing the constitutive behavior of composite structures (i.e., fiber-reinforced materials, laminate composites, masonry, etc.), a number of interface models have been developed in the literature which take into account different interface behaviors, like elasticity, plasticity, viscoelasticity, thermal effects, residual stresses and damaging (e.g., [3-6, 8, 10, 14, 15, 17-24, 33, 35, 39-42, 45]).

Constitutive models for interphases/interfaces in masonry structures are usually formulated in terms of the contact traction and relative displacements of the two surfaces interacting through the mortar joints [START_REF] Alfano | Combining interface damage and friction in a cohesive-zone model[END_REF][START_REF] Fileccia Scimeni | The interphase model applied to the analysis of masonry structures[END_REF][START_REF] Giambanco | The interphase finite element[END_REF]. Thus, the mortar joints are generally approximated with material interphases in which material nonlinearities eventually occur. From a numerical point of view, these constitutive laws being an internal spring-like boundary condition can be implemented within the finite element formulation via zero-thickness interface elements [START_REF] Giambanco | Numerical analysis of masonry structures via interface models[END_REF].

In this paper, a regular masonry made of blocks connected together by regular but not too thin layers of mortar is considered. The original contribution of the paper consists in the introduction of a third thin interphase interposed between bricks and mortar joints. This material interphase represents a thin layer which has the physico-chemical properties allowing the bonding between bricks and mortar joints.

A novelty with respect the existing literature about micromodelling approach for masonry structures is thus represented by the introduction of a third component which has different mechanical properties from brick and mortar and which undergoes geometrical nonlinearities (i.e., large displacements) together with material nonlinearities (i.e., damaging). To account for large displacements, the Saint Venant-Kirchhoff constitutive model is adopted, which is the simplest among nonlinear elastic material models. Damaging behavior is modeled using a procedure similar to the approach introduced by Bonetti et al. [START_REF] Bonetti | A model of imperfect interface with damage[END_REF]: the thin interphase interposed between bricks and mortar is assumed to be a microcracked material undergoing a degradation process ruled by a pseudopotential of dissipation. The microcracked material is a generalized Kachanov-type material. In particular, within the Kachanov theory [START_REF] Sevostianov | On some controversial issues in effective field approaches to the problem of the overall elastic properties[END_REF][START_REF] Tsukrov | Effective moduli of an anisotropic material with elliptical holes of arbitrary orientational distribution[END_REF], constitutive equations are obtained via the homogenization of a microcracked material characterized by k families of cracks. The elastic coefficients depend on the lengths l k of these cracks. In the present study, only one family of microcracks is considered and l is taken to denote their equivalent length. Moreover, the interaction between microcracks is neglected (non-interacting approximation, cf. [START_REF] Sevostianov | On some controversial issues in effective field approaches to the problem of the overall elastic properties[END_REF]). As a result, the equivalent microcracks length l represents our damage parameter. Thus, according to the Kachanov theory for microcracked media, the elasticity coefficients of the Saint Venant-Kirchhoff material are derived and they explicitly depend on the evolving damage parameter l.

This novel approach used to derive the imperfect interface laws coupling homogenization for microcracked media and matched asymptotic techniques is called Imperfect Interface Approach (IIA) [START_REF] Raffa | Micromechanical modeling of imperfect interfaces and applications[END_REF][START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part I: theory[END_REF][START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part II: in silico analysis[END_REF]. Within the large displacement framework, the IIA allows to derive a nonlinear spring-like model which can be easily implemented in a finite element analysis. In analogy to former studies [START_REF] Giambanco | Numerical analysis of masonry structures via interface models[END_REF], the proposed interface model is implemented in a commercial FEM-based software via zero-thickness interface elements including damage evolution.

The present paper is organized as follows. Section 2 is devoted to the formulation of the constitutive aspects of the damaged interphase undergoing large displacements. Then, the problem of a composite body constituted of three deformable solids, a block, a thin microcracked damaging interphase and a mortar joint bonded together, is investigated in Section 3. In details, an asymptotic expansion method is applied to the elastic problem of the composite and a model of nonlinear imperfect interface for the thin interphase is derived. In such a model, surface damage effects turn out to be included in the limit as the thickness of the thin interphase vanishes. The numerical validation is detailed in Section 4. Firstly, an analytical analysis of the brick/mortar interface is carried out, in which I and II pure loading modes in the sense of fracture mechanics are investigated. Then, two finite element analyses based on the extensive experimental campaign carried out in the nineties by Vermeltfoort and Raijmakers [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] concerning shear tests on masonry walls, are proposed. Numerical results as well as comparisons with the experimental data are discussed. Finally, conclusions are presented in Section 5.

A model of interphase

Generalities

In this study a regular masonry is considered, i.e., the masonry structure is made of parallelepipedic bricks linked together by regular but not too thin layers of mortar. The novelty of this approach lies in considering as a third deformable body, a thin interphase between the brick and the mortar. The physico-chemical properties of this interphase allow the bonding between the masonry principal constituents. In the following, and without loss of generality, the interphase thickness will be considered as a constant. The interphase mechanical model is detailed in the next paragraph.

A model of damaged material

In this section, a deformable solid occupying a smooth bounded domain B of IR 3 is considered. An orthonormal Cartesian frame (O, e 1 , e 2 , e 3 ) is introduced and let (x 1 , x 2 , x 3 ) be taken to denote the three coordinates of a particle. The solid is assumed to comprise a material:

hyperelastic of a Saint Venant-Kirchhoff-type [START_REF] Rizzoni | On Saint Venant -Kirchhoff Imperfect Interfaces[END_REF][START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part I: theory[END_REF] microcracked of a generalized Kachanov-type [START_REF] Mauge | Effective elastic properties of an anisotropic material with arbitrarily oriented interacting cracks[END_REF][START_REF] Tsukrov | Effective moduli of an anisotropic material with elliptical holes of arbitrary orientational distribution[END_REF] As a result, the elastic coefficients depend on the length l, which is taken to denote the equivalent cracks length, therefore, it has the meaning of a damage parameter. As an example, for an isotropic material weakened by a cracks family lying orthogonal to e 3 , the Young modulus in the third direction E 3 is equal to

E 3 = E 0 1 + 2 ρ C E 0 , (1)
with

C = π 2 1 √ E 0 1 µ 0 -2 ν 0 E 0 + 2 E 0 1/2 (2)
where E 0 (resp. µ 0 , resp. ν 0 ) is the Young modulus (resp. the shear modulus, resp. the Poisson ratio) of the undamaged material, and ρ is the cracks density, which form is ρ = l 3 V in 3D problems and ρ = l 2 S in 2D, being V (resp. S) the volume (resp. the surface) of the representative elementary domain [START_REF] Mauge | Effective elastic properties of an anisotropic material with arbitrarily oriented interacting cracks[END_REF][START_REF] Tsukrov | Effective moduli of an anisotropic material with elliptical holes of arbitrary orientational distribution[END_REF].

Let the evolution of crack length l be introduced. A pseudo-potential of dissipation φ is considered [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF]. As proposed in former papers [START_REF] Bonetti | A model of imperfect interface with damage[END_REF][START_REF] Orefice | An experimental/numerical study on the interfacial damage of bonded joints for Fibre-Reinforced Polymer profiles at service conditions[END_REF], the chosen pseudo-potential is given by the sum of a quadratic term and a positively 1homogeneous functional, as a result the dissipative character will be directly related to a rate-dependent and a rate-independent contributions:

φ( l) = 1 2 η l2 + I [0,+∞[ ( l), (3)
where η is a positive viscosity parameter and I A denotes the indicator function of the set A, i.e., I A (x) = 0 if x ∈ A and I A (x) = +∞ otherwise. The term I [0,+∞[ ( l) forces l to assume non-negative values (i.e., the crack length cannot decrease) and it renders the irreversible character of the degradation process of the bond. Moreover, the free energy of the solid is chosen as follows:

ψ(E, l) = 1 2 B(l)E(u) : E(u) -ωl + I [l0,∞[ (l) (4)
where E(u) is the Green-Lagrange tensor, E(u) = 1 2 ∇u + ∇u T + ∇u T ∇u , B(l) is the stiffness tensor of the body, with usual properties of symmetry and positivity, ω is a (negative) parameter similar to the Dupré's energy [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF][START_REF] Fouchal | Contribution to the modelling of interfaces in masonry construction[END_REF] and l 0 is a given initial crack length. It is worth to highlight that, due to the weak thickness of the bond, it seems natural to introduce finite strains in the interphase. The constitutive equation is given by:

S = F -1 P = B(l)E(u) (5)
where P (resp. S) is the first (resp. second) Piola stress tensor and F = I +∇u is the gradient of transformation tensor.

Remark. It is possible to limit the crack length by assuming a vanishing stiffness if l ≤ l max ; this can be done by introducing the characteristic function of [l 0 , l max ], χ [l0,lmax] , and by assuming for the free energy the following form:

ψ(E, l) = χ [l0,lmax] ( 1 2 B(l)E(u) : E(u) -ωl) + I [l0,∞[ (l).
Therefore, the equilibrium of the body is written as:

DivP + f = 0, in B (6)
where f is a given body force. Finally, from the chosen form of the pseudo-potential of dissipation φ (Eq. ( 3)) and because of the positivity of the parameter η > 0, the following equation governing the damage parameter l results:

η l = ω - 1 2 B ,l (l)E(u ε ) : E(u) + , (7)
where (•) + is taken to denote the positive part of a function.

A model of rescaled damaged interphase

The thin interphase domain B ε is assumed to have a constant thickness ε in the third direction e 3 (see Fig. 2), i.e.,

B ε = x ∈ IR 3 , |x 3 | < ε/2 . ( 8
)
In the following, all the parameters and fields are indexed by ε. The upper and lower surfaces of the thin interphase are denoted as

S ε ± = x ∈ IR 3 , x 3 = ±ε/2 and S ε L = ∂B ε /S ε ± , (9)
where ∂B ε indicates the boundary of B ε . Body forces in B ε and lateral forces on S ε L are neglected. Moreover, the material body is assumed to be soft, i.e., B ijhk (l) ≈ εb ijhk (l).

Thus, the equilibrium problem of the body is described by the following system of equations:

               P ε ij,j = 0 in B ε P ε ij n j = 0 on S ε L P ε ij n j = P ± i on S ε ± P ε ij = εF ε ir b rjhk (l)E hk (u ε ) in B ε η ε l = ω ε - 1 2 εb ,l (l)E(u ε ) : E(u ε ) + in B ε (10)
where P + = P -are a priori given vectors.

Since the weak thickness of the interphase, it is natural to seek the solution of problem (61) using asymptotic expansions with respect to the small parameter ε [START_REF] Lebon | Asymptotic study on a soft thin layer: The non-convex case[END_REF][START_REF] Lebon | Asymptotic analysis of a thin interface: The case involving similar rigidity[END_REF][START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: An energy approach[END_REF][START_REF] Lebon | Asymptotic analysis of Drucker-Prager and Mohr-Coulomb soft thin interfaces[END_REF][START_REF] Lebon | Asymptotic modelling of interfaces taking contact conditions into account: Asymptotic expansions and numerical implementation[END_REF][START_REF] Rizzoni | Higher order model for soft and hard interfaces[END_REF]. In particular, the following asymptotic series with fractional powers are assumed [START_REF] Rizzoni | On Saint Venant -Kirchhoff Imperfect Interfaces[END_REF][START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part I: theory[END_REF]:

u ε = u 0 + ε 1/3 u 1 + ε 2/3 u 2 + εu 3 + o(ε) P ε = P 0 + ε 1/3 P 1 + ε 2/3 P 2 + εP 3 + o(ε) (11)
The following change of variable is introduced [START_REF] Ciarlet | Mathematical Elasticity. Volume I: Three-Dimensional Elasticity[END_REF]:

(x 1 , x 2 , x 3 ) ∈ B ε → (z 1 , z 2 , z 3 ) ∈ B, with (z 1 , z 2 , z 3 ) = (x 1 , x 2 , x 3 ε ) and it is set ûε (z 1 , z 2 , z 3 ) = u ε (x 1 , x 2 , x 3 ) and P ε (z 1 , z 2 , z 3 ) = P ε (x 1 , x 2 , x 3 ), where B = {(x 1 , x 2 , x 3 ) ∈ Ω : |x 3 | < 1 2 }.
The superscript (∧) denotes the fields in the rescaled configuration.

The governing equations of the rescaled problem result in:

                 P ε ij,j = 0 in B P ε ij n j = 0 on S L P ε ij n j = P ± i on S ± P ε ij = ε F ε ir brjhk (l) Êhk (û ε ) in B ηε l = ωε - 1 2 ε b,l (l) Ê(û ε ) : Ê(û ε ) + in B (12)
where

S ± = {(x 1 , x 2 , x 3 ) ∈ Ω : x 3 = ± 1 2
} and S L = ∂B/S ± , ∂B are the boundaries of B. In Eqs. [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF], Ê(û ε ) is the rescaled Green-Lagrange tensor, and its asymptotic expansion reads as:

Ê(û ε ) = ε -2 Ê-2 + ε -1 Ê-1 + Ê0 + O(ε), (13)
where

     Ê-2 = 1 2 | ûε ,3 | 2 e 3 ⊗ e 3 , Ê-1 = 1 2 u ε ,3 ⊗ e 3 + 1 2 (u ε ,3 • u ε ,α )(e α ⊗ e 3 ), Ê0 = 1 2 (u ε ,α ⊗ e α ) + 1 2 u ε 3,α (e α ⊗ e 3 + e 3 ⊗ e α ) + 1 2 (u ε ,α • u ε ,β )(e α ⊗ e β ). (14)
In view of Eqs. [START_REF] Fouchal | Contribution to the modelling of interfaces in masonry construction[END_REF] the displacement and stress fields are written as asymptotic expansions in the rescaled adhesive, as follows:

P ε = P 0 + ε 1/3 P 1 + ε 2/3 P 2 + ε P 1 + o(ε) ûε = û0 + ε 1/3 û1 + ε 2/3 û2 + εû 1 + o(ε). (15)
Substituting Eqs. [START_REF] Giambanco | Numerical analysis of masonry structures via interface models[END_REF] into the first equation of system [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF], the following conditions are found to hold in B:

P k i3,3 = 0, k = 0, 1, 2, i = 1, 2, 3 (16)
The latter result implies that P k i3 is independent of z 3 and, additionally, that the relative jumps are equal to zero:

P k i3 = 0, k = 0, 1, 2, i = 1, 2, 3 (17)
where the notation

[f ] = f (x 1 , x 2 , 1 2 ) -f (x 1 , x 2 , - 1
2 ) indicates the jumps of function f in the rescaled model.

Substituting Eqs. [START_REF] Giambanco | Numerical analysis of masonry structures via interface models[END_REF] and Eqs. ( 14) into the fourth equation of system [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF] it is deduced that the following conditions hold in B:

ûk i,3 = 0, k = 0, 1, i = 1, 2, 3 (18)
and û2 i,33 = 0, i = 1, 2, 3 (19)
i.e. ûk i and û2 i,3 are independent of z 3 , leading to the conditions:

ûk = 0, k = 0, 1, (20)
û2 ,3 = û2 . (21)
Moreover, the following relationship is obtained:

1 2 b3333 (l) û2 i,3 2 û2 i,3 ,3 = 0 ( 22
)
and thus 1 2 b3333 (l) û2 i,3

2 û2 i,3 = P 0 ( 23
)
where P 0 is the first term in the expansion of P . By introducing Eq. ( 21) in Eq. ( 23) the following equation is obtained:

P 0 = 1 2 b3333 (l) û2 2 û2 (24)
Let consider the fifth equation in system [START_REF] Frémond | Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications[END_REF]. The term

1 2 ε 1 2 - 1 2 b,l (l) Ê(û ε ) : Ê(û ε )dz 3 ( 25
)
has to be estimated. Using the asymptotic expansions of ûε and Eqs. ( 18)-( 19), the first non-vanishing term in the expansion (the second order term) is found to be governed by the expansion of Ê33 (û ε ). Eqs. ( 18) and ( 19) leads to

Ê33 (û ε ) ≈ 1 2 ε -2/3 (û 2 1,3 + û2 2,3 + û2 3,3 ) (26)
As a results, Eq. ( 25) can be rewritten as follows:

1 2 ε 1 2 - 1 2 b,l (l) Ê(û ε ) : Ê(û ε )dz 3 ≈ ≈ 1 2 ε b3333,l (l) 1 2 - 1 2
1 2 ε -2/3 (û 2 1,3 ) 2 + (û 2 2,3 ) 2 + (û 2 3,3 ) 2 2 dz 3 (27)
It is worth recalling that the second order terms in the displacements expansion û2 i,3 , i = 1, 2, 3 are found to be constant with respect to z 3 (cf. Eq. ( 19)), accordingly, the expression (û 2 1,3 ) 2 + (û 2 2,3 ) 2 + (û 2 3,3 ) 2 leads to û2 2 after integration.

Finally, the damage evolution equation results in:

η-1 l = ω-1 - 1 8 b3333,l (l) û2 4 + ( 28
)
where η -1 and ω -1 are the terms in the expansion of η ε and ω ε respectively of power -1/3,

ηε = ... + ε -1/3 η-1 + ... ωε = ... + ε -1/3 ω-1 + ..., (29)
One remarks that η ε and ω ε are assumed to be independent of the z 3 coordinate.

Study of a composite

The body B ε studied in the two previous sections is, in what follows, sandwiched between two hyperelastic bodies occupying respectively the domains Ω ε ± ⊂ IR 3 , denoted as adherents. It is assumed that Ω ε ± and B ε are perfectly bonded along the surfaces S ε ± . We assume that an external load g is applied on a boundary part S ± g ∈ ∂Ω ε ± /S ε ± , and that the boundary condition u = 0 holds on S ± u ∈ ∂Ω ε ± /S ε ± , with S ± u having strictly positive measure and such that S ± g ∩ S ± u = ∅. Finally, a body force f ± is applied to Ω ε ± .

Fig. 3 An example of the geometry of the composite

The equilibrium problem of the composite structure is described by the following system of equations:

                               P ε ij,j + f ± i = 0 in Ω ε ± P ε ij n j = g ± i on S g P ε ij,j = 0 in B ε [[P ε i3 ]] = 0 on S ε ± [[u ε i ]] = 0 on S ε ± u ε i = 0 on S ± u S ε ij = a ± ijhk E hk (u ε ) in Ω ε ± S ε ij = εb ijhk (l)E hk (u ε ) in B ε η ε l = ω ε - 1 2 εb ε ,l (l)E(u ε ) : E(u ε ) + in B ε (30)
where a ± are the elasticity tensors in Ω ε ± . The following change of variable is introduced in the adherents:

(x 1 , x 2 , x 3 ) ∈ Ω ε ± → (z 1 , z 2 , z 3 ) ∈ Ω ± , with (z 1 , z 2 , z 3 ) = (x 1 , x 2 , x 3 ± 1/2 ∓ ε/2), and it is set ūε (z 1 , z 2 , z 3 ) = u ε (x 1 , x 2 , x 3 ) and P ε (z 1 , z 2 , z 3 ) = P ε (x 1 , x 2 , x 3 ), where Ω ± = {(x 1 , x 2 , x 3 ) ∈ Ω : ±x 3 > 1 2 }.
The external forces are assumed to be independent of ε. As a consequence, it is set f ± (z 1 , z 2 , z 3 ) = f ± (x 1 , x 2 , x 3 ) and ḡ± (z 1 , z 2 , z 3 ) = g ± (x 1 , x 2 , x 3 ). Note that the superscript -(respectively ∧) denotes the fields in the adherents Ω ε ± (respectively in the interphase B ε ) in the rescaled configuration.

The governing equations of the rescaled problem are the following:

                                 P ε ij,j + f ± i = 0 in Ω ± P ε ij n j = ḡ± i on S± g P ε ij,j = 0 in B P ε i3 = P ε i3 on S ± ūε i = ûε i on S ± ūε i = 0 on Su Sε ij = ā± ijhk Ēhk (ū ε ) in Ω ± Ŝε ij = ε bijhk (l) Êhk (û ε ) in B ηε l = ωε - 1 2 ε b,l (l) Ê(û ε ) : Ê(û ε ) + in B (31)
In view of Eqs. [START_REF] Fouchal | Contribution to the modelling of interfaces in masonry construction[END_REF], the displacement and stress fields, written as asymptotic expansions, read as:

P ε = P 0 + ε 1/3 P 1 + ε 2/3 P 2 + ε P 3 + o(ε) ūε = ū0 + ε 1/3 ū1 + ε 2/3 ū2 + εū 3 + o(ε), (32)
Substituting Eqs. [START_REF] Orefice | An experimental/numerical study on the interfacial damage of bonded joints for Fibre-Reinforced Polymer profiles at service conditions[END_REF] in system [START_REF] Milani | Homogenised limit analysis of masonry walls, part I: failure surfaces[END_REF], the first order of expansions (power 0) provides the following conditions:

       P 0 ij,j + f ± i = 0 in Ω ± P 0 ij n j = ḡ± i on S± g ū0 i = 0 on S± u S0 ij = ā± ijhk Ēhk (ū 0 ) in Ω ± (33)
Using the results obtained in Section 2, we have that the following conditions

P 0 i3 = 1 2 b3333 (l) û2 2 û2 i (34)
and

η-1 l = ω-1 - 1 8 b 3333,l (l) û2 4 + ( 35
)
hold on the surfaces S ± . The latter conditions are complemented by the continuity conditions of the stress vectors along the surfaces S ± .

To summarize, the following system of equations is obtained:

                           P ij,j + f ± i = 0 in Ω ± P ij n j = g ± i on S± g u i = 0 on S± u S ij = a ± ijhk E hk (u) in Ω ± [P i3 ] = 0 on S P i3 = 1 2ε 2 b 3333 (l) |[u]| 2 [u i ] on S η l = ω - 1 8ε 3 b 3333,l (l) |[u]| 4 + on S (36)
The system [START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part II: in silico analysis[END_REF] governs the equilibrium problem of an adhesive interface characterized by a nonlinear damaging behavior, which is sandwiched between two adherents. In the next section, some numerical and analytical applications of the proposed interface model concerning modeling of the brick/mortar interfaces in masonry structures are investigated.

Mode I: uniaxial tension/compression condition

The proposed example is based on the former paper by R. Rizzoni et al. [START_REF] Rizzoni | On Saint Venant -Kirchhoff Imperfect Interfaces[END_REF] (Section 7). In the reference configuration, the composite structure is subjected to a tensile (compressive load) Q > 0 (Q < 0) directed along i 3 and acting on the upper and the lower boundaries Γ + and Γ -, as shown in Fig. 4. On the remaining part of the boundary ∂Ω 0 , the surface forces are taken to vanish. The load intensity Q is assumed to not depend on ε and the body forces are neglected. For the sake of simplicity, the adherents are supposed to comprise the same Saint Venant-Kirchhoff material with Lamé constants λ and µ, however the engineering elastic constants E and ν are used in what follows.

By combining the formulation proposed in [START_REF] Rizzoni | On Saint Venant -Kirchhoff Imperfect Interfaces[END_REF] with the one obtained in Sections 2-3 (cf. Eqs. ( 36)), the equilibrium problem of the joined structure reads as:

                                   DivP = 0 in Ω 0 P = (I + ∇u) ( E 1+ν )E(u) + νE (1+ν)(1-2ν) (I • E(u))I in Ω 0 E(u) = 1 2 ∇u + ∇u T + ∇u T ∇u in Ω 0 [P i 3 ] = 0 on S 0 P i 3 = 1 2ε 2 b 3333 (l) |[u]| 2 [u] on S 0 P i α = 0 with α = 1, 2 on ∂Ω 0 P i 3 = Q i 3 on Γ ± η l = ω - 1 8ε 3 b 3333,l (l) |[u]| 4 + on S 0 (37)
We seek the following solution for the displacement field u : Ω 0 → IR 3 :

u = (λ 1 -1)(x 1 i 1 + x 2 i 2 ) + (λ 3 -1)x 3 ± 1 2 [u 3 ] i 3 in Ω 0 ± (38)
with λ 1 ,λ 3 constants to be determined representing the stretches along x 1 and x 3 axes, respectively, and [u 3 ] the i 3 -component of the constant vector [u] representing the displacement jump at the interface S 0 . The first Piola-Kirchhoff stress tensor corresponding to [START_REF] Rekik | Homogenization methods for interface modeling in damaged masonry[END_REF] is:

P = P 11 (i 1 ⊗ i 1 + i 2 ⊗ i 2 ) + P 33 (i 3 ⊗ i 3 ) (39)
with

P 11 = E λ 1 2(1 + ν)(1 -2ν) (λ 2 1 -1) + ν E λ 1 2(1 + ν)(1 -2ν) (λ 2 3 -1) (40)
P 33 = ν E λ 3 (1 + ν)(1 -2ν) (λ 2 1 -1) + E (1 -ν)λ 3 2(1 + ν)(1 -2ν) (λ 2 3 -1) (41)
The stretch λ 1 , obtained by imposing the vanishing of P 11 has the following form (cf. [START_REF] Rizzoni | On Saint Venant -Kirchhoff Imperfect Interfaces[END_REF]):

λ 1 = 1 + ν -νλ 2 3 1 2 (42)
By substituting Eq. ( 42) in Eq. ( 41) and coupling the boundary condition on Γ ± (cf. Eqs. ( 37)), the equation which determines λ 3 as a function of the load Q is obtained:

Q = 2 E λ 3 λ 2 3 -1 (43)
Moreover, by recasting the fifth equation of the system (37) and taking into account the boundary condition on Γ ± , the displacement jump [u 3 ] is derived:

[u 3 ] = ε 2 3 Q |Q| 2 3 2 b 3333 (l) 1 3 (44)
which is positive if Q > 0 and negative otherwise.

In case of compression (Q < 0), ω is a negative-valued constant and therefore the last equation of the system [START_REF] Rekik | Identification of the representative crack length evolution in a multi-level interface model for quasi-brittle masonry[END_REF] gives η l = 0, implying no damage evolution. We thus consider the case of traction (Q > 0). After assuming b 3333 (l) = L 2 C l 2 , we have b 3333,l (l) = -L C l 3 , in which C is a positive-defined elastic constant (see Eq. ( 2)) and L is a constant playing the role of a characteristic interface length. By using this last expression as well as Eq. ( 44) in the last equation of the system (37), the damage evolution law as a function of the load Q is obtained:

η l ,Q Q = ω - C 2 ε L l(Q) 1 3 Q 4 3 + (45)
Because ω < 0, we find no damage evolution for 0 ≤ Q ≤ Q 0 with:

Q 0 = -ω 4 3 2 ε L l 0 C 1 4 with l(Q 0 ) = l 0 ( 46
)
Finally we reduce to the following differential problem:

     η l ,Q Q = ω - C 2 ε L l(Q) 1 3 Q 4 3 + l(Q 0 ) = l 0 (47)
whose solution is represented in Fig. 5 for the following values of the constants: C = 0.0014 MPa -1 , L = 210 mm, ε = 10 -2 L, l 0 = 0.05 mm, η = 1.7 e 9 MPa s, ω = -3.6 e 6 N/mm and Q = 1 MPa s -1 . Moreover, the displacement jump [u 3 ], shown in Fig. 6, is expressed as follows:

         [u 3 ] = ε 2 3 Q 1 3 4 C l 2 0 L 1 3 0 ≤ Q ≤ Q 0 [u 3 ] = ε 2 3 Q 1 3 4 C l(Q) 2 L 1 3 Q > Q 0 ( 48
)
with l(Q) solution of the differential problem (47).

Mode II: simple shear condition

In the reference configuration, the composite structure is subjected to a pure shear load q > 0 directed along i 1 and acting on the upper and the lower boundaries Γ + and Γ -, as shown in Fig. 4. On the remaining part of the boundary ∂Ω 0 , the surface forces are taken to vanish. The load intensity q is assumed to not depend on ε and the body forces are neglected. In analogy with the previous case, the adherents are supposed to comprise the same Saint The equilibrium problem of the composite structure is analogous to the system (37), except for replacing the boundary condition on Γ ± with

P i 3 = q i 1 . ( 49
)
We seek the following solution for the displacement field u : Ω 0 → IR 3 :

u = γ x 3 ± 1 2 [u 1 ] i 1 in Ω 0 ± , (50)
with γ a constant to be determined, representing the shear strain in the plane {x 1 x 3 } and [u 1 ] the i 1 -component of the displacement-jump vector [u] at the interface S 0 . The first Piola-Kirchhoff stress tensor corresponding to (50) is:

P = P 11 (i 1 ⊗ i 1 ) + P 22 (i 2 ⊗ i 2 ) + P 33 (i 3 ⊗ i 3 ) + P 13 (i 1 ⊗ i 3 ) + P 31 (i 3 ⊗ i 1 ) (51)
with

P 11 = E γ 2 4(1 + ν)(1 -2ν) , (52)
P 22 = E ν γ 2 2(1 + ν)(1 -2ν) , (53)
P 33 = E (1 -ν) γ 2 2(1 + ν)(1 -2ν) , (54)
P 13 = E 2(1 + ν)(1 -2ν) γ 1 -2ν 2 + γ 3 (1 -ν) , (55)
P 31 = E γ 4(1 + ν) . ( 56
)
If the boundary condition ( 49) is imposed, then the following equation which determines γ as a function of the load q is derived:

q = γG γ 2 1 -ν 1 -2ν + 1 2 , ( 57
)
being G = E 2(1+ν) the shear modulus. In analogy with the previous example, the following relationship holds:

[u 1 ] = ε 2 3 q |q| 2 3 2 b 3333 (l) 1 3 (58)
which can be rewritten as follows:

         [u 1 ] = ε 2 3 q 1 3 4 C l 2 0 L 1 3 0 ≤ q ≤ q 0 , [u 1 ] = ε 2 3 q 1 3 4 C l(q) 2 L 1 3 q > q 0 , (59)
with l(q) solution of the differential problem:

     η l ,q q = ω - C 2 ε L l(q) 1 3 q 4 3 + l(q 0 ) = l 0 (60) in which q 0 = -ω 4 3 2 ε L l 0 C 1 4
.

From a qualitative point of view, the damage parameter l(q) and the displacement jump [u 1 ] show a behavior similar to the one shown in Figs. 5 and 6.

Validation at the macroscale: Vermeltfoort and Raijmakers walls

The proposed numerical simulations made reference to the extensive experimental campaign carried out in the nineties by Vermeltfoort and Raijmakers [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] concerning shear tests on masonry walls. Vermeltfoort and Raijmakers results are well-known in structural analysis of masonry structures and they are usually taken as benchmarks for validating models [START_REF] Lourenço | Computational strategies for masonry structures[END_REF][START_REF] Alfano | Combining interface damage and friction in a cohesive-zone model[END_REF][START_REF] Zucchini | A micro-mechanical homogenisation model for masonry: Application to shear walls[END_REF].

Experimental tests [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] concern two one-meter-square and 100 mm thick walls, either without and with an opening, as sketched in Fig. 7. The central opening is one brick wide and six layers high and it is situated five layers from the bottom and one and half brick lengths from the loaded side of the specimen (see Fig. 7, left side). The walls comprise of full bricks of 210 * 100 * 52 mm 3 and head and bed mortar joints of 10 mm thickness. The up and the bottom sides of the specimens are connected via a high-performance mortar with two steel beams. In particular, the beam at the bottom side is connected to the test ring. Moreover, the up side beam is prevented to move in vertical direction by three vertical jacks in displacement controlled mode, and it is forced to move horizontally by a jack placed on the left side of the stiff beam. The horizontal load acts at the upper edge of the specimen as an uniform distributed shear load [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF]. Vermeltfoort and Raijmakers tested many specimens comprising units and mortars of different material properties, nevertheless, only the JD-type walls [START_REF] Lourenço | Computational strategies for masonry structures[END_REF] Both specimens (with and without opening) are subjected to a vertical preload of 0.3 MPa, then the shear force is imposed increasingly and monotonically via the horizontal jack. Concerning the wall without opening, some horizontal cracks caused by traction stresses and located at the application corner of the force and at its opposite in the bottom side, propagate rapidly when the horizontal load is applied. Successively, cracks propagate following the mortar joints path along the compressed diagonal with the load increasing, until the rupture. Moreover, the application of the horizontal load to the wall with opening causes first cracks nucleation and propagation localized at the corner of the central opening. Consecutively, after load increasing, a predominant diagonal crack is formed along the compressed diagonal. Finally, the rupture mechanism consists in some localized compressive failures at the top and at the bottom of the wall.

Numerical modeling of the walls

The proposed numerical simulations have been carried out with the FEMbased software Comsol Multiphysics. The vertical pre-loading is not taken into account for this study, accordingly, only the applied shear force is simulated. The problem is modeled in 2D plane stress conditions. The same boundary conditions are assumed for both walls, as sketched in Fig. 8. In particular, the bottom side beam is assumed to be fixed to a rigid plane. Moreover, the top side beam is constrained in the vertical direction and a horizontal displacement is applied on its left side. The simulations comprise a quasi-static loading in displacement controlled mode, accordingly the applied horizontal displacement is linearly increased until its maximal value.

The three principal domains, i. e. bricks, mortar and steel beams, are modeled as linear elastic isotropic materials whose material properties are taken in according with experiments by Vermeltfoort and Raijmakers [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF], as detailed in the previous paragraph.

It is worth noting that the proposed finite element analyses are carried out within the small strain theory framework. This fact is not incompatible with the implementation of an interface law obtained via the finite strain hypothesis, as demonstrated in a previous work by authors [START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part I: theory[END_REF]. A free triangular mesh of quadratic Lagrangian order is used in whole domains of both models, steel beams included (Fig. 8). The simulations are solved within a quasi-static framework via a temporal discrete solver included in Comsol. This choice has been made in order to limit solution convergence problems, which should be due to the strong nonlinearities of coupling brick/mortar interface laws with the evolution law of l, sixth and seventh equations of system [START_REF] Raffa | On modelling brick/mortar interface via a St. Venant-Kirchhoff orthotropic soft interface. Part II: in silico analysis[END_REF], respectively. Consecutively, a fixed time step is chosen as ∆t = 0.1 s and the shear test simulations are carried out within a range varying between 0 and 5 s. Accordingly, the imposed displacement increases linearly during this time range and it varies between 0 at the initial step and its maximal value d max = 5.5 mm at the final step.

Brick/mortar interfaces modeling

Within the present study the brick/mortar interfaces, either in horizontal and vertical principal directions of masonry walls, are assumed to be adhesive interfaces characterized by a nonlinear damaging behavior. These interfaces are governed, as a recall, by the following equations:

       P i3 = 1 2ε 2 b 3333 (l) |[u]| 2 [u i ] on S η l = ω - 1 8ε 3 b 3333,l (l) |[u]| 4 + on S (61)
where the first equation represents the soft nonlinear interface law and the second one the nonlinear evolution of the damage parameter l. These equations are specialized in the case of horizontal and vertical interfaces. The interface stiffness b 3333 (l) is obtained via a micromechanical homogenization of Kachanov-type [START_REF] Mauge | Effective elastic properties of an anisotropic material with arbitrarily oriented interacting cracks[END_REF][START_REF] Tsukrov | Effective moduli of an anisotropic material with elliptical holes of arbitrary orientational distribution[END_REF][START_REF] Rekik | Identification of the representative crack length evolution in a multi-level interface model for quasi-brittle masonry[END_REF][START_REF] Rekik | Homogenization methods for interface modeling in damaged masonry[END_REF][START_REF] Fouchal | An interface model including cracks and roughness applied to masonry[END_REF], its closed-form expression reads as:

b 3333 (l j ) = L j 2 C l 2 j (62)
where the subscript j = h, v indicates the horizontal and vertical interface cases, respectively. The compliance C is calculated from Eq. ( 2) starting from the undamaged interphase properties, and in the proposed models it results in C = 0.0014 MPa -1 (refers to [START_REF] Rekik | Identification of the representative crack length evolution in a multi-level interface model for quasi-brittle masonry[END_REF][START_REF] Fouchal | An interface model including cracks and roughness applied to masonry[END_REF] for further details). Parameter L is the characteristic length of the interfaces, resulting in L h = 210 mm and L v = 52 mm for horizontal and vertical interfaces, respectively. The first derivative of Eq. ( 62) with respect the damage parameter l results in:

b 3333,l (l j ) = - L j C l 3 j (63)
The interphase thickness ε appearing in Eqs. ( 61) is specialized for horizontal and vertical interfaces and it is assumed to be a fraction of the characteristic length L, as follows: ε h = 10 -2 L h , ε v = 10 -2 L v . This last choice has been made in according with [START_REF] Lebon | Numerical study of soft adhesively bonded joints in finite elasticity[END_REF] in order to assure a good approximation of the interphase modeling problem with the physically-based interface one. Let indicate l(t k ) = l k and [u(t k )] = [u k ] with t k = k∆t, the damage parameter and the displacement jump vector calculated at the k-ieme time step, respectively. Similarly to the numerical approach in [START_REF] Orefice | An experimental/numerical study on the interfacial damage of bonded joints for Fibre-Reinforced Polymer profiles at service conditions[END_REF], by approximating the first derivative of the damage parameter with respect to time as l ≈ l k+1 -l k ∆t , an explicit form of the damaging equation (second equation of Eqs. (61)) is derived:

l k+1 j = l k j + ∆t η ω + L j 8 C ε 3 j (l k j ) 3 u k j 4 + j = h, v (64)
Eq. ( 64), where Eq. ( 63) has been taken into account, is used for the numerical implementation of damaging in both horizontal and vertical interfaces of masonry walls. Moreover, parameter l has as an inferior limit bound, its initial value l 0 , so the following condition holds: l k ≥ l 0 = l(0). Model parameters l 0 = 0.05 mm, η = 1.7 e 9 and ω = -3 e 6 have been identified from experimental data by [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] concerning solid walls of JD-type, and they are used in both proposed simulations, i. e. solid wall and wall with a central opening.

As a result, the evolution of damage parameter l h in the solid wall is shown in Fig. 9. The l curve (Fig. 9) highlights a first path in which l is constant and equal to the initial value l 0 , then while load increases the damaging phenomenon begins. It is worth to remark that for the proposed numerical investigations a superior limit bound for the damage parameters is not taken into account, nevertheless it may be easily implemented. In the following section main results are discussed.

Results and discussions

Concerning the solid wall an identification process of the three model parameters l 0 , η and ω has been made. The obtained values allow to compare the global behavior of the simulated shear test with the experimental one [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF]. Particularly, Fig. 10 shows the horizontal imposed force against the horizontal displacement of the top side of the wall, either in numerical and experimental case. A good correlation is found, specially after the elastic limit. It is worth to remark that the fracture phenomenon is not object of the present study, which focuses on damaging and the consequent loss of global stiffness.

The initial damage parameter l 0 = 0.05 mm is assumed to be uniform in whole masonry wall at the initial step, this means assuming to have a small initial damage uniformly distributed at the brick/mortar interfaces, either horizontal and vertical. This pre-damaging can be caused, in the masonry manufacturing, by creep and shrinkage phenomena in the mortar. The parameter l 0 influences the initial elastic behavior of the wall. The model parameters η = 1.7 e 9 and ω = -3 e 6 govern the post-peak masonry behavior, i. e. after the elastic limit, by acting on the damage evolution (Eqs. (61)). Particularly, ω represents an energy threshold beyond which damage triggers.

The distributions of σ yy stresses shown in Fig. 11 are in agreement with the experimental findings by Vermeltfoort and Raijmakers [START_REF] Vermeltfoort | Shear tests on masonry panels of 1x1 m2[END_REF][START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF]. In details, these contour plots put in evidence, for solid wall as well as for the wall with the opening, the compressive and tensile zones at the top and bottom wall corner, and at the corner of the central opening. These stresses localization zones are the zones where experimentally cracks nucleate and begin propagating following the brick/mortar interfaces and/or the mortar joints. Within the numerical simulations, this last phenomenon is highlighted in Figs. 12 and13, where the distribution of damage parameters l h and l v is shown.

Particularly, highest values of damage are localized, accordingly with the experimental evidences, at the corner of the wall in the solid wall, and at the opening corner in the wall with the opening. The latter finding confirms that the proposed interface model is able to make a good prediction of damaging in masonry structures. Therefore, for both numerical models, parameter l v is smaller than l h . This result is related to the loading condition. In fact, by applying a shear loading in the horizontal direction the displacement jumps Fig. 11 Contour plots of yy-component of stress tensor with deformed shape at the final simulation step. Undeformed shape is traced also Fig. 12 Contour plots of the distribution of damage parameters l h and lv in solid wall in horizontal interfaces are higher than in the vertical interfaces, providing higher concentrations of damage (see Eqs. 61).

Conclusions

In the present study a model of soft imperfect interface comprising large displacements and evolving microcracking is derived via asymptotic techniques. The proposed interface model does not depend on loading conditions, on geometry and it is not mesh-dependent. Modeling the brick/mortar interfaces via the proposed imperfect interface model is a suitable and less computationally expensive alternative to existing detailed micro-models aimed at describing damage in masonry structures. Moreover, it allows to use linear elastic constitutive laws for the masonry's principal constituents, i.e., bricks and mortar joints.

A limitation of the proposed study is to not make a distinction between the damage level of compressive and tensile zones, as also highlighted in Figs. 12 and 13. As a future enhancement, an unilateral condition preventing damaging in compression will be accounted within the pseudo-potential of dissipation in order to better describe the general behavior of quasi-brittle materials.

Fig. 1

 1 Fig. 1 Sketch of the three body modeling approach

Fig. 2

 2 Fig. 2 An example of the geometry of the interphase

4 Numerical validations 4 . 1

 41 Analytical analysis of the brick/mortar interface model In what follows, two simple analytical applications of the proposed soft interface model are presented. The presence of I and II pure loading modes in the sense of fracture mechanics is investigated. The geometrical configurations are shown in Fig. 4.

Fig. 4

 4 Fig. 4 Sketch of the investigated systems composed of two nonlinear elastic isotropic adherents Ω 0 + and Ω 0 -joined by a nonlinear soft interface condition along a common face S 0 . I and II pure loading modes are represented on the left and on the right side, respectively

Fig. 5

 5 Fig. 5 Evolution of the damage variable l (Q) as a function of the load Q, obtained as the solution of the differential problem (47).

Fig. 6

 6 Fig. 6 Evolution of the displacement jump [u 3 ] as a function of the load Q.

  are considered in the proposed numerical simulations, accordingly bricks have a Young modulus E b = 16.7 GPa and a Poisson ratio ν b = 0.15 and mortar has a Young modulus E m = 0.782 GPa and a Poisson ratio ν m = 0.13.

Fig. 7

 7 Fig. 7 Geometry of the walls used in experimental tests by Vermeltfoort and Raijmakers.

Fig. 8

 8 Fig. 8 Schematic sketch of mesh and boundary conditions for the wall with opening

Fig. 9

 9 Fig. 9 Damaging in horizontal interfaces of the solid wall

Fig. 10

 10 Fig. 10 Horizontal force-displacement graphs for solid wall. Comparison between the proposed model and the experimental data by [48, 49] (black dotted curve)

Fig. 13

 13 Fig. 13 Contour plots of the distribution of damage parameters l h and lv in wall with opening

Acknowledgements This research was partially supported by Vinci Program 2013 (n. C2-73) of Italo-Francese University.

Conflict of interest

The authors declare that they have no conflicts of interest.


Ehsan Miandji

Jonas Unger

Christine Guillemot

Multi-Shot Single Sensor Light Field Camera Using a Color Coded Mask

Keywords:

We present a compressed sensing framework for reconstructing the full light field of a scene captured using a single-sensor consumer camera. To achieve this, we use a color coded mask in front of the camera sensor. To further enhance the reconstruction quality, we propose to utilize multiple shots by moving the mask or the sensor randomly. The compressed sensing framework relies on a training based dictionary over a light field data set. Numerical simulations show significant improvements in reconstruction quality over a similar coded aperture system for light field capture.

I. INTRODUCTION

Light field imaging has emerged as a promising technology for a variety of applications including photorealistic rendering, computational photography and computer vision.Many acquisition devices have been designed to capture light fields, ranging from camera arrays, [START_REF] Wilburn | High performance imaging using large camera arrays[END_REF], to single cameras mounted on moving gantries, [START_REF] Unger | Free Form Incident Light Fields[END_REF], and plenoptic cameras [START_REF] Adelson | Single lens stereo with a plenoptic camera[END_REF], [START_REF] Ng | Light field photography with a handheld plenoptic camera[END_REF]. Plenoptic cameras rely on multiplexing and use an array of microlenses placed in front of the sensor to capture multiple low resolution views in one 2D sensor image [START_REF] Ng | Light field photography with a handheld plenoptic camera[END_REF], [START_REF] Lumsdaine | The focused plenoptic camera[END_REF]. This is an efficient and easy to implement way of capturing multiple view-points, but reduces the spatial resolution by orders of magnitude compared to the a raw sensor image.

Recent research efforts have been directed towards overcoming the spatio-angular trade-off inherent to plenoptic cameras. In [START_REF] Liang | Programmable aperture photography: multiplexed light field acquisition[END_REF], Liang et al. used a programmable non-refractive mask placed at the aperture to sequentially capture the light field as a linear combination subsets of light rays. This series of 2D measurements, however, usually requires a large number of captured images and a long exposure times. Two architectures are instead described in [START_REF] Ashok | Compressive light field imaging[END_REF] for compressive acquisition of light fields, one exploiting correlations along the spatial dimensions and the other one along the angular dimensions, allowing for capture using short exposure times.

In [START_REF] Xu | A high-resolution lightfield camera with dual-mask design[END_REF], an architecture is proposed to compute optically coded projections on a fixed 2D sensor using two attenuation masks separately placed at the aperture plane and in front of the sensor. Given the measurements recorded on the sensor, the light field is then reconstructed using a least square minimization with a total variation regularization constraint. Similarly, the authors in [START_REF] Babacan | Compressive light field sensing[END_REF] place a randomly coded mask in front of the aperture to obtain incoherent measurements of the light field. Multiple shots are captured as random linear combinations of angular images by separately opening one region of the aperture and blocking light in the others.

A camera architecture is also proposed in [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] based on optically coded projections on a single image sensor. The light field is reconstructed using a compressive sensing framework, assuming that the light field is sparse in a domain defined by an overcomplete dictionary. In [START_REF] Miandji | Compressive image reconstruction in reduced union of subspaces[END_REF], incoherent spatioangular measurements are extracted using a random binary mask and the light field is reconstructed using an ensemble of 2D separable dictionaries. While the results significantly outperform [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] for general purpose light field sensing, the special structure of the sensing matrix required for a single sensor light field camera is not considered.

In this paper, we propose a camera architecture which captures incoherent measurements of the light field via a controllable color mask placed in front of the sensor. In order to increase the incoherence, we take multiple shots, where for each shot the mask configuration is changed to create a new random pattern. To reduce computations and increase the incoherence, we optionally also perform a random sampling of the spatial domain. The reconstruction algorithm is computationally efficient and existing 1 or 0 solvers [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF]- [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] can be used. In the analysis and evaluation we use RGB color filters, but the theory also applies to multi-spectral filters.

Throughout the paper we use the following notation. Vectors and matrices are denoted by boldface lower-case (a) and boldface upper-case (A) letters, respectively. Elements of a vector or matrix is denoted using subscripts, e.g. a i or A i,j . The jth column of a matrix is denoted A j . Finally, we denote a set of algebraic objects by superscripts; e.g. Φ i,j , i ∈ {1, . . . , m}, i ∈ {1, . . . , n}, denotes a set of mn matrices.

II. LIGHT FIELD COMPRESSIVE ACQUISITION

A popular approach for light field photography is coded aperture [START_REF] Liang | Programmable aperture photography: multiplexed light field acquisition[END_REF], [START_REF] Liang | Light field analysis for modeling image formation[END_REF]. By employing a geometric optics perspective, this method constructs several views of a scene using different regions of the aperture [START_REF] Babacan | Compressive light field sensing[END_REF]. Let the function l(r i , t j , u α , v β ) describe the two plane light field parametrization [START_REF] Gortler | The lumigraph[END_REF] defined by a pair of (discrete) locations on the sensor plane (r i , t j ) and the aperture plane (u α , v β ). The resolution of the sensor and the aperture plane are |r| × |t| and |u| × |v|, respectively. The coding of the aperture is achieved by using a mask of resolution |u|×|v| on the aperture plane. Each region of this mask produces an angular image on the sensor, which we denote by y i , where i ∈ {1, . . . , ν}, and

ν = |u||v|. If                  y R,1 y G,1 y B,1 y R,2 y G,2 y B,2 . . . y R,s y G,s y B,s                  y∈R ωλs =                  Φ 1,R,1 . . . Φ ν,R,1 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,G,1 . . . Φ ν,G,1 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,B,1 . . . Φ ν,B,1 Φ 1,R,2 . . . Φ ν,R,2 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,G,2 . . . Φ ν,G,2 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,B,2 . . . Φ ν,B,2 . . . . . . . . . . . . Φ 1,R,s . . . Φ ν,R,s 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,G,s . . . Φ ν,G,s 0 . . . 0 0 . . . 0 0 . . . 0 Φ 1,B,s . . . Φ ν,B,s                  Φ∈R ωλs×ωνλ                  x 1,R . . . x ν,R x 1,G . . . x ν,G x 1,B . . . x ν,B                  x∈R ωνλ , (1)
the mask is fully transparent, then the image formed on the sensor is the integrated radiance from all aperture regions, i.e. y = ν i=1 y i . Assuming the mask has values a i , the sensing model for a coded aperture can be formulated as follows

y = a 1 I a 2 I . . . a ν I      x 1 x 2
. . .

x ν      , (2)
where x i ∈ R ω , with ω = |r||t|, is the incoming radiance corresponding to the aperture region i, and I ∈ R ω×ω is the identity matrix. To maximize incoherence of the sensing matrix and a deterministic basis function, we would like the elements of Φ to be independent and identically distributed random variables [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF]- [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF]. Unfortunately, the sensing matrix a 1 I . . . a ν I in (2) is highly structured, and hence leads to coherent measurements. To alleviate this problem, Babacan et al. [START_REF] Babacan | Compressive light field sensing[END_REF] propose to use multiple shots, where for each shot a different set of mask values is used. Assuming s shots are taken, equation (2) becomes

     y 1 y 2 . . . y s      =      A 1,1 I A 1,2 I . . . A 1,ν I A 2,1 I A 2,2 I . . . A 2,ν I . . . . . . . . . . . . A s,1 I A s,2 I . . . A s,ν I           x 1 x 2 . . . x ν      . (3)
To further improve the incoherence, Marwah et al. [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] propose to place the mask at a small distance d m from the sensor, such that d m < d a , where d a is the distance of the aperture plane from the sensor. In this setup, the pixel value at location (r i , t j ) over the sensor is calculated by integrating the incoming radiance from different directions while modulating them by the mask:

y(r i , t j ) = |u| α=1 |v| β=1 f (r i + γ(u α -r i ), t j + γ(v β -t j )) l(r i , t j , u α , v β ), (4)
where γ = d m /d a . The function f (.) describes the mask value given a spatial location. Indeed there is a one to one correspondence between the pixels on the sensor and on the mask. The sensing model according to ( 4) is

y = Φ 1 Φ 2 . . . Φ ν      x 1 x 2
. . .

x ν      , (5)
where Φ i ∈ R ω×ω are diagonal matrices, each containing sheared mask values obtained by evaluating the function f (.) in ( 4) for the corresponding view direction (u α , v β ) and all the spatial locations over the sensor. The nonzero elements of the mask are typically drawn from an appropriate distribution (e.g. a nonzero Gaussian) or optimized for maximal incoherence with the sparsifying dictionary [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF]. In all the aforementioned acquisition frameworks, a monochrome mask is used. As a result, the same sensing matrix Φ is applied to all the color channels. This indeed leads to coherent measurements. In the following section, a new acquisition framework will be described that increases the measurement incoherence by randomizing the sensing of different spectral bands of the light field.

III. MULTI-SHOT COLOR CODED ACQUISITION

Unlike [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] that uses a monochrome mask, we utilize a random color mask. Recent advances in micro-lithography techniques [START_REF] Lapray | Multispectral filter arrays: Recent advances and practical implementation[END_REF]- [START_REF] Yi | Fabrication of multispectral imaging technology driven MEMS-based micro-arrayed multichannel optical filter mosaic[END_REF] has enabled the utilization of colored microfilters in many imaging applications such as multi-spectral imaging [START_REF] Arguello | Colored coded aperture design by concentration of measure in compressive spectral imaging[END_REF], depth extraction [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF], and image enhancement [START_REF] Lee | Color shift model-based image enhancement for digital multifocusing based on a multiple color-filter aperture camera[END_REF]. To further increase the measurement matrix incoherence, we use multiple shots, each with a different random colored mask. To achieve this, similar to [START_REF] Kittle | Multiframe image estimation for coded aperture snapshot spectral imagers[END_REF], we propose to use a piezo system for rapid mask movement. Depending on the spatial patch size (see Section IV) a movement of few pixels is enough for a new random mask pattern in each shot. Such systems are capable of capturing more than 60 frames per second, hence moving objects can be faithfully reconstructed.

The new light field sensing model incorporating per-shot colored mask patterns is formulated in [START_REF] Wilburn | High performance imaging using large camera arrays[END_REF]. Without loss of generality, and to simplify our notation, we consider three color channels (RGB). The entries 0 denote zero-valued matrices of appropriate size. Additionally, the sensing matrix Φ takes into account all the shots. Hence, the reconstruction TABLE I: Effect of sub-sampling on reconstruction quality. The angular resolution of the light field is 5 × 5 and we have created 9×9 non-overlapping patches over the spatial domain. The data set used here is described in section V. method, see [START_REF] Ashok | Compressive light field imaging[END_REF], can recover the light field x in a single pass. Alternatively, one can recover each shot individually and compute the final result as the average of all the recovered shots (see section V for a comparison of these approaches).

As it can be seen in ( 1), the total number of samples, i.e. the number of rows in Φ, is ωλs, where λ is the number of color channels. Since ω and λ are determined by the hardware implementation, the number of samples is mainly determined by s, which in turn affects the computational complexity of the reconstruction algorithm (see section IV). To provide flexibility over the trade off between computational complexity and reconstruction quality, we use spatial subsampling in the measurement model ( 1). This can be done by a sampling matrix P ∈ R rωλs×ωλs as follows

y = PΦx, (6)
where r ∈ (0, 1] is the sampling ratio. The matrix P can be constructed by sampling rωλs rows from I ∈ R ωλν×ωλν uniformly at random without replacement. Using [START_REF] Liang | Programmable aperture photography: multiplexed light field acquisition[END_REF], the total number of samples is proportional to sr. When s is fixed, e.g. due to hardware limitations, one can use r to reduce the reconstruction time (see Table II). On the other hand, spatial sub-sampling can be used to increase the reconstruction quality when the total number of samples, or the reconstruction time, is constant. To illustrate this, we consider two cases: 1. We take one shot without sub-sampling 2. We take two shots and set r = 0.5 so that the total number of samples is equal to the previous case. Results of this experiment are summarized in Table I. Using this simple approach we gain 0.24dB in Peak Signal to Noise Ratio (PSNR) without a significant increase in reconstruction time.

The nonzero values of the matrix Φ in (1) can be drawn from various distributions. A common approach is to use a nonzero centered Gaussian distribution |N (0, σ)|. Subgaussian distributions, e.g. Bernoulli distribution with probability 0.5, can also be used. We will compare the reconstruction quality for different distributions in section V.

IV. RECONSTRUCTION ALGORITHM

The reconstruction algorithm takes as input the measured light field (y) and the sensing matrix (PΦ) in order to obtain the original light field x, see [START_REF] Liang | Programmable aperture photography: multiplexed light field acquisition[END_REF]. Since in practice ν s, the sensing matrix is wide and there exists infinitely many solution to [START_REF] Wilburn | High performance imaging using large camera arrays[END_REF]. Using compressed sensing [START_REF] Candes | An introduction to compressive sampling[END_REF], the measurement model in (1) becomes x = PΦDα, where D ∈ R ωνλ×ρωνλ is an overcomplete dictionary with ρ as the overcompleteness factor,

we can recover or estimate the original light field as x = Dĉ. Equation ( 7) is known as Basis Pursuit Denoising (BPDN) and there exists several approaches for solving this problem [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF], [START_REF] Van Den | Probing the pareto frontier for basis pursuit solutions[END_REF], [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF]. In this paper we use the SL0 algorithm [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] since we observed superior reconstruction quality. Moreover, similar to [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF], we use a training based method [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF] to obtain an overcomplete dictionary. Due to the large size of the vectorized light field, x, it is a common practice to divide the light field over the spatial domain [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF], [START_REF] Miandji | Compressive image reconstruction in reduced union of subspaces[END_REF], while including all the angular and spectral information. As a result, the dimensionality of each vectorized light field patch will be ωνλ, where ω is the spatial patch size over the sensor. In addition, one can create overlapping patches [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] from the light field. After the reconstruction, the final light field elements are calculated by averaging the corresponding values in overlapped patches.

V. EXPERIMENTAL RESULTS

In this section we will report simulation results for our light field capture method. For the purpose of comparison, the synthetic light field data set of [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] is used. This data set consists of five light fields, which we have divided into training and testing sets as shown in Fig. 1. The angular resolution of each light field is 5 × 5, covering a field of view of 20 degrees; i.e. the angle separation is 4 degrees. The training set is used to learn a 1.5 times overcomplete, 10sparse dictionary. The spatial patch size is set to 9 × 9. Hence the signal dimensionality is 9 × 9 × 5 × 5 × 3 = 6075, where the last component is for color.

Figure 2 compares our method with that of Marwah et al. [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF]. Parameter values, along with PSNR and timing results, are summarized in Table II. We trained one dictionary for both methods. It should be noted that in [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF], a dictionary of size ων × ρων is trained and used for all the color channels independently. In our implementation of their method, we trained a dictionary of size ωνλ × ρωνλ. This dictionary improves the results of [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF] by about 3dB. We compared our results with the improved implementation of [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF]. Our method achieves about 3.6dB to 4dB improvement over [START_REF] Marwah | Compressive light field photography using overcomplete dictionaries and optimized projections[END_REF], II. using non-overlapping and overlapping patches, respectively.

Comparing figures 2b and 2f, we see that the grayscale mask leads to noise-like artifacts in the reconstruction. In contrast, our color coded mask produces results that are almost indistinguishable from the original light field.

When spatial sub-sampling is used, while keeping the number of shots fixed, we observe a drop of 1.5dB in PSNR. However, the computation time is reduced almost in half. Given the small reduction in visual quality (compare 2d and 2e), sub-sampling can be used to accelerate the reconstruction process. On the other hand, increasing the number of shots and decreasing the sub-sampling ratio while keeping the total number of samples constant will lead to a slight increase in reconstruction quality (see Table I). Moreover, we observed that independent reconstruction of shots followed by linear interpolation (Fig. 2c), performs inferior to simultaneous reconstruction (Fig. 2d) as described in section III.

The effect of the number of shots (s) on reconstruction quality is plotted in Fig. 4, along with a visual quality comparison in Fig. 3. It can be seen that when the number of shots is more than three, the PSNR grows linearly; note that a linear growth in PSNR implies exponential reduction in error. As a result, it is beneficial to increase s as far as the camera's ability to handle fast moving objects is not hindered. Results in figures 3 and 4 were obtained using non-overlapping patches and a Gaussian color mask.

Finally, we compare the reconstruction quality for different mask distributions in Table III. The following distributions are considered: uniform in the interval [0, 1], Gaussian and nonnegative Gaussian, Bernoulli distribution with probability 0.5, Rademacher distribution, and a quasi random sequence, specifically the Halton sequence [START_REF] Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF]. Note that the distributions that take negative values like the Gaussian cannot be implemented as a mask. However, we include these results for comparison since they are commonly used in compressed sensing. Indeed the Halton sequence is inferior to random masks. Moreover, nonnegative distributions are performing very similarly. Here we used non-overlapping patches.

VI. CONCLUSIONS

We introduced a compressive light field camera using a color coded mask. Moreover, a multi-shot capture by micro movements of the mask shows significant improvements in reconstruction quality. For future work, we are interested in determining the conditions for exact recovery that take into account the sparsity of the light field and the number of shots.

( a )Fig. 1 :

 a1 Fig. 1: Data sets used in our experiments: (a) Four light fields used for training an overcomplete dictionary, and (b) the light field we used in our simulations. and α defines the sparse coefficients. By solving the following optimization problem ĉ = min b b 1 s.t. p = ΦDb,(7)

Fig. 2 :

 2 Fig. 2: Visual quality comparison: (a) and (b) are obtained by the method of Marwah et al. [10] using non-overlapping and overlapping patches, respectively. And (c)-(f) show the results of our method using the parameters listed in TableII.

Fig. 3 :

 3 Fig. 3: Visual quality comparison for the number of shots (s). See Fig. 4 for the corresponding PSNR values.

Fig. 4 :

 4 Fig. 4: The Effect of number of shots (s) on PSNR. See Fig. 3 for a visual comparison. Non-overlapping patches are used.

TABLE II :

 II Reconstruction quality (PSNR) for various compressive coded aperture methods. For these results we used uniform distribution for nonzero values of (1).

		Mask	Method	s	r	overlap PSNR (dB) Time (min.)	Figure
	Marwah et al. [10]	Grayscale Grayscale	Simultaneous 5 Simultaneous 5	1 1	no yes	27.91 33.59	43 380	2a 2b
		RGB	Independent	5	1	no	29.23	30	2c
	Our method	RGB RGB	Simultaneous 5 Simultaneous 5 0.5 1	no no	31.50 30.14	43 24	2d 2e
		RGB	Simultaneous 5	1	yes	37.63	381	2f

TABLE III :

 III Reconstruction quality (PSNR) for various distributions used for non-zero values of (1). The zero-mean Gaussian distribution with variance σ is denoted N (0, σ).

	Uniform N (0, 1)	|N (0, 1)|	Bernoulli Rademacher	Halton
	30.72	31.50	30.80	30.87	31.50	27.70

ACKNOWLEDGMENT This work has been funded by Swedish Foundation for Strategic Research through grant IIS11-0081 and in part by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).