▶ examine and compare alternative methods of incorporating competitive interactions in plant and animal breeding programs ▶ theory, simulations, and a biological experimentation

Outline

Introduction

Theory

Experimental validation

Discussion Griffing (1967) Group of n unrelated individuals:

P i = P D,i + n-1 j̸ =i P S,j = A D,i + E D,i + n-1 j̸ =i
(A S,j + E S,j ) ▶ P i : observed phenotypic value of individual i ▶ A D,i : heritable direct effect (breeding value) of i ▶ E D,i : nonheritable direct effect of i ▶ A S,j : heritable social/indirect effect (breeding value) of j ▶ E S,j : nonheritable indirect effect of j T. Flutre GQE Muir (2005) 16/05/2023 6 / 30 Griffing (1967) Expected change in the mean from selection on individual records:

∆µ = ī σ P (σ 2 A D + (n -1)σ A DS )
where ī is the standardized selection differential (selection intensity)

Competition for a limited resource: Griffing (1967) Expected change in the mean from selection on groups:

▶ σ A,DS < 0 ▶ positive selection
∆µ = 1 n īgr σ Pgr (σ 2 A D + 2(n -1)σ A DS + (n -1) 2 σ 2 A S )
▶ always positive

Needed improvements

Group selection is based entirely on between-group variation and ignores within-group variation. The optimal index would separate the direct and associative effects and weight each according to the variance-covariance structure of the genetic parameters.

▶ ρ D,A = -0.54

▶ σ 2 D /(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.38 ▶ (σ 2 D + σ 2 A )/(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.62 T. Flutre GQE Muir (2005) 16/05/2023 14 / 30 ▶ ρ D,A = -0.53 ▶ σ 2 D /(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.495 ▶ (σ 2 D + σ 2 A )/(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.505

Impact of group size

Simulations with groups of size 16 versus 4 (for the same progeny size: 256).

Correlations between estimated and true genetic effects.

▶ both direct and associative effects were estimated more precisely in small groups than in large, but the difference was small Outline Introduction Theory

Experimental validation

Discussion

Competition for ressources

Restrict access to food:

▶ 15.2 cm per cage ▶ once a day ⇒ 240 g / day ▶ adequate to meet all nutritional requirements provided the birds did not waste feed ▶ in some cages, some birds dominating the social order restricted access to the feeder Parameter estimation

▶ ρ D,A = -0.56 ▶ σ 2 D /(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.21 ▶ (σ 2 D + σ 2 A )/(σ 2 D + σ 2 A + σ 2 ϵ ) = 0.

  can reduce rather than increase the mean

	T. Flutre GQE	Muir (2005)	16/05/2023	6 / 30

  23 BLUP feed conversion efficiency can be improved only if feed efficiency is included in the selection program▶ because of the biology of most traits, it is difficult to envision a situation in which the genetic correlation between direct and associative effects would be positive

	Outline
	▶ even if variation due to associative effects is small relative to
	the environmental variance, the impact of associative effects
	Introduction on group performance can be dramatic, particularly if group
	size is large
	Theory ▶ with D-
	Experimental validation
	Discussion

▶ cases of complementarity, altruism?

Influence of tree i on j: c i→j ∝ 1

▶ phenotypic associative effect of tree i: p,i = a g ,i + a e,i


Franc ¸ois ROUX • Florian Dupriez-Robin
email: florian.dupriez-robin@cea.fr

Franc ¸ois Auger
email: francois.auger@univ-nantes.fr

Real-time simulation of an electric ship in normal and faulty conditions.

Franc ¸ois ROUX • Florian Dupriez-Robin • Guénaël Le Solliec • Franc ¸ois Auger Abstract Today's requirement for improved electrical systems increase needs for accurate, fast and versatile models. MATLAB TM /SIMULINK is a good candidate in terms of average simulation for accuracy, versatility and speed, but for detailed simulation (i.e. designed with the SimScape toolbox) it can be very slow, especially when multiple power converters are simulated. This issue can be solved by real time simulators based on FPGAs like OP-5607 from Opal-RT. With these real-time simulators, complex systems can be simulated with a time step of a few hundreds of nanoseconds. This article presents a new methodology in simulation fields which let researchers to interface fast response power electronics converters and slower physicals parts. This new methodology allows designers to evaluate systems at multiple time scales: minutes, seconds and milliseconds. This article emphasizes today's possibilities of simulating a system as complex as the propulsion line of a cruise ship during an entire 30 min trip. Finally, this kind of simulation allows engineers to increase the accuracy of modeling, even to work on realistic faulty conditions.

Introduction

Today's concern about environmental issues is pushing public institutions to take action about maritime pollution [START_REF] Lin | Compliance with international emission regulations: Reducing the air pollution from merchant vessels[END_REF]. Nowadays, many harbors have stringent environmental regulations (i.e. limitation of CO 2 emission, fine particules matters, noise...) [START_REF] López-Aparicio | Shipping emissions in a nordic port: Assessment of mitigation strategies[END_REF]. All Electrical Ships (AES) are commonly described as a good way to address these issues. Late trends in ship energy management are about Medium Voltage Direct Current (MVDC) AES [START_REF] Mccoy | Trends in ship electric propulsion[END_REF] because of their high integration possibilities of multiples energy source, capabilities of reconfiguration... Last decades evolution on power electronics, especially high voltage/current switches (6.5 kV, 1 kA [START_REF] Abb | HiPak IGBT module[END_REF]) allows designers to reuse regular/well known medium power topologies in high power converters. This drives the MVDC ship [5], which is the best fit for integration of renewable source/storage energy inherently Direct Current (DC). Finally, new high power Fuel Cells (FC), (e.g. Ballard 200 to 1000 kW marine FC), can be easily integrated as the main power supply for a cruise ship.

The choice of modeling type in power electronics is essential, and is driven by the simulation goals. From our point of view, we have 3 main types of simulations :

-Power simulations showing the power exchanges between all subsystems in a complex system without the complexity of control management. It is useful for testing the global behavior, -Average simulation is closer to real systems. The linear behavior of converters can be implemented, stability of regulators can be studied [START_REF] Bacha | Power Electronic Converters Modeling and Control[END_REF]. -Detailed simulation is about switches behavior, it is the most accurate way of studying converters and their interactions (high frequency noise, harmonics...) The main issue is computation time: for a dozen of switches in a simulation, time ratio can reach 1/60 (i.e. for 1 s of simulation, computation time will take 1 min).

The last evolution in field-programmable gate arrays (FP-GAs), especially about simulating power electronics [START_REF] Yamane | Multi-FPGA solution for large power systems and microgrids real time simulation[END_REF] on FPGAs, allows researchers to simulate complex systems with a lot of switches in real-time (i.e. unity time ratio) [START_REF] Genić | Hardware in the loop testing of shipboard power system's management, control and protection[END_REF]. So, in this article we are going to simulate an AES with MVDC bus supplied by a FC and a battery on FPGAs. The following system is applied to a travel with a fault on DC bus. Simulation on FPGAs allows us to study the faulty period and impacts on the overall system. By means of new FP-GAs resolution, this study can be done on multiple times scales: minutes, seconds and milliseconds. The following part of this paper is organized into tree sections: section 2 introduces the methodology of this work and how information is exchanged between each part of the system. Section 3 presents used model and interactions between converters and control. Finally, section 4 shows some results of this work. running together under the SIMULINK environment, electronic part is implemented under OPAL-RT Shematic editor and compiled on FPGAs. Control is intentionally based on very common/well known strategies in order to evaluate the impact of faulty conditions on power electronics without advanced regulation strategies. In order to be closer to the physical behavior, most of the exchanged information between control and physical simulation is done by wired signals. Finally in this paper, we perform a simulation of the propulsion line of an AES, from production to propeller.

Mode

The overall system was applied to a cruise ship 30 min long travel (from Fromentine to île d'Yeu. e.i. from a continental french harbor to a near french island). Eventually FPGAs programming allows us to explore faulty configurations like creating short circuits on the DC bus.

Model and parameters

In this section, we will present as exhaustively as possible the model and the parameters of this work. We will start by the Power Management System (PMS), which manages elements in this simulation, then the electrical model of power converters and their close control. Finally, we will show how the mechanical response of a ship impacts the electrical part. Models and control are always managed by a close safety supervisor. Its first purpose is to manage starting/faulty conditions, high level implementation is presented in Fig. 1. It is in these close supervisors that safety and respect of standards are managed. The main practice that we are following is IEEE Std 1709 TM 2010 about faulty tolerance in MVDC, which mainly consists of 10% tolerance on DC bus and 10 ms to clear faulty condition [5]. The second purpose is the application at converters level of PMS control signals and limitations.

General view of the case study

In ship, symmetry/redundancy is a key for safety, in case of failure of a part of propulsion line the ship can return to the nearest harbor with the remaining power. Because of the that, we are only studying a half ship, the global architecture is presented in Fig. 2.

PMS

The PMS is the supervisor of the entire ship. Its goal is to give orders to every converters supervisors :

start and stop, -power limitation, -high level control (i.g. DC bus level, battery charge, etc), -dispatching between power sources etc. The PMS is based on information provided by converter supervisors, state (starting, steady and fault) and few other information (power, Stat Of Charge (SOC)...). The PMS algorithm is based on the algorithm of supervisors.

FC and DCDC converter

The FC used in this work is based on [START_REF] Njoya | A generic fuel cell model for the simulation of fuel cell vehicles[END_REF] (MATLAB TM Sim-Scape Model). The associated converter is an Interleaved Boost Converter (IBC) (Fig. 2), because of inherent operation of FC, this converter is not reversible. This converter is current controlled by a single Proportional Integral (PI) controller. Design of IBC is based on boost equations [START_REF] Hasaneen | Design and simulation of DC/DC boost converter[END_REF].

Because of the unavailability of data sheet of 800 kW FC we use data sheet of a 50 kW FC. In this purpose, in order to preserve dynamics and behavior, current provided to the FC model is divided by 16.

D = 1 - V in min V out max (1)
L Boost ≥ R ch * D * (1 -D) 2 2 * F Scw (2)
C Boost ≥ V DC * V out * D F Sw * ∆ V out * R ch (3)

Battery and DCDC reversible converter

The battery used in this work is based on [START_REF] Tremblay | Experimental validation of a battery dynamic model for EV applications[END_REF] (MATLAB TM SimScape Model). In this topology, the only converter which can act as an Active Front End (AFE) is the battery DCDC converter. Its main goal is to control the bus voltage. As an AFE, it must be reversible in order exchange energy with the bus. We choose a well known boost medium power topology presented in Fig 2 . This boost is controlled by two cascaded PI controllers, the outer loop controls bus voltage and inner loop controls current [START_REF] Bacha | Power Electronic Converters Modeling and Control[END_REF]. The design of the boost converter is based on Eqs. 1 to 3.

DC bus

A capacitor has been introduced between DCDC converters and the DCAC converter, to ensure a better stability, for starting stage a preload resistor and its switch were added (Fig. 2).

DCAC converter and motor

Finally, a two level DCAC converter (Fig. 2) converts from DC bus to the Alternative Current (AC) Permanent-Magnets Synchronous Motor (PMSM). The control is based on Direct Torque Control (DTC) [START_REF] Bimal | Modern power electronics and ac drives[END_REF]. Pulse Width Modulation (PWM) generation is managed by specialized OPAL-RT digital outputs.

Ship and propeller

The ship model is based on [START_REF] Peng | Robust controller design for marine electric propulsion system over controller area network[END_REF] and experimental data. Table 2 shows power over velocity of the studied ship, total resistance to ship advance is computed by R ship

[N] = P [W ] /v [m/s] .
Eq. 4 shows the result of a second order linear regression between resistance and velocity, which is used as model for the ship resistance.

R ship = 371 V 2 a + 24 V a (4)
The estimation of external characteristics of the propeller is based on OpenProp [START_REF] Epps | OpenProp: An Open-source Design Tool for Propellers and Turbines[END_REF]. Multiple analysis were performed in order to meet ship linear velocity, required thrust and rotational motor speed. Propeller characteristics are presented in Table 3. Torque and thrust of propellers are usually computed from an adimensional advance (5), thrust [START_REF] Bacha | Power Electronic Converters Modeling and Control[END_REF] and torque [START_REF] Yamane | Multi-FPGA solution for large power systems and microgrids real time simulation[END_REF] coefficients [START_REF] Bernitsas | Kt, kq and efficiency curves for the wageningen b-series propellers[END_REF].

J s = V a ω D (5)
T Prop = K T (J s ).ρ ω 2 D 4 (6)
Q Prop = K Q (J s ).ρ ω 2 D 5 (7) J s : advance coefficient Q Prop : torque [Nm] K Q : torque coefficient T Prop : thrust [N] K T : thrust coefficient V a : ship speed [m/s] ρ : water density [kg/m 3 ] ω : rotor speed [rad/s] D : propeller diameter [m]
OpenProp computes thrust and torque coefficients from external characteristic propeller. Table 3 summarizes parameters which feed OpenProp. These K T , K Q versus J s coefficient are feeding the propeller model ((5) to ( 7)). At last, because unsymmetrical behavior of propeller for positive and negative advance coefficients, negative thrust is decreased by 50%. Finally fundamental dynamics relations are applied to compute the ship and rotor speed [START_REF] Genić | Hardware in the loop testing of shipboard power system's management, control and protection[END_REF] to [START_REF] Njoya | A generic fuel cell model for the simulation of fuel cell vehicles[END_REF].

m Va = ∑ f orce = T Prop -R Ship (8) J ωrotor = ∑ torque = Γ rotor -Q Prop
Start V d : V desired V s : V ship P : P ower consign V d -V s V d < -0.2 > 0.5 P = P max P = -P max P = 600 * V 3 d -209 * V 2 d + 316 * V d

Pilot

Pilot determines the level of power applied to the propulsion engine. Basic rules are presented in Fig. 3, in order to achieve the ship travel (speed over time) defined in Fig. 4 (Speed envelope) 5 kn into harbor and channel, 20 kn at sea.

Short circuit

As an illustration of advanced simulation possibilities on FPGAs, we implement a short circuit between DC bus and DCAC converter. To this aim, line impedance has been added at relevant locations (see Fig. 2). To ensure compliance with IEEE standard, a short circuit of a duration < 0.01s was performed. In theses conditions, the system should stay in steady state.

Results

Evaluation of the model on 30 min travel

Fig. 4 shows the modeled speed of the ship during the 30 min case study. On this figure, we can see the evolution of speed (rotor and linear) over speed shape. The simulated ship is within specifications. The third part emphasizes implemented pilot commands (Fig. 3). Fig. 5 shows power source (FC), storage (Battery) and engine. We can see that the battery is charging for the first 20 minutes until its reaches 80% of SOC (limit which has been chosen in the PMS as the end of the force battery recharge). So on FC curve a drop occurs when battery is fully charged. As expected, we can see that the FC supplies an average power to the system (engine power and charging the battery), battery is continuously working as an AFE (i.e. keeps DC voltage bus constant) whether its charging or discharging. On the engine section Fig. 5, the pilot commands are displayed.

Faulty study

Fig. 6 shows the impact of the short circuit on the DC bus. As we can see, the time counter of the boost (used as a watchdog) is counting until its end. If it reaches 10 ms the system switches in safety mode (all switches are forced open) and is waiting for a reset. But in this case, we want to cause a fault that does not prevent the system from continuing to work. On the engine session we can see the reaction of the engine : on magnified views, the engine speed decreases during short circuit, as expected, then the speed increases because of the PI controller response during fault, but everything is back to normal within a second. At time t = 0.2 s, a second fault occurs, because the motor restarts, and thus starts to draw energy from the DC bus again, whereas during the first fault it was supplying energy. This phenomenon is amplified by the saturation state of the controller. Finally we can see one limitation of this set up, engine currents are limited at 1600 A because of acquisition line saturation.

Conclusion

As described in this paper, a new bit-stream approach dedicated to power electronics in FPGAs brings electronics simulation to the next level in terms of accuracy, simulation time and versatility. Engineers with the same model can explore multiple aspects of a system with a high level of fidelity with an improved in simulation time. This article reveals a lack about control strategy : PI control is enough for quasi steady state but not relevant when faulty conditions occur, but it was not the objective. From our point of view, this article is a good baseline, future development could explore different control strategies. The objective was to develop a method which allows engineers to simulate a system with different dynamics (power electronics [µs], engine [ms], hydrodynamic and inertia [s]) with high accuracy and long term impacts assessment. Moreover, with this kind of simulation, research on power, even energy management system can be pursued with detailed model rather than less accurate average models. These kind of results can only be obtained using FPGAs, CPU simulation would be very time consuming. Finally, this method is a step forward of HIL toward Power HIL and rapid control prototyping.

1 Fig. 1

 11 Fig. 1 General view of a supervisor algorithm

Fig. 2

 2 Fig. 2 General view of the studied electrical architecture

( 9 )

 9 m : ship mass [kg] Γ rotor : rotor torque [Nm] J : inertia (rotor + Propeller) [kg.m 2 ]

Fig. 3

 3 Fig. 3 Pilot logigram

Fig. 5

 5 Fig. 4 full travel

Fig. 6

 6 Fig. 6 Faulty and recovery

Table 2

 2 (half) ship power and velocity

	Engine power [kW] Ship velocity [kn] Resistance [kN]
	485	21.2	44
	340.5	18.8	35
	280.5	17.6	31
	243.5	14.1	20
	83	11.8	13
	26.5	8.2	8
	0	0	0
	Table 3 ship power (half ship) and velocity
	Ship	Propeller
	Maximum thrust : 44 kN	Diameter : 2.3 m
	Speed : 20 kn	Number of blades : 3
	Ship mass : 107 t Inertia : 1440 kgm 2	Hub diameter : 0.23 m Rotor speed : 200 rpm
		Propeller type : NACA 65A010

This work was supported by the project "Monitoring and management of marine renewable energies" granted by the French "Pays de la Loire" region.


T Allanche

R Clerc

P Paillet

O Duhamel

V Goiffon

J Rousson

J-P Baudu

J-R Macé

H Desjonqueres

C Monsanglant Louvet

Y Ouerdane

A Boukenter

S Girard

Impact of high ionizing dose on high-power white LEDs

Keywords: LED, radiation, ABC model, camera

Visible cameras operating at 1 MGy of  radiation require radiation tolerant illumination system. Here, commercially available high-power white LEDs have been tested under X and  radiations. The light power of white commercial LED was shown to slightly decrease (less than 10 %) in post irradiation measurements. A closer analysis, based on the ABC model for efficiency droop, reveals that irradiation essentially increases the amount of trap assisted recombination, impacting the LED quantum efficiency mostly at low supply voltage. This suggests that favoring the LED operation at higher supply voltage could be a good strategy to enhance the radiation hardness of LED-based illumination systems.

INTRODUCTION

The CAMRAD project [START_REF] Andra | CAMRAD[END_REF] aims at designing a visible camera able to operate autonomously in harsh (ionizing radiations up to 1 MGy) and eventually in completely dark environment. It is thus necessary to equip such camera with a proper illumination system. In particular, in order to deal with the issue of radiation induced dark current, it needs to operate at low integration times (typically lower than 2 ms).

In this context, LED based lighting systems offer a very good trade off in terms of portability, power consumption and overall lighting power. However, the question of the impact of -radiation on LED remains open. Radiation effects on LEDs have already been studied in the literature, particularly for low-dose space applications [3] [4]. It is admitted that the degradation of LEDs by irradiation is mainly due to atomic displacements caused by particles such as protons or neutrons [5]- [7]. -rays are mainly associated with ionization processes but can possibly cause rare atomic displacements. But little is known about the impact of high level of -rays on LED performance, which is the aim of this paper.

EXPERIMENTAL PROCEDURE

Irradiations conditions

Gamma radiation irradiation was performed with 60 Co at IRMA [START_REF] Irsn | Installation IRMA[END_REF] (IRSN, Saclay, France). X-rays irradiation were performed at the laboratoire Hubert Curien (Saint-Etienne, France) thank to the LabHC MOPERIX machine (100 kV Xray tube with W-target filtered by 4 mm thick of beryllium) that delivers X-rays beam with a maximum at 40 keV. Irradiations were performed at room temperature.

Sample Characterizations

LEDs emission spectra in W/nm were performed using an integrating sphere and a calibrated lamp for several supply current. In order to ensure that the LED thermalized, a delay of 5 min was let between each current step. high energies. The LED spectrum is not modified by the radiation. Table 1 summarizes our results for 5 Commercial Off-The-Shelf (COTS) high-power white LEDs  irradiated at 1 MGy. For all LEDs, the power decrease remains below 10%.

RESULTS

Post-irradiation

Brand

LUMI LED

CREE

Result analysis using the ABC Model

The ABC model is a simple and extensively used model to explain the InGaN LED efficiency droop, by the competition between radiative recombination on one side and Auger and trap assisted recombination on the other side. Following [START_REF] Schubert | Light-emitting diodes[END_REF], in strong injection regime, the external quantum efficiency (EQE) is computed as a function of the diode current density j according to:

𝑬𝑸𝑬 ∝ 𝑩𝒏 𝟐 𝑨𝒏 + 𝑩𝒏 𝟐 + 𝑪𝒏 𝟑 𝒋 ∝ (𝑨𝒏 + 𝑩𝒏 𝟐 + 𝑪𝒏 𝟑 )
where n is the electron/hole concentration in the emission zone, and An the trap assisted, Bn 2 the radiative recombination, Cn 3 the Auger recombination rate. This model has been used to fit the EQE versus current curves of two different devices, for different levels of irradiation (Fig. 3). Interestingly, it was possible to fit all the curves by tuning only the value of the parameter A (which depends on the traps concentration). B and C, which are material dependent parameters, were kept constant, and their value found in good agreement with the literature [START_REF] Lin | [END_REF]. Moreover, Fig. 3 shows that 𝐴 has the same order of magnitude in the 2 LEDs, and tends to increase linearly with irradiation dose. These results suggest that losses in LED are caused by a growing concentration of recombination of traps, bringing insight on the impact of irradiation of the performance of InGaN LEDs.

Moreover, this conclusion has also an impact on the LED selection for radiation tolerant system. Indeed, LED systems are typically composed of several mid-power LED operating at low current, close to the maximum of quantum efficiency, in order to maximize the overall luminous efficiency and to minimize the power consumption. According to previous conclusions, such system is expected to be more impacted by gamma radiation, as an increase of A impacts mostly the low current regime of the EQE versus current curve. An alternative to produce the same amount of luminous flux consists in using a reduced number of power LED operating at high current. If this second solution is usually avoided (it has a higher power consumption), these results suggest that for irradiation purpose, this second option would be less sensitive to damage produced by gamma radiation.

Fig. 1 :

 1 Fig. 1 : Schematic representation of the set up for measuring LED spectral power.

Fig. 2 :

 2 Fig. 2: MKRBWT LED before/after X and  radiation a) Spectral measurements at the nominal current (400 mA) for different irradiation conditions, b) Power vs current before/after radiation.

Fig. 2

 2 Fig.2presents one of the results obtained with one LED type (MKRBWT by CREE) with the comparison between its output powers before and after 1MGy TID. The fact that the LED is biased or not during the exposure does not change the good tolerance of the device. -rays cause more damage than X-rays because Compton electrons are more probable at

Fig. 3 :

 3 Fig. 3 : EQE versus current for pristine and gamma irradiated a) MKRAWT, b) XPEBWT LED device.

Fig. 4 :

 4 Fig.4: Trap assisted parameter A extracted by fitting data of Fig.3versus irradiation dose for both LEDs.

Table 1 :

 1 Power decrease for 5 commercial LEDs at 1 MGy -rays (nominal voltage and current)

	Product	LCH1	MKR	MKR	SDW8	XPEB
	Number	3090	AWT	BWT	4F1C	WT
	1 MGy	-8.2 %	+-1 %	-8.8 %	-4.6 +2.6 %

CONCLUSIONCommercial high-power white LEDs were irradiated with  or X-rays at the MGy level. In both cases, the emitted power decrease was measured to be below 10%. Thanks to the ABC model for efficiency droop, gamma radiation were found responsible for an increase of recombination traps, without impacting the LED radiation and auger recombination properties. In consequence, LEDs operating at high current at the expense of power consumption were found less sensitive to radiation effects.


Kerkyacharian Gérard

Kyriazis George

Pennec Le

Petrushev Erwan

Picard Pencho

Dominique

Inversion of noisy Radon transform by SVD based needlets

A linear method for inverting noisy observations of the Radon transform is developed based on decomposition systems (needlets) with rapidly decaying elements induced by the Radon transform SVD basis. Upper bounds of the risk of the estimator are established in L p (1 ≤ p ≤ ∞) norms for functions with Besov space smoothness. A practical implementation of the method is given and several examples are discussed.

Introduction

Reconstructing images (functions) from their Radon transforms is a fundamental problem in medical imaging and more generally in tomography. The problem is to find an accurate and efficient algorithm for approximation of the function to be recovered from its Radon projections. In this paper, we consider the problem of inverting noisy observations of the Radon transform. As in many other inverse problems, there exists a basis which is fully adapted to the problem, in particular, the inversion in this basis is very stable; this is the Singular Value Decomposition (SVD) basis. The Radon transform SVD basis, however, is not quite suitable for decomposition of functions with regularities in other than L 2 -related spaces. In particular, the SVD basis is not quite capable of representing local features of images, which are especially important to recover.

The problem requires a special construction adapted to the sphere and the Radon SVD, since usual tensorized wavelets will never reflect the manifold structure of the sphere and will necessarily create unwanted artifacts, or will concentrate on special features (such as ridgelets...).

Our idea is to design an estimation method for inverting the Radon transform which has the advantages of maximum localization of wavelet based methods combined with the stability and computability of the SVD methods. To this end we utilize the construction from [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF] (see also [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]) of localized frames based on orthogonal polynomials on the ball, which are closely related to the Radon transform SVD basis. As shown in the simulation section the results obtained are quite promising.

To investigate the properties of this method, we perform two different studies. The first study is of theoretical kind and investigates the possible losses (in expectation) of the method in the 'minimax framework'. This principle, fairly standard in statistics, consists in analyzing the mathematical properties of estimation algorithms via optimization of their worst case performances over large ensembles of parameters. We carry out this study in a random model which is also well known in statistics, the white noise model. This random model is a toy model well admitted in statistics since the 80's as an approximation of the 'real' model on scattered data. It is proved, for instance in [START_REF] Brown | Asymptotic equivalence of nonparametric regression and white noise[END_REF] that the regression model with uniform design and the white noise model are close in the sense of Le Cam's deficiency -which roughly means that any procedure can be transferred from one model to the other, with the same order of risk-. This model has the main advantage of avoiding unnecessary technicalities. In this context we prove that over large classes of functions (described later), our method has optimal rates of convergence, for all the L p losses. To our knowledge, most of the parallel results are generally stated for L 2 losses, as in [START_REF] Madych | Polynomial based algorithms for computed tomography[END_REF] for example, very few (if any) consider L p losses while it is a warrant for instance that the procedure will be able to detect small features. Again, the problem of choosing appropriated spaces of regularity in this context is a serious question, and it is important to consider the spaces which may be the closest to our natural intuition: those which generalize to the present case the approximation properties shared by standard Besov and Sobolev spaces. We can also prove that our results apply for ordinary Besov spaces.

In the case p ≥ 4 we exhibit here new minimax rates of convergence, related to the ill posedness coefficient of the inverse problem d- 1 2 along with edge effects induced by the geometry of the ball. These rates are interesting from a statistical point of view and have to be compared with similar phenomena occurring in other inverse problems involving Jacobi polynomials (e.g. Wicksell problem), see [START_REF] Kerkyacharian | Needlet algorithms for estimation in inverse problems[END_REF].

Our second study of the performances of our procedure is performed on simulations. Since in practical situations scattered data are generally observed, we carried out our simulation study in the scattered data model. We basically compared our method to the SVD procedure -since it is the most commonly studied method in statistics-and the simulation study consistently predicts quite good performances of our procedure and a comparison extensively in favor of our algorithm. One could object that it is a rather common opinion that 'one should smooth the SVD'. However, there are many ways to do so (for instance, we men-tion a parallel method, employing a similar idea for smoothing out the projection operator but without using the needlet construction and in a no-noise framework, which has been developed by Yuan Xu and his co-authors in [START_REF] Xu | Reconstruction from radon projections and orthogonal expansion on ball[END_REF][START_REF] Xu | Fast oped algorithm for reconstruction of images from radon data[END_REF][START_REF] Xu | Image reconstruction by OPED algorithm with averaging[END_REF]). Ours has the advantage of being optimal for at least one point of view since we are able to obtain the right rates of convergence in L p norms.

The paper is structured as follows. In Section 2 we introduce the model and the Radon transform Singular Value Decomposition. In Section 3 we give the class of linear estimators built upon the SVD. We also give the needlet construction and introduce the needlet estimation algorithm. In Section 4 we establish bounds for the risk of this estimate over large classes of regularity spaces. Section 5 is devoted to the practical implementation and results of our method. Section 6 is an appendix where the proofs of some claims from Section 3 are given.

2 Radon transform and white noise model

Radon transform

Here we recall the definition and some basic facts about the Radon transform (cf. [START_REF] Helgason | The Radon transform[END_REF], [START_REF] Natterer | The mathematics of computerized tomography[END_REF], [START_REF] Logan | Optimal reconstruction of a function from its projections[END_REF]). Denote by

B d the unit ball in R d , i.e. B d = {x = (x 1 , . . . , x d ) ∈ R d : |x| ≤ 1} with |x| = ( d i=1 x 2 i ) 1/2
and by S d-1 the unit sphere in R d . The Lebesgue measure on B d will be denoted by dx and the usual surface measure on S d-1 by dσ(x) (sometimes we will also deal with the surface measure on S d which will be denoted by dσ d ). We let |A| denote the measure

|A| = A dx if A ⊂ B d as well as |A| = A dσ(x) if A ⊂ S d-1 .
The Radon transform of a function f is defined by

Rf (θ, s) = y∈θ ⊥ sθ+y∈B d f (sθ + y)dy, θ ∈ S d-1 , s ∈ [-1, 1],
where dy is the Lebesgue measure of dimension d -1 and θ ⊥ = {x ∈ R d : x, θ = 0}. With a slight abuse of notation, we will rewrite this integral as

Rf (θ, s) = y,θ =s f (y)dy.
It is easy to see (cf. e.g. [START_REF] Natterer | The mathematics of computerized tomography[END_REF]) that the Radon transform is a bounded linear operator mapping

L 2 (B d , dx) into L 2 S d-1 × [-1, 1], dµ(θ, s) , where dµ(θ, s) = dσ(θ) ds (1 -s 2 ) (d-1)/2 .

Noisy observation of the Radon transform

We consider observations of the form

dY (θ, s) = Rf (θ, s)dµ(θ, s) + εdW (θ, s), (2.1)
where the unknown function f belongs to L 2 (B d , dx). The meaning of this equation is that for any φ(θ, s)

in L 2 (S d-1 × [-1, 1], dµ(θ, s)) one can observe Y φ = φ(θ, s)dY (θ, s) = S d-1 ×[-1,1] Rf (θ, s)φ(θ, s)dµ(θ, s) + ε φ(θ, s)dW (θ, s) = Rf, φ µ + εW φ .
Here W φ = φ(θ, s)dW (θ, s) is a Gaussian field of zero mean and covariance

E(W φ , W ψ ) = S d-1 ×[-1,1] φ(θ, s)ψ(θ, s)dσ(θ) ds (1 -s 2 ) (d-1)/2 = φ, ψ µ .
The goal is to recover the unknown function f from the observation of Y . As explained in the introduction, this model is a toy model, fairly accepted in statistics as an approximation of the 'real model' of scattered data. The study is carried out in this setting to avoid unnecessary technicalities.

Our idea is to devise an estimation scheme which combines the stability and computability of SVD decompositions with the superb localization and multiscale structure of wavelets. To this end we utilize a frame (essentially following the construction from [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]) with elements of nearly exponential localization which is compatible with the SVD basis of the Radon transform. This procedure is also to be considered as a first step towards a nonlinear procedure especially suitable to handle spatial adaptivity since real objects frequently exhibit a variety of shapes and spatial inhomogeneity.

Polynomials and Singular Value Decomposition of the Radon transform

The SVD of the Radon transform was first established in [START_REF] Cormack | Representation of a function by its line integrals with some radiological applications II[END_REF][START_REF] Davison | A singular value decomposition for the radon transform in ndimensional eucledian space[END_REF][START_REF] Louis | Orthogonal function series expansions and the null space of the radon transform[END_REF]. In this regard we also refer the reader to [START_REF] Natterer | The mathematics of computerized tomography[END_REF][START_REF] Xu | Reconstruction from radon projections and orthogonal expansion on ball[END_REF]. In this section we record some basic facts related to the Radon SVD and recall some standard definitions which will be used in the sequel.

Jacobi and Gegenbauer polynomials

The Radon SVD bases are defined in terms of Jacobi and Gegenbauer polynomials. The Jacobi polynomials P (α,β) n

, n ≥ 0, constitute an orthogonal basis for the space

L 2 ([-1, 1], w α,β (t)dt) with weight w α,β (t) = (1 -t) α (1 + t) β , α, β > -1.
They are standardly normalized by P (α,β) n

(1) = n+α n and then [START_REF] Andew | Special Functions[END_REF][START_REF] Erdélyi | Higher Transcendental Functions[END_REF][START_REF] Szegő | Orthogonal polynomials[END_REF] 1

-1 P (α,β) n (t)P (α,β) m (t)w α,β (t)dt = δ n,m h (α,β) n , where
h (α,β) n = 2 α+β+1 (2n + α + β + 1) Γ(n + α + 1)Γ(n + β + 1) Γ(n + 1)Γ(n + α + β + 1) . (2.
2)

The Gegenbauer polynomials C λ n are a particular case of Jacobi polynomials, traditionally defined by

C λ n (t) = (2λ) n (λ + 1/2) n P (λ-1/2, λ-1/2) n (t), λ > -1/2,
where by definition (a) n = a(a + 1) . . .

(a + n -1) = Γ(a+n) Γ(a) . It is readily seen that C λ n (1) = n+2λ-1 n = Γ(n+2λ) n!Γ(2λ) and 1 -1 C λ n (t)C λ m (t)(1 -t 2 ) λ-1 2 dt = δ n,m h (λ) n with h (λ) n = 2 1-2λ π Γ(λ) 2 Γ(n + 2λ) (n + λ)Γ(n + 1)
.

(2.3)

Polynomials on B d and S d-1

We detail the following well known notations which will be used in the sequel. Let Π n (R d ) be the space of all polynomials in d variables of degree ≤ n. We denote by P n (R d ) the space of all homogeneous polynomials of degree n and by V n (R d ) the space of all polynomials of degree n which are orthogonal to lower degree polynomials with respect to the Lebesgue measure on B d . V 0 is the set of constants. We have the following orthogonal decomposition:

Π n (R d ) = n k=0 V k (R d ).
Also, denote by H n (R d ) the subspace of all harmonic homogeneous polynomials of degree n and by H n (S d-1 ) the restriction of the polynomials from H n (R d ) to S d-1 . Let Π n (S d-1 ) be the space of restrictions to S d-1 of polynomials of degree ≤ n on R d . As is well known

Π n (S d-1 ) = n m=0 H m (S d-1 )
(the orthogonality is with respect of the surface measure dσ on S d-1 ).

Let Y l,i , 1 ≤ i ≤ N d-1 (l), be an orthonormal basis of H l (S d-1 ), i.e.

S d-1 Y l,i (ξ)Y l,i ′ (ξ)dσ(ξ) = δ i,i ′ .
Then the natural extensions of Y l,i on B d are defined by Y l,i (x) = |x| l Y l,i x |x| and satisfy

B d Y l,i (x)Y l,i ′ (x)dx = 1 0 r d-1 S d-1 Y l,i (rξ)Y l,i ′ (rξ)dσ(ξ)dr = 1 0 r d+2l-1 S d-1 Y l,i (ξ)Y l,i ′ (ξ)dσ(ξ)dr = δ i,i ′ 1 2l + d .
For more details we refer the reader to [START_REF] Dunkl | Orthogonal polynomials of several variables[END_REF].

The spherical harmonics on S d-1 and orthogonal polynomials on B d are naturally related to Gegenbauer polynomials. The kernel of the orthogonal projector onto H n (S d-1 ) can be written as (see e.g. [START_REF] Stein | Introduction to Fourier Analysis on Euclidian spaces[END_REF]) if

N d-1 (n) is the dimension of H n (S d-1 ): N d-1 (n) i=1 Y l,i (ξ)Y l,i (θ) = 2n + d -2 (d -2)|S d-1 | C d-2 2 n ( ξ, θ ). (2.4)
The "ridge" Gegenbauer polynomials

C d/2 n ( x, ξ ) are orthogonal to Π n-1 (B d ) in L 2 (B d
) and the kernel L n (x, y) of the orthogonal projector onto V n (B d ) can be written in the form (see e.g. [START_REF] Petrushev | Approximation by ridge functions and neural networks[END_REF][START_REF] Xu | Reconstruction from radon projections and orthogonal expansion on ball[END_REF])

L n (x, y) = 2n + d |S d-1 | 2 S d-1 C d/2 n ( x, ξ )C d/2 n ( y, ξ )dσ(ξ) (2.5) = (n + 1) d-1 2 d π d-1 S d-1 C d/2 n ( x, ξ )C d/2 n ( y, ξ ) C d/2 n 2 dσ(ξ).
The following important identities are valid for "ridge" Gegenbauer polynomials:

B d C d/2 n ( ξ, x )C d/2 n ( η, x )dx = h (d/2) n C d/2 n (1)
C d/2 n ( ξ, η ), ξ, η ∈ S d-1 , (2.6)
and, for

x ∈ B d , η ∈ S d-1 , S d-1 C d/2 n ( ξ, x )C d/2 n ( ξ, η )dσ(ξ) = |S d-1 |C d/2 n ( η, x ), (2.7)
see e.g. [START_REF] Petrushev | Approximation by ridge functions and neural networks[END_REF]. By (2.5) and (2.7) L n (x, ξ)L n (y, ξ)dσ(ξ) = (2n + d)L n (x, y).

L n (x, ξ) = (2n + d) |S d-1 | C d/2 n ( x, ξ ), ξ ∈ S d-1 , f 3,1,0 f 4,4,1 f 14,4,1 f 14,12,0

The SVD of the Radon transform

Assume that {Y l,i : 1 ≤ i ≤ N d-1 (l)} is an orthonormal basis for H l (S d-1 ). Then it is standard and easy to see that the family of polynomials

f k,l,i (x) = (2k+d) 1/2 P (0, l+d/2-1) j (2|x| 2 -1)Y l,i (x), 0 ≤ l ≤ k, k-l = 2j, 1 ≤ i ≤ N d-1 (l),
form an orthonormal basis of V k (B d ), see e.g. [START_REF] Dunkl | Orthogonal polynomials of several variables[END_REF]. On the other hand the collection

g k,l,i (θ, s) = [h (d/2) k ] -1/2 (1-s 2 ) (d-1)/2 C d/2 k (s)Y l,i (θ), k ≥ 0, l ≥ 0, 1 ≤ i ≤ N d-1 (l), is obviously an orthonormal basis of L 2 (S d-1 × [-1, 1], dµ(θ, s)).
2.1. For any f ∈ L 2 (B d ) Rf = k≥0 λ k 0≤l≤k, k-l≡0 (mod 2) 1≤i≤N d-1 (l) f, f k,l,i g k,l,i (2.8)
and for any g

∈ L 2 (S d-1 × [-1, 1], dµ(θ, s)) R * g = k≥0 λ k 0≤l≤k, k-l≡0 (mod 2) 1≤i≤N d-1 (l)
g, g k,l,i µ f k,l,i .

(2.9)

Furthermore, for f ∈ L 2 (B d ) f = k≥0 λ -1 k 0≤l≤k, k-l≡0 (mod 2) 1≤i≤N d-1 (l)
Rf, g k,l,i µ f k,l,i .

(2.10)

In the above identities the convergence is in L 2 and

λ 2 k = 2 d π d-1 (k + 1)(k + 2) . . . (k + d -1) = 2 d π d-1 (k + 1) d-1 ∼ k -d+1 .
(2.11)

Remark : Observe that if k ≥ 0, 0 ≤ l ≤ k, k -l ≡ 0 (mod 2)
, and

1 ≤ i ≤ N d-1 (l), then R * f k,l,i = 0. ⋆
For the proof of this result, we refer the reader to [START_REF] Natterer | The mathematics of computerized tomography[END_REF] and [START_REF] Xu | Reconstruction from radon projections and orthogonal expansion on ball[END_REF]. We will only show in the appendix that λ 2 k has the value given in (2.11).

3 Linear estimators built upon the SVD

The general idea

In a general noisy inverse model

dY t = Kf dt + εdW t ,
where K is a linear operator mapping f ∈ H → Kf ∈ K, and H and K are two Hilbert spaces, the SVD yields a family of linear estimators via the following classical scheme. Suppose K has an SVD

Kf = m σ m f, e m e * m , f ∈ H,
where {e m } and {e * m } are orthonormal bases for H and K, respectively, and Ke m = σ m e * m and K * e * m = σ m e m with K * being the adjoint operator of K. We also assume that

σ m → 0. Then, if σ m = 0, e * m dY t = Kf • e * m dt + ε e * m dW t = f • K * e * m dt + ε e * m dW t = σ m f e m dt + ε e * m dW t and hence 1 σ m e * m dY t = f, e m + ε σ m e * m dW t . (3.1)
In going further, suppose that {φ l } is a tight frame for H. Therefore, for any

f ∈ H f = l α l φ l , α l = f, φ l .
We can represent φ l in the basis {e m }:

φ l = m φ l , e m e m = m γ m l e m
and hence

α l = m γ m l f, e m .
On account of (3.1) this leads to the estimator

fN = l≤N αl φ l with αl = m γ m l 1 σ m e * m dY t , (3.2)
where N is a parameter. By (3.1) we have

αl = m γ m l f, e m + m γ m l ε σ m e * m dW t = α l + Z l ,
where

Z l has a normal distribution N(0, m (γ m l ) 2 ε 2 σ 2 m
). In this scheme the factors

1 σ 2 m
, which are inherent to the inverse model, bring instability by inflating the variance.

The selection of the frame {φ l } is critical for the method described above. The standard SVD method corresponds to the choice φ l = e l . This SVD method is very attractive theoretically and can be shown to be asymptotically optimal in many situations (see Dicken and Maass [START_REF] Dicken | Wavelet-Galerkin methods for ill-posed problems[END_REF], Mathé and Pereverzev [START_REF] Mathé | Geometry of linear ill-posed problems in variable Hilbert scales[END_REF] together with their nonlinear counterparts Cavalier and Tsybakov [START_REF] Cavalier | Sharp adaptation for inverse problems with random noise[END_REF], Cavalier et al [START_REF] Cavalier | Oracle inequalities for inverse problems[END_REF], Tsybakov [START_REF] Tsybakov | On the best rate of adaptive estimation in some inverse problems[END_REF], Goldenschluger and Pereverzev [START_REF] Goldenshluger | On adaptive inverse estimation of linear functionals in Hilbert scales[END_REF], Efromovich and Kolchinskii [START_REF] Efromovich | On inverse problems with unknown operators[END_REF]). It also has the big advantage of performing a quick and stable inversion of the operator K. However, while the SVD bases are fully adapted to describe the operator K, they are usually not quite appropriate for accurate description of the solution of the problem with a small number of parameters. Although the SVD method is suitable for estimating the unknown function f with an L 2 -loss, it is also rather inappropriate for other losses. It is also restricted to functions which are well represented in terms of the Sobolev space associated to the SVD basis.

Switching to an arbitrary frame {φ m }, however, may yield additional instability through the factors (γ m l ) 2 's. Our idea is to utilize a frame {φ l } which is compatible with the SVD basis {e m }, allowing to keep the variance within reasonable bounds, and has elements with superb space localization and smoothness, guaranteeing excellent approximation of the unknown function f . In the following we implement the above described method to the inversion of the Radon transform. We shall build upon the frames constructed in [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF] and called "needlets".

Construction of needlets on the ball

In this part we construct the building blocks of our estimator. We will essentially follow the construction from [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF].

The orthogonal projector

L k on V k (B d ) . Let {f k,l,i } be the orthonormal basis of V k (B d ) defined in §2.3.3. Denote by T k the index set of this basis, i.e. T k = {(l, i) : 0 ≤ l ≤ k, l ≡ k(mod 2), 0 ≤ i ≤ N d-1 (l)}. Also, set ν = d/2 -1. Then the orthogonal projector of L 2 (B d ) onto V k (B d ) can be written in the form L k f = B d f (y)L k (x, y)dy with L k (x, y) = l,i∈T k f k,l,i (x)f k,l,i (y).
Using (2.4) L k (x, y) can be written in the form

L k (x, y) (3.3) = (2k + d) l≤k, k-l≡0(2) P (0,l+ν) j (2|x| 2 -1)|x| l P (0,l+ν) j (2|y| 2 -1)|y| l i Y l,i x |x| Y l,i y |y| = (2k + d) |S d-1 | l≤k, k-l≡0(2) P (0,l+ν) j (2|x| 2 -1)|x| l P (0,l+ν) j (2|y| 2 -1)|y| l 1 + l ν C ν l x |x| , y |y| .
Another representation of L k (x, y) has already be given in (2.5). Clearly

B d L m (x, z)L k (z, y)dz = δ m,k L m (x, y) (3.4) and for f ∈ L 2 (B d ) f = k≥0 L k f and f 2 2 = k L k f 2 2 = k L k f, f . (3.5) 3.2.2 Smoothing Let a ∈ C ∞ [0, ∞) be a cut-off function such that 0 ≤ a ≤ 1, a(t) = 1 for t ∈ [0, 1/2]
and supp a ⊂ [0, 1]. We next use this function to introduce a sequence of operators on L 2 (B d ). For j ≥ 0 write

A j f (x) = k≥0 a k 2 j L k f (x) = B d A j (x, y)f (y)dy with A j (x, y) = k a k 2 j L k (x, y).
Also, we define

B j f = A j+1 f -A j f . Then setting b(t) = a(t/2) -a(t) we have B j f (x) = k b k 2 j L k f (x) = B d B j (x, y)f (y)dy with B j (x, y) = k b k 2 j L k (x, y).
Evidently, for f ∈ L 2 (B d )

A j f, f = k a k 2 j L k f, f ≤ f 2 2 (3.6) and lim j→∞ A j f -f 2 = lim j→∞ (A 0 + j-1 m=0 B m )f -f 2 = 0. (3.7)
An important result from [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF] (see also [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]) asserts that the kernels A j (x, y), B j (x, y) have nearly exponential localization, namely, for any M > 0 there exists a constant c M > 0 such that

|A j (x, y)|, |B j (x, y)| ≤ C M 2 jd (1 + 2 j d(x, y)) M W j (x) W j (y) , x, y ∈ B d , (3.8)
where

W j (x) = 2 -j + 1 -|x| 2 , |x| 2 = |x| 2 d = d i=1 x 2 i , and
d(x, y) = Arccos ( x, y + 1 -|x| 2 1 -|y| 2 ), x, y = d i=1 x i y i . (3.9)
The left part (1) of Figure 2 illustrates this concentration: it displays the influence of a point x to the value of B j f at a second point y 0 , namely the values of B j (x, y 0 ) for a fixed y 0 and j = 4. This influence peaks at y 0 and vanishes exponentially fast to 0 as soon as one goes away from y 0 . The central part (2) of Figure 2 shows the modification of the concentration when j is set to a large value (j = 6). The right part (3) of Figure 2 shows the lack of concentration of B j when the cut-off function a used is far from being C ∞ . The resulting kernel still peaks at y 0 but the value of B j f at y 0 is strongly influenced by values far away from y 0 . Remark : At this point it is important to notice the following correspondence which will be used in the sequel. For

B j (x, y 0 ) j 4 6 4 a C ∞ C ∞ non smooth (1) (2) (3)
S d + = {(x, z) ∈ R d × R + , |x| 2 d + z 2 = 1}, we have the natural bijection x ∈ B d → x = (x, 1 -|x| 2 ) and d(x, y) = d S d + (x, ỹ), where d S d
+ is the geodesic distance on S d + . ⋆

Approximation

Here we discuss the approximation properties of the operators {A j }. We will show that in a sense they are operators of "near best" polynomial L p -approximation.

Denote by

E n (f, p) the best L p -approximation of f ∈ L p (B d ) from Π n , i.e. E n (f, p) = inf P ∈Πn f -P p . ( 3
B d |A j (x, y)|dy ≤ c * , x ∈ B d , j ≥ 0,
where c * is a constant depending only on d. Therefore, the operators A j are (uniformly) bounded on L 1 (B d ) and L ∞ (B d ), and hence, by interpolation, on

L p (B d ), 1 ≤ p ≤ ∞, i.e. A j f p ≤ c * f p , f ∈ L p (B d ). (3.11)
On the other hand, since a(t) = 1 on [0, 1/2] we have A j P = P for P ∈ Π 2 j-1 .

We use this and (3.11) to obtain, for f ∈ L p (B d ) and an arbitrary polynomial

P ∈ Π 2 j-1 , f -A j f p = f -P +P -A j f p ≤ f -P p + A j (P -f ) p ≤ (1+c * ) f -P p = K f -P p .
Consequently, f -A j f p ≤ KE 2 j-1 (f, p). In the opposite direction, evidently,

A j f ∈ Π 2 j and hence E 2 j (f, p) ≤ f -A j f p . Therefore, for f ∈ L p (B d ), 1 ≤ p ≤ ∞, E 2 j (f, p) ≤ f -A j f p ≤ KE 2 j-1 (f, p).
(3.12)

These estimates do not tell the whole truth about the approximation power of A j .

It is rather obvious that because of the superb localization of the kernel A j (x, y) the operator A j provides far better rates of approximation than E 2 j-1 (f, ∞) away from the singularities of f . In contrast, the kernel S j (x, y) = 0≤k≤2 j L k (x, y) of the orthogonal projector S j onto Π 2 j is poorly localized and hence S j is useless for approximation in L p , p = 2. This partially explains the fact that the traditional SVD estimators perform poorly in L p -norms when p = 2.

Splitting procedure

Let us define

C j (x, y) = m a m 2 j L m (x, z)) and D j (x, y) = m b m 2 j L m (x, z).
Note that C j and D j have the same localization as the localization of A j , B j in (3.8) (cf. [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF]). Using (3.4), we get the desired splitting

A j (x, y) = B d C j (x, z)C j (z, y)dz (3.13) and B j (x, y) = B d D j (x, z)D j (z, y)dz. (3.14) Obviously z → C j (x, z)C j (z, y) is a polynomial of degree < 2 j+1 and z → D j (x, z)D j (z, y) is a polynomial of degree < 2 j+2
. The next step is to discretize the kernels A j (x, y) and B j (x, y).

Cubature formula and discretization

To construct the needlets on B d we need one more ingredient -a cubature formula on B d exact for polynomials of a given degree.

Recall first the bijection between the ball B d (equipped with the usual Lebesgue measure) and the unit upper hemisphere in R d+1 :

S d + = {(x, y), x ∈ R d , 0 ≤ y ≤ 1, |x| 2 d + y 2 = 1}
equipped with dσ the usual surface measure.

T : (x, y) ∈ S d + → x ∈ R d and T -1 : x ∈ R d → x = (x, 1 -|x| 2 d ) ∈ S d +
Applying the substitution T one has (see e.g. [START_REF] Xu | Orthogonal polynomials and cubature formulae on spheres and on balls[END_REF])

S d + F (x, y)dσ(x, y) = B d F (x, 1 -|x| 2 ) dx 1 -|x| 2 (3.15)
and hence for f :

R d → R (x,y)∈S d + f (x)ydσ d (x, y) = x∈B d f (x) 1 -|x| 2 d dx 1 -|x| 2 d = B d f (x)dx. (3.16)
Therefore, given a cubature formula on S d one can easily derive a cubature formula on B d . Indeed, suppose we have a cubature formula on S d + exact for all polynomials of degree n + 1, i.e, there exist χn ⊂ S d + and coefficients ω ξ > 0, ξ ∈ χn , such that

S d + P (u)dσ(u) = ξ∈ χn ω ξ P ( ξ) ∀P ∈ Π n+1 (R d+1 ). If P ∈ Π n (R d ) then P (x)y ∈ Π n+1 (R d+1
) and hence ξ∈ χn

ω ξP (ξ) 1 -ξ 2 = S d + P (x)ydσ = B d P (x)dx.
Thus the projection χ n of χn on B d is the set of nodes and the associated coefficients given by ω ξ = 1ξ 2 ω ξ induce a cubature formula on B d exact for Π n (R d ).

The following proposition follows from results in [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF] and [START_REF] Xu | Orthogonal polynomials and cubature formulae on spheres and on balls[END_REF].

Proposition 3.1. Let {B( ξi , ρ) : i ∈ I} be a maximal family of disjoint spherical caps of radius ρ = τ 2 -j with centers on the hemisphere S d + . Then for sufficiently small 0 < τ ≤ 1 the set of points χ j = {ξ i : i ∈ I} obtained by projecting the set { ξi : i ∈ I} on B d is a set of nodes of a cubature formula which is exact for Π 2 j+2 (B d ). Moreover, the coefficients ω ξ i of this cubature formula are positive and ω ξ i ∼ W j (ξ i )2 -jd . Also, the cardinality #χ j ∼ 2 jd .

Needlets

Going back to identities (3.13) and (3.14) and applying the cubature formula described in Proposition 3.1, we get

A j (x, y) = B d C j (x, z)C j (z, y)dz = ξ∈χ j ω ξ C j (x, ξ)C j (y, ξ) and B j (x, y) = B d D j (x, z)D j (z, y)dz = ξ∈χ j ω ξ D j (x, ξ)D j (y, ξ).
We define the father needlets φ j,ξ and the mother needlets ψ j,ξ by

φ j,ξ (x) = √ ω ξ C j (x, ξ) and ψ j,ξ (x) = √ ω ξ D j (x, ξ), ξ ∈ χ j , j ≥ 0.
We also set ψ -1,0 =

½ B d |B d | and χ -1 = {0}.
From above it follows that

A j (x, y) = ξ∈χ j φ j,ξ (x)φ j,ξ (y), B j (x, y) = ξ∈χ j ψ j,ξ (x)ψ j,ξ (y).
Therefore,

A j f (x) = B d A j (x, y)f (y)dy = ξ∈χ j f, φ j,ξ φ j,ξ = ξ∈χ j α j,ξ φ j,ξ , α j,ξ = f, φ j,ξ .
(3.17) and

B j f (x) = B d B j (x, y)f (y)dy = ξ∈χ j f, ψ j,ξ ψ j,ξ = ξ∈χ j β j,ξ ψ j,ξ , β j,ξ = f, ψ j,ξ .
(3.18) By (3.17) and (3.6) we have

φ j,ξ 2 2 ≥ A j φ j,ξ , φ j,ξ = ξ ′ ∈χ j φ j,ξ , φ j,ξ ′ φ j,ξ ′ , φ j,ξ = ξ ′ ∈χ j | φ j,ξ , φ j,ξ ′ | 2 ≥ φ j,ξ 4 2
and hence

φ j,ξ 2 ≤ 1. (3.19)
From (3.5) and the fact that j≥0 b(t2

-j ) = 1 for t ∈ [1, ∞), it readily follows that f = j≥-1 ξ∈χ j f, ψ j,ξ ψ j,ξ , f ∈ Ä 2 (B d ),
and taking inner product with f this leads to

f 2 2 = j ξ∈χ j | f, ψ j,ξ | 2 ,
which in turn shows that the family {ψ j,ξ } is a tight frame for Ä 2 (B d ) and conse- quently

ψ j,ξ 2 2 ≥ ψ j,ξ 4
2 , i.e. ψ j,ξ 2 ≤ 1.

(3.20)

Observe that using the properties of the cubature formula from Proposition 3.1, estimate (3.8) leads to the localization estimate (cf. [START_REF] Petrushev | Localized polynomials frames on the ball[END_REF]):

|φ j,ξ (x)|, |ψ j,ξ (x)| ≤ C M 2 jd/2 W j (ξ)(1 + 2 j d(x, ξ)) M ∀M > 0. (3.21)
Nontrivial lower bounds for the norms of the needlets are obtained in [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]. More precisely, in [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF] it is shown that for 0 < p ≤ ∞

ψ j,ξ p ∼ φ j,ξ p ∼ 2 jd W j (ξ) 1/2-1/p , ξ ∈ χ j . (3.22)
We next record some properties of needlets which will be needed later on. For convenience we will denote in the following by h j,ξ either φ j,ξ or ψ j,ξ . Theorem 3.2. Let 1 ≤ p ≤ ∞ and j ≥ 1. The following inequalities hold

ξ∈χ j h j,ξ p p ≤ c2 j(dp/2+(p/2-2) + ) if p = 4, (3.23
)

ξ∈χ j h j,ξ p p ≤ cj2 jdp/2 if p = 4, (3.24)
and for any collection of complex numbers {d ξ } ξ∈χ j ξ∈χ j

d ξ h j,ξ p ≤ c ξ∈χ j |d ξ | p h j,ξ p p 1/p . ( 3

.25)

Here c > 0 is a constant depending only on d, p, and τ .

To make our presentation more fluid we relegate the proof of this theorem to the appendix.

Linear needlet estimator

Our motivation for introducing the estimator described below is the excellent approximation power of the operators A j defined in §3.2.2 and its compatibility with the Radon SVD. We begin with the following representation of the unknown func-

tion f f = k,l,i f, f k,l,i f k,l,i ,
where the sum is over the index set {(k, l, i)

: k ≥ 0, 0 ≤ l ≤ k, l ≡ k(mod 2), 1 ≤ i ≤ N d-1 (l)}.
Combining this with the definition of A j we get

A j f = ξ∈χ j f, φ j,ξ (y) φ j,ξ = ξ∈χ j α j,ξ φ j,ξ ,
where

α j,ξ = f, φ j,ξ = k,l,i γ j,ξ k,l,i f, f k,l,i = k,l,i γ j,ξ k,l,i 1 λ k R(f ), g k,l,i µ = k,l,i γ j,ξ k,l,i 1 λ k g k,l,i R(f )dµ.
Here γ j,ξ k,l,i = f k,l,i , φ j,ξ (y) can be precomputed. It seems natural to us to define an estimator f j of the unknown function f by

f j = ξ∈χ j α j,ξ φ j,ξ , (3.26)
where

α j,ξ = k,l,i γ j,ξ k,l,i 1 λ k g k,l,i dY. (3.27)
Here the summation is over {(k, l, i) : 0 ≤ k < 2 j , 0 ≤ l ≤ k, l ≡ k(mod 2), 1 ≤ i ≤ N d-1 (l)} and j is a parameter. Some clarification is needed here. The father and mother needlets, introduced in §3.2.6, are closely related but play different roles. Both φ j,ξ and ψ j,ξ have superb localization, however, the mother needlets {ψ j,ξ } have multilevel structure and, therefore, are an excellent tool for nonlinear n-term approximation of functions on the ball, whereas the father needlets are perfectly well suited for linear approximation. So, there should be no surprise that we use the father needlets for our linear estimator.

Furthermore, even if the needlets are central in the analysis of the estimator, the estimator f j can be defined without them. Indeed,

f j = ξ∈χ j k,l,i γ j,ξ k,l,i 1 λ k g k,l,i dY φ j,ξ
as all the sum are finite, their order can be interchanged, yielding

f j = k,l,i 1 λ k g k,l,i dY ξ∈χ j γ j,ξ k,l,i φ j,ξ = k,l,i 1 λ k g k,l,i dY A j f k,l,i
and thus the estimator is obtained by a simple componentwise multiplication on the SVD coefficients

f j = k,l,i a k 2 j λ k g k,l,i dY f k,l,i .
However, as will be shown in the sequel, the precise choice of this smoothing allows to consider L p losses and precisely because of the localization properties of the atoms this approach will be extended using a thresholding procedure in a further work, using this time the mother wavelet.

The risk of the needlet estimator

In this section we estimate the risk of the needlet estimator introduced above in terms of the Besov smoothness of the unknown function.

Besov spaces

We introduce the Besov spaces of positive smoothness on the ball as spaces of L p -approximation from algebraic polynomials. As in §3.2.3 we will denote by E n (f, p) the best L p -approximation of f ∈ Ä p (B d ) from Π n . We will mainly use the notations from [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]. Definition 4.1. [15] Let 0 < s < ∞, 1 ≤ p ≤ ∞, and 0 < q ≤ ∞. The space B s,0 p,q on the ball is defined as the space of all functions f ∈ Ä

p (B d ) such that |f | B s,0 p,q = n≥1 (n s E n (f, p)) q 1 n 1/q < ∞ if q < ∞,
and |f | B s,0 p,q = sup n≥1 n s E n (f, p) < ∞ if q = ∞.
The norm on B s,0 p,q is defined by

f B s,0 p,q = f p + |f | B s,0 p,q . Remark : From the monotonicity of {E n (f, p)} it readily follows that f B s,0 p,q ∼ f p + j≥0 (2 js E 2 j (f, p)) q 1/q
with the obvious modification when q = ∞.

⋆

There are several different equivalent norms on the Besov space B s,0 p,q . Theorem 4.2. With indexes s, p, q as in the above definition the following norms are equivalent to the Besov norm f B s,0 p,q :

(i) N 1 (f ) = f p + (2 js f -A j f p ) j≥0 l q , (ii) N 2 (f ) = f p + (2 js B j f p ) j≥1 l q , (iii) N 3 (f ) = f p + (2 js ξ∈χ j | f, ψ j,ξ | p ψ j,ξ p p ) j≥-1 l q .
Proof. The equivalence N 1 (f ) ∼ f B s,0 p,q is immediate from (3.12). To prove that N 2 (f ) ∼ N 1 (f ), we recall that

B j = (A j+1 -A j ) (see §3.2.2) and hence B j f p ≤ f -A j+1 f p + f -A j f p which readily implies N 2 (f ) ≤ cN 1 (f ).
In the other direction, we have

f -A j f p = ∞ l=j B l f p ≤ ∞ l=j B l f p . Assuming that N 2 (f ) < ∞ we have B l (f ) p = α l 2 -ls with {α l } ∈ l q . Hence ∞ l=j B l (f ) p = ∞ l=j α l 2 -ls = 2 -js ∞ l=j α l 2 -(l-j)s =: 2 -js β j
and by the convolution inequality {β j } ∈ l q . Therefore, N 1 (f ) ≤ cN 2 (f ).

For the equivalence N 3 (f ) ∼ f B s,0 p,q , see [15, Theorem 5.4].

Comparison with the "standard" Besov spaces

The classical Besov space B s p,q (B d ) is defined through the L p -norm of the finite differences:

∆ h f (x) = (f (x + h) -f (x))½ x∈B d ½ x+h∈B d and in general ∆ N h f (x) = ½ x∈B d ½ x+N h∈B d N k=0 (-1) N -k N k f (x + kh).
Then the Nth modulus of smoothness in L p is defined by

ω N p (f, t) = sup |h|≤t ∆ N h f p , t > 0.
For 0 < s < N, 1 ≤ p ≤ ∞, and 0 < q ≤ ∞, the classical Besov space B s p,q is defined by the norm

f B s p,q = f p + 1 0 [t s ω N p (f, t)] q dt t 1/q ∼ f p + ∞ j=0 [2 js ω N p (f, 2 -j )] q 1/q
with the usual modification for q = ∞. It is well known that the definition of B s p,q

does not depend on N as long as s < N [START_REF] Johnen | On the equivalence of the k functional and moduli of continuity and some applications[END_REF]. Moreover, the embedding

B s p,q ⊆ B s,0 p,q , (4.1) is immediate from the estimate E n (f, p) ≤ cω N p (f, 1/n) [13].
4.2 Upper bound for the risk of the needlet estimator

Theorem 4.3. Let 1 ≤ p ≤ ∞, 0 < s < ∞, and assume that f ∈ B s,0 p,∞ with f B s,0 p,∞ ≤ M. Let f J = ξ∈χ J α J,ξ φ j,ξ
be the needlet estimator introduced in §3.3, where J is selected depending on the parameters as described below.

1. If M2 -J(s+d) ∼ ε when p = ∞, then f -f J ∞ ≤ c ∞ M d s+d ε s s+d log M/ε. 2. If M2 -Js ∼ ε2 J(d-2/p) when 4 ≤ p < ∞, then f -f J p p ≤ c p M (d-2/p)p s+d-2/p ε sp s+d-2/p ,
where when p = 4 there is an additional factor ln(M/ε) on the right.

If M2

-Js ∼ ε2 J(d-1/2) when 1 ≤ p < 4, then f -f j p p ≤ c p M (d-1/2)p s+d-1/2 ε sp s+d-1/2 .

Remarks :

• It will be shown in a forthcoming paper that the following rates of convergence are, in fact, minimax, i.e. there exist positive constants c 1 and c 2 such that sup

f B s,0 p,∞ ≤M inf f estimator f -f p p ≥ c 1 max{ε sp s+d-2/p , ε sp s+d-1/2 }, sup f B s,0 ∞,∞ ≤M inf f estimator f -f ∞ ≥ c 2 ε s s+d log 1/ε.
• The case p = 2 above corresponds to the standard SVD method which involves Sobolev spaces. In this setting, minimax rates have already been established (cf. [START_REF] Dicken | Wavelet-Galerkin methods for ill-posed problems[END_REF], [START_REF] Mathé | Geometry of linear ill-posed problems in variable Hilbert scales[END_REF] [4], [START_REF] Cavalier | Oracle inequalities for inverse problems[END_REF], [START_REF] Tsybakov | On the best rate of adaptive estimation in some inverse problems[END_REF], [START_REF] Goldenshluger | On adaptive inverse estimation of linear functionals in Hilbert scales[END_REF], [START_REF] Efromovich | On inverse problems with unknown operators[END_REF]); these rates are ε 2s s+d-1/2 . Also, it has been shown that the SVD algorithms yield minimax rates. These results extend (using straightforward comparisons of norms) to L p losses for p < 4, but still considering the Sobolev ball { f B s,0 2,∞ ≤ M} rather than the Besov ball { f B s,0 p,∞ ≤ M}. Therefore, our results can be viewed as an extension of the above results, allowing a much wider variety of regularity spaces.

• The Besov spaces involved in our bounds are in a sense well adapted to our method. However, the embedding results from Section 4.1.1 shows that the bounds from Theorem 4.3 hold in terms of the standard Besov spaces as well. This means that in using the Besov spaces described above, our results are but stronger.

• In the case p ≥ 4 we exhibit here new minimax rates of convergence, related to the ill posedness coefficient of the inverse problem d-1 2 along with edge effects induced by the geometry of the ball. These rates have to be compared with similar phenomena occurring in other inverse problems involving Jacobi polynomials (e.g. Wicksell problem), see [START_REF] Kerkyacharian | Needlet algorithms for estimation in inverse problems[END_REF].

A j f -f p ≤ c f B s,0 p,∞ 2 -js ≤ cM2 -js . (4.2)
Now from dY = Rf dµ + εdW we have

g k,l,i dY = Z Rf g k,l,i dµ + ε g k,l,i dW = B d f R * g k,l,i dx + ε Z k,l,i = λ k B d f f k,l,i dx + ε Z k,l,i and hence 1 λ k g k,l,i dY = B d f f k,l,i dx + ε λ k Z k,l,i .
On account of (3.27) this leads to

α j,ξ = k,l,i γ j,ξ k,l,i B d f f k,l,i dx + k,l,i γ j,ξ k,l,i ε λ k Z k,l,i = α j,ξ + Z j,ξ .
Here the summation is over

{(k, l, i) : 0 ≤ k < 2 j , 0 ≤ l ≤ k, l ≡ k(mod 2), 1 ≤ i ≤ N d-1 (l)}.
Since Z k,l,i are independent N(0, 1) random variables, Z j,ξ ∼ N(0, σ 2 j,ξ ) with

σ 2 j,ξ = ε 2 k,l,i |γ j,ξ k,l,i | 2 (k) d π d-1 2 d k ≤ (2 j ) d-1 π d-1 2 d ≤ c2 j(d-1) ε 2 (4.3) with c = (d/2π) d-1 .
Here we used that {f k,l,i } is an orthonormal basis for Ä 2 and

hence k,l,i |γ j,ξ k,l,i | 2 = φ j,ξ 2
2 ≤ 1. From (3.26) f j = ξ∈χ j α j,ξ φ j,ξ and using (4.2) we have, whenever 1 ≤ p < ∞,

f -f j p p ≤ 2 p-1 { f -A j f p p + A j f -f j p p } ≤ 2 p-1 {cM p 2 -jsp + A j f -f j p p } (4.4)
and, for p = ∞,

f -f j ∞ ≤ f -A j f ∞ + A j f -f j ∞ ≤ cM2 -js + A j f -f j ∞ . (4.5)
On the other hand, using inequality (3.25) of Theorem 3.2 we obtain, if 1 ≤ p < ∞,

A j f -f j p p = ξ∈χ j (α j,ξ -α j,ξ )φ j,ξ p p ≤ c ξ∈χ j |α j,ξ -α j,ξ | p φ j,ξ p p
and hence

A j f -f j p p ≤ c ξ∈χ j |Z j,ξ | p φ j,ξ p p ≤ c(ε2 j(d-1)/2 ) p ξ∈χ j φ j,ξ p p , (4.6)
where we used that |Z j,ξ | p ≤ c(ε2 j(d-1)/2 ) p . Similarly, for p = ∞,

A j f -f j ∞ = ξ∈χ j (α j,ξ -α j,ξ )φ j,ξ ∞ ≤ c max ξ∈χ j |α j,ξ -α j,ξ | φ j,ξ ∞
and hence

A j f -f j ∞ ≤ c {max ξ∈χ j |Z j,ξ | φ j,ξ ∞ } ≤ cε2 j(d-1)/2 max ξ∈χ j φ j,ξ ∞ 2 log 2 2 jd ≤ cε2 jd j. (4.7)
For the second inequality above we used Pisier's lemma:

If Z j ∼ N(0, σ 2 j ), σ j ≤ σ, then ( sup 1≤j≤N |Z j |) ≤ σ 2 log 2 N.
We also used that max ξ∈χ j φ j,ξ ∞ ≤ c2 j(d+1)/2 , which follows by inequality (3.23) of Theorem 3.2. Combining (4.5) and (4.7) we obtain, for p = ∞,

f -f j ∞ ≤ c{M2 -js + ε2 jd j} and if M2 -j(s+d) ∼ ε, then f -f j ∞ ≤ cM d s+d ε s s+d log M/ε.
Similarly, combining estimate (3.23) of Theorem 3.2 with (4.6) and inserting the resulting estimate in (4.4) we obtain in the case 4

≤ p < ∞ f -f j p p ≤ c{M2 -jsp + (ε2 j(d-1)/2 ) p 2 jdp/2+p/2-2 } = c{M p 2 -jsp + ε p 2 j(dp-2) }. If M2 -js ∼ ε2 j(d-2/p) this yields f -f j p p ≤ cM (d-2/p)p s+d-2/p ε sp s+d-2/p .
Accordingly, for p = 4 we combine inequality (3.24) with (4.6) and insert the result in (4.4) to obtain

f -f j p p ≤ c{M p 2 -jsp + (ε2 j(d-1)/2 ) p j2 jdp/2 } = c{M p 2 -jsp + j(ε2 j(d-1/2 ) p } and if M2 -js ∼ ε2 j(d-1/2) this yields f -f j p p ≤ cM (d-2/p)p s+d-2/p ε sp s+d-1/2 log M/ε.
Finally, if 1 ≤ p < 4 as above we obtain using (3.23), (4.6), and (4.4)

f -f j p p ≤ c{M p 2 -jsp + (ε2 j(d-1)/2 ) p 2 jdp/2 } = c{M p 2 -jsp + (ε2 j(d-1/2 ) p }. So, if M2 -js ∼ ε2 j(d-1/p) , then f -f j p p ≤ cM (d-1/2)p s+d-1/2 ε sp s+d-1/2 .
This completes the proof of the theorem. We have implemented this scheme for d = 2 in the radiological setting of Cormack [START_REF] Cormack | Representation of a function by its line integrals with some radiological applications II[END_REF]. This case corresponds to the fan beam Radon transform used in Computed Axial Tomography (CAT). As shown if Figure 3, an object is positioned in the middle of the device. X rays are sent from a pointwise source S(θ 1 ) located on the boundary and making an angle θ 1 with the horizontal. They go through the object and are received on the other side on uniformly sampled array of receptors R(θ 1 , θ 2 ). The log decay of the energy from the source to a receptor is proportional to the integral of the density f of the object along the ray and thus one finally measures

Rf (θ 1 , θ 2 ) = e θ 1 +λe θ 1 -θ 2 ∈B 2 f (x)dλ
with e θ = (cos θ, sin θ) or equivalently the classical Radon transform

Rf (θ, s) = y∈θ ⊥ sθ+y∈B 1 f (sθ + y)dy, θ ∈ Ë 1 , s ∈ [-1, 1],
for θ = θ 1θ 2 and s = sin θ 2 . The device is then rotated to a different angle θ 1 and the process is repeated. Note that dθ ds (1-s 2 ) is nothing but the measure corresponding to the uniform dθ 1 dθ 2 by the change of variable that maps (θ 1 , θ 2 ) into (θ, s).

The Fan Beam Radon SVD basis of the disk is tensorial in polar coordinates:

f k,l,i (r, θ) = (2k + 2) 1/2 P (0, l) j (2|r| 2 -1)|r| l Y l,i (θ), 0 ≤ l ≤ k, k -l = 2j, 1 ≤ i ≤ 2,
where P 0,l j is the corresponding Jacobi polynomial, and Y l,1 (θ) = c l cos(lθ) and

Y l,2 (θ) = c l sin(θ) with c 0 = 1 √ 2π and c l = 1 √ π otherwise. The basis of S 2 × [-1, 1]
has a similar tensorial structure as it is given by

g k,l,i (θ, t) = [h k ] -1/2 (1 -t 2 ) 1/2 C 1 k (t)Y l,i (θ), k ≥ 0, l ≥ 0, 1 ≤ i ≤ 2,
where C 1 k is the Gegenbauer of parameter 1 and degree k. The corresponding eigenvalues are

λ k = 2 √ π √ k + 1 .
In this paper, we have considered the theoretical framework of the white noise model. In this model, we assume that we have access to the noisy "scalar product" g k,l,i dY , that is to the scalar product of Rf with the SVD basis g kl,i up to a i.i.d. centered Gaussian perturbation of known variance ǫ 2 . This white noise model is a convenient statistical framework closely related to a more classical regression problem with a uniform design on θ 1 and θ 2 , which is closer to the implementation in real devices. In this regression design, one observe

Y i 1 ,i 2 = Rf 2π i 1 N 1 - i 2 N 2 , sin 2π i 2 N 2 + ǫ i 1 ,i 2 , i 1 ≤ N 1 , i 2 ≤ N 2
where N 1 and N 2 gives the discretization level of the angles θ 1 and θ 2 and ǫ i 1 ,i 2 is an i.i.d. centered Gaussian sequence of known variance σ 2 . Note that this points are not cubature points for the SVD coefficients. The correspondence between the two model is obtain by replacing the noisy scalar product g k,l,i dY with by the corresponding Riemann sum

Rf, g k,l,i = 1 N 1 × N 2 N 1 -1 i 1 =0 N 2 -1 i 2 =0 g k,l,i 2π i 1 N 1 - i 2 N 2 , sin 2π i 2 N 2 Y i 1 ,i 2
and using the calibration ǫ 2 = σ 2 /(N 1 × N 2 ). It is proved, for instance in [START_REF] Brown | Asymptotic equivalence of nonparametric regression and white noise[END_REF] that the regression model with uniform design and the white noise model are close in the sense of Le Cam's deficiency -which roughly means that any procedure can be transferred from one model to the other, with the same order of risk-. The estimator f j defined in the white noise model by

f j = k,l,i a k 2 j λ k g k,l,i dY f k,l,i = k,l,i a k 2 j √ k + 1 2 √ π g k,l,i dY f k,l,i
is thus replaced in the regression model by

f j = k,l,i a k 2 j λ k Rf, g k,l,i f k,l,i = k,l,i a k 2 j √ k + 1 2 √ π
Rf, g k,l,i f k,l,i .

Numerical results

To illustrate the advantages of the linear needlet estimator over the linear SVD estimator, we have compared their performances on a synthetic example, the classical Logan Shepp phantom [START_REF] Shepp | The fourier reconstruction of a head section[END_REF], for different L p norm and different noise level, and for both the white noise model and the regression model. The Logan Shepp phantom is a synthetic image used as a benchmark in the tomography community. It is a simple toy model for human body structures simplified as a piecewise constant function with discontinuities along ellipsoids (see Figure 6). This example is not regular in a classical sense. Indeed, it belongs to B 1,0 1,1 but not to any B s,0 p,q with s > 1.

To conduct the experiments, we have adopted the following scheme. Denote by f of the Logan Shepp function presented above, its decomposition in the SVD basis f k,l,i up to degree k = 512 has been approximated with an initial numerical quadrature χ valid for polynomial of degree 4 × k = 2048,

f, f k,l,i ≃ (r i ,θ i )∈χ ω (r i ,θ i ) f (r i , θ i )f k,l,i (r i , θ i ) = c k,l,i
and used this value to approximate the original SVD coefficients of R(f ), the noiseless Radon transform of f ,

R(f ), f k,l,i ≃ λ k c k,l,i .
In the white noise setting, for all k ≤ k 0 = 256, a noisy observation g k,l,i dY is generated by

g k,l,i dY ≃ λ k c k,l,i + ǫW k,l,i
where ǫ is the noise level and W k,l,i a iid sequence of standard Gaussian random variables. Our linear needlet estimator f J of level J = log 2 (k N ) is then computed as

f J = k≤k 0 ,l,i a k 2 J (c k,l,i + ǫ λ k W k,l,i )f k,l,i
while the linear SVD estimator f S k S of degree k S is defined as

f S k S = k≤k S ,l,i (c k,l,i + ǫ λ k W k,l,i )f k,l,i .
We also consider the naive inversion up to degree k 0 f I which is equal to f S k 0 . The L p estimation error is measured by reusing the initial quadrature formula,

f -f p ≃ (r i ,θ i )∈χ ω (r i ,θ i ) |f (r i , θ i ) -f(r i , θ i )| p .
In the regression setting, we have computed the values of the Radon transform Rf of f on a equispaced grid for the angles θ 1 and θ 2 specified by its sizes N 1 and N 2 using its SVD decomposition up to k = k = 512. We have then defined the noisy observation as

Y i 1 ,i 2 = (Rf 2π i 1 N 1 - i 2 N 2 , sin 2π i 2 N 2 + ǫ i 1 ,i 2
with ǫ i 1 ,i 2 an i.i.d. centered Gaussian sequence of known variance σ 2 . The estimated SVD coefficients are obtained through the Riemann sums

Rf, g k,l,i = 1 N 1 × N 2 N 1 -1 i 1 =0 N 2 -1 i 2 =0 g k,l,i 2π i 1 N 1 - i 2 N 2 , sin 2π i 2 N 2 Y i 1 ,i 2 .
We plug then these values instead of the g k,l,i dY in the previous estimators.

For each noise level and each norm, the best level J and the best degree K has been selected as the one minimizing the average error over 50 realizations of the noise. Figure 4 displays, in a logarithmic scale, the estimation errors f -f p in the white noise model plotted against the logarithm of the noise level ǫ. It shows that, except for the very low noise case, both the linear SVD estimator and the linear Needlet estimators reduce the error over a naive inversion linear SVD estimate up to the maximal available degree k 0 . They also show that the Needlet estimator outperforms the SVD estimator in a large majority of cases from the norm point of view and almost always from the visual point of view as shown in Figure 6. The localization of the needlet also 'localizes' the errors and thus the "simple" smooth regions are much better restored with the needlet estimate than with the SVD because the errors are essentially concentrated along the edges for the needlet. Remark that the results obtained for the regression model in Figure 5 are similar. We have plotted, in a logarithmic scale, the estimation errors against the logarithm of the equivalent of the noise ǫ 2 in the regression σ 2 /(N 1 × N 2 ) with N 1 = N 2 = 64 and various σ 2 . Observe that the curves are similar as long as σ 2 is not too small, i.e. as long as the error due to the noise dominate the error due to the discretization. As can be seen both analysis do agree. This confirms the fact that the white noise model analysis is relevant for the corresponding fixed design.

A fine tuning for the choice of the maximum degree is very important to obtain a good estimator. In our proposed scheme, and in the Theorem, this parameter is set by the user according to some expected properties of the unknown function or using some oracle. Nevertheless, an adaptive estimator, which does not require this input, can already be obtained from this family, for example, by using some aggregation technique. A different way to obtain an adaptive estimator based on thresholding is under investigation by some of the authors. From [20, p. 99] with some adjustment of notation, we have

λ 2 k = |Ë d-2 |π (d-1)/2 Γ d+1 2 C d/2 k (1)C (d-2)/2 l (1) 1 -1 C d/2 k (t)C (d-2)/2 l (t)(1 -t 2 ) (d-3)/2 dt,
where 0 ≤ l ≤ k and l ≡ k (mod 2). As will be seen shortly λ k is independent of l.

We will only consider the case d > 2 (the case d = 2 is simpler, see [20, p. 99]). To compute the above integral we will use the well known identity (cf. [25, (4 Γ(l + 2λ) (l + λ)Γ(l + 1) .

We use this and that C λ n (1) = Γ(n+2λ n!Γ(2λ) (see §2. We insert this in (6.2) and then a little algebra shows that λ 2 k = 2 d π d-1 (k+1) d-1 .

Proof of Theorem 3.2

For the proof of estimates This yields estimates (3.23)- (3.24).

We now turn to the proof of estimate (3.25). We will employ the maximal operator M t (t > 0), defined by

M t f (x) := sup B∋x 1 |B| B |f (y)| t dy 1/t , x ∈ B d , (6.3)
where the sup is over all balls B ⊂ B d with respect to the distance d(•, •) from (3.9) containing x. It is easy to show that (see §2.3 in [START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]) the Lebesgue measure on B d is a doubling measure with respect to the distance d(•, •). Hence the general theory of maximal operators applies. In particular, the Fefferman-Stein vector-valued maximal inequality is valid: If 0 < p < ∞, 0 < q ≤ ∞, and 0 < t < min{p, q} then for any sequence of functions {f ν } ν on B d This completes the proof of (3.25).

( ∞ ν=1 |M t f ν (•)| q ) 1/q p ≤ c ( ∞ ν=1 |f ν (•)| q ) 1/q p . ( 6

Figure 1 :

 1 Figure 1: A few radon SVD basis elements (Low quality figure due to arXiv constraint)

Figure 1

 1 displays a few f k,l,i and illustrates their lack of localization. The following theorem gives the SVD decomposition of the Radon transform.

Theorem

Figure 2 :

 2 Figure 2: Smoothing kernel B j (x, y 0 ) for a fixed y 0 for (1) a C ∞ cut-off function a with j = 4 , (2) for the same a with j = 6 and (3) for a non smooth cut-off function a with j = 4 (Low quality figure due to arXiv constraint)

⋆ 4 . 3

 43 Proof of Theorem 4.3 Assume f ∈ B s,0 p,∞ and f B s,0 p,∞ ≤ M. Then by Theorem 4.2,

Figure 3 :

 3 Figure 3: Simplified CAT device

Figure 4 :

 4 Figure 4: Error decay in the white noise model: the red curve corresponds to the needlet estimator and the black one to the SVD estimator.

Figure 5 :Figure 6 : 6 Appendix 6 . 1

 56661 Figure 5: Error decay in the regression model: the red curve corresponds to the needlet estimator and the black one to the SVD estimator.

(6. 1 ) 2 ( 1 -t 2 )

 1212 This with λ = (d -2)/2 and the orthogonality of the polynomials C

3 . 1 ) 2 Γ 2 )

 3122 and |Ë d-2 | = 2π (d-1)/2 Γ((d-1)/2) to obtain (d)l!Γ(d -2) Γ(k + d)Γ(l + d -2) 2 4-d π (d -2)Γ d-2 2The doubling formula for Gamma-function says: Γ(2z) = 2 2z-1 √ π Γ(z)Γ(z + 1) (see e.g.[START_REF] Szegő | Orthogonal polynomials[END_REF]) and hence Γ(d)Γ(d-2) = (d-1)(d-2)Γ(d-2)

1 2 1 2 -j + √ 1 -|ξ| 2 r 1 dπ d- 1 with 2 0( 1 ( 2

 111211212 (3.23)-(3.24) we first note that by (3.3) -j + 1 -|ξ| 2 p/2-1 and we need an upper bound for Ω r := 2 -jd ξ∈χ j . To this end, we will use the natural bijection between B d and Ë d + considered in the remark in §3.2.2. Thus for x ∈ B d we write x = (x, 1 -|x| 2 ) ∈ Ë d + . Let p = (0, 1) be the "north pole" of Ë d . For ξ ∈ χ j we denote byB Ë d ( ξ, ρ) is the geodesic ball on Ë d of radius ρ centered at ξ, i.e. B Ë d ( ξ, ρ) := {x ∈ Ë d : d Ë d (x, p) < ρ}, where d Ë d (x, p) = Arccos ( 1 -|x| 2 ) = Arccos x, p is the geodesic distance between x, p. Using that ||u| -|ξ|| ≤ | ũ, p -ξ, p | ≤ d Ë d ( ξ, ũ) ≤ ρ for ũ ∈ B Ë d ( ξ, ρ), and ρ = τ 2 -j ≤ 2 -j (see Proposition 3.1), it follows that 1 2 -j + 1 -|ξ| 2 ≤ 2 2 -j + 1 -|u| 2 = 2 2 -j + ũ, p ∀ũ ∈ B Ë d ( ξ, ρ).On the other hand, we have|B Ë d ( ξ, ρ)| = |Ë d-1 | ρ 0 (sin θ) d-1 dθ ≥ ρ d |Ë d-1 | 2 d-|Ë d-1 | = 2π d/2 Γ(d/2) . We use the above and the fact that the balls {B Ë d ( ξ, ρ)} ξ∈χ j are disjoint to obtainΩ r ≤ 2 -jd ξ∈χ j 1 |B Ë d ( ξ, ρ)| B Ë d ( ξ,ρ) 1 (2 -j + ũ, p ) r dσ(ũ) j + ũ, p ) r dσ(ũ) ≤ c|Ë d-1 | π/sin θ) d-1(2 -j + cos θ) r dθ ≤ c π/2 0 sin θ (2 -j + cos θ) r dθ = c 1 0 -j + t) r dt ≤ c(d, τ, r) 2 2 -j t -r dt.

. 4 )

 4 Denote by B(ξ, r) the projection of B Ë d ( ξ, r) onto B d , i.e. B(ξ, r) := {x ∈ B d : d(x, ξ) < r}. By [15, Lemma 2.5], we have(M t ½ B(ξ,r) )(x) ≥ c 1 + d(ξ, x) r -(d+1)/t , ξ ∈ B d , 0 < r ≤ π. (6.5)It is easy to see (cf.[START_REF] Kyriazis | Decomposition of weighted triebellizorkin and besov spaces on the ball[END_REF]) that|B(ξ, ρ)| ∼ 2 -jd (2 -j + 1 -|ξ| 2 ) ∼ 2 -jd W j (ξ), ξ ∈ χ j . (6.6) Also, we let 1E := 1 |E| ½ E denote the L 2 -normalized characteristic function of E ⊂ B d . Then (3.21) and (6.6) imply h j,ξ p ∼ 1B(ξ,ρ) p , ξ ∈ χ j . (6.7) Now, pick 0 < t < 1 and M > (d + 1)/t. From (3.21) and (6.5) it follows that |h j,ξ (x)| ≤ c(M t 1B(ξ,ρ) )(x), x ∈ B d . (6.8) Using this, the maximal inequality (6.4), and (6.7) we obtain ξ∈χ j d ξ h j,ξ p ≤ c ξ∈χ j M t (d ξ 1B(ξ,ρ) ) p ≤ c ξ∈χ j d ξ 1B(ξ,ρ) p ≤ c ξ∈χ j d ξ 1B(ξ,ρ)

  .7.29)]) (n + λ)C λ n (t) = λ(C λ+1 n (t) -C λ+1 n-2 (t)). Summing up these identities (with indices n, n -2, . . . ) and taking into account that C λ 0 (t) = 1, C λ 1 (t) = 2λ(t), we get

	C λ+1 n (t) =	⌊n/2⌋ j=0	n -2j + λ λ	C λ n-2j (t).