Singular Value Decomposition - Applications In Image .

Singular Value Decomposition Applications in Image ProcessingIveta HnětynkováKatedra numerické matematiky, MFF UKÚstav informatiky, AV ČR1

Outline1. Singular value decomposition2. Application 1 - image compression3. Application 2 - image deblurring2

1. Singular value decompositionConsider a (real) matrixA Rn m, r rank (A) min {n, m} .A hasmnrcolumns of length n ,rows of lenght m ,is the maximal number of linearly independentcolumns (rows) of A .3

There exists an SVD decomposition of A in the formA U ΣV T,where U [u1, . . . , un] Rn n, V [v1, . . . , vm] Rm m are orthogonal matrices, and"Σ Σr 00 0 # Rn m, Σr σ1 .σr r r , Rσ1 σ2 . . . σr 0 .4

Singular value decomposition – the matrices:A{ui} i 1,.,n{vi} i 1,.,m{σi} i 1,.,rUSVTare left singular vectors (columns of U ),are right singular vectors (columns of V ),are singular values of A.5

The SVD gives us:span (u1, . . . , ur ) range (A) Rn,span (vr 1, . . . , vm) ker (A) Rm,span (v1, . . . , vr ) range (AT) Rm,span (ur 1, . . . , un) ker (AT) Rn,spectral and Frobenius norm of A, rank of A, .6

Singular value decomposition – the subspaces:ATARnRmTker(A ).uns1s2sr}00}u1range(A) u.2urur 1v1v2Trange(A).vrvr 1. ker(A)vm7

The outer product (dyadic) form:We can rewrite A as a sum of rank-one matrices in the dyadic formA U ΣV T Tv1σ1 . . [u1, . . . , ur ] . σrvrT u1σ1v1T . . . ur σr vrTrX σiuiviTi 1rXAi . i 1 Moreover kAik2 σi gives kA1k2 kA2k2 . . . kAr k2 .8

Matrix A as a sum of rank-one matrices: AA1A2.Ar -1ArrASAii 1SVD reveals the dominating information encoded in a matrix. Thefirst terms are the “most” important.9

Optimal approximation of A with a rank-k:The sum of the first k dyadic termskXAi i 1kXσiuiviTi 1is the best rank-k approximation of the matrix A in the sense ofminimizing the 2-norm of the approximation error, tj.kXuiσiviT argminX Rn m , rank (X) ki 1{kA Xk2}This allows to approximate A with a lower-rank matrixA kXi 1Ai kXσiuiviT .i 110

Different possible distributions of singular values:Singular values of the BCSPWR06 matrix (from Matrix Market).2Singular values of the ZENIOS matrix (from Matrix Market).101010010010 210 410 1010 610 8 201010 1010 3010 1210 1410 4010 1610 1810 50050010001500101002003004005006007008009001000The hardly (left) and the easily (right) approximable matrices(BCSPWR06 and ZENIOS from the Harwell-Boeing Collection).11

2. Application 1 - image compressionGrayscale image matrix, each entry represents a pixel brightness.12

Grayscale image: scale 0, . . . , 255 from black to white 5255255.255255255 Colored image: 3 matrices for Red, Green and Blue brightness values13

MATLAB DEMO:Approximate a grayscale image A using the SVD by ki 1 Ai. Comparestorage requirements and quality of approximation for different k.PMemory required to store:an uncompressed image of size m n: mn valuesrank k SVD approximation: k(m n 1) values14

3. Application 2 - image deblurring15

Sources of noise and blurring: physical sources (moving objects, lens out of focus, .), measurement errors, discretization errors, rounding errors (finite precision arithmetics in computer), .Challenge: Given only the blurred noisy image B and having information about how it was obtained try to find (or approximate) exactimage X.16

Model of blurring process:Blurred photo: Barcode scanning: X(exact image) A(blurring operator) B(blurred noisy image)17

Obtaining a linear model:Using some discretization techniques, it is possible to transform thisproblem to a linear problemAx b,A Rn n,x, b Rn,where A is a discretization of A, b vec(B), x vec(X).Size of the problem: n number of pixels in the image, e.g., evenfor a low resolution 456 x 684 px we get 723 432 equations.18

Image vectorization B b vec (B):B [b1,b2,.,bw][]b b1b2.bw19

Solution back reshaping x vec (X) X:x x1x2.xw[]X [x1,x2,.,xw]20

Solution of the linear problem:Let A be nonsingular (which is usually the case). ThenAx bhas the unique solutionx A 1b.Can we simply take it?21

Image deblurring problem: Original image x, noisy blurredimage b and the “naive” solution xnaive A 1b:x true imageb blurred, noisy imagex inverse solutionWhy it does not work? Because of the properties of our problem.22

Consider that bnoise is noise and bexact is the exact part in our imageb. Then our linear model isAx b,b bexact bnoise.Usual properties: the problem is ill-posed (i.e. sensitive to small changes in b); often A is nonsingular, but its singular values σj decay quickly; bexact is smooth, and satisfies the discrete Picard condition (DPC); bnoise is often random - white (does not have dominant frequencies), and does not satisfy DPC; kbexactk kbnoisek, but kA 1bexactk kA 1bnoisek .23

SVD components of the naive solution:From the SVD A U ΣV T , we haveA 1 V Σ 1U T nXnX1 T1vj uj vj uTj .j 1 σjj 1 σjThusxnaive A 1b nX1j 1 σjn uT bexactXjσjvj uTjb vj j 1 {z}xexact A 1 bexactn uT bXjj 1 σjn uT bnoiseXjσjj 1 {zA 1 bnoisevjvj .}24

xnaive n uT bexactXjvj σjj 1 {z}exact 1exactx A bn uT bnoiseXjσjj 1 {z 1A bnoisevj}We want the left sum xexact. What is the size of the right sum(inverted noise) in comparison to the left one?Exact data satisfy DPC:exact decay faster than the singular values σ of A.On average, uTbjjWhite noise does not:noise do not exhibit any trend.The values uTjb25

T exact uT bnoise are for small indexes j dominated byThus uTjj b uj bthe exact part, but for large j by the noisy part.singular values of A and projections uTbi410right hand side projections on left singular subspaces uTbisingular values σi210noise level010 210 410 610 810 1010 12100123456784x 1026

Because of the division by the singular value, the components of thenaive solutionT bexactT bnoiseuuXkXkjjxnaive A 1b v vjjj 1j 1σjσj {z} {z}inverted noisexexactnoiseexactuTuTXNXNjbjbvj j k 1vj j k 1σjσj {z} {z}inverted noisexexactcorresponding to small σj ’s are dominated by inverted noise.Remember, there are usually many small singular values.27

Basic regularization method - Truncated SVD:Using the dyadic formA U ΣV T nXuiσiviT ,i 1we can approximate A with a rank k matrixA Sk kXi 1Ai kXui σi viT .i 1Replacing A by Sk gives an TSVD approximate solutionx(k) k uT bXjj 1 σjvj .28

TSVD regularization: removing of troublesome componentssingular values of A and TSDV filtered projections uTbi410filtered projections φ(i) uTbisingular values σi210noise levelfilter function φ(i)010 210 410 610 810 1010 12100123456784x 1029

Here the smallest σj ’s are not present. However, we removed alsosome components of xexact.An optimal k has to balance between removing noise and not loosingtoo many components of the exact solution. It depends on the matrixproperties and on the amount of noise in the considered image.MATLAB DEMO: Compute TSVD regularized solutions for different values of k. Compare quality of the obtained image.30

Other regularization methods: direct regularization; stationary regularization; projection (including iterative) regularization; hybrid methods combining the previous ones.31

Other applications: computer tomography (CT); magnetic resonance; seismology; crystallography; material sciences; .32

References:Textbooks: Hansen, Nagy, O’Leary: Deblurring Images,Spectra, Matrices, and Filtering, SIAM, 2006. Hansen: Discrete Inverse Problems,Insight and Algorithms, SIAM, 2010.Software (MatLab toolboxes): on the homepage of P. C. Hansen HNO package,Regularization tools,AIRtools,.33

Image deblurring problem: Original image x, noisy blurred image band the\naive" solution xnaive A 1b: x true image b blurred, noisy image x inverse solution Why it does not work? Because ofthe properties of our problem. 22. Consider that bnoise is noise and bexact is the exact partin our image b. Then our linear model is