Ahmed Amari

Nadjib Achir

Paul Muhlethaler

Anis Laouiti

Optimized Broadcast Scheme for Mobile Ad hoc Networks

In this paper we propose an optimized broadcasting mechanism which uses very limited signaling overhead. The main objective is to select the most appropriate relay nodes according to a given cost function. Basically, after receiving a broadcast packet each potential relay node computes a binary code according to a given cost function. Then, each node starts a sequence of transmit/listen intervals following this code. In other words, each 0 corresponds to a listening interval and each 1 to a transmit interval. During this active acknowledgment signaling period, each receiver applies the following rule: if it detects a signal during any of its listening intervals, it quits the selection process, since a better relay has also captured the packet. Finally, we split the transmission range into several sectors and we propose that all the nodes within the same sector use the same CDMA orthogonal spreading codes to transmit their signals. The CDMA codes used in two different sectors are orthogonal, which guarantees that the packet is broadcast in all possible directions.

Introduction

Mobile ad hoc networks (MANETs) enable mobile nodes to communicate with each other without the need for a preexisting infrastructure. Since the radio range of each node is limited, multi-hop routing protocols are required to allow nodes that are not within range of each other to communicate. In this case, all the nodes need to act not only as simple nodes but also as routers to forward packets to the next hop towards the destination. Pointto-point communication is not the only communication mode, broadcast communication is also very important and is prominent in many applications. Moreover, broadcast is widely used for a large number of routing protocols proposed in the literature. Indeed, many of existing routing protocols use the easiest broadcasting scheme, the simple flooding algorithm, to construct the routing tables, such as in OLSR protocol. Unfortunately, many studies pointed out the fact that pure flooding may lead to many contentions and collisions, and thus serious performance degradations. This situation is known as the "broadcast-storm" problem.

In order to cope with these issues many challenges need to be faced to provide efficient broadcasts over mobile ad hoc networks. In this paper we focus on one of the main challenges, which is how to select the most appropriate forwarding nodes. Basically, classical flooding scheme can be improved if only a few nodes forward the messages sent by the initiator of the broadcast. More precisely, the main challenge is to find and select the fewest number of nodes to forward a message while maximizing the coverage. This problem has been defined in the literature as the problem of finding the minimum Connected Dominating Set (CDS). Unfortunately, finding this solution is NP-hard [START_REF] Khan | Application-driven cross-layer optimization for video streaming over wireless networks[END_REF].

Here we propose an optimized broadcasting mechanism which uses a very limited overhead. Our broadcasting algorithm uses an advanced acknowledgment scheme to select the appropriate relay nodes from all the potential relays within the transmission range. We propose that each node that captured a packet is required to transmit a short acknowledgment made up of signaling bursts. This sequence is computed by each node according to a given cost function. The acknowledgment is used to select the best potential relay nodes and consists of a sequence of intervals of the same length in which a given receiver can either transmit or listen. During this active acknowledgment signaling period, if a receiver detects a signal during any of its listening intervals, it quits the selection process, since it indicates that a better relay has also captured the packet. We thus obtain a powerful selection process to determine the most appropriate relay nodes. Finally, to distribute the packet in all possible directions, we propose that the transmission area is split into several sectors. Moreover, to guarantee that one relay is selected in each sector, a different CDMA spreading code is used in each sector. In two different sectors CDMA spreading codes are orthogonal.

The rest of this paper is organized as follows: In section 2 we provide the state of the art. In Section 3 we detail the protocol we propose. Simulation parameters, metrics, and results obtained with our broadcasting scheme are presented in Section 4. Finally, in Section 5 we conclude the paper, and suggest future research directions.

State of the Art

Several mechanisms have been proposed in the literature in order to optimize the broadcast transmissions in MANETs, in terms of reliability, efficiency, energy consumption and ability to satisfy the application QoS constraints. Assessment of the proposed schemes have been conducted for several applications like video or audio streaming, file sharing, communications for military applications or for civilians through wireless sensor networks [START_REF] Tseng | The broadcast storm problem in a mobile ad hoc network[END_REF] [3] [4] [5] [START_REF] Peng | Ahbp: An efficient broadcast protocol for mobile ad hoc networks[END_REF]. Williams and Camp [START_REF] Williams | Comparison of broadcasting techniques for mobile ad hoc networks[END_REF] categorize the broadcasting schemes into four families: simple flooding, probability based methods, area based methods and neighborhood based methods.

The first family, simple flooding, is considered as the default solution for packet broadcasting in ad hoc networks. In this scheme, each node needs to retransmit each received packet to its neighbors once. Since packets are duplicated at each node simple flooding ensures that every node in the network receives each broadcast packet at least once. Moreover no signaling is needed. Unfortunately, when considering high density networks, Ni et al. [START_REF] Tseng | The broadcast storm problem in a mobile ad hoc network[END_REF] show that simple flooding may lead to serious message redundancy, a high level of contention and many collisions. This situation is often referred to the "broadcast storm" problem.

The second family, probability based methods, includes two broadcasting schemes [START_REF] Tseng | The broadcast storm problem in a mobile ad hoc network[END_REF]: probability based and counter based. In probability based scheme, each node is affected a predefined packet re-broadcasting probability: p. Whenever a node receives a packet, it will rebroadcast it with probability p, or drop it with a probability: 1 -p. Note, that if p=1 then this scheme leads to simple flooding. In the second scheme, counter based, each node rebroadcasts a given packet if the same copies of this packet have already been received less than a given number of times. Thus, compared to previous broadcasting schemes, the reduction in the number of redundant broadcast packets that could be obtained by the counter based scheme will rely on the network density to ensure the packet delivery. In a dense MANET, few nodes will rebroadcast packets, while with a sparse topology, a majority of nodes will have to rebroadcast the packets.

In the third family, area based methods, the retransmission of a packet is decided according to the estimation of the relative physical distance between the transmitter and its relay. Within this family, we can cite two schemes [START_REF] Tseng | The broadcast storm problem in a mobile ad hoc network[END_REF]: distance based and location based. In the first scheme, each node is given a predefined distance threshold (d ). When a node receives a broadcast packet it estimates its physical distance from the last relaying node. If the distance is less than the threshold, the packet is then drooped. Otherwise, the packet is retransmitted. The main problem with this scheme concerns accurately assessing this distance. For instance, the distance can be estimated from the measured signal strength of the received packet. Unfortunately, this technique does not operate well in a realistic ad hoc network. In location based schemes, each node adds its location information in the header of each packet sent. Thus, packet retransmission is done only if the additional coverage area added by the receiver is greater than a given threshold.

Finally, the fourth family contains neighborhood based methods which exploit information on the neighborhood of the nodes. The main schemes belonging to this family are: self pruning [START_REF] Lim | Multicast tree construction and flooding in wireless ad hoc networks[END_REF], dominant pruning [START_REF] Lim | Multicast tree construction and flooding in wireless ad hoc networks[END_REF], scalable broadcasting [START_REF] Peng | On the reduction of broadcast redundancy in mobile ad hoc networks[END_REF], multipoint relay [START_REF] Qayyum | Multipoint relaying: An efficient technique for flooding in mobile wireless networks[END_REF], Ad Hoc broadcast protocol [START_REF] Peng | Ahbp: An efficient broadcast protocol for mobile ad hoc networks[END_REF] and Simplified Multicast Forwarding for MANETs [START_REF] Macker | Simplified multicast forwarding for manet[END_REF]. In self pruning, each node piggybacks, into each retransmitted packet, its adjacent nodes. Thus, nodes that received this packet can check if their own adjacent nodes are the same as those piggybacked in the received packet. If so, the packet is dropped. Otherwise, it is re-transmitted. While self pruning exploits only one-hop node information, the remaining schemes dominant pruning, scalable broadcasting, multipoint relay, Ad Hoc broadcast protocol and Simplified Multicast Forwarding for MANETs require that all nodes have knowledge of their two-hop neighbors. Unfortunately, all neighborhood based methods require extra transmission overhead, especially in dense MANETs.

Recently Blaszczyszyn et al. [START_REF] Blaszczyszyn | Opportunistic broadcast in vanets (ob-van) using active signaling for relays selection[END_REF] propose an opportunistic broadcast protocol for VANET. This protocol is very similar to our approach since it uses an active signaling acknowledgment to select the best relay among all the vehicles that have received the packet. Unfortunately, this algorithm was proposed only in the case of one dimension corresponding to roads.

Proposed Algorithm

As introduced earlier, the problem we study in this paper is, given an input MANET, how to determine the most adequate relays with minimum/maximum cost, according to a given cost function, in order to provide an efficient broadcast strategy. To meet this requirement, our idea is, as for unicast communications, to use an acknowledgment scheme to select the appropriate relay nodes from all potential relays within the transmission area. More precisely, we propose that each node that captured a packet is required to transmit a short acknowledgment made up of active signaling bursts. The idea of active signaling bursts was initially introduced in the HiPERLAN type 1 standard [START_REF] Jacquet | Priority and collision detection with active signaling -the channel access mechanism of hiperlan[END_REF], where the objective is to merge priority contention and collision detection into an efficient channel access mechanism. Basically, each node willing to access the medium has to transmit access pattern of non-data signals in a preamble to the packet transmission attempts .

In HiPERLAN active signaling, each node must access the medium by transmitting an "access pattern" of non-data signals in the preamble to the packet. This access pattern basically consists of listening and signaling periods. The signal transmitted in the signaling period contains no data at all. The procedure of collision detection and conflict resolution is performed during the packet access preamble where the node alternates between listening and transmitting. The access rule in the access pattern is that when a node is in a listening period, if it detects some energy from another node, then it immediately stops the access scheme and defers its packet transmission for a later attempt. Otherwise the nodes transmit their packet at the end of the access pattern transmission. In this case, when a transmission attempt starts, each node in competition selects a new access pattern. Since this access pattern has a random component, the probability of two nodes having the same access pattern is reduced. The access patterns also depend on the access priority assigned to the packet by the MAC level data transfer service. Thus, it is clear that the node which has priority to transmit its packet is the one holding the pattern which presents the largest digital number expanded in binary, if we consider that one codes the pattern time slots with a 1 corresponding to a transmit period and 0 corresponding to a listening period.

Our proposed acknowledgment active signaling burst is quite similar to the HiPERLAN active signaling scheme, except that we propose to use it as an acknowledgment. We propose that each node receiving a packet has to transmit a short acknowledgment. In this case, several acknowledgments could be sent at the same time by all the nodes capturing the packet. In classical communication, this operation leads to several collisions between the acknowledgments, but not in our case. In fact each acknowledgment, denoted as a burst in the following, is not carrying any data at all but it consists of a sequence of intervals of the same length in which a given receiver can either transmit or listen. During the transmitting period each node sends a signal on the medium and during a listening period each node listens on the medium. These signaling bursts can be represented by binary sequences denoted in the following as burst codes. 0 denotes a listening interval and 1 denotes a transmitting interval. In the rest of this paper, we fix the size of the burst code to L bits We propose to use the following relay selection rule: if a receiver detects a signal (energy) during any of its listening intervals, it quits the selection process. This means that the receiver stops transmitting during the entire remaining part of the burst. We consider that the detection of a transmission during a listening interval implies that a better relay has also captured the data packet. When a node receives a data packet it computes its own burst code as a function of the criterion that we wish to optimize to select our "best" relay. The only condition that we need to satisfy is that the better the relay is, the higher the burst code is. For instance, if we consider that the best relays are the farthest nodes, we should code the distance separating the source node and the relay nodes in base 2; this distance in base 2 will be the burst code. Thus, we can easily check that the selection mechanism will always select the relay nodes having the highest burst code (i.e. the greatest distance), since there will always be an interval in which another relay with a smaller burst code (i.e. a smaller distance) listens when the relay with the highest burst code is transmitting. Finally, in order to discriminate between nodes using the same burst code (e.g. at the same distance from the transmitter) we add r bits randomly selected.

In Figure 1, we represent an example of the selection process of three potential relays, a, b and c. We can notice that node a, with a higher burst code wins the selection process. In this case, only node a will relay the packet. Unfortunately, considering only this selection mechanism does not guarantee that the packet will be propagated in all possible directions. In addition to the acknowledgment signaling mechanism, we propose that the transmission area of each node is split into a given number of sectors, noted as k. The objective is that at least one node in each sector has to relay the packet. On the one hand, if we consider the whole transmission range as one unique sector (k=1), the packet is relayed by only a few nodes. This may negatively impact the performance. On the other hand, having a large number of sectors will also affect the performance by increasing the number of packet collisions. The number of sectors should not be too small nor too large. Moreover, since we want at least one node to act as a relay node in each sector, we also need to guarantee that nodes belonging to adjacent sectors will not eliminate each other. In other words, the relay node selection process should act only between nodes belonging to the same sector. To achieve this requirement, one possible approach is to perform a sequential selection process phase, starting from the first sector up to the last sector. This approach will however increase the delay of the acknowledgment phase and thus reduce the network throughput. The solution that we propose is that during the transmitting period (corresponding to 1 in the burst code) each node must send a CDMA spreading code. The code that must be used by each potential relay is selected according to its sector.

CDMA orthogonal spreading codes are one of the major elements within the whole CDMA system. Essentially, the concept of CDMA is based around the fact that a data sequence is multiplied by a spreading code or sequence, which increases the bandwidth of the signal. Then, in the receiver, the same spreading code or sequence is used to extract the required data. Only when the required code is used, does the required data appear from the signal. In our case, we use the same principle, since each potential relay node has to multiply the burst code by the spreading code associated to its sector (denoted in the following as the sector code). If a potential relay is in listening mode, it must compute the autocorrelation between the received signal and its sector code. If the autocorrelation is equal to 1, then this means that another node, belonging to the same sector, is in its transmitting period. This implies that a better potential relay has also captured the data packet and the node must quit the selection process. Here, this node has to listen to the medium during the first interval and transmit during the second interval. We also assume that this node is in a sector having [1 1 1 -1 1 -1 -1 -1] as its sector code. As illustrated, the node multiplies the burst code to the sector code and transmits the result on the channel. At the same time, we assume that a second node (node b) having a burst code equal to [0 0] belongs to the same sector as a. In this case, b is listening during the two intervals. As we can see, node b receives a signal during the second interval and computes the autocorrelation between the received signal and the code of its sector. If the result obtained is equal to 1, which is the case in this example, than the node quits the selection process.

1 1 1 1 -1 -1 -1 -1 data 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 data Burst
In the second example, Figure 3, we consider the case where two nodes belonging to two adjacent sectors receive the acknowledgment signals. As in the previous example, we assume that node a has a burst code equal to [0 1] and a sector code equal to [1 1 1 -1 1 -1 -1 -1]. We consider here when we split the communication range in more then 2 sectors our approach gives the same delivery ratio as flooding which is close to 100%. Given that, in Figure 4.d we plot the broadcast delay. By broadcast delay, we mean the delay required for one packet to reach all the network nodes. As expected, and for the same reasons as above, we notice that in the case of 1 or 2 sectors the flooding scheme exhibits smaller delays and outperforms our approach. However, when we increase the number of sectors our approach obtains better delays than flooding. In order to evaluate the impact of the broadcast algorithm on battery lifetime, we plot the consumed energy for each node for our broadcast and flooding schemes in Figure 4.b. We can see that the average energy consumed by the flooding scheme is much greater than the average energy consumed by our optimized broadcast scheme. Even if we increase the number of sectors to 32, our approach outperforms the flooding scheme. This is confirmed in Figure 4.c, where we plot the number of relays required.

1 1 1 1 -1 -1 -1 -1 data 1 1 1 1 -1 -1 -1 -1 1 1 1 -1 1 1 1 1 -1 -1 -1 -1
In addition, we analyzed the performances of our scheme when we increase the traffic load. In this case, we fix the number of broadcast sources to 1%, 2%, 4%, 6%, 8% and finally 10% of the total number of deployed nodes. Also, we consider that each source is generating a CBR packet traffic with an average packet inter arrival time varying from 100ms to 400ms, with 5, the delivery ratio achieved by our proposal consistently outperforms the delivery ratio obtained using flooding. Finally, in Figure 6 we plot the per packet energy consumption. As we can see, the optimized broadcast scheme consumes less energy than the flooding scheme.

Conclusion

In this paper, we present an efficient broadcast strategy for MANETs. We consider a cost function for packet transmissions and we we study how to determine the best relays for this cost function. Our idea is, as for unicast communications, to use an acknowledgment scheme to select the appropriate relay nodes from all the potential relays within the transmission area. We split the transmission range into several sectors and within each sector, all the nodes that have captured the broadcast packet are required to transmit a short acknowledgment made up of active signaling bursts. This burst code consists of a sequence of intervals of the same length in which a given receiver can either transmit or listen. This acknowledgment phase guarantees that the potential relays that have the highest burst code are selected as relay nodes. Moreover, in order to guarantee that nodes belonging to adjacent sectors do not eliminate each other, we propose that during the transmitting interval (corresponding to 1 in the burst code) each node has to send CDMA orthogonal spreading codes. The code that must be used by each potential relay is selected according to its sector. Finally, we evaluated our approach and we compared it to the classical flooding scheme. The results obtained demonstrate that our proposal outperforms flooding by increasing the delivery ratio and decreasing the number of required relays and thus the energy-cost.

Figure 1 :

 1 Figure 1: Selection process

Figure 2 :

 2 Figure 2: Selection process: two nodes in the same sector

Figure 4 :

 4 Figure 4: Delivery ratio, delay, number of relays and energy consumption when transmitting only one packet

Figure 6 :

 6 Figure 5: Delivery ratio

Code Burst Code Channel output Sector Code Sector Code Node a Node b Code Burst Code Channel output Sector Code Sector Code Node a

that node b has a burst code equal to [0 0] but a sector code equal to [1 -1 1 1 1 -1 1 1] since it is not in the same sector as a. We can clearly notice, from Figure 3, that node b does not receive the data sent by node a since the two nodes are using different sector codes. The autocorrelation between the received signal and the sector code of b is equal to zero. In this case, node b continues the selection process.

Performance Evaluation

In order to evaluate the performances of our proposed broadcast scheme we carried out a set of simulations. We implemented and integrated our broadcast schemes in the SENSE 3.0 simulator. In our simulations setup, we considered nodes following the IEEE 802.11 standard with a channel capacity of 1 Mb/s. We assumed an area of 1000 meters by 1000 meters where we deployed 100 nodes uniformly. We assumed that all the nodes are static. Moreover, we compared our approach with the classical flooding scheme. Finally, as the relay selection criterion, we considered a distancebased cost function, where the main objective was to force the farthest nodes to act as relay nodes.

We start our evaluation by analyzing the performances of our proposed scheme when varying the number of sectors from 1 to 32. In Figures 4.a,4.b, 4.c and 4.d, we plot the delivery ratio, the broadcast delay, the number of packet duplications and the energy broadcasting cost respectively. In Figure 4.a, we denote the packet delivery ratio (i.e. the ratio of the number of nodes taht received the packet) for both our optimized broadcast and pure flooding schemes. We study different numbers of sectors: 1, 2, 4, 8, 16 and 32 sectors. It is clear that, with only 1 or 2 sectors, the flooding scheme outperforms our approach. This due to the fact that we do not allow the packet to be propagated in all directions, which greatly limits the number of nodes relaying the packet and thus reduces the performances. However,


Fabien Momey
email: fabien.momey@cea.fr

Éric Thiébaut
email: eric.thiebaut@univ-lyon1.fr

Catherine Burnier
email: catherine.burnier@univ-st-etienne.fr

Loïc Denis
email: loic.denis@univ-st-etienne.fr

Jean-Marie Becker
email: jean-marie.becker@univ-st-etienne.fr

Laurent Desbat
email: lau-rent.desbat@imag.fr

Regularized 4D-CT reconstruction from a single dataset with a spatio-temporal prior

Keywords: Dynamic tomography, Reconstruction, Inverse Problems, Regularization, Signal processing. I

X-ray Computerized Tomography (CT) reconstructions can be severely impaired by the patient's respiratory motion and cardiac beating. Motion must thus be recovered in addition to the 3D reconstruction problem. The approach generally followed to reconstruct dynamic volumes consists of largely increasing the number of projections so that independent reconstructions be possible using only subsets of projections from the same phase of the cyclic movement. Apart from this major trend, motion compensation (MC) aims at recovering the object of interest and its motion by accurately modeling its deformation over time, allowing to use the whole dataset for 4D reconstruction in a coherent way.

We consider a different approach for dynamic reconstruction based on inverse problems, without any additional measurements, nor explicit knowledge of the motion. The dynamic sequence is reconstructed out of a single data set, only assuming the motion's continuity and periodicity. This inverse problem is solved by the minimization of the sum of a data-fidelity term, consistent with the dynamic nature of the data, and a regularization term which implements an efficient spatio-temporal version of the total variation (TV). We demonstrate the potential of this approach and its practical feasibility on 2D and 3D+t reconstructions of a mechanical phantom and patient data.

object". Reconstructing the patient's anatomy ignoring motion causes noticeable artifacts [START_REF] Balter | Uncertainties in CT-based radiation therapy treatment planning associated with patient breathing[END_REF], [START_REF] Shimizu | Impact of respiratory movement on the computed tomographic images of small lung tumors in three-dimensional (3D) radiotherapy[END_REF], [START_REF] Chen | Artifacts in computed tomography scanning of moving objects[END_REF]. In other words, the inverse problem of reconstruction has to be addressed in 4D; more specifically in 3D+t, where the motion and deformation of anatomical structures as a function of time are recovered. In image guided radiotherapy (IGRT), this is a critical matter because a X-ray treatment protocol requires a precise localization of the lung tumor, in order to prevent at best healthy tissues from irradiation. As a result, the "object" to be reconstructed is modelized as a spatio-temporal signal f (x, t), defined for spatial coordinates x ∈ R 3 and time t ∈ R.

The so-called 4D-CT reconstruction problem has benefited from active research for the past twenty years. Some investigations were made on the scanning protocol itself [START_REF] Goldberg | Evaluation of ultrafast CT scanning of the adult abdomen[END_REF], [START_REF] Ritchie | Minimum scan speeds for suppression of motion artifacts in CT[END_REF], [START_REF] Willis | Optimal scan for time-varying tomography. I. Theoretical analysis and fundamental limitations[END_REF], [START_REF] Willis | Optimal scan for time-varying tomography. II. Efficient design and experimental validation[END_REF]. Ritchie et al. [START_REF] Ritchie | Minimum scan speeds for suppression of motion artifacts in CT[END_REF] established early that artifacts were still present even with ultrafast scanning. Such an observation still holds twenty years later. A standard solution is based on a correlation between the 3D-CT data and a 1D temporal record of the patient's pseudo-periodical movements. The acquired projections are then sorted, according to the phase of the cycle to which they are related. Data subsets are extracted, and independent static reconstructions are then performed using each of these subsets. Such a method, called gated 4D-CT or phase-correlated 4D-CT [START_REF] Seppenwoolde | Precise and real-time measurement of 3D tumor motion in lung due to breathing and heartbeat, measured during radiotherapy[END_REF], [START_REF] Low | A method for the reconstruction of four-dimensional synchronized CT scans acquired during free breathing[END_REF], [START_REF] Ford | Respiration-correlated spiral CT: a method of measuring respiratory-induced anatomic motion for radiation treatment planning[END_REF], [START_REF] Vedam | Acquiring a four-dimensional computed tomography dataset using an external respiratory signal[END_REF], [START_REF] Keall | Acquiring 4D thoracic CT scans using a multislice helical method[END_REF], [START_REF] Berbeco | Towards fluoroscopic respiratory gating for lung tumours without radiopaque markers[END_REF], [START_REF] Sonke | Respiratory correlated cone beam CT[END_REF], [START_REF] Rietzel | Four-dimensional computed tomography: image formation and clinical protocol[END_REF], requires a sufficient number of projections for each reconstructed phase to avoid motion artifacts; as a consequence it increases the amount of X-ray dose delivered to the patient for the reconstruction of the whole sequence.

Another widely investigated class of methods aims at recovering the motion itself, i.e. models the deformation of the object as a function of time, and uses it to reproduce the motion. The motion is generally estimated in the form of a deformation vector field Γ t (x) that maps the volume at any time t back to its state at a reference time t 0 :

f (x, t) = f (Γ t (x), t 0 ) . (1)
Such a motion model can be incorporated into the tomographic projector in order to use all the available projections to reconstruct the patient's anatomy at this reference state. This class of methods is called motion compensation (MC). Research in this area has been active in the last decade, yielding very efficient algorithms. Such approaches rely on the quality of the estimation of the motion, a challenging inverse problem, especially because of the complex modeling of deformations such as breathing or cardiac beating. They generally require a separate high resolution 3D+t dataset for motion estimation.

In the IGRT context, 4D Cone-Beam CT (4D-CBCT) constitutes a particular modality because phase correlation is made difficult due to the poor number of remaining projections for the reconstruction of a single phase. One must turn to more robust approaches such as motion compensation. In 4D-CBCT, motion estimation is often done on an additional 4D-CT reconstruction [START_REF] Brankov | Motion-compensated 4D processing of gated SPECT perfusion studies[END_REF], [START_REF] Blondel | 3D tomographic reconstruction of coronary arteries using a precomputed 4D motion field[END_REF], [START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF], [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF], [START_REF] Rit | On-the-fly motioncompensated cone-beam CT using an a priori model of the respiratory motion[END_REF]. IGRT makes this prior information naturally available thanks to the necessary 4D planning CT prior to treatment. This so-called planning CT dataset is reconstructed by a phase-correlated method using a sufficient, i.e. enhanced, number of projections. Hence the total amount of data required for MC methods in this context is very high, meaning higher X-ray doses delivered to the patient. However, recent works have established that a very efficient MC reconstruction can be obtained from a single dataset, by a direct estimation of the motion model directly on the current measurements. In [START_REF] Brehm | Selfadapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy[END_REF], [START_REF] Brehm | Artifact-resistant motion estimation with a patient-specific artifact model for motioncompensated cone-beam CT[END_REF], the authors estimate the motion from a phase-correlated 4D-CT reconstruction on CBCT data, correcting the motion artifacts to recover accurately the deformation. Then they use the motion model to get the 4D MC reconstruction in a second step. In [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF], [START_REF] Yan | A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging[END_REF], the two inverse problems (motion estimation and reconstruction) are solved in a single algorithm using alternating schemes.

In parallel with the development of MC methods, phasecorrelated reconstruction approaches have also benefited from very interesting improvements. Following the emergence of the theory of compressive sensing, robust reconstruction algorithms were proposed, based on sparsity constraints, e.g. piecewise continuity, to allow efficient motion artifacts-free 3D reconstructions of each phase from a limited number of projections [START_REF] Chen | Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[END_REF]. However such methods still reconstruct each phase independently from each other. Only the most recent approaches propose to take into account the temporal correlation between the cyclic phases of the object under motion, extending the spatial prior knowledge to the temporal dimension. Their basic asumption is the continuity -or piecewise continuity -in both spatial and temporal dimensions, of the 4D image to be reconstructed. In particular, the approach proposed in [START_REF] Tian | Low-dose 4DCT reconstruction via temporal nonlocal means[END_REF] regularizes a 4D-CT reconstruction by temporal correlation of the different phases thanks to the minimization of nonlocal means between patches taken over the temporal dimension. Other methods [START_REF] Waibo | Spatial-temporal Total Variation Regularization (STTVR) for 4D-CT reconstruction[END_REF], [START_REF] Ritschl | Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatiotemporal sparsity prior[END_REF] use the Rudin et al. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] total variation prior in 4D, applying it on numerical 2D XCAT phantom CBCT-like data [START_REF] Waibo | Spatial-temporal Total Variation Regularization (STTVR) for 4D-CT reconstruction[END_REF] and in vivo cardiac micro-CT mouse data [START_REF] Ritschl | Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatiotemporal sparsity prior[END_REF].

Our approach follows the same principles as these spatiotemporally regularized 4D phase-correlated reconstruction methods: we use only the current dataset for the reconstruction of the 3D+t object, suppressing in this way any need for additional projections. No explicit knowledge of the motion is required and the whole 3D+t sequence is directly reconstructed, globally addressing the inverse problem of reconstruction from dynamic projections as the joint minimization of a data-fidelity term and a regularization term, as detailed in [START_REF] Balter | Uncertainties in CT-based radiation therapy treatment planning associated with patient breathing[END_REF]. The method is therefore very general and can readily be applied to many types of 4D-CT modalities. One of the novel contribution in our work is that we developed a global 4D tomographic projection model, based on classical static projection approach coupled with a temporal interpolation of the 4D object, allowing to calculate a projection at any given date t. Thus we overcome the classical phase correlation approach which only sorts the projections phase by phase with a given temporal gap. In addition, the specific continuity of the spatio-temporal object is enhanced by a 3D+t edge-preserving smoothness regularization module based on relaxed Rudin et al. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] total variation (TV).

In this introduction, we have presented the state-of-the-art of 4D-CT reconstruction, from the phase-correlated approaches to motion compensation, and presented the method we describe in this paper (cf. Section III). As explained in Section II, the demonstration of our method is made in the context of 4D-CBCT. In Section IV, we demonstrate its efficiency on 2D+t reconstructions from numerically simulated 2D+t data. Then we apply our method to real data, performing 3D+t reconstructions of a mechanical phantom and patient data, for which we confront our results to a standard MC reconstruction presented in [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF]. Applying our spatio-temporally regularized 4D phase-correlated reconstruction method to real human data acquired in the context of IGRT scanning protocol, is the other major contribution of our work, proposing to explore the feasibility of such an approach, in contrast to MC methods.

II. CONE-BEAM CT: DEMONSTRATION FRAMEWORK

A. Motion is critical in 4D-CBCT 4D-CBCT is an especially interesting modality for the study of 4D respiratory motion CT. Indeed, CBCT scanners have been implanted on linear accelerators for radiotherapy in the early 2000's [START_REF] Jaffray | A radiographic and tomographic imaging system integrated into a medical linear accelerator for localization of bone and soft-tissue targets[END_REF], [START_REF] Jaffray | Cone-beam computed tomography with a flat-panel imager: initial performance characterization[END_REF], [START_REF] Jaffray | Flat-panel cone-beam computed tomography for image-guided radiation therapy[END_REF]. The CT scan acquired during the treatment is for example used to control a posteriori the delivered therapeutic X-ray dose. This type of system has 2 specificities: the slow period of rotation of its gantry (from 1 to 2 minutes), and the fact that all the projections (about 600-700) are acquired on a single rotation of its flat detector. On the one hand, this slow rotation speed allows several respiratory cycles to be completed during the acquisition, giving a satisfactory angular coverage for each reconstructed phase. But, on the other hand, the lack of projections makes 4D-CT reconstructions difficult.

Many MC reconstruction methods have been proposed specifically for this type of data [START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF], [START_REF] Zeng | Estimating 3-D respiratory motion from orbiting views by tomographic image registration[END_REF], [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF], [START_REF] Rit | On-the-fly motioncompensated cone-beam CT using an a priori model of the respiratory motion[END_REF]. Li et al. [START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF] have incorporated the deformation field in the analytic reconstruction algorithm FDK [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. Rit et al. [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF], [START_REF] Rit | On-the-fly motioncompensated cone-beam CT using an a priori model of the respiratory motion[END_REF] have compared the Li et al.'s approach with an iterative reconstruction method based on the SART algorithm [START_REF] Andersen | Simultaneous Algebraic Reconstruction Technique (SART): a superior implementation of the ART algorithm[END_REF], inserting the deformation field into the tomographic projection model. These approaches assume that motion is unchanged between the acquisition on which motion is estimated and the CBCT acquisition for which motion is compensated. Zeng et al. [START_REF] Zeng | Estimating 3-D respiratory motion from orbiting views by tomographic image registration[END_REF] have proposed to avoid a preliminary 4D-CT reconstruction step by a direct estimation of the deformation model on the CBCT projections, using the complementary information given by an a priori 3D static reconstruction of the patient. The 3D anatomical model still has to be obtained from a previous reconstruction step, generally from breath-hold acquired data, a constraining step that generates additional irradiation. However, as already mentioned in Section I, very efficient MC reconstruction methods have been recently developed that perform the motion estimation and the 4D reconstruction by using only the CBCT dataset [START_REF] Brehm | Selfadapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy[END_REF], [START_REF] Brehm | Artifact-resistant motion estimation with a patient-specific artifact model for motioncompensated cone-beam CT[END_REF], [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF], [START_REF] Yan | A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging[END_REF]. Thus they overcome the approximations and possible errors induced by the hypothesis of reproducibility of the motion between two scans. Moreover they reduce the effective X-ray dose required for the reconstruction, a common objective with our approach.

B. "2D CBCT" dynamic data simulation

For testing purposes, we developed a simulated dynamic 2D phantom based on the well known Shepp-Logan model [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF], [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. The components of this phantom are ellipses whose position and shape may vary following a temporal sinusoidal signal with a given period. More precisely, these ellipses undergo changes (translations, rotations, dilations) that simulate the anatomical variations induced by the respiratory and cardiac motions at any instant t. Fig. 1 shows the phantom and the positions of its ellipses at the phase t = 0 of the motion cycle. The details of the motion features of each ellipse are summarized in Tab. I. The column Mot indicates the motion type : T for a translation, R for a rotation and D for a dilation/contraction. The shape parameters are : the center of the ellipse (x 0 , y 0 ) (in cm), the half-axis lengths a and b (in cm), the tilt from the horizontal axis of the major half-axis α (in • ), and the attenuation coefficient µ given as a factor of the attenuation coefficient of water (in cm -1 ). The amplitude A of the motion type is attached to the corresponding shape parameter of the moving ellipse.

Note that some ellipses are kept static (motion type "-").

The effective motion is given by the variation of the shape parameter following a temporal periodic sine or cosine signal (column Sig). The respiratory motion is simulated by giving to some of the ellipses a movement period of t a = 5 seconds.

The cardiac cycle is simulated by giving ellipse E8 a 1 second period. Hence the phantom can be calculated at any date t of a cycle. The attenuation value at a given spatial point is computed as the sum of all the ellipses levels of gray passing on it. Fig. 2 illustrates, on 3D views, the anatomical evolution of one motion cycle of this phantom. The observed waves in the temporal direction show the anatomical variations of the ellipses.

Fan beam projections are computed in close form under the hypothesis of an instantaneous measurement. This yields CBCT-like dynamic data in 2D. We simulated 600 such projections for a 360 • angle coverage corresponding to a total acquisition time of 120 seconds. As a result, 24 complete motion cycles are done during the acquisition, with 25 projections per cycle.

Fig. 3 displays a static reconstruction using the generated set of dynamic projections. The strong motion artifacts due to the inappropriate "static hypothesis" highlight the need for a 2D+t sequence reconstruction.

III. METHOD

A. Reconstruction criterion

We solve the inverse problem [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] of 3D+t reconstruction by minimizing the joint criterion:

f + = arg min f {J data (f ) + J prior (f )} , (2)
where J data and J prior are respectively the data-fidelity and regularization terms (detailed in the following paragraphs) and f denotes the parameters of the 3D+t sequence to be reconstructed, i.e. a 4D image of "spatio-temporal" voxels, corresponding to one period of motion.

B. Data fidelity and direct model

Assuming independent projections and Gaussian noise, the data fidelity can be written:

J data (f ) = θ∈Θ y θ -R θ •S t θ •f 2 W θ , (3)
where y θ denotes the measured projections at orientation θ of the detector and time t θ , Θ is the set of projection angles. The term R θ •S t θ •f is the model of the data, operator S t θ interpolates its arguments at time t θ while operator R θ performs the projection for orientation θ (as further detailed below). The discrepancy between the data and its model is measured by the squared Mahalanobis distance:

u 2 W θ = u • W θ • u
with W θ a weighting matrix equal to the inverse of the covariance of the noise.

As already mentioned, we assume that the acquisition of a projection is instantaneous. This hypothesis motivates the decomposition of the dynamic projection of the 3D+t object f in two successive operations: 1) The operator S t θ performs a temporal interpolation of the 3D+t voxels in order to extract the 3D object f (x, t θ ) at the projection time t θ . Our choice is justified by the fact that the temporal variation of the object is a continuous phenomenon. Hence it is possible to approximate the anatomical state of the object at any date by interpolating the frames sampled by the 4D voxels. Of course, the quality of the interpolation depends on the temporal sampling rate, as well as on the type of interpolation used. In this way we go further from classical 4D phase-correlated methods that simply divide the period of a cycle into equally spaced intervals and sort the projections. Fig. 4 shows the evolution of the root mean square (RMS) error between the data simulated in Section II-B and the output of our dynamic projection model as a function of the temporal sampling rate, for different interpolation methods.

As expected, the error decreases when the number of frames increases. In phase correlated 4D-CT, the selection of projections can be compared to using a nearest neighbor interpolation, which generates much more modeling errors than linear or B-spline-based interpolations as can be seen on Fig. 4. Oscillations in the RMS error curves (notably for the nearest neighbor interpolation) correspond to the sampling rates where some frames are exactly synchronized with the actual projections. In these particular cases, S t θ is just the identity operator and the remaining errors are due to the spatial sampling of f and to the approximations made by the projector R θ . In our implementation, we use linear temporal interpolation as it gives the least approximation errors (see Fig. 4). The coefficients of S t θ are determined from a 1D temporal signal giving the periodic evolution of the motion, similar to the signals used for 4D-CT reconstructions. It is easy from this signal to extract the period of a cycle, and to identify the phase of each temporal frame. Then the projection dates t θ give the locations of the interpolation points on the temporal axis, from which the coefficients of S t θ are deduced. Since the period of the respiratory or cardiac motion may slightly vary from one cycle to another, the sequence f has to be normalized to a mean period. Hence the cycles in the 1D signal are dilated or contracted in order to coincide with the normalized cycle. This permits to achieve a finer and more realistic calculation of the interpolation coefficients. A similar temporal registration method is used by Blondel et al. [START_REF] Blondel | 4D deformation field of coronary arteries from monoplane rotational X-ray angiography[END_REF], [START_REF] Blondel | 3D tomographic reconstruction of coronary arteries using a precomputed 4D motion field[END_REF] for cardiac SPECT.

In our tests, we made the same simplification as Blondel et al. [START_REF] Blondel | 4D deformation field of coronary arteries from monoplane rotational X-ray angiography[END_REF], [START_REF] Blondel | 3D tomographic reconstruction of coronary arteries using a precomputed 4D motion field[END_REF], i.e. we have neglected the variation of the spatial amplitude of the motion from one cycle to another. Our results on empirical data discussed in sections IV-B and IV-C support the use of this simplifying hypothesis.

2) The operator R θ is the tomographic projector which performs a static projection of the 3D object interpolated at t θ . R θ can be any numerical model implemented for iterative CT reconstruction. As already mentioned, projections must be separated according to the motion phase in order to reconstruct the temporal sequence. This implies a drastic reduction of the number of projections available per reconstructed temporal frame. In this context of reconstruction from a projections, we have shown [START_REF] Momey | A new representation and projection model for tomography, based on separable B-splines[END_REF], [START_REF] Momey | Spline driven: high accuracy projectors for tomographic reconstruction[END_REF] that the accuracy of the projector is a critical issue to make Fig. 5. A spatio-temporal representation of a situation where a same voxel P is "seen" at different times in 2 different neighbouring tissues, due to the motion of the interface between tissues A and B, point P is in tissue A at time t 1 , and in tissue B at time t 2 . The attenuation at point P then quickly changes between t 1 and t 2 when the boundary reaches P .

the best of available data. In order to cope with such cases, we have developed an accurate and fast model. It exploits B-spline basis functions to represent the 3D object and to approximate their projections. This model will be designated henceforth under the name spline driven [START_REF] Momey | Spline driven: high accuracy projectors for tomographic reconstruction[END_REF]. When reconstructing objects from a limited number of projections, we have shown spline driven to be more accurate than classical projectors such as distance driven

[13], with a satisfactory computional time. In order to use dynamic data in a better way, we therefore used the spline driven projector to implement our dynamic reconstruction approach.

C. Spatio-temporal regularization

The term J prior in equation ( 2) is essential to regularize the otherwise ill-conditioned inverse problem of dynamic reconstruction. In contrast to phase-correlated CT reconstructions that are performed independently, the regularization J prior can enforce a temporal coherence in the reconstructions. The spatial regularization must be chosen so as to favor smooth areas while preserving discontinuities at the boundary of anatomical structures. An edge-preserving smoothness prior such as TV [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is typically chosen in this context. Having a look at the 3D representation of our 2D+t simulated phantom displayed on Fig. 2, we observe that the variations of shapes are exclusively noticeable at the interfaces between the different parts of the object. For such an object, the high amplitude spatial gradients give the tissues' boundaries. They are therefore intrinsically correlated to the temporal gradients from which moving boundaries can be traced back. Fig. 5 illustrates the case where a given spatial position "sees" 2 different tissues at time t 1 and time t 2 . The temporal profile of attenuation sketched in Fig. 5 corresponds to a large motion where the change of the intensity value at position P is sharp. On the contrary, for sub-voxel motions, the sharp boundary will be more progressive and thus will give a smoother transition. Cardiac and respiratory movements also cause smooth spatial and temporal changes because they induce fluctuations of the absorption, particularly at the interfaces between different tissues, e.g. the lung walls or the myocardium's tissues. Consequently, we advocate the use of an 2 -1 edge preserving regularization [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF], [START_REF] Elbakri | Statistical image reconstruction for polyenergetic X-ray computed tomography[END_REF], [START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF], [START_REF] Wang | Iterative image reconstruction for CBCT using edge-preserving prior[END_REF] which is known to favor smoothness for small amplitude changes while preserving larger discontinuities. Being concerned with spatio-temporal (dis)continuity, we propose to use the following generalized form of relaxed TV:

J prior (f ) = k, 2 + µ 2 space ∇ space k, • f 2 + µ 2 time (∂ time k, • f ) 2 1/2 . ( 4
)
Here operators

∇ space k, et ∂ time k,
are respectively the 3D spatial gradient and the temporal derivative of the 3D+t object at the voxel (k, ) with k and the spatial index and temporal index (i.e., the frame) of the voxel. In order to tune the strength of the regularization (relatively to the data fidelity term) and to account for the heterogeneity of the physical dimensions of the spatial gradient and of the temporal derivative, we introduce two different hyperparameters µ space ≥ 0 and µ time ≥ 0 in Eq. ( 4). Note that the ratio of these two hyperparameters is homogeneous to a velocity. In Eq. ( 4), the parameter > 0 is not to be mistaken with a relaxation parameter introduced to avoid a singularity: it has to be appropriately tuned to set the trade-off between sharp and smooth changes. Note that our 4D regularizer is similar to those used in [START_REF] Waibo | Spatial-temporal Total Variation Regularization (STTVR) for 4D-CT reconstruction[END_REF], [START_REF] Ritschl | Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatiotemporal sparsity prior[END_REF].

Having a convex differentiable criterion (2), we perform the numerical minimization with the VMLM algorithm [27] which is a limited memory quasi-Newton optimization method with BFGS updates.

IV. RESULTS

A. Numerical data

To test our approach, we reconstructed our dynamic Shepp-Logan phantom presented in section II-B, from a simulated set of 600 projections. We chose to sample the 2D+t sequence with 25 frames. In this configuration, the frames are exactly synchronized with the projections which, as explained below, lets us mimic a phase-correlated reconstruction. Fig. 6 displays some frames of two different reconstructions from this dataset: 1) The sequence in the second column was reconstructed while imposing only the spatial regularization for each frame. This was achieved by setting µ time = 0 in Eq. ( 4). Thanks to the perfect synchronisation of the frames with the projections, the result corresponds to an ideal spatially TV regularized gated-like reconstruction, very close to the PICCS approach [START_REF] Chen | Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[END_REF]. The spatial hyperparameter µ space > 0 was tuned to achieve the best visual quality for the result.

2) The sequence in the third column was reconstructed using the spatio-temporal regularization in Eq. ( 4) with µ time > 0 and µ space > 0 both tuned to improve the visual quality of the result. To accurately compare the quality of the 2 reconstructions, Fig. 7 zooms into 2 regions of interest (ROI) centered at small and fine moving structures, which are very close from each other.

In the first case, each frame has been reconstructed independently from the others because there was no temporal correlation imposed between frames in this configuration. A set of only 24 projections is available per frame, which explains the rather poor quality of the reconstruction. This is particularly noticeable on the fine moving structures. Indeed we can observe on Fig. 7 a strong blurring effect due to the motion which makes the discrimination between the structures very difficult all along the cycle.

The second case, µ time = 0, shows the noticeable improvement brought by the temporal regularization of the 2D+t sequence in the quality of the reconstruction. The imposed continuity between frames helps to reduce artifacts due to the lack of available projections per frame. As a result, structures are better recovered, as well as their motion, in a nonambiguous way. The gain is particularly clear on the ROIs in Fig. 7 where the fine structures are well resolved and can be discriminated without ambiguity at each frame. As a result their respective motion is very well recovered.

B. Mechanical phantom data

We have reconstructed a mechanical phantom from a dataset acquired at the Centre Léon Bérard, on a scanner Elekta Synergy Cone-Beam CT. The flat panel is composed of a grid of 512 × 512 detector pixels of size 0.08 × 0.08 cm 2 . The acquisition process is similar to that described for CBCT scanners in section II-A.

The phantom has the features of a human thorax. Several zones of various densities mimic the lungs, the muscles and the spine. A small sphere of 2 cm in diameter is inserted in the right lung. It can be mechanically animated with a periodic circular motion in the transversal plane, and a periodic translation in the cranio-caudal direction. The two motions can be combined, reproducing the trajectory of a tumor during the respiratory cycle.

1) 2D+t reconstructions: A set of 630 projections has been acquired regularly in time and angle on 360 • during 116 seconds. The spherical insert was animated with a periodical motion in the transversal plane. The period of this "respiratory" cycle was 4 seconds. Thus, about 29 cycles are done during the acquisition, corresponding to about 22 projections per cycle.

The phantom was positioned so that the isocenter of the The corresponding reconstruction. The disk's interior corresponds to the part of the object that is projected, hence "visible", on all the frames. A remarkable feature of our inverse approach is its ability to produce a rather satisfactory reconstruction outside the "visible" part where there are some missing data.

scanner is contained in the same plane as the trajectory of the spherical insert. Extracting the median detector lines corresponding to this plane from the CB projections, we obtained a fan beam sinogram from which we reconstructed a 2D+t sequence of the corresponding phantom's slice. A single temporal frame of this reconstruction is shown in Fig. 8. The projections are clearly truncated by the detector. With iterative reconstruction methods, this problem can be partially addressed by reconstructing a larger volume of interest in order to encompass the whole object, the regularization helping to extend the reconstruction to regions with missing data. The static reconstruction displayed by Fig. 8 shows that, despite the truncation of the field of view, the object support is correctly restored even though there are some artifacts and a loss of details outside the fully seen region. However note that these defects have no impact on the quality of the reconstruction inside the fully seen region.

Our dynamic 2D+t reconstructions are presented on Fig. 9. As for the dynamic Shepp-Logan phantom, we have performed two reconstructions: a plain spatially regularized one (Fig. 9(d)), and a fully spatial and time regularized one (Fig. 9(e)). We chose to reconstruct a 2D+t sequence of 22 frames. The operator S t θ of our direct model (see Section III-A) is implemented by a linear temporal interpolation. Hence, as the direct model links each object frame to two projections, about 29×2 projections per frame are available for the reconstruction. For comparison, Fig. 9 also displays static reconstructions from a varying number of projections with the mechanical insert at still (a-c). Thanks to the spatio-temporal regularization, the reconstructed frames of the dynamic reconstruction in Fig. 9(e) have almost the same visual quality as the reconstruction from 630 projections in Fig. 9(a). This visual quality is in neat improvement over the static reconstruction from 64 projections shown in Fig. 9(b) or from the gatedlike reconstruction in Fig. 9(d). The motion of the insert is correctly recovered and the anatomical structures are better reconstructed, in particular the spine.

2) 3D+t reconstructions: Fig. 10 displays a 3D+t reconstruction of the phantom from another set of CB projections, acquired in the same conditions as previously. This time, the spherical insert was periodically animated with two kinds of motion: circular in the transversal plane and translational in the cranio-caudal direction. Hence the insert had a realistic 3D movement. We reconstructed a sequence of 22 frames made of 225×225×75 voxels of size 2×2×2 mm 3 . Outside the field of view, the reconstruction appears corrupted by motion and truncation artifacts. However, inside the field of view, the quality of the reconstruction is comparable to that of our 2D+t reconstructions, thanks to the spatio-temporal regularization.

To focus on the quality of the restored images, Fig. 10 only shows the non-truncated field of view. Transition from 2D+t to 3D+t involves no additional difficulties, but a bigger computional burden (many more variables, computations of 4D instead of 3D gradients in the regularization process).

C. Patient data

We have processed the same CBCT dataset of a patient's thorax as the one used in the article of Rit et al. [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF]. The CB projections were acquired at the Centre Léon Bérard, (Lyon, France), on a scanner Elekta Synergy Cone-Beam CT. Thus the acquisition process is the same as for the mechanical phantom in section IV-B. In [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF], the reconstructions were performed with a MC method developed by the authors and with a gated method. These reconstructions have been made available to us, in order to make a comparison with our method. We also obtained the recorded 1D temporal signal giving the periodic linear phase of the respiratory cycles, from which we have calibrated the temporal interpolator S t θ . The mean period of a cycle was evaluated to 2.4 s. We only kept the full cycles, reducing the number of required CB projections to 625. The duration of acquisition was 116 seconds, involving about 48 cycles, and 13 projections per cycle. We chose to reconstruct a 3D+t sequence of 13 frames. Using again linear interpolation for S t θ , about 48×2 projections were available per reconstructed frame.

Figure 11 displays the full 3D+t reconstruction. Each frame consists in 275×200×135 voxels of 2×2×2 mm 3 . Inside the field of view, the thorax is correctly reconstructed, and one can identify the different anatomical structures: the right lung containing a tumor, its bronchi, and also the heart and the ribs. The respiratory motion is well recovered, particularly that of the diaphragm and of the tumor, without ambiguities nor noticeable artifacts. For comparison, Fig. 12 shows side by side the MC and gated reconstructions by Rit et al. in [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF] and the result with our method. The end-exhale phase reconstructed by the MC method has a size of 262×261×132, i.e. a voxel's size of 0.98×0.98×2 mm 3 . The gated 3D+t sequence was reconstructed with the same parameters. We compare the reconstructed frame corresponding to the end-exhale phase. We have re-interpolated our reconstruction on the grid of finer voxels used by Rit et al., using a cubic B-spline interpolation kernel.

Regarding the gated reconstruction, our method has eliminated the motion artifacts and better recovered the anatomical structures. Indeed gated methods have to be applied on extended datasets (4D-CT) to be efficient, which is not the case for CBCT and explains the poor quality of the gated reconstruction. Note that this gated is basic, and that it could be improved using spatial regularization of each phase. However it evidences the gain of our approach which achieves a good quality of reconstruction using only the CBCT dataset.

The MC reconstruction reveals a finer resolution of recovered structures than ours, which looks smoother. Note that this MC method requires an accurate motion model computed from a high resolution 4D-CT data set. Our reconstruction has only needed the current set of CB projections to reconstruct the patient's thorax and the respiratory motion, recovering the motion as well as the static structures without ambiguity or motion artifact. As a result, a lower effective X-ray dose is necessary with our method than with this particular MC method. Indeed the latter method does not show the best result that could be obtained with the MC approach, since efficient MC reconstructions can be performed on the single CBCT dataset [START_REF] Brehm | Selfadapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy[END_REF], [START_REF] Brehm | Artifact-resistant motion estimation with a patient-specific artifact model for motioncompensated cone-beam CT[END_REF], [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF], [START_REF] Yan | A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging[END_REF]. However these results also demonstrate the potential of the spatio-temporally regularized 4D phase-correlated reconstruction approach, which is able to recover a motion artifacts-free 3D+t sequence with very simple and general hypotheses.

V. CONCLUSION This paper described an approach for phase-correlated 4D-CT reconstruction with a spatio-temporal regularization. We address the 4D reconstruction problem with very simple and general hypothesis. Our motivation was to perform 3D+t reconstructions from a limited number of projections, such as provided by a CBCT system, without resorting to additional 4D-CT data. Our results on CBCT data show that satisfying reconstructions can be obtained from limited data sets using a simple prior model based on spatial and temporal coherence.

Reconstruction of a 3D+t sequence from limited data is possible by considering jointly all projections in a global inverse problem. We carefully designed the model of pro- jections, based on spatial and temporal interpolation. Temporal frames are linearly interpolated to match phases of the cyclic movement under reconstruction. A joint spatio-temporal regularization exploits spatial and temporal coherence of the dynamic volume to provide images with smooth regions and motion. In order to preserve sharp boundaries and prevent from blurring motion, we proposed a relaxed 3D+t total-variation regularization. Such spatio-temporally regularized phase-correlated 4D-CT reconstruction method has already been explored in the context of numerical 2D XCAT phantom CBCT-like data [START_REF] Waibo | Spatial-temporal Total Variation Regularization (STTVR) for 4D-CT reconstruction[END_REF] and in vivo cardiac micro-CT mouse data [START_REF] Ritschl | Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatiotemporal sparsity prior[END_REF]. We applied our method to 4D-CBCT reconstruction problems from a mechanical phantom and real patient data. Our results both confirmed that our method can produce 3D+t reconstructions using only CBCT data, with a significantly improved quality compared to pure and even spatially regularized gated reconstructions. The motion of a lung tumor can be clearly identified in our reconstruction, which is confirmed by comparison with a reconstruction obtained using a standard MC method based on an accurate estimation of motion from high-resolution 4D-CT images. The reconstruction quality does not match that achieved using accurate motion models but may be a simple and general alternative when high-resolution 4D-CT datasets are not available.

By avoiding the use of additional 4D-CT data in dynamic CBCT reconstructions, the approach we followed should reduce the effective X-ray dose required for a given 4D reconstruction.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Temporal phase at t = 0 of our Shepp-Logan-like 2D dynamic phantom. It simulates a 2D object with 256 × 256 spatial voxels of size 0,1 × 0,1cm 2 , i.e. an global virtual size of 25,6 × 25,6cm 2 . It is composed of 15 ellipses Ei identified on the image. The values are shown in Hounsfield units, with µeau = 0,1928cm -1 .

5 Fig. 3 .

 53 Fig. 3. Static reconstruction from the set of dynamic projections generated by the simulated dynamic Shepp-Logan phantom.

Fig. 4 .

 4 Fig.4. RMS errors between simulated data of our dynamic Shepp-Logan phantom and the output of our dynamic projection model, R θ •S t θ , as a function of the temporal sampling rate. Each curve is representative of a given temporal interpolation method implemented by the operator S t θ .

Fig. 6 .

 6 Fig. 6. 2D+t reconstructions (25 frames) of the dynamic Shepp-Logan phantom (values are given in Hounsfield units). 24 projections are available per frame of the reconstructed cycle, with which they are synchronized (perfect spatially TV regularized gated-like case). 5 frames have been selected on a period-6 base (numbers 1, 7, 13, 19, 25). The first column shows the reference images. The second column shows the reconstructed images without temporal regularization (µ time = 0), i.e. the frames are independently reconstructed with a spatial regularization. The third column shows the reconstruction with a 3D (2D+t) spatio-temporal regularization (µ time = 0).

Fig. 7 .

 7 Fig. 7. Zoom into 2 regions of interest (ROI) of the 2D+t reconstructions of Fig. 6, showing small and fine moving structures of the phantom. Each of the 25 frames are shown for each ROI and compared to the truth with the root mean square error (RMS) map.

Fig. 8 .

 8 Fig. 8. 2D+t reconstruction of a mechanical phantom. Left: A sinogram extracted from the CB projections of the mechanical phantom. Right:The corresponding reconstruction. The disk's interior corresponds to the part of the object that is projected, hence "visible", on all the frames. A remarkable feature of our inverse approach is its ability to produce a rather satisfactory reconstruction outside the "visible" part where there are some missing data.

22 Fig. 9 .

 229 Fig. 9. 2D+t reconstructions of the median transversal slice of the mechanical phantom. Each restored sequence consists in 22 temporal frames with 500×500 pixels of size 1×1 mm 2 (the region of interest is zoomed in). (a,b,c): Static reconstructions from respectively 630, 64 and 29 projections of the static phantom (i.e., the insert was motionless). (d): Only spatially regularized reconstruction. (e): Spatio-temporally regularized reconstruction. The reference position of the spherical insert is identified on each frame by a small circle.

Fig. 10 .

 10 Fig. 10. 3D+t reconstruction of the mechanical phantom from the CB projections. A sliced view of each frame in the axial (center), sagittal (right), and coronal (top) directions is displayed. Only the non-truncated field of view is shown.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Our 3D+t reconstruction of the thorax of a patient (Elekta Synergy Cone-Beam CT data from [29], used by courtesy of Centre Léon Bérard). A sliced view of each frame in the axial (center), sagittal (right), and coronal (top) directions is shown.

TABLE I TABLE OF ATTRIBUTES

 IOF OF THE 15 ELLIPSES OF THE SHEPP-LOGAN-LIKE 2D DYNAMIC PHANTOM OF FIG. 1.

ACKNOWLEDGMENTS

The empirical data were kindly provided by Simon Rit from the Centre Léon Bérard (Lyon, France). Our algorithms were developped in Yorick [http://yorick.sourceforge. net/index.php], a freely available data processing language provided by David Munro. Computations for the 3D+t reconstructions were carried out on the Horizon cluster [http: //www.projet-horizon.fr].

This work was supported by the MiTiV project (Méthodes Inverses pour le Traitement en Imagerie du Vivant), funded by the French ANR (N • ANR-09-EMER-008).