Yohann Pitrey
email: ypitreyh@insa-rennes.fr

Yann Serrand
email: yserrand@insa-rennes.fr

Marie Babel
email: mbabel@insa-rennes.fr

Olivier Déforges
email: odeforge@insa-rennes.fr

ρ-domain for low-complexity rate control on MPEG-4 Scalable Video Coding

Scalable Video Coding was designed in response to the growing need for flexibility in video transmission over networks and channels. MPEG-4 Scalable Video Coding (SVC) is a recently finalized standard which introduces new coding tools such as spatial, temporal and quality scalability, to produce a layer-based scalable video stream. Additionally, inter-layer prediction allows a layer to use information from other layers as a basis for motion and texture prediction, improving the overall coding efficiency. Rate control is a capital issue in video coding, as it is designed to regulate the bitrate at the output of the encoder and keep it close to a specified constraint. Whereas rate control has been extensively studied for non-scalable video coding, only few propositions were made for scalable video coding. In this paper, we adapt an attractive rate control approach, based on a bitrate modeling framework called ρ-domain, for scalable video coding. We show that this model performs well on all spatial, temporal and quality scalabilities, and handles inter-layer prediction quite accurately. After validating the approach in MPEG-4 SVC, we use the ρ-domain model to build a simple-accurate rate control scheme. Results show that the mean frame bitrate error is below 7% on a representative set of configurations, while the impact on the complexity of the encoder is very low.

Introduction

In a world where video is omnipresent, with television, internet broadcasting and video-capable mobile devices, efficient video coding becomes more and more crucial. More than fifteen years of research have led to the release of the MPEG-4 AVC/H.264 standard in 2003 [START_REF] Wiegand | Overview of the h.264/avc video coding standard[END_REF]. Designed for video communications over networks, it overpasses the compression performance of previous standards such as MPEG-2 by nearly 50%.

Although very effective when addressing some given setup requirements, this standard is not suited to deal with heterogeneous networks and devices. Current video transmitting services have to face various decoding contexts, with various device screen sizes, computing abilities and network bandwidths. Existing solutions re-encode one specific video stream for each target set-up, causing storage and bandwidth wastage. The new MPEG-4 Scalable Video Coding (SVC) standard [START_REF] Reichel | Scalable video coding -joint draft 4[END_REF] provides more flexible video streams, which contain several layers. A base layer is first encoded, then improved using spatial, temporal and fidelity/quality enhancement layers. Spatial enhancement layers increase the frame resolution and address variable screen sizes. Temporal layers increase the number of frames per second and add smoothness to the video motion. Quality layers refine motion and texture information and increase the signal-to-noise ratio in the reconstructed video stream. The standard also introduces a new tool called interlayer prediction, which allows a layer to use information from the lower layers to predict motion and texture information. This way the redundancy between layers is removed, and the overall coding efficiency is improved.

Rate control is a capital issue when dealing with video transmission. It is intended to control the bitrate at the output of the encoder, to fit with a given bandwidth constraint while maximizing the quality of the reconstructed video stream. The quantization parameter (QP) determines the amount of information lost during the encoding process and has a direct impact on the bitrate and the quality of the reconstructed data. The main challenge concerning rate control is to find a way to predict the output bitrate from the value of the QP. Some approaches have been proposed for different video coding standards. In [5], a simple rate control scheme is proposed for the MPEG-2 standard. This early approach is based on a simplified rate-distorsion model that does not allow fine tunning of the encoder output. In [START_REF] Ma | Proposed draft description of rate control on JVT standard[END_REF], the authors extend the algorithm from [5] to the H.264 standard. In particular, following the original technique, a second optimization step is performed to fur-ther improve the rate-distorsion tradeoff. The method in [START_REF] Li | Adaptive basic unit layer rate control for JVT[END_REF] is more elaborate and enables better performance at the expense of a higher complexity. An interesting approach called ρ-domain was introduced in [START_REF] He | ρ-domain bit allocation and rate control for real time video coding[END_REF]. ρ-domain is a bitrate modeling framework, where ρ is defined as the percentage of zero-valued coefficients in a frame after quantization. It has been demonstrated that the relationship between ρ and the output bitrate is quite linear. This enables good bitrate prediction with very low complexity using ρ as an intermediate between the QP and the output bitrate.

Although rate control has been widely studied on single layer video, only few propositions were made for scalable video. In [START_REF] Leontaris | Rate control reorganization in the joint model soft reference[END_REF], the authors present some early work and focus on the regulation of the bitrate in the base layer. In [START_REF] Xu | Rate control for scalable video model[END_REF], the method from [START_REF] Ma | Proposed draft description of rate control on JVT standard[END_REF] is adapted for spatial and temporal scalabilities only. These two approaches fail to deal with all the features provided by MPEG-4 SVC. A complete rate control approach is proposed in [START_REF] Liu | Rate control of h.264/AVC scalable extension[END_REF]. The method deals with the three types of scalability and handles inter-layer prediction. However, the good performance of the method is at the expense of an increase of the complexity.

In this contribution, we adapt a ρ-domain based approach for scalable video coding. We first validate the approach in the context of MPEG-4 SVC by studying the relationship between ρ and the bitrate on the three types of scalability and on inter-layer prediction. Then, this approach is used to build a simple effective rate control scheme that performs well on spatial, temporal and quality scalabilities, while being able to handle inter-layer prediction. The effectiveness of the proposed algorithm is assessed on a representative set of set-ups proposed by MPEG-4 SVC. Results show that it controls the bitrate accurately, while the complexity increase is very low.

The paper is organized as follows. Section 2 introduces the ρ-domain bitrate modeling framework. In section 3, we describe the low-complexity rate control scheme derived from the ρ-domain model. In section 4, we illustrate the validity of ρ-domain for MPEG-4 SVC and discuss the performance of the proposed rate control scheme.

ρ-domain based rate model

Aside from classical rate control approaches, an original rate control framework called ρ-domain has shown interesting performances on single-layer DCT-based coders [START_REF] He | ρ-domain bit allocation and rate control for real time video coding[END_REF]. We will now describe this framework and extend it to scalable video coding.

ρ-domain for single-layer video coding

Let ρ be the percentage of zero coefficients in a frame after quantization. Due to the entropy coding that ends the encoding process, this parameter is directly related to the output bitrate of the encoder [START_REF] He | ρ-domain bit allocation and rate control for real time video coding[END_REF][START_REF] He | Linear rate control for JVT video coding[END_REF]. Whether a transform coefficient is coded as a zero or conserved after quantization depends on its position in the macroblock (MB) and the QP value. Let c m ij be the coefficient at position (i, j) in the macroblock m. This coefficient is coded as a zero if its value is below a certain dead zone threshold. Let z(c m ij , i, j, q) be an indicator function defined as follows:

z(c m ij , i, j, q) = 1 if c m ij < t(i, j, q) 0 otherwise, (1)
where t(i, j, q) is the value of the dead zone threshold for a coefficient at position (i, j) and QP value q [START_REF] Shin | Rate control using linear rate-ρ model for h.264[END_REF]. From 1, the value of ρ for a frame is computed as:

ρ(q) = 1 M m,i,j z(c m ij , i, j, q). ( 2
)
where M is the total number of coefficients in the frame.

An interesting point about ρ is that its relation with the bitrate can be approximated quite accurately by means of a linear model. This assumption has been validated in various standards such as MPEG-2, H.263 [START_REF] He | Optimum bit allocation and accurate rate control for video coding via ρ-domain source modeling. Circuits and Systems for Video Tech[END_REF] and MPEG-4 AVC/H.264 [START_REF] Shin | Rate control using linear rate-ρ model for h.264[END_REF]. In [START_REF] Shin | Rate control using linear rate-ρ model for h.264[END_REF], this relation is such as:

ρ(R) = R 0 -R * (1 -ρ 0 ) R 0 , (3)
where R 0 and ρ 0 are two parameters of the model which can be easily evaluated by encoding the frame once with an initial value of QP. Equations 2 and 3 allow us to find a relation between the bitrate and the quantization parameter, which is crucial for the design of a rate control scheme, as we show in the next section.

Validation of the ρ-domain for scalable video coding

The linearity of the relationship between ρ and the output bitrate has been demonstrated for MPEG-4 AVC/H.264 in [START_REF] Shin | Rate control using linear rate-ρ model for h.264[END_REF]. Although, the behaviour of ρ in the MPEG-4 SVC scalable context has never been studied. The new features such as the three types of scalability and inter-layer prediction imply some changes in the very nature of the bitstream. To be able to use a ρ-domain based approach in MPEG-4 SVC, we first need to verify that the relationship between ρ and the output bitrate is still linear.

We tested this assumption on each spatial, temporal and quality scalabilities. Each type of scalability was tested separately, first without inter-layer prediction, then with interlayer prediction. We measured the values of ρ and output bitrate for each frame and each value of QP between 0 and 51. Figure 1 shows a brief view of the typical relationship between ρ and the bitrate. It is quite visible that the relationship is still linear. the output bitrate for P and B frames in three layers for each type of scalability. There is a strong correlation between the two variables, which allows us to assert that their relationship is quite linear.

Our tests show that the linearity of the relationship between ρ and the output bitrate is still valid in the scalable context of MPEG-4 SVC. In next section, we use this observation to build a simple rate control scheme and show the interest of a ρ-domain based approach for MPEG-4 SVC.

Rate control scheme for SVC

The main objective of rate control is to regulate the output bitrate of the encoder. A rate control scheme classically contains two modules. First, a bit allocation policy is used to dispatch the available bitrate among groups of pictures (GOPs), frames or macroblocks. Then, achieving the target bitrate can be accomplished by acting on the quantization parameter (QP).

In our scheme, a target bitrate is specified for each layer separately. Each frame in a layer l is allocated the same target bitrate R t :

R t = T l F l , (4)
where T l is the target bitrate for the layer l and F l is the number of frames per second in the layer l. Once bit alloca-tion is performed, the optimal value of QP is computed for each frame. Assuming that the model in 3 is valid, there is a univoque relation between ρ and the bitrate. Therefore, for a target bitrate R t , there exists a unique value of ρ, say ρ t . From equation 2, there is a one-to-one mapping between ρ and the QP. Thus, finding the QP that generates the closest bitrate to R t is equivalent to finding the QP that generates the closest value of ρ to ρ t . This search is performed by computing equation 2 for each value of QP between 0 and 51. Computing the value of ρ t first requires to evaluate the parameters R 0 and ρ 0 in equation 3. To do this, we encode the frame with an initial value of QP, as mentioned in section 2. This rate control scheme was designed to be as simple as possible, to evaluate the performances of the linear rate model itself. It can still be greatly improved by including some more advanced features such as a wise bit allocation strategy and joint rate-distorsion optimisation.

Our rate model and rate control scheme both operate at frame-level. As displayed by figure 2, reaching ρ t would imply to consider a real QP value. Then, using only integer QP values introduces a restriction in terms of achievable bitrates. Refining the rate control at macroblock-level such as in [START_REF] He | ρ-domain bit allocation and rate control for real time video coding[END_REF] could alleviate this problem. Although, we wanted to build a low-complexity scheme, and controling the QP at macroblock level this would add large complexity and make it inapplicable, especially in a scalable video encoding context with several layers such as MPEG-4 SVC.

QP

ρ ρ R R t    q +1 t q t ρ t  Figure 2.
Frame-level rate control threshold effect. x the value of ρ t is processed, which corresponds to R t ; y ρ t corresponds to a float value of QP, which lies between two integer QP values; z choose q t as the value of QP which corresponds to the closest value of ρ from ρ t ; { the QP thresholding operation also yields a thresholding on the achievable bitrate.

Experimental results

We tested the rate control scheme proposed in Section 3 for each type of scalability enabled by MPEG-4 SVC. Our tests were performed on the JSVM Reference Software 8.6.

First, inter layer prediction is disabled, while a second set of tests evaluates the impact of this coding tool on the accuracy of the ρ-domain rate model. To evaluate the performance of the scheme, we measure the error between the target bitrate for a frame and the real number of bits produced by the encoder after rate control. As stated in Section 3, our frame-level rate-control scheme only allows us to reach a close value to the specified target bitrate for each frame. From equations 2 and 3, we compute the bitrate generated by encoding a frame with the chosen QP, say q t :

Rt = R 0 * (1 -ρ(q t )) (1 -ρ 0 ) . (5)
Then, the bitrate error percentage δ is computed at frame level from Rt and the real achieved bitrate R f :

δ = R f -Rt Rt * 100. ( 6
)
We also calculate the error percentage between the achieved bitrate and the original bitrate constraint:

∆ = R f -R t R t * 100. (7)

Spatial scalability

The generated bitstreams contain three spatial layers (QCIF, CIF and 4CIF) at 30 frames per second. Only dyadic spatial scalability is enabled. GOPs contain 4 frames (i.e.: PBBB). Each enhancement layer uses the lower layer as a dependency ID, and all quality IDs were all set to 0. The base layer is encoded at 100 kbits per second, while the two spatial enhancement layers are encoded at 300 and 900 kbits per second, respectively. The bitrate per second evolution is presented in figure 3. It can be seen that the bitrate constraint is respected over the three layers with very small error. Table 2 presents the frame error between achieved and target bitrates over two test sequences. µ δ and µ ∆ denote the mean value of δ and ∆ as defined in equations 6 and 7, while σ δ and σ ∆ denote their standard deviation. l 0 stands for base layer, l 1 and l 2 stand for the enhancement layers. Figure 4 shows achieved bitrate for each frame. The variations of the frame target bitrate Rt (dashed lines) we observe around the constraint R t (dotted lines) are due to the frame level of our rate control scheme, as stated in Section 3. The achieved bitrates are actually close to their targets, especially for P frames which have been marked with a dot on the figure.

Quality scalability

Quality scalability is quite similar to spatial scalability. Instead of changing the dimensions of pictures, residual information is refined. Our tests were run on three CIF layers at 30 frames per second with Coarse Grain Scalability (CGS). The same GOP structure as for spatial scalability was used. Dependency IDs are set to 0 and each enhancement layer use the previous layer as a quality ID.

Figure 5 presents the bitrate per second obtained for quality scalability. Our rate control scheme has a similar behaviour with CGS and spatial scalabilities, as confirmed by the results in table 3. Figure 6 shows achieved bitrates at frame level with quality scalability.

Temporal scalability

Temporal scalability is based on a hierarchical GOP structure. Each enhancement layer increases the number of frames per second of its base layer by two. This is performed by adding a B frame between each two frames from the base layer. Our tests on temporal scalability were performed on three CIF layers at 15 Hz, 30 Hz and 60 Hz. In the base layer, GOPs only contain P frames. In the first enhancement layer they contain one P frame and one B frame (PB-PB...), in the highest layer they contain one P frame and three B frames (PBBB-PBBB...). Results are shown in figure 7 andtable 4.

On temporal scalability, our ρ-domain rate control is quite accurate. The error is particularly small on the base layer and on the first enhancement layer. Figure 8 shows achieved frame bitrates for temporal scalability. It confirms the results from table 4, as the achieved bitrate is very close to the target for the two lower layers.

Inter Layer Prediction

We reproduced the previous tests over spatial, quality and temporal scalabilities, this time using inter layer prediction. Adaptive mode was enabled, which means that the mode decision and motion estimation choose the best prediction among macroblocks from either inside the frame, other frames in the same layer or in the base layer. We used the same configurations as before : the output streams contain three layers, encoded at 100, 300 and 900 kbits per second.

Figure 9 displays the achieved bitrate per second for spatial scalability with inter layer prediction. Frame allocation error is reported in table 5, while figure 10 presents the achieved frame bitrate. We observe that the error variation is made higher by the use of inter layer prediction. However, the mean bitrate per second remains relatively close to the constraint for each layer. Figures 11 to 14 and tables 6 to 7 display the achieved bitrates and frame errors for quality and temporal scalabilities with inter layer prediction. Our rate control scheme shows the same behavior for these configurations.

Results interpretation

Our rate control scheme achieves a mean frame error below 7% of the allocated bitrate on each type of scalabil- ity, including when inter layer prediction is enabled. This makes it quite suitable to fit with a bitrate per second constraint. We would like to report the particular effectiveness on P frames, for which the mean frame error below 3% of the target bitrate. On B frames though, the performances are not as good, with a mean frame error around 4% to 8%. Inter layer prediction causes the error variation to increase, but the mean frame error is not significantly affected. At frame level though, the error variation is relatively important as we can see in figures 4 and 10. This is caused by a well-known problem called the chicken and egg dilemma [START_REF] Li | Adaptive rate control for h.264. Image Processing[END_REF]. In MPEG-4 AVC/H.264 and SVC, the QP is involved in the choice of prediction modes and motion vectors. It means that the contents of the residual frame will vary from one QP to another, not only because of the quantization, but also because the QP affects the motion and texture prediction. As a result, if the target QP used for the second pass is too different from the QP used for the first pass, the bitrate prediction loses its accuracy. In figure 15, we can see that the linear rate model from 3 fits with the actual R(ρ) function with very small error. This further confirms that the linear ρ-domain based rate model from [START_REF] Shin | Rate control using linear rate-ρ model for h.264[END_REF] is quite efficient on MPEG-4 SVC. Although, there is a slight shift from the ρ and rate predicted by the model to their actual values. As an example, two points were marked on the figure. The first one (inside the circle) was predicted by the linear rate model. The first pass used a QP of 41, and a target QP of 37 was chosen as the best QP to fit with the target bitrate. The second point (inside a square), corresponds to the actual ρ and rate generated after encoding the frame with a QP equal to 37. The gap between these two points show the loss of precision of the rate model when the QPs are too different. Actually, the expected bitrate for this frame was 156 kbits, whereas the final output bitrate was 117 kbits, which is about 25% wrong.

This problem justifies that our rate control scheme is less performant on B frames and with inter layer prediction. Indeed, B frames can use more prediction modes and motion vectors than P frames, which means that changing the QP in a B frame is more likely to have a great impact on the coefficients. As well, inter layer prediction provides additional prediction modes, thus augmenting the possible combinations and perturbing the efficiency of bitrate prediction.

In this paper we studied the relevance of a ρ-domain based linear rate model in MPEG-4 Scalable Video Coding. After checking that the relationship between ρ and the bitrate was remains valid for scalable video, the model was incorporated into a simple and efficient rate control scheme.

Our tests on spatial, temporal and quality scalabilities showed that ρ-domain is an interesting approach for rate control in MPEG-4 SVC. We stated that the ρ-rate model manages to fit with the actual relationship with little imprecision on small QP variations. The mean frame error was below 7% of the allocated bitrate, and even lower on P frames. We also observed that inter layer prediction had a limited influence on the performance of this model.

At the moment, our rate control scheme is very simple. Each frame is allocated the same number if bits, whereas taking into account its type (I, P or B) and residual activity can help improving the smoothness of the reconstructed stream. Wise QP initialization should also be added to alleviate the precision loss issue we observed for large QP variations.

In our future work, the presented linear rate model will be incorporated into a more elaborate rate control scheme with frame-type dependant bitrate allocation, error beedback at both frame and GOP levels and distorsion considerations. Cross layer rate control will also be widely investigated. First, we will focus on a wise bitrate dispatching among the scalability layers, then the possibilities offered by inter layer prediction will be fully studied in terms of bitrate reusability and flexibility. We will also be interested in perceptual video quality measures over scalable video coding, which are better user-felt quality indicators.

Figure 3 .

 3 Figure 3. Bitrate per second for three spatial layers. From bottom to top : QCIF base layer; CIF spatial layer; 4CIF spatial layer. Constraint : • • • , achieved bitrate : -.

Figure 4 .

 4 Figure 4. Frame bitrate on HARBOUR for three spatial layers. R t : • • • , Rt : ---, R f : -.

Figure 5 .

 5 Figure 5. Bitrate per second for three CGS layers. From bottom to top: base layer encoded at 300 kbps; CGS layer encoded at 600 kbps; CGS layer encoded at 900 kbps. Constraint : • • • , achieved bitrate : -.

Figure 6 .

 6 Figure 6. Frame bitrate on HARBOUR for three CGS layers. R t : • • • , Rt : ---, R f : -.

Figure 7 .

 7 Figure 7. Bitrate per second for three temporal layers. From bottom to top : base layer at 15 fps; temporal layer at 30 fps; temporal layer at 60 fps. Constraint : • • • , achieved bitrate : -.

Figure 8 .

 8 Figure 8. Frame bitrate on HARBOUR for three temp. layers. R t : • • • , Rt : ---, R f : -.

Figure 9 .

 9 Figure 9. Bitrate per second for three spatial layers with inter layer prediction. From bottom to top : QCIF base layer; CIF layer; 4CIF layer. Constraint : • • • , achieved bitrate : -.

Figure 10 .

 10 Figure 10. Frame bitrate on HARBOUR for three spatial layers with inter layer prediction. R t : • • • , Rt : ---, R f : -.

Figure 12 .

 12 Figure 12. Bitrate per second for three temporal layers with inter-layer prediction. From bottom to top : base layer at 15 fps; temporal layer at 30 fps; temporal layer at 60 fps. Constraint : • • • , achieved bitrate : -.

Figure 13 .

 13 Figure 13. Frame bitrate on HARBOUR for three CGS layers with inter layer prediction. R t : • • • , Rt : ---, R f : -.

Figure 14 .

 14 Figure 14. Frame bitrate on HARBOUR for three temporal layers with inter layer prediction. R t : • • • , Rt : ---, R f : -.

  Figure 15. Linear approximation of the R(ρ) function.

Table 1 . Correlation coefficients between ρ and the bitrate.

 1 Table 1 further confirms this observation and shows the correlation coefficients between ρ and

		18					1.4
	Bitrate (bpp)	2 6 10 14				Bitrate (bpp)	0.2 0.6 1
		0.2	0.4	0.6	0.8	1	0.2	0.4	0.6	0.8	1
				ρ
			SPATIAL		TEMPORAL		QUALITY
			P		B	P	B		P		B
	layer 0	-0.9992	-0.9950	-0.9992	-0.9955	-0.9987	-0.9959
	layer 1	-0.9960	-0.9942	-0.9964	-0.9945	-0.9988	-0.9942
	layer 2	-0.9957	-0.9955	-0.9974	-0.9949	-0.9991	-0.9954
				USING INTER-LAYER PREDICTION
	layer 0	-0.9987	-0.9949	-0.9993	-0.9923	-0.9992	-0.9919
	layer 1	-0.9964	-0.9947	-0.9987	-0.9935	-0.9977	-0.9934
	layer 2	-0.9961	-0.9945	-0.9978	-0.9941	-0.9982	-0.9943

ρ Figure 1. Relationship between ρ and the bitrate (Left: P frame in the base layer for each value of QP. Right: spatial B frame in the highest enhancement layer with inter-layer prediction).

Table 2 . Target and constraint frame error for three spatial layers.

 2

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	-0.38%	-3.96%	-5.87%	-3.44%	-4.26%	-5.67%
	σ δ	12.53	11.19	7.37	6.79	7.42	8.80
	µ∆	-1.20%	-5.50%	-6.80%	-3.87%	-5.76%	-6.37%
	σ∆	9.58	8.95	6.02	10.01	10.46	12.24

Table 3 . Target and constraint frame error for three CGS layers.

 3

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	-0.29%	-1.48%	-2.61%	-2.99%	-4.98%	-4.10%
	σ δ	10.74	9.07	4.30	8.52	6.01	5.17
	µ∆	-0.90%	-0.84%	-1.00%	-4.51%	-5.02%	-4.90%
	σ∆	8.53	10.90	11.04	10.57	9.98	10.43

Table 4 . Target and constraint frame error for three temporal layers.

 4

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	2.05%	0.59%	-2.98%	2.77%	0.60%	-3.99%
	σ δ	3.21	4.31	8.56	9.21	9.15	6.06
	µ∆	-0.30%	-0.06%	-4.33%	2.27%	-0.04%	-4.40%
	σ∆	9.48	7.74	10.14	10.57	7.93	11.13

Table 5 . Target and constraint frame error for three spatial layers with inter-layer predic- tion.

 5

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	-0.38%	-0.03%	0.91%	-3.44%	-4.73%	-4.78%
	σ δ	8.53	13.11	14.19	6.79	9.71	11.29
	µ∆	-1.20%	-0.21%	-4.95%	-3.87%	-4.86%	-4.96%
	σ∆	12.58	14.52	15.99	10.01	12.30	14.63

Table 6 . Target and constraint frame error for three CGS layers with inter layer prediction.

 6

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	0.63%	-1.91%	-4.77%	-2.97%	-4.99%	-3.16%
	σ δ	15.71	14.83	6.23	8.27	11.27	5.60
	µ∆	0.57%	4.30%	-5.23%	-4.48%	-5.40%	-3.85%
	σ∆	16.12	16.66	13.70	10.63	13.27	9.77
		HARBOUR			SOCCER

Figure 11. Bitrate per second for three CGS layers with inter-layer prediction. From

  bottom to top : base layer encoded at 300 kbps; CGS layer encoded at 600 kbps; CGS layer encoded at 900 kbps.

Constraint : • • • , achieved bitrate : -.

Table 7 . Target and constraint frame error for three temporal layers with inter layer predic- tion.

 7

			HARBOUR			SOCCER
		l0	l1	l2	l0	l1	l2
	µ δ	2.05%	1.22%	1.58%	2.77%	0.08%	1.04%
	σ δ	3.21	13.31	16.17	11.21	6.85	16.69
	µ∆	-0.30%	5.05%	2.90%	2.27%	0.88%	1.52%
	σ∆	9.10	14.77	15.69	10.57	6.89	16.20


Angelo Alessandri

Patrizia Bagnerini

Christophe Prieur

Anna Rossi

Control of Normal Flow PDEs with ISS Properties

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Input-to-State Stability (ISS) was developed by E. Sontag in the late eighties [START_REF] Sontag | Smooth stabilization implies coprime factorization[END_REF], [START_REF] Sontag | Remarks on stabilization and input-to-state stability[END_REF], for systems described by ordinary differential equations, and has enabled the solution of numerous robust nonlinear control problems that were previously neither known nor addressed. No other foundational approach has reached comparable significance in the analysis of forced and interconnected nonlinear systems. Recently, a lot of efforts have been devoted to the extension of ISS to infinite-dimensional systems described by partial differential equations (PDEs). See recent textbook [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF] for an introduction of this notion for infinite-dimensional systems, and [START_REF] Mironchenko | Input-to-state stability of infinitedimensional systems: recent results and open questions[END_REF] for a survey with some results on nonlinear ones, mainly semilinear parabolic or nonlinear delay systems. Despite the large literature, most of such results concern systems described by linear hyperbolic and parabolic PDEs, while very little is known about nonlinear PDEs as considered in this paper.

In this work, the robustness of a control scheme proposed in [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF] for a class of nonlinear PDEs known as normal flow equations (NF, for short) is presented by showing its ISS properties in a suitable sense. NF equations are pretty popular in the research area of the so-called level set methods [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]- [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], which have been developed to describe the moving fronts for multi-dimensional applications such as fluid mechanics, image processing, and material science. Among other, the optimal control of heat equations is addressed in [START_REF] Hinze | Optimal control of the free boundary in a two-phase Stefan problem[END_REF], [START_REF] Bernauer | Optimal control of the classical twophase Stefan problem in level set formulation[END_REF], where the free boundaries are modelled by level sets. In [START_REF] Alessandri | Optimal control of propagating fronts by using level set methods and neural approximations[END_REF], the use of approximation is pursued for the purpose of optimal control involving level sets. Feedback controllers for noise-free systems described by NF equations are presented in [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF]. Modeling and identification of fire fronts by using level methods based on NF equations are considered in [START_REF] Alessandri | Parameter estimation of fire propagation models using level set methods[END_REF].

The stability of this class of nonlinear systems is in general nonrobust with respect to external vanishing disturbances, A. Alessandri, P. Bagnerini. and A. Rossi are with the University of Genoa (DIME), Via all'Opera Pia 15, I-16145 Genoa, Italy (emails: angelo.alessandri@unige.it; patrizia.bagnerini@unige.it; anna.rossi@unige.it).

C. Prieur is with the University of Grenoble Alpes, CNRS, Grenoble-INP, GIPSA-lab, F-38000, Grenoble, France (christophe.prieur@gipsa-lab.fr). His work has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003) as illustrated by a specific example in this paper (see the numerical results section). Therefore a dedicated control theory needs to be developed for such nonlinear PDEs. In this context, the main objective of this work is to derive a suitable control design method ensuring an ISS property of the closed-loop system, in presence of external disturbance and output measurement noise. To do that, we rely on the wellposedness framework of [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF] and extend the stability results presented therein by explicitly accounting for the presence of disturbances. It gives output feedback controllers so that the nonlinear closed-loop NF system is ISS for inputs defined by external disturbances and output measurement noise.

The importance of ISS in order to devise robust controllers for one-phase and two-phase Stefan models is highlighted in [START_REF] Koga | Materials Phase Change PDE Control and Estimation: From Additive Manufacturing to Polar Ice, ser. Systems & Control: Foundations & Applications[END_REF], where a number of real-world applications to processes with phase changes are presented. PDEs with moving boundaries have recently attracted a lot of attention for the capability to describe complex phenomena, which have been poorly investigated up to now (see also [START_REF] Koga | Towards implementation of PDE control for Stefan system: Input-to-state stability and sampled-data design[END_REF]). Clearly, the main difficulty is to find the appropriate ISS Lyapunov function under the most general assumptions [START_REF] Prieur | ISS-Lyapunov functions for time-varying hyperbolic systems of balance laws[END_REF].

The paper is organized as follows. Notation and preliminaries are described in Section II. The main result is presented in Section III. Section IV provides an example of asymptotically stable NF equation that is unstable in presence of a vanishing external disturbance. Moreover, applying the control design method of Section III yields an ISS closedloop system, as shown also by means of numerical simulations. Conclusion and prospect of future work are discussed in Section V.

II. NOTATION AND PRELIMINARIES

Given a set Ω ⊂ R q , Ω, ∂Ω, µ(Ω) denote its closure, boundary, and Lebesgue measure, respectively. The set C 0 (Ω) is the class of continuous functions in Ω, while C 1 (Ω) is the class of continuous functions in Ω with continuous first-order derivatives. BU C(Ω) denotes the set of bounded uniformly continuous functions on Ω. For any T > 0, consider the Hamilton-Jacobi equation

φ t (x, t) + H(x, t, φ(x, t), φ x (x, t)) = 0 in Ω×(0, T ) (1)
where Ω is a bounded, smooth (e.g., with Lipschitz boundary), open set; H : Ω × [0, T ) × R × R q → R is the Hamiltonian function; φ x denotes the gradient of φ w.r.t.

x; φ 0 (x) ∈ BU C(Ω) denotes the initial condition in Ω; boundary conditions are needed as well. The function φ 1 ∈ C 0 (Ω × (0, T )) is a viscosity subsolution of (1) if and only if there exists v ∈ C 1 (Ω × (0, T )) such that

v t (x 1 , t 1 ) + H(x 1 , t 1 , φ 1 (x 1 , t 1 ), v x (x 1 , t 1 )) ≤ 0 holds in (x 1 , t 1 ) ∈ Ω×(0, T ) local maximum point of φ 1 -v.
The function φ 2 ∈ C 0 (Ω×(0, T )) is a viscosity supersolution of (1) if and only if there exists v ∈ C 1 (Ω × (0, T )) such that

v t (x 2 , t 2 ) + H(x 2 , t 2 , φ 2 (x 2 , t 2 ), v x (x 2 , t 2 )) ≥ 0 holds in (x 2 , t 2 ) ∈ Ω×(0, T ) local minimum point of φ 2 -v.
If φ is a viscosity subsolution and a viscosity supersolution of (1), then φ is a viscosity solution of (1) (see [START_REF] Barles | Existence results for first order Hamilton Jacobi equations[END_REF]Definition II.1]). For any integer n ≥ 1 and

x ∈ R n , let |x| := n i=1 x 2 i . L 2 (Ω) denotes the Hilbert space of square integrable functions γ(•, t) : Ω → R n with norm |γ(•, t)| L2 := Ω |γ(x, t)| 2 dx 1/2 < ∞ for all t ≥ 0. H 1 (Ω)
is the Sobolev space of square integrable functions with square integrable first derivatives, i.e., λt) holds instead of the condition in (iii), i.e., without explicit dependence on |φ 0 -φ| L2 in the r.h.s., we simply say that φ converges exponentially to φ in the L 2 sense. Finally, the Young inequality is the upper bound of the cross product of any couple of real numbers a and b, as follows: 2ab ≤ a 2 + b 2 .

H 1 (Ω) := {γ ∈ L 2 (Ω) : γ x ∈ L 2 (Ω)}. Let φ ∈ H 1 (Ω) be an equilibrium of (1) if H(x, t, φ(x), φx (x)) = 0 for any t ≥ 0. Then, φ is said to be: (i) L 2 stable if for all ε > 0 there exists δ ε > 0 such that |φ 0 -φ| L2 < δ ε ⇒ |φ(•, t) -φ| L2 < ε for all t ≥ 0; (ii) L 2 asymptotically stable if it is stable and lim t→+∞ |φ(•, t) -φ| L2 = 0; (iii) L 2 exponentially stable if there exists λ > 0 such that |φ(•, t) -φ| L2 ≤ c |φ 0 -ϕ| L2 exp(-λt) for some c > 0 and all t ≥ 0. If |φ(•, t) -φ| L2 ≤ c exp(-

III. MAIN RESULTS

In this section, we consider the control of a multidimensional NF equation affected by disturbances. More specifically, we deal with

φ t (x, t) + f (x, t)|φ x (x, t)| = u(x, t) + w(x, t) , y(x, t) = φ(x, t) + v(x, t) , (2)
in Ω × (0, T ), where Ω ⊂ R q , T > 0, f is a fixed velocity field, u is the in-domain input, y is the output, w denotes the in-domain disturbance, and v accounts for output measurement noises. As pointed out in [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF], the knowledge of the gradient of φ would facilitate the design of a stabilizing controller for (2) in the absence of noises, given by u

(x, t) = -k φ(x, t) + f (x, t) |φ x (x, t)|.
Likewise in [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF], we do not assume the knowledge of the spatial derivatives of the solution and rely on a Luenberger observer for the second term in the left-hand side of (2) (i.e., f (x, t) |φ x (x, t)|); this observer provides an estimate of f (x, t) |φ x (x, t)| in such a way to apply the control given by

u(x, t) = -k y(x, t) + f (x, t) φx (x, t) , (3)
where k ∈ R is the gain to be chosen and φ(x, t) is the estimate of φ(x, t) given by the Luenberger observer

φt (x, t) + f (x, t)| φx (x, t)| + k( φ(x, t) -y(x, t)) = u(x, t) . (4)
The system based on the NF equation ( 2) entails the solution of the following boundary-value problem upon the application of a feedback law (this motivates the little abuse of notation concerning u):

               φ t (x, t)+f (x, t)|φ x (x, t)| = u(x, t, φ(x, t))+w(x, t) in Ω×(0, T ), y(x, t) = φ(x, t) + v(x, t) in Ω×(0, T ), φ(x, t) = ϕ(x) in ∂ Ω×(0, T ), φ(x, 0) = φ 0 (x) in Ω, ( 5
) where ϕ ∈ BU C(∂ Ω), φ 0 ∈ BU C(Ω) and f ∈ BU C(Ω × [0, T )
) is a known function acting as velocity field. Compared to (1), the Hamiltonian is given by H(x, t, s, p) = f (x, t)|p|-u(x, t, s)-w(x, t). The existence and uniqueness of viscosity solutions to (5) are guaranteed by the following assumptions:

w ∈ BU C(Ω × [0, T )), u ∈ BU C(Ω × [0, T ) × R), u nonincreasing w.r.t. s for each x ∈ Ω, t ∈ [0, T ), v ∈ BU C(Ω × [0, T ))
, as well as f, u , w and v Lipschitz w.r.t. x ∈ Ω. Thus, we can apply [17, Theorem 7.1, p. 82] since assumption (H1) is verified for f, u, w, v Lipschitz w.r.t. x ∈ Ω, and (H2) holds with γ R = 0, as u is nonincreasing w.r.t. s for each x ∈ Ω, t ∈ [0, T ). Concerning the observer (4), we have to deal with the problem

         φt (x, t) + f (x, t) | φx (x, t)| + k ( φ(x, t) -φ(x, t)) -k v(x, t) = u(x, t, φ(x, t)) in Ω × (0, T ), φ(x, t) = ϕ(x) in ∂ Ω × (0, T ), φ(x, 0) = φ0 (x) in Ω, (6)
for any T > 0 and with φ0 ∈ BU C(Ω), f and u satisfying the same assumptions considered for (5), φ 0 Lipschitz in Ω, and φ solution of [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF]. The problem (6) has a unique solution φ ∈ BU C(Ω × [0, T )) since (i) φ 0 is Lipschitz in Ω, the unique solution φ of ( 5) is Lipschitz in Ω × [0, T ) (see, e.g., [START_REF] Lions | Generalized Solutions of Hamilton-Jacobi Equations[END_REF]), and (ii) the Hamiltonian given by H(x, t, s, p) = f (x, t)|p| + ks -u(x, t, s) -k φ(x, t) -k v(x, t) satisfies the conditions of existence stated in [START_REF] Barles | An Introduction to the Theory of Viscosity Solutions for First-Order Hamilton-Jacobi Equations and Applications[END_REF].

The following result connects the solution of (5) to the sub and supersolutions of (6).

Lemma 1: For any T > 0, let φ be the viscosity solution of (5) with φ 0 Lipschitz in Ω. Then, for any α ≥ 1, φ1 (x, t) := φ(x, t) + c e -α t -sup

Ω×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| ∈ BU C(Ω × [0, T )) (7)
with c < 0 and φ2 (x, t) := φ(x, t) + d e -α t + sup + k sup

Ω×[0,T ) |w(x, t)|
Ω×[0,T ) |v(x, t)| ∈ BU C(Ω × [0, T )) (8)
with d > 0 are a viscosity subsolution and a viscosity supersolution of ( 6) with k ≥ α, respectively.

The proof of the previous lemma is given in Appendix. For the proposed control (shown in Fig. 1) the following ISS property holds. This is the main result of this paper (proven in appendix section).

Theorem 1: For the closed-loop system given by ( 2), (3), and (4) under the assumptions stated in the formulations of the problems ( 5) and ( 6)

, if k ≥ 1 the solutions φ, φ ∈ BU C(Ω × [0, T )) for all T > 0 satisfy the following bound |φ(•, t)| L2 ≤ √ 2 µ(Ω) 1/2 β + | φ(•, 0)| L2 e -kt + √ 2 µ(Ω) 1/2 sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| , (9)
for all t ∈ [0, T ), where β > 0 is such that |φ 0 (x)-φ0 (x)| ≤ β for all x ∈ Ω.

Remark. Some observations are in order.

• It is worth noting that the above result is more stringent as compared with the stability condition for the noisefree case presented in [5, Theorem 2, p. 4], for which it is required k > 0 instead of k ≥ 1. • As for nonlinear finite-dimensional control systems, there exist globally asymptotically stable NF equations that are unstable in presence of vanishing disturbance.

Next section will concern the presentation of the results obtained with a simple case study.

IV. A SIMPLE ONE-DIMENSIONAL EXAMPLE

In this section, we evaluate the effectiveness of the proposed approach by dealing with a simple one-dimensional example, which is analytically proved not to be ISS by showing an unstable behavior when subject to a converging disturbance in the absence of control. Next, we will illustrate how to stabilize this system by using the proposed control scheme and show simulations performed in Matlab with the level set methods toolbox developed by Mitchell [START_REF] Mitchell | The flexible, extensible and efficient toolbox of level set methods[END_REF].

Note that ISS property does not hold in general for globally asymptotically stable NF equations. Consider, for example, the one-dimensional NF equation

φ t (x, t) + f (x, t) |φ x (x, t)| = 0 in (0, 1) × [0, +∞) (10) where f (x, t) := λ(x) (t + 2) log(t + 2) λ(x) := min(x, 1 -x) (11)
with boundary conditions φ(0, t) = φ(1, t) = 0 for all t ≥ 0.

One can prove that ( 10) is globally stable in the L 2 sense by using a simple quadratic Lyapunov functional. In addition, it is straightforward to verify that [START_REF] Bernauer | Optimal control of the classical twophase Stefan problem in level set formulation[END_REF] with initial condition φ(x, 0) = λ(x) over [0, 1] admits the viscosity solution

φ(x, t) = log(2) λ(x) log(t + 2) , (12)
having an L 2 norm that converges asymptotically to zero. However, in the presence of an external in-domain disturbance as described by

φ t (x, t) + f (x, t) |φ x (x, t)| = w(x, t)
with initial condition φ(x, 0) = λ(x) and

w(x, t) := 2λ(x) log(2)(t + 2) (13)
tending to zero asymptotically, the viscosity solution

φ(x, t) = log(t + 2) λ(x) log(2)
is asymptotically divergent.

Let us address the problem to design the proposed control to [START_REF] Bernauer | Optimal control of the classical twophase Stefan problem in level set formulation[END_REF] and analyze it by means of simulations. The solution of the problems ( 5) and ( 6) together with (3) were obtained on the interval [0, 1] discretized via a regular grid composed of 100 discretization points. Both ( 5) and (6) were approximated with numerical schemes for nonlinear hyperbolic equations. In particular, a total variation diminishing Runge-Kutta scheme of second order was used for time discretization and an upwind second-order essentially non-oscillatory scheme [START_REF] Kimmel | Numerical Geometry of Images[END_REF] was utilized for space approximation. In all simulations we chose initial condition φ 0 (x) = φ0 (x) = λ(x) and boundary conditions ϕ(0) = ϕ(1) = 0 for all t ≥ 0.

In the first collection of tests, w was given by ( 13), whereas we set v ≡ 0. The graphs of φ(•, t) for different values of t ∈ [0, 1.5] are reported in Fig. 2 for k = 0, 0.4, 1, and 5. As expected, φ diverges for k = 0. Note that, even if for k = 0.4 the stability condition k ≥ 1 is not fulfilled, a stable behavior is observed since φ increases for a certain time and then decreases. This result is not surprising since the stability condition we found is proved to be only sufficient and not necessary, thus presumably conservative to some extent. For k = 1 and k = 5, φ shows an asymptotically stable behavior to the origin as expected. In the second test illustrated in Fig. 3, random noises uniformly distributed in the ranges [-0.08, 0.08] were considered for both w and v with the control parameter k set to 2.

V. CONCLUSION

The paper extends some of the results presented in [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF] by demonstrating that ISS bounds hold in a suitable sense for the considered observer-based control scheme for NF equations upon the choice of a sufficiently high gain. Such theoretical findings have been confirmed by simulation results. The relevance of this work concerns the potential application to the control of shape dynamics or moving fronts/interfaces, usually modelled by means of normal flow equations in two or higher dimensions, for which the present results hold. Thus, future efforts will be devoted to the control of more complicated systems such as airlift bioreactors (see, e.g., [START_REF] Bagnerini | Simulation of the rising of gas bubbles in a pilot-scale external loop airlift photobioreactor[END_REF], [START_REF] Neviani | Gas bubble dynamics in airlift photo-bioreactors for microalgae cultivation by level set methods[END_REF]) or in general multi-phase flows described by Navier-Stokes equations combined with NF equations to account for moving interfaces (e.g., bubble dynamics for bioreactors), while providing robustness guarantees in terms of ISS.

APPENDIX

In this appendix we first prove Lemma 1, and then Theorem 1.

A. Proof of Lemma 1

First, we show that φ1 is a viscosity subsolution of [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF].

Let θ ∈ C 1 (Ω × (0, T )). Since θ(x, t) := θ(x, t)-c e -α t + sup Ω×[0,T ) |w(x, t)|+k sup Ω×[0,T ) |v(x, t)| is a C 1 function in Ω × (0, T ) and φ is the viscosity solution of (5), it follows that θt (x 1 , t 1 ) + f (x 1 , t 1 )| θx (x 1 , t 1 )| -u(x 1 , t 1 , φ(x 1 , t 1 )) -w(x 1 , t 1 ) ≤ 0 if (x 1 , t 1 ) ∈ Ω × (0, T ) is a local maximum point of φ -θ.
By definition of θ, we obtain

θ t (x 1 , t 1 ) + c α e -α t1 + f (x 1 , t 1 ) |θ x (x 1 , t 1 )| -u(x 1 , t 1 , φ(x 1 , t 1
)) -w(x 1 , t 1 ) ≤ 0 and hence, using [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF],

θ t (x 1 , t 1 ) + f (x 1 , t 1 ) |θ x (x 1 , t 1 )| + α ( φ1 (x 1 , t 1 ) -φ(x 1 , t 1 )) -u(x 1 , t 1 , φ(x 1 , t 1 )) + α sup Ω×[0,T ) |w(x, t)| -w(x 1 , t 1 ) + k α sup Ω×[0,T ) |v(x, t)| ≤ 0.
If k ≥ α, as φ1 (x 1 , t 1 ) -φ(x 1 , t 1 ) < 0, we obtain

θ t (x 1 , t 1 ) + f (x 1 , t 1 ) |θ x (x 1 , t 1 )| + k ( φ1 (x 1 , t 1 ) -φ(x 1 , t 1 )) -k v(x 1 , t 1 ) -u(x 1 , t 1 , φ(x 1 , t 1 ) ≤ θ t (x 1 , t 1 ) + f (x 1 , t 1 ) |θ x (x 1 , t 1 )| + α ( φ1 (x 1 , t 1 ) -φ(x 1 , t 1 )) -u(x 1 , t 1 , φ(x 1 , t 1 )) -k v(x 1 , t 1 ) + α sup Ω×[0,T ) |w(x, t)| -w(x 1 , t 1 ) -α sup Ω×[0,T ) |w(x, t)| + w(x 1 , t 1 ) + k α sup Ω×[0,T ) |v(x, t)| -k α sup Ω×[0,T ) |v(x, t)| ≤ -α sup Ω×[0,T ) |w(x, t)| + w(x 1 , t 1 ) -k v(x 1 , t 1 ) + α sup Ω×[0,T ) |v(x, t)| ≤ 0
where the last inequality holds since α ≥ 1. Therefore, φ1 is a viscosity subsolution of ( 6) since (x 1 , t 1 ) is a local maximum also for

φ1 -θ = φ + c e -α t -sup Ω×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| -θ = φ -θ .
To prove that φ2 is a viscosity supersolution of (6), let (x 2 , t 2 ) ∈ Ω × (0, ) be a local minimum point of φ2 -θ, for all θ ∈ C 1 (Ω × (0, T )). Since θ(x, t) := θ(x, t)-d e -α t -sup

Ω×[0,T ) |w(x, t)|-k sup Ω×[0,T ) |v(x, t)|
is a C 1 function in Ω × (0, T ), φ is a viscosity supersolution of ( 5) and (x 2 , t 2 ) is also a local minimum point of φ -θ = φ2 -θ, we have

θt (x 2 , t 2 ) + f (x 2 , t 2 ) | θx (x 2 , t 2 )|- u(x 2 , t 2 , φ(x 2 , t 2 )) -w(x 2 , t 2 ) ≥ 0.
Thus, by definition of θ, it follows that

θ t (x 2 , t 2 ) + d α e -α t2 + f (x 2 , t 2 ) |θ x (x 2 , t 2 )| -u(x 2 , t 2 , φ(x 2 , t 2 )) -w(x 2 , t 2 ) ≥ 0
and, by [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], reads

θ t (x 2 , t 2 ) + f (x 2 , t 2 ) |θ x (x 2 , t 2 )| + α ( φ2 (x 2 , t 2 ) -φ(x 2 , t 2 )) -α sup Ω×[0,T ) |w(x, t)| -w(x 2 , t 2 ) -u(x 2 , t 2 , φ(x 2 , t 2 )) -k α sup Ω×[0,T ) |v(x, t)| ≥ 0. If k ≥ α, as φ2 (x 2 , t 2 ) -φ(x 2 , t 2 ) > 0, we obtain θ t (x 2 , t 2 ) + f (x 2 , t 2 ) |θ x (x 2 , t 2 )| + k( φ2 (x 2 , t 2 ) -φ(x 2 , t 2 )) -k v(x 2 , t 2 ) -u(x 2 , t 2 , φ(x 2 , t 2 )) ≥ θ t (x 2 , t 2 ) + f (x 2 , t 2 ) |θ x (x 2 , t 2 )| + α( φ2 (x 2 , t 2 ) -φ(x 2 , t 2 )) -u(x 2 , t 2 , φ(x 2 , t 2 )) -α sup Ω×[0,T ) |w(x, t)| -w(x 2 , t 2 ) + α sup Ω×[0,T ) |w(x, t)| + w(x 2 , t 2 ) -k v(x 2 , t 2 ) + k α sup Ω×[0,T ) |v(x, t)| -k α sup Ω×[0,T ) |v(x, t)| ≥ α sup Ω×[0,T ) |w(x, t)| + w(x 2 , t 2 ) + k α sup Ω×[0,T ) |v(x, t)| -v(x 2 , t 2 ) ≥ 0 ,
namely, φ2 is a viscosity supersolution of ( 6). This concludes the proof of Lemma 1.

B. Proof of Theorem 1

Before proving the main result, let us first state and prove the following intermediate result based on Lemma 1.

Lemma 2: For any α ≥ 1, there exists β > 0 such that

| φ(x, t)| ≤ β e -α t + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)|, t ∈ [0, T ) (14)
for all T > 0, where φ(x, t) := φ(x, t) -φ(x, t) and φ, φ are viscosity solutions of ( 5) and ( 6), respectively, with k ≥ α and φ 0 , φ0 Lipschitz in Ω.

Proof. Since φ 0 , φ0 ∈ C 0 (Ω), there exists β > 0 such that |φ 0 (x) -φ0 (x)| ≤ β. If k ≥ α, by Lemma 1 we can write that φ1 (x, t) = φ(x, t) -β e -α t -sup

Ω×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| and φ2 (x, t) = φ(x, t) + β e -α t + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)|
are sub and supersolution of (4), respectively. Moreover, Using the feedback control (3), i.e., u(x, t) = -k φ(x, t) -k v(x, t) + f (x, t) φx (x, t) , in (6) we get φt (x, t) = -k φ(x, t) and hence it is straightforward to prove the L 2 exponential stability of φ to zero (see [START_REF] Alessandri | State and observer-based feedback control of normal flow equations[END_REF]Theorem 2,p. 4]). Using such a property, the Young inequality φ(x, t) 2 = φ(x, t) -φ(x, t) + φ(x, t) for any k ≥ α with α ≥ 1. Thus, the choice α = k ≥ 1 provides the tightest bound given by ( 9). This concludes the proof of Theorem 1.

φ1 (x, 0) = φ 0 (x) -β -sup Ω×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| ≤ φ0 (x) ≤ φ 0 (x) + β + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| = φ2 (x,

Fig. 1 .

 1 Fig. 1. Sketch of the control loop.

Fig. 2 .Fig.

 2 Fig. 2. Simulation results with in-domain disturbances (13) and control (3) for different values of k = 0, 0.4, 1, 5. The initial condition is depicted in red. For k = 0 the solution diverges in accordance with what analytically shown, whereas a small overshoot occurs for k = 0.4 since φ increases for a certain time and then decreases.

  0)for x ∈ Ω and φ1 (x, t) = ϕ(x) -βe -αt -supΩ×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| ≤ φ(x, t) = ϕ(x) ≤ φ2 (x, t) = ϕ(x) + βe -αt + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| on ∂Ω × (0, T ). Hence, using [17, Theorem 5.1], it follows that, on Ω × [0, T ), φ1 (x, t) = φ(x, t) -β e -α t -sup Ω×[0,T ) |w(x, t)| -k sup Ω×[0,T ) |v(x, t)| ≤ φ(x, t) ≤ φ(x, t) + β e -α t + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| = φ2 (x, t) and | φ(x, t) -φ(x, t)| ≤ β e -α t + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)|, t ∈ [0, T )for all T > 0, thus obtaining[START_REF] Koga | Towards implementation of PDE control for Stefan system: Input-to-state stability and sampled-data design[END_REF] and concluding the proof of Lemma 2.Using Lemma 2, we are now in position to prove the main result.Proof of Theorem 1. Since φ, φ ∈ BU C(Ω × [0, T )), from (14) it follows that | φ(•, t)| 2 L2 = Ω | φ(x, t)| 2 dx ≤ β e -α t •, t)| L2 ≤ µ(Ω)1/2 β e -α t + sup Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| .

2

  ≤2 φ(x, t) 2 + 2 φ(x, t) 2 ,

										+	√	2 µ(Ω) 1/2	sup	|w(x, t)| + k sup	|v(x, t)|
										Ω×[0,T )	Ω×[0,T )
	and the inequality	√	a + b ≤	√	a+	√	b for a, b ≥ 0, it follows
	that
	|φ(•, t)| L2 ≤	√	2		φ(x, t) 2 dx	1/2 +	√	2	φ(x, t) 2 dx	1/2
	≤	√	Ω 2 µ(Ω) 1/2 β e -α t + sup	Ω

Ω×[0,T ) |w(x, t)| + k sup Ω×[0,T ) |v(x, t)| + √ 2 Ω φ(x, 0) 2 e -2kt dx 1/2 = √ 2 µ(Ω) 1/2 β e -α t + | φ(•, 0)| L2 e -kt


Lei Zhu

Fernando Soldevila

Claudio Moretti

Alexandra D'arco

Antoine Boniface

Xiaopeng Shao

Hilton B De Aguiar

Sylvain Gigan
email: sylvain.gigan@lkb.ens.fr

Large field-of-view non-invasive imaging through scattering layers using fluctuating random illumination

Non-invasive optical imaging techniques are essential diagnostic tools in many fields. Although various recent methods have been proposed to utilize and control light in multiple scattering media, non-invasive optical imaging through and inside scattering layers across a large field of view remains elusive due to the physical limits set by the optical memory effect, especially without wavefront shaping techniques. Here, we demonstrate an approach that enables non-invasive fluorescence imaging behind scattering layers with field-of-views extending well beyond the optical memory effect. The method consists in demixing the speckle patterns emitted by a fluorescent object under variable unknown random illumination, using matrix factorization and a novel fingerprintbased reconstruction. Experimental validation shows the efficiency and robustness of the method with various fluorescent samples, covering a field of view up to three times the optical memory effect range. Our non-invasive imaging technique is simple, neither requires a spatial light modulator nor a guide star, and can be generalized to a wide range of incoherent contrast mechanisms and illumination schemes.

INTRODUCTION

Non-invasive optical imaging has important applications in various fields ranging from biotechnology [START_REF] Zhao | Non-invasive detection of apoptosis using magnetic resonance imaging and a targeted contrast agent[END_REF], [START_REF] Artzi | In vivo and in vitro tracking of erosion in biodegradable materials using non-invasive fluorescence imaging[END_REF] to optical detection [START_REF] Kozloff | Non-invasive optical detection of cathepsin Kmediated fluorescence reveals osteoclast activity in vitro and in vivo[END_REF]. However, inhomogeneous samples, such as biological tissues, scatter light, which results in a complex speckle pattern on the detector [START_REF] Goodman | Some fundamental properties of speckle[END_REF], [START_REF] Bender | Customizing speckle intensity statistics[END_REF]. With increasing depth, separating the low amount of ballistic light from the scattered light becomes a big challenge [START_REF] Abramson | Light-in-flight recording by holography[END_REF], [START_REF] Huang | Optical coherence tomography[END_REF]. Over the years, many approaches have been put forward to overcome this problem by exploiting or suppressing the scattered light. With the development of spatial light modulators (SLMs), multiple ways to control and manipulate scattered light have been developed [START_REF] Mosk | Controlling waves in space and time for imaging and focusing in complex media[END_REF], [START_REF] Rotter | Light fields in complex media: Mesoscopic scattering meets wave control[END_REF]. Several techniques have been proposed to focus light by making use of feedback signals to optimize the incident wavefront to recreate a focus that is then used for raster-scanning microscopy [START_REF] Vellekoop | Focusing coherent light through opaque strongly scattering media[END_REF], [START_REF] Horstmeyer | Guidestar-assisted wavefrontshaping methods for focusing light into biological tissue[END_REF]. These techniques require access to both sides of the scattering layer to optimize the wavefront, which strongly limits their application in realcase scenarios. To overcome this, other strategies have been proposed based on wavefront shaping and various feedback signals such as fluorescence or ultrasound signals [START_REF] Horstmeyer | Guidestar-assisted wavefrontshaping methods for focusing light into biological tissue[END_REF], [START_REF] Katz | Controlling light in complex media beyond the acoustic diffraction-limit using the acoustooptic transmission matrix[END_REF], [START_REF] Popoff | Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media[END_REF], [START_REF] Hofer | Manipulating the transmission matrix of scattering media for nonlinear imaging beyond the memory effect[END_REF]. However, these approaches either require long acquisition times or are limited to small fields of view (FoV). On the other hand, several techniques exploiting the angular speckle correlations, known as the optical memory effect (ME) [START_REF] Freund | Memory effects in propagation of optical waves through disordered media[END_REF], [START_REF] Yılmaz | Angular memory effect of transmission eigenchannels[END_REF], [START_REF] Osnabrugge | Generalized optical memory effect[END_REF], have also been proposed for imaging objects hidden behind scattering media [START_REF] Bertolotti | Non-invasive imaging through opaque scattering layers[END_REF], [START_REF] Katz | Non-invasive singleshot imaging through scattering layers and around corners via speckle correlations[END_REF]. While these approaches are fast, their FoV is still limited by the ME range.

Linear fluorescence is widely used in biology and biomedical sciences [START_REF] Ruan | Fluorescence imaging through dynamic scattering media with speckle-encoded ultrasoundmodulated light correlation[END_REF], [START_REF] Lichtman | Fluorescence microscopy[END_REF], [START_REF] Mangeat | Super-resolved live-cell imaging using random illumination microscopy[END_REF]. It enables imaging of cellular, subcellular, or molecular components, and has the advantages of high spatial resolution, contrast, and speed. Recent advances have allowed both focusing and imaging through scattering media using fluorescent light. Even so, these methods either rely on the use of guide stars [START_REF] Horstmeyer | Guidestar-assisted wavefrontshaping methods for focusing light into biological tissue[END_REF], are limited to the ME range [START_REF] Hofer | Wide field fluorescence epi-microscopy behind a scattering medium enabled by speckle correlations[END_REF], or need to characterize the scattering medium [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF].

Here, we present a robust approach that allows to noninvasively image through scattering layers at depths far beyond the ME range, using fluctuating random speckle illumination simply generated with a rotating diffuser. Once excited, each fluorescent emitter generates a unique speckle pattern on the detector, which constitutes its fingerprint. Each image captured by the camera is an incoherent sum of the fingerprints from the emitters, with different relative weights due to the variable random illumination patterns and sample structure. To retrieve each individual fingerprint, we capture a set of images while randomly changing the illumination, and use a Non-negative Matrix Factorization (NMF) algorithm to demix the set of acquired frames. After that, the fingerprints are used to reconstruct the final image by exploring the correlations between them. To validate the technique, we experimentally demonstrate our non-invasive approach both on fluorescent beads and on continuous fluorescent objects.

RESULTS

The experimental setup is depicted in Fig. 1a. A rotating holographic diffuser modulates the incident light coming from a laser by adding a random phase when the light propagates through it. Then, the modulated light travels through the scattering medium and generates a random unknown speckle arXiv:2107.08158v1 [physics.optics] 17 Jul 2021 Fluo. Obj.

Scat.

Ground truth Reconstruction

Fluo. Obj.

Scat.

I speckle. which illuminates the object. The excitation speckle induces a fluorescent response from the object, which propagates back through the scattering medium, generating an incoherent sum of speckles that is measured by the camera. Although the captured images are low-contrast, random, and seemingly information-less, they contain all the fingerprints from the independent emitters of the object, but with time-varying weights. Furthermore, independent emitters within the ME range will produce correlated but shifted fingerprints on the camera [START_REF] Freund | Memory effects in propagation of optical waves through disordered media[END_REF], while emitters outside the ME range will produce totally uncorrelated fingerprints. For a given speckle illumination, the captured image, I fluo , can be expressed as a linear superposition of those fingerprints with different weights. Thus, the camera image is given by:

CW

I f luo (r, t) = P k=1 w k (r)h k (t), (1)
where I f luo (r, t) corresponds to a low contrast speckle for the t th illumination, r is the spatial coordinate, w k (r) represents the fingerprint of the k th independent emitter of the object, h k (t) stands for the amount of excitation light at the k th emitter during the t th illumination, and P is the number of independent emitters. Given enough different random illuminations, a collection of frames can be used to retrieve each individual fingerprint, w k (r), by using an NMF algorithm, that we will now explain in detail.

Fingerprint demixing procedure. After randomly exciting the object with a variety of t random speckles, a series of camera images, I f luo (r, t), are collected. It is possible to retrieve the fingerprints, w k , corresponding to each independent emitter, from the measurements I f luo (r, t) by finding the solution to the minimization problem [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]:

min W >0,H>0 I -W H 2 , (2)
This minimization problem can be formulated as a low rank factorization, where the matrix I ∈ R r×t + contains all the I f luo (r, t), can be approximated with two real positive matrices W ∈ R r×ρ + (the fingerprints) and H ∈ R ρ×t + (the temporal evolutions), where r are the pixels, ρ is the estimated rank of I and t indicates the frames. Since the collected images and the demixed fingerprints are positive, this problem corresponds exactly to the family of NMF problems. The NMF framework has been employed in demixing scenarios, both in structural imaging [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF] and functional imaging [START_REF] Moretti | Readout of fluorescence functional signals through highly scattering tissue[END_REF], [START_REF] Pégard | Compressive light-field microscopy for 3D neural activity recording[END_REF], [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF]. In our case, the estimated rank ρ approximately corresponds to the number of independent emitters P and it can be estimated from the data by minimizing the root mean squared residual of NMF as a function of the rank (see Supplementary I).

Fingerprint-based reconstruction. After the demixing step, the fingerprints are retrieved. Due to the ME, emitters close to each other will produce highly correlated fingerprints, with a spatial shift that is directly determined by their relative position [START_REF] Katz | Non-invasive singleshot imaging through scattering layers and around corners via speckle correlations[END_REF]. By exploring these correlations, a position map of the emitters can be recovered, thus yielding an image of the object. Several approaches can be used to explore the correlations and calculate the shifts between fingerprints. Usually, this process is performed by doing a cross-correlation between the fingerprints and locating the position of the maximum [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF]. However, here we introduce a novel approach based on deconvolution, that we denote as Fingerprint-based Reconstruction (FBR). Compared to the cross-correlation procedure, we found that this approach allows to suppress noise and strongly improves the quality of the reconstruction.

No matter if two emitters are within the same ME patch or not, one can perform the pairwise deconvolution of the i-th emitter by the k-th emitter, which can be written as:

arg min o i,k µ 2 ||w i -o i,k w k || 2 2 + ||o i,k || T V , (3)
where µ is a regularization parameter, denotes the convolution operator, ||•|| 2 indicates the L 2 norm and ||•|| T V represents the Total Variation (TV) norm, and the two fingerprints are w i (considered as the "image") and w k (considered as the "point spread function" (PSF)). When the two emitters lay within one ME range, the pairwise deconvolution yields a uniform image with a narrow delta-like peak, which is located at a distance from the center given by the relative position of the two emitters ( r i,k = r i -r k ). If the two emitters are located beyond the ME range, the deconvolution yields noise.

For a given emitter k, it is possible to obtain, O k , the partial image of the object in the vicinity of the emitter, by simply adding the result of all the pairwise deconvolutions related to that emitter, o i,k . (see Fig. 1b-d).

O k = ρ i=1 o i,k , (4)
Even if the ensemble of emitters expands well beyond the ME range, the full spatial distribution can be recovered if the different isoplanatic patches are "connected" by emitters (see Fig. 1d). For example, if emitters i and k are beyond the ME range but emitter j is between them, we can always calculate the shift between them as r i,k = r i,j + r j,k . The global reconstruction O Global can be obtained by composing all the partial images, O k , into one image, taking into account their relative positions with respect to the first emitter, r k,1 :

O Global = ρ k=1 O k ( r -r k,1 ), (5)
Experimentally, we prove that our technique can be used to recover very sparse objects by using 2D distributions of beads with a diameter of 1 µm. As shown in Fig. 2, our approach can reconstruct the objects without constraints from the ME range. Non-sparse, continuous objects are common in scattering biological samples, and this often poses a difficult challenge for non-invasive imaging through scattering media approaches. To demonstrate that our technique can also work with non-sparse and continuous objects, we use fluorescencestained pollen grains, whose reconstructed images are shown in Fig. 3. The size of fluorescence-stained pollen grains in Fig. 3 is smaller than the ME range and the process of reconstruction used for them is the same as what is presented in Fig. 1.

In both Fig. 2 and Fig. 3 we report the reconstruction of hidden beads distributions and pollen seeds. However, the more complex the objects are, the larger the required number of independent illuminations. Indeed, given more speckle patterns, NMF is able to provide more reliable fingerprints, thus a more reliable reconstruction (See supplementary Fig. 10). We note that a few hundred illuminations are sufficient to recover the object reasonably well in our case. Importantly, our technique is not limited by the number of independent illuminations that we can generate with the rotating diffuser, as it is possible to produce a very large number of independent illumination with different rotating diffusers by tunning their scattering angle

DISCUSSION

Here, we report on successfully recovering the hidden objects by exploiting the correlation between the fingerprints. However, we believe, based on the spectral ME or the 3D ME, that the technique could be used to recover multi spectral or 3D objects. This is a significant advantage as it does not require the scattering medium to be static. Another important future direction will be to explore the ME of dynamic scattering media to recover the hidden object through it. Interestingly, our technique, contrarily to the techniques based on the optical transmission matrix, does not require focusing light on the object, and relies only upon the random video frames generated from the random illumination. Another important point is the simplicity of the technique, which does not require calibration of the PSFs [START_REF] Antipa | DiffuserCam: lensless single-exposure 3D imaging[END_REF] of the imaging system, and can be implemented without the need of expensive SLMs. Thus, our technique can be easily implemented in various scattering media and imaging scenarios. It is necessary to stress that our approach reconstructs both non-sparse and continuous objects over previous autocorrelation approaches [START_REF] Bertolotti | Non-invasive imaging through opaque scattering layers[END_REF], [START_REF] Katz | Non-invasive singleshot imaging through scattering layers and around corners via speckle correlations[END_REF].

In conclusion, we have shown a non-invasive technique to computationally retrieve images of objects hidden behind a scattering medium from low-contrast fluorescent speckles using random illumination. We have demonstrated that our approach works with both sparse objects, even beyond the ME range, and continuous objects. Importantly, the proposed approach neither relies upon ballistic light nor uses wavefront shaping, and it is adaptable to various scattering media and objects. Our technique is flexible, robust, and opens a promising avenue towards deep fluorescence imaging in highly scattering media. Finally, it can be generalized to a wide range of incoherent contrast mechanisms and illumination schemes.

METHODS

Experimental setup. A continuous-wave laser (λ = 532 nm, Coherent Sapphire) is expanded and illuminates the rotating holographic diffuser (Edmund, DG10). Then the modulated light is delivered onto the fluorescent sample through a 200 mm lens (LA1708-A, Thorlabs) and objective (Zeiss W "Plan-Apochromat" ×20, NA 1.0). After excitation, the fluorescence is scattered by the medium and collected with a 150 mm tube lens (L, AC254-150-A, Thorlabs), which is employed to produce an image onto the detector, a sCMOS camera (Hamamatsu ORCA Flash). Two dichroic filters (short pass 532 nm, Thorlabs and 533 nm notch MF525-39, Thorlabs) are used to block any signal that does not come from the fluorescence emission. The fluorescent objects, which are made of orange beads (540/560 nm, Invitrogen FluoSpheres, size 1.0 µm) or pollen seeds (Carolina, Mixed Pollen Grains Slide, w.m.), are placed below the scattering medium. The distance between the scattering medium and the fluorescent objects is 0.2 mm. A transmission pathway that consists of a microscope objective (Olympus "MPlan N" ×20, NA 0.4), a 150 mm tube lens (L, AC254 -150-A, Thorlabs), and CCD camera (Allied Vision, Manta), is used as a passive control only. This control part is used to correctly select the position of fluorescent object with a white light source (Moritex, MHAB 150W) and it also allows us to align the experimental setup. For the scattering medium, we either use a single holographic diffuser (Newport, 10DKIT-C1,10°) or a combination of two (Newport, 10DKIT-C1,10°and 10DKIT-C1,1°) in order to get different memory effect ranges.

The exposure time has been set from 10 ms to 20 ms, depending on the scattering medium and the fluorescent sample. Once captured, the speckle images, which contain few tens of speckle grains, are cropped from the raw images. Then, a high-pass Fourier Gaussian filter is employed to remove the background from the cropped images and the processed data set is analyzed with the NMF algorithm to obtain the fingerprint of each emitter. The experimental setup is shown in Supplementary V.

In Fig. 2a-b, the size of each cropped image is 70 × 72 pixels, the number of patterns t is 5120, the scattering medium is a holographic diffuser (Newport, 10DKIT-C1, 10), and the exposure time is 15 ms. In Fig. 2c-d, the size of each cropped image is 64 × 66 pixels, the number of patterns t is 5120, the scattering medium are two holographic diffusers (Newport, 10DKIT-C1, 10 + 1), and the exposure time is 20 ms. In Fig. 3, the size of each cropped image is 70×74 pixels, the number of patterns t is 5120, the holographic diffuser (Newport, 10DKIT-C1, 10) is used as the scattering medium, and the exposure time is 10 ms.

NMF+fingerprint-based reconstruction algorithm. For the NMF, knowing the rank of the system is necessary. The rank ρ is estimated by looking at the root mean square residual ||I f luo -W H|| F ro as a function of rank ρ (detailed in Supplementary I), where ||•|| F ro indicates the Frobenius norm, and minimizing it. For the NMF, a random initialization is employed. The retrieved fingerprints are used as the input data of the FBR.

I. ESTIMATING THE RANK OF THE SYSTEM

In order to operate the NMF algorithm, the rank, ρ, is a unique parameter which needs to be pre-determined. Several methods can be used to estimate this rank, such as looking at clusters using a k -means algorithm [START_REF] Moretti | Readout of fluorescence functional signals through highly scattering tissue[END_REF], or by estimating the 'elbow' of the loss function [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF], [START_REF] Hutchins | Positiondependent motif characterization using non-negative matrix factorization[END_REF]. In our case, we use a simple procedure based on the study of the evolution of the root mean square residual (RMSR) of the NMF for different ranks to estimate the true rank of the system. Once we acquire a dataset, it is possible to demix it by setting any desired rank ρ. For each ρ, we can estimate the quality of the NMF solution by the RMSR: ||I f luo -W H|| F ro . As we change the rank, the RMSR varies, reaching a minimum when the estimated rank equals the true rank of the system. Also, when overestimating the rank, it is possible to have a big variance on the result of the NMF for multiple realizations. Both these behaviors can be used to delimit the true value of the rank of the system. As an example of the procedure, we show both a simulation with different numbers of beads and the same procedure with experimental datasets in Fig. 4. As can be seen in Fig. 4a, the minimum value of RMSR is ρ = 5, 10, and 15, respectively, which is in excellent agreement with the ground truth number of sources (5, 10, and 15). The rank estimation procedure is performed 12 times with different random initializations in order to see the variance of the results. Furthermore, we investigate the effect of incorrectly estimating the ranks ρ of the system on the quality of our reconstruction. We experimentally choose a fluorescent object which contains 10 beads. Then, we recover the object using our approach with different rank estimations, ranging from 6 to 15. The structural similarity index metric (SSIM) [START_REF] Wang | Image Quality Assessment: From Error Visibility to Structural Similarity[END_REF] is used to quantify the performance of reconstruction. As shown in Fig. 5, our technique is capable of achieving reconstruction with different ranks even if the proposed rank is not fully accurate. However, underestimating the rank may lead to losing some beads in the global reconstruction. On the other hand, overestimating the rank on the NMF will produce mixed patterns that hinder the capability to obtain the relative position between different emitters.

II. OPTICAL MEMORY EFFECT

With the purpose of proving that our technique can be applied in imaging beyond the ME, different diffusers with different scattering properties are tested. We estimate the ME range of each diffuser by exploring the correlation between the speckle patterns generated while laterally moving an emitter (a small fluorescent bead). As the source moves, the correlation between speckle patterns gradually decreases, reaching a minimum correlation when the distance is longer than the ME range. Measuring this distance allows estimating the ME range for each diffuser. As reported in Fig. 6, the ME range is approximately 50 µm and 20 µm corresponding to diffuser #1 and diffuser #2, respectively. Our non-invasive technique not only can retrieve objects within ME range as shown in Fig. 2a of the manuscript, which is possible to be recovered via autocorrelation approaches, but also can reconstruct an extended object, which could not be retrieved with autocorrelation approaches. As shown in Fig. 2b of the manuscript, the size of the reconstructed extended object is around 58 µm × 22 µm with the diffuser #2.

III. FINGERPRINT-BASED RECONSTRUCTION

As described in the main text, it is possible to reconstruct the whole object as long as the ME patches have some overlaps. To present the details of FBR, we experimentally choose a fluorescent object which contains 3 fluorescent beads and display the various results from the pairwise deconvolution, o i,k , and the different partial images, O k . As shown in Fig. 7, the O k of the emitter k can be recovered by choosing w k as the PSF. By looking at the maximum value of those o i,k , the shift r i,k between fingerprints w i and w k can be retrieved, as shown in Fig. 1c of the manuscript. In practice, the fingerprints coming from two emitters which are beyond the ME range will not provide the relative position information of their emitters (as they will be totally uncorrelated). Thus, it is necessary to infer whether the fingerprints w i and w k are within or beyond the ME range. In our method, we study α = max{o i,k } max{o k,k } as a function of relative distance, where max{o i,k } stands for maximum value of o i,k . A given threshold, α tres of α is introduced to evaluate it. For example, if α is greater than α tres , w i and w k belong to the same ME range. Otherwise, they belong to different ME ranges, and their relative position cannot be retrieved. As shown in Fig. 8, we experimentally investigate the maximum value of various results of pairwise deconvolutions as a function of distance between emitters. We set α tres to 0.01. The full workflow of technique is depicted in Alg.1.

Algorithm 1: Image retrieval procedure

Input: Series of camera images, I f luo (r, t).

Output: Image of the object, O Global . Estimate the rank (ρ) of the system from I f luo (r, t) (see Supplementary I). Retrieve the spatial fingerprints (w i ) by using NMF with the estimated rank.

IV. PRODUCING RANDOM ILLUMINATION USING A SLM

Our technique retains the potential of employing a SLM for producing random illumination. The SLM can produce a large number of independent speckle illuminations. Experimentally, we replace the rotating diffuser with the SLM and perform experiments on fluorescent point-like objects and continuous volumetric objects by generating random illuminations with the SLM. The reconstruction is presented in Fig. 9. To prove that it is possible to recover a more reliable image with more patterns, we show the reconstruction with the different number of patterns in Fig. 10. Note that the results as shown in Fig. 9 and Fig. 10 are just tests to recover with the illumination patterns produced by SLM that allows generating a higher diversity of patterns than our current rotating diffuser, but it is not a limitation of our technique, as we could use different ro-tating diffusers to get more independent illumination patterns.

Figure 1 .

 1 Figure 1. Schematic of the experimental setup and reconstruction principle. (a) Schematic view of experimental setup. A coherent light source illuminates a rotating diffuser in order to excite the fluorescent object through a scattering medium with a random modulated speckle pattern. Once excited, the emitted signal from the fluorescent objects is recorded with a camera. I f luo is a series of t fluorescent speckles corresponding to different random speckle illuminations. The fingerprints can be recovered from I f luo by using NMF. (b-e) Fingerprint-based reconstruction. (b) Pairwise deconvolution (labelled as -1 ) between all the possible pairs of emitter fingerprints is performed. (c) The result of each deconvolution provides the relative position between one emitter and its neighbors. (d) By adding the resulting images for each emitter, it is possible to recover a partial image of the object centered at that emitter (see Eq.4). (e) All the partial images can be merged into the final reconstruction according to the relative position between neighbouring emitters. Dashed circle indicates the optical memory range. Scale bar sizes are 10 µm. RD: rotating diffuser, DM: dichroic mirror, OB: objective, Scat.: scattering medium, Fluo. Obj.: fluorescent object, SF: spectral filter, TL: tube lens.

Figure 2 .

 2 Figure 2. Experimental results of imaging through a scattering medium with fluorescent beads. (a)-(b) Fluorescent images of beads recorded without scattering medium. (c)-(d) Reconstruction of the object using NMF+FBR approach. The estimated rank of NMF is ρ = 26 for (c) and ρ = 16 for (d). In both cases, t = 5120 fluorescent speckle patterns are captured. The exposure time of (c) and (d) is set to 15 ms and 20 ms, respectively. Dashed circles indicate the optical memory effect range.

Figure 3 .

 3 Figure 3. Experimental results of imaging through scattering media with continuous objects. (a)-(b) Fluorescent images of different pollen seed structures recorded without scattering medium. (c)-(d) Reconstruction of the object with the NMF+FBR approach. The estimated rank for the NMF is ρ = 68 for (c) and ρ = 85 for (d), respectively. In both cases, t = 5120 fluorescent speckle patterns are recorded with an exposure time of 10 ms. Dashed circle indicates the optical memory effect range.

Figure 4 .

 4 Figure 4. NMF rank estimation procedure. (a) Simulation results with P fluorescent sources. The triangle symbol marks the minimum mean value of the RMSR, the solid line stands for the mean value of the RMSR over 12 times with the random initialization value of NMF, and the asterisk is the error bar that indicates its standard deviation in both positive and negative direction. The minimum mean RMSR value is achieved when the estimated rank ρ equals the true rank P of the system. (b) Experimental results with a fluorescent object which contains P = 11 beads. (c) Estimating the rank of a fluorescent object shown in (b), the solid line stands for the mean value of the RMSR over 12 repetitions, the triangle symbol marks the minimum mean value of the RMSR, and the error bars indicate the standard deviation of the RMSR.

Figure 5 .

 5 Figure 5. The plot of SSIM coefficient between different reconstructions and reference image. The reconstruction corresponding to ρ = 10 is treated as a reference image (the number of beads in the sample is 10).

Figure 6 .

 6 Figure 6. Plot of maximum value of cross-correlations between different patterns corresponding to different fluorescent bead positions for diffuser #1 (a) and #2 (b), where the colorcode indicates the different data sets.

  for k = 1, ..., ρ do Perform the pairwise deconvolution between w k and all the other fingerprints (w i =k ). Retrieve the relative position between emitter k and its neighbours inside the ME range. Calculate the partial image of the object in the vicinity of the emitter (O k ) by adding the result of all the pairwise deconvolutions related to that emitter (o i,k ). end Merge all the partial images (O k ) into the final reconstruction (O Global ) using the relative position between emitters.

Figure 7 .Figure 8 .

 78 Figure 7. Detail of FBR. (a) Ground truth taken without diffuser. (b) The detail of pairwise deconvolution. The estimated rank ρ of this data set is 3 and the exposure time is 10 ms.

Figure 9 .

 9 Figure 9. Experimental results of imaging through scattering media by producing random illumination with a SLM. (a,c) are ground truths of fluorescent beads object and fluorescent continuous volumetric object. (c,d) are reconstructions corresponding to (a,b), respectively. The estimated rank ρ of (c) and (d) is 14 and 53, respectively. t = 10240 for (c) and t = 5120 for (d) fluorescent speckle patterns are recorded. The dataset (c) is recorded with an exposure of 50 ms and the exposure time of (d) is 20 ms. Dashed circle indicates the optical memory effect range.

Figure 10 .Figure 11 .

 1011 Figure 10. Performance of reconstruction with different number of speckle patterns by producing random illumination with a SLM. (a) Ground truth of fluorescent continuous volumetric object taken without scattering medium. (c)-(d) Reconstructions with different number of illumination patterns. The rank of the NMF is estimated at ρ = 27 for (b), at ρ = 31 for (c), at ρ = 33 for (d), at ρ = 36 for (e), at ρ = 41 for (f), and at ρ = 53 for (g). The dataset is recorded with an exposure of 20 ms.

Acknowledgements. This research has been funded by the FET-Open (Dynamic-863203), European Research Council ERC Consolidator Grant (SMARTIES-724473), China Scholarship Council (CSC) (201906960055).