Page 1 Inverse Problems - Examples Inverse Problems - Theory - PDF document

Outline � Theory � example 1D deconvolution � Fourier method � Algebraic method Inverse Problems � discretization � matrix properties � regularization � solution methods Ivo Ihrke � Computed Tomography (CT) � Radon transform � Filtered Back-Projection � natural phenomena � glass objects Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Inverse Problem - Definition Inverse Problems - Examples � forward problem � forward problem – volume rendering � given a mathematical model M and its � given voxel data and image formation model, parameters m, compute (predict) observations o compute a view of the object � m are the volume coefficents, c is o = M ( m ) the ray that determines the pixel's � inverse problem value o (observation) � � given observations o and a mathematical model o = c m ( c ( s )) ds M, compute the model's parameters m = M − � ( o ) Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Inverse Problems - Examples Inverse Problems - Examples � inverse problem – CT � forward problem – convolution � given the pixel values o, the ray geometry c and � example blur filter the image formation model, compute the volume � given an image m and a filter kernel k, compute densities m the blurred image � invert 3D � o = m ⊗ k o = m ( c ( s )) ds + n c � n is a noise component � we will later see how to do this Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Page 1

Inverse Problems - Examples Inverse Problems - Theory � inverse problem – deconvolution � deconvolution in Fourier space � convolution theorem ( F is the Fourier transform ): � example blur filter � given a blurred image o and a filter kernel k, compute the sharp image o = m ⊗ k ⇔ F ( o ) = F ( m ) � F ( k ) � need to invert F ( m ) = F ( o ) o = m ⊗ k + n � deconvolution: F ( k ) � problems � n is again noise � division by zero � Gibbs phenomenon (ringing artifacts) Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 A One-Dimensional Example – Deconvolution A One-Dimensional Example – Deconvolution Spectral Spectral F ( k ) � most common: is a low pass filter reconstruction is noisy even if data is perfect ! � � � � � , the inverse filter, is a high pass filter F � k � � amplifies noise and numerical errors � � � Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 A One-Dimensional Example – Deconvolution A One-Dimensional Example - Deconvolution Spectral Spectral spectral view of signal, filter and inverse filter � solution: restrict frequency response of high pass filter (clamping) � � � , if | F � k � | < γ F � k � G = γ F � k � , else | F � k � | M = O � G Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Page 2

A One-Dimensional Example - Deconvolution A One-Dimensional Example- Deconvolution Spectral Algebraic � alternative: algebraic reconstruction reconstruction with clamped inverse filter � convolution � ∞ o ( x ) = m ( t ) k ( x − t ) dt −∞ � discretization: linear combination of basis functions N � m ( t ) = m i φ i ( t ) i Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 A One-Dimensional Example – Deconvolution A One-Dimensional Example – Deconvolution Algebraic Algebraic � discretization: � discretization: o = m ⊗ k o = m ⊗ k � observations are � observations are � ∞ � ∞ linear combinations linear combinations = m ( t ) k ( x − t ) dt = m ( t ) k ( x − t ) dt of convolved basis of convolved basis −∞ −∞ functions functions � ∞ � ∞ � N � N ≈ m i φ i ( t ) k ( x − t ) dt ≈ m i φ i ( t ) k ( x − t ) dt � linear system with � linear system with m i m i −∞ −∞ unknowns unknowns i i � ∞ � ∞ � N � N � often over- � often over- = m i φ i ( t ) k ( x − t ) dt = m i φ i ( t ) k ( x − t ) dt determined, i.e. determined, i.e. −∞ −∞ i i more observations o more observations o unknown N N � � than degrees of than degrees of = m i ( φ i ⊗ k ) = m i ( φ i ⊗ k ) freedom ( basis freedom ( basis i i functions ) functions ) � = �� � = �� linear system linear system Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 A One-Dimensional Example – Deconvolution A One-Dimensional Example – Deconvolution Algebraic Algebraic � normal equations – in case you forgot � Why ? � analyze distribution of eigenvalues � ( �� − � ) T ( �� − � ) = min � || �� − � || � min � = min � f ( � ) � remember ▽ f = 2 � T �� − 2 � T � = 0 det( � ) = Π N i �� λ i , det( � ) = 0 ⇒ matrix � is under�determined � T �� = � T � � solve to obtain solution in a least � actually we will check the singular values squares sense � T � (square root of eigenvalues of ) � apply to deconvolution problem solution is completely broken ! Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Page 3

A One-Dimensional Example – Deconvolution A One-Dimensional Example – Deconvolution Algebraic Algebraic � T � � � Why is this bad ? matrix has a very wide range of singular values! � more than half of the singular values are smaller than machine � Singular Value Decomposition: U, V are ≈ 10 − �� epsilon for double precision orthonormal, D is diagonal 10 � � = ��� T ( ��� T ) − � � − � � Inverse of M: = � − T � − � � − � = 10 − � �� − � � T = Log-Plot ! � singular values are diagonal elements of D � − � = diag( � inversion: � � �� ) 10 − �� Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 A One-Dimensional Example – Deconvolution A One-Dimensional Example – Deconvolution Algebraic Algebraic � computing model parameters from observations: � inverse problems are often ill-conditioned � = � − � � = �� − � � T � (have a numerical null-space) � inversion causes amplification of noise � again: amplification of noise 10 � � potential division by zero 10 �� 10 − � Log-Plot ! 10 � � − � = diag( � � �� ) numerical null space 10 − �� 10 � Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Well-Posed and Ill-Posed Problems Condition Number � Definition [Hadamard1902] � measure of ill-conditionedness: condition number � measure of stability for numerical inversion � a problem is well-posed if � ratio between largest and smallest singular value a solution exists 1. the solution is unique 2. κ ( � ) = σ � σ � , σ � > . . . > σ N are the singular values of � the solution continually depends on the data 3. � a problem is ill-posed if it is not well-posed � smaller condition number � � less problems when � � � most often condition (3) is violated inverting linear system if model has a (numerical) null space, parameter � � condition number close to one implies near choice influences the data in the null-space of the data orthogonal matrix very slightly, if at all � noise takes over and is amplified when inverting the model Computational Photography Ivo Ihrke, Summer 2007 Computational Photography Ivo Ihrke, Summer 2007 Page 4