J Marco

Marco Mendes
email: marco.mendes@fe.up.pt

Francisco Restivo

Paulo Leitão
email: pleitao@ipb.pt

Armando W Colombo
email: armando.colombo@de.schneider-electric.com

J Marco Mendes

Petri Net Based Engineering and Software Methodology for Service-Oriented Industrial Automation

Keywords: Service-oriented Architecture, Service-oriented Computing, Industrial Automation, Software engineering, Petri Nets

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

Introduction

In Service-oriented Architectures (SoA) for industrial automation [START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF] several research directions and solutions have been presented, but overwhelmingly directed to a specific part of the whole engineering problem that SoA reflects in industrial automation. Thus, an inroad is required for integrated solutions of engineering, respecting users, developers, available hardware and software. Petri nets (PN) were identified as being part of the solution by presenting a set of useful characteristics supporting the life-cycle of service-oriented systems. With the choice over Petri net based formalism, the research question and motivation of this work follows as: "How does Petri Nets based methodology's features contribute as a unified tool to the integration and collaboration of autonomous service-oriented automation control components and support their life-cycle management?" This work focuses on the actual results in terms of research, development and application of the integrated methodology for the engineering of service-oriented automation systems using Petri nets. This paper is organized in the following sections: Section 2 describes the contribution to technological innovation, section 3 presents an overview on related literature and state of the art, section 4 explains the overall engineering approach and software for service-oriented automation systems using as core methodology Petri nets, the application is resumed and discussed in section 5 and finally conclusions and future work are exposed in section 6.

Contribution to Technological Innovation

The present research work by using Petri nets for the engineering of service-oriented automation systems offer several novel aspects that already were reviewed by international committees (see e.g. the SOCRADES project at http://www.socrades.eu) and recent disseminations. Note that the aspects discussed here are part of the individual research work used in the context of the project. From the SoA side, it is well known that it represent a future trend in the automation and production world, but from the software engineering and also end-user perspective, integrated solutions are missing. As such, the research work focus on a Petri net based solution for those systems, for the design, analysis and also operation. Contributions are in a form of an open methodology for developing custom and feature-full Petri net applications and control solutions, an engineering process based on Petri nets and several runtime features such as distributed orchestration, composition and active conflict resolution.

Overview on Related Literature and State of the Art

In SoA, the main objective is the interoperability between different participants in the system by providing and using resources in form of services via an open protocol based communication platform (e.g. Web Services). In automation and manufacturing, service providers and requesters include also industrial devices, manufacturing equipment and products (besides the typical software agents that run on normal computers and servers). Several efforts are being done in integrating and managing these new systems. The first visible application resulted from the SIRENA project [START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF] with the objective to develop a service infrastructure for real-time embedded networked applications. Other projects in the same research area are currently in advance research and development state, such as the SOCRADES project.

Important subjects are still related to engineering solutions of service-oriented systems from the software/hardware and user point of view, and also to stimulate the industrial adoption. In service-oriented automation systems, a formal and unified method is required to provide design facilities, with the ability to validate models and to be also used as an integration middleware at runtime with enough flexibility and features. Previous works in the domain report some advances with engineering methodologies. Main research directions are related to SoA device integration and Device Profile for Web Services (DPWS) [START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF] [2] [START_REF] Li | Service-Oriented Embedded Device Model in Industrial Automation[END_REF], orchestration and semantics of devices and systems [START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF] [4], the use of Business Process Execution Language (BPEL) and other formalisms [START_REF] Puttonen | An application of BPEL for service orchestration in an industrial environment[END_REF] [6], process optimization [START_REF] Rahman | Supply chain optimization towards personalizing web services[END_REF], business integration (for example [START_REF] Nguyen | Exploiting SOA for adaptive and distributed manufacturing with cross enterprise shop floor commerce[END_REF], integration of SoA and Multi-agent Systems (MAS) [START_REF] Ribeiro | MAS and SOA: Complementary Automation Paradigms[END_REF] and virtual SoA enabled production environments [START_REF] Cachapa | An approach for integrating real and virtual production automation devices applying the service-oriented architecture paradigm[END_REF].

Petri net is a language for design, specification, simulation and verification of systems. Petri nets are a graphical and mathematical modeling tool applicable to any systems. They are a promising tool for describing and studying information processing systems that are characterized as being concurrent, asynchronous, distributed, parallel, nondeterministic, and/or stochastic [START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF]. Petri Nets can be applied to many areas, including the industrial automation and control. Since the formal Petri nets are very basic in the application point of view, numerous extensions nets have been introduced over the years, normally designated by high-level Petri net. Logic controllers of Petri nets have also been developed (see the sequence controllers developed in the early 80's by Hitachi Ltd. [START_REF] Murata | A Petri Net-Based Controller for Flexible and Maintainable Sequence Control and its Applications in Factory Automation[END_REF] and the application of fuzzy Petri net controllers [START_REF] Gomes | Petri net based Programmable Fuzzy Controller targeted for distributed control environments[END_REF]) or converted into other formalisms (such as the industrial standard languages of IEC 61131) that can be understood by the devices [START_REF] Uzam | Conversion of Petri net controllers for manufacturing systems into ladder logic diagrams[END_REF].

The traditional application of Petri nets in SoA environments seems to be orchestration and choreography, to define sequences, conditions, interactions and compositions of services. Petri nets and high-level Petri nets are used in orchestration and choreography for modeling processes/composition [START_REF] Hamadi | A Petri net-based model for web service composition[END_REF] [16], analysis purposes [START_REF] Yang | Verifying web services composition based on hierarchical colored petri nets[END_REF] [18] and negotiation [START_REF] Jiang | Modeling of Web-based collaborative negotiation systems using colored Petri net[END_REF]. However, the OASIS standard of Web Services Business Process Execution Language (WS-BPEL) is the most prominent standard for orchestration of services. The application of BPEL directly in automation environment can be discussed. From one side, it has already a well defined syntax in XML that can be used directly with Web services. From the other side, BPEL is a specification mainly targeting business requirements, it is unknown to automation system engineers, too complex and descriptive to be interpreted by resource constrained devices (typically used in automation), not suitable for internal service process description based on device/software capabilities, technology dependent on Web services and missing direct analysis and validation support. However, some efforts have been done in terms of the application of BPEL and orchestration in general for industrial automation with results (see for example [START_REF] Puttonen | An application of BPEL for service orchestration in an industrial environment[END_REF]).

Research Contribution and Innovation

Especially for service-oriented industrial automation, the application of Petri nets must be open to several specifications, requirements and methodologies for engineering that will possibly come in future. Considering the several aspects, an approach was specified for a concrete Petri nets methodology that is feasible and customized for the studied needs. The resulting open methodology for the definition of Petri net based applications permits the construction of custom software using a ground basis made of several elementary packages. This includes the Petri net formalism (according to [START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF]), analysis routines and conflict management. Besides the required openness, another requirement is the introduction of the time factor when evolving the Petri net. Consequently, the basis of the Petri nets is guaranteed, except that the doors are open in terms of delays that can be used for customized operations. Over this, a modular property system was specified and a novel approach for the creation of token games (an "engine" that runs a Petri net) based on a customized template. Please consult the documentation of [START_REF] Mendes | Customizable Service-oriented Petri Net Controllers[END_REF] for more information.

Having the basis for Petri nets applications and engineering, the next step would be identify what are the characteristics that can be extracted from them to be used in service-oriented automation systems. Most of these required or welcome features were already researched and documented. The remaining ones can be considered as future research work and possibly will also generate new characteristics. Table 1 lists the features and applications of Petri nets (and extensions) that are considered and also references to some publications with more detail over them. [START_REF] Mendes | Customizable Service-oriented Petri Net Controllers[END_REF] Application software tools for PC [START_REF] Mendes | Software Methodologies for the Engineering of Service-Oriented Industrial Automation: The Continuum Project[END_REF] Composition of Petri nets [START_REF] Mendes | High-Level Petri Nets Control Modules for Service-Oriented Devices: A Case Study[END_REF] Petri net controller for devices [START_REF] Mendes | Software Methodologies for the Engineering of Service-Oriented Industrial Automation: The Continuum Project[END_REF] Decision support system for PN [START_REF] Leitão | Decision Support System in a Service-oriented Control Architecture for Industrial Automation[END_REF] The use of Petri nets can be presented in the life-cycle of processes as modeling, analysis (simulation) and execution (control). For the development of the required software applications, there are several needs since an integrated engineering approach and resulting software. Based on the requirements of physical hardware and users, it was decided that the basic building blocks that compose the distributed system should be configurable software components assuming different tasks. Therefore, the software components were designated as "bots" (that have the so called "orchestration engine" embedded inside) and are able (in a service-oriented fashion) coordinate their activity and proceed also to collaboration processes with other components in the system. The central part of Fig. 1 has represented the concept, where the domain of autonomy of the bot is, for example, a transfer unit (conveyor). Communication to other software components in the system (e.g. other bots) is done via the exposition of services as well as the requesting of necessary services from others. To design, configure and maintain bots, there is a need of specific tools, that are user-friendly and speed-up the development, using a high-level programming approach (visual languages, such as Petri nets).

Since SoA-based automation systems are lacking in dedicated software applications and methodologies, effort where done in the development of several software packages for the demonstrator. In any case, system engineering and associated tools are required to facilitate the developer's intervention. From the Petri nets side, the practical usage is limited by the lack of computer tools which would allow handling large and complex nets in a comfortable way. The project was baptized Continuum Development Tools (CDT), named after the continuum concept used in physics and philosophy. First developments were started by integrating already developed software components, in special the old editor called PndK (Petri nets development toolKit), under the same umbrella. Along with the integration, it was identified that several software packages are needed, namely: a framework for developing bots, engineering tools for the design and managing of bots and several utilities (mainly libraries) for supporting activities (e.g. such as communication and interface for devices). The main component is the Continuum Bot Framework (CBF) for the development of bots and their functional modules, inspired in the anatomy of living beings. Another component, the Continuum Development Studio (CDS) that is based on an extensible Document/View framework, provides an engineering tool for service-oriented bots, for example, supporting the visual description, analysis and simulation of their behavior (for now, in Petri nets formalism). The utilities package includes several reused software libraries and tools, some developed internally others adopted from the outside, such as the SOA4D DPWS library (available at https://forge.soa4d.org) providing facilities for the development of Web services and the Qt toolkit (see http://qt.nokia.com/products/), used mainly as a graphical toolkit for human interaction in the CDS. Fig. 1. Concept of the automation bot and engineering process (adapted from [START_REF] Mendes | Software Methodologies for the Engineering of Service-Oriented Industrial Automation: The Continuum Project[END_REF]).

The engineering is overviewed in the external boxes of Fig. 1. The CDS is employed for designing and analyzing the Petri nets template models for describing the behavior of the bots (Fig. 1.a). When importing the device/connection information, several steps are done by the CDS: 1) instantiate Petri nets models for each bot based on the designed template models, 2) create the necessary properties of the Petri nets models so that several parameters of the given device/connection information and 3) based on the "enrichment" of the Petri nets models, composition of models can be done for creating connection logic and for the overall system analysis (Fig. 1.b). After the analysis and simulation (that can be done with the CDS), bots must now be configured. The process of deploying a service that encapsulates its logic as a Petri nets model to a bot that provides an embedded Petri Nets Kernel Module is depicted in Fig. 1.c. The deployment functionality is a standard feature of the DPWS and is exposed as a dynamic deployment service. The target and the deployment service can be discovered by stacks built-in discovery service. After deployment, a new service endpoint has been added and the execution of the services logic has been initiated. Deployment information includes the Petri nets behavior model, connection information of neighbors (required services), provided services by the bot and also extra configuration information for the other modules of the bot. The bot will configure itself (and its modules) and is then ready for operation (Fig. 1.d).

Application and Discussion

The implementation and evaluation of the current methodology is done using a virtual and real hardware demonstrator. The virtual one was modeled using the DELMIA solution to provide a simulation environment based on the real system. A service framework was added to DELMIA so that the environment can be controlled by the exposition and request of services, in the same way as pretended by the physical counterpart. The assembly equipment compromises a set of conveying modules, two lifters that make the interface between the two levels and two assembly workstations. The objective is to route the incoming pallets to the correct workstation, according to their production plan. The equipment is controlled by automation devices where the bots (including the Petri net engine, DPWS framework and other modules) run. These devices are then connected to the network, so that they can be configured, and their services can be used.

First results show the applicability and correct behavior of the system and reduced efforts in development. However and since it represents a first prototype application, there are some tasks of improvement on both methodology and application. One requirement is the previous knowledge about Petri nets, SoA and the engineering approach. This information has to be transmitted to the users and developers that are normally used to the traditional automation systems. The software part has to be improved in terms of stability, performance (e.g. improving response form embedded devices) and user-friendliness (mostly the engineering procedures using the CDS). Not less important is the automatic integration into higher level IT-systems, especially dealing with production orders and business needs.

This approach is used actively in the SOCRADES project for the specified demonstrator and has been evaluated by experts in the field. More information can be seen on the public available video at YouTube ("IP SOCRADES Demonstration of Service-Oriented Architecture integrating Real and Virtual Devices in the Electronic Assembly Scenario") that shows some of the application results of this approach.

Conclusions and Further Work

The current contributions of this work can be resumed by the specification of unified service-oriented engineering method that can be used in a collaborative automation environment and prove features of Petri net-based solution in service-oriented automation and production systems as a unified methodology and integration middleware. In the end, it should contribute to the reduction of the design, operational and reconfiguration phases in service-oriented automation systems. New features of this work can be addressed by the combination of the methodology itself, especially the token game template for customized implementations, conflict management, property system and also some aspects of engineering in service-oriented industrial automation. Worth of mention is that it concludes also in a new form of engineering that is different form the traditional used in automation systems.

Table 1 .

 1 Features and applications of Petri nets and reference publications.

	Features and applications	Ref.	Features and applications	Ref.
	Open methods for PN token games	[20]	Migration and virtual environments	[23]
	Property system	[20]	Knowledge extraction from PN	-
	Analysis and validation	[21]	PN interchange format (XML)	-
	Petri net modularization/interfaces	[21]	Petri net-based engineering	[24]
	Petri net service association

Acknowledgements. The authors would like to thank the European Commission and the partners of the EU IST FP6 project "Service-Oriented Cross-layer infrastructure for Distributed smart Embedded devices" (SOCRADES), the EU FP6 "Network of Excellence for Innovative Production Machines and Systems" (I*PROMS), and the European ICT FP7 project "Cooperating Objects Network of Excellence" (CONET) for their support.


Miguel Pereira

Luís Bernardo

Rui Dinis

Rodolfo Oliveira

Paulo Carvalho

Paulo Pinto

Delay Analysis for TDMA Schemes with Packet Recombining

Keywords: Delay analysis, ARQ schemes, Mac protocol and Markov Chains

This paper considers the use of SC modulations (Single-Carrier) with FDE (Frequency-Domain Equalization) with low-complexity soft combining ARQ schemes (Automatic Repeat reQuest). With our technique, packets associated to different transmission attempts are combined in a soft way, allowing improved performances. Its low complexity makes it particularly interesting for the uplink of wireless systems. This paper proposes an accurate analytical model for a TDMA (Time Division Multiple Access) scheme where packet combining ARQ is applied. It evaluates the uplink non-saturated packet delay for a generic message arrival process. Our analytical results are validated using physical and MAC layer simulations.

Introduction

Packet errors can be significant in wireless systems due to fading and shadowing effects, which may lead to significant decrease in the receiving power. The traditional approach to cope with an erroneous packet is to discard it, and ask for its retransmission, which corresponds to the conventional ARQ techniques (Automatic Repeat reQuest). The major problem with conventional ARQ techniques is that persistent unfavorable propagation conditions may originate a very high PER (Packet Error Rate). The individual packet error probability is not affected by the number of packet retransmissions because the information contained in signal of the discarded packets is not used. Hybrid ARQ/FEC (Forward Error Correction) strategies were proposed [START_REF] Hagenauer | Rate-compatible Punctured Convolutional Codes (RPCP Codes) and Their Applications[END_REF][START_REF] Gusmão | Adaptive HARQ Schemes Using Punctured RR Codes for ATM-compatible Broadband Wireless Communications[END_REF] to cope with these scenarios. They retain the signal associated to an erroneous packet and they may ask for additional redundancy. The traditional H-ARQ (Hybrid ARQ/FEC) approach relies on a powerful error correcting code that is punctured, to increase the data rate. Successive retransmissions reduce the puncturing, and therefore increase the error correction capacities of the code. This limits the error correction capability to the base code capabilities. We proposed on [START_REF] Dinis | Soft Combining ARQ Techniques for Wireless Systems Employing SC-FDE Schemes[END_REF] the use of soft packet combining techniques to implement H-ARQ on SC modulations (Single Carrier) with FDE (Frequency-Domain Equalization). It is based on repetition codes with soft decision, which is not bounded by the performance of the basic code. The use of H-ARQ techniques improve the network throughput compared to conventional ARQ technique because the packet error probability after retransmission i ( i p ) is usually lower than 1 i p -. It also reduces the average packet delay. In this paper we study the performance of H-ARQ techniques for SC-FDE schemes. SC-FDE schemes are generally accepted as one of the best candidates for the uplink of future broadband wireless systems [START_REF] Falconer | Frequency Domain Equalization for Single-Carrier Broadband Wireless Communications[END_REF]. We consider a TDMA (Time Division Multiple Access) scheme, where a station shares a multiple access communications channel by transmitting its messages during its dedicated time slots, and a generic packet arrival distribution function.

TDMA is used in several current wireless network systems (802.16 , 802.11s, etc.) when hard QoS (Quality of Service) guarantees are needed. TDMA was the focus of an extensive list of past works (e.g. [START_REF] Saeki | An Analysis of a TDMA Channel Using Stop-and-Wait, Block, and Select-and Repeat ARQ Error Control[END_REF][START_REF] Rubin | Message Delay Analysis for TDMA Schemes Using Continous-Slot Assignments[END_REF][START_REF] Chen | Per connection delay analysis of a frame-based TDMA/CDMA MAC protocol[END_REF]), which assumed a constant packet error probability for the different retransmissions. ARQ error control was analyzed in [START_REF] Saeki | An Analysis of a TDMA Channel Using Stop-and-Wait, Block, and Select-and Repeat ARQ Error Control[END_REF] for a Poisson distributed message arrival process. Rubin and Zhang [START_REF] Rubin | Message Delay Analysis for TDMA Schemes Using Continous-Slot Assignments[END_REF] analyzed the message delay for TDMA using multiple contiguous-slot assignments, considering a generic message arrival process. Chen and Chang [START_REF] Chen | Per connection delay analysis of a frame-based TDMA/CDMA MAC protocol[END_REF] proposed a delay analysis model for TDMA model where slot assignment is characterized by a random variable, for a Poisson distributed message arrival process.

Contribution to Technological Innovation

This paper is the first to propose an exact analytical model concerning the queue-size and the message delay analysis of the TDMA system with an H-ARQ technique, for a generic message arrival distribution. This model is validated for Poisson traffic using ns-2, showing a significant packet delay reduction compared to the conventional ARQ technique.

System Overview

In this paper we consider the uplink transmission in an SC-FDE system. The timedomain block associated to a given user (i.e., the corresponding packet) is

{ ; = 0,1, , 1} n a n N - K
, where n a is selected from a given constellation and N is the DFT size (Discrete Fourier Transform). As with other block transmission techniques, a suitable cyclic prefix is added to each time-domain block.

If we detect errors in the packet we ask for its retransmission, but we store the signal associated to each transmission attempt. Although we could keep trying to transmit the packet until there were no errors, in practice there is a maximum number of transmission attempts N R . If we fail after N R attempts we need to change the transmission parameters (transmit power, carrier frequency, base station, etc.) since the channel is too bad.

The packet associated to r th attempt to transmit { ; = 0,1, , 1} The received signal associated to the r th transmission attempt is sampled and the cyclic prefix is removed, leading to the time-domain block ( ) { ; = 0,1, , 1}

n a n N - K is ( ) { ; = 0,1, , 1}
r n y n N -K . If the cyclic prefix is longer than the overall channel impulse response then the corresponding frequency-domain block is ( )

{ ; = 0,1, , 1} r k Y k N - K , where ( ) ( ) ( ) ( ) ( ) = , r r r r r k k k k k Y A H N I + + (1) with (
{ ; = 0,1, , 1} r k A k N - K is the DFT of ( ) { ; = 0,1, , 1} r n a n N - K and ( ) r k H
is the overall channel frequency response for the r th transmission attempt.

Let us assume that we have R versions of the packet (i.e., there were R transmission attempts). Our receiver, which is based on the IB-DFE receivers proposed in [START_REF] Benvenuto | Block Iterative DFE for Single Carrier Modulation[END_REF][START_REF] Dinis | On Broadband Block Transmission over Strongly Frequency-Selective Fading Channels[END_REF], is depicted in fig 1 . We have an iterative frequency-domain receiver where, for a given iteration i , the frequency-domain samples at the output are given by

( 1) ( ) ( , ) ( ) ( ) =1 = R i i r i r i k k k k k r A F Y B A - - ∑ % (2) where ( , )
{ ; = 0,1, , 1}

r i k F k N - K ( = 1,2, , r R K ) and ( ) { ; = 0,1, , 1} i k B k N - K can be
regarded as the feedforward and the feedback coefficients, respectively. , where n a denotes the average symbol values conditioned to the FDE output, which can be obtain as described in [START_REF] Gusmão | A Class of Iterative FDE Techniques for Reduced-CP SC-Based Block Transmission[END_REF].

The optimum feedforward coefficients for a given iteration, can be written as (see

[3]), ( , ) ( , ) ( ) = , r i r i k i k F F γ ( with 1 ( ) ( , ) ( ) = 0 =1 1 = N R i r i r k k k r F H N γ - ∑∑ ( and ( )* ( , ) ( 1) 2 ( ) 2 =1 = , (1 ( ) ) | | r r i k k R i r k r H F H α ρ ′ - ′ + - ∑ ( = 1,2, , , r R K (3) where ( ) 2 2 [| | ] = [| | ] r k k E N E A α (i.e.
, α is the inverse of the SNR (Signal-to-Noise Ratio)) and the correlation coefficient ( )

i ρ is given by 1 ( ) ( ) ( ) = 0 1 = (| { } | | { } |). 2 N i i i n n n Re a Im a N ρ - + ∑ (4)
The optimum feedback coefficients (also to minimize the signal-to-noise plus interference ratio) are given by

( ) ( , ) ( ) =1 = 1. R i r i r k k k r B F H ′ ′ ′ - ∑ (5)

The System Model

In this section we present the exact message delay and queue-size analysis for the system.

System Characterization

We consider a communication channel on a TDMA basis. Time is divided into equal length slots, each of duration equal to 1 block data of 32 bytes. The start of a message transmission across the channel must coincide with the beginning of a slot. The success of an arriving message is given by l q , 0 < < 1 l q , with l , 1 l R ≤ ≤ , representing the transmission number of each message and R the maximum number of retransmissions allowed for each message. Each packet contains 32 bytes, so that a packet transmission time is equal to 1 slot. n N represents the number of messages arrivals at the station during the n th slot, according to a stochastic process. We assume that { } , and characterizes the system dynamics for a generic N value.

The Generating Function of the Queue Size

Let k X represent the number of messages in the station's queue at the start of slot k (including message arrivals that will occur during slot 1 k ), and let k Z represent the number of transmissions of one message at the same slot. So, at the start of slot k we have k X messages in queue, and we will transmit the corresponding message for the k Z th time. The process {( , ) 1} )

, , = = | = , = lim j i l k k k k P P X i Y j Z l →∞ (6) ( ) ( ) =1 =1 =1 1 = = | = , = lim J N L R k k k J k j l i P X i Y j Z l J π + →∞ ∑∑∑ (7) Next, we define a set of z -transforms, | | 1 z ≤ . ( ) , , , =0 = i j l j i l i U z P z ∞ ∑ 1 , 1 j N L l R ≤ ≤ + ≤ ≤ (8) ( ) ( ) =0 = i i z i z π ∞ Π ∑ (9) ( ) ( ) ( ) 2 = 0 = 1 = , 1 = i i i A z a z A A α α ∞ ′ ′′ ∑ ( 10
)
The steady-state probabilities referred previously satisfy the following equilibrium equations. For

0 ≥ i 1, , , , = 0 = i i l m N L i m l m P a P + - ∑ , , 1, , = 0 = , 2 i j i l m j i m l m P a P N j N L -- + ≤ ≤ + ∑ (11) (12) 1 , ,1 1,0,1 1, 1 , 1, 1 , =1 = 0 = 0 = , 2 1
R i i j i i j l m j i m l m j i m R l m m P a P q a P a P j N - - -+ - -+ - + + ≤ ≤ + ∑ ∑ ∑ (13) ( ) 1 , , 1 1, , 1 =0 = 1 , 2 ,2 1 i j i l l m j i m l m P q a P l R j N - - -- - - ≤ ≤ ≤ ≤ + ∑ ( 14
)
Looking the fact that ( ) ( )

1 = = k P Y j N L - +
, for 1 j N L ≤ ≤ + , we have ( )

, , =1 =1 1 = N L R j i l j l i P N L π + + ∑∑ (15)
Considering ( 8) and ( 9) we then obtain ( )

, =1 =1 1 = N L R j l j l z U N L + Π + ∑∑ (16)
The objective is to derive an expression for ( ) z Π . Multiplying the i th equation of (11) to (14) by i z and summing over i , we obtain the following set of equations

( ) ( ) ( ) 1, , = l N L l U z A z U z + ( ) ( ) ( ) , 1, = , 2 j l j l U z A z U z N j N L - + ≤ ≤ + ( ) ( ) ( ) ( ) ( ) , 1 1, 1 1,0, 1 = 1 , 2 1 j l l j l j l U z q A z U z P j N - -- - - - - ≤ ≤ + (17) (18) (19)
        - + +         ∑ (22) ) ( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 0 1 =2 = 1 R L R L l L R l l l g z z q A z Q A z q Q A z - + + + -     - + +     ∑ (23) and ( ) 1 =1
= 1 X X m m Q q - - ∏ (24)
By the normalization condition, ( )

2, =1 1 = 1 R l l U ∑ , we then obtain, 1 1,0,1 =2 =2 = 1 R R z l l z l l D P Q Q L -         + -         ∑ ∑ (25) where ( ) 2,1 1 U
is determined applying the L'Hôspital rule and is equal to

( ) 1 2,1 1 1 = f U g (26)
and where

1 1,0,1 = z f P D (27) ( )( ) ( )( ) 1 1 =2 = 1 1 1 1 1 1 1 1 R z R l l l D q Q R L q Q l L α -         + -- - + + - - + +         ∑ (28)
( )

1 1 1 =2 = 1 1 R R l l l g L q Q R q Q l α -     -+ + +     ∑ (29)
From ( 16) we obtain the generating function fully transmitted from the system. For the general batch arrival case, this last component of the delay can be obtained using Little's formula assuming a first-comefirst-served (FCFS) service discipline. That is,

[ ] 1 [ ] = 2 E X E D α + (36)
where [ ] = (1) E X ′ Π is given by (33). This equation can be used with any message arrival distribution, as long as the message arrival process, ( ) A Z , is defined. The next section validates the proposed model, assuming arrivals to be governed by a Poisson process, so ( ) ( ) ( )

1 = , 1 = = . z A z e A λ α λ -- ′ (37)

Simulations

In this section, we present a set of performance results concerning the proposed packet combining ARQ scheme for a TDMA uplink channel. Each FFT (Fast Fourier Transform) block has N=256 data symbols and are selected from a QPSK constellation, with Gray mapping totalling 512 bits per data block. The data packets correspond to one FFT block, with 32 bytes. We consider an uncoded scenario where the channels for each packet retransmission remains fixed. To allow an efficient packet separation we consider the cyclic-shifted versions of the FFT-blocks for each retransmission attempt which formally is equivalent to have cyclic-shifted versions of the channel in each retransmission attempt [START_REF] Dinis | Soft Combining ARQ Techniques for Wireless Systems Employing SC-FDE Schemes[END_REF]. We assumed perfect synchronization and channel estimation conditions. We consider, as referred earlier, one slot per station ( ) = 1 N and 8 stations ( ) = 8 L generating traffic following a Poisson process with λ packets/slot. Fig. 1 shows the average PER of our interactive receiver when we have l transmission attempts ( l q ). It shows a significant PER reduction with a l 's increment, especially when the number of R is not too large (the ). This is not surprising, since the total transmitted power grows with l and we take full advantage of all transmitted power.

We implemented the proposed packet combining approach in the ns-2 simulator [12]. We assumed that l q for packets is a mean-ergodic process, and we modeled the PER variation with l using the average PER values measured in the first set of experiments, and represented in Fig. 1. The average l q in function of 0 / b E N is introduced as a matrix parameter for the simulations. q ≈ for 3dB), clearly showing its advantage compared to a classical TDMA scheme.

Conclusions

In this paper, we have carried out exact message delay and queue-size analysis for a single slot per frame packet-switched TDMA applying a packet combining ARQ scheme. We propose a generic analytical model that applies to any message arrival distribution. The proposed analysis was validated using simulations for a SC-FDE scheme, showing its correctness. The results show the significant gains that can be obtained using a packet combining ARQ scheme compared to a classical TDMA scheme.

Future work includes the testing of the model with other message arrival processes, including long-tail processes. At model level extended models are required to handle multiple slots per frame, and other queueing disciplines.

  performance in the presence of strong inband interference.

α, with 1 N, with 1 L

 11 is a sequence of i.i.d. random variables, and we set where α is the mean number of messages arrivals per slot. The station is allocated a possible multiple contiguous slots per frame N ≥ . A frame contains N L + slots≥ , representing the group of slots allocated to other stations. Messages can only be transmitted in those N dedicated slots during a frame. This section proposes a complete solution for the delay with = 1 N

  is not Markovian because it is time dependent. To make it Markovian we use one supplementary variable, k Y , to represent the slot index in the frame, 1 k k ≥ .It can be shown that the system is stable if the average number of messages arrivals per slot α is less than the average number of messages transmitted per slot. Under the stability condition, the conditional steady-state distribution { }

(

Fig. 2

 2 depicts the simulation and the analytical model message delay in function of 0 / b E N for different Poisson traffic load (when 1 / λ equal to 10, 20 and 30 packets), for 4 = R (it attempts up to 4 transmissions per packet). The analytical results were computed using 43 and 47. It shows that the measured delay follows the analytical model values, thus validating it. It also shows that message delay decreases as 0 / b E N increases. Notice that packet combining allows a finite delay even for a initial very high PER ( 1 0

Fig. 1 .

 1 Fig. 1. PER performance for a single packet transmission.

Fig. 2 .

 2 Fig. 2. Message delay with

( )

Henceforth this paper considers only the resolution of { } k X 's steady-state distribution for = 1

N

. Using (17) to (20) iteratively, and due to the fact of 1,0, = 0 l P for 2 l ≥ we obtain after some algebric manipulations

where

where

The Mean Queue Size

The mean the number of messages in the station's queue can be computed by taking the derivative of ( ) z Π and letting z equal to 1. From (30) we obtain

To compute ( )

we apply twice the L'Hôpital's rule to obtain

Delay Analysis

The message delay [ ] E D is defined as the time interval elapsed between the message arrival and departure times. It can be calculated as the sum of the time elapsed between the message arrival and the beginning of the next slot (equal to one half-slot) plus the number of slots elapsed between the slot of its arrival and the instant it is


André Mora

Pedro Vieira

José Fonseca

Advances in Image Processing Techniques for Drusens Detection and Quantification in Fundus Images

Keywords: image processing, retina imaging, automatic segmentation

Age-Related Macular Degeneration (ARMD) is considered the leading cause of irreversible blindness in developed countries. One of its risk factors is the presence of drusens, which are retina abnormalities appearing as yellowish spots in fundus images. In this article a methodology using image processing techniques for the quantification of drusens is presented. The method uses splines combined with a contrast normalization to correct uneven illumination, followed by a drusen detection and modelling algorithm. The detection uses a gradient based segmentation algorithm that isolates drusens. They are then fitted by Gaussian functions, producing a model that is used to compute the area affected. To validate the methodology, 22 images were marked by three ophthalmologists and compared to the automated method. The sensitivity and specificity for the automated process (0.664 and 0.963) were comparable to that obtained among the specialists (0.656 and 0.971). Also, the Intraclass Correlation Coefficient showed an agreement of 74.9% between the processed images and the specialists' analysis.

Introduction

Some ophthalmologic lesions have visible symptoms in the retina surface. Drusens, accumulations of extra-cellular materials beneath the retina surface, are one of those mentioned above. They are a risk factor for Age-Related Macular Degeneration, the leading cause of irreversible blindness in developed countries [START_REF]World Health Organization, International Agency for the Prevention of Blindness: State of World's Sight -Vision 2020: the right to sight[END_REF]. They appear in fundus images as speckles of higher intensity (see Fig. 1) and are qualitatively evaluated by the ophthalmologists by using fuzzy variables like: shape, size, location, etc. For this examination to be accurate, by quantifying the abnormalities, it requires a fastidious work from the expert that must be skilled for this task.

An accurate quantification of drusens in retina images using automatic procedures is one of the challenges for researchers who study medical imaging. This is a subjective analysis and its criteria can differ among specialists and among analysis performed by one specialist at different moments. If we also consider that the design of an automatic diagnosis system is complex, then the development of such accurate procedure is a difficult task.

Contribution to Technological Innovation

Several studies have been published in the last twenty-years for drusens detection using image enhancement techniques. They included image pre-processing for illumination compensation and contrast improvement, and for the drusens segmentation have been proposed algorithms using adaptive local thresholds [START_REF] Peli | Drusen Measurement from Fundus Photographs Using Computer Image Analysis[END_REF][START_REF] Kirkpatrick | Quantitative image analysis of macular drusen from fundus photographs and scanning laser ophthalmoscope images[END_REF][START_REF] Rapantzikos | Detection and segmentation of drusen deposits on human retina: potential in the diagnosis of age-related macular degeneration[END_REF][START_REF] Morgan | Automated extraction and quantification of macular drusen from fundal photographs[END_REF][START_REF] Phillips | Quantification of diabetic maculopathy by digital imaging of the fundus[END_REF][START_REF] Shin | Computer-assisted, interactive fundus image processing for macular drusen quantitation[END_REF], global thresholds [START_REF] Smith | A method of drusen measurement based on the geometry of fundus reflectance[END_REF] or fuzzy logic thresholds [START_REF] Thdibaoui | A fuzzy logic approach to drusen detection in retinal angiographic images[END_REF]. Most of these techniques proved that good noise removal algorithm is required for getting reliable results. We have also contributed to this subject by introducing new image processing algorithms for illumination correction [START_REF] Mora | Drusen Deposits Modelling with Illumination Correction[END_REF] and for drusens detection and modelling [START_REF] Mora | Modeling of Drusen Deposits Based on Retina Image Tridimensional Information[END_REF].

In this article it is proposed a methodology for automated drusens detection and quantification. The main goal is to create an algorithm which is as non-parametric as possible. Parameters should be calculated based on the image information, requiring less user intervention. The proposed methodology includes an image pre-processing step to compensate non-uniform illumination and contrast normalization; a detection algorithm that determines the number of drusens and an image modelling to quantify the area affected. The methodology was validated on a set of 22 images with the collaboration of three ophthalmologists.

Illumination Correction

In ophthalmologic imaging it is usually difficult to obtain good quality images. Illumination is not always uniform due to numerous factors: the retina surface is convex; the macula reflectance is lower than others parts of the retina; the patient eye is not stable, etc. The consequence is that the image has different contrast areas increasing the complexity of the quantification methods. We propose to correct the non-uniformity of illumination and contrast before further processing.

The proposed technique is designed for gray scale images and since the retina images used are mostly in RGB format, a RGB to Gray conversion is applied to those. This conversion is obtained using the green channel information, which is a procedure widely used in other retina image works [START_REF] Rapantzikos | Detection and segmentation of drusen deposits on human retina: potential in the diagnosis of age-related macular degeneration[END_REF][START_REF] Shin | Computer-assisted, interactive fundus image processing for macular drusen quantitation[END_REF][START_REF] Soliz | Computer-aided methods for quantitative assessment of longitudinal changes in retinal images presenting with maculopathy[END_REF]. It is chosen by its better dynamic range, better contrast and less influence by the illumination differences, while it maintains all retina structures visible.

Using the gray scale image, and to correct the effects of non-uniform illumination, we then determine the image background intensity and correct the image accordingly. Several methods were tested (FFT, Gaussian fitting, Gaussian blurring, Splines) [START_REF] Mora | Drusen Deposits Modelling with Illumination Correction[END_REF] and the smoothing spline fitting was the one to produce better contrasted images. This method fits the image background f(x,y) with a smoothing cubic spline (Fig. 2.b). This is a cubic spline S(x,y) that minimizes the function presented in equation [START_REF]World Health Organization, International Agency for the Prevention of Blindness: State of World's Sight -Vision 2020: the right to sight[END_REF], where p is a smoothing factor and ) , ( ( )

dy dx y x S D p y x S y x f p m y n x ∫∫ ∑∑ - + - ⋅ = = 2 2 1 1 2 ) , ( 1 ) , ( ) , ( (1)
The smoothing factor p, defined in the interval [0..1], controls the relative weight between fitting all data points (p=1) and having a smooth spline (p=0). For p=1 the spline is an interpolating spline that crosses every data point, while p=0 it maximizes the smoothness generating a straight line. The reference value for p was determined empirically within our set of test images and obtained p = 1e -6 , for an image resolution of 12.5 µm/pixel.

The corrected image is obtained by dividing the original image by the illumination pattern, generating a unitary average intensity. To recreate a grey scale image with 256 levels a constant value, 85 (1/3 scale), is multiplied by the image. However, depending on the image contrast this value can generate out of scale values or low contrast images. To overcome this problem a contrast normalization procedure was introduced. It involves the determination of the standard deviation (σ) of the image intensity and its normalization to a predefined value. The retinal vessels are one of the structures that are present in all images and provide good contrast with the background. We propose a method that evaluates locally the contrast between the vessels and the background and normalizes it to a predefined value.

Using a sliding window-based analysis the σ value and the mean intensity are calculated for each window. A window width of 1/8 of the optic disc diameter was considered to be adequate, as it can contain small vessel branches and background. Then from the 50 darker windows, that are typically windows containing vessels over a uniform background, is determined the median of the σ values ( median σ

). The image intensity is finally updated using the equation ( 2), where predefined σ is the predefined σ value (0.15) and A the constant background value (85).

  

     + - * ⋅ = 1 ) 1 ) , ( ( ) , ( _ y x image A y x image norm median predefined σ σ (2)

Drusens Detection

Our Drusens detection and quantification algorithm is based on the modelling of Drusen spots to ensure a shape consistent image segmentation. For getting a better performance from the modelling algorithm it is recommended that an image preprocessing is applied to calculate the approximate number of Drusens, their location and optionally their size.

We propose a maximum detection algorithm based on labelling the maximum gradient path. In medium resolution images there are more then just one pixel pointing to each intensity maximum. Therefore, it is possible to follow one or more ascending paths that inevitably reach an intensity maximum.

The algorithm begins by determining the image gradient using a 3x3 Sobel operator which evaluates intensity changes in the horizontal and vertical axis. The pixel gradient is a vector which is oriented to the ascending pixel. The two- dimensional gradient analyses the 3x3 neighbourhood of each pixel what eliminates several noise related artefacts. But when signal/noise ratio is low the algorithm can be significantly improved by changing every pixel gradient by the mean (vectors sum) of its neighbour's gradients. Fig. 3.a contains an example of an image after applying a 3x3 gradient mean. The original image represents two small Drusens.

In the first stage of the labelling procedure, initial label propagation, every pixel is examined in a top-left to bottom-right direction and is assigned a new label to all pixels that are not marked yet. Whenever a label is assigned to a pixel it is propagated to the neighbour pixel that is in the gradient azimuth direction. This labelling continues on the propagated pixels until the next pixel is already marked. Fig. 3.b presents the result of the labelling first stage. Every time the label propagation finishes on the same label that is being propagated, that pixel is marked as an intensity maximum. In the same situation when there is a different label the two labels are defined as being compatibles. This means that they belong to the same intensity maximum, i.e., the same Drusen spot.

The second stage of the labelling procedure is to apply the label compatibilities resulting on an image with as many labels as possible Drusens. The result is a segmented image where all pixels that contribute to a spot have the same label. It is important to notice that the algorithm is so far a non-parameterized algorithm.

But a problem arises when flat valleys or flat hills occur. In these cases not all gradient paths end on the same intensity maximum pixel since there might exist more then one pixel for that intensity maximum. Consequently, the algorithm will generate more intensity maximums then it really should. For solving this problem a merging strategy was followed. Using a connected components approach, the local intensity maximums are merged to other neighbour intensity maximums if they can be connected by a path that doesn't goes lower than a predefined amplitude.

At this stage drusens are detected and characterized by centre coordinates and size. This information will be used as the initial parameters for the modelling algorithm.

Drusens Modelling

The final phase is to analytically characterize the shape and intensity of the drusens, providing information to evaluate them quantitatively. The elevations shown on drusen areas on the tri-dimensional representation motivated us to create a model of the image intensity and from it extract the drusens dimensions. From the model it is computed the total affected area and drusens dimensions. This information can be valuable to an alongside comparison that evaluates sequence of images during a long term treatment.

In a three-dimensional view of the retina surface, as in Fig. 4.a, it can be seen a drusen that has a waveform similar to a Gaussian function superimposed to a flat background with low amplitude white noise. These observations and the possibility to describe the drusens by means of elementary functions motivated the use of Gaussian functions to model drusens. This approach allows the drusen modelling and its isolation from the background.

The modified Gaussian function used is given by:

0 2 2 2 * ) , ( Z d s Y s X y x e A y x G +         + - = .
(

) 3
where

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) x F y s s s y y x x Y y y x x X ⋅ = ⋅ - + ⋅ - - = ⋅ - + ⋅ - = θ θ θ θ cos sin sin cos 0 0 0 0 (4)
The Gaussian function used represented on equation ( 3) allows: translations in the xy plane (x 0 , y 0 ) and in the background plane (Z 0 ), amplitude scaling (A), rotation (θ) and shape modifications (s x , s F , d). These latter defines the width in the x-plane (s x ), the shape factor (s F ) used to obtain the width in y-plane, allowing stretched shapes in the xy plane, and the shape factor (d) that controls the amplitude profile between square shape, round shape and thin shape.

For searching the function parameters that best fits the image we adopted the Levenberg-Marquardt non-linear least squares optimization algorithm with predefined parameters constraints ( [START_REF] Marquardt | An algorithm for least-squares estimation of non-linear parameters[END_REF]). Because some Drusens are overlapped, a multiple function (one or more for each spot) simultaneous optimization must be used for studying the sum of all the Gaussian functions. In Fig. 4.b there is an example of superimposed Drusens modeled through multiple Gaussian functions.

Finally, for highlighting Drusen areas in the image it is shown the contour or the area of each detected spot at predefined percentage of its amplitude (without background).

Results

For evaluation of the drusen detection and quantification process a set of 22 images containing colour and greyscale fundus images was marked by three independent ophthalmologists. Our study was focused on the macula, as suggested by the Wisconsin Grading System [START_REF] Klein | The Wisconsin age-related maculopathy grading system[END_REF] a circular region of interest centred in the macula with 3000µm diameter. The specialists marked this set of images using the computer assisted procedure with the software MD3RI [START_REF] Mora | MD3RI a Tool for Computer-Aided Drusens Contour Drawing[END_REF] which were then used to compare the results produced by the proposed method (AD3RI) over the same images.

For the results analysis the percentage of drusen area, the false positive and false negative pixels, were computed. The results were then evaluated using the algorithm's sensitivity and specificity as well as using the Intraclass Correlation Coefficient (ICC) to evaluate the agreement between all the specialists and the AD3RI method and among each specialist.

The sensitivity and specificity indicators were obtained by comparing the processed image with each specialist marking. The average sensitivity and specificity were 0.664 and 0.963, respectively. The same indicators were evaluated among the specialists as a comparative indicator and the average sensitivity and specificity were comparable to the previous ones, 0.656 and 0.971, respectively.

As a comparative analysis method, the ICC with 95% confidence levels [START_REF] Mcgraw | Forming inferences about some intraclass correlation coefficients[END_REF] evaluated the areas' values between specialists and the AD3RI (ICC AD3RI ) and among

AD3RI

Specialist #1 Specialist #2 Specialist #3 specialists (ICC Spec ). These indicators were all within the 95% confidence interval and obtained an average agreement of 74.9% between the specialists and the AD3RI and an average agreement of 81.1% among the specialists.

Discussion

The development of methods to quantitatively measure drusens in a reproducible and accurate procedure will improve the quality of the follow up of this disease and potentiate epidemiologic studies and clinical trials. In this article we presented a new method to quantitatively measure drusens and its' comparison with 22 images marked by three independent ophthalmologists. The algorithm is based on the detection and modelling of drusens to automatically calculate the affected areas including an image pre-processing for non-uniform illumination correction. The presented non-uniform illumination correction algorithm was an important step to obtain a less parameterized methodology, since it is capable to deliver an image with normalized illumination and contrast to the following steps. The detection and modelling of drusens with Gaussian functions demonstrated to be capable to detect drusens even with low contrast and to segment the image maintaining the drusens typical shape. It also provides an analytical model which allows the determination of drusen indicators, such as confluence, average size, type of drusens and others.

From the statistical analysis of the algorithm with three specialists related to the 22 images was observed that the algorithm performs as do specialists. On the one hand, the sensitivity (0.664) was slightly higher than the one observed among the specialists (0.656). This is due to the algorithm detection of a higher overall area, increasing the probability of agreement with the specialists. On the other hand, the high specificity value (0.963) demonstrates that the automated process is very accurate and is also comparable to results obtained among the specialists (0.971). With the algorithm, the value is lower due to the increase of false positives, which results of the detection of more drusens. We consider that the lower quantification of drusens by the specialists is mainly due to the marking process being fastidious.

The other measure of agreement used, the Intraclass Correlation Coefficient (ICC), demonstrated that the agreement among specialists on the overall areas values is approximately 81% and that the automatically processed images have 74.9% of agreement when compared to specialists. This analysis showed a good concordance between the automated process and the specialists, verifying that the analysis is statistically significant, as all the obtained ICC values were within the 95% confidence levels.

From the above, the AD3RI method showed promising results to be tested by a wider panel of specialists and with new sets of images. It detected the most important drusen as did the specialists, having always the same level of detail independently of the specialists' attention and accuracy. Also, the qualitative feedback from the specialists confirmed that the results were accurate and that the software can play an important role in the follow-up of retina diseases.

Fig. 1 .

 1 Fig. 1. Examples of retina images obtained by fundus photography.

  derivative of the cubic spline:

Fig. 2 .

 2 Fig. 2. Illumination correction: (a) Original image f(x,y); (b) Estimated illumination model i(x,y); (c) Corrected image.

Fig. 3

 3 Fig. 3 Example of Drusens center detection algorithm. (a) Original Image and gradient; (b) Labelling 1 st stageinitial label propagation -(with label propagation on the two upper rows); (c) Labelling 2 nd stageapply compatibilities -(two Drusens centered on the highlighted pixels).

Fig. 4

 4 Fig. 4 (a) 3D original image view; (b) Drusens modeling using multiple Gaussian functions.

Fig. 5

 5 Fig. 5 Two examples of drusens markings.

Acknowledgments. This work was financed by the Fundação para a Ciência e Tecnologia, through the POCTI and POSI Research Programs (project nº POCI/SAU-ESP/57592/2004). The authors acknowledge the University of Aberdeen, Rudolfstiftung Hospital, Hospital Santa Maria and Faculdade de Ciências Médicas for supplying the retina images used in this work, for their valuable feedback on the software testing and for marking the sample images.


Hadi Aliakbarpour

Joao Filipe Ferreira

Kamrad Khoshhal

Jorge Dias

A Novel Framework for Data Registration and Data Fusion in Presence of Multi-modal Sensors

Keywords: Multi-modality, data registration, data fusion, Occupancy grid and Homography

This article presents a novel framework to register and fuse heterogeneous sensory data. Our approach is based on geometrically registration of sensory data onto a set of virtual parallel planes and then applying an occupancy grid for each layer. This framework is useful in surveillance applications in presence of multi-modal sensors and can be used specially in tracking and human behavior understanding areas. The multi-modal sensors set in this work comprises of some cameras, inertial measurement sensors (IMU), laser range finders (LRF) and a binaural sensing system. For registering data from each one of these sensors an individual approach is proposed. After registering multi-modal sensory data on various geometrically parallel planes, a two-dimensional occupancy grid (as a layer) is applied for each plane.

Introduction

Data registration and data fusion are two crucial issues in a multi sensory environment since heterogeneous sensors are substantially different. In this work we propose a novel framework for data registration and fusion in presence of various modalities such as image, range, sound and inertial data. This paper is structured as following. In Sec. 2 the main contribution is described. Sec. 3 is for related work. Sensors models are introduced in Sec.4. Sec. 5 is dedicated to registering data from different modalities and then applying a data fusion strategy. A preliminary result of an ongoing experiment is shown in Sec. [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Then in Sec. 7 conclusion and future work are discussed.

Contribution to Technological Innovation

The use of surveillance systems is growing in different areas such as airports, banks and human behavior interpretation. Each particular sensor has its own drawbacks. Although to overcome that nowadays researchers are trying to use several sensors from different modalities instead of using just a single type, however data registration and data fusion are still two critical and determinant phases in surveillance systems. Our main contribution is to work on these subjects and make a framework based on geometric data registration and probabilistic data fusion.

Related Works

Fusion of the sensors outputs is a crucial challenge related in such applications. A model namely JDL is proposed in [START_REF] Smith | Approaches to multisensor data fusion in target tracking: A survey[END_REF] for data fusion. Armesto et al. in [START_REF] Armesto | On multi-rate fusion for non-linear sampled-data systems: Application to a 6d tracking system[END_REF] presented an approach to model and fuse of non-linear data from multi-rate systems which have visual and inertial sensors. Bellotto and Hu also presented an approach in [START_REF] Bellotto | Vision and laser data fusion for tracking people with a mobile robot[END_REF] to fuse the LRF data and vision data by PTZ camera. Chakravarty and Jarvis in [START_REF] Chakravarty | Panoramic vision and laser range finder fusion for multiple person tracking[END_REF] discussed a fusion method for LRF and panoramic vision data to track people from a stationary robot simultaneously. Khan et al. in [START_REF] Khan | A homographic framework for the fusion of multi-view silhouettes[END_REF] presented a new approach which fuse different views silhouette data form some.

Sensors Models

A pinhole camera model [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] is used in this work. The 3D LRF in this work is built by moving a 2D LRF along one of its axes. A detailed model of this configuration can be found in [START_REF] Aliakbarpour | An efficient algorithm for extrinsic calibration between a 3d laser range finder and a stereo camera for surveillance[END_REF]. A detailed version of IMU model used in this paper can be seen in [START_REF] Lobo | Relative pose calibration between visual and inertial sensors[END_REF]. And eventually for the sound system, a Bayesian binaural system is considered for imitating human binaural system (see [START_REF] Ferreira | Implementation and calibration of a Bayesian binaural system for 3d localisation[END_REF]). The Bayesian binaural system (described in [START_REF] Pinho | A Bayesian binaural system for 3d soundsource localisation[END_REF]) composes three distinct and consecutive processors (Fig. 3-a): the monaural cochlear unit, which processes the pair of monaural signals {X 1 ,X 2 } coming from the binaural audio transducer system by simulating the human cochlea, so as to achieve a tonotopic representation (i.e. a frequency band decomposition) of the left and right audio streams; the binaural unit, which correlates these signals and consequently estimates the binaural cues and segments each sound-source; and, finally, the Bayesian 3D sound-source localization unit, which applies a Bayesian sensor model so as to perform localization of sound-sources in 3D space, using the egocentric frame of reference } {ε .

Multi-modal Data Registration and Data Fusion

In a multi-modal sensory environment, data registration and data fusion are two crucial issues. In this section we talk about these issues where there are four heterogeneous sensors. Prior to that, a coordinate reference } {W and also a 3D plane π ref inside this space will be introduced which are universal and common for all sensors inside the framework (see fig. 1-left). Let the origin of this coordinate frame be on our reference plane and also consider two coordinate vectors X and Y extending across (aligned) the plane (XY-plane of } {W corresponds to π ref ). Moreover assume that the z direction of this plane is parallel to the earth gravity vector but in an apposite direction. Having these definitions, π ref becomes a virtual horizontal plane [START_REF] Mirisola | Tracking from a moving camera with attitude estimates[END_REF] which is considered a common geometrical plane for all sensors in this setup. The idea is to geometrically register data observed by sensors with different modalities (cameras, LRFs and microphones) onto this plane (in section 5.4 it is extended to some parallel planes). After having the registration structure, a BOF [START_REF] Mekhnacha | Bayesian occupancy filter based "fast clustering-tracking" algorithm[END_REF] will be applied to achieve the final framework. Right: Considering a virtual camera for each real camera. The principle axis of such a virtual camera is coincided to the gravity vector. The rotation matrix between real camera } {C and virtual camera } C { ' is given by the IMU [START_REF] Lobo | Relative pose calibration between visual and inertial sensors[END_REF].

Image Data Registration

The idea of this work is to have all data coming from different sensors registered on a common reference plane, namely π ref .

= | { = c i n i I I . A virtual image plane, } 1... = | { ' c i n i I
, is considered for each camera. Such a virtual plane is assumed to be a horizontal plane at a distance f below the camera sensor, f being the focal length [START_REF] Mirisola | Tracking from a moving camera with attitude estimates[END_REF]. In other words, it can be said that beside of each real camera C in the setup, a virtual camera C' is also considered whose center, {C'}, is coincided to the center of the real camera {C}. So that the transformation matrix between {C'} and {C} will just have a rotation part and the translation part is a zero vector. In order to have In this case H is called infinite homography since there is just a pure rotation between real camera and virtual camera centers [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In this case the equation for H will become as [START_REF] Pinho | A Bayesian binaural system for 3d soundsource localisation[END_REF][START_REF] Criminisi | Accurate visual metrology from single and multiple uncalibrated images[END_REF] being R as the rotation matrix between } {C and } { ' C (real and virtual camera centers). Reminding that normal vector of the virtual camera's plane is coincided to the earth gravity vector, then the rotation matrix R can be achieved by using the IMU coupled to the camera. The Camera Inertial Calibration Toolbox is used in order to calibrate a rigid couple of a IMU and camera [START_REF] Lobo | Relative pose calibration between visual and inertial sensors[END_REF].

1 ' ' = - K R K H I I
Homography Between Virtual Image Plane and π π π π ref : The homography matrix between the virtual image plane I' (the image of virtual camera) and π ref needs to be computed. A homography matrix H can be represented by its axis and vertex and cross-ratio parameters:

a v a v I H T T 1) ( = - + µ (3)
in which v is the vertex coordinates and a is the axis line [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Mirisola | Tracking from a moving camera with attitude estimates[END_REF][START_REF] Criminisi | Accurate visual metrology from single and multiple uncalibrated images[END_REF]. For the case of a 3D plane being parallel to the image plane, like in our case, the Eq. 3 becomes a simple matrix by considering The intersection of the arrival sound plane and π ref will be registered as a line.

Range Data Registration

LRF is one of the interesting sensors in our multi-modal setup. The model of our LRF is described in [START_REF] Aliakbarpour | An efficient algorithm for extrinsic calibration between a 3d laser range finder and a stereo camera for surveillance[END_REF]. Geometrical registration of range data is the goal of this section. The general intention is to be able geometrically registering data coming from LRF on the reference plane π ref (see Fig. , can be expressed by

x T x L L W ref           1 0 0 0 0 0 1 0 0 0 0 1 = π (5)
where L W T is the transformation matrix between LRF's local reference and can be acomputed by the method described in [START_REF] Aliakbarpour | An efficient algorithm for extrinsic calibration between a 3d laser range finder and a stereo camera for surveillance[END_REF].

Registration of Azimuth Plane of Arrival of Sound

In this work, we use a simplified version of the sensor model described in Sec. 4, taken the decomposition equation ),

θ θ θ S P S P P S P ( 6
)
where τ is the interaural time-difference binaural cue [START_REF] Pinho | A Bayesian binaural system for 3d soundsource localisation[END_REF],

} , { 0 n θ θ θ L ∈
denotes an azimuth angle taken from a discretised span of n angles, θ S is a binary variable signaling the existence of a sound-source within the azimuthal plane at θ extending from the frame of reference {ξ} (the so-called azimuth plane of arrival for the respective sound-waves),

) (θ P is an uninformative uniform distribution, ) | ( θ θ S P
is the prior on the existence of sound-sources for a given θ, also chosen to be a uniform distribution, and

) | ( θ τ θ S P
is the binaural sensor model, a set of normal distributions obtained through calibration (see [START_REF] Ferreira | Implementation and calibration of a Bayesian binaural system for 3d localisation[END_REF]), which indicate the probability of the measurement of θ (i.e. the angle for the azimuth plane of arrival) knowing whether or not a sound-source is present in that plane. Using Bayes rule it is possible to invert the sensor model so as to obtain

) | S ( P θ τ θ
; an auditory saliency map can then be constructed, and a MAP (maximum a posteriori) method can be used to extract the azimuth angle estimate for the most salient sound-source (see Fig. 3-b). Based on these definitions it is possible to also register sound information on the reference plane π ref (see Fig 2). The idea for sound registering is to calculate the intersection between the π ref and the plane of arriving sound. The result will be registered as a line on the π ref . Such an extension needs to be done for cameras, LRFs and microphones. For LRF it becomes simple since a general method to find the homogeneous transformation between LRF and } {W has been defined in Sec. 5.2. For arrival sound it is exactly the same method described in Sec. 5.3. Here we continue to describe this extension for images. The intention is to calculate homography for any of these layers having the homography of the reference layer

I ref H π .
For doing so we use Khan's approach described in [START_REF] Khan | A homographic framework for the fusion of multi-view silhouettes[END_REF].

Multi-layer Homography-based Occupancy Grid

Fig. 3-c shows a schematic of the proposed multi-layer framework for data registration. The idea is to make an occupancy grid layer for each plane of the registration framework. In order to make the final data fusion a BOF grid [START_REF] Mekhnacha | Bayesian occupancy filter based "fast clustering-tracking" algorithm[END_REF] per each plane is proposed as an occupancy layer (see Fig. 3-(c)). It can be used as input of a classification algorithm in order to classify the objects inside the scene and then be used by a tracking algorithm. We use the concept and also the model described in [START_REF] Mekhnacha | Bayesian occupancy filter based "fast clustering-tracking" algorithm[END_REF] for using BOF in our case. Here the variables definition is the same as [START_REF] Mekhnacha | Bayesian occupancy filter based "fast clustering-tracking" algorithm[END_REF]. Using that the probability joint distribution for the model becomes also the same as what is defined in [START_REF] Mekhnacha | Bayesian occupancy filter based "fast clustering-tracking" algorithm[END_REF] (l being the layer index):

) | ( ) 1 | ( ) 1 | ( ( ) 1 ( = ) ... 1 ( 1 = 1 O t A Z P t A t O P t A t A P P t A P Z Z t O t A t A P l c l t i S i l c l l c l c l l c l l c l t S t l c l c l l c l ∏ - - - - (7)

Experiments

As a preliminary experiment result the image data registration part is performed for just one layer. We used the CVLab -EPFL dataset [START_REF] Fleuret | Multi-camera people tracking with a probabilistic occupancy map[END_REF]. Fig. 4:(a)-(d) show four images from four views and their projection on the ground plane using homography concept. Then after background subtraction, they are combined and intersected in order to obtain the feet positions on the ground plane (Fig. 4-e).

Conclusion and Future work

Data registration and data fusion have been investigated in this paper where there are four heterogeneous sensors such as camera, microphone, laser range finder and inertial sensor. Then a multi-modal multi-layer framework is porposed. In such a framework firstly sensory data coming from different sensors are geometrically registered on different parallel planes, then an occupancy grid is applied for each layer. A prelimanry experiment result has been shown for image data registration in just one layer in Sec. 6. Data registration for other layers, range and sound data and also constructing an occupancy grid for each layer will be remained as our future work.

Fig. 1 .

 1 Fig. 1. Left: World coordinate system, W, and the reference plane, ref π . Middle: Registering

Fig. 1 -

 1 middle shows a setup with c reference plane π ref . In this setup each camera is rigidly fixed with an IMU (Inertial Measurement Unit). The intention is to register a 3D point ) camera C onto the reference plane π ref as x ref π using the IMU. The image plane of these cameras can be expressed as } 1...

  are considered: Firstly, a 3D point X W is projected on the camera reference frame ({C}) using projection matrix). Secondly, x I (the imaged point on the camera's image plane) is projected to the corresponding virtual image plane as x I ' . This can be done by having the related homography transformation, namely point on the virtual image plane, is reprojected to the main reference plane π ref by having the related homography transformation matrix, called ' I ref H π . For the first step it is by the camera model described in the Sec. 4. The second and third steps are described in the following two sub-section. Assuming to already have Inertial Compensated Homography: Lets consider I'' X as re-projection of a point I X from image plane I of camera C onto the virtual image plane I' of camera C' (see figure 1-right). Then the following general equation can be used: between real image plane and its virtual plane.

  homography matrix between the reference plane π ref and virtual image plane I'.

Fig. 2 .

 2 Fig. 2. Left: Range data registration. Right: Registration of azimuth plane of arrival of sound: The intersection of the arrival sound plane and π ref will be registered as a line.

  2-left). More specifically, for a 3D point X in the scene which is observed in the LRF local reference frame as homogeneous form), its projection on π ref , in homogeneous form ,

Fig. 3 . 5 . 4

 354 Fig. 3. (a): Bayesian binaural system. (b): Auditory saliency map and azimuth plane of arrival extracted for the most salient sound source. (c): Multi-modal multi-layer data registration and data fusion framework.

Fig. 4 .

 4 Fig. 4. (a) to (d): Images from different views and their projection on the ground plane using homography concept. e: combining four ground-projected images in order to have the feet intersection.

Acknowledgment: Hadi Ali Akbarpour and Joao Filipe Ferreira are supported by the FCT (Portuguese Fundation for Sceince and Technology). This work is partially supported by the European Union within the FP7 Project PROMETHEUS, www.prometheus-FP7.eu.


Fernanda Da
email: fernanda@isr.uc.pt

C Faria

Jorge Batista
email: batista@isr.uc.pt

Helder Araújo
email: helder@isr.uc.pt

Bio-inspired Binocular Disparity with Position-Shift Receptive Field

Keywords: Depth perception, Binocular disparity, Bio-inspired stereo vision, Energy model

Depth perception starts with the binocular interaction receptive fields of simple cells modeled by two Gabor functions followed by a half-squaring function. Simple cells do not have reliable disparity computation. Complex cells combine two simple cells in quadrature. They are better adapted to encode disparity information. The image disparity can be determined by fixing the receptive field in one eye and varying it in the other. Pooling information of spatial frequency and orientation is very important to improve the quality of results of real world stereograms. In this work, a bio-inspired method to calculate binocular disparity based on the energy model for depth perception of real world images is described and implemented. The performance of the proposed method is also evaluated in what concerns the algorithm computational cost, root-mean-square error and percentage of wrongly matched pixels.

Introduction

In the literature several methods to compute the binocular disparity are described such as statistical approaches, probabilistic approaches and artificial neural networks [START_REF] Scharstein | A Taxonomy and Evaluation of Dense Two-frame Stereo Correspondence Algorithms[END_REF]. The focus of this work is on the application of bio-inspired models.

The brain combines information from both eyes and constructs a visual world that is perceived as single and three dimensional. Primary visual cortex (V1) is the first site at which single neurons can be activated by stimuli in both eyes. Primary visual cortex neurons encode information specifically about the relationship between the images in the two eyes [START_REF] Cumming | The Physiology of Stereopsis[END_REF]- [START_REF] Qian | Computing Stereo Disparity and Motion with Known Binocular Cell Properties[END_REF].

Biologically plausible models for stereo vision inspired in the visual cortex employ receptive field (RF) profiles of binocular cells for disparity computation. A region of a scene that is visible by left and right eyes is named binocular view field. All points in the binocular view field have different positions in the left and right eyes [START_REF] Anzai | Neural Mechanisms for Processing Binocular Information. I. Simple Cells[END_REF], [START_REF] Qian | Computing Stereo Disparity and Motion with Known Binocular Cell Properties[END_REF], [START_REF] Ohzawa | Stereoscopic Depth Discrimination in the Visual Cortex: Neurons Ideally Suited as Disparity Detectors[END_REF]. Binocular disparities are the small positional differences between the images from the left and right eyes. The slightly different viewpoint is due the horizontal separation of the eyes in the head. Depth perception is deduced through binocular disparity [START_REF] Anzai | Neural Mechanisms for Processing Binocular Information. I. Simple Cells[END_REF]- [START_REF] Qian | Computing Stereo Disparity and Motion with Known Binocular Cell Properties[END_REF].

Stereoscopic vision at the neuronal level based on events within the primary visual cortex can be studied through two subdivisions of cortical cells, named simple cells and complex cells. Simple cells are composed of orientation selective RFs that have spatially separate regions. These regions are mutually antagonistic, i.e., they respond to either ON (excitatory) or OFF (inhibitory) to a stimulus [START_REF] Parker | Binocular Depth Perception and the Cerebral Cortex[END_REF], [START_REF] Ohzawa | Stereoscopic Depth Discrimination in the Visual Cortex: Neurons Ideally Suited as Disparity Detectors[END_REF]- [START_REF] Read | Early Computational Processing in Binocular Vision and Depth Perception[END_REF]. The receptive fields of binocular simple cells are selective to binocular disparity and they are also affected by several stimulus parameters, such as [START_REF] Anzai | Neural Mechanisms for Processing Binocular Information. I. Simple Cells[END_REF], [START_REF] Qian | Computing Stereo Disparity and Motion with Known Binocular Cell Properties[END_REF], [START_REF] Ohzawa | Stereoscopic Depth Discrimination in the Visual Cortex: Neurons Ideally Suited as Disparity Detectors[END_REF]: position, contrast polarity, spatial frequency and orientation. The other type of V1 cells, complex cells, are composed of orientation selective RFs and respond to a stimulus anywhere within the RF, i.e., there is not spatially separate regions. Complex cells are not sensitive to changes in stimulus position and contrast polarity. The fact that they do not suffer from the influence of irrelevant parameters makes them well adapted to encode disparity information [START_REF] Parker | Binocular Depth Perception and the Cerebral Cortex[END_REF]- [START_REF] Hubel | Receptive Fields, Binocular Interaction, and Functional Architecture in the Cat's Visual Cortex[END_REF].

The binocular energy model is a model for disparity-selective complex cells based on a combination of simple-cell subunits. Binocular interaction RFs of a simple cell is described as the product of left and right eye RFs. This interaction can be demonstrated by a linear binocular filter followed by a half-power function [START_REF] Anzai | Neural Mechanisms for Processing Binocular Information. I. Simple Cells[END_REF], [START_REF] Qian | Computing Stereo Disparity and Motion with Known Binocular Cell Properties[END_REF], [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF], [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF].

This work describes a bio-inspired method to calculate binocular disparity based on energy model for depth perception of real world images [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF]- [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF].

Contribution to Technological Innovation

Visual perception is essential for many applications, especially in Robotics. For example, in applications related to human-robot interaction, the machine needs to cope with and share the human environment in a friendly way. For applications such as assistants in surgery, service robots, caretakers for the elderly and disabled people, a technological solution is necessary requiring advanced perception capabilities as well as advanced skills that allow a robot to perform the desired tasks.

The properties of a scene analyzed by a robot head can involve tracking movements, focusing by vergence and 3D perception of the scene structure. The relative orientation of the eyes, the target location and depth are important for surveillance activities [START_REF] Batista | A Focusing-by-Vergence System Controlled by Retinal Motion Disparity[END_REF] [START_REF] Batista | Binocular Tracking and Accommodation controlled by Retinal Motion Flow[END_REF].

The algorithm presented in this paper is a preliminary work related to the development of bio-inspired visual perception modules for robotic heads. The model created is based on two algorithms [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF], [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF]. The principal idea is to develop a faster algorithm to compute disparity maps employing the concepts of energy model for real world stereograms.

Stereoscopic Model

Binocular disparity forms the basis of stereoscopic depth perception. The binocular energy model is based on a combination of simple cells. Each complex cell is the linear combination of the outputs of four simple type subunits. A typical binocular simple cell is identical to two Gabor functions (equation 1), one for each of its receptive fields in the left and the right retinas [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF].

) ' cos( 2 ' ' exp ) , , ( 2 2 2 φ ω σ φ -         - - = x y x y x g (1)
where ω is the cell preferred spatial frequency (2πf), θ is its preferred orientation, φ is its phase, x'=xcos(θ)+ysin(θ) and y'=ycos(θ)-xsin(θ). The standard deviation σ is computed as follows:

1 2 1 2 2 2 ln 5 . 1 5 . 1 - + × = f π σ (2)
where the half-maximum bandwidth is 1.5 octaves [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF].

The left and right RFs are shifted in position (equations 3 and 4 -the receptive field is fixed in one eye and varying it in the other) [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF].

) , , ( ) , (

φ y x g y x g l = (3)
) , , ( ) , (

φ y d x g y x g r + = ( 4
)
where d is a positional difference between the left and right eyes Gabor functions.

The response of a simple cell can be determined by first filtering, for each eye, the retinal image by the corresponding receptive field profile (equation 5), and then adding the two contributions from the two eyes (equation 6) [START_REF] Qian | Physiological Computation of Binocular Disparity[END_REF].

∫ +∞ ∞ - = dy dx y x g y x I L l l ) , ( ) , ( (5)
R L S + = (6)
The response of a complex cell can be modeled by summing the squared outputs of a quadrature pair of simple cells, i.e., their spatial phases differ by 90 degrees while all the other parameters of the cells are identical (equation 7). The combination of simple cells in quadrature eliminates the Fourier phase dependence. Therefore, complex cells encode reliable disparity information [START_REF] Qian | Binocular Disparity and the Perception of Depth[END_REF], [START_REF] Read | Early Computational Processing in Binocular Vision and Depth Perception[END_REF], [START_REF] Qian | Physiological Computation of Binocular Disparity[END_REF].

S S C + = (7)
Briefly, the process of depth perception begins with the binocular interaction RFs of simple cells. They can be modeled by two Gabor functions followed by a halfsquaring function. Simple cells depend on the Fourier phase of the stimulus. As a result they do not perform a reliable disparity computation. Complex cells combine two simple cells in quadrature and eliminate the Fourier phase dependence. They are better adapted to encode disparity information than simple cells.

Computing Disparity

To compute the image disparity, the most responsive complex cell at each spatial location is selected at zero phase disparity for the first spatial frequency (f). We use 6 spatial frequencies (0.5, 1.0, 2.0, 4.0, 8.0 and 16.0 cycles per degree) and 6 RF orientations (30, 60, 90, 120, 150 and 180 degrees). Next the disparity responses of the complex cells population with the fixed f are computed, for the full range of the positional differences (d) for each θ. Complex cells at a local maximum or a local minimum are selected and if it is not the first spatial frequency f, then the selected complex cell is the cell with the closest filtered d value computed. If spatial frequency f = 0.5, then the selected complex cell is the one with the biggest d value. Another situation is when a local maximum or a local minimum does not exist then the d value is -99. The value -99 means that the d value will be ignored in the pooling final disparity map. The result images and the filtered images are saved. A 5x5 median filter is used. The size of the receptive field decreases when f increases. The range of values for d is split into several steps. The steps have a small and fixed size. Fig. 1 shows a block diagram of the algorithm. The outputs of the algorithm are 36 images (6 spatial frequencies and 6 orientations) that will be pooled for the computation of the final disparity map. The images employed in the pooling are the images saved without filtering. The values of the filtered images are used only during the process of selection of d (Fig. 1).

Pooling

A disparity map created by pooling the cells with multiple spatial frequencies and orientations produces results that are more reliable and accurate than those that would be obtained without the pooling step. The pooling is the same method described in [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF] where each spatial position of the 36 disparity maps with different spatial frequencies and orientations are combined to compute the average. The resulting average is compared with each spatial position of all maps and the position whose value is farthest from the mean is removed. The process is repeated until the total of spatial positions is reduced by half.

Evaluation Methodology and Results

The algorithm evaluation was done employing two quantitative measures of fit quality: root-mean-square error (R) and percentage of wrongly matched pixels (B) [13]. R is computed using equation 8,

∑ - = ) , ( 2 | ) , ( ) , ( | 1 y x T C y x d y x d N R ( 8
)
where d C is the computed disparity map, d T is the ground truth map and N is the total number of pixels.

Equation 9 is used to calculated B with a disparity error δ d = 1.

∑ > - = ) , ( ), | ) , ( ) , ( (| 1 y x d T C y x d y x d N B δ (9)
Fig. 2 and Fig. 3 show the resulting disparity maps for the images from the Middlebury stereo repository [START_REF] Scharstein | A Taxonomy and Evaluation of Dense Two-frame Stereo Correspondence Algorithms[END_REF] (http://cat.middlebury.edu/stereo). Fig. 4 also presents an output disparity map for the method described in section 3. Table 1 shows the results of root-mean-square error (R) and the percentage of wrongly matched pixels (B) to Figures 2 (d) and 3 (d). Errors of table 1 are only evaluated in non-occluded regions [START_REF] Scharstein | A Taxonomy and Evaluation of Dense Two-frame Stereo Correspondence Algorithms[END_REF] (http://cat.middlebury.edu/stereo). Table 1. Evaluation results. R B Fig. 2 (d) 2.0 20 Fig. 3 (d) 1.5

Discussion of the Results

Several tests were performed with the algorithm described, using a Pentium 4 with 3.00GHz and 2.00GB of RAM memory. The operating system used was Suse Linux 10.

The algorithm proposed by Chen and Qian [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF] (using a coarse-to-fine approach) was implemented in C++ and tested with the goal of analyzing its performance with real world images and to check whether it would be possible to run it in real-time. All functions that are used to perform convolution were implemented in the frequency domain to allow a decrease in the processing time. However, in spite of several optimizations, it was not possible to obtain real-time performance.

The coarse-to-fine algorithm spent approximately twenty seven minutes to process an image with 256x256 pixels (Fig. 4. (a) and (b)). The coarse-to-fine algorithm can be improved using multiple processors and a library with multithreading support. The algorithm proposed by Read and Cumming [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF], was also implemented in C++. However its processing time was higher than the processing time of the coarse-to-fine algorithm. The test with the images presented in Fig. 4. (a) and (b) required five days and nineteen hours of running time.

The proposed algorithm in this work required 7 minutes and 33 seconds to obtain the disparity map of Fig. 4 (c). For Fig. 2 (d) the algorithm spent 8 min and 45 sec. and for Fig. 3. (d) the disparity map was computed in 6 min. and 20 sec. The method described method is based partially on the algorithms presented in [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF] and [START_REF] Chen | A Coarse-to-Fine Disparity Energy Model with Both Phase-Shift and Position-Shift Receptive Field Mechanisms[END_REF].

The root-mean-square error and the percentage of bad matching pixels presented in table 1 are similar to the results shown in [START_REF] Read | Sensors for Impossible Stimuli May Solve the Stereo Correspondence Problem[END_REF] but do not have yet the precision of stereo algorithms in 13 (http://cat.middlebury.edu/stereo).

Conclusion and Future Work

The estimation of disparity maps using the energy model is slow. Real images are inherently complex due to the non-uniformity of the disparity values. This problem can be handled by pooling images with different orientations and spatial frequencies. On the other hand, pooling requires additional processing time. One alternative to develop a real time system (using models of populations of cortical neurons) is to employ DSPs and FPGA chips [START_REF] Tsang | Neuromorphic Implementation of Active Gaze and Vergence Control[END_REF].

As future work, we intend to improve the precision of the algorithm and its suitability for tracking applications. Depth accuracy is crucial for binocular vergence control.

Fig. 1 .

 1 Fig. 1. Block diagram of the algorithm.

Fig. 2 .

 2 Fig. 2. (a) Left image; (b) right image; (c) the ground truth disparity map; (d) the estimated disparity map.

Fig. 3 .

 3 Fig. 3. (a) Left image; (b) right image; (c) the ground truth disparity map; (d) the estimated disparity map.

Fig. 4 .

 4 Fig. 4. (a) Left image; (b) right image; (c) the estimated disparity map.