Removing Flicker from Old Movies

François PITIÉ

SIGMEDIA lab Dept. of Electronic & Electrical Engineering Trinity College, Dublin

Supervisor Dr Anil Kokaram Dr Rozenn Dahyot head of the SIGMEDIA lab [email protected]. [email protected] Submitted as part fullfilment for Diplôme d’ingénieur for ENST Bretagne/Eurecom, FRANCE DEA Image-vision for the Université de Nice Sophia-Antipolis, FRANCE

Contents 1

Introduction

3

2

Flicker Modeling

4

3

Detection of Intensity Flicker in video sequences 3.1 Flicker Analysis . . . . . . . . . . . . . . . . . . . . 3.1.1 Using histograms . . . . . . . . . . . . . . . 3.1.2 Using a mapping function f . . . . . . . . . 3.2 Detection Measures . . . . . . . . . . . . . . . . . . 3.2.1 1st Measure: difference of the mean intensity 3.2.2 2nd Measure: Detection of main shift lobe . . 3.3 Localization of the flicker . . . . . . . . . . . . . . . 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

6 6 6 7 7 7 8 9 9

. . . . . . . . . . . . . .

10 10 10 12 12 13 13 14 14 15 15 16 18 19 20

4

Flicker Estimation 4.1 Overview of the existing methods . . . . . . . . 4.1.1 Parameters estimation . . . . . . . . . . 4.1.2 Flicker Localization . . . . . . . . . . . 4.1.3 Ohuchi’s method [12] . . . . . . . . . . 4.2 Robust Estimation . . . . . . . . . . . . . . . . . 4.2.1 Matrix formulation . . . . . . . . . . . . 4.2.2 Energy of the residuals . . . . . . . . . . 4.2.3 Iterative Reweighted Least Square (IRLS) 4.2.4 Scale parameter . . . . . . . . . . . . . . 4.3 Results: evaluation of the estimator robustness . . 4.4 Cosine Basis vs. Polynomial Basis . . . . . . . . 4.5 Energy and Symmetry . . . . . . . . . . . . . . 4.6 Regularization . . . . . . . . . . . . . . . . . . . 4.7 A new Estimator . . . . . . . . . . . . . . . . . 2

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . .

4.8 5

6

7

4.7.1 A new model for the flicker . . . . . . . . . . . . . . . . . . . . . . . 4.7.2 The restoration process . . . . . . . . . . . . . . . . . . . . . . . . . . Results of the new Estimator . . . . . . . . . . . . . . . . . . . . . . . . . . .

Flicker Compensation 5.1 Overview of the existing methods . . . . . . . 5.1.1 A Simple Method . . . . . . . . . . . . 5.1.2 Roosmalen’s Restoration . . . . . . . . 5.2 Smoothing the integrated parameters . . . . . . 5.2.1 The problem . . . . . . . . . . . . . . 5.2.2 Formulation of the problem. . . . . . . 5.2.3 Removing the noise . . . . . . . . . . . 5.2.4 Nature of the noise: choice of the filter 5.2.5 Conclusion . . . . . . . . . . . . . . . IFA Movie 6.1 Irish Film Archive . . . 6.2 The movie . . . . . . . 6.3 The restoration process 6.4 Some results . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

20 20 21

. . . . . . . . .

22 22 22 22 23 23 23 25 25 26

. . . .

27 27 27 28 30

Conclusion and future work.

31

Bibliography

34

A Energy & Symmetry A.1 The regression problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Looking at the residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35 35 35

B Rory O’More: the rising of 1641 in Ireland

38

3

Abstract The video archives of the world contain many important historic, artistic and cultural records that are stored in a fragile state, on a volatile medium. They require conservation and restoration. Automated tools for video restoration will be crucial in preserving our cultural heritage, since manual image restoration is a tedious and time-consuming process. This work concerns digital flicker removal in old corrupted footage and was undertaken as part of the European Union project BRAVA. The effect of flicker is a common degradation in old movies. It is perceived as a fluctuation in the brightness of the frames. The aim of this work was to develop and implement algorithms for the detection, estimation and then removal of flicker. We proposed some important improvements to the classical algorithms and we tested our new methods on a real movie of 12500 heavily degraded frames.

Les archives vidéo contiennent de nombreux chefs-d’oeuvre artistiques qui sont actuellement stockés sur de fragiles supports et qui doivent être impérativement conservés et restaurés. Pour préserver ce patrimoine culturel qu’il est impensable de réparer à la main, il est crucial d’élaborer des outils automatiques de restauration vidéo. Ce travail - qui s’inscrit dans le cadre du projet européen BRAVA - porte sur la restauration numérique dans les vieux films de la dégradation de flicker (ou bruit de battement ou pompage). Il s’agit d’une dégradation classique qui se traduit par une fluctuation de l’intensité lumineuse au cours du temps. Le but de ce projet est de développer et d’implémenter les algorithmes de détection, d’estimation et de compensation du flicker. Nous avons apporté des contributions importantes aux algorithmes actuels, ce qui nous a permis de restaurer un vieux film d’environ 12500 images particulièrement endommagées.

Acknowledgements I would like to thank at first my supervisor Dr Anil Kokaram for his guidance and his enthusiastic dynamism. I am also very grateful to Dr Rozenn Dahyot for her pertinent advices and for the time she dedicated to me. Finally I would like to thank Hugh Denman and all the other people from the lab, especially Justine Grave, Francis Kelly, Elena Ranguelova, Denis Sidorov, Ludvick Tesar, Yanming Yu.

1

Notations • In the intensity level of observed frame n • hn intensity histogram / probability density functions of In • f (u) = v mapping function between gray levels of two pictures • IR n restored frame n • IC n compensated frame n • gn,m (Im ) = In differences observed between frame n and m R • µn,m (IR n ) = Im differences wanted between frame n and m

• αn (In ) = IR n compensation function (inverse of the degradation) of frame n • θαn parameters associated with the function αn • θgn,m parameters associated with the function gn,m

2

Chapter 1 Introduction The video archives of the world contain many important historic, artistic and cultural records that are stored in a fragile state, on a volatile medium, and in need of conservation and restoration. Automated tools for video restoration will be crucial in preserving our cultural heritage, since manual image restoration is a tedious and time-consuming process. The Broadcast Archives Restoration Through Video Analysis (BRAVA) Project, initiated in 1999, is supported by the European Commission. It is coordinated by INA, in partnership with the company Snell&Wilcox, other broadcast archive owners (RTP, the BBC), Universities (TU Delft, TCD Dublin), and SGT. The project aims to ensure the widest possible access to valuable archive material by significantly enhancing the efficiency of the video and film restoration process. This report concerns digital flicker removal in old corrupted footage. This work was undertaken as part of the European Union project BRAVA. The effect of flicker is a common degradation in old movies. It is perceived as a fluctuation in the brightness of the frames. The aim of this work is to develop and implement algorithms for the detection and removal of flicker. The flicker removal framework has been divided into three parts. The first part consists of a diagnosis module capable of detecting frames affected by flicker in a footage (chapter 3). The second part is the robust characterization of the flicker through the estimation of parameters (chapter 4). The last part is the restoration of the degraded footage, given these estimated parameters (chapter 5). The proposed algorithms have been used in a practical full movie restoration process, which will be described at the end of the report (chapter 6). This work was undertaken during a 6 month internship as Master student in the Electronic and Electrical Engineering Department, Trinity College Dublin, as part of my studies at Eurecom/ENST Bretagne and toward the DEA Image-Vision at the University of Nice-Sophia Antipolis. The work was performed under the supervision of Dr Rozenn Dahyot and Dr Anil Kokaram in the team led by Dr Anil Kokaram.

3

Chapter 2 Flicker Modeling The effect of flicker is perceived as a temporally random, and spatially varying luminance change in images. This is a typical degradation observed in archived film and video. For movies shot in the silent era, this effect was due to the different exposure times of the film as it was hand cranked through the film camera. More modern material can also be affected by luminance changes due to poor lighting. The figure 2.1 shows this intensity fluctuation in a sequence. The flicker displayed on this figure is really strong and the visual perception is even worse when the movie is played in real time. For more examples see http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/shots/

Figure 2.1: Flicker effect on a sequence of the Irish movie Rory O’More (1911). Previous work [4], [15] consider that flicker can be modeled as a spatially varying gain and 4

offset function. In+1 (x, y) = a(x, y) · In (x, y) + b(x, y)

(2.1)

This equation establishes the relationship between the intensity levels of two consecutive degraded frames: In (x, y) and In+1 (x, y). This model is not based on any physical properties, but only deduced from observation and it will be improved in section 4.7.1. Sequences affected by flicker show low frequency spatial degradation. Therefore the functions of gain a(x, y) and offset b(x, y) should be low frequency. Previous work by [4], [15] employs a 2-dimensional polynomial basis. For example: a(x, y) = a0,0 + a1,0 x + ... + ai,j xi y j

(2.2)

b(x, y) = b0,0 + b1,0 x + ... + bi,j xi y j

(2.3)

The figure 2.2 shows an example of flicker following this model.

a(x, y)

b(x, y)

Original Image

Altered Image

Figure 2.2: Example of original, artificially altered frame with the corresponding representations of the functions a(x, y) and b(x, y).

5

Chapter 3 Detection of Intensity Flicker in video sequences Within a long programme of video, there may be only short segments affected by flicker. Typically an operator would view the sequence and decide where certain jobs need to be done. The aim of this part is to present a diagnosis module that should be able to measure automatically as quickly as possible the strength of the flicker in the frames. The detection has to be fast and robust: 1. Fast enough to make the process useful, 2. Robust enough to detect all the corrupted frames and limit the number of false alarms. Different measure methods are reviewed in this chapter.

3.1

Flicker Analysis

3.1.1

Using histograms

Histograms have been used in context analysis for some time. Histograms are an approximation of probability density functions (p.d.f.) of the image intensity. Working with intensity p.d.f. avoids problems linked with the position of the information. So the p.d.f. will not change in case of such noise as jitter, shaking and all other noise based on position information, as well as the p.d.f. will not change dramatically in a same shot. It is certain that flicker can affect the shape of image histograms and therefore it seems useful to consider them as a potential fast diagnostic. We decided therefore to work on p.d.f., even if we drop besides the spatial information of the frames, which is not important for the detection module.

6

3.1.2 Using a mapping function f Consider the problem at order zero in x and y. The relationship between successive frames is then: In+1 (x, y) = a · In (x, y) + b The parameters a and b are spatially constant. And so this function is applied globally to the image. For each gray level u in In we can associate the gray level v = f (u) in In+1 . v = f (u) = a · u + b Come the p.d.f. of the image h(u). We have the general equation [8]: hn (u) du = hn+1 (v) dv

(3.1)

hn (u) du = hn+1 (f (u)) df (u) hn (u) = hn+1 (f (u)) f 0 (u)

(3.2)

Considering flicker case of zero order gain and offset: hn (u) = a · hn+1 (a · u + b)

(3.3)

The figure 3.1 shows this relation between the two histograms. Note in particular the shift between both curves (due to the offset b).

3.2 3.2.1

Detection Measures 1st Measure: difference of the mean intensity

This first method to measure the flicker strength is based on the fact that the strength of the flicker is directly linked to the intensity mean of the image. This method has been often used in literature [4, 15, 12, 11] to check the quality of the results of a flicker-removal process in a sequence since the mean intensity of a sequence without flicker is generally steady. If we note: f (u) = u + δ(u) δ(u) represents the intensity shift observed for the gray level u. Then we can write: Eu [f (u)] = Eu [u] + Eu [δ(u)] From (3.1) yields Ev [v] = Eu [f (u)] and then : Eu [δ(u)] = Ev [v] − Eu [u] 7

(3.4)

Figure 3.1: Example of histogram distortion due to flicker effect. Histograms of frames In and In+1 are plotted in the same graph.

So the mean of the intensity shift can be easily computed. Other noises than flicker usually occur in a movie. For instance dark or bright blotches are typical degradations (see section 6.3) in old movies. This method is sensible to large blotches. Furthermore, if there is a spatial flicker fluctuation (higher order model in x and y), δ(u) may be positive and negative in the same image. Therefore Eu [δ(u)] may be null and so the frame may be not detected.

3.2.2

2nd Measure: Detection of main shift lobe

In order to be more robust we can look at the content of the picture. We know that the main lobe of an intensity histogram of a function corresponds to the main objects of the image (the background for instance). The idea is so to look at the position of the main lobe of the histograms. Its position do not change too much except in presence of flicker. From (3.3) we know that the positions of the main lobe of hn and hn+1 are linked by this relation: arg max(hn+1 (u)) = a · arg max(hn (u)) + b (3.5) u

u

So if we measure the distance between the lobes positions, we will have a measure of the flicker intensity. The main problem is to detect the lobes. This can be done by taking the zeros of the derivative of the histogram function. But sometimes histograms can have several lobes of same importance, so that one lobe is not always the most important in a shot. To avoid this problem 8

we can keep not only the main lobe but also several other important lobes and try to match correspondences between histograms of consecutive frames. Nevertheless this method causes some problems of precision since histograms may be flat. In such cases, the lobe detection will not be accurate and false alarms will be generated.

3.3

Localization of the flicker

We have described two methods to detect flicker. These method were designed for zero order of gain and offset functions. But we know that the flicker may not be global. Localized flicker can be perceived by the viewer but detections processes described above may not detect these corrupted frames. To improve the results, we can split the frame into 5 blocks as described in figure 3.2. The process is then done separately on the 5 blocks. For each frame, the final measure of flicker can be taken as the maximum of the 5 measures.

Figure 3.2: To deal with the localization of the flicker, the frame is split into 5 overlapping blocks.

3.4

Results

Unfortunately we had not enough time to really test further the methods described. The code of the methods is written and can be reused. However a video showing results of the maximum lobe shift detector is available on http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/detection/demo. avi.

9

Chapter 4 Flicker Estimation Having detected flicker we can now concentrate on estimating the characteristics of the degradation to be able to remove it. The aim of the chapter is to estimate the parameters involved in the flicker model. We will present at first in section 4.1.1 the existing methods to estimate the flicker. Then in section 4.2 we will show how a robust estimation of the flicker parameters can be done. Then in section 4.5 and 4.4 we will present some major improvements to the current methods. Eventually in section 4.7 we will present a new model for the flicker as well as a new estimator.

4.1 Overview of the existing methods This section presents a brief overview of the estimation methods that have been proposed in the literature.

4.1.1

Parameters estimation

Decencière’s method [4]. The method proposed by Decencière in his thesis is based on histogram matching. He uses a zeroth order model for the flicker: In+1 = a · In + b with a and b spatially constants for each frame. Assuming that the maximum and minimum gray values have the same values for the frames

10

of a same shot, we have:

    max(In+1 ) = a · max(In ) + b min(In+1 ) = a · min(In ) + b    E[In+1 ] = a · E[In ] + b

Parameters a and b can be easily retrieved from these equations. Roosmalen’s method [15]. The zeroth order parameters are found by comparing the moments (variance and mean) of the 2 images. s var[In+1 (x, y)] a= ; b = E[In+1 (x, y)] − a · E[In (x, y)] var[In (x, y)] This method corresponds to the least mean square solution of the problem. This method is a popular method. It has been used by [17, 12] and the company Snell&Wilcox had for instance based their flicker removal on this algorithm. Roosmalen had also proposed some other improvements for the flicker removal (cf. sections 4.1.2 and 5.1.2) Naranjo’s method [11] can be used at zeroth order in x and y. His method is based on histogram matching. The histogram of the degraded frame is transformed into a target histogram which corresponds to the expected histogram of the original frame. The idea is to estimate the function f such that: f (In ) = In+1 Finding f can be performed by using standard method of histogram matching. This method uses the cumulated histograms of both images [11, 8] (Cn and Cn+1 ). Z

u0

hn (u) du = hn+1 (v) dv Z v0 hn (u) du = hn+1 (v) dv

u=0

0

Cn (u) = Cn+1 (f (u)) Since the cumulated histograms functions are monotonous it is possible to compute their inverse. Then we have: −1 f (u) = Cn+1 ◦ Cn (u) (4.1) The figure 4.1 shows an example of mapping function between two degraded frames. We can see for this particular example that f is actually not affine but curved. But this method is quite sensitive to outliers. A big dark blotch can for instance produce a lobe in the histogram and make other parts of the curve smaller. Then applying 4.1 will shift all the correspondences in the inversion of the cumulated histograms. 11

Figure 4.1: Example of intensity level matching function f .

4.1.2 Flicker Localization All the previous methods consider only global flicker (zero order in x and y). But we know that the flicker may have some spatial fluctuations. To deal with this aspect, Roosmalen [15] split the image into small overlapping blocks. As the flicker is assumed to be spatially low frequency, it is possible to use his zeroth order estimator on these blocks. But the estimation of parameters a and b may be not reliable on some blocks, if the variance of the signal is null for instance or if some events appeared in blocks such as blotches, local motion or outliers. In this case, the parameters are rejected and replaced by an interpolation of more reliable neighbor blocks parameters. Once parameters a and b have been found for all blocks, to avoid block effect, parameters are interpolated in each pixel by a bilinear interpolation.

4.1.3

Ohuchi’s method [12]

Local motion and blotches cause problems to the estimator proposed by Roosmalen. Robust estimators are a particular class of estimators suitable to solve these kind problems. From the work of Roosmalen, Ohuchi proposed a method based on robust estimation. He established a hierarchy of flicker-correction model: 1. In+1 (x, y) = a · In (x, y) + b 2. In+1 (x, y) = a(x, y) · In (x, y) + b 3. In+1 (x, y) = a · In (x, y) + b(x, y) 4. In+1 (x, y) = a(x, y) · In (x, y) + b(x, y) with a(x, y) and b(x, y) second order polynomials. 12

The parameters are solved hierarchically. The zeroth order problem is solved with the estimation of Roosmalen[15]. Higher order models are solved by using M-estimators (4.2.2) and by introducing the solutions of lower order model as first estimations. According to him, this method gives better results than Roosmalen’s one.

4.2 Robust Estimation Following Ohuchi, we decided to work with robust cost functions. The flicker estimation presented here relies on the following image sequence model: In+1 (x) = a(x, θ) · In (x) + b(x, θ) + ε(x)

(4.2)

Where In (x) and In+1 (x) are the gray levels of the pixel at the location given by position vector x in frames In and In+1 . The vector θ contains the polynomial coefficients of the gain function a(x, θ) and the offset function b(x, θ).

4.2.1

Matrix formulation

The functions of gain a(x) and offset b(x) are low frequency spatially. They can be expressed in a 2-dimensional polynomial basis. Let Φ be the matrix of the functional basis. 

 φ1 (x1 ) φ2 (x1 ) . . .   Φ =  φ1 (x2 ) φ2 (x2 ) . . .  .. .. .. . . . with φi (x) = φi (x, y) = xk y l for instance. The functions a(x, θ) and b(x, θ) can now be written as follows :  a(x, θ) = Φa b(x, θ) = Φb where a and b are the vectors of parameters: θ = (a, b) This formulation gives for each x: (In+1 )x = (Φa)x · (In )x + (Φb)x + ε(x) The expression (Φa)x · (In )x can be expressed for the whole image by diag (In ) · Φa where In is the column vector that contains all the intensity levels of the picture and diag (In ) is the diagonal matrix whose diagonal elements are the elements of In . e = In+1 − diag (In ) · Φa − Φb 13

(4.3)



    e1 In+1 (x1 ) φ1 (x1 ) In (x1 )  ..     .. .. −  .  =  . . eM

In+1 (xM )  φ1 (x1 )  .. − .

...

...

 φp (x1 ) In (x1 )  ..  .

φ1 (xM ) In (xM ) . . . φp (xM ) In (xM )   φp (x1 ) b1   ..  ..  .  .

φ1 (xM ) . . . φp (xM )

 a1 ..  .  ap

bp

with M the number of pixels in In .

4.2.2 Energy of the residuals The problem is now to find the parameters which minimize some robust cost function E(ε) of the error made on the model. The minimization of this energy was done here through M-estimation [9]. The M-estimation has the advantages of being a robust estimation method, relatively easy to implement and fast. With Least Mean Square formulation, the energy minimized is quadratic in the error. Here a function ρ is introduced to reduce the effect of outliers by replacing the squared residuals by another function of the residuals. The idea is that the large errors are assumed to be most likely due to local motion or blotches and not to flicker. From an implementation point of view, these regions will be weighted out by some weight mask.

E(ε) =

X

ρ(ε(x))

(4.4)

ρ(In+1 − diag (In ) · Φa + Φb)

(4.5)

x

=

X x



 X ˆ ˆ, b = arg min a ρ(In+1 − diag (In ) · Φa + Φb) (a,b)

(4.6)

x

Several functions ρ can be used ([9, 7, 3]). We arbitrarily decided to choose the Cauchy’s one: ρ(ε(x)) =

4.2.3


1 log 1 + ε(x)2 2

Iterative Reweighted Least Square (IRLS)

The minimization of the energy can be done through the classical IRLS algorithm [9, 13]. We need at first the expression of the Jacobian of the residual ε : h i J = diag (In ) · Φ , Φ (4.7) 14

 

  φp (x1 ) In (x1 ) φ1 (x1 )   .. ..  , . . φ1 (xM ) In (xM ) . . . φp (xM ) In (xM ) φ1 (xM )

φ1 (x1 ) In (x1 )   .. J=  .

...

...

  φp (x1 )   ..   .

. . . φp (xM ) (4.8)

The IRLS algorithm.  begin  initialization    do  " # 
−1 T i ai+1  update: = JT W i J J W In+1  i+1  b   update: Wi+1 = ρ(ε)  ε   until convergence: ai+1 − ai ≈ 0, bi+1 − bi ≈ 0 end Where Wi is the diagonal matrix collecting the weights over all the positions x. Wi is estimated at each step given ai and bi .

4.2.4

Scale parameter

The robust estimator depends generally on a scale parameter σρ that controls the value after which errors are considered as outliers. X  ε(x)  E(ε) = ρ σρ x

(4.9)

The IRLS algorithm is working the same way. Only the weights update is modified. The scale parameter is changed at each step of the estimation. We decided to start at σρ = 30 since the flicker can occur intensity shifts of more than 20 levels of gray intensity (over 256 gray levels). And the estimation loop finishes with σρ = 3, which corresponds to the noise level in our case. The scale parameter is divided by 1.5 at each step.

4.3 Results: evaluation of the estimator robustness The robustness of the flicker estimator has been tested on a sequence of pictures of Lenna, on which artificial flicker has been added as well as gaussian noise (N (0, 3)) and uniform distribution of black pixels. The black pixels aim to model outliers (especially blotches). The 15

figure 4.2 shows the absolute difference between original and estimated parameters values in function of the percentage of black pixels in the frame. The estimator gives correct results (a difference of 0.04 in the parameters is not visually perceptible) up to 60% of black pixels. After 60% of black pixels, results are wrong since the estimator simply gives the black picture as a solution.

Figure 4.2: Pictures from the Lenna sequence with artificial flicker, gaussian noise (N (0, 3)) and uniform distribution of black pixels.

Figure 4.3: Absolute difference between estimated and true parameters vs. percentage of black pixels in the image. (because of time, only one estimation has been done for each level of noise)

4.4

Cosine Basis vs. Polynomial Basis

It is interesting to notice that we can choose the nature of the functional basis. It can be a polynomial one or it can be for instance a cosine basis. 16

One advantage of the cosine basis is its orthogonality. We have ΦT Φ = id. And if (but only if!) Φ is squared, ΦΦT = id. This introduces some simplifications in the IRLS algorithm and thus reduces the computational complexity. Optimizations can be even better if we use the LEGEND algorithm [2] instead of the classical IRLS. This algorithm was also described in [9] using iterative modified residuals [3]. With this algorithm we need only one matrix inversion: (ΦΦT )−1 = id. But unfortunately, the computation complexity was not our first priority, we had not enough time to fully compare both algorithms. But this is not enough, we need to compare also the qualities of the results. This can be done by comparing the values of the error energy (section 4.2.2). We have done this experience on a sequence from Rory O’More of 100 frames. Both estimations had been carried out for a fixed number of iterations and at different orders of the basis (up to order 12). The figure 4.4 shows the results of this experience. The blue line corresponds to the mean of the ratio between energies using cosine and polynomials. We can see that using a cosine basis gives lower error energy than using a polynomial basis. That does not mean that the underlying model of flicker is rather cosine than polynomial. But this experience shows at least that cosine basis give comparable results to polynomial basis one. Results are even visually better with cosine basis. So considering the computational advantages we decided to work with a cosine basis. We had not enough time to carry out any studies about the choice of the order of the basis. Usually order 4 or 6 is enough but the movie Rory O’More appeared to need higher order (at least 12 on some shots). The problem is of course to find automatically the needed order.

Figure 4.4: mean of the ratio of the error energy (E(cosine)/E(poly)) vs. order of the basis.

17

4.5

Energy and Symmetry

Consider the problem at zeroth order in x and y, we have basically to do a linear regression between random variables In and In+1 . We could naively expect that the two following symmetric problems (at zeroth order in x and y) will give symmetric solutions:  ε (x) 1 ε2 (x)

= In (x) − a1 · In+1 (x) − b1

(4.10)

= In+1 (x) − a2 · In (x) − b2

Unfortunately this expectation is not correct. As Galton [14, 6] showed when establishing the principle of regression, if the correlation between the random variables In and In+1 is not 1, then we cannot have:  a ˆ2 ˆb2

=

1 a ˆ1 ˆ

= − aˆb11

(4.11)

For more informations about this problem (sometimes called the regression ‘fallacy’ or regression ‘trap’) see the Appendix A. What does this mean for the estimation of the flicker? This means that the estimation will be biased. Especially the gain a will be usually less than 1. To make the problem symmetric, the model can be changed into:  ε (x) = (I (x) − I (x)) − u · 1 (I (x) + I (x)) − b 1 n+1 n 1 2 n+1 n 1 ε2 (x) = (In (x) − In+1 (x)) − u2 · 1 (In+1 (x) + In (x)) − b2 2

(4.12)

with u1 = a1 − 1, u2 = a2 − 1. With such a model and assuming that u ≈ 0 the symmetric estimations give effectively symmetric results. See the appendix A for more informations. We processed on a sequence of 100 frames the both estimations. We found with the classical √ model that the correlation coefficient [16] is not 1 but rather a ˆ1 · a ˆ2 = 0.99. The bias observed is quite important since the gain of each estimation is under evaluated of about 0.01, which means that usually the frames will look darker of about 0.01 ∗ 100 = 1 gray level. √ By adding the symmetry aspect, we have a ˆ1 · a ˆ2 = 0.99999 (for 100 frames), which is visually not perceptible.

18

4.6

Regularization

This part had been implemented but not yet tested in a full restoration process. This section is therefore presented for information. Why do we need a Regularization? If the image is itself zeroth order then there are then too many solutions for a and b: the problem is ill-posed and need therefore a regularization. The regularization can be implemented by adding a constraint on the parameter a(x). The parameter a is usually close to 1. We decided accordingly to base the constraint on this fact. With a bayesian approach, the M-estimation is actually equivalent to the Maximum Likelihood estimation [3]. ˆ = arg max P (In , In+1 |a, b) (ˆ a, b) (a,b) (   ) Y 1 ε(x) = arg max exp − ρ (a,b) 2 σρ x ( ) X  ε(x)  = arg min ρ (a,b) σρ x If we look at the Maximum A Posteriori Estimation: ˆ = arg max P (In , In+1 , a, b) (ˆ a, b) (a,b)

= arg max P (In , In+1 |a, b) · P (a, b) (a,b)

We can now introduce the constraint on a(x) in P (a, b). # " ((Φ · a) (x) − 1)2 1 exp − P (a, b) = P (Φa) = p 2σa2 2πσa2 σa2 is the variance associated to Φa. The regularized energy is: X  ε(x)  1 X E(ε) = ρ + 2 ((Φ · a) (x) − 1)2 σ σ ρ a x x For the IRLS implementation, we then have to solve: " #! " # " # T T Φ Φ 0 a Φ Φa 1 1 0 JT W i J + 2 = JT Wi In + 2 σa σa 0 0 b 0 19

(4.13)

(4.14)

4.7

A new Estimator

4.7.1 A new model for the flicker The previous method (IRLS+DCT basis at order 2) has been successfully tested on a test sequence (see http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/shots/tunnel.avi. But the results appeared to be disappointing for the Rory O’More movie. The main impression of brightness flicker were removed, but some weird fluctuations remained. These fluctuations seemed to appear because the real underlying model for the flicker is not an affine model as it is presented in our previous model: In+1 (x) = a · In (x) + b but a more complex model: In+1 (x) = f (In (x))

(4.15)

where the mapping function f is not affine. To take into account this aspect, we propose a new model for the flicker: In+1 (x) = f (In (x)) · a(x) The global non-linearity is expressed by the function f and the localization aspect by the gain a. The offset function b does not seem here to be useful. We noticed that results are good enough without using b and that parameter b are usually too small to be relevant. And, to take into account the problem of symmetry the model becomes: 1 (In+1 (x) − f (In (x))) = (In+1 (x) + f (In (x))) · u(x) 2

4.7.2

The restoration process

Since we had not enough time to find a way to estimate at the same time f and a, we decided to treat these aspects separately. The model at order zero in x (equation 4.15) can be solved with the Naranjo’s estimator (section 4.1.1). The idea is here to remove the non-linearity by preprocessing, as a first step, a flicker compensation based on Naranjo’s estimator. (see the chapter 5 for more details about the compensation procedure itself). After this first step, another restoration process is needed to deal with the localization aspect of the flicker. The second process is only based on the estimation of the gain a(x). This is done as explained in previous section 4.2 but since the offset function b(x) appear to be not anymore useful, we gain a lot in term of computational complexity:
1 p
Ipn+1 (x) − Ipn (x) = In (x) + Ipn+1 (x) · u(x) 2 20

Where u(x) = Φd and Ipn (x) are the frames obtained after the first step. But we know that Naranjo’s estimator is not robust to outliers such as blotches. So that after the preprossing global flicker still remain (see http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/ shots/shot_08.naranjo.avi for an example of Naranjo’s estimator based restoration of the Rory O’More sequence). Fortunately the error done on the matching of the histograms at this first step can be compensated by the gain factor a(x) so that global flicker is corrected and local flicker is compensated as well.

4.8

Results of the new Estimator

The figure 4.5 show are some results for the Rory O’More sequence. Especially we can see on this sequence how even strong localization of the flicker can be well estimated and then removed (restoration at order 12). Once again, see http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/shots/ to see more results.

Figure 4.5: Frames from the Irish movie Rory O’More (1911). On the top degraded original frames. On the bottom the restored frames. The black diagonal band has been correctly removed

21

Chapter 5 Flicker Compensation We know how to estimate flicker parameters between frames. We need now to compensate the flicker in each frame. The existing methods presented in section 5.1 have the main disadvantages to be partly manual and not able to compensate completely the flicker degradations. We propose therefore in section 5.2 a compensation algorithm capable of removing flicker and preserving the wanted effects, like fading or shadows.

5.1 Overview of the existing methods 5.1.1

A Simple Method

To restore a shot, a simple idea is take one frame as reference frame and then to compensate all the other frames in reference to this frame. But choosing the frame is not done automatically and usually we have to take not one but several frames as reference frame, especially if the shot is long. Moreover, this method is unable to preserve wanted effects like fading or camera flash: all frames will look like the reference frame.

5.1.2

Roosmalen’s Restoration

Roosmalen proposed another method [15]. A reference frame has still to be chosen as reference frame. But here all the other frames are restored recursively: the second frame to the reference frame, the third frame to the restored second frame . . . . But due to some bias in Roosmalen’s estimators, the sequence are darker in time. To correct this fading, Roosmalen introduced a bias in his restoration. The compensated frame IC n is found by estimating the flicker between the current frame In and the previous restored frame IR n−1 . But R the restored current frame In is here a mixture of the compensated frame and the original frame

22

In . C IR n = k · In + (1 − k) · In

The tuning constant k is generally set to k = 0.85.

5.2

Smoothing the integrated parameters

Employing Roosmalen’s process on Rory O’More shows disappointing results: flicker degradations still remained. In particular, the intensity mean of the frames was sightly fluctuating over a period of 20 frames. To improve the results, we developped another restoration process that aims to suppress the flicker noise and to preserve the effects wished by the movie director (fading, camera flashes). This is a general process that we used also for the deshaking (see section 6.3).

5.2.1

The problem

Using the previous estimators we can estimate the coefficients involved in flicker. These estimations give the functions (gn,m )n,m that represent the estimations of the flicker/shake between frame n and m. The functions gn,m can be obtained either directly by estimating the flicker between frames n and m, either by using the estimations between consecutive frames (gn,m (Im ) = gn,n+1 ◦ . . . ◦ gm−1,m (Im ) = gn,n+1 (. . . gm−1,m (Im ))). gn,m (Im ) = In Let n be the current frame. Define this as a reference frame: n is fixed. The figure 5.1 shows the ‘position’ of the neighbor frames according to the reference frame (m 7−→ θgn,m ). The noisy part of the plot is due to the flicker/shake noise and the global shape of the curve corresponds to the effects that we want to preserve.

5.2.2 Formulation of the problem. Our formulation of the problem is partly based on the restoration process that Decencière proposed to remove translational shake [4]. But with our modeling it is possible to deal with much more general cases and not only with transitional shake but also with the two steps of the flicker restoration or with our affine shake removal. We know the degraded frames (In )n and we want
to find the deflickered frames IR n n . The flicker/shake compensation, which is the inverse of the degradation) can modeled by some functions αn such as for each n we have: αn (In ) = IR n 23

Figure 5.1: Example of parameters θgn,m (m) Figure 5.2: Example of parameters θgn,m (m) and θalpha (n).

and smooth corrected parameters.

Actually the estimation is not only measuring the flicker/shake parameters but other effects of same nature. This can be a pan or a zoom effect in case of shake estimation, it can be a fading, a camera flash, a shadow,. . . in case of flicker estimation. We want to preserve these effects in the restored sequence. The functions (µn,m )n,m represent these effects: R µn,m (IR m ) = In

We don’t have access directly to this information. The only available data are the functions (gn,m )n,m that represent the estimations of the flicker/shake between frame n and m. The problem is now to find the compensation functions αn and to apply them on the degraded frames. µn,n+1

IR ←−−−− IR ... n+1 xn x αn αn+1   gn,n+1

In ←−−−− In+1

...

µm−1,m

IR ←−−−− xm−1 αm−1 

IR xm αm 

gm−1,m

Im−1 ←−−−− Im

Figure 5.3: Scheme of the restoration process (inspired from [4]). From the diagram 5.3 yields: gn,n+1 ◦ · · · ◦ gm−1,m = αn−1 ◦ µn,n+1 ◦ · · · ◦ µm−1,m ◦ αm

(5.1)

gn,m = αn−1 ◦ µn,m ◦ αm

(5.2)

θgn,m (m) = θα−1 (m) + εn,m (m) n ◦µn,m

(5.3)

Consider the parameter space:

24

Where θgn,m (m) corresponds to the parameters of the function gn,m (m), θα−1 to the parameters n −1 of αn and εn,m (m) is the perturbations on parameters due to αm .

5.2.3 Removing the noise To suppress the noise εn,m (m), we decided to apply a simple hamming filter on the observed parameters. The noise which is most likely high frequency will be removed by the hamming filter. h i   (5.4) G ∗ θgn,m (m) (k) = G ∗ θα−1 (k) + [G ∗ εn,m (m)] (k) n ◦µn,m (m) {z } | =0

An hypothesis can also be made on the signal signal θα−1 (m). We assume here that the n ◦µn,m filter preserve the signal. This can be expressed by: h i (n) = θα−1 = θα−1 (5.5) G ∗ θα−1 n n ◦µn,m (m) n ◦µn,n (n) Considering that the noise is removed and the signal preserved, we have now:   G ∗ θgn,m (m) (n) = θα−1 n

(5.6)

That means we only need to smooth the curve of the integrated parameters (θgn,m (m)) and to take the smoothed value in n to have the parameters of the degradation (θα−1 ). Knowing these n parameters we can apply the compensation functions αn on all the frames of the sequence. On the figure 5.1 it is possible to see the integrated parameters θgn,m (m) and the corresponding smoothed parameters.

5.2.4

Nature of the noise: choice of the filter

The original curve is noisy, but is it really a high frequency noise? If we have a look at figure 5.1. We can see that it is rather an impulsive noise. For the shake case for instance, some big shakes can occur. Then applying a hamming filter will modify all the other frames (see 5.4). More robust filters can then be used. We decided to use at first a median filter (to be insensible against big shake) and then to apply a hamming filtering for the smoothness. Figure 5.5 shows the results obtained with the combination of median and hamming filter. Usually we used a size for the filter window of 3 frames for the median pass of 30 frames (more than 1 second of movie) for the lowpass filter. Depending on the nature of the degradation other solutions are possible.

25

Figure 5.4: Example of integrated param- Figure 5.5: Example of integrated parameeters smoothed by a hamming filter.

ters smoothed by a median/hamming filter.

5.2.5 Conclusion The method proposed above was successfully used to restore the Rory O’More movie, as well for the deshaken process as for the flicker removal process. Pans were protected and the steadiness restored. But should a movie be completely restored? Some people asked us if our restoration was not too smooth. An old movie should be sightly shaking! We can not answer this question but the method proposed allows the user to define the kind of restoration that he wishes. By changing the nature of the filter he can keep some slightly shaking or flicker.

26

Chapter 6 IFA Movie 6.1

Irish Film Archive

The Irish Film Archive (IFA) acquires, preserves and makes permanently accessible Ireland’s film and audiovisual heritage. Its holdings include an extensive collection of film and magnetic tape, film stills and posters, a paper archive, film equipment and memorabilia. Material from the Archive collection is screened on a monthly basis in the Irish Film Centre cinemas. To establish a link between the Irish Film Archive and the SIGMEDIA lab, Dr Anil Kokaram proposed to restore an old Irish movie from the beginning of the 20th century. The movie, named Rory O’More and produced by Sydney Olcott in 1911, needed a strong restoration procedure. Different recent automatic restoration processes were employed. A public presentation of the results was performed by Dr Anil Kokaram the 4th of September 2002 in the Irish Film Centre of Dublin.

6.2

The movie

Rory O’More was an Irish hero, leader of the rising of 1641 in Ireland against King of England Charles I. The movie relates the events of his capture and his escape (see Appendix). The movie contains some 12291 frames in poor state of preservation. The figure 6.1 shows some of these frames. The frames are 720x576 pixel each. Severe degradations appear in the movie and especially shake, flicker, blotches, lines and grain noise (these degradations are described below).

27

Figure 6.1: Frames from the Irish movie Rory O’More (1911).

6.3

The restoration process

Several processes have been implemented and used to carry out the full restoration of the movie. Digitalization. The first step was the digitalization of the movie. The physical movie was scanned by telecine machine, put onto BestaSP and then captured to Abekas files. Detection of the shotcuts. For easier manipulations, we decided to work on the shots of the movie. There is 38 shots of a few hundred mega bytes each. Motion estimation/removal. To the viewer, the presence of shake is the most disconcerting example of movie degradation. Instead of steady camera motion, the picture appears to shake. There can be many reasons for this phenomenon. One is that the scanner or the digitizing device (téléciné) is not properly matched to the position of the film: scanned frames are displaced. Another is that the camera was held by hand. Similar to the flicker, we need at first to estimate the parameters of the shake before restoring the frames. The estimation of the shake concerns the estimation of the global motion of the camera. What is tracked is the background and not the local motions of the objects. An affine model was chosen to describe the displacement of the frames (translation, zoom and rotation). The affine model is needed since some frames are rotated. The estimation itself is performed 28

through M-estimators. Part of the algorithm is described in [5, 1]. We used then the restoration process described in chapter 5 to compensate the frames. The estimation process takes around 3 seconds per frame on a PC and about 10 hours for the whole movie. The compensation itself is very fast. Flicker estimation/removal.

We used the process described in section 4.7.

The restoration process takes around 40 seconds per frame on a computer and about 10 hours on 10 PC’s for the whole movie. But programs are here running on matlab. More clever programming should definitely reduces the computation time. Blotch removal/Noise reduction [10]. Blotches are one of the most common degradations appearing in archived footage. Large pieces of dirt, holes in the film emulsion, and dust on the film are all perceived as blotches. Because blotches are only present in one frame (i.e. for 1/25th of a second), they appear and disappear very rapidly, and so are particularly distracting. Noise in film refers to small, random fluctuations in individual pixels. It looks like a grainy kind of flicker. Noise is quite common in old footage, and has a strong effect on the perceived quality of the film. Furthermore, the presence of noise in the film makes blotch removal harder. However, noise reduction is very difficult in blotchy film! Our deblotching process therefore treats noise and blotches simultaneously, using a Markov Chain Monte Carlo (MCMC) technique. This is the JOMBADI algorithm described in [10]. The restoration process takes around 40 seconds per frame on a computer and about 10 hours on 10 PC’s for the whole movie. But some frames where too badly damaged and needed more attention. Line Removal.

The algorithm used to remove line scratches is described in [10].

Image enhancement.

An unsharp masking has been applied to enhancement the pictures.

Crop and zoom. The deshaken process causes a problem: since pictures are displaced, parts of the compensated image are not defined and so black. The deblotch process can inpaint some of these regions but not all. So we crop the visible part of the movie and we zoom it to have the same size. Why this order? it is interesting to see which defect should be removed at first. Each process will of course works better if all other problems have been fixed before. For instance flicker es29

timation works better if differences between images are only due flicker. So how did we choose this order? Problems can actually be separated in global problems (flicker and shake) and local problems (noise reduction, lines, blotches, . . . ). The restoration of global problems is typically more robust than for the local problems. Flicker and shake should be therefore fixed prior to the other problems. Deshaking process was chosen as first process since shake was very strong in the movie. Actually the chosen order corresponds also the human perception. The viewer is at first disturbed by the shake, then, when the shake is fixed, by he will noticed the flicker and eventually the blotches and the grain noise.

6.4

Some results

It is hard to show on a paper results of a movie restoration because the temporal aspect is very important. So you should go to http://www.mee.tcd.ie/∼sigmedia/ifc/flicker/shots/ifc/ too see video results. The figure 6.2 shows some of the original frames and restored frames.

Figure 6.2: Frames from the Irish movie Rory O’More (1911). On the top the degraded frames whom blotches has been circled in red. On the bottom the restored frames

30

Chapter 7 Conclusion and future work. The goals of our work were to establish a diagnosis module for flicker detection, to improve the existing method of flicker estimation, and to eventually to provide a full flicker removal algorithm. We proposed some ideas and built a prototype for the flicker diagnosis module. This part needs some further testing and certainly some improvement to increase the robustness. We implemented a robust estimator of flicker severity based on Ohuchi’s method [12]. However, the classical restoration procedures failed [12, 11] to restore correctly the heavily degraded movie Rory O’More. Therefore, we improved several important aspects of the classical estimators, such as removing the bias due to regression and changing the functional basis. Eventually we proposed a new model of flicker. This new model now takes both the non-linearity of the histogram distortion and the localization of the flicker into account. From this model we implemented a new flicker-estimator and the corresponding restoration procedure. Concerning the restoration procedure, the method we proposed can now remove the flicker degradation while preserving desired effects such as fading. The new estimators and restoration procedures have been successfully used to restore the 12500 frames of Rory O’More. What are the possible improvements? Since our algorithms succeeded in restoring Rory O’More, which was degraded to an exceptional extent, we can expect that this restoration process should work for almost all movies. But some improvements are still possible, especially in term of the computational complexity of the algorithm. The quality may be also improved through the use of a different basis or by using a more robust histogram matching technique. In the longer term, it will be interesting to consider the interactions between the different modules involved in the digital movie restoration. Flicker and shake are for instance interdependant: a big shake engenders a wrong flicker estimation and flicker degradation engenders a wrong global motion estimation. Another research direction could be to consider the digital restoration in a feature based approach. For the flicker case, the residual energies we used

31

are pixel based, but we could imagine other ways to express these energies by looking other basic features. Finally, it is remarkable to see how important it was to test our algorithms on 12500 heavily degraded frames and not only on test sequences of 25 frames. Working on a such amount of frames encouraged us to improve our algorithms, sometimes even to consider new models. It was also a great experience to apprehend how all these algorithms should be linked.

32

Bibliography [1] P. Bouthémy, M. Gelgon, and F. Ganansia. A unified approach to shot change detection and camera motion characterization. IEEE Transactions on Circuits and Systems for Video Technology, 9:1030–1044, 1999. [2] P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud. Two deterministic half quadratic regularization algorithms for computed imaging. In International Conference on Image Processing, pages 168–172, Austin, USA, 1994. [3] R. Dahyot. Appearance based road scene video analysis for the management of the road network. PhD thesis, University of Strasbourg I, France, November 2001. [4] E. Decencière. Restauration automatique de films anciens. PhD thesis, Ecole Nationale Supérieure des Mines de Paris, December 1997. [5] P. Delacourt, A. Kokaram, and R. Dahyot. Comparison of global motion estimators. In proceedings of Irish Signals and Systems Conference, Cork, Ireland, June 2002. [6] F. Galton. Regression towards mediocrity in hereditary stature. Journal of the Anthropological Institute, pages 246–263, 1886. [7] S. Geman and G. Reynolds. Constrained restoration and the recovery of discontinuities. IEEE Transactions on Pattern Analysis Machine Intelligence, 14(3):367–383, March 1992. [8] R. Gonzalez and P. Wintz. Digital Image Processing. 2nd edition, 1987. [9] P.J. Huber. Robust Statistics. John Wiley and Sons, New York, 1981. [10] A. Kokaram. Motion Picture Restoration: Digital Algorithms for Artefact Suppression in Degraded Motion Picture Film and Video. Springer Verlag, 1998. [11] V. Naranjo and A. Albiol. Flicker reduction in old films. In Proc. of the 2000 International Conference on Image Processing (ICIP-2000), September 2000. 33

[12] T. Ohuchi, T. Seto, T. Komatsu, and T. Saito. A robust method of image flicker correction for heavily-corrupted old film sequences. In Proc. of the 2000 International Conference on Image Processing (ICIP-2000), September 2000. [13] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C - The Art of Scientific Computing. Cambridge University Press, 1995. [14] S.M Stigler. The History of Statistics. Cambridge: Belknap Press of Harvard Press, 1986. [15] P.M.B. van Roosmalen. Restoration of archived film and video. PhD thesis, Delft University of Technology, October 1999. [16] E. Weisstein. The world of mathematics. Available at http://mathworld.wolfram.com/. [17] X. Yang and N. Chong. Enhanced approach to film flicker removal. Society of PhotoInstrumentation Engineers, 2000.

34

Appendix A Energy & Symmetry A.1

The regression problem

This is a mathematical expression of the regression problem exposed by Stigler [14]. Let X have a normal distribution N (0, c2 ), and let Y = X + Z, where Z corresponds to the noise and is N (0, b2 ), independent of X. Then Y is N (0, b2 + c2 ) and clearly the conditional expectation of Y given X = x is found to be E(Y |X = x) = x. But what is the conditional expectation of X given Y = y? After some calculations[14] we find that E(X|Y = y) is not y 2 but rather it is E(X|Y = y) = b2c+c2 y.

A.2

Looking at the residuals

We expose here another formulation which is compatible with the M-estimator theory. Assume that the gain a is closed to 1: a = 1 + u ≈ 1. Instead of solving the following problems:  0 (x, a , b ) = a · I (x) + b − I (x) 1 1 1 n 1 n+1 1 0 (x, a2 , b2 ) = a2 · In+1 (x) + b2 − In (x)

(A.1)

2

We propose by mixing the previous equations to solve:   (x, u , b ) = (I (x) − I (x)) − u · 1 (I (x) + I (x)) − b 1 1 1 n+1 n 1 2 n+1 n 1 2 (x, u2 , b2 ) = (In (x) − In+1 (x)) − u2 · 1 (In+1 (x) + In (x)) − b2 2

35

(A.2)

The hypothesis done here are only a = 1 + u ≈ 1. We want to show that the new problems are symmetric. The symmetry is expressed in the estimated parameters by:  a ˆ2 = 1 + uˆ2 = aˆ11 ≈ 1 − uˆ1 (A.3) ˆ ˆb2 = − b1 ≈ −ˆb1 a ˆ1

Both problems will be equivalent if we can express one problem from the other one just by changing the name of the variables. Let define u2 = −u1 and b2 = −b1 . Consider the expression of 2 : 2 (xi , u2 , b2 )2 = 2 (xi , −u1 , −b1 )2 = ((In − In+1 ) + u1 · (In + In+1 ) + b1 ))2 = ((In+1 − In ) − u1 · (In + In+1 ) − b1 ))2 = 1 (xi , u1 , b1 )2

The residuals are equivalent. So both problems are equivalent: (ˆ u2 , ˆb2 ) = arg min E2 (u2 , b2 ) u2 ,b2 X = arg min 2 (x, u2 , b2 )2 u2 ,b2

x

X

= arg min

−u1 ,−b1

2 (x, −u1 , −b1 )2

x

= arg min E1 (u1 , b1 ) −u1 ,−b1

= (−ˆ u1 , −ˆb1 ) We have (ˆ u2 , ˆb2 ) = (−ˆ u1 , −ˆb1 ) and we know that we should have (ˆ u2 , ˆb2 ) ≈ (−ˆ u1 , −ˆb1 ). Depending on the validity of the assumption u ≈ 0 (cf. experiment developed in section 4.5) our model is correct. Note that if we don’t change the modeling of the problems, then the expressions of the residual are not equivalent: 2 (xi , a2 , b2 )2 =

1 1 (xi , a1 , b1 )2 a21

36

And so: (ˆ a2 , ˆb2 ) = arg min E2 (a2 , b2 ) a2 ,b2 X = arg min 2 (x, a2 , b2 )2 a2 ,b2

= arg

x

X 1  (x, 1/a1 , −b1 /a1 )2 2 2 1/a1 ,−b1 /a1 a 1 x min

= arg

min 1/a1 ,−b1 /a1

6= (

1 E1 (a1 , b1 ) a21

ˆb1 1 ,− ) a ˆ1 a ˆ1

The problems are here not symmetric.

37

Appendix B Rory O’More: the rising of 1641 in Ireland Rising of 1641 in Ireland: Charles I’s policies and confiscations cause insurrection in Ulster and this rebellion triggered the English Civil War, which put an end to King Charles I’s attempt to create an absolutist state (represented in Ireland by the policies of his lord deputy, Thomas Wentworth, 1st earl of Strafford). Great Catholic-Gaelic rebellion for return of lands, later joined by Old English Catholics in Ireland. Under leadership of Irish chieftain, Rory O’More, conspiracy was formed to seize Dublin and expel the English. English settlers were driven out of Ulster. Catholics hold 59% of land in Ireland.. The native Irish assisted by the Hiberno-English rose in 1641 rebellion and attacked the settlements of Anglicans and Presbyterians. They drowned, murdered, and burned alive, men, women and children. While the stories say that they killed 200,000 people, this is probably grossly overestimated.( from http://www.celticcrossroads.com/celt1600.htm)

38