Michael Felsberg

Fredrik Larsson

Learning Bayesian Tracking for Motion Estimation

A common computer vision problem is to track a physical object through an image sequence. In general, the observations that are made in a single image determine the actual state only partially and information from several views has to be merged. A principled and wellestablished way of fusing information is the Bayesian framework. In this paper, we propose a novel way of doing Bayesian tracking called channelbased tracking. The method is related to grid-based tracking methods, but differs in two aspects: The applied sampling functions, i.e., the bins, are smooth and overlapping and the system and measurement models are learned from a training set. The results from the channel-based tracker are compared to state-of-the-art tracking methods based on particle filters, using a standard dataset from the literature. A simple computer vision experiment is shown to illustrate possible applications.

Introduction

Due to the projection into the image plane, the 3D motion of an object cannot be determined uniquely from monocular views. However, if the object motion is non-degenerate, i.e., if the motion is not restricted to certain subspaces, 3D tracking of the object allows one to estimate its full state, i.e., all considered dimensions of the motion state.

Problem Formulation

Assume a 3D object moving in space and we want to estimate a subspace of its configuration, in practice mostly its 3D position. The observation data consists of consecutive 2D views of the object, taken from a single camera. For simplicity the camera will not move here, but the case where the camera moves is believed to be solvable within the same framework. The relative configuration of image plane and motion trajectory are assumed to be such that all state space dimensions of interest can be estimated, i.e., the system is observable [START_REF] Dahl | On observer error linearization for perspective dynamic systems[END_REF]. It is further assumed that the observation data is already extracted from the sequence, i.e., a sequence of observations (image coordinates) is the input. The output of the method consists of a sequence of 3D coordinates. All models are assumed to be non-linear and the noise terms are assumed as non-Gaussian. Under these circumstances, we believe that only learning systems can provide stable and robust solutions. Therefore, the proposed approach makes no use of manually engineered system models or measurement models. Both models are learned from training data, where during the learning phase samples for system state sequences are required. Note that both models are multi-modal and thus non-trivial to learn.

Learning the measurement model implies that no camera calibration is required in contrast to state-of-the-art methods that rely on accurate calibration. Coordinate systems can be chosen fairly arbitrarily and according to the requirements of the application in mind. Learning the system model, i.e., learning the dynamics of objects, allows one to change to suitable, application dependent state space coordinate systems. For instance, when it comes to car safety systems one might consider the distance in terms of field of safe travel [START_REF] Gibson | A theoretical field-analysis of automobile-driving[END_REF] which is more relevant than the Euclidean distance to the car. In contrast to state-of-the-art methods using predefined system models, the same implementation can be used on different state spaces.

Related Work

The most relevant work from the literature can be found in the area of non-linear, non-Gaussian Bayesian tracking [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Isard | CONDENSATION -conditional density propagation for visual tracking[END_REF]. In Bayesian tracking, the current state of the system is represented as a probability density function of the system's state space. In the prediction step, this density is modified according to the system model and a new, typically smoother density is obtained as the prediction of the next system state. In the observation step, measurements of the system are used to update the predicted density. Typically, the density is sharpened by inference and it serves as the next system state.

Non-linear, non-Gaussian Bayesian tracking excludes (extended) Kalman filtering. Common approaches are particle filters and grid-based methods [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. Whereas particle filters apply Monte Carlo methods for approximating the relevant density function, grid based methods discretize the state-space, i.e., apply histogram methods for the approximation. In the case of particle filters, densities are propagated through the models by computing the output for individual particles. Grid-based methods use discretized transition maps to propagate the histograms and are closely related to Bayesian occupancy filtering [START_REF] Coué | Using Bayesian programming for multi-sensor multitarget tracking in automotive applications[END_REF].

An extension to grid based methods is to replace the rectangular histogram bins with overlapping, smooth kernel functions. This leads to the recently proposed channel-based tracking [START_REF] Felsberg | Fusing dynamic percepts and symbols in cognitive systems[END_REF]. Channel-based tracking implements Bayesian tracking using channel representations [START_REF] Granlund | An Associative Perception-Action Structure Using a Localized Space Variant Information Representation[END_REF] and linear mappings on channel representations, so-called associative networks [START_REF] Johansson | The application of an oblique-projected landweber method to a model of supervised learning[END_REF]. The term channel is well established in vision literature for a representation using a band pass tuning function [START_REF] Howard | Binocular Vision and Stereopsis[END_REF], while some may at times have a reason to view it as a population coding [START_REF] Zemel | Probabilistic interpretation of population codes[END_REF][START_REF] Pouget | Information processing with population codes[END_REF]. The main advantage compared to grid-based methods is the reduction of quantization effects and reduced computational effort.

Channel representations can be considered as regularly sampled kernel density estimators [START_REF] Felsberg | Channel smoothing: Efficient robust smoothing of low-level signal features[END_REF], implying that channel-based tracking is related to kernelbased prediction for Markov sequences [START_REF] Georgiev | Nonparamtetric system identification by kernel methods[END_REF][START_REF] Yakowitz | Nonparametric density estimation, prediction, and regression for markov sequences[END_REF]. In the cited work, system models are estimated in a similar way as described here, but the difference is that we consider sampled densities, making the algorithm much faster than the cited work. Another way to represent densities in tracking are Gaussian mixtures (e.g. [START_REF] Han | Probabilistic fusion tracking using mixture kernelbased Bayesian filtering[END_REF]) and models based on mixtures can be learned using the EM algorithm, cf. [START_REF] North | Learning dynamical models using expectation-maximisation[END_REF], although the latter method is restricted to uni-modal cases (Kalman filter).

The main novelty of this paper compared to [START_REF] Felsberg | Fusing dynamic percepts and symbols in cognitive systems[END_REF] is the generalization of channel-based tracking to multi-dimensional state vectors, higher-order system models, and generic measurement models. Higher-order system models are required for tracking partially observable states and generic measurement models are obtained through learning. Furthermore, this paper contains a quantitative evaluation of results in comparison to particle filters and grid-based methods.

Organization of the Paper

The paper is organized as follows. After the introduction, the methods required for further reading are introduced: Bayesian tracking and channel representations of densities. As novelties of this paper, the channel-based tracking algorithm is formulated for the multi-dimensional case and the learning of system and observation models is described. In section 4 channel-based tracking is compared to particle filters and grid-based methods in a standard performance test for Bayesian tracking. A simple vision-based experiment is shown to illustrate applicability in practical problems. In section 5 we discuss the achieved results.

Bayesian Tracking and Channel Representations

Channel-based tracking can be considered as a generalization of grid-based methods for implementing non-linear, non-Gaussian Bayesian tracking. Hence we give a brief overview on Bayesian tracking and channel representations.

Bayesian Tracking

For the introduction of concepts from Bayesian tracking we adopt the notation from [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. Bayesian tracking is commonly defined in terms of a process model f and a measurement model h, distorted by i.i.d. noise v and n

x k = f k (x k-1 , v k-1 ) (1) z k = h k (x k , n k ) . (2)
The symbol x k denotes the system state at time k and z k denotes the observation made at time k. Both models are in general non-linear and time-dependent.

The current state can be estimated, given that the previous state and all previous observations are known, by using the prediction equation. Assuming a Markov process of order one allows us to consider the conditional density of the novel state as an integral over its conditional density given the previous state

p(x k |z 1:k-1 ) = p(x k |x k-1 )p(x k-1 |z 1:k-1 ) dx k-1 . (3)
Taking into account the new measurement, the prediction is updated through

p(x k |z 1:k ) = p(z k |x k )p(x k |z 1:k-1 ) p(z k |x k )p(x k |z 1:k-1 ) dx k . (4)
In the case of non-linear problems with multi-modal densities, basically two relevant approaches for implementing (3) and ( 4) are known: The particle filter and grid-based methods. Particle filters are considered superior to grid-based methods concerning accuracy [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF].

Particle filter methods apply Monte Carlo simulation for approximating the relevant densities, and the crucial step is the re-sampling of particles from the previous estimates in order to avoid degeneracy and loss of diversity. Cumulative histograms are required for the re-sampling step [START_REF] Isard | CONDENSATION -conditional density propagation for visual tracking[END_REF], and in case of many particles and high-dimensional spaces, the method suffers from the computational burden.

Grid-based methods assume a discrete state space such that the densities are approximated with histograms. Thus, conditional probabilities of state transitions are replaced with linear mappings. In contrast to [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] where densities were formulated using Dirac distributions weighted with discrete probabilities, we assume band-limited densities and apply sampling theory. However, the final equations for grid-based tracking and the model assumptions remain the same.

For the remainder of this section, sampling is meant in the signal processing sense and not in its meaning as a stochastic sampling. Sampling the posterior of the previous time step gives us

w i k-1|k-1 p(x k-1 |z 1:k-1 ) * δ(x i -x k-1 ) (5)
where * denotes convolution and δ(x ix) is the Dirac impulse at x i . Note that the sampling itself is time independent. Accordingly, the sampled prior pdf at time k and the sampled posterior at time k are obtained as

w i k|k-1 p(x k |z 1:k-1 ) * δ(x i -x k ) (6)
w i k|k p(x k |z 1:k ) * δ(x i -x k ) . (7)
Plugging ( 5) and ( 3) into (6) and applying the power theorem gives us

w i k|k-1 = p(x k |x k-1 )p(x k-1 |z 1:k-1 ) dx k-1 * δ(x i -x k ) = j f ij k w j k-1|k-1 (8)
where

f ij k = p(x k |x k-1 ) * δ(x i -x k ) * δ(x j -x k-1
). ( 9) Accordingly, plugging ( 6) and ( 4) into [START_REF] Granlund | An Associative Perception-Action Structure Using a Localized Space Variant Information Representation[END_REF] gives us

w i k|k = p(z k |x k )p(x k |z 1:k-1 ) * δ(x i -x k ) p(z k |x k )p(x k |z 1:k-1 ) dx k = h i k (z k )w i k|k-1 j h j k (z k )w j k|k-1 (10)
where

h i k (z k ) = p(z k |x k ) * δ(x i -x k ). ( 11
)
Note that ( 8) and ( 10) are exactly the same equations as in the formulation of grid-based methods in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], except for change of notation. Re-writing grid-based methods in terms of the sampling theorem allows us to analyze approximation errors more systematically. Grid-based methods require the more samples the higher the upper band limit of the pdf, i.e., the wider the characteristic function ϕ x (t) = E{exp(it T x)}. Furthermore, the sampling-based formulation makes it easier to switch to other sampling schemes than the ordinary impulse sampling.

Channel Representations of Densities

The channel representation [START_REF] Granlund | An Associative Perception-Action Structure Using a Localized Space Variant Information Representation[END_REF][START_REF] Snippe | Discrimination thresholds for channel-coded systems[END_REF] can be considered as a way of sampling continuous functions or, alternatively, as histograms where the bins are replaced with smooth, overlapping basis functions, see e.g. [START_REF] Pampalk | Using Smoothed Data Histograms for Cluster Visualization in Self-Organizing Maps[END_REF]. Consider a density function p(x) as a continuous signal that is sampled with a smooth basis function, e.g., a B-spline. It is important to realize here that the sampling takes place in the dimensions of the stochastic variables, not along the time axis k. Sampling is meant in the signal processing sense, not in its meaning as a stochastic sampling. It has been shown in the literature that an averaging of a stochastic variable in channel representation is equivalent to the sampled kernel density estimator with the channel function as kernel function [START_REF] Felsberg | Channel smoothing: Efficient robust smoothing of low-level signal features[END_REF]. For the purpose of Bayesian tracking, we replace the sampled densities ( 5) -( 7) with

w i k1|k2 p(x k1 |z 1:k2 ) * b(x i -x k1 ) (12)
where b(x) is a channel basis function. For the remainder of this paper it is chosen as [START_REF] Forssén | Low and Medium Level Vision using Channel Representations[END_REF] b(x) 2a π

n cos 2 (ax n ) |x n | < π 2a 0 otherwise. ( 13
)
Here a determines the relative width, i.e., the sampling density. According to [START_REF] Felsberg | Channel smoothing: Efficient robust smoothing of low-level signal features[END_REF], the channel representation reduces the quantization effect compared to ordinary histograms by a factor of up to 20. Switching from histograms to channels thus allows us to reduce computational load by using fewer bins, to increase the accuracy for the same number of bins, or both at the same time. In theory, the reduction of bins is limited by the upper band limit of the density function, but in practice the number of drawn samples (in a stochastic sense now) is often much too low anyway and regularization becomes necessary.

For performing maximum likelihood or MAP estimation using channels, a suitable algorithm for extracting the maximum of the represented distribution is required. For cos 2 -channels with a spacing of π 3a , an optimal algorithm in least-squares sense is obtained in the one-dimensional case as [START_REF] Forssén | Low and Medium Level Vision using Channel Representations[END_REF]

xk1 = l + 1 2a arg   l+2 j=l w j k1|k1 exp(i2a(j -l))   . (14)
N -dimensional decoding is obtained by local marginalization in a window of size 3 N and subsequent decoding of the N marginals. The index l of the decoding window is chosen using the maximum sum of a consecutive triplet of coefficients.

Channel-Based Tracking

This section contains the multi-dimensional channel-based tracking algorithm and the learning method for the system and the measurement models.

Channel-Based Tracking Algorithm

When plugging channel representations into ( 5) -( 7), the two conditional densities p(x k |x k-1 ) and p(z k |x k ) have to be considered more closely. The power theorem which has been used to derive ( 8) and [START_REF] Zemel | Probabilistic interpretation of population codes[END_REF] does not hold if we sample with channels instead of impulses, because some high-frequency content is removed and thus the scalar product between sampled densities will always be less than the integral product of the corresponding continuous densities. However, if the densities are band-limited from the start, the regularization by the channel basis functions removes no or only little high-frequency content. Hence, the power theorem holds approximately, i.e., the difference between the scalar product and the integral product can be neglected. This has been confirmed in experiments with synthetically generated densities. Henceforth, ( 8) and ( 10) can be applied for the channel-based density representations as well.

For what follows, the coefficients of ( 8) are summarized in the matrix F k = {f ij k } and the coefficients of ( 10) are summarized in the vector-valued function h k (z k ) = {h j k (z k )}. In the following section, both operators will be learned from a set of training data. This requires that both remain stationary and we remove the time index k (not from z k though): F and h(z k ). This removes absolute time reference from the model. The posterior density is now obtained by

w k|k-1 = Fw k-1|k-1 (15)
w k|k = h(z k ) • w k|k-1 h T (z k )w k|k-1 , ( 16
)
where • is the element-wise product, i.e., the enumerator remains a vector. We have not yet addressed the state space and the prediction and measurement model. The models will be treated separately in the subsequent section, for the moment it is assumed that they are known. The state space should contain sufficiently many degrees of freedom to describe the observed phenomena. Typical choices are first or second order Euclidean motion states, since the Markov model of order one applied in (3) only keeps track of the previous state.

The choice of such a complex state space requires prior knowledge about the problem which is not available in the considered task, as mentioned in Sect. 1.1. Applying the Markov theorem the other way round, an equivalent setting is to use only positions as states and to apply a higher-order Markov model. Hence, more than just the previous state is considered in the prediction step:

w k|k-1,k-2,...,k-n = [F 1 F 2 . . . F n ][w T k-1|k-1 w T k-2|k-2 . . . w T k-n|k-n ] T (17)
In the observation step [START_REF] North | Learning dynamical models using expectation-maximisation[END_REF] the prior pdf has to be changed accordingly.

inria-00321934, version 1 -16 Sep 2008

In [START_REF] Snippe | Discrimination thresholds for channel-coded systems[END_REF], the state space is represented by the concatenation of channel representations. Using concatenation of channels instead of an outer product of channel vectors leads to linear asymptotic growth of computational complexity instead of exponential growth. Using a concatenation of channel vectors corresponds to a marginalization of the joint density of previous states. In theory this could lead to predictions based on non-corresponding previous states, but following the line of arguments in [START_REF] Jonsson | Correspondence-free associative learning[END_REF], this is very unlikely to happen for channel based linear mappings. The whole channel-based tracking algorithm is summarized in Algorithm 1. The mentioned prior is treated in the following section.

Learning of System and Measurement Models

A particular feature of channel-based tracking is that the system model f and the measurement model h can easily be learned -which is different from most particle-filter based methods which need pre-specified models. The system model is trained from a sequence of states, i.e., at some time the system needs access to the entire state space. This can be thought of as a calibration phase in the case of computer vision applications or as a bootstrapping phase in the case of cognitive robotics. Another option is to observe the entire state space of another, accessible system, e.g., observing the own car, and using the model to predict an inaccessible system, e.g., another car.

In all subsequent experiments, the system model is trained by estimating the matrix [F 1 F 2 . . . F n ] from the covariance of the state channel vector at time k and the n previous instances. Since the model matrix corresponds to the conditional pdf and not to the joint pdf, the covariance is normalized with the marginal distribution for the n previous states (see also [START_REF] Yakowitz | Nonparametric density estimation, prediction, and regression for markov sequences[END_REF], plugging (3.3) into (2.7))

[ F1 F2 . . . Fn ] = Kmax k=n w k|k [w T k-1|k-1 w T k-2|k-2 . . . w T k-n|k-n ] 1 Kmax k=n [w T k-1|k-1 w T k-2|k-2 . . . w T k-n|k-n ] (18)
where 1 denotes a one-vector of suitable size and the quotient is evaluated pointwise. An important advantage of the covariance-based method for estimating the model matrix is that one can easily incrementally update the matrix by adding the covariance of the updated new state and the n previous states to the enumerator and the marginal to the denominator.

Algorithm 1 channel-based tracking algorithm.

Require: [F1F2 . . . Fn], h(z k ) and prior are known For the first n steps, no complete previous state sequence is available and the prediction equation cannot be computed using the model matrix above. Instead, the empirical distribution for the first n states is stored and used as a prior in line 1 of Algorithm 1.

1: set [w T n-1|n-1 w T n-2|n-2 . . . w T 0|0 ]
The measurement model is the most difficult part in the processing chain. Since no analytic formulations of the measurement equation are available, h(z k ) cannot be a continuous function in z k and has to be approximated by a suitable scheme. For the sake of simplicity, the observation data is also represented in a channel representation such that the measurement model becomes a product of the current observation channel vector v k and a matrix H

h(z k ) ≈ recode(v T k H). (19)
Here, recode(•) denotes the subsequent decoding ( 14) of modes, pruning, and weighted re-encoding of modes, which corresponds to an inhibition of sidemaxima. The matrix H is estimated in a similar way as the system model

Ĥ = Kmax k=1 v k w T k|k 1 Kmax k=1 w T k|k ( 20
)

Experiments

In this section, experiments validating the concept of channel-based Bayesian tracking are discussed. In Sect. 4.1 we evaluate our presented method on the classical Carlin's experiment. We compare our result to state-of-the-art methods such as the SIR particle filter and the likelihood particle filter and show that we can achieve competitive results. We demonstrate the validity of our method on a real world visual experiment in Sect. 4.2.

Carlin's Experiment

In the first experiment, the following system is considered

x k = x k-1 2 + 25x k-1 1+x 2 k-1 + 8 cos(1.2k) + v k-1 z k = x 2 k 20 + n k (21)
where v k-1 and n k is zero mean Gaussian white noise with variances 10.0 and 1.0 respectively. This system is highly nonlinear regarding both the system and measurement equation and the symmetric nature of the measurement equation poses a challenging task. This example has been used for evaluation in several publications e.g. [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF][START_REF] Carlin | A Monte Carlo approach to nonnormal and nonlinear state-space modeling[END_REF][START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. We use the Root Mean Squared Error (RMSE) as a measure of performance to be able to compare our result to the ones reported in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. In our setup the initial state x 0 is set to 8 according to x k = 0 for k < 0. The true state and the estimated result for one evaluation can be seen in Fig. 1. Note that for most of the time the results are remarkable good. However for a few instances, e.g. k = 23, the estimated state seems to be mirrored in x = 0. This can be explained by the measurement function. From the measurement function alone, there is no way to tell if we are in -x or +x.

A comparison with the results in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] is shown in table 1. The RMSE presented for our method was obtained with a second order model using 12 cos 2 -channels for both the observation and state space. We trained our method on 100 sets and performed 100 evaluation runs. It should be noted that we did model the additive time dependent component of the system, i.e. the cos(1.2k) term, since it does not comply with the stationary assumption needed for learning. However, we did learn the remaining components of the system model as well as the full observation model. The likelihood particle filter performs best. It should be noted that the test scenario is heavily biased toward the particle filters, since the only unknown components for the particle filters are the different noise components, while the system and the measurement model are fully known, which in a real world scenario must be considered unrealistic due to model errors. Our results incorporate system model and measurement model errors as well as noise errors. Despite these facts, we still manage to produce competitive result. The scenario consists of a spherical object, an orange, attached to the roof by a string. We captured images of the object by an uncalibrated stereo camera setup where external and internal parameters are unknown and lens distortions are not compensated. 1 The observations z k are the 2D position of the object in the right image and the states x k are the position in the left image. We used a third order model with 17 cos 2 -channels for both the observation and state space and used 37 frames for training. The resulting F matrix is visualized in Fig. 3.

At frame 38 we simulate a sensor failure and from here on rely on channel-based tracking. The tracking result 14-16 frames after the sensor failure can be seen in Fig. 2. The entire sequence is available at http://www.diplecs.eu/publications. The continuous formulation requires re-evaluation of all kernels for each new sample. 3 Nothing is said in the paper about computational complexity / performance.

Conclusion

In this paper, a novel variant of Bayesian tracking has been proposed: channelbased tracking with learned models. The approach is related to grid-based methods, but uses smooth, overlapping bins and the system and measurement models are acquired through learning. A number of advantages have been postulated and a standard experiment and a vision experiment have been presented.

The most important advantage of channel-based tracking compared to particle filters is that competitive results are achieved while the system model and measurement model are learned from given state and observation data. In Carlin's experiment, channel-based tracking performs similarly well as particle filters and significantly better than grid-based methods, both concerning the number of bins and the accuracy. However, the accuracy is slightly lower than for stateof-the-art particle filter methods, which had to be expected, since the particle filter methods make use of analytic system and measurement models.

In a second, qualitative computer vision experiment it has been shown that channel-based tracking can be applied to learn the mapping from uncalibrated cameras to dynamic object states. This type of experiment can be considered as a prototype for estimation problems under partial occlusion or sensor failure. Other properties in relation to cited work is summarized in table 2.

Fig. 1 .

 1 Fig.1. The true (solid line) and estimated state (dashed line) for the Carlin's experiment.

Fig. 3 .

 3 Fig. 3. Visualization of the F-matrix that has been learned during the first 37 frames. Note the similarity to Lissajous figures of order (1,1).

Table 1 .

 1 RMSE obtained on Carlin's experiment. Our method is third in performance even though we do not model the system or measurement model. All results except for CBT (12 channels) were taken from[START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] (50 particles and 50 grid points).

	Algorithm	RMSE
	Extended Kalman filter	23.19
	Approximate grid-Based Filter	6.09
	Regularized Particle filter	5.55
	SIR Particle filter	5.54
	Channel Based Tracking	5.43
	Auxiliary Particle filter	5.35
	Likelihood Particle filter	5.30

Table 2 .

 2 Comparison to cited work. CBT refers to channel-based tracking and [3] to the four particle filter methods from table 1.

	Method	[13]	[14]	[15]	[16]	[3] CBT
	learned models	system system	no	system no	both
	multi-modal densities	no	yes	yes	no	yes	yes
	fast implementation	no 2	no 2	no 3	yes	yes	yes
	2

Actually, offline data has been used consisting of two image sequences only. One of the sequences contained two frame-drops, which were automatically discovered by the channel-based tracking.

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n • 215078 (DIPLECS) and from the CENIIT project CAIRIS.


Jia-Bin Huang

Sih-Ying Wu

Chu-Song Chen

Enhancing Color Representation for the Color Vision Impaired

In this paper, we propose a fast re-coloring algorithm to improve the accessibility for the color vision impaired. Compared to people with normal color vision, people with color vision impairment have difficulty in distinguishing between certain combinations of colors. This may hinder visual communication owing to the increasing use of colors in recent years. To address this problem, we re-map the hue components in the HSV color space based on the statistics of local characteristics of the original color image. We enhance the color contrast through generalized histogram equalization. A control parameter is provided for various users to specify the degree of enhancement to meet their needs. Experimental results are illustrated to demonstrate the effectiveness and efficiency of the proposed re-coloring algorithm.

Introduction

Due to the availability of color printers and color display devices, the use of colors in multimedia contents to convey rich visual information has dramatically increased. It becomes more important to perceive colors for effective visual communication. However, people with certain types of color vision impairment have difficulty in distinguishing between some colors. In Fig. 1 we show examples of how color vision impaired people perceive color images. The two images in the left are the original color images perceived by people with normal vision, and the visual information can be easily interpreted. On the other hand, important color information in the original images may disappear or become indistinct in the six images in the right, which are the simulation results for people with different types of color vision deficiency (CVD).

In this paper, our aim is to resolve the potential inconsistency of visual perception between people with normal vision and people with CVD. We develop a fast algorithm to re-color images so that the original color contrast can be well preserved for color vision impaired viewers. Also, a user-specified parameter is provided for the trade-off between the enhancement degree and the naturalness in the re-colored images. We propose to re-map the hue components in the original color image through a spatially invariant global color transformation. The contrast is enhanced by forcing the confusing colors to have wider dynamic ranges in the hue space. The luminance and the saturation components remain unaltered to ensure the naturalness of the re-colored image.

The remainder of this paper is structured as follows. Section 2 describes human visual color perception of colors and various types of CVD. In Section 3, we review previous works on re-coloring images for accommodating the color blindness. The proposed re-coloring algorithm is introduced in Section 4. Section 5 presents the experimental results, and Section 6 concludes this paper. The energy received by the L, M, and S cones can be computed by a numerical integration over the wavelength λ:

[L, M, S] = E(λ)[ l, m, s]dλ, (1)
where E(λ) is the power spectral density of the light and l, m, s are the fundamental spectral sensitivity functions for L-, M-, and S-cones.

Color vision impairment, or color vision deficiency, results from partial or complete loss of function of one or more types of cone cells. There are three major types of CVD: anomalous trichromacy, dichromacy, and monochromacy, which are summarized in Table 1.

Anomalous trichromacy, a mild color deficiency, is often characterized by a shift of one of the three fundamental spectral sensitivity functions so that the pigments in one type of cones are not sufficiently distinct from the pigments in another type. The three types belonging to this category are protanomaly, deuteranomaly, and Tritanomaly, depending on deficiency (shift) in L-, M-, or S-cones.

A more severe color deficiency, dichromacy, is present when one of the three cone types is absent. Protanopia and deuteranopia, lacking L-and M-cones respectively, are unable to distinguish between colors in the green-yellow-red section of the spectrum, while those without S-cones are tritanopia that have difficulty in discriminating blue from yellow.

The last class is the severest but rarest type of CVD, called Monochromacy. Monochromats are also known as "total color blindness", lacking all types of cone cells. They are unable to distinguish any colors and can perceive brightness variations only.

There have been several works that attempt to simulate color deficient vision [START_REF] Brettel | Computerized simulation of color appearance for dichromats[END_REF][START_REF] Kondo | A computer simulation of anomalous color vision[END_REF][START_REF] Meyer | Color-defective vision and computer graphics displays[END_REF]. In [START_REF] Brettel | Computerized simulation of color appearance for dichromats[END_REF], a computational model of CVD was formulated in the three dimensional LMS space, where three orthogonal axes L, M, and S represent the responses of the three different cones. The simulation algorithm computes the LMS tristimulus values from the RGB color space and projects colors in the LMS color space onto a pair of planes to represent the perceived responses by CVD viewers. The projected LMS values are then transformed to the RGB color space, allowing people with normal vision to experience how the CVD viewers perceive colors.

Fig. 1 shows examples of the simulation results of CVD. The top row in the two sets of images are the simulation results of anomalous trichromacy. Images listed from left to right are protanomaly, deuteranomaly, and tritanomaly. Similarly, the bottom row are results of dichromacy. Images listed from left to right are the results of protanopia, deuteranopia, and tritanopia. The first image set is an image from Ishihara test chart. We can see that the number "45" disappears for protanopia and deuteranopia and become hard to recognize for protanomaly and deuteranomaly. The second image set contains blue-yellow content, which can not be distinguished by tritanopia.

With the ability to simulate color perception of CVD viewers, we can perform color transformation such that the color information in the original color images can be preserved for color-deficient users.

Related Work

Many previous works have been devoted to address CVD accessibility. We classify them into two main categories: 1) tools that provide guidelines for designers to avoid ambiguous color combinations, and 2) methods that (semi-)automatically reproduce colors that are suitable for CVD viewers. Methods in the first category assist designers in color reproduction by providing guidelines [4], using a restricted CVD palette [START_REF] Rigden | The Eye of the Beholder-Designing for Colour-Blind Users[END_REF] [6], or verifying color schemes [START_REF] Walraven | Color displays for the color blind[END_REF] [START_REF] Jenny | Designing Maps for the Colour-Vision Impaired[END_REF]. However, it still takes a lot of effort for designers to select colors that are friendly for CVD viewers. Also, these methods can not be applied to existing natural images, which may contain tens of thousands of colors.

Based on the observation that red-green colorblind people (those with protanopic and deuteranopic deficiencies) are unable to distinguish between red and green, Daltonization [START_REF] Dougherty | Daltonize[END_REF] presents a procedure to allow users to specify three parameters for re-coloring images. The first parameter is for stretching contrast between red and green hues, the second one is for blue and yellow contrast modulation, and the last one is for modulating luminance component. The results are sensitive to the selection of the parameters. Improper parameters may result in unnatural images. Visual content adaptation for CVD viewers in the MPEG-21 digital item adaptation framework had been proposed in [START_REF] Song | Digital item adaptation for color vision variations[END_REF]. The adaptation for dichromats aims to provide better information by adjusting hue and saturation in the HSI color space. Ichikawa et al. [11] [12] proposed a more general method. First they used color quantization to select a subset of colors from the images, and then build an objective function that maintains the color distances and the extent of re-coloring. They solved the optimization problem by a genetic search algorithm. The final result was obtained using an interpolation method. However, the lack of considering luminance consistency and the gamut constraint may raise some problems. Recently, Rasche et al. proposed automatic methods for re-coloring images using linear transformation [START_REF] Rasche | Detail Preserving Reproduction of Color Images for Monochromats and Dichromats[END_REF] or constrained multidimensional scaling [START_REF] Rasche | Re-coloring Images for Gamuts of Lower Dimension[END_REF]. To alleviate the computational burden, 256 landmark colors are selected and used to minimize a pre-defined objective function, which preserves the proportional color differences. Unlike previous approaches focusing on re-coloring images only for CVD users, Huang et al. [START_REF] Huang | Information Preserving Color Transformation for Protanopia and Deuteranopia[END_REF] aimed to maintain the naturalness of the original images as well as preserving the color contrasts . They constrained the color mapping as a rotation operation in the a*b* plane in the CIELab color space and introduced two objective functions: a detail term for CVD viewers and a naturalness term for normal views. A parameter was provided to control the trade-off between these two terms. In [START_REF] Wakita | SmartColor: Disambiguation Framework for the Colorblind[END_REF], Wakita et al. introduced three mapping errors to preserve the color information that the authors expect to convey in a document. They presented that color effects can be modeled as three types: color contrast, distinguishing ability, and natural coloring. However, the optimization process is very computational demanding and can only handle a small number of key colors (usually no more than 10). This prohibits the application for re-coloring natural images. Jefferson and Harvey [START_REF] Jefferson | Accommodating color blind computer users[END_REF] presented a framework for re-coloring color documents similar to [START_REF] Ichikawa | Preliminary study on color modification for still images to realize barrier-free color vision[END_REF], [START_REF] Rasche | Re-coloring Images for Gamuts of Lower Dimension[END_REF], and [START_REF] Huang | Information Preserving Color Transformation for Protanopia and Deuteranopia[END_REF]. The novelties of their method lie in using the W3C color evaluation criteria, mapping without the gamut problem, and developing a new method for selecting key colors. Later in [START_REF] Jefferson | An interface to support color blind computer users[END_REF], they proposed an interface to assist CVD users to access digital images. The main idea is to transfer the color information of the defective cones to the working ones One major disadvantage of the previous re-coloring methods through optimization approaches is the expensive computational cost. Moreover, the optimal results for the pre-defined objective functions may not necessarily yield visual appealing images for CVD viewers due to the luminance consistency or gamut mapping problems.

4 The Proposed Re-coloring Algorithm

Design Objectives

To address the problems stated in Section 3, we propose a fast approach that automatically enhances the color contrast for CVD viewers. Listed below are the three design objectives of our method:

1. Maintaining luminance and saturation consistency:

Our method aims not to enhance the contrast in the original image, but to maintain the contrast that may not be perceived by CVD users. The consistency of luminance and saturation is important to generate a natural image for people with normal vision.

Preserving the order of the hue:

By "the order of hue" we mean the ordering in the hue component in the HSV color space. In other words, we wish not to corrupt the natural hue ordering of the original image, but to stretch or compress the distance between them.

Fast computation:

We aim to design a fast algorithm that is applicable to real-time applications.

Since the luminance and saturation components are often not the major factors that lead to confusion for CVD users, we choose the simplest way to preserve them: leave them unchanged. An RGB color images is first converted to the HSV color space and then we apply a color transfer function to the hue component only. To preserve the natural ordering of hue, the hue transfer function is constrained to be a non-decreasing function. In the next subsection, we show how the hue transfer function can be obtained without performing a time-consuming optimization procedure. After transferring the hue values, the re-colored image is converted back to the RGB color space for display.

Re-coloring using Generalized Histogram Equalization

We enhance the contrast in the hue channel through a histogram transformation approach. Histogram equalization (HE) is one of the most well-known techniques for image enhancement. HE enhances the contrast of an image by expanding the dynamic range of the original image such that the histogram of the resultant image is better distributed. To this end, the cumulative density function is used for intensity mapping.

The intensity transfer function is found by first generating the histogram of the image, and then the mapping function can be written as:

T (g) = g min + (g max -g min ) g g min hist(g) dg, (2)
where the hist(g) is the normalized histogram, g is the original intensity value, and g max and g min are the maximum and minimum intensity values of g, respectively. The normalized histogram hist(g), that is, the probability distribution of the grey levels in an image, can also be interpreted as an expansion function, where a larger value of hist(g) tends to stretch the range of grey levels around g. Since the histogram only contains the global statistics (the occurrence frequency of each grey level) of an image, unexpected effects like over-enhancement may occur in some regions because local characteristics are not considered. Therefore, we adapt the generalized histogram equalization (GHE) [START_REF] Jen | Generalized Histogram Equalization Based on Local Characteristics[END_REF] to our re-coloring algorithm.

The generation of the histogram in HE is a masking-and-accumulating operation with 1 × 1 window. In the GHE, the mask is generalized from 1 × 1 to n × n to encode the local information into the histogram. We use a 3 × 3 mask for simplicity, and in practice larger masks like 5 × 5 or 7 × 7 do not affect the results significantly.

We measure three local characteristics α, β, and γ within the neighborhood system N (x, y) centered at point (x, y). The first measurement α is the hue value at point (x, y), which is the only measurement in HE: α(x, y) = hue(x, y).

(

) 3
The second one β measures the maximum local hue difference within the neighborhood system N (x, y):

β(x, y) = max i,j {hue(i, j)} -min i,j {hue(i, j)}, (i, j) ∈ N (x, y). ( 4
)
Obviously, more complicated measurement like local hue variance can be used as well. The last one evaluates the local color information loss due to CVD in the CIELab color space, where the Euclidian distance corresponds to human visual perception of color difference:

γ(x, y) = (i,j)∈N (x,y) (||(C(x, y) -C(i, j))|| -||(Sim(C(x, y)) -Sim(C(i, j)))||) 2 ,
(5) where C(x, y) is the CIELab color at point (x, y), || • || is the Euclidean norm, and Sim(•) is the simulation CIELab color for CVD viewers [START_REF] Brettel | Computerized simulation of color appearance for dichromats[END_REF].

The third measurement γ evaluates how much local color contrast is lost when the image is perceived by CVD viewers. One can expect that γ has a large value in regions where color confusion occurs. We show in Fig. 2(a) a natural image with red and green colors. In Fig. 2(b), the color contrast between red and green is lost when perceived by protanopia. The expansion function from a pixel in HE is a dirac delta function δ(x -g), while in GHE, the expansion function can be conditional on the measured local characteristics, and common kernel functions such as Gaussian, Rectangle and Epanechnikov can be used. With the computational efficiency in mind, we choose the rectangle kernel as our expansion function:

S(h|α, β, γ) = γ × Rect( h -α β ), (6)
where the Rect(•) is the rectangle kernel which can be expressed as:

Rect(x) ≡ 1, if -0.5 ≤ x ≤ 0.5, 0, otherwise. (7)
After scanning over the input color image, calculating the local characteristics, the generalized histogram GH(h) can be obtained:

GH(h) = S(h|α(x, y), β(x, y), γ(x, y))dxdy, (8)
Similar to HE, we can then construct the hue transfer function T (h):

T (h) = h min + (h max -h min ) h hmin GH(h)dh, (9)
where h denotes the input hue value. Note that the transfer function is always non-negative and non-decreasing, and thus can yield more visual pleasing images. The h min is not altered after the hue transfer function h min = T (h min ). Since the hue axis is circular, users can specify a pivot hue value so that this hue value will not be changed after enhancing. In this paper, h min is set to be zero.

Controlling the Degree of Enhancement

Although the hue transfer function can be automatically computed from the given image, we provide a parameter for users to control the degree of enhancement to satisfy various needs of different viewers. This is achieved by introducing a magnitude mapping function:

M (x) = x p , ( 10
)
where p is a user-specified control parameter. The transfer function is further modified as:

T (h) = h min + (h max -h min ) h hmin M (GH(h))dh h max hmin M (GH(h))dh . ( 11
)
If p is larger than one, the normalized expansion function is emphasized, resulting in more significant enhancement. If users prefer to generate a more natural image, p < 1 can be specified so that the dynamic range of the normalized expansion function is compressed. Note that if p = 0, the transfer function maps every input value to itself, generating the original image.

We have implemented and tested our algorithm using C++ on a Pentium 4 3.4GHz PC. The complexity of our algorithm scales linearly with the size of the input image. Re-coloring an image with 200x200 pixels takes less than 0.5 second, which is significantly faster than previous published works (a few minutes or more [START_REF] Jefferson | Accommodating color blind computer users[END_REF] [13] [START_REF] Rasche | Re-coloring Images for Gamuts of Lower Dimension[END_REF]). Since it is difficult to gather various types of color-deficient viewers to evaluate our method, we use the computational model to simulate color perception of people with CVD. We use the images from the Ishihara test chart to demonstrate the effectiveness of the proposed re-coloring method. Ishihara test chart is a book that consists of a number of cards for discovering congenital color blindness and red-green blindness. Two samples from Ishihara test chart are shown in the left of Fig. 3. People with normal vision can easily identify the number "6" and "8" in the original images. However, people with dichromacy are unable to see the number from the test chart (simulated in Fig. 3(a)(b)), while the numbers are very unclear to those with anomalous trichromacy (simulated in Fig. 3(e)(f)). We apply our re-coloring algorithm to these two test charts. The re-colored and simulated results for different types of CVD are shown in Fig. 3 (c)(d)(g)(h). After re-coloring, we can see that the numbers appear (for dichromacy) or become clearer (for anomalous trichromacy) for CVD viewers.

We compare our algorithm with that proposed by Rasche et al. in [START_REF] Rasche | Re-coloring Images for Gamuts of Lower Dimension[END_REF]. The left part of Fig. 4 is the simulation result for protanopia and the right part are for deuteranopia. Fig. 4 (a) (d) show the simulation results for protanopia and deuteranopia, respectively. The re-colored and simulated images by our method are shown in Fig. 4 (b) (e). Fig. 4 (c) (f) show the results by Rasche et al. Their method did generate images with high contrast, but may yield over-enhanced images because they did not take hue order preserving into consideration. Note also that, their method maps the green to blue in the upper image "berries" and maps the red to blue in the lower image "flower". This causes no problems in still images, but may result in corruption of temporal coherence in videos.

Fig. 5 shows the comparison with method proposed by Jefferson et al. in [START_REF] Jefferson | Accommodating color blind computer users[END_REF]. Fig. 5(a) are the two original color image, which contain coloring pencils. The simulated images for deuteranopia (in the left) and tritanopia (in the right) are illustrated in Fig. 5(b). We can see that the color contrast in the original image disappears when perceived by CVD users. Fig. 5(c) (e) show the re-colored images by our approach and by Jefferson et al.'s method, respectively. The corresponding simulation results of the re-colored images are shown in Fig. 5(d) (f). The images re-colored by Jefferson et al.'s method are unnatural, and the luminance of the re-colored images are inconsistent with the original ones. Compared with their method, our re-coloring algorithm produces more perceptually pleasing images.

We also show the effects caused by the control parameter p. A natural image "flower" and an image taken from Ishihara test chart are re-colored using different parameter p. Fig. 6(a) (b) show the original images and their simulation results for deuteranopia. From Fig. 6(c) to (h), the control parameter p is set inria-00321936, version 1 -16 Sep 2008 from 0.2 to 1.2, with an increasing step 0.2. As p increases, the simulated image exhibits higher contrast, which is particularly obvious in the test chart images.

(a) (c) (e) (g) (b) (d) (f) (h) (a) (c) (e) (g) (b) (d) (f) (h)

Conclusions

In this paper, we have proposed a fast and effective re-coloring algorithm for people with color vision impairment. We propose to enhance the contrast in the hue channel only, which is the major factor resulting in confusion in visual communication. We employ the generalized histogram technique to encode local information in the original color image. We also provide a control parameter for users to specify the degree of enhancement. Our algorithm can run at real-time, and thus can be easily extended to video processing applications. Computational simulation results have demonstrated the effectiveness and efficiency of the proposed algorithm.

Fig. 1 .

 1 Fig. 1. Original color images and the simulation results for various types of color blindness. Images in the left are the original color images perceived by people with normal vision, while the six images in the right are the simulation results. The top row are the results for anomalous trichromacy and the bottom row are for dichromacy. From left to right, the first column is for protanopic viewers, the second is for deuteranopic viewers, and the third is for Tritanopic deficiency.

Fig. 2 (

 2 c) illustrates the map of γ values in logarithmic scale. The brighter regions correspond to larger γ values.

Fig. 2 .

 2 Fig. 2. The information loss map in image "berries". (a) The original color image. (b) The simulation result for protanopia. (c)The local color contrast loss map in logarithmic scale.

Fig. 3 .

 3 Fig. 3. Enhancing color contrast for red-green types of CVD on images from Ishihara test chart. (a)(b) The simulation results for protanopia and deuteranopia, respectively. (c)(d) The simulation results of the re-colored images for protanopia and deuteranopia. (e)(f) The simulation results for protanomaly and deuteranomaly, respectively. (g)(h) The simulation results of the re-colored images for protanomaly and deuteranomaly.

Fig. 4 .Fig. 5 . 2 Fig. 6 .

 4526 Fig. 4. Comparison with Rasche et al.'s method using two images "berries" and "flower". The left part are the simulation results for protanopia and the right image set is the simulated images for deuteranopia. (a)(d) Simulated images. (b)(e) Simulation results of the re-colored images by our method (with p = 0.6). (c)(f) Simulation results of the re-colored images by our method by Rasche et al.'s method.

Table 1 .

 1 Major genetic color deficiencies and prevalence for western races.

	Type	Name	Cause of defect	Prevalence
		Protanomaly	L-cone defect	1.3%
	Anomalous trichromacy	Deuteranomaly	M-cone defect	4.9%
		Tritanomaly	S-cone defect	0.01%
		Protanopia	L-cone absent	1%
	Dichromacy	Deuteranopia	M-cone absent	1.1%
		Tritanopia	S-cone absent	0.002%
	Monochromacy	Rod Monochromacy no functioning cones very rare

Acknowledgement

This research was supported by the National Science Council of Taiwan under Grant No. NSC 96-3113-H-001-011.