Finaritra Raoelijaona

Stéphane Thore
email: stephane.thore@inserm.fr

Sébastien Fribourg

Domain definition and interaction mapping for the endonuclease complex hNob1/hPno1

Keywords: ribosome biogenesis, endonuclease, pre-rRNA, KH domain, PIN domain, zinc ribbon

Ribosome biogenesis requires a variety of trans-acting factors in order to produce functional ribosomal subunits. In human cells, the complex formed by the proteins hNob1 and hPno1 is crucial to the site 3 cleavage occurring at the 3'-end of 18S pre-rRNA. However, the properties and activity of this complex are still poorly understood. We present here a detailed characterization of hNob1 organization and its interaction with hPno1. We redefine the boundaries of the endonuclease PIN domain present in hNob1 and we further delineate the precise interacting modules required for complex formation in hNob1 and hPno1. Altogether, our data contributes to a better understanding of the complex biology required during the site 3 cleavage step in ribosome biogenesis.

Introduction

Ribosome biogenesis is one of the most complex pathway existing in eukaryotic cells. This is due to the high number of factors involved and to the variety of the enzymatic reactions occurring along the entire pathway [START_REF] Barandun | Assembly and structure of the SSU processome-a nucleolar precursor of the small ribosomal subunit[END_REF]. In addition, the pathway is physically segregated between the nucleus and the cytoplasm by the nuclear membrane, which acts as a natural checkpoint. In recent years, a number of studies have revealed with atomic details the architecture of pre-ribosomal particles such as the 90S, the pre-60S and the pre-40S ribosomes prior or after their export to the cytoplasm [START_REF] Scaiola | Structure of a eukaryotic cytoplasmic pre-40S ribosomal subunit[END_REF][3][4] [START_REF] Heuer | Cryo-EM structure of a late pre-40S ribosomal subunit from Saccharomyces cerevisiae[END_REF][6] [START_REF] Sun | Molecular architecture of the 90S small subunit preribosome[END_REF][8] [START_REF] Sanghai | Modular assembly of the nucleolar pre-60S ribosomal subunit[END_REF]. It is noteworthy that most studies have been performed using fungal preribosomal particles, and therefore may differ from the situation encountered in human cells. The study of human ribosome biogenesis is directly relevant for understanding the set of mammalian diseases affecting this pathway and collectively referred to as ribosomopathies [START_REF] Narla | Ribosomopathies: human disorders of ribosome dysfunction[END_REF]. In this context, dedicated studies of the structure and function of the mammalian factors are highly relevant.

The recent structural studies on fungal ribosome biogenesis revealed the location of numerous trans-acting factors on the pre-ribosomal particles allowing for a comprehensive understanding of their biological role [START_REF] Peña | Eukaryotic ribosome assembly, transport and quality control[END_REF]. However, some critical factors are only partially visible or even absent from these snapshot images of the ribosome assembly process. For example, comprehensive structural data is lacking for Rrp5, a protein involved in both the pre-40S and the pre-60S synthesis pathway and for human Nob1 (hereafter hNob1), the endonuclease required for the cleavage at site 3 during the last step of pre-40S maturation. hNob1 is often associated with its natural protein partner hPno1/hDim2, which in turn regulates Nob1's catalytic activity [START_REF] Tone | Nob1p is required for biogenesis of the 26S proteasome and degraded upon its maturation in Saccharomyces cerevisiae[END_REF]. The position of the yeast Pno1 protein on late pre-ribosome particles has been recently identified [START_REF] Heuer | Cryo-EM structure of a late pre-40S ribosomal subunit from Saccharomyces cerevisiae[END_REF], but Nob1's location on the ribosomal particle in any organism was until very recently unknown [START_REF] Ameismeier | Visualizing late states of human 40S ribosomal subunit maturation[END_REF]. In addition, it appears that functional differences exist for Nob1 depending on the organism and thus studies on orthologues are of specific interest. As an example, the protein Fap7 has been shown to inhibit Nob1's endonuclease activity in Pyrococcus horikoshii [START_REF] Hellmich | Essential ribosome assembly factor Fap7 regulates a hierarchy of RNA-protein interactions during small ribosomal subunit biogenesis[END_REF], while the Fap7 protein promotes the same endonucleolytic activity in human [START_REF] Bai | The ATPase hCINAP regulates 18S rRNA processing and is essential for embryogenesis and tumour growth[END_REF]. The hNob1 is therefore a paradigm of an essential trans-acting protein for which limited molecular information is available.

The structure of the Nob1 protein from P. horikoshii is available, but does not allow for full modelling of the mammalian protein due to various large sequence insertions, which do not exist in the archaeal protein [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF]. hNob1 harbors a similar organization to the archaeal orthologue and consists of PIN (PilT N terminus) endonuclease and zinc ribbon domains, separated by linker region that is poorly conserved with the exception of its last C-terminal residues [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF], [START_REF] Fatica | Nob1p is required for cleavage of the 3' end of 18S rRNA[END_REF]. In lower eukaryotic species, the linker has a minimal length of ~15 residues while it exceeds ~100 residues in human. The conserved C-terminal linker region was shown to be critical for the Nob1/Pno1 complex formation in Chaetomium thermophilum [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF]. The presence of an active PIN domain required to cleave the pre-rRNA was demonstrated via structure prediction and mutational studies more than a decade ago while the precise function of the zinc ribbon is unknown [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF]. A recent review classifying the PIN domain-like superfamily [START_REF] Matelska | Comprehensive classification of the PIN domain-like superfamily[END_REF] defined the folding landscape of PIN domains as always adopting a prototypical β1-α1-α2-β2-α3-β3-α4-β4-α5-β5 topology [START_REF] Matelska | Comprehensive classification of the PIN domain-like superfamily[END_REF]. Their overall structure forms a five-stranded β-sheet packed by five α-helices with possible additional secondary structure elements although these extra elements do not change the global organization [START_REF] Matelska | Comprehensive classification of the PIN domain-like superfamily[END_REF]. Moreover, the four active site residues in the majority of the cases are located in the secondary structure element β1, α3, α4 and β4. Surprisingly, one of the amino acid predicted to be essential for the yeast Nob1 PIN domain activity, Asp110, appeared dispensable in yeast complementation experiments [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF]. In addition, PIN domains generally display low affinity for their target RNA, leading to the hypothesis that hNob1 requires additional stabilization upon interacting with its cleavage site on ribosomal rRNA. This hypothesis would explain the inability to reconstitute the cleavage reaction in vitro with hNob1 only [START_REF] Lamanna | Nob1 binds the single-stranded cleavage site D at the 3'-end of 18S rRNA with its PIN domain[END_REF]. Even if the zinc ribbon may be involved in stabilizing pre-rRNA contacts, one likely candidate for this role is the hNob1 partner hPno1. hPno1 contains two hnRNP K-homology (KH) domains and these domains are known to possess RNAbinding ability [START_REF] Valverde | Structure and function of KH domains[END_REF]. Additional proteins such as CINAP/Fap7 could also be implicated in this multi-factor reaction [START_REF] Hellmich | Essential ribosome assembly factor Fap7 regulates a hierarchy of RNA-protein interactions during small ribosomal subunit biogenesis[END_REF].

Here, we present a detailed molecular analysis of the interaction between the human transacting factor hNob1 and its binding partner hPno1. We revisit and propose a new delineation for the hNob1 PIN domain on the basis of biochemical data. We further demonstrate that the interaction between hNob1 and hPno1 requires the conserved sequence in hNob1 that precedes the second half of the PIN domain, along with the KH-like domain of hPno1. A single amino acid substitution abrogates the complex formation, supporting our finding that we have described the minimal interaction surface between the two proteins. We also observe that the non-conserved extensions found in hNob1 C-terminus and hPno1 N-terminus are, in our in vitro conditions, dispensable.

Results

In the last two decades, a number of studies have contributed to our understanding of Nob1's function in ribosome biogenesis. Most of our knowledge comes from seminal research performed by the group of D. Tollervey using the yeast protein [START_REF] Fatica | Nob1p is required for cleavage of the 3' end of 18S rRNA[END_REF]. Based on sequence alignment of the Nob1 protein family, it was proposed that Nob1 would contain a PIN domain and a zinc ribbon separated by a linker of variable length [START_REF] Fatica | Nob1p is required for cleavage of the 3' end of 18S rRNA[END_REF]. Almost a decade later, the NMR structure of the orthologous Nob1 protein from P. horikoshii was solved and the atomic data confirmed the presence of a PIN domain and of a zinc ribbon domain separated by a unstructured linker [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF].

Careful analysis of the Nob1 sequence alignment, 3D-structural modelling using the web server PHYRE2 [START_REF] Kelley | The Phyre2 web portal for protein modeling, prediction and analysis[END_REF] and recent PIN superfamily definition [START_REF] Matelska | Comprehensive classification of the PIN domain-like superfamily[END_REF] consistently suggest that the PIN domain present in the hNob1 protein (and likely in yeast) must be redefined. The domain should rather be divided in two pieces separated by a long linker (Figure 1A) [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF]. The first "half" would contain the required secondary structure elements β1-α1-α2-β2-α3-β3-α4 (β1-α1-α2-β2-α3-α4-β3-α5-α6 according to the predicted secondary structure elements present in the human sequence). The second part would cover the last elements β4-α5-β5, the missing elements required to form a functional PIN domain according to the recent classification (β4-α7-β5 according to our prediction; Figure 1A).

Interestingly, such a split organization changes the identity of the putative active site residues. When the PIN domain was defined with the yeast sequence, residues D15, E43, D84, D92 and D110 (residues D10, E36, D75, D83 and S101 in the human sequence) were proposed as conserved acidic residues possibly important for the enzymatic activity [START_REF] Fatica | Nob1p is required for cleavage of the 3' end of 18S rRNA[END_REF] [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF]. Since then, several studies have experimentally validated Asp15, Glu43 and Asp92 as critical for the yeast Nob1 activity, while Asp84 and Asp110 were apparently not required [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF], [START_REF] Pertschy | RNA Helicase Prp43 and Its Co-factor Pfa1 Promote 20 to 18 S rRNA Processing Catalyzed by the Endonuclease Nob1[END_REF]. Our predicted hNob1 PIN domain model fits the secondary structure requirements for a functional PIN domain of the VapC PIN superfamily (Figure 1B). According to our new definition, the residues D75 and S101 (corresponding to residues D84 and D110 in the yeast protein) are not predicted to be located in the active site (Figure 1B). This is likely the reason why the growth of S. cerevisiae cells was unaffected by the Nob1 D110N mutation [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF]. Based on our secondary structure prediction and sequence alignment, we predicted the β4-α5-β5 secondary structure elements between residues 234 and 255 (Figure 1A, labelled β4-α7-β5 in hNob1). With this revised arrangement, the fourth active site residues would correspond to D238 (D271 in the yeast protein, Figure 1A andB). Therefore, this model explains the published biochemical data (Figure 1B).

We next validated this computational prediction using co-expression experiments with multiple hNob1 constructs (Figure 1C). We first demonstrate that hNob1 fragment encompassing residues 1-134 can form a stable complex with hNob1 (206-314; Figure 1D). However, such a complex cannot be formed if we remove the stretch from residues 206 to 249, strongly suggesting that the interaction surface between the two proposed "halves" are located between residues 1 to 134 and the protein fragment 206 to 249 (Figure 1D). In order to challenge our hypothesis, we generated point mutations to destabilize the bi-partite PIN domain. Hence, D83K, L87S and Y91A point mutants targeting the interface between the two proposed PIN domain fragments were introduced based on the predicted model (Figure 1B and1E). hNob1 (1-134) carrying the mutation L87S is unable to form a stable complex with the hNob1 fragment 206-314 (Figure 1E) whereas D83K or Y91A mutations do not affect the complex formation. We thus show that the PIN domain of hNob1 is formed by two "halves" separated by a dispensable loop. Such a modular organization has been very recently confirmed by the cryo-EM based structural analysis of human pre-40S particles [START_REF] Ameismeier | Visualizing late states of human 40S ribosomal subunit maturation[END_REF].

To gain further insight into the human Nob1 and Pno1 interaction, we decided to precisely map the recognition surface. We performed co-expression experiments using full-length hPno1 with various truncations of hNob1. We initiated the experiment using the full-length proteins and showed an efficient and stable interaction despite limited degradation of hNob1 (Figure 2A). As a negative control, we also verified that untagged hPno1 alone did not show nonspecific binding to Ni-NTA beads (data not shown). Since the C-terminal part of hNob1 is predicted to be disordered and is not conserved through evolution (Figure 1A), we supposed that the observed degradation of the full-length (FL) protein occurred within this region. We next deleted the entire region that follows the zinc ribbon and reassessed the impact of this deletion on the stability of the complex.

The binding to hPno1 was unaffected by the truncation (construct 1-314, Figure 2A). Next, we further shortened the protein by removing the zinc ribbon domain (resulting in a construct from residues 1 to 255) and showed that it had no effect on hPno1 association (Figure 2A). Thus, the entire C-terminal region including the zinc ribbon is dispensable in vitro for complex formation.

We attempted to use the construct encompassing the first half of hNob1 PIN domain but we could not obtain the protein (data not shown). Instead, we used a construct containing only the linker located inside the PIN domain (construct 104-230) and co-expressed it with hPno1. Co-purification was observed for this construct, indicating that the interaction between hNob1 and hPno1 is at least partially mediated by the linker (Figure 2A).

We then turned our interest towards the binding site of hNob1 on hPno1. Pno1 is a conserved protein, which contains two KH domains [START_REF] Nicastro | RNA interactions: back in the groove[END_REF]. The first KH domain of hPno1 is termed "KH-like" in the literature because it does not contain the critical GXXG motif known to mediate ssRNA association [START_REF] Woolls | Roles of Dim2 in ribosome assembly[END_REF]. Over the years, several studies performed in yeast have looked into the association of Nob1 and Pno1 [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF], [START_REF] Woolls | Roles of Dim2 in ribosome assembly[END_REF]. In particular, Nob1 binding sites on the Pno1 protein were identified on the KH-like domain [START_REF] Woolls | Roles of Dim2 in ribosome assembly[END_REF]. Several mutations located at different places on the KHlike domain including the 310-helix fragment normally devoted to RNA association have been shown to modulate or abolish Nob1 association [START_REF] Woolls | Roles of Dim2 in ribosome assembly[END_REF]. We performed similar experiments with the human proteins by testing the capacity of the full-length hPno1 protein (hPno1 construct 1-252) to bind to hNob1 fragment encompassing the PIN domain (hNob1 construct 1-255; Figure 2B). In absence of pre-rRNA, the two proteins can form a stable complex (Figure 2B). We further tested the importance of the N-terminal non-conserved extension present in the hPno1 protein (residues 1 to 73) and showed that its deletion does not impair complex formation (Figure 2B). Similarly to hNob1, the evolutionary non-conserved region of hPno1 is dispensable for complex formation.

Finally, we truncated hPno1 to the KH-like domain (residues 73-155) and showed that it is sufficient to observe near WT-level of interaction with hNob1 (Figure 2B).

To map as precisely as possible the site of interaction between the two proteins, we decided to find the responsible residues on each protein. We previously observed that the internal loop present in the PIN domain is necessary and sufficient for a stable association of hPno1 (construct 104-230, Figure 2A). The loop was artificially divided in two fragments based on sequence composition (Figure 1A) and used to test hPno1 association by pull-down experiments (Figure 3).

The fragment from 104 to 179 is highly enriched in charged residues and poorly conserved in evolution, while the sequence of the fragment 180 to 220 is highly conserved through evolution (Figure 1A). To facilitate the stable expression of the linker sequences, we fused these two sequences to a glutathione S-transferase tag. After purification of the GST-fused sequences, we incubated them with the N-terminal deleted hPno1 construct (Figure 3B). The C-terminal part of the linker (residues 180-220) is able to retain hPno1 on the beads while the fragment 104-179 fails to interact (Figure 3B). These data therefore suggest that the PIN domain, the zinc ribbon and the large unstructured C-terminal region are dispensable for the interaction with hPno1. Based on sequence conservation (Figure 1A), we mutated the W208 in the core of the highly-conserved linker region and reassessed complex formation by co-expression of the mutated hNob1 together with hPno1. Mutant protein was readily expressed but did not stably bind to hPno1 (Figure 3C).

We thus concluded that the formation of the human complex relies on a short sequence motif found in the loop present between the bipartite PIN domain, a situation reminiscent to the one recently identified in the fungal complex (Figure 3) [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF].

Finally, we turned to hPno1 and using a sequence conservation mapping performed by the Consurf server [START_REF] Ashkenazy | an improved methodology to estimate and visualize evolutionary conservation in macromolecules[END_REF] identified a hot spot of conservation located between residues 105 to 110 (data not shown). Unfortunately, point mutations in this region led to poor stability of solubility of hPno1, which prevented us from testing the capacity of the mutated hPno1 protein to bind to hNob1. Recently, the same peptide was shown to mediate numerous contacts with hNob1 [START_REF] Ameismeier | Visualizing late states of human 40S ribosomal subunit maturation[END_REF].

Discussion

Nob1/Pno1 is the key complex responsible for the final maturation step of the 18S pre-rRNA. Despite numerous studies on the yeast proteins, information is still lacking on their exact function and on the formation of the complex especially within the context of human ribosome biogenesis. With the aim to fill this lack of knowledge, we first observed that the boundaries of the PIN domain defined a decade ago for hNob1 should be modified to reflect the published biochemical data [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF], [START_REF] Fatica | PIN domain of Nob1p is required for D-site cleavage in 20S pre-rRNA[END_REF]. We demonstrate that the PIN domain extends further into the hNob1 primary sequence by using bioinformatics and biochemical methods (Figure 1). With this new organization of the PIN domain, a long partially conserved sequence has now become an internal loop of the PIN domain corresponding to approximately residues 96 to 232. However, the unstructured loop may not be a simple "separator" between the newly defined bi-partite PIN domain.

We then set out to define the determinant of the interaction within the hNob1 / hPno1 complex by using co-expression and pull-down assays. Several deletion or point mutation constructs were generated and used to characterize in vitro the minimal protein fragments compulsory to observe a stable complex. Using this approach, the region of contact between the two proteins could be reduced to a short linear fragment located between residues 206 and 220 for hNob1 and to the KH-like domain of hPno1. Further fine mapping of hPno1 was precluded by an inability to obtain stable mutant protein using our experimental setup.

By redefining the limits of the PIN domain of hNob1 and the location of the interaction surface with hPno1, we now postulate that hPno1 may regulate the endonuclease activity of hNob1.

It is possible that the binding of Pno1 sterically hinders the PIN domain active site or modulates the accessibility of the PIN domain for its substrate 18S pre-rRNA. Recent data obtained with the C. thermophilum Nob1 and Pno1 proteins demonstrate that co-expressing ctPno1 with two ctNob1 constructs encompassing residues 1-169 and 250-354 can form a stable complex [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF]. The two ctNob1 constructs can be found in the same complex either by sandwiching ctPno1 or via a direct interaction between them. These constructs correspond to residues 1-115 and 188-291 in hNob1.

According to our experiments, the construct 1-115 covers the first half of the PIN domain, thus insufficient to bind stably to hPno1 in vitro (Figure 3B). The second construct, i.e. residues 188-291, contains the other half of the PIN domain and the zinc ribbon. Taken together, our data help explain why ctNob1/ctPno1 tri-partite complex has been reconstituted by Sturm and colleagues [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF]. Their ctNob1 constructs effectively covered the two halves of the PIN domain, an observation undetected at that time.

In our hands, the zinc ribbon does not interact stably with the first half of the PIN domain (Figure 1D). Although we did not observe any decrease in complex formation when deleting this region from hNob1 in vitro (Figure 2A), it is possible that in vivo, the zinc ribbon and/or the Cterminal region participate in complex formation with Pno1. In addition to potentially helping the interaction with hPno1, the close proximity between the zinc ribbon and the active site of the PIN domain makes us suggest that the zinc ribbon could be important to regulate the hNob1 endonucleolytic activity.

In conclusion, the human protein hPno1 binds to hNob1 using a similar surface as compared to the yeast counterpart [START_REF] Woolls | Roles of Dim2 in ribosome assembly[END_REF]. hNob1 contacts hPno1 using residues between 206 and 220, where the key amino acid W208 is located [START_REF] Sturm | Interdependent action of KH domain proteins Krr1 and Dim2 drive the 40S platform assembly[END_REF], suggesting highly conserved mechanism between species even though each organism presents some specificities due to sequence evolution. During the preparation of this manuscript, Ameismeier and colleagues reported the cryo-EM based model of several late stage pre-40S particles either purified from native source or reconstituted where they validate the present conclusion on the interaction between hNob1 and hPno1 [START_REF] Ameismeier | Visualizing late states of human 40S ribosomal subunit maturation[END_REF].

A key question remains regarding hNob1's activity on the 18S pre-rRNA. Indeed, information on hNob1 recognition and cleavage of 18S pre-rRNA are missing since hNob1's active site is located at more than 30 Å away from the newly defined PIN domain active site [START_REF] Ameismeier | Visualizing late states of human 40S ribosomal subunit maturation[END_REF]. The role of individual components in 18S pre-rRNA cleavage may in future be provided by additional structures of pre-40S intermediates or by individual protein complexes together with their target RNA. Alternatively, one cannot exclude that additional players such as the CINAP/Fap7 protein, may be required to induce the obligatory local rearrangements leading to the close proximity between hNob1 PIN domain active site and pre-rRNA allowing the direct observation of the site 3 cleavage configuration. More data will therefore be necessary to fully comprehend the role of the Nob1/Pno1 complex during this critical step of ribosome production.

Materials and Methods

PCR amplification and cloning

Plasmids expressing hNob1 and hPno1 were built by inserting the PCR-amplified sequences of hNob1 or hPno1 into modified pET vectors allowing for coexpression and with or without an N-terminal hexa-Histidine or Glutathione S-Transferase (GST) Tag. Cloning was performed with a modified sequence ligation independent cloning (LIC) protocol and plasmids were verified by DNA sequencing.

Protein expression and purification

hNob1 and hPno1 constructs were individually transformed or cotransformed into Escherichia coli BL21 Rosetta 2 cells. Bacteria were cultured in TB medium at 37 °C and protein expression was induced with 0.2 mM isopropyl 1-thio-β-D-galactopyranoside overnight at 15 °C. Cells were pelleted down by centrifugation and resuspended with 1 mL of ice cold Lysis buffer (50 mM Tris-HCl pH 8, 300 mM NaCl, 0.1 % Triton X100, 5 mM β-mercaptoethanol and 10 % glycerol w/v). Cells were lysed by sonication. The soluble fraction was separated from the crude extract by centrifugation at 14,000 rpm for 30 min at 4 °C. Cleared lysate was incubated with 50 µL Ni-NTA beads SIGMA (H0537) for 2 h at 4 °C on a rotating wheel. Beads were spun down to remove the unbound material and washed five times with Wash buffer (50 mM Tris-HCl pH 8, 300 mM NaCl, 20 mM imidazole, 5 mM β-mercaptoethanol and 10 % glycerol w/v). His-tagged proteins were eluted by competition with Elution buffer (Wash buffer plus 250 mM imidazole). Eluted samples were analyzed by SDS PAGE and stained with Coomassie blue. GST-fused proteins were purified on GST beads SIGMA (G4510) following same purification protocol with the exception that the Wash buffer did not contain imidazole. Protein elution was carried out with Wash buffer supplemented with 20 mM reduced glutathione.

Pull down assay

Proteins were purified as indicated above. Histidine-tagged full length hPno1 and a construct from residues 73-252 were subjected to a cleavage step with TEV protease during an overnight incubation at 16 °C. Cleaved tag, TEV proteases and undigested proteins were removed by reloading the sample on Ni-NTA and collecting the flow-through fraction. Tagged hNob1 constructs were added to hPno1 in a 1:1 molar ratio. The mixture volume was adjusted to 1 mL with binding buffer (50 mM Na2HPO4 / NaH2PO4, pH 7.5, 500 mM NaCl, 5 mM βmercaptoethanol and 10 % glycerol w/v). Samples were incubated at 4 °C for 1 h and afterwards 30 µL of NiNTA beads was added to the mixture and further incubated on a rotating wheel for 1 h at 4 °C. Unbound material (flow-through) was collected after centrifugation for 30 s at 2,000 g.

Beads were washed five times with Wash buffer 2 (50 mM Na2HPO4 / NaH2PO4, pH 7.5, 500 mM NaCl, 20 mM imidazole, 5 mM β-Mercaptoethanol and 10 % glycerol w/v) before elution with 30 µL of Elution buffer 2 (Wash buffer 2 supplemented with 250 mM imidazole). Input and bound fractions were analyzed by SDS PAGE and stained with Coomassie blue. GST pull-downs were performed similarly except that buffers did not contain imidazole and the elution buffer was supplemented with 20 mM of reduced glutathione.

Bioinformatics analysis of hNob1 or hPno1

Protein sequences for hNob1 (UniProt entry Q9ULX3) and hPno1 (UniProt entry Q9NRX1) were uploaded to the online server PHYRE2 for structure prediction [START_REF] Kelley | The Phyre2 web portal for protein modeling, prediction and analysis[END_REF] and to the ConSurf server for sequence conservation prediction [START_REF] Ashkenazy | an improved methodology to estimate and visualize evolutionary conservation in macromolecules[END_REF].

The structural model obtained for hNob1 was used to prepare the panel B from Figure 1 and subsequently design the truncation/point mutation boundaries. The most closely related model found by PHYRE is the atomic structure of Nob1 from P. horikoshii (PDB code 2LCQ) [START_REF] Veith | Structural and functional analysis of the archaeal endonuclease Nob1[END_REF].

The sequence conservation score within the Pno1 protein family was calculated using 150 sequences from homologous proteins taken in the UniRef90. The conservation score was then mapped onto the surface of the archaeal Pno1 structure from P. horikoshii (PDB code 2E3U) [START_REF] Jia | Crystal structure of Dim2p: a preribosomal RNA processing factor, from Pyrococcus horikoshii OT3 at 2.30 A[END_REF].

  Figure legends.

Figure 1 .

 1 Figure 1. Human Nob1 domain organization. A) Sequence alignment and secondary structure

Figure 2 .

 2 Figure 2. Definition of the minimal fragment of hNob1 or hPno1 required for complex formation.

Figure 3 .

 3 Figure 3. hNob1 interaction with hPno1 depends on the tryptophan at position 208 in hNob1 and

Acknowledgements

We thank Cameron Mackereth and Fabien Darfeuille for critical reading of the manuscript. We further acknowledge Lucile Maleuvre and Manon Coursières for participation to the experiments. This work was funded by ANR grant (ANR Blanc 2010 grant RIBOPRE40S (ANR-10-BLAN-1224) and RIBOMAN ANR-16-CE11-0029 to S.F.), the INSERM and the University of Bordeaux.


Lotfi Belkoura
email: lotfi.belkoura@univ-lille1.fr

Jean-Pierre Richard
email: jean-pierre.richard@ec-lille.fr

A DISTRIBUTION FRAMEWORK FOR THE FAST IDENTIFICATION OF LINEAR SYSTEMS WITH DELAYS

Keywords:

extended to delays identification.

INTRODUCTION

The delay phenomenon constitutes one of the major complexity components of networked control, since actuators, sensors, computers, field networks and wireless communications that are involved in feedback loops unavoidably introduce variable deadtimes. Despite numerous advances in this field, delay remains a theoretical and practical challenge (see, e.g. the survey [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]). Among these open problems, the online delay identification has a particular part: On the one hand, the delay knowledge can benefit many of the control techniques suited for such systems (predictors, finite spectrum assignments, observers...). On the other hand, the existing identification techniques for time delay systems (TDS) (see [START_REF] Ren | On line identification of continuoustime systems with unknown time delay[END_REF] for a modified least squares technique and a survey in [START_REF] Belkoura | Identifiability and identification of linear systems with delays[END_REF])) generally suffer from poor speed performance. The authors regard this paper as a first step to the TDSadaptation of the fast identification techniques that were recently proposed in (Fliess M., 2003) in the framework of nonlinear, finite-dimensional models. Let us recall that those techniques are not asymptotic, and do not need any statistical knowledge of the noises corrupting the data. (See, e.g., (Fliess and Sira-Ramirez, 2004) for applications to nonlinear state estimation, (Fliess et al., 2004), [START_REF] Fliess | Closed-loop fault-tolerant control for uncertain nonlinear systems. Control and Observer Design for Nonlinear Finite and Infinite[END_REF] for linear and nonlinear diagnosis, and (Fliess et al., 2003) for signal processing.

We adopt in this paper a distributional formulation from which the parameters as well as the input delays can be easily identified. The plant is assumed to be linear with respect to its parameters, but not necessarily with respect to its inputs and outputs. The identification procedure mainly consists in the following three steps: (1) Differentiations, (2) Multiplication with an appropriate C ∞ function and (3) Integrations.

We first recall here some standard definitions and results from distribution theory and fix the notations we shall use in the sequel. Let Ω be an open subset of R. The space of C ∞ functions having compact support in Ω is denoted by D(Ω), and D ′ (Ω) is the space of distributions on Ω, i.e. the space of continuous linear forms on D(Ω).

The complement of the largest open subset of Ω in which a distribution T vanishes is called the support of T and is denoted by supp T . D ′ + (resp. E ′ ) is the space of distributions with support contained in [0, ∞) (resp. compact support). It is an algebra with respect to convolution with identity δ, the Dirac distribution. When concentrated at a point {τ }, the latter distribution δ(t -τ ) is denoted δ τ . A distribution is said to be of order r if it acts continuously on C r -functions but not on C r-1 -functions. Measures and functions are of order 0. Functions are considered through the distributions they define and are therefore indefinitely differentiable. If y is a continuous function except at a point a with a finite jump σ a , its derivative dy/dt writes

dy/dt = ẏ + σ a δ a ( 1
)
where ẏ is the distribution defined from the usual derivative of y. Derivation, integration and translation can be formed from the convolution products

ẏ = δ (1) * y, y = H * y, y(t-τ ) = δ τ * y (2)
where δ (1) is the derivative of the Dirac distribution, and H denotes the Heaviside function. With a slight abuse of notations, we shall write H k y the iterated integration of y and more generally T k the iterated convolution product of order k.

For S, T in D ′ + , one has: supp S * T ⊂ supp S + supp T

(3)

where the sum in the right hand side is defined by {x + y ; x ∈ supp S, y ∈ supp T }. Finally, when there is no danger of confusion, we shall sometimes denote T (s), s ∈ C, the Laplace transform of T .

The paper is organized as follows. Section 2 gives the main theoretical tools we shall need for our identification problem. Section 3 is devoted to specific cases involving piecewise constants inputs, and Section 4 considers the general case with arbitrary inputs.

MULTIPLICATION AND CONVOLUTION PRODUCTS

In the general case, the multiplication of two distributions (say α and T ) is not always defined. However, this operation always make sense when one of the two terms is a smooth function. Several properties can be derived from such product and the most important for our developments are given bellow. The next Theorem is the key result from which most of the parameters (including the delays) can be identified from step input responses.

Theorem 1. [START_REF] Schwartz | Thorie des distributions[END_REF] If T has a compact support K and is of order m (necessarily finite), α T = 0 whenever α and its derivatives of order ≤ m vanish on K.

The following examples illustrate this statement in case α is a polynomial and T a singular distribution. Note that, in forming the product α T , the delay τ involved in the argument T (t -τ ) now appears also as a coefficient.

t δ = 0, (t -τ ) δ τ = 0, t 2 (t -τ ) (a δ (1) + b δ τ ) = 0.
The usual derivation rule (α T ) ′ = α ′ T +α T ′ also applies. Thanks to the previous Theorem and the Leibniz's formula, the statement for the product t δ can be generalized to t l δ (n) = 0 for l > n, and

t l δ (n) = (-1) l (n!/(n -l)!) δ (n-l) , l ≤ n. (4)
We shall make use of another property involving both multiplication with t n and the convolution product, in case one of the two distributions (S or T ) has a compact support.

t n (S * T ) = n k=0 C k n (t k S) * (t n-k T ). (5)
The combination of the rules ( 4) and ( 5) with S = δ (p) and T = y allows us to transform terms of the form t n y (p) into linear combinations of derivatives of products t k y. Denoting z i = t i y, one has for example,

t 3 y (2) = t 3 (δ (2) * y) = -6 z 1 + 6 z (1) 2 -z (2)
3 . (6) Note that integrating twice this expression by considering H 2 t 3 y (2) results in nothing but the integration by parts formula with available data z i . When facing delayed terms, relation (5) also applies and leads to expressions involving the delay not only as an argument, but also as a unknown coefficient to be identified. For example, one easily gets

t 2 ẏ(t -τ ) = t 2 (δ (1) τ * y) = δ τ * [τ 2 z (1) 0 -2τ (z 0 -z (1) 1 ) -2z 1 + z (1) 2 ](7)
Here again, integration with H yields a causal relation with available data z i .

IDENTIFICATION FROM STEP INPUT RESPONSES

Using simple examples, this section shows hows Theorem 1 can be used for both parameters and delays identification in case of piecewise constant inputs.

3.1 First order system with a single input delay

Let us consider a first order system with a delayed input governed by:

ẏ + a y = y(0) δ + γ 0 H + b u(t -τ ), (8)
where γ 0 is a constant perturbation, a, b, and τ are constant parameters, and the coefficient a is assumed to be known (for the moment). Consider also a step input u = u 0 H. A first order derivation yields

ÿ + a ẏ = y(0) δ (1) + γ 0 δ + b u 0 δ τ . (9)
By virtue of Theorem 1, the right hand side of equation ( 9) can be canceled by means of a multiplication with a function α such that

α(0) = α ′ (0) = 0, α(τ ) = 0, ( 10
)
and the choice of the polynomial α(t) = t 3 -τ t 2 results in

t 3 (ÿ + a ẏ) = τ t 2 (ÿ + a ẏ). ( 11
)
As an equality of singular distributions, this relation doesn't make sense for any t (otherwise we would have τ = t). However, k ≥ 1 successive integrations (or a convolution with H k ) result in functions equality from which the delay τ becomes available. More precisely, since supp H k δ τ ⊂ (τ, ∞), we can easily show that all the obtained functions will vanish on (0, τ ) and the delay is consequently not identifiable on this interval. Conversely, being nonzero for all t > τ , the delay is everywhere identifiable on (τ, ∞). We therefore get from (11):

τ = H k (t 3 y (2) + a t 3 y (1) ) H k (t 2 y (2) + a t 2 y (1) ) , t > τ, (12)
where, by virtue of (5) and recalling the notation z i = t i y of the previous section,

t 3 y (2) = -6 z 1 + 6 z (1) 2 -z (2) 3 , t 2 y (2) = -2 z 0 + 4 z (1) 1 -z (2) 2 , t 2 y (1) = 2 z 1 -z (1) 2 , t 3 y (1) = 3 z 2 -z (1) 3 .
These coefficients show that k ≥ 2 integrations can avoid any derivation in the delay identification algorithm, and the next figure show a partial realization scheme (limited to H k t 3 y (2) ) of the terms involved in (12).

Note that the algorithm given in (12) only requires the values of a and the output y. In case the coefficient a is also unknown, the same approach can be used for a simultaneous identification. The following relation is easily derived from ( 11)

τ (t 2 y (2) ) + a τ (t 2 y (1) ) -a (t 3 y (1) ) = t 3 y (2) (13) + 1 s 1 s 1 s 1 s 1 s y -6 6 -1 H 2 t 3 y (2) H 3 t 3 y (2) H 4 t 3 y (2) t t 2 t 3 Fig. 1. Realization scheme of H k t 3 y (2)
and a linear system with unknown parameters (τ, a τ, a) is obtained by using different integration orders

  H 2 w 1 H 2 w 2 H 2 w 3 H 3 w 1 H 3 w 2 H 3 w 3 H 4 w 1 H 4 w 2 H 4 w 3     τ aτ -a   =   H 2 w 0 H 3 w 0 H 4 w 0   (14)
where we have denoted w 0 = t 3 y (2) , w 1 = t 2 y (2) , w 2 = t 2 y (1) , and w 3 = t 3 y (1) . A simulation result from equation ( 14) is given in Figure 2 for k = 2, 3, 4 integrations and the parameters

y(0) = 0.3, a = 2, τ = 0.6, γ 0 = 2, b = 1, u 0 = 1.
For the previous identifiability reason, the obtained linear system may be not consistent for t < τ . Moreover, and unlike the single delay case, a local loss of identifiability may occur for t > τ as suggested in Figure 2 for t ≈ 1.5s. For a complete identifiability study, the reader may refer to (Fliess M., 2003) for nonlinear and free delay cases, and [START_REF] References Belkoura | Identifiability of systems described by convolution equations[END_REF] Note that this simultaneous identification of a and τ can also be formulated as a generalized eigenvalue problem as

[ H 2 w 3 H 2 w 0 H 3 w 3 H 3 w 0 -τ H 2 w 2 H 2 w 1 H 3 w 2 H 3 w 1 ] a 1 = 0 (15)
Finally, others combinations of parameters can be obtained and identified by means of a multiplication of (9) with different C ∞ functions. For instance, considering α(t) = t(t -τ ) and using the same previous technique results in

H k t(ÿ + a ẏ) τ + H k τ y(0) = H k t 2 (ÿ + a ẏ) (16)
from which both delay and initial condition can be identified. The only coefficient for which the explicit value of τ is required is its associated parameter b. Due to the fast convergence of the algorithms, we may also consider a separate procedure in which the undelayed terms are firstly identified and reused for the delay identification.

Second order system with multiple input delays

The previous technique can be extended to multiple input delays by means of one of the three different following approaches: (a) A multiplication with the

C ∞ function t 2 (t -τ 1 ) • • • (t -τ K ) (if K is finite)
. This may however lead to a large size linear system for which all the delays are not identifiable until t > τ K , (b) A recursive identification whose main drawback comes from propagation errors, and (c) A local (in time) identification if one assumes an incompressible delay such that

τ k+1 -τ k > ∆. ( 17
)
We shall focus in this section on delays identification based on the third case and consider the example of a second order system governed by

ÿ + a 1 ẏ + a 0 y = ϕ 0 + γ 0 H + 4 k=0 b k u(t -τ k ) ( 18
)
where τ 0 = 0, γ 0 is a constant perturbation, and ϕ 0 (of order 1 and support {0}) contains the initial condition terms. A first order derivation, a multiplication by t 3 (t -λ), and a "local integration" using H(s) = (1-e -∆s/3 )/s (instead of 1/s) result in

0 = H 3 t 3 (t -λ)(y (3) + a 1 y (2) + a 0 y (1) ) = N 1 -λ N 0 . ( 19
)
Terms of the form H 3 t p y (q) (and hence N 0 and N 1 ) are formed as in the previous section, while from the right hand side of ( 18) we also have,

N i = H 3 4 k=1 b k τ 3+i k δ τ k , i = 0, 1. ( 20
)
Since supp H 3 ∈ (0, ∆), and by virtue of (3), we get that N 0 = N 1 = 0 on R\(τ k , τ k + ∆) while on the remaining subsets, the delays are obtained from

λ = τ k = N 1 /N 0 on (τ k , τ k + ∆). ( 21
)
A simulation result is shown in Figure 3 for the parameters , 4, 6, 9). 21) is set to zero until N 0 reaches a significant nonzero value.

a 0 = 2, a 1 = 1, y(0) = 1.3, ẏ(0) = -2.3, γ 0 = 0.5, ∆ = 0.25, (b 0 , • • • , b 4 ) = (0.5, -5, 2, 1), (τ 1 , • • • , τ 4 ) = (2

A particular case with state delay

For some particular cases, a local state delay identification based on step input response is possible as well. The following example illustrates the proposed approach.

ẏ(t) + a y(t -τ ) = b u(t), (22) y
(θ) = 0 θ ∈ (-τ, 0).
Denoting H τ = H(t -τ ), the left hand side member writes (δ + a H τ ) * ẏ in which δ + a H τ admits the convolution inverse [START_REF] Hirsh | Elements of functional analysis[END_REF],

(δ + a H τ ) * -1 = δ -a H τ + • • • + (-1) n-1 a n-1 H n-1 τ + • • • . (23)
Since, by virtue of (3), supp H k τ ⊂ (k τ, ∞), a second order derivation of (22) combined with (23) results in

y (3) = b δ (1) -a b δ τ + ζ (24)
where supp ζ ⊂ (2 τ, ∞).

Following the line of the previous sections, a cancellation of the singularities using the smooth function t 2 (t -τ ) yields

t 2 (t -τ )(y (3) -ζ) = 0. (25)
Therefore, when restricted to the interval (τ, 2 τ ) the state delay is obtained from

τ = (H k t 3 y (3) )/(H k t 2 y (3) ), t ∈ (τ, 2 τ ). ( 26
)
The next simulation shows the realization of ( 26) with a = 2 and τ = 0.6.

As expected, the delay is identified only on the interval (τ, 2τ 2τ . Note also that unlike the input delay case, a simultaneous identification including the state coefficient a can not be considered in this case.

IDENTIFICATION WITH ARBITRARY INPUTS

For systems with input/state delays and general inputs, there are mainly two limitations of different nature. The first one concerns the online aspects while the second one is more fundamental. More precisely, the delays can be generally identified by means of convolution products and this may constitute a drawback if online results are expected. However, this constraint can be relaxed if there is no state delay and if the input is not prescribed. On the other hand, for a state delay case and unless particular situations, the initial condition (involving a function) is by definition unknown and depends, via its support within (-τ, 0), on the delay to be identified.

These observations as well as the identification procedure are illustrated in the following first order example ẏ +a 0 y +a 1 y(t-τ ) = ϕ 0 +b 0 u+b 1 u(t-τ ), (27) in which the initial condition ϕ 0 writes

ϕ 0 = y(0) δ -a 1 ϕ, supp ϕ ∈ (0, τ ),
for some function ϕ. For simplicity reasons and although a simultaneous identification of both the parameters and the delay can be considered, we shall restrict our study to the delay identification.

For notational convenience, let us define from ( 27) the delayed and undelayed terms as

η ẏ + a 0 y -b 0 u, β b 1 u -a 1 y. ( 28
)
Using the developments given in Section 2-eq.( 7), a first order derivation of ( 27) followed by a multiplication by t and t 2 results in

t(η -ϕ 0 ) = t(δ τ * β) = δ τ * (τ β + tβ), (29)
t 2 (η -ϕ 0 ) = t 2 (δ τ * β) = δ τ * (τ 2 + 2τ tβ + t 2 β).
Now, a cross convolution product of these two relation yields the following convolution equation involving τ , τ 2 , and in which the shift operator (δ τ * ) has been removed.

τ 2 (tβ * tη) + τ (2tβ * tη -β * t 2 η) = (t 2 β * tη -tβ * t 2 η) (30) +(t 2 β + 2τ tβ -τ 2 tβ) * tϕ 0 + (τ β -tβ) * t 2 ϕ 0 .
It remains to consider the initial condition term ϕ 0 = y(0) δ -a 1 ϕ that appears in this equation.

If ϕ ≡ 0 (the plant is initially at rest) or a 1 = 0 (there is no state delay), then tϕ 0 = t 2 ϕ 0 = 0 and successive integrations lead to the linear system

Hw 1 Hw 2 H 2 w 1 H 2 w 2 τ 2 τ = Hw 0 H 2 w 0 , (31)
where we have denoted

w 1 = tβ * tη, w 2 = 2tβ * tη -β * t 2 η, w 0 = t 2 β * tη -tβ * t 2 η. (32)
Note that by means of successive derivations, the same technique can be used in case of constant or piecewise polynomial function ϕ. On the other hand, in the free state delay case (a 1 = 0 ⇒ β = b 1 u), and if the input is free, the above convolution products can be easily realized online by considering simple expressions for u(s). For instance, with u(t) = He -at , and recalling that du/ds is the Laplace transform of tu, the coefficient w 1 in (32) also writes

w 1 (s) = b 1 (tη)(s) (s + a) 2 . ( 33
)
Figure 5 shows a simulation result with τ = 2, a 1 = 0, a 0 = 2, y(0) = 0 b 0 = b 1 = 1, and the discontinuous input u = H(t) -3H(t -4) + cos(t)(0.2+sin(2t)). This correspond to a single input delay which does not require a multiplication by t 2 . Note that a singularity occurred at t ≈ 8.5s. In Figure 6, the input/state delay is considered with a 1 = 0.4, a 0 = 2, b 0 = b 1 = 1, τ = 2, and the input u = cos(t)(0.2 + sin(2t)). The initial condition ϕ is set to zero although y(0) = 0.5 = 0. Here again, the small deviations observed for t ≈ 4.5s may result from a local ill-conditioned linear system. This note has presented a simple method for both parameter and delay identification of dynamical systems. As in (Fliess M., 2003), the high speed convergence of the proposed algorithms allows us to consider the problem of on line identification and control of time delay systems. Robustness issues including the noise effects and a deeper study of the possible singularities, multivariable systems with partial state measurements as well as the extension to discrete time processes are under active investigation.

Fig. 3 .

 3 Fig. 3. Multiple delays identificationDue to the non identifiability of the delays on R\(τ k , τ k + ∆), the value of λ in (21) is set to zero until N 0 reaches a significant nonzero value.

Fig. 4 .

 4 Fig. 4. Local identification of a state delay

Fig. 5 .

 5 Fig. 5. Input delay identification

  Fig. 6. Input/State delay identification

This work was supported by the Rgion Nord Pas-de-Calais and the FEDER Fonds Europen de Dveloppement Rgional (European Funds of Regional Development) under the ARCIR Robocoop and the AUTORIS-TAT T31 project.


Chengfang Ren

Alexandre Renaux

Sylvain Azarian

Performance Bounds for Change Point and Time Delay Estimation

Keywords: Performance bounds, Change point estimation, Time delay estimation

 with an additional change on the transmitted message. A theoretical derivation of the Barankin bound and a simplified version of this bound are proposed in order to predict the Maximum Likelihood Estimator (MLE) behavior. Simulations illustrate the validity of both bounds but it is pointed out that the MLE seems asymptotically efficient for the normalized time shift estimation contrary to the estimation of message's starting time.

I. INTRODUCTION

Time delay estimation between two spatially separated sensors is fundamental for localization and tracking problems in radar, acoustics and ultrasonic applications [START_REF]A comparative study of two discrete-time phase delay estimators[END_REF], [START_REF] Carter | Coherence and time delay estimation[END_REF]. It turns out that this problem could also be interesting for physically separated sensors synchronization [START_REF] Zhong | Time delay estimation in the presence of clock frequency error[END_REF]. Previous studies on this problem always assume that the shape of reference signal did not change over time. But this hypothesis might not be verified for an opportunity source and only a portion of the reference signal can be exploited. To solve this problem, our method will be based on retrospective change point estimation with an additional parameter (time delay) to estimate. Change point estimation/detection problems have been widely investigated in the literature e.g. [START_REF] Basseville | Detection of Abrupt Changes, Theory and Application[END_REF] for a general introduction and [START_REF] Chen | Parametric Statistical Change Point Analysis[END_REF] for the specific offline change point estimation. In terms of performance analysis, the first asymptotic results concerning the MLE for very specific problems have been proposed in [START_REF] Hinkley | Inference about the change-point in a sequence of random variables[END_REF] and recently reanalyzed in [START_REF] Fotopoulos | Exact asymptotic distribution of change-point MLE for change in the mean of Gaussian sequences[END_REF] and [START_REF] Fotopoulos | Maximum likelihood estimation of a change-point for exponentially distributed random variables[END_REF]. The term "asymptotic" generally refers to the case where a large number of observations are collected before and after the change and then the asymptotic behavior of MLE in terms of probability density function (pdf) can be derived for most of the time. However, from a practical point of view, these results are currently hard to be useful. This is why the signal processing community has tried to bypass this drawback by focusing on lower bound on the Mean Square Error (MSE) and not directly on the pdf of the MLE.

In the context of change point estimation, the discrete nature of parameters to estimate implies that the regularity conditions are not satisfied in order to apply the classical Cramér-Rao bound [START_REF] Cramér | Mathematical Methods of Statistics, ser[END_REF], [START_REF] Rao | Information and accuracy attainable in the estimation of statistical parameters[END_REF] similarly to the case of cyclic parameters estimation [START_REF] Routtenberg | Non-Bayesian periodic Cramér-Rao bound[END_REF] or when the support of data pdf depends on the parameters to estimate [START_REF] Lu | The multidimensional Cramér-Rao-Leibniz lower bound for likelihood functions with parameter-dependent support[END_REF]. Then, other bounds requiring less regularity conditions have been proposed to overcome this difficulty. The Chapman-Robbins bound has been derived in [START_REF] Ferrari | Barankin lower bound for change points in independent sequences[END_REF] in the context of one change point and extended to the multiple change point problem in [START_REF] La Rosa | Barankintype lower bound on multiple change-point estimation[END_REF]. In the Bayesian context, the Weiss-Weinstein bound has been studied in [START_REF] Bacharach | Weiss-Weinstein bound for change-point estimation[END_REF]- [START_REF] Bacharach | A Bayesian Lower Bound for Parameter Estimation of Poisson Data Including Multiple Changes[END_REF]. In this paper, we stay in the non-Bayesian context to study the Barankin (or McAulay-Seidman bound) [START_REF] Barankin | Locally best unbiased estimates[END_REF], [START_REF] Mcaulay | A useful form of the Barankin lower bound and its application to PPM threshold analysis[END_REF] for a change point estimation problem when (contrary to previous works) two sets of non synchronized data are available. Consequently, one has to additionally estimate the time delay between both receivers. We first start by presenting the mathematical framework of our analysis. Then, we derive the Barankin bound and propose also a simplified hybrid bound which is found to be surprisingly simple to compute. Finally, some simulation results are presented in order to show the interest of the proposed bounds.

II. MODEL SETUP

The application context of our paper is trying to synchronize multiple receivers on a well known message m (t) (sinus, chirp, etc...) sent by a well located opportunity source e (t) (could be a radio station, communication channel, etc...). For readability, we restrict our analysis to two receivers denoted r 1 and r 2 . The geometry of both emitter and receivers are known, so we can easily compensate the radial distance difference between e to r 1 and e to r 2 by introducing a time delay on the received signals. Nevertheless, in order to simplify the analytic expression of the proposed bound, these two sensors are assumed to be equidistantly located from the opportunity source. The main goal here is to take advantage of this well known message m (t) in order to synchronize our receivers. Therefore, the synchronization tasks consist in three steps: detect m (t) then estimate exactly the starting time of this message. This time will be denoted by τ in the following. And finally, use this message to synchronize r 1 and r 2 by estimating the normalized time shift denoted by ∆ between the received signals which could be viewed as the time delay between receivers. Since our goal is trying to establish performance bounds on the estimation of τ and ∆, we will skip the detection part. Mathematically, we can formulate the estimation problem as following:

Assume that we have already detected a single message m (t) received by r 1 and r 2 between the time interval [0; N T e ] where T e is the common sampling interval of r 1 and r 2 . We only assume phase jump between receivers' clock and no frequency drift, otherwise it will be useless to synchronize them. Additionally, the sampling interval should respect the Nyquist-Shannon sampling theorem, i.e. 2BT e < 1, where B is the bandwidth of m(t). This condition is sufficient for not loosing information on m(t) when sampling. Therefore, we can model the received signal as

r 1 [k] = s (kT e ) + n 1 [k], k = 1, . . . , N (1)
for the first receiver and

r 2 [k] = s (kT e + ∆T e ) + n 2 [k], k = 1, . . . , N (2)
for the second receiver where s (t) is the noiseless signal of form

s (t) = 0, if t ≤ τ m (t) , if t > τ , , (3)
N is the total number of observations, ∆T e represents the fact that these two receivers are not synchronized (0 < ∆ < 1) and the receivers noise are represented by n i [k], i = 1, 2, assumed to be independent and identically distributed (i.i.d.). This model is similar to the classical time delay estimation problem proposed in [START_REF] So | Time-delay estimation for sinusoidal signals[END_REF]- [START_REF] Zhong | Time delay estimation in the presence of clock frequency error[END_REF] but with an additional parameter on the reference signal change. Indeed, the shape of the reference signal is assumed to be not static. Therefore the receivers should be synchronized on the portion of reference signal that provides useful information on the time shift between receivers signal. Additionally, the time delay considered is on the receivers sampling time and not on the (noiseless version of) received signal itself, this is why we can restrict the estimation of ∆ between 0 and 1.

Consequently, the pdfs of both receivers signal can be formulated as:

r i [k] ∼ f 1 (r i [k]) if k ≤ ν f 2 (r i [k]) if k > ν (4)
with ν = τ Te ( . denotes the floor function). In order to simplify notations, we will work with normalized frequency which is equivalent to take the sampling period T e = 1. One can also regroup equations ( 1) and ( 2) into a single formula

r (t) = s (t) + n (t) such that r (k) = r 1 [k] and r (k + ∆) = r 2 [k] with n (t) a strict stationary signal with the same pdf than the pdf of n 1 [k] d = n 2 [k] (" d =" means equality in distribution).
An illustration of the aforementioned problem is shown in figure [START_REF] So | Time-delay estimation for sinusoidal signals[END_REF] where the message m(t) is sinusoidal with a signal change at τ = 0.00015 and a sampling shift of ∆ = 0.5 between both receivers. The goal will be to retrieve parameters τ and ∆ from collected data r 1 [k] and r[k].

Consequently, the received signals are modeled by: Finally, we can concatenate all information into a single observation vector given by:

r 1 [k] ∼ f 1 (r (k)) if k = 1, . . . , ν f 2 (r (k)) if k = ν + 1, . . . , N r 2 [k] ∼ f 1 (r (k + ∆)) if k = 1, . . . , ν f 2 (r (k + ∆)) if k = ν + 1, . . . , N (5)
r = [r 1 [1] r 2 [1] • • • r 1 [N ] r 2 [N ]] T = [r (1) r (1 + ∆) • • • r (N ) r (N + ∆)] T . ( 6
)
We denote θ = [ν ∆] T the vector of unknown parameters.

Contrary to the classical performance bound analysis, the unknown parameter ν only takes integer value and does not parameterize pdfs f 1 and f 2 . This information should be taking into account in order to derive a relevant bounds. This is why in the following Section, we present the Barankin (or McAulay-Seidman) bound and a simplified hybrid version in the aforementioned context. A first general expression is given whatever the pdf of the noise is and we then detail in the Gaussian noise case.

III. PROPOSED BOUNDS A. Relation to prior work

On the one hand, performance bounds have recently been developed for change point estimation [START_REF] Ferrari | Barankin lower bound for change points in independent sequences[END_REF], [START_REF] La Rosa | Barankintype lower bound on multiple change-point estimation[END_REF]. However, this bound did not take into account the estimation of an additional continuous parameter (time delay between receivers). On the other hand, time delay estimation with a change on the reference signal has not been investigated. Therefore, we establish in this paper a lower bound on the MSE in order to predict the accuracy of both parameters: time when signal pattern changes called "change point" and time delay between receivers. Since the time delay is estimated over a compact interval, we treat here an hybrid estimation problem. Consequently, the proposed bound differs from existing lower bounds in the literature.

B. General bound

The Barankin bound is the tightest lower bound on the variance of any unbiased estimator. In this paper, the bound is denoted BB (θ) and is given by [20, Eqn. 2]:

BB (θ) = sup hν ,h∆ [h 1 h 2 ] D -1 (θ, h) [h 1 h 2 ] T (7) with h 1 = [h ν 0]
T and h 2 = h ∆ ] T . The so-called test-points h ν and h ∆ are chosen such that ν + h ν ∈ {2, . . . , T -1} and ∆ + h ∆ ∈ ]0; 1[. The elements of the matrix D (θ, h) are defined as follows:

D i,i (θ, h) = E f (r; θ + h i ) f (r; θ) 2 -1, i = 1, 2 (8)
and

D 1,2 (θ, h) = D 2,1 (θ, h) = E f (r; θ + h 1 ) f (r; θ + h 2 ) f 2 (r; θ) -1. (9)
Since all component of observations vector r are independent, thus its pdf is given by:

f (r; θ) = ν k=1 f 1 (r (k)) f 1 (r (k + ∆)) (10)
× N k=ν+1 f 2 (r (k)) f 2 (r (k + ∆))
We first start by analysing the expression of

D 1,2 (θ, h) = D 2,1 (θ, h). If one assumes h ν > 0 then f (r; θ + h 1 ) f (r; θ + h 2 ) f 2 (r; θ) = ν+hν k=ν+1 f 1 (r (k)) f 1 (r (k + ∆)) f 2 (r (k)) f 2 (r (k + ∆)) × ν k=1 f 1 (r (k + ∆ + h ∆ )) f 1 (r (k + ∆)) N k=ν+1 f 2 (r (k + ∆ + h ∆ )) f 2 (r (k + ∆)) (11)
After some cumbersome calculus, one obtains

D 1,2 (θ, h) + 1 = (12) ν+|hν | k=ν+1 f 1 (r(k + ∆)) f 2 (r(k + ∆)) f 2 (r(k + ∆ + h ∆ ))dr(k + ∆) Lets continue with D 1,1 (θ, h) = E f 2 (r;θ+h 1 ) f 2 (r;θ) -1, it is
almost the previous one with h 2 = h 1 , therefore one has

D 1,1 (θ, h) + 1 = (13) ν+|hν | k=ν+1 f 2 1 (r (k)) f 2 1 (r (k + ∆)) f 2 (r (k)) f 2 (r (k + ∆)) dr (k) dr (k + ∆) Finally, D 2,2 (θ, h) = E f (r;θ+h 2 ) f (r;θ) 2 -1 is obtained similar by previous calculus D 2,2 (θ, h) + 1 = ν k=1 f 2 1 (r (k + ∆ + h ∆ )) f 1 (r (k + ∆)) dr (k + ∆) × N k=ν+1 f 2 2 (r (k + ∆ + h ∆ )) f 2 (r (k + ∆)) dr (k + ∆) . ( 14
)
This expression can be applied to any problem by specifying the probability density functions f 1 and f 2 . We next study the Gaussian case.

C. Gaussian case

Let us assume that n i [k], i = 1, 2 are complex circular Gaussian with zero mean and variance σ 2 . Then the pdf f 1 (r (k)) and f 2 (r (k)) are

f 1 (r(k)) = 1 πσ 2 e -1 σ 2 r(k) 2 , f 2 (r(k)) = f 1 (r(k) -m(k))
By plugging the above expressions in D 1,1 (θ, h), one obtains

D 1,1 (θ, h) + 1 = ν+|hν | k=ν+1 1 π 2 σ 4 e -2 σ 2 ( r(k) 2 + r(k+∆) 2 ) e 1 σ 2 ( r(k)-m(k) 2 + r(k+∆)-m(k+∆) 2 ) dr(k)dr(k + ∆). ( 15
)
By noticing that ∀t

1 πσ 2 e -2 σ 2 r(t) 2 + 1 σ 2 r(t)-m(t) 2 dr (t) = e 2 m(t) 2 σ 2
, one can obtain a simple expression for D 1,1 (θ, h)

D 1,1 (θ, h) = e 2 σ 2 ν+|hν | k=ν+1 ( m(k) 2 + m(k+∆) 2 ) -1. ( 16
)
Concerning D 1,2 (θ, h), by noticing that ∀t

f 1 (r(t)) f 2 (r(t)) f 2 (r(t + h ∆ ))dr(t) = 1 πσ 2 e m(t) 2 -r(t) 2 σ 2 × e -1 σ 2 (Re(r * (t)(m(t)-m(t+h∆)))+ m(t+h∆) 2 ) dr(t). ( 17
)
Let us set

r (t) 2 + 2 Re (r * (t) (m (t) -m (t + h ∆ ))) (18) = r(t) + m(t) -m(t + h ∆ ) 2 -m(t) -m(t + h ∆ ) 2 , then f 1 (r (t)) f 2 (r (t)) f 2 (r (t + h ∆ )) dr (t) = e 1 σ 2 ( m(t)-m(t+h∆) 2 -m(t+h∆) 2 + m(t) 2 ) ,
and finally

D 1,2 (θ, h) + 1 (19) = e 1 σ 2 ν+|hν | k=ν+1   m (k + ∆) -m (k + ∆ + h ∆ ) 2 -m (k + ∆ + h ∆ ) 2 + m (k + ∆) 2   .
The last expression of D 2,2 (θ, h) is obtained by noticing that

f 2 1 (r (t + h ∆ )) f 1 (r (t)) dr (t) = 1,
and that

f 2 2 (r(t + h ∆ )) f 2 (r(t)) dr(t) = 1 πσ 2 e -r(t) 2 -m(t) 2 σ 2 × e -2 σ 2 (Re(r * (t)(m(t)-2m(t+h∆)))+ m(t+h∆) 2 ) dr(t). (20)
By letting

r(t) 2 + 2 Re (r * (t)(m(t) -2m(t + h ∆ ))) (21) = r(t) + m(t) -2m(t + h ∆ ) 2 --2m(t + h ∆ ) 2 , one obtains f 2 2 (r(t + h ∆ )) f 2 (r(t)) dr(t) = e 1 σ 2 ( m(t)-2m(t+h∆) 2 -2 m(t+h∆) 2 + m(t) 2 ) , (22)
and finally

D 2,2 (θ, h) + 1 = e 1 σ 2 N k=ν+1 m(k+∆)-2m(k+∆+h∆) 2 × e 1 σ 2 N k=ν+1 m(k+∆) 2 -2 m(k+∆+h∆) 2 . (23)

D. Simplified bound

Since the synchronization shift ∆ lies in a compact set, the Cramér-Rao bound can be calculated for it, but one has to keep the Barankin bound for the change point location ν. This can be done from the aforementioned proposed bound by noticing that the Barankin bound tends to the Cramér-Rao bound when the test-points tend to zero. Here, one consequently needs to study the behavior of the Barankin bound when h ∆ → 0 [START_REF] Van Trees | Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking[END_REF].

By noticing that around h ∆ → 0,

m (k + ∆ + h ∆ ) = m (k + ∆) + h ∆ ∂m (k + ∆) ∂∆ + o (h ∆ ) (24)
one can obtain (after simple calculus)

ν+|hν | k=ν+1 m(k + ∆) -m(k + ∆ + h ∆ ) 2 (25) -m(k + ∆ + h ∆ ) 2 + m(k + ∆) 2 = ν+|hν | k=ν+1 -2h ∆ Re m(k + ∆) ∂m(k + ∆) ∂∆ + o (h ∆ ) ,
and

N k=ν+1 m(k + ∆) -2m(k + ∆ + h ∆ ) 2 (26)
-2 m(k + ∆ + h ∆ ) 2 + m(k + ∆) 2 = N k=ν+1 2h 2 ∆ ∂m(k + ∆) ∂∆ 2 + o h 2 ∆ .
By using these expressions into D (θ, h), one obtains a "hybrid" Cramér-Rao-Barankin bound, denoted BB (ν, ∆) , which is very simple. This bound will be given in the next Section for a particular function m (t) .

IV. SIMULATIONS

Let's consider two receivers for which their receiving data respectively follow the model given in equations ( 1) and ( 2). This model assumes that receivers are equidistant from the opportunity source. Otherwise, a preprocessing step is necessary to compensate the time delay introduced by the radial distance gap from the opportunity source to each receiver.

The sent message is of sinus form given by m [k] m (kT e ) = e j2πf kTe = e j2πf0k where f 0 = f T e is the normalized frequency. In order to estimate message arrival and the normalized time shift, i.e., θ = [ν ∆]

T , we perform the MLE given by

θ = arg max ν,∆ f (r;ν, ∆)
where r is given in equation ( 6) with

f (r; ν, ∆) = ν k=1 f 1 (r(k))f 1 (r(k + ∆)) N l=ν+1 f 2 (r(l))f 2 (r(l + ∆)).
Additionally, if we assume that both receiver noise follow a i.i.d. zero mean Gaussian distribution i.e.

n 1 [k] = n 2 [k] ∼ N 0, σ 2
, then, after some calculus, the ML criterion is given by the following expression:

θ = arg max ν,∆ ν + N k=ν+1
Re e -j2πf0k r(k) + r(k + ∆)e -j2πf0∆ .

Since parameter ν is discrete and ∆ belongs to [0; By using the formula developed in the Gaussian noise case (see Section III-C), the Barankin bound ( 7) is obtained with the following information matrix

D (θ, h) =   e 4|hν | σ 2 -1 e 4|hν | sin 2 (πf0h ∆ ) σ 2 -1 e 4|hν | sin 2 (πf0h ∆ ) σ 2 -1 e 8 N -ν σ 2 sin 2 (πf0h∆) -1   .
However, the estimation problem is linear w.r.t normalized time ∆, so we can expect no threshold effect for its MLEs MSE. Therefore, we can simplify this bound by mixing with the Cramér-Rao bound for the estimation of ∆. Then, by using the formula developed in Section (III-D), our simplified bound is given by:

BB (ν, ∆) =    sup hν h 2 ν σ e 4|hν | σ 2 -1 0 0 σ 2 8π 2 f 2 0 (N -ν)    . (27)
We compare the MSE of the MLE to the Barankin bound and its simplified form for estimation of ∆ and ν in the following setup: f 0 = 0.1, N = 45, ν = 25 and ∆ = 0.4. The Barankin bound is the maximization of the criterion in [START_REF] Hinkley | Inference about the change-point in a sequence of random variables[END_REF] over a grid search with steps δh ν = 1 and δh ∆ = 10 -3 and the simplified bound is obtained with the same step for ν but does not depend on ∆. The MSE of the MLE is estimated with 1000 Monte-Carlo trials. In the figure (2), the Barankin bound and its simplified form give the same result since the estimation problem w.r.t. ∆ is linear. The performance of MLE estimate for ∆ is well predicted by both bounds. In figure [START_REF] Carter | Coherence and time delay estimation[END_REF], the MSE of the MLE for ν is lower bounded by both Barankin bound and its simplified form. The simplified form gives almost the same performance than the Barankin with a lower computational cost. However, the prediction is not accurate when the SNR becomes higher than 15dB. This difference is probably explained by the fact that ν is a discrete parameter. Indeed, it is proved in [START_REF] Choirat | Estimation in discrete parameter models[END_REF] that there is no asymptotic efficiency of the MLE when estimating discrete parameters from an i.i.d. sequence. Here we are in the more constrained case where the observations are independent but not identically distributed and consequently, there is no reason that the MLE becomes efficient.

V. CONCLUSION

In this paper, we established the Barankin bound and a simplified version of this bound for a synchronization problem where an opportunity source is used. The MLE's accuracy for the normalized time shift is well predicted by our proposed bounds. Therefore, the MLE seems asymptotically efficient for estimation of ∆. However, the prediction of MLEs MSE still need to be improved for parameter ν because of its asymptotic inefficiency for discrete parameter estimation.

Fig. 1 .

 1 Fig. 1. Example of data collected on the receivers from a sinusoidal source with τ = 0.00015 and ∆ = 0.5

Fig. 2 .

 2 Fig. 2. Comparison of MLE MSE, BB and BB for estimation of ∆. The yaxis is the amplitude of the MSE expressed in the square unit of the normalized time.

Fig. 3 .

 3 Fig. 3. Comparison of MLE MSE, BB and BB for estimation of ν. The yaxis is the amplitude of the MSE expressed in the square unit of the normalized time.

  1[, the implementation of the MLE is performed by maximizing the above criterion over a grid search [|1; N |] × [0; 1[ with a step δ ∆ = 10 -3 for parameter ∆.

	10 -1
	10 -2
			BB for ∆
			BB for ∆
		MLE MSE for ∆
	16 10 -3	18	20	22	24	26	28	30
				SNR [dB]


Julien Fayer

Géraldine Morin

Simone Gasparini

Maxime Daisy

Benjamin Coudrin

Radiometric confidence criterion for patch-based inpainting

come

I. INTRODUCTION

Contrary to Augmented Reality, Diminished Reality (DR) [START_REF] Mori | A survey of diminished reality: Techniques for visually concealing, eliminating, and seeing through real objects[END_REF] aims at virtually removing one or more objects from the image of a scene. The image regions containing the object must then be filled seamlessly with a texture coherent with the surrounding regions in order to provide a diminished, photo-realistic view. To this end, completion methods such as inpainting [START_REF] Bertalmio | Image inpainting[END_REF] can be used to fill the region hidden by the object according to the photometric data of their surrounding regions. The results of inpainting approaches fit well for the expectations of Diminished Reality. In order to improve inpainting results in the context of DR, previous works [START_REF] Kawai | Diminished reality based on image inpainting considering background geometry[END_REF], [START_REF] Silberman | Indoor segmentation and support inference from RGBD images[END_REF] exploit the indoor setting: as the background consists of planes (floor, walls, ceiling), instead of applying inpainting techniques directly on the image, they take advantage of the planar geometry of the scene to better guide the inpainting process. For each plane, they compute the homography between the camera plane and the 3D plane. Perspective rectification is applied to the plane region to allow inpainting on distortionfree images. In the undistorted plane, the texture structure recovers its regularity in size and pattern, which is a great advantage for inpainting. However, inpainting methods do not take into account the quality variation in terms of image resolution of these rectified images: this can lead to degraded results, as discussed in Section III. We propose an approach which first computes a trust criterion on the rectified image and uses it to improves the inpainting process in the rectified image.

The paper is organized as follows: §II reviews the stateof-the-art inpainting techniques, §III highlights the issues in inpainted images due to resolution variation and §IV introduces the proposed method to guide and improve the inpainting process. §V shows how the proposed criterion can be integrated in two classical inpainting algorithms, and the corresponding results are given in §VI. Finally, §VII concludes the paper, proposes perspective and future extensions of our approach.

II. RELATED WORKS

Inpainting finds its natural application in image restoration to remove, e.g., scratches and cracks from the frames of old movies or photos, or in special effects to remove entire objects from an image. There are two main approaches for inpainting techniques [START_REF] Bertalmío | Inpainting[END_REF]. Diffusion-based methods [START_REF] Chan | Nontexture inpainting by curvaturedriven diffusions[END_REF] are used in image restoration to fill or correct small regions of the images for which a mask is provided by the user. These methods are generally based on Partial Differential Equations (PDEs) and a diffusion model that iteratively propagate the information from the outside of the mask along the isophotes, i.e. the level lines perpendicular to the gradient vectors of the pixels on the contour. These methods perform well when filling small and smooth regions but are not adapted if a structure or texture needs to be propagated. Moreover, these methods, being iterative, have a high computational cost. Patch-based methods [START_REF] Criminisi | Object removal by exemplarbased inpainting[END_REF], [START_REF] Buyssens | Exemplarbased inpainting: Technical review and new heuristics for better geometric reconstructions[END_REF] are instead used to fill larger portions of the image by copying either single pixels (sparsitybased [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF]), entire patches or a mixture of those from other parts of the image (exemplar-based [START_REF] Criminisi | Object removal by exemplarbased inpainting[END_REF]). For each pixel p of the mask, they search the most similar patch in the image to the one centered in p, and they copy it. The search for this similar patch is the most important but also the most expensive step of the algorithm. Many variants and optimizations have been proposed over the last decade. One of the most effective approaches is PatchMatch [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF], which efficiently finds for every patch the approximate nearest-neighbour in the image using a randomized cooperative hill climbing strategy. In [START_REF] Chican | Constrained patchmatch for image completion[END_REF] the search is restricted to the most likely offsets, reducing both the complexity and also enhancing the propagation of the geometric structures of the image. The other critical step in patch-based methods is the selection of p and the order of filling. Onion-peel order fills the missing data starting from the pixels on the border and proceeding layer after layer towards the region's center. This sometimes leads to unexpected results at the center of the region and, in general, structures are not propagated inside the region. Data-aware methods, instead, give priority to pixels lying on borders of objects, thus favouring the preservation of structures. On the other hand, a known issue of the PatchMatch approach is that it cannot handle well regular textures, i.e., textures embedding regular patterns or structures. Other methods have been proposed to handle regular textures by performing a statistical analysis of the texture that allows to find the map of the dominant directions (or offsets) [START_REF] He | Statistics of patch offsets for image completion[END_REF], [START_REF] Köppel | Fast image completion method using patch offset statistics[END_REF]: the inpainting problem is then cast as a global minimization of an energy function written in terms of the offset map that enforces the structure and texture consistency. In [START_REF] Liu | Exemplar-Based Image Inpainting Using Multiscale Graph Cuts[END_REF], the minimization problem is solved via graph cuts in order to reduce the computational complexity. The latest works, however, uses information such as the camera position and the 3D structure of the scene to improve the final rendering. In the context of DR applications [START_REF] Siltanen | Diminished reality for augmented reality interior design[END_REF], [START_REF] Kawai | Diminished reality based on image inpainting considering background geometry[END_REF], [START_REF] Pavic | Interactive image completion with perspective correction[END_REF], the objects to be removed usually appear lying on a plane or occluding other planes in the background (eg. floor or walls). In this case, the inpainting is not applied directly to the image. First the perspective distortion is removed and the original structure of the texture is recovered by computing the rectifying homography. An inpainting algorithm is then applied to the undistorted view, and the completed image is re-projected into the original image through the inverse homography. However, if the area to be filled belongs to a plane having a skew angle w.r.t. camera optical axis (e.g., the floor of a room), the photometric information of the undistorted image has usually a large variation of resolution. Applying inpainting independently of this resolution variation generates a blur in the reconstructed data, including in the final rendering after reprojection as explained in the next section. Eisenacher et al. [START_REF] Eisenacher | Texture synthesis from photographs[END_REF] proposed a user-driven texture synthesis algorithm: the user can describe the local geometry supporting the texture which can be then re-synthesized it on new surfaces. Although it is not strictly speaking an inpainting method, this is, to the best of our knowledge, the only work that, like ours, tries to took into account the different resolution of the image according to the 3D scene geometry.

III. PROBLEM STATEMENT

In our indoor DR application, we apply inpainting techniques to the perspectively undistorted view of a plane. The undistorted view has a non-uniform quality of data and possible interpolation artefacts, which affects the final inpainting result. We designed a simple experiment to illustrate such effect. Consider a textured plane Π seen by a camera C under an angle γ w.r.t. the plane normal n (see Fig. 1): we rendered two synthetic images, each with a different value of γ (60 • and 85 • , respectively). Fig. 2 shows the two images of the textured plane Π (column (a)), along with the outline of the region to inpaint, and the undistorted views of Π (col. (b)). Once PatchMatch is applied to the undistorted view (col. (c)) and the result is projected back on the original image (col. (d)), we can clearly see that low-resolution data corresponding to the most perspectively distorted part of the image has been used to fill the mask region: the result is visually unsatisfying with a noticeable blur effect in the inpainted region. This simple example shows the influence of the perspective distortion on the result of a classic inpainting technique. In this work we propose to introduce a novel criterion in order to better guide the completion process and thus avoid the propagation of lowresolution data.

IV. RADIOMETRIC CONFIDENCE CRITERION

We first seek a criterion that associates to each image pixel a score characterizing the quality of its projection when it is mapped by the rectifying homography. In particular, our score relies on Bouguer's law, which measures the light flux dF of a light point source I over a surface dS seen under a solid angle dΩ (see Fig. 3):

dF = I dΩ = I cos θ r 2 dS. (1)
Thus the light flux dF over the surface unit dS is

dF dS = I cos θ r 2 . (2)
We deduce, for each point P of the image plane, a trust score as the ratio between the point light flux dF over an infinitesimal area around P on the image plane and the flux around its projected point Q on the plane Π. Using (2), we compute at the points P and Q :

dF dS (P ) = I cos φ r 2 P = I cos 3 φ f 2 (3)
and

dF dS (Q) = I cos θ r 2 Q = I cos 3 θ d 2 , ( 4
)
where f is the camera focal length, and d is the distance between C and Π. Note that r P = f cos φ and r Q = d cos θ . We define our criterion as the ratio between (3) and ( 4)

trust(P ) = dF dS (Q) dF dS (P ) = f d 2 cos θ cos φ 3 . (5)
whose value decreases as the distance d increases and for increasing angles θ. Higher values of trust are pixels of good quality, whereas lower values denote points that are likely to represent regions of the plane far from the camera and/or seen under a very skew angle. We consider P as the center of the pixel p of the image and we define the confidence map C for each pixel of the rectified image domain U × V ⊂ N 2 : with C(p) = 0 if p is not a projection of a point of the plane Π. Fig. 4 shows an example of the map computed for a synthetic image. In the next section we show how we integrate our criterion into two classic inpainting algorithms. Note that our criterion is not restricted to these two methods, but it can be easily integrated for filtering (as shown in §VI) or in other inpainting approaches.

C : U × V -→ R p -→ C(p) = trust(P ) (6)

V. RESOLUTION-AWARE INPAINTING

We integrate the proposed criterion in two state-of-the-art inpainting methods, PatchMatch and statistic-Graphcut. First, let us define a predicate we use in both approaches to handle the selection of a new candidate through the map C and avoid 4, we used a colormap to identify each pixel position and we show inside the mask the final NNF map, i.e. each color inside the mask indicate from which zone of I the pixel has been taken. On the left, the original PatchMatch NNF contains many reddish pixels coming from zones of I with "lower quality data". On the right, our computed NNF uses pixels from "higher data quality" regions of I (greenish pixels are closer to the camera).

the selection of "low resolution" data. The predicate returns false (reject new candidate) or true (accept) and it is defined as:

validation(p, q, C, α) = C(q) ≤ α C(p), (7)
where p is the pixel to complete, q is the candidate proposed by the approach, and α ∈ ]0, 1] is a tolerance parameter.

A. Adapted PatchMatch approach

PatchMatch [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF] efficiently finds matches between patches of an image as detailed in Alg. 1. It defines a nearest neighbour field (NNF) as a function of offsets map : R 2 → R 2 that binds any patch a of an image A to a patch b of an image B according to a distance D (typically, the Sum of Squared Differences (SSD)). The NNF is initialized randomly and updated with a multi-resolution pyramid. Two iterative updates are performed: (a) the offsets with best scores are propagated to adjacent pixels (1.12-15), and (b) a random search in the neighbourhood of all offsets with best distance scores (1.16). for i ← 1 : N do 5:

if i = 1 then 6:
map ← RAND( ) First floor, initialisation of NNF map In PatchMatch the patches are treated identically regardless of their quality, in particular in the low resolution levels of the pyramid. This leads to a spread of low quality data throughout the mask. To counter that, in the propagation step, beside verifying that the candidate patch has a better similarity than the current considered patch, we also use validation to verify that the center pixel q of the candidate patch has a larger confidence than the mask pixel p. Note we are using α = 1 in validation for this case. Our modifications are shown in red in Alg. 1. Fig. 5 visually illustrates the effect of our modifications.

B. Adapted statistic and graphcut approach

As mentioned in Section II, statistic analysis followed by graphcut is more adapted than PatchMatch for inpainting textures with a regular pattern. Alg. 2 details the steps of the considered algorithm [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF], [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. Statistic analysis consists in calculating for each patch of the known zone I |M C its associated offset (2.4). The offset of a patch centered in p is the translation vector t such that:

I(p + t) = arg min u, u >τ sim I(p), I(p + u)
where sim is a similarity function (based on a Sum of Squared Differences) and τ is a threshold to avoid selecting patches which are too close to the considered patch. Then, for each offset, its occurrence is computed and the most expensive K Iinpainted ← APPLYLABELS(labels) 8: end procedure for all p : (x, y) ∈ Ie do In the case of regular textures, it may happen that we do not find an offset that allows to complete a pixel because the only offsets considered point to pixels in the mask. In this case, we could allow a pixel p with a C(p) value to use an offset pointing to a pixel q with a smaller C(q) value but close to C(p). For that, we reset the tolerance parameter α in [START_REF] Criminisi | Object removal by exemplarbased inpainting[END_REF]: at the beginning of the Graphcut step, α is set to 1. If we lack candidates for finding the right label for a mask pixel, we decrease α by a fixed value in order to loosen the constraint and use data of lower slightly quality.

We show in the next section some results based on our adapted algorithm.

VI. EXPERIMENTAL RESULTS

We implemented the proposed method in C++, using CImg [START_REF] Tschumperlé | The CImg Library[END_REF] for image manipulation and its PatchMatch implementation on Android devices. For perspective removal we use the interpolation implementation by [START_REF] Thevenaz | Interpolation revisited [medical images application[END_REF], while for pixel labeling in the offsets approach we use the implementation of Graphcut [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF], [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF]. We compute the plane homography by using its corners which are automatically computed in synthesis cases or with ARCore [23].

Fig. 6 shows the results for two examples of undistorted views with synthetic textures, a stochastic and a regular texture. The mask, in red, has been defined by the removal of objects (see the first row of Fig. 6). We use the state of the art PatchMatch algorithm for the first case and the Graphcut algorithm for the second case. The second row of the Fig. 6 shows that the original result is very blurred because low quality data has been propagated. The third row (with the validation criterion) shows the propagation has been handled without degrading the inpainting process. Fig. 7 shows results of the final output of our DR pipeline on a synthetic image (top row) and two real images (second and third row). For the rectified images shown in the last section, we applied a luminosity adjustment using a diffusion based inpainting method in order to create a more photorealistic output. For each case, we have a propagation of the "low" resolution data to the "high" resolution data zone in the classic approaches, whereas our validation criterion added in the inpainting methods handles better the propagation.

Our implementation gives improved results with stochastic textures, even in complex situations. In the regular texture case, we adjust the parameter α as explained in §V. However, our method supports the propagation of the high quality data, possibly creating a resolution discontinuity at the boundaries of the mask. Thus, after the inpainting process, we apply a Gaussian blur whose kernel varies for each pixel p according to their confidence C(p). At each pixel p the kernel width σ is chosen inversely proportional to the square root of C(p).

VII. CONCLUSION AND FUTURE WORK

We proposed a criterion to better guide the selection of patch in inpainting methods that takes advantage of the planar geometry of the scene to better guide the inpainting process and avoid the propagation of low resolution data. We have shown these improvements handle the blur effect without degrading the quality of the final output. We have also used the criterion a posteriori to recover a continuous quality in the inpainted data.

As future work, we plan to extend our approach to the multiview setting. In such scenario the global rectified image of the plane is generated from all the views. The trust map is not any more continuous, as the C(p) values at the boundaries between the projections of views are not continuous. Pixels of the rectified plane may come from multiple original images, thus leading to the problem of properly defining the confidence when multiple sources are available.

Fig. 1 .

 1 Fig. 1. Camera positions w.r.t. the plane Π according to the angle γ.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Perspective influence on inpainting processing on two synthetic images rendered with γ = 60 • (first row) and γ = 85 • (second row). From right to left: (a) input image with the outline of the zone that will be removed and completed using inpainting, (b) rectified image of the floor, (c) rectified image completed using classical inpainting [10], (d) remapped output, (e) our proposed solution.

Fig. 4 .

 4 Fig. 4. Left: input image of a floor plane (green texture). Center: rectified image of the plane. Right: confidence map C for the rectified image, the Jet colormap (high values in red and low values in blue) shows that the value of trust decreases for the pixels that have a lower quality, i.e. those on the top right part of the image.

Fig. 5 .

 5 Fig. 5. The effects of our modifications of PatchMatch. In a rectified image I like the one in Fig.4, we used a colormap to identify each pixel position and we show inside the mask the final NNF map, i.e. each color inside the mask indicate from which zone of I the pixel has been taken. On the left, the original PatchMatch NNF contains many reddish pixels coming from zones of I with "lower quality data". On the right, our computed NNF uses pixels from "higher data quality" regions of I (greenish pixels are closer to the camera).

Algorithm 1 2 : 3 :

 123 PatchMatch algorithm with trust criterion Require: Rectified image I, rectified binary mask M , trust map C Ensure: Rectified inpainted image Iinpainted 1: procedure PATCHMATCH(I, M, C) Apply diffusion-inpainting in I |M I |M is the restriction of I to M Create multi-resolution pyramid for I and M with N floors 4:

Algorithm 2

 2 Statistic+Graphcut algorithm with trust criterionRequire: Rectified image I, rectified binary mask M , trust map C Ensure: Rectified inpainted image Iinpainted 1: procedure STATISTIC+GRAPHCUT(I, M,C)

9 : 10 :

 910 function COMPUTELABELS(I, M, O,C) Ie ← I |M ∪∂M ∂M is the boundary of M 11:

Fig. 6 .

 6 Fig. 6. Results from PatchMatch (left) and Statistic analysis and Graphcut (right) on an empty area coloured in red (first row) without validation criterion (second row) and with validation criterion (third row).

  Graphcut is performed (2.5) to create the labels map by minimizing an energy function in the formE = E data (p, t) + E reg (p, q, t p , t q ),where E data measures the cost of associating the offset t to pixel p, and E reg is the regularization term used to express the visual consistency of neighbour pixel values pointed by the offset map. Like PatchMatch, this approach is multi-resolution and there is no distinction between the "low" and "high" quality data during the energy minimisation. We integrate our criterion into E data : for a pixel p and a offset t, E data (p, t) is set to 0 if p + t does not belong to the mask AND if validation(p, p + t, C, α) = 1, otherwise it is set to +∞. It allows to penalize the couple (p, t) if trust(p + t) < trust(p).

	12:	if p ∈ ∂M then
	13:	labels(p) ← 0	Match adjacent pixels to null vector
	14:	else
	15:	labels(p) ← RAND(1, |O|)	Random initialisation
	16:	end if
	17:	labels ← GRAPHCUT(Ie, O, C, Edata, Ereg)
	18:	end for
	19:	return labels
	20: end function
	21: function EDATA(p, t, I, M,C)
	22:	if p ∈ ∂M OR (p + t) / ∈ M AND validation(p, p + t, C, α) then
	23:	return 0
	24:	end if
	25:	return +∞
	26: end function
	offsets are selected (they should characterize the regularity of
	a structured texture). Finally,

ACKNOWLEDGMENTS

This research was supported by the CIFRE ANRT 2016/0139 and the regional project FEDER-FSE Midi-Pyrénées et Garonne REALISM n.15056690.