International audienceWe present a new blind deconvolution method for video sequence. It is derived following an inverse problem approach in a Bayesian framework. This method exploits the temporal continuity of both object and PSF. Combined with edge-preserving spatial regularization, a temporal regularization constrains the blind deconvolution problem, improving its effectiveness and its robustness. We demonstrate these improvements by processing various real video sequences obtained by different imaging techniques

Dauphin, Raphaël

Gressard, Alain

Soulez, Ferréol

Thiébaut, Éric

Tourneur, Yves

HAL-INSU

BLIND DECONVOLUTION OF VIDEO SEQUENCESFerre´ol Soulez, E´ric Thie´baut, Yves Tourneur, Alain Gressard, Raphae¨lDauphinTo cite this version:Ferre´ol Soulez, E´ric Thie´baut, Yves Tourneur, Alain Gressard, Raphae¨l Dauphin. BLINDDECONVOLUTION OF VIDEO SEQUENCES. 15th International Confe´rence on Image Pro-cessing, Oct 2008, San Diego, United States. pp.0, 2008. <ujm-00293656>HAL Id: ujm-00293656https://hal-ujm.archives-ouvertes.fr/ujm-00293656Submitted on 7 Jul 2008HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.BLIND DECONVOLUTION OF VIDEO SEQUENCESFerre´ol Soulez1,2, E´ric Thie´baut1, Yves Tourneur3, Alain Gressard4, Raphae¨l Dauphin41Centre de RechercheAstrophysique de Lyon,CNRS-UMR 5574, Universite´Lyon 1, Ecole NormaleSupe´rieure de Lyon, Observatoirede Lyon, 9 avenue CharlesAndre´, Saint-Genis Laval cedex,F-69561, France.2 Laboratoire Hubert Curien(ex-LTSI), CNRS-UMR 5516,Universite´ Jean Monnet, 18 ruePr Benoıˆt Lauras, F-42000Saint-Etienne, France.3Centre Commun deQuantimetrie, Universite´ Lyon1,8 avenue Rockefeller - 69373Lyon cedex 08, France4 Service de Cardiologie del’hoˆpital de la Croix Rousse,Hospices Civiles de Lyon, FranceABSTRACTWe present a new blind deconvolution method for video sequence.It is derived following an inverse problem approach in a Bayesianframework. This method exploits the temporal continuity of both ob-ject and PSF. Combined with edge-preserving spatial regularization,a temporal regularization constrains the blind deconvolution prob-lem, improving its effectiveness and its robustness. We demonstratethese improvements by processing various real video sequences ob-tained by different imaging techniques.Index Terms— blind deconvolution, denoising, image recon-struction, video signal processing.1. INTRODUCTIONThe problem of reconstructing a signal using observation blurred byunknown process is called blind deconvolution. Nowadays, blindimage deconvolution receives increasing attention from the aca-demic world (see [1] for a review). Although, many authors usevideo sequences in a super-resolution framework (a multi-frameblurred observation of the same scene (see e.g. [2]), to our knowl-edge, few of them have studied the specific problem of blind videosequence deconvolution [3]. In this paper, we propose to use thetemporal correlations in video sequence to design an effective blinddeconvolution algorithm.2. MODEL DESCRIPTIONThe video sequence g(s, t) observed in a direction s is composedof T successive images gt(s) taken at different instant t between itsbeginning t1 and its end tT . It is given by:g(s, t) ="h(s, t|s′, t′)x(s′, t′)ds′dt′ + n(s, t), (1)where x(s′, t′) is the object brightness at instant t′, h(s, t|s′, t′) is thepoint spread function (PSF), and n(s, t) account for the noise (sourceand detector). The PSF h(s, t|s′, t′) is the observed brightness distri-bution at instant t in the direction s for a point source located at theinstant t′ in direction s′.For data sampled on N pixels, Eq. (1) can be written in a matrixform:g = H · x + n, (2)where g = (gTt1 , . . . , gTtT)T and x = (xTt1 , . . . , xTtT)T are vectors of sizeN · T , n is a N · T random vector and H is a N · T × N · T matrix.In this study, we make two major assumptions. First, there is notemporal spread (gt = Ht ·xt+nt ,∀t) and H becomes block diagonal.Then, the PSF is shift invariant (isoplanatic) and Ht can be expressedusing spatial convolution product ∗ and its first row ht:gt = ht ∗ xt + nt. (3)Under circulant approximation, this can be evaluated rapidly usingFFTs.3. BAYESIAN APPROACHIn blind deconvolution both the object vector and the PSF h (definedby h = (hT1 , . . . , hTT )T), of size N · T , must be estimated using (N · T )measurements in g.Following [4], the “maximum a posteriori” solution is:{x, h}MAP = argmax{x,h}Pr{x, h|g}. (4)From Bayes’ theorem, assuming x and h are independent:Pr{x, h|g} =Pr{g|x, h}Pr{x}Pr{h}Pr{g}(5)and since Pr{g} does not depend on the model, we can write:{x, h}MAP = argmax{x,h}(Pr{g|x, h}Pr{x}Pr{h}) . (6)Maximizing Pr{x, h|g} is equivalent to minimizingε(x, h) = − log Pr{x, h|g} = Φlkl(x, h; g) + Φobj(x) + Φpsf(h), (7)The penalizing function to minimize defined in Eq. (7) is the sumof three terms : a likelihood penalty Φlkl(x, h; g) ensuring the agree-ment between the model x∗ h and the data g, and two regularizationpenalty Φobj(x) and Φpsf(h) introducing subjective a priori knowl-edge about the object and the PSF respectively.3.1. The likelihood penalty termFor Gaussian noise, the likelihood penalty reads:Φlkl(x, h; g) = [g − m(x, h)]T · C−1noise · [g − m(x, h)], (8)where Cnoise is the covariance matrix of the noise and m(x, h) themodel. Using Eq. (2) the model is defined as:m(x, h) = H · x, (9)and is evaluated using Eq. (3) for all t. For uncorelated noise, Cnoiseis diagonal and Eq. (8) simplifies to:Φlkl(x, h, g) =∑t∑kwk,t((Ht · xt)k − gt,k)2, (10)where 1/wr,t is the noise variance for pixel k of frame t. This modelcan cope with non-stationary noise and can be used to express con-fidence on measurements on each pixel of the data. Thus it can dealwith unmeasured pixels on the sensor (1/σ2 = 0 for such pixels).3.2. a priori on the objectAs the different dimensions of the object are not homogeneous, wetake a separable regularization term:Φobj(x) = µobjΨobj(x) +T∑t=1λobjt Θobj(xt) , (11)where Ψobj(x) is a temporal regularization and Θobj(xt) is a spatialregularization at each instant t.3.2.1. Spatial regularization on the objectThere are many different kinds of spatial regularization used in im-age reconstruction. As the noise mostly contaminates high frequen-cies, smoothness is the most effective regularization constraint toavoid amplification of noise. To avoid oversmoothing of sharp fea-tures caused by quadratic regularization we choose an edge preserv-ing regularization. This is achieved by taking, at each instant t:Θobj(xt) =∑k∑v∈Vkϕ(xt,k − xvdv), (12)where xv is the value of a pixel v in the neighborhood Vk of pixelk (here a V8 neighborhood) and dv its distance to this pixel. In thiswork, we choose a l1–l2 norm. This norm is asymptotically quadratic(resp. linear) for small (resp. large) pixel differences compared to thethreshold η. It is defined by:ϕ(u; η) = 2 η2[|u|/η − log (1 + |u|/η)], (13)The value of the parameter η is not critical for the reconstruction andcan be approximately fixed once for all to the value of one quan-tization level. In that case, this regularisation is close to a pure ℓ1regularization that promote sparsity of gradient[5], but it can be min-imized faster by our optimization algorithm.3.2.2. Temporal regularization on the objectIf the scene evolution is slow compared to the time sampling, asmoothing temporal regularization, similarly to the spatial regular-ization but along the time axis, can be defined.3.3. a priori on the PSFAs for the object, our PSF regularization term is split in:Φpsf(h) = µpsfΨpsf(h) +T∑t=1λpsft Θpsf(ht) . (14)3.3.1. Spatial PSF regularizationAs there are many different cause of blur (defocus, motion, diffrac-tion, diffusion...), there are a lot of totally different shapes of PSF. Inthis context we choose a parametric function p(θ) (Gaussian, gener-alized Gaussian, Lorentzian...) to set the prior shape of the PSF, thenthe regularization becomes [6]:Θpsf(h) = (h − p(θ))TW(h − p(θ)), (15)where W is a weight matrix. In our work to constraint the PSF tobe increasingly close to the prior shape as it moves away from itscenter, we consider a diagonal W with a power law on the diagonal.The functions p(θ), chosen according to the experimental conditionsshall have few parameters (width, orientation).3.3.2. Temporal PSF regularizationAs for the object, if the PSF evolution is slow compared to thetime sampling, a quadratic smoothing temporal regularization canbe used, e.g. :Ψpsf(h) =T−1∑t=2(2 ht − ht−1 − ht+1)2 , (16)4. ALGORITHM SUMMARYIn this Bayesian framework, reconstructing the de-blurred video se-quence amounts to determine the couple {x+, h+} that minimizes thecriterion defined in Eq. (7), which writes:ε =Φlkl(x, h; y) + µobjΨobj(x) + µpsfΨpsf(h)+T∑t=1(λobjt Θobj(xt) + λpsft Θpsf(ht)),(17)The optimal reconstructed image x+ and PSF h+ depend on the valueof each parameter (four hyper-parameters per frame).4.1. Hyper-parameters SettingTo simplify the determination of hyper-parameters, we have madeseveral simplification hypothesis. We suppose the noise statisticalproperty, the image dynamic and the PSF shape identical in eachframe. As a consequence: (i) the object spatial hyper-parametersin each frame are identical ; (ii) the PSF spatial hyper-parametersin each frame are identical. Thus we only have now four hyper-parameters (λobj, λpsf , µobj and µpsf).Despite this simplification, choosing the optimal values of thehyper-parameters is cumbersome and difficult. Whether methodssuch as generalized cross-validation (GCV)[7] or the L-curve [8]are suitable for this task deserves an extensive study which is outof the scope of this paper. In the present work, we simply choosehyper-parameter values by visual inspection of the resulting image.4.2. Minimization MethodAn alternating minimization[9, 10] scheme is used to minimize thecriterion.1. initialize the PSF with its prior shape h(0) = p(θ) ,2. estimate the optimal object x(k+1) given the PSF h(k),3. estimate the optimal PSF h(k+1) given the object x(k+1),4. repeat steps 2 and 3 until convergence or after a defined num-ber (k) of iterations.In order to determine the optimal image x+ and PSF h+ in ourinverse problem approach, one has to minimize a criterion with re-spect to a very large number of variables (all the pixel values forevery frames). To that end, we used the VMLM-B algorithm [11]which is a limited memory variant of the variable metric method withBFGS updates [12]. This algorithm, which can further accounts forbound constraints on the parameters. We make use of these boundconstraints to enforce PSF positivity. This algorithm has proven ef-fectiveness for image reconstruction and only requires the computa-tion of the penalty function to be minimized and its gradient. Thememory requirement is a few times the size of the problem.5. EXPERIMENTAL RESULTSThis algorithm was used with different experimental data sets fromboth medical and biological fields, and with different dynamicalimaging techniques.5.1. Coronarography(a) Raw image (b) Reconstruction of 1(a)Fig. 1. One frame of a coronarography sequenceThe coronarography is a radiological exam to observe coronaryarteries (heart arteries). This gives short (several seconds) video se-quences displaying the motion of an impervious to X-ray product inthe coronary arteries. These sequences are examined by cardiologistto localized obstacle in the arteries, possibly responsible of an heartattack.Ten sequences were acquired with the same system, and werethen processed with the same assumptions. Several item can be no-ticed by a simple observation of the raw data (see Fig. 1(a)). (i) Themovement of the heart is too important between successive framesto use a temporal constraint on the object. So we set µobj = 0. (ii)As the blur seems to be caused by tissue in the axe of projection, andas neither the patient nor the acquisition system was moving duringthe recording, the PSF is supposed to be constant. As consequence,µpsf = ∞ or equivalently ht = h1 ∀t. In that case the number of un-known parameters on the PSF is divided by T , increasing the compu-tation speed. (iii) The measured area is not square and is cropped bysome parts of the apparatus which remains the same in every framesof every sequences. This is taken in account in the a priori weightmap W estimated by a basic thresholding and shown in Fig. 2(b) :wk,t = wk ={1 if k-th pixel is measured,0 otherwise.(18)The prior PSF p(θ) has a Lorentzian shape with only one parame-ters: the full width at half maximum θ. Finally only three hyper-parameters have to be determined to perform this blind deconvolu-tion: λpsf , λobj and θ.(a) Estimated PSF (b) Weight mapFig. 2. Estimated PSF and weight map used for coronarographysequence shown in Fig. 1These sequences of about 75 frames of 512 × 512 pixels wereprocessed in about 80 seconds per frames on a Pentium IV CPU at3.60 GHz. The assumption of an identical PSF in every frames ofa same sequence considerably constrains the PSF and improves theconvergence of the h+ determination stage. A further relaxing ofthis temporal constraint on the PSF, does not show significant im-provements and thus the assumption of a constant PSF is verified. Aframe of one of the studied sequences is shown Fig. 1(a), along withthe coresponding restored frame and the estimated PSF Fig. 2(a). Al-though the working quantization was very small (at most 10 levelsbetween pixels in the arteries and the background), a visual assess-ment shows the effectiveness of our technique. Motion perceptionin the deconvolved video sequences is greatly improved as the sep-aration of the coronaries and the background is enhanced. Let uspoint out that the cropped area does not disturbed the deconvolution,even for pixels close to the border of this area. These unmeasuredparts were just filled with uninformative smooth background in ac-cordance with the smoothing prior. These restored sequences wereshown to cardiologists who confirmed the potential utility of the pre-sented technique. If this method proves efficiency with with noisiersequences, they can consider a decrease of the X-ray beam to de-crease the irradiation of the patient.5.2. Confocal MicroscopyFig. 3. Mitochondria: frames from the raw sequence (right) and theblind deconvolution (left).Confocal microscopy is an optical technique used for imagingwith short depth of field, eliminating out of focus images. In thisexperiment, a fluorescence product was used to mark mitochondriain a “non beating HL-1” cell [13]. Temporal evolution of these mi-tochondria was observed by confocal microscopy. As in the coro-narography case and for the same reasons, the assumptions µobj = 0and µpsf = ∞ are made. Both raw image and its coresponding framein reconstructed sequence are presented in Fig. 3. The raw sequenceis very noisy and the results clearly demonstrates similar ability ofregularized blind deconvolution to not only enhance the resolutionbut also to reduce the noise.5.3. Conventional MicroscopyFig. 4. Epithelial Cell: frames from the raw sequence (right) and theblind deconvolution (left).A sequence of epithelial hair cells taken with conventional trans-mission microscope was processed by our method. It is composed of400 frames of 253×178 square pixels. In this sequences, the depth offield is relatively thick and, as cilia move in the spatial three dimen-sions, successive focalisations and defocalisations of the same ciliumcan be observed. The blur on these cilia can be divided in blur due tothe apparatus, which is isoplanatic and out of focus blur, which is notisoplanatic. As our method cannot cope with none isoplanatic PSFwe tried to remove only the blur due to the apparatus which wereconsidered isoplanatic and identical on every frames (µpsf = ∞). Asonly the hairs are moving rapidly, a segmentation based on temporalvariation of pixels is made. A temporal constraint µobj is set for thepixels of the background and the cell, but it is relaxed for pixels nearcilia (µobj = 0). A frame of both raw and deconvolved sequences areshown in Fig. 4. On several pixels, disturbing diffraction figures dueto dust on a glass of the apparatus can be seen. This phenomenon istaken into account by setting the weight these pixels at w2k= 0 as inSec. 5.1. Even in this difficult case, our method achieves to improveresolution of both hairs and inner structures of the cell.6. CONCLUSIONThis paper presents a new method for blind deconvolution of videosequences. We exploit both spatial and temporal continuity toachieve a good characterization of the PSF and then a good videosequence blind deconvolution. This method had proven its capabilityand robustness with various experimental data.ACKNOWLEDGMENTSThe authors are grateful to those who have developed biologicalpreparations and acquired the video: Sophie Pelloux for the mito-chondria and Brigitte Chhin for the epithelial cells.7. REFERENCES[1] P. Campizi and K. Egiazarian, Eds., blind image deconvolu-tion: theory and application, CRC Press, 2007.[2] T. J. Schulz, “Multiframe blind deconvolution of astronomicalimages,” JOSA-A, vol. 10, no. 5, pp. 1064–1073, may 1993.[3] D. Li, R.M. Mersereau, D.H. Frakes, and M.J.T. Smith, “Newmethod for suppressing optical turbulence in video,” in Proc.European Signal Processing Conference (EUSIPCO), 2005.[4] E. Thie´baut and J.-M. Conan, “Strict a priori constraints formaximum likelihood blind deconvolution,” JOSA-A, vol. 12,no. 3, pp. 485–492, March 1995.[5] David L Donoho and Michael Elad, “Optimally sparse rep-resentation in general (nonorthogonal) dictionaries via l min-imization.,” Proc Natl Acad Sci U S A, vol. 100, no. 5, pp.2197–2202, Mar 2003.[6] A. Tarantola, Inverse Problem Theory and Methods for ModelParameter Estimation, Society for Industrial Mathematics,2005.[7] Gene H. Golub, Michael Heath, and Grace Wahba, “Gener-alized cross-validation as a method for choosing a good ridgeparameter,” Technometrics, vol. 21, pp. 215–223, 1979.[8] P. C. Hansen, “Analysis of discrete ill-posed problems bymeans of the l-curve,” SIAM Rev., vol. 34, no. 4, pp. 561–580,1992.[9] Y.L. You and M. Kaveh, “A regularization approach to jointblur identification and imagerestoration,” Image Processing,IEEE Transactions on, vol. 5, no. 3, pp. 416–428, 1996.[10] Tony F. Chan and C.K. Wong, “Convergence of the alternat-ing minimization algorithm for blind deconvolution,” LinearAlgebra and its Applications, vol. 316, pp. 259–285, apr 2000.[11] E. Thie´baut, “Optimization issues in blind deconvolution algo-rithms,” in Astronomical Data Analysis II., Jean-Luc Starck,Ed., dec 2002, vol. 4847, pp. 174–183.[12] Jorge Nocedal, “Theory of algorithms for unconstrained opti-mization,” Acta Numerica, vol. 1, pp. 199–242, 1992.[13] S. Pelloux, J. Robillard, R. Ferrera, A. Bilbaut, C. Ojeda,V. Saks, M. Ovize, and Y. Tourneur, “Non-beating hl-1 cellsfor confocal microscopy: Application to mitochondrial func-tions during cardiac preconditioning,” Progress in Biophysicsand Molecular Biology, vol. 90, no. 1-3, pp. 270–298, 2006.

BLIND DECONVOLUTION OF VIDEO SEQUENCES

Abstract

Similar works

Full text

Available Versions

HAL-INSU

HAL-UJM

HAL-ENS-LYON

HAL-HCL

Hal-Diderot