Samuel Garcia
email: samuel.garcia@ensea.fr

Bertrand Granado
email: bertrand.granado@ensea.fr

Task model and online operating system API for hardware tasks in OLLAF platform

This article present an original hardware task model and the corresponding online API for Fine Grained Dynamically Reconfigurable Architecture. We cover the integration of this API in the OLLAF platform and more specifically its application to memory access management in a dynamically reconfigurable environment. Methods offered by this platform are compared to existing software and hardware solutions. We also discuss of the design complexity of an application using difference solutions. We demonstrate that our solution cans give application developer the same flexibility than with a software implementation, with a very close design complexity while ensuring the same performance gain a common FPGA based IP would permit.

I. INTRODUCTION

More than moore philosophy tries to focus more on how can we get the most out of a given technology as Moore's law tend to focus on bringing a new technology to achieve more. Within this philosophy an application developer can usually enhance his system performance by designing a custom designed computing architecture for his specific application algorithm instead of using a generic architecture such as a general purpose micro-processor(more than Moore option). Fact is that most often application designers simply choose a newer technology and/or bigger micro-processor when it comes to performance enhancement(Moore's law). Custom computing architecture implemented on a fine grained dynamically reconfigurable architecture (FGDRA) can offer the same degree of flexibility than a micro-processor does while offering the performance benefit of a custom architecture. The dynamic nature of both application and platform imply a run time management in the form of an Operating System (OS) that propose an abstracted and simplified model of the design possibilities offered by the host architecture in the form of an Application Programmer Interface (API). This article present an online API for hardware task implemented on the OLLAF platform [START_REF] Garcia | Ollaf: A fine grained dynamically reconfigurable architecture for os support[END_REF]. We will here discuss of an innovative vision of data management architecture for dynamically reconfigurable platform as an example of OLLAF online API use. This paper will be organized as follows. Section II covers related works on this subject, it also gives an overview of the OLLAF platform. Section III present our work on online API for dynamically reconfigurable hardware task and more specifically memory access services. Section IV is a comparison of the several example solutions presented. Conclusions are then drawn in section V, as well as perspectives on this work.

II. RELATED WORK

A. bibliography

Several researches have been led in the field of OS for FGDRA [START_REF] Simmler | Multitasking on FPGA Coprocessors[END_REF], [START_REF] Chen | Configuration-Sensitive Process Scheduling for FPGA-Based Computing Platforms[END_REF], [START_REF] Walder | Reconfigurable Hardware Operating Systems: From Design Concepts to Realizations[END_REF], [START_REF] Wigley | Introducing reconfigme: An operating system for reconfigurable computing[END_REF]. All those studies present an OS more or less customized to enable specific FGDRA related services. Example of such services are partial reconfiguration management, hardware task preemption or hardware task migration. They are all designed on top of a commercial FPGA coupled with a micro-processor. This micro-processor may be a soft-core processor, an embedded hardwired core or even an external processor.

Those works tend to consider hardware task as subtask of a software task. Operating system kernel being usually adaptation of a classical software operating system, a simple solution is to give the software mother-task the responsibility to send system call. In [START_REF] Segard | Soc and rtos: Managing ips and tasks communications[END_REF], authors propose a method for both software to hardware and hardware to hardware intertask communication service. This method rely on a IP wrapper that encapsulate the task architecture and communicate with a particular software task being called the software alter-ego of the hardware task. While this method can work for small message passing exchanges or simple system call, it would provide poor performance for high volume data transfers. The system being purely software, thus sequential, it would not allow to perform more than one operating system service at a time.

The simplest solution, recommended in [START_REF] Xilinx | Two flows for partial reconfiguration: Module based or difference based," Xilinx, Application Note, 2004, application Note: Virtex, Virtex-E, Virtex-II, Virtex-II Pro Families XAPP[END_REF] is to use a simple bus as communication media and a hardware task model like a simple micro-processor peripheral. This model has the advantage of being very straitforward, but it bound hardware tasks as a control software's slave.

In [START_REF] Gantel | Dataflow programming model for reconfigurable computing[END_REF] authors propose an actor task model on a data-flow computing scheme. The API consist on a control and a data channel given by a wrapper that take care of communication and synchronization work. This method differs with our approach as it imply strong constraint on architectural choice in the task implementation. This restrictive task model allow to simplify this API as only few specific services will ever be called by a task. This architecture would be bound to few particular classes of applications.

Our approach is to consider a hardware task as a full fledged, independent system task. We also want a task model being as generic as possible. The application designer can chose the architecture that is most optimal to implement it.

In the state of the art are widely accepted API such as those of POSIX norm [START_REF]ieee standard for information technologyportable operating system interface (posix(r))[END_REF] or µC OS-II [START_REF] Labrosse | MicroC/OS-II The Real-Time Kernel[END_REF]. The work presented here aim at porting those API fundamental principles to a dynamically reconfigurable hardware system. This API is to be integrated as a part of the OLLAF platform presented in the next paragraph.

B. OLLAF architecture overview

OLLAF architecture stand for Operating system enabled Low LAtency FGDRA. It is a Fine Grained Dynamically Reconfigurable Architecture (FGDRA) platform specifically designed to support operating system services to manage hardware task. In this overview we will first explain the specification we follow in the design of this platform, this should give a good understanding of the original philosophy of this design. We will present the mechanisms used in OLLAF to answer to those specifications and then show the contribution of those architectural solutions in terms of operating system support.

1) Specifications of a FGDRA with OS support: We have designed a FGDRA with OS support following those specifications.

It should first address the problem of the configuration speed of a task. This is one of the primary concerns because if the system spend more time configuring itself than actually running tasks its efficiency will be poor. The configuration speed will thus have a big impact on the scheduling strategy.

In order to enable more choice on scheduling scheme, and to match some real time requirements, our FGDRA platform must also include preemption facilities. For the same reasons than configuration, the speed of context saving and restoring process will be one of our primary concerns. On this particular point, previous work that have been discussed in [START_REF] Garcia | Hardware task context management for fine grained dynamically reconfigurable architecture[END_REF] will be adapted and reused.

Scheduling on a classical micro-processor is just a matter of time. The problem is to distribute the computation time between different tasks. In the case of a FGDRA the system must distribute both computation time and computation resources. Scheduling in such a system is then no more a one dimensional problem, but a three dimensional one. One dimension is the time and the two others represent the surface of reconfigurable resources. Performing an efficient scheduling at run time for minimizing processing time is then a very hard problem that the FGDRA should help getting close to solve. The primary concern on this subject is to ensure an easy task relocation. For that, the reconfigurable logic core should be spited into several equivalent blocks. This will allow to move a task from one block to any another block, or from a group of blocks to another group of blocks of the same size and the same form factor, without any change on the configuration data. The size of those blocks would be a tradeoff between flexibility and scheduling efficiency.

Another aspect of an operating system is to provide inter task communication services. In our case we will distinguish two cases. First the case of a task running on top of our FGDRA and communicating with another task running on a different computing unit. This last case will not be covered here as this problem concern a whole heterogeneous platform, not only the particular FGDRA computing units. The second case is when two, or more, tasks run on top of the same FGDRA communicate together. This communication channel should remain the same wherever the task is placed on the FGDRA reconfigurable core and whatever state those tasks are (running, pending, waiting, ...). That mean that the FGDRA platform must provide a rationalized communication medium including exchange memories.

The same arguments could also be applied to inputs/outputs. Here again two cases exist. First the case of I/O being a global resource of the whole platform. Secondly the case of special I/O directly bound to the FGDRA. In the center stand the reconfigurable logic core of the FGDRA. This core is a dual plane, an active plane and a hidden one, organized in slices. Each slice can be reconfigured separately and offers the same set of services. A task is mapped on an integer number of slices. This topology as been chosen for two reasons. First, using a partial reconfiguration by slice transforms the scheduling problem into a two dimensional problem (time + 1D surface) which will be easier to handle for minimizing the processing time. Secondly as every slice are the same and offers the same set of services, tasks can be moved from one slice to another without any change on the configuration data.

In the figure, at the bottom of each slice you can notice two hardware blocks called CMU1 and HCM2 . The CMU is an IP able to manage automatically task's context saving and restoring. The HCM standing for Hardware Configuration Manager is pretty much the same but to handle configuration data also called bitstream. More details about this controller can be found in [START_REF] Garcia | Hardware task context management for fine grained dynamically reconfigurable architecture[END_REF]. On each slice a local cache memory named LCM is added. This memory is a first level of cache memory to store contexts and configurations close to the slice where it might most probably be required. The internal architecture of the core provides adequate materials to work with CMU and HCM. More about this will be discussed in the next section.

On the right of the figure stands a big block called "HW Sup + HW RTK + CCR". This block contains a hardware supervisor running a real time kernel specially adapted to handle FG-DRA related OS services and platform level communication services. In our first prototype presented here, this hardware supervisor is a classical 32 bits micro-processor. Along with this hardware supervisor a central memory is provided for OS use only. This memory will store configurations and contexts of every task that may run on the FGDRA. This supervisor communicates with all slices using a dedicated control bus. The hardware supervisor can initiate context transfers, from and to the hidden plane, by writing in CMU's and HCM's registers through this control bus.

Finally, on top of the figure 1 you can see the application communication media. This communication media provides a communication port to each slice. Those communication ports are bound to the reconfigurable interconnection matrix of the core through API endpoint presented here. Smart Memory Block (SMB) and Smart IO Block are also connected to this communication media.

This architecture has been developed as a VHDL model in which the size and number of slices are generic parameters. On this article we will consider slices composed of 512LE with 32 input line and 32 output line to/from either the two adjacent slices and the communication media, this is illustrated on figure 2 3) OLLAF platform and OS interaction: In previous sections an architectural view of our FGDRA has been exposed. In this section, we discuss about the impact of this architecture on OS services. We will here consider the four services most specifically related to the FGDRA :

• First, the configuration management service : on the hardware side, each slice provides a HCM and a LCM. That means that configurations have to be prefetched in the LCM. The associated service running on the hardware supervisor will thus need to take that into account. This service must manage an intelligent cache to prefetch task configuration on the slices where it might most probably be mapped.

• Secondly, the preemption service: The context management service must ensure that it never exists more than one valid context for each task in the entire FGDRA. Contexts must thus be transferred as soon as possible from LCM to the centralized global memory of the hardware supervisor.

• Scheduling service, and in particular the space management part of the scheduling : it takes advantage of the slice topology and the centralized communication scheme. The reconfigurable resource could then be managed as a virtual infinite space containing an undetermined number of slices. The job is to dynamically map the virtual space into the real space (the actual reconfigurable logic core of the FGDRA).

• Communication service and system management of I/O and memory access : it abstract all the communications of a task and the complexity due to dynamical relocation of task necessary for logic resource management.

III. PROPOSED SOLUTIONS

A. Task model and API in OLLAF

Our goal is to transpose to hardware task the same management facilities that an application developer can dispose when designing a software task. Operating Systems for reconfigurable hardware task management is nothing new, lot of works have been and still are led on the subject, but very few of them define clearly what is the Task model in such a system. For a full software system, a task is basically a C language function. This simple assessment tells us how we can exchange data with it, how it can be ran and how it acts within the system. Moreover, it also makes clear of what a system call is and how it can be programmed by the application designer, with just an other C function. This define the system API as a collection of C functions to call when required, in other terms a C library. To call a system service, all you have to do is include the corresponding API's library and call the corresponding function with the right parameters. Not only this is simple, this is also natural to any software developer. The code sample in Listing 1 show a typical software task performing a simple four taps FIR filter on a set of data, using the µC OS-II API. We can see that input data are passed as an argument of the task function while the output data are passed back using an OS provided message passing service.

Figure3 show a typical functional architecture of a hardware version of this task. To make this architecture a hardware task to be ran on an OS managed FGDRA we have to tell where the task stop and where the system begin. Figure4 show this division as we propose it to be. The system part will have to

* ( i n p u t D a t a + i ) * P1 + * ( i n p u t D a t a + i -1 ) * P2 + * ( i n p u t D a t a + i -2 ) * P3 + * ( i n p u t D a t a + i -3 ) * P4 ) / 4 ; OSMboxPost ( * mBoxID , o u t p u t B u f f e r ) ;
OSTaskDel ( OS PRIO SELF ) ; } Listing 1: A four taps FIR filter software implementation using µC OS-II API be included into the OLLAF architecture. Then the point we want to discuss in more depth in this paper is the interface between the task part and the system part of the functional architecture. This interface will constitute a hardware API for the operating system. It will also define the task model in our system.

To answer to that question we will first focus on what we already have in the OLLAF architecture. In the previous section we have seen that OLLAF's reconfigurable logic core provide in each slice 32 inputs signals and 32 output signals to the outer world (ins[31..0] and outs[31..0]). We also know that each slice is composed of 512 Logic Elements (LE). Thus we want an interface that use only 32 input lines and 32 output lines, and for the use of witch the required logic on the task side can be made using far less than 512 LEs.

The solution we chose in OLLAF is to use a 32 bit access register which act as a command register and/or as a data register with sequential behavior. The first time a given task instance access the API, it must always provide a command word. Some simple access can be performed directly inserting eventual data in the same 32 bits word. Access requiring more data as for example a standard memory acces that require an Fig. 3. Functional architecture of a hardware implementation of a four taps FIR filter Fig. 4. Four tap FIR filter implemented as hardware task, a) and b) are two possible division between task and operating system address word and a 32 bits data word, the same register will be used sequentially.

Figure5's chart shows the the command word API. We can see that the two-first bits determine the type of destination with which the task want to communicate. System (00) is for system call. Memory (01) is for memory access. Task [START_REF] Labrosse | MicroC/OS-II The Real-Time Kernel[END_REF] is for intertask communication. And I/O (11) is for external input/output.

If we conciders the 4 taps FIR filter, we can see that we need to write output data to the memory. This can be done by initiating a full 32 bits transfer. Using our API it then consist on first writing the control word "01 0 1 1 xxx xxxx xxxx xxxx xxxx xxxx", then writing the 32 bits address and finally the 32 bits data. This imply to add to the actual task architecture a multiplexor and a small state machine, the whole represent only 34 Logic elements.

B. Data streams

The communication media used in OLLAF is a circuit switching NoC3 rather than the more popular packet switching NoC. It is based on a dynamic adaptation of the work presented in [START_REF] Pujol | Rseau dinterconnexion haut dbit pour les architectures parallles connexionnistes[END_REF]. This choice allows to dynamically create a point to point connection between two network nodes on which the full bandwidth can then dedicated to data transfers. We call this data stream transaction. Task to memory communication is one particular case where such a data stream transaction can give a real performance gain while still allowing a simple dynamic management. One can figure such a data stream transaction as a hardware DMA transfer with a dynamically created dedicated bus. The transfer initialization is performed as it would for a common DMA, using the upper mentioned API. Once the request is accepted by the global communication arbiter, which in this case is part of the operating system, a dedicated communication channel is created for the transaction. If the task is to be preempted and moved elsewhere, the system will automatically recreate this channel, with the new hostplace of the task, as a part of the context restoration process. The dedicated channel automatically follow the task so that application developer doesn't have to bother of the dynamical management of the task. As long as it is designed in respect with the API, the system guaranty that it runs as if it was the only task in a static system.

In the previous section we have seen that we could perform standard memory access using the OLLAF runtime hardware API. We will show here another type of access that can be performed using this API still based on our simple four tap FIR filter example.

As said earlier, task to memory data stream transaction are very similar to DMA transfer. In OLLAF the actual management of the transfer is done by the memory itself. Embedded data memory is decomposed into several separated smart memory blocs (SMB). Each can serve one or more transactions including data stream transaction. We will not here get into more details about the smart memory engines used in OLLAF as this article want to focus on the API mechanism and associated principles.

In our example we can use a data stream transaction for both input and output data as they are both continuous blocs of data.

As we would do using a traditional DMA we can initiate a transfer from the input memory to our task, initiate a writing transaction on a second API endpoint and then compute the data as a stream allowing an optimal computing power use. Note that this requires two API endpoint so the task should be at least two slice big on the reconfigurable logic core. Figure6 show the actual functional architecture of the new implementation. According to the API tree chart in Fig5 we see that to initiate a 32 bits simple block read transaction, we need to write the word "01 1 0000010 DATA SIZE" where DATA SIZE is a 22 bits integer number representing the size of the input data block to read, then we must write the beginning address of this data block in the memory space. Data will then be sent one 32 bits word each clock cycle. The same principle stand for the writing transaction. We will now discuss of some more detail about internal mechanism of the API. The communication media offers a 100% dedicated channel to system communication. All system calls including data transfer initialization use this channel and are addressed to the appropriate node if required (in this case the smart memory bloc). In fact every nodes that can serve system call listen to this channel. The main system "router" only have to determine which of these nodes will be concerned by the current call. This mechanism paired with a fully parallelized hardware system router permits a zero delay system call serving.

As looking at Fig6, we can see that the task now consist only of its own functionality and a little state machine handling API communication for transactions initialization. This little state machine can easily be automatically generated at design time using simple macro in the HDL source code of the task. In the same way, if we want the input to be now a data flow from external I/O port, the OLLAF API permit a similar data flow transaction for external I/O to task communication involving smart I/O engines. Not only an application developer does not have to bother about dynamical management, but even a simple, possibly static, task becomes easier to design compared with a standard FPGA design flow. This comes to the cost of a dedicated platform architecture specially designed to support operating system managed dynamical system on a chip.

IV. RESULTS

In this section we summarize the results of the four implementation of the four tap FIR filter presented in this article. The first case is the C code running on a general purpose micro-processor (GPP), second case is a traditional implementation of the functional architecture presented on the previous section targeting commonly available FPGA, third case is the OLLAF implementation of this architecture using API's standard memory access, and the fourth case is the OLLAF implementation using API's transactional access. For each case we prompt the execution time for a 10,000 words input file at a fixed frequency of 100 MHz. For the three hardware implementation the logic size of task is given in terms of both logic elements count and slices count. A slice being 512 LE and providing one API endpoint. We also compare for each case the size of the execution context as a clue of the preemption cost. At last we will also discuss of the ease of development. Table I display those results.

In terms of performance, hardware implementation show up-to 16 times improvement compared to the software solution. This result was easily expectable, so it is that this gain could be much greater with a more complex algorithm. The same comparison with a 16 taps FIR filter would have lead to a gain up-to 64. Concerning the OLLAF task implementation, using standard access leads to a doubled computation time compared with the classic FPGA implementation. This is due to the fact that addresses and data have to be transfered sequentially. Transactional access in the fourth implementation case withdraw this problem giving just a 2 clock cycle overhead for transactions initialization.

Considering the logical size of the three hardware implementation, we can see that they are roughly equivalent. Stan-dard access OLLAF implementation present a 34LE overhead as mentioned earlier due to the state machine and multiplexor needed to cope with the API endpoint protocol. The transactional implementation is slightly smaller as the system now take care of all the memories address processing. Note this size reduction is due to a simplification of the task architecture. In OLLAF, the reconfiguration atom is the slice. Slices consist in 512 LE and provide one API endpoint. Any task is constituted of an integer number of slices. Thus the little size difference will not have any real effect here as both implementation will occupy two full slices.

The execution context represents the data that must be saved prior to preemption. The time cost of a task preemption will be directly affected by the size of this context. Execution context size then give a good clue of the preemptability of a task. On a classic 32 bit general purpose microprocessor the execution context is basically the set of registers of the processor, this represent approximately 1 Kbit of data. Two slices of reconfigurable logic in OLLAF represent 1 Kbit of data. But in a classical FPGA, as memories are bound to the reconfigurable logic core, ensuring a correct link between a memory bloc and the task architecture can be difficult as the task is moved. It will require to create a fixed communication structure, such as the communication media in OLLAF, using reconfigurable logic thus limiting its complexity and performance compared with hardwired resources. Not only this solution offers suboptimal performances, but it also loose most of the advantage of having memory blocs inserted into the reconfigurable logic core. Without this fixed structure, the only way to ensure the link between the task and the memory is to move memory data with the task and to keep the same relative position of the task prior to the memory bloc. This make the three dimensional scheduling a real nightmare, and makes execution context size more than 300 times bigger. This is where a specifically designed platform such as OLLAF take all it's meaning compared with a common FPGA primarily designed for a static use.

Ease of development might be a little trickier to discuss as it is a very subjective question. It is then hard to give any real scientific proof about it. However it is commonly accepted that designing a software task is simpler than developing its hardware equivalent. According to [START_REF] Bdti | Bdti certified(tm) results for the autoesl autopilot highlevel synthesis tool[END_REF], using available C to HDL high level synthesis tools does not overcome this statement. We can also easily consider that implementing a dynamically managed hardware task using common FPGA platforms such as Virtex family using available tools is more difficult than designing a static hardware architecture without operating system. Designing a hardware task for OLLAF platform, thanks to this API, require no other operations than to develop the HDL code of the task, using either macro blocs of macro-functions to address system calls. All the complexity due to dynamical management are abstracted. Moreover, using operating system services can also abstract parts of the complexity of the task itself. While this article takes the example of memory access, other services can permit to abstract inter task communications, I/O access, task synchronization, (...), using the API presented here. We mentioned earlier that the software task model and system call have a major advantage of being natural to any developer. The model proposed here should also feel natural to any HDL designer as it relies on straightforward design model using traditional widespread tools.

V. CONCLUSION AND PERSPECTIVES In this article we presented an original solution for online hardware task system calls. This solution relies on a task model as close as possible of any random hardware IP as it would be designed targeting existing fine grained reconfigurable platforms such as FPGA. All system services can be called via hardware API endpoints bound to the reconfigurable switching matrix of each reconfigurable logic slice. This simple interface can cover any communication needs from simple system call to more complex data transfer, I/O access or memory access. In this article we demonstrated the contribution of this API for memory access. More generally we presented the bases of a complete online operating system API supported by a platform architecture which permit to support all required services with both low overhead. This complete solution offers to application developers a highly abstracted view of the system, freeing them to the complexity of system level resources management including computational resources, data storage and exchange resources and I/O.

Perspectives of this work are numerous. A VHDL model of the platform is being developed, at this day all the reconfigurable logic core, configuration and context management units as well as a hardware supervisor are developed. The application communication media is under development as well as API endpoint and Smart Memory Blocks.

Along with the development of the platform, a complete design flow is being developed based on Mentor Graphics HDL designer is being developed, including synthesis, place and route as well as HDL system call. This custom design suite is composed of TCL scripts around Mentor Graphics tools and a HDL library mentioned in this article. As this design flow takes a standard HDL description as input, it is very easy to add existing higher level synthesis tool on top of it and/or to build custom ones.

In order to validate such computing solution, we would also need a benchmarking method adapted to dynamically reconfigurable computing.

Fig. 1 .

 1 Fig. 1. Global view of OLLAF

Fig. 2 .

 2 Fig. 2. external view of a slice in OLLAF

#

  i n c l u d e <OSMainAPI . h> / / C o e f f i t i e n t s # d e f i n e P1 1 / 2 # d e f i n e P2 1 / 2 # d e f i n e P3 1 / 8 # d e f i n e P4 1 / 4 v o i d f i r F i l t e r T a s k ( i n t * i n p u t D a t a , i n t i n p u t T a b l e S i z e , v o i d * p d a t a ) { i n t i ; i n t * o u t p u t B u f f e r = m a l l o c ( s i z e o f ( i n t ) * i n p u t T a b l e S i z e ) ; f o r ( i = 3 ; i <i n p u t T a b l e S i z e ; i ++) * ( o u t p u t B u f f e r + i ) = (

Fig. 5 .

 5 Fig. 5. Tree view of OLLAF online API, each arrow means a word to write or read on the API endpoint, following a branch gives the actual word to write.

Fig. 6 .

 6 Fig.6. Functional view of the FIR filter task architecture using transactional block read and write memory access in OLLAF

CMU : Contexte Management Unit

HCM : Hardware Configuration Manager

NoC : Network on Chip


Gaël Rigaud
email: gael.rigaud@ensea.fr

Maï K Nguyen

Alfred K Louis
email: louis@num.uni-sb.de

Circular Harmonic Decomposition Approach for Numerical Inversion of Circular Radon Transforms

Keywords: Radon transforms (RT), Circular Radon transform (CRT), Circular-arc Radon transform (CART), image reconstruction, Circular Harmonic Decomposition (CHD), Compton

Numerical inversions via circular harmonic decomposition for two classes of circular Radon transforms are established. The first class deals with the Radon transform (RT) defined on circular arcs having a chord of fixed length rotating around its middle point (CART) and the second is the RT defined on a set of circles passing through a fixed point of the plane (CRT). These circular Radon transforms arise from the modeling of different modalities in Compton scattering tomography (CST). Inversions via circular harmonic decomposition are used for image reconstructions in CST. Simulation results show the efficiency and interest of this inversion method in imaging science.

INTRODUCTION

Since the seminal work of J Radon [START_REF] Radon | Uber die bestimmung von funktionnen durch ihre integralwerte längs gewisser mannigfaltikeiten[END_REF], which finds numerous applications (Computed Tomography, Single Photon Emission Computed Tomography (SPECT) and Positron Emission Tomography (PET), etc.) many extensions of this integral transform have been widely discussed, in particular in the literature of imaging science. This is the case when results of measurements appear under the form of integrals of a physical quantity over lower dimensional manifolds. The relevant problem to solve is the recovery of the physical quantity of interest as a function in R 2 . The field of such problems is known in mathematics as integral geometry and in image processing as image reconstruction.

As the circle is the simplest curve in the plane next to the straight line, it becomes the natural object on which a new Radon transform can be defined. Circular Radon transforms along paths that are not on the zero sets of harmonic polynomials were proved to be invertible. Quinto [START_REF] Quinto | Radon transforms of curves in the Problems[END_REF] has discussed the case of all translations of a circle of fixed radius as well as circles centered on a circle. Finally Compton scattering tomography, a two-dimensional imaging process based on scattered radiation, has two modalities for which image formation is built on Radon transforms defined on a set of circles passing through a fixed point of the plane and on circular arcs having a chord of fixed length rotating around its middle point.

In image reconstruction all these modelings of physical phenomenon lead to an inversion of the established forward transform in order to reconstruct the original function. But generally the obtained inversion formula is not adapted for computation because of singularities. A regularization has to be applied and computational methods can be derived.

In this paper we describe two novel computational methods via circular harmonic decomposition (CHD) for circular Radon transforms. The CHD was introduced by Cormack [3, 4] and is very adapted to a polar approach of the studied forward transforms. Indeed Cormack has shown several properties of the circular harmonic components of the data for a given family of curves in the plane which lead to an inversion and to a regularization procedure. This way permits to fulfill the consistency criterion of the data and so to reduce the number of artifacts, which is not the case of the well-known "Filtered Back-Projection" algorithm. Our work is based on these two procedures and on a computational approach of Chapman and Cary [5] who have discussed an alternative algorithm to the "Filtered Back-Projection" algorithm for the inverse Radon transform. We first use the same algorithm in the case of a circular-arc Radon transform then we will establish a algorithm for the inverse circular Radon transform defined on a set of circles passing through a fixed point. Finally, we present numerical results of these two algorithms used in emerging CST and compare them with those obtained for conventional tomography modeled by the classical Radon transform and computed by the filtered backprojection algorithm in the first case and with the inversion formula established by S.J. Norton [6] in the second case.

COMPTON SCATTERING TOMOGRA-PHIES BASED ON CIRCULAR RADON TRANSFORMS

In this part, we will present the working principles of two modalities in Compton scattering tomography and the associated Radon's problems. We will see how to model the physical phenomenon as an integral of the studied function on the corresponding curve and how to inverse these transforms.

Novel CST modality based on Circular-Arc Radon transform

Consider a two-dimensional object represented by a nonnegative continuous function f (M) with bounded support in R 2 . Fig. 1 shows how Compton scattering tomography works. An emitting point source S is placed at a distance 2p from a point detector D. We consider only the upper part of space. This is possible because an angle collimator is placed at D. The segment SD rotates around its middle O and its angular position is given by ϕ.

Emitted photons are scattered at site M and some of them are detected by the detector D at an energy Eω. Therefore the detector can record scattered photons according to scattered energy which is related to the scattering angle ω by the Compton formula. Thus, for a fixed ϕ, to each energy Eω corresponds a set of scattering sites on a circular arc C(ϕ, ω). Finally the detected radiation flux density g(ϕ, ω) is proportional to the integral of the electron density f (M) with M ∈ C(ϕ, ω) which can be written as the Radon transform extended on the corresponding circular arcs C(ϕ, ω)

C2 f (ϕ, ω) = (r,θ)∈C(ϕ,ω) f (r, θ)ds ( 1
)
where C2 is the circular-arc Radon transform (CART) es-tablished in (ϕ, ω)-space and ds is the elementary length of circular arc to be computed from the circular arc equation

r = p( 1 + τ 2 cos 2 γ -τ cos γ) , (2)
where τ = cot ω and γ = θϕ and with (ω, ϕ) ∈]0, π 2 ] × [0, 2π]. Putting τ in equation (1), we obtain :

C 2 f (ϕ, τ ) = π 2 -π 2 f (r(γ), γ + ϕ) r(γ) √ 1 + τ 2 1 + τ 2 cos 2 γ dγ .
(3) where C 2 is the circular-arc Radon transform (CART) established in (ϕ, τ )-space. Thus equation (3) is the image formation.

The inverse transform can be worked out using the Fourier angular components of f and C 2 f :

             f (r, θ) = l f l (r)e ilθ with f l (r) = 1 2π 2π 0 f (r, θ)e -ilθ dθ and              C 2 f (ϕ, τ ) = l C 2 f l (τ )e ilϕ with C 2 f l (τ ) = 1 2π 2π 0 C 2 f (ϕ, τ )e -ilϕ dϕ (4)
Following [7] we give a new integral equation linking C 2 f l (τ ) to f l (r), the circular components of C 2 f (ϕ, τ ) and f (r, θ).

Since γ = θϕ and accounting for the invariance of the integrand under γ ←→ -γ, equation (3) takes the form

C 2 f l (τ ) = 2 π 2 0 r(γ) √ 1 + τ 2 1 + τ 2 cos 2 (γ) f l (γ) cos(|l|γ)dγ .
(5)

The absolute value of l is due to the invariance :

f -l (r) = f * l (r).
Because of the equation of the circular arc (equation ( 2)), we can show that dγ r(cos γ)

1 + τ 2 cos 2 γ = dr τ 2 -1 4 p r -r p 2 . ( 6
)
Now using this relationship between the differential elements, we change back to the r-variable equation (5) to obtain:

τ C 2 f l (τ ) √ 1 + τ 2 = 2 p p( √ 1+τ 2 -τ ) dr 1 -1 4τ 2 p r -r p 2 f l (r) × cos |l| cos -1 1 2τ p r - r p . (7)
This new form is adapted to Cormack's inversion procedure [3]. In the term cosh(l cosh -1 (h(r, τ ))), the function h(r, τ ) is a product of functions of r and of τ , i.e h(r, τ ) = h1(r) × h2(τ ).

We can choose the following change of variables [7],

q = 1 τ = tan ω, and t -1 = 1 2 p r - r p . (8)
Then we apply the Cormack's procedure [3] to invert equation (7) and obtain the inverse transform of the CART through the circular harmonic decomposition :

f l (r) = (-) 2p(p 2 + r 2 ) π(p 2 -r 2 ) 2 d dt ∞ t cosh(|l| cosh -1 ( q t )) q ( q t ) 2 -1 C 2 f l ( 1 q ) 1 + q 2 dq t= 2pr p 2 -r 2 . (9)
An alternative form of this result can be obtained with a further change of variables in the integration. Equation ( 9) may be now recast as

f l (r) = (-) 2p(p 2 + r 2 ) π(p 2 -r 2 ) 2 ∞ t cosh(|l| cosh -1 ( q t )) q 2 -t 2 d dq C 2 f l ( 1 q ) 1 + q 2 dq t= 2pr p 2 -r 2 . ( 10
)
Finally f (r, θ) is reconstructed through its Fourier expansion with the circular harmonic components f l (r). We have established the forward and inverse transforms in the case of this new Compton scattering tomography. In the next section, we will discuss how to use these transforms for numerical image reconstruction.

Norton's CST and Circular Radon transform

In 1994, S.J. Norton [6] worked out a CST modality which is based on a Radon transform on circles having a fixed common point. The working principle is given by Fig. 2. A point source S emits primary radiation towards an object defined by its electron density function ne(r, θ), of which M is a scattering site (running point).

A point detector D moves along an Ox -axis and collects, at given energy E, scattered radiation from the object. The physics of Compton scattering demands that the registered radiation flux density g at site D is due to the contribution of all scattering sites M lying on an arc of circle from S to D subtending an angle (πω), where ω is the scattering angle corresponding to the outgoing energy E, as given by the Compton formula.

Ω

S D ρ M (r 2 ,θ 2 ) 2 Source ω ω Detector ϕ θ2 θ1 M (r 1 ,θ 1 ) 1 r1 r2 Object Figure 2: Principle of Norton's CST
Norton gave the expression of the projections g as :

g(ρ, ϕ) = π 0 dθ ∞ 0 dr ρ ne(r, θ) w(r, θ; ρ, ϕ) × δ [r -2ρ cos(θ -ϕ)] (11)
where δ(.) is the 1-D Dirac delta function and w(.) is defined by:

w(r, θ; ρ, ϕ) = a r s(θ) P (ϕ) 4π (2ρ) 3 sin 2 θ . (12)
In the above equation, a represents the area of an element of detection, s(θ) expresses any angular dependance of the γray source distribution, and P (ω) (where ω = π/2+ϕ) is the Klein-Nishina differential cross section. Mathematically, g is essentially the Radon transform of the object electron density ne(M ) on arcs of circle, when radiation attenuation and photometric effects on radiation propagation are neglected.

Norton proposed an inverse formula given by:

ne(r, θ) = 1 π 2 2π 0 dϕ ∞ 0 dρ g(ρ, ϕ) w(r, θ; ρ, ϕ) × h [r -2ρ cos(θ -ϕ)] (13)
where

h(x) = ∞ -∞ e -iζx |ζ| dζ . ( 14
)
This expression is the same convolution kernel employed in the filtered Back-Projection algorithm used in x-ray transmission CT. The difference is that the Back-Projection is performed along straight lines in transmission CT, whereas here the Back-Projection is performed around the circles r = 2ρ cos(θϕ).

Now let :

     C 1 f (ρ, ϕ) = g(ρ, ϕ) (2ρ) 3 P (ϕ) f (r, θ) = ne(r, θ) as(θ)r 4π sin 2 θ p = 2ρ (15)
Substituting equation ( 15) in ( 11) we obtain:

C 1 f (p, ϕ) = π 0 dθ +∞ 0 dr p f (r, θ) δ {r -p cos(θ -ϕ)} . (16)
From ( 16) we use the Fourier angular components of f and C 1 f (like previously) and we can show that :

C 1 f l (p) = 2 p 0 dr f l (r) cos |l| cos -1 r p 1 -r p 2 dr . (17)
Equation ( 17) is precisely the integral equation for f l in the case of a (β = 1)-curve of Cormack [3,4]. Thanks to a property between cos and cosh, he could derive the following inverse formula for this circular Radon transform

f l (r) = 1 πr r 0 cosh(|l| cosh -1 ( r p )) ( r p ) 2 -1 d dp C 1 f l (p) dp . ( 18
)
Finally f (r, θ) is reconstructed through its Fourier expansion with the circular harmonic components f l (r). We have thus established the forward and inverse transforms (eq. ( 17) and ( 18)) in the case of the Norton's Compton scattering tomography.

In the next section, we will discuss how to use these transforms for numerical image reconstruction and how to derive novel algorithms based on circular harmonic decomposition.

NUMERICAL INVERSIONS

In this part, we will work out the numerical inversions that we propose for the both Circular Radon transforms presented above.

Computation of inverse CART

In principle, one can use equation [START_REF] Nguyen | Novel approach to stationary transmission scanning using Compton scattered radiation[END_REF] to perform numerical computations. However, assuming that the original function is bounded, a close inspection of the integral kernel of equation [START_REF] Nguyen | Novel approach to stationary transmission scanning using Compton scattered radiation[END_REF] shows that it behaves as

lim q→+∞ cosh(|l| cosh -1 ( q t )) q ( q t ) 2 -1 ≈ lim u→+∞ 2 e (|l|-2)u , (19)
where u = cosh -1 ( q t ). Equation ( 19) presents an apparent divergence when |l| > 2, since for q → ∞ the integrand grows very rapidly. This explains that when |l| > 2, this apparent divergence makes the q-integral unstable, and as such the presence of noise in the data C 2 f l (1/q) for large q is badly propagated into the calculation of f l (r) which prevents simulation studies. It becomes then obvious that equation [START_REF] Nguyen | Novel approach to stationary transmission scanning using Compton scattered radiation[END_REF] needs to be regularized.

Even if the circular arc doesn't belong to the α-curves family (defined by A.M. Cormack [3]), we will see that there is a frequential link between them. Indeed by introducing :

   F l (t = 2pr p 2 -r 2 ) = (p 2 -r 2 ) 2 2p(p 2 +r 2 ) × f l (r) G l (q) = C 2 f l ( 1 q ) √ 1+q 2 .
(

) 20
we obtain :

F l (t) = (-) 1 π ∞ t dq cosh(|l| cosh -1 (q/t)) q 2 -t 2 dG l (q) dq t= 2pr p 2 -r 2 (21)
Equation ( 21) is precisely the straight-line Radon transform inversion formula (α = 1) given by Cormack's regularization procedure. Thus with three changes of function, we pass from the original space of CART definition to a new space in which the CART becomes the ordinary Radon transform. Cormack proposed a regularization procedure based on consistency conditions of the circular harmonic components of the data ∞ 0 dq dG l (q) dq q α l = 0 for α l = (l -1), (l -3), ..., > 0 .

(

) 22
Thanks to these conditions and a property of the Tchebychev polynomial linking the polynomial of first kind (T l (.)) to the second kind one (U l (.))

T l (x) √ x 2 -1 = x - √ x 2 -1 l √ x 2 -1 + U l-1 (x) , (23)
he could derive the following regularised inverse formula for the CART:

F l (t) = 1 πt t 0 dq G l (q) U |l|-1 q t - 1 πt ∞ t dq G l (q) (q/t) -(q/t) 2 -1 |l| (q/t) 2 -1 . (24)
Then we follow the Chapman and Cary computional approach in the space where the Circular-Arc Radon transform becomes the standard Radon transform. In this algorithm we will consider that (q, t) ∈]0; Q] 2 with the same sampling ∆t and where Q = tan(max{ω}). Now let q = t cos χ, (respectively q = t cosh χ) in the first (respectively the second) integral of equation ( 24). Thus :

F l (t) = 1 π π 2 0 dχ G l (t cos χ) sin(lχ) - 1 π cosh -1 (Q/t) 0 dχ G l (t cosh χ) e -lχ (25)
Now we can define the discretized forms of F l (t) and G l (q) as :

G lk = G l (k∆t) F lj = F l (j∆t), (26) where (j
, k) ∈ [0, K] 2 with K = Q ∆t .
The derivative G l (q) can be approximated by coefficients a lkm where:

a lk = G l(k+1) -G lk ∆t . (27)
Equation (25) becomes so :

F lj = 1 π j-1 k=0 a lk (I l (χj(k + 1)) -I l (χjk)) + K-1 k=j a lk (J l (χj(k + 1)) -J l (χjk))   (28)
in terms of the primitive

I l (χ jk ) = χ jk sin(|l|x)dx = - cos(|l|χ jk ) |l| for l = 0 (29)
and

J l (χ jk ) = χ jk e -|l|x dx = -e -|l|χ jk /|l| if l = 0 χ jk if l = 0 ( 30
)
where

χ jk =    cos -1 k j for 0 ≤ k ≤ j cosh -1 k j for j ≤ k ≤ K. (31)
Finally the discretized form of the reconstruction equation appears as

F lj = 1 |l|π j-1 k=0 a lk cos |l|χ j(k+1) -cos |l|χ jk + 1 |l|π K-1 k=j a lk e -|l|χ j(k+1) -e -|l|χ jk (32)
And for l = 0 a separate expression exists

F0j = - 1 π K-1 k=j a 0k χ j(k+1) -χ jk . (33)
So having obtained all the F lj , which are the discretized circular components of the associated Radon transform, we recover the theoretical circular components of our circulararc problem f l (r) by using formula (20) given by

f l (r) = 2p(p 2 + r 2 ) (p 2 -r 2 ) 2 F l 2pr p 2 -r 2 , (34)
The final step consists in working out the summation in discrete form.

f (r, θ) = l f l (r) e ilθ . ( 35
)
We have established an alternative algorithm of image reconstruction based on circular harmonic decomposition for this new Compton scattering tomography modality.

Computation of inverse CRT

To establish our numerical inversion in the case of the Circular Radon transform presented above, we start from its analytical inverse:

f l (r) = 1 πr r 0 cosh(|l| cosh -1 ( r p )) ( r p ) 2 -1 d dp C 1 f l (p) dp (36)
This formula has the same kind of singularities that in the Circular-Arc Radon transform case since this transform belongs to the family of β-curves. So we can use the same method to regularize it and obtain :

f l (r) = 1 πr r 0 r p -( r p ) 2 -1 |l| ( r p ) 2 -1 C 1 f l (p)dp - 1 πr +∞ r U |l|-1 (r/p) C 1 f l (p)dp (37)
Then we will adapt the Chapman and Cary approach to this case. Indeed Chapman and Cary have studied the case of the (α = 1)-curve and so the case of the standard Radon transform. Now from the same idea we will propose a computional approach to reconstruct the original function f (r, θ) from the data C 1 f (p, ϕ). As we cannot consider an infinite size, we shall define the maximum value of p, pmax, and a minimum value, pmin, because a scattering site cannot be the source. Since,

cosh -1 (x) = ln(x + x 2 -1) = -ln(x -x 2 -1) (38)
we obtain

f l (r) = 1 πr r p min e -|l| cosh -1 (r/p) ( r p ) 2 -1 C 1 f l (p)dp - 1 πr pmax r U |l|-1 (r/p) C 1 f l (p)dp (39)
We make the change of variable r = p cosh(χ) in the first integral and r = p cos(χ) in the second one. The transform (39) becomes :

f l (r) = 1 π cosh -1 r p min 0 e -|l|χ cosh 2 χ C 1 f l r cosh χ dχ - 1 π cos -1 r pmax 0 sin(|l|χ) cos 2 χ C 1 f l r cos χ dχ (40)
To simplify the algebra and notation we consider that p and r have the same sampling ∆r. So we can define the discretized forms of f l (r) and C 1 f l (p) as :

g lk = C 1 f l (k∆r) and f lj = f l (j∆r) (41)
where (j, k) ∈ [1, K] 2 with K = ∆r pmax . Moreover we use linear interpolation to simplify the algorithm and the calculation of the derivative of the data g l . So in the interval k∆r < p < (k + 1)∆r, we can write:

ḡlk = g l(k+1) -g lk ∆r . ( 42
)
Finally we obtain a discretized form of our transform (eq (40))

f lj = - 1 π j-1 k=1 ḡlk J l (χ j(k+1) ) -J l (χ jk ) + K-1 k=j ḡlk I l (χ j(k+1) ) -I l (χ jk )   (43)
in terms of the primitive integrals :

I l (χ jk ) = χ jk sin(|l|x) cos 2 x dx and J l (χ jk ) = χ jk e -|l|x cosh 2 x dx . ( 44
)
The variables χ jk corresponds to the discrete radii p = k∆r

χ jk = cosh -1 j k for 1 ≤ k ≤ j cos -1 j k for j ≤ k ≤ K (45)
We evaluate the primitive integrals (44) by recurrence :

In(x) = 2 n -2 (tan x sin((n -2)x) -cos((n -2)x)) - n n -2 In-2(x) (46)
with the initial conditions :

I1(x) = 1 cos x and I2(x) = -2 ln(cos x) (47)
In the same way, we make the change of variable u = e χ jk in the primitive and we let n = |l| -1:

Jn(u) = du u n (u 2 + 1) 2 . ( 48
)
So we obtain the following recurrence relation :

Jn+1(u) = - 1 nu n (u 2 + 1) - n + 2 n Jn-1(u) (49)
with the initial conditions :

J0(u) = u 2(u 2 + 1) + 1 2 arctan(u) and
J1(u) = - 1 2 ln(u 2 + 1) + ln(u) + 1 2(u 2 + 1) (50)
We have also the special case l = 0. In this case, we evaluate the primitive integral J0(χ jk )

J0(χ jk ) = χ jk 1 cosh 2 x dx = sinh χ jk cosh χ jk = 1 - k j 2 (51) hence, f0j = - 1 π j-1 k=1 ḡ0k   1 - k + 1 j 2 -1 - k j 2   . ( 52
)
So having obtained all the f lj , which are the discretized circular components of our circular problem by constructing the theoretical f l (r). Finally we work out the summation in discrete form to obtain the original function.

f (r, θ) = l f l (r) e ilθ . ( 53
)
We have established an alternative algorithm of image reconstruction based on circular harmonic decomposition for the Norton's Compton scattering tomography. Moreover this algorithm is an alternative to which proposed Norton.

In the last section, we present the numerical results in order to prove the efficiency of our algorithms.

SIMULATION RESULTS

In this section, we present the numerical simulations of both circular Radon transforms (CART-CHD and CRT-CHD) and compare the reconstructions obtained by the ordinary Radon transform filtered back-projection (RT FBP) in the case of C 2 and by the Norton's inverse formula, which can be apparented to a "filtered back-projection" on circles (NRT FBP), in the case of C 1 .

Moreover the main aim of this work is to establish a new imaging principle by Compton scattered radiation. This is why other factors in realistic imaging systems are not treated here such as medium attenuation, inhomogeneous electron density and Poisson emission noise.

To illustrate and compare the quality of the reconstructions, we define the normalized mean square error (N M SE) and the normalized mean absolute error(N M AE):

N M SE = Ir -Io 2 2 N 2 , (54)
N M AE = Ir -Io 1 N 2 , ( 55
)
where Ir is the reconstructed image, Io is the original image (ground truth) and N 2 is the size of the studied image.

Image reconstructions by CART-CHD and RT-FBP

As an illstration of the feasibility of our new algorithm, we present numerical simulations applied on the medical phantom (Fig. 3). The scattering medium is discretized with 512 × 512 of length units (pixels). We consider the number of rotational positions Nϕ and the number of energy levels Nω. These numbers define the corresponding angular sampling steps dϕ and dω by :

dϕ = 2 π Nϕ and dω = π 2Nω . (56)
In order to have a "well-conditioned" problem, the number of projections (Nϕ × Nω) must be larger than the number of image pixels (512 × 512). This is why we take Nϕ = 720 and Nω = 500. Moreover we take the constant term p = 512 and min{ω} = dω.

Let us recall the classical Radon transform which is defined as integral of object function on straight lines. The forward Radon transform is:

g(u, ϕ) = R 2 dxdy f (x, y) δ (u -x cos ϕ -y sin ϕ) , (57)
and its inverse transform is:

f (x, y) = π 0 dϕ +∞ -∞ du +∞ -∞ dν |ν| × e -2iπν(u-x cos ϕ-y sin ϕ) g(u, ϕ) = - 1 2π 2 π 0 dϕ P.V. +∞ -∞ du ∂g(u, ϕ) ∂u 1 u -x cos ϕ -y sin ϕ . (58)
Equation ( 58) is called the filtered backprojection method (FBP). In this case, FBP is an exact inversion formula obtained by combining the action of the ramp filter (|ν|)and the backprojection operation of the Radon transform. This is the most popular inversion method for the ordinary Radon transform owing to its rapid algorithmic implementation.

The CART-CHD approach gives in general a reasonable image quality in the reconstruction of the medical phantom (Fig. 8). Indeed, contours and small structures are well recovered. Moreover, we can observe that there are a lot fewer artifacts compared to FBP reconstruction (Fig. 7) around the studied object. This is due to the fact that the CART-CHD reconstruction is consistent with the Radon data. The numerical error measurements (Table 1) obtained using the CART-CHD are very close to those of the ordinary Radon transform and proves the feasibility and the interest of this algorithm.

Image reconstructions by CRT-CHD and by NRT-FBP

Now we will present the numerical simulations in the case of our algorithm of the circular Radon transform applied on the Shepp-Logan phantom (Fig. 4). But this case we put the phantom farther from the source in order to avoid artifacts generated by the 1/r factor. Moreover the source is placed below on the left of the image and the detector moves along the axis Sx. The scattering medium is discretized with 256×256 of length units (pixels). We consider the number of rotational positions Nϕ and the number of radii Np. These numbers define the corresponding sampling steps dϕ and dp by :

dϕ = 2 π Nϕ and dp = 4 * 256 Np . ( 59
)
We take Nϕ = Np = 4 * 256. As the system of the corresponding CST is to not rotate around the object, the (ϕ, p)space is very large in front of the studied image. So to "well-observe" the object we have to take a large maximum value of p (4 times the image size).

The CRT-CHD approach gives a very interesting image quality in the reconstruction of the Shepp-Logan phantom (Fig. 11). Indeed, contours and small structures are better recovered than with the Norton's inverse formula (Fig. 10).

Nevertheless we have more artifacts generated around the studied object than in the Norton's case. This is due to the fact that the p ∈ R + and that we have to fix a maximum value pmax. Therefore this loss of information generated artifacts. To reduce this artifacts, we have to increase pmax and so to increase the length of the detector in the related CST. The numerical error measurements (Table 1) obtained using the CRT-CHD are smaller as compared to those of the Norton's inverse formula and prove the feasibility of this algorithm, particularly in the domain of the non-destructive testing. Table 1: "NMSE and NMAE for various reconstructions of the medical phantoms (Figs. 3 and4)"

Method

CONCLUSION

Alternative algorithms for the numerical inversion of two circular Radon transforms are established. Based on the circular harmonic decomposition, the advantage of these algorithms is to produce consistent image with the data (forward transforms) as opposed to analytical conventional methods like the "filtered back-projection" algorithm and to have the same computational complexity than the "filtered back-projection" method.

In the circular-arc Radon transform case, the simulation results prove the feasibility of the associated Compton scattering tomography which is suitable for biomedical imaging or non-destructive testing for example could present an alternative to X-ray tomography. In the second case we have compared our algorithm to the Norton's inverse formula.

Even if it seems hard to choose the best reconstruction way, we can say that both approaches are complementary since Norton's method gives a general reconstruction whereas ours reconstructs with more details of the small structures and the boundaries. The associated Compton scattering tomography is more suitable for the industrial nondestructive evaluation since we can reconstruct the original medium from one side of the object and without rotation.

Figure 1 :

 1 Figure 1: Principle of the new CST based on the CAR transform.

  × 10 -5 NRT-FBP Fig. 10 0.029 7.48 × 10 -4 CRT-CHD Fig. 11 0.028 5.9 × 10 -4

Figure 3 :Figure 4 :

 34 Figure 3: Illustration of a medical phantom.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 : 4 xFigure 10 :

 56789410 Figure 5: Radon transform of the medical phantom shown in figure 3

Figure 11 :

 11 Figure 11: Reconstruction of the Shepp-Logan phantom shown in figure 4 using CRT-CHD and data in figure 9.