PROBABILISTIC SURFACE CHANGE DETECTION AND MEASUREMENT FROM DIGITAL - PowerPoint PPT Presentation

PROBABILISTIC SURFACE CHANGE DETECTION AND MEASUREMENT FROM DIGITAL AERIAL STEREO IMAGES André Jalobeanu, Cristina Gama Geophysics Center of Évora - University of Évora, Portugal José A. Gonçalves Geospatial Sciences Research Center - University of Porto, Portugal

Outline Goals Bayesian inference Forward model The Bayesian network Data terms Inference (inversion) Probabilistic DSM DSM comparison Future work & conclusions

Goals and objectives ๏ Digital Surface Model (DSM) with error map ‣ Predict non-stationary uncertainties ➔ error propagation ‣ Compute consistent and physically meaningful error bars ‣ Dense model with sub-pixel accuracy ๏ Why we use stereo optical images ‣ Availability, coverage, redundancy, price ๏ Requirements ‣ Raw and well-sampled images, metadata or GCP or ref. DSM ๏ Necessary tools ‣ Probability theory, signal processing, computer vision, applied math, and some Physics!

Bayesian Inference prior model likelihood OBJECTIVE : (a priori knowledge image formation model about the observed object) posterior probability density function (pdf) p ( θ | observations ) = p ( observations | θ ) × p ( θ ) p ( observations ) parameters of interest (unknown solution) evidence (useful for model comparison) • Integrate out the unwanted parameters • Obtain the a posteriori probability distribution function (pdf) • This pdf contains both the optimum and the predictive uncertainties of all the parameters of interest

Basic ingredients & mathematical tools ‣ Forward modeling: • All parameters are random variables • Likelihood - image formation or forward model • Prior pdf - object modeling • Graphical models convenient for design and understanding ‣ The proposed Bayesian inference scheme: • Marginalization (integration) of nuisance variables • Approximations - otherwise intractable! • Deterministic optimization for speed requirements • Uncertainty prediction

Graphical models / Bayesian networks unknown, unknown known pdf unwanted � A B Y known (observed data) ‣ Converging arrows: conditional pdfs ‣ Terminal nodes: prior pdfs ‣ Joint pdf = Product (Priors) x Product (Conditionals) P( Θ ,A,B,Y) = P( Θ ) P(A) P(B) P(Y | Θ ,A,B) ‣ Inference, joint pdf: P( Θ ,A,B | Y) ∝ P( Θ ,A,B,Y) ‣ Inference, marginal pdf : P( Θ | Y) = ∫ P( Θ ,A,B | Y) dAdB ‣ Use the graph structure for an efficient, hierarchical inversion

Forward model 0. Finite resolution surface model ๏ Target = band-limited DSM ‣ No information beyond the cut-off frequency of the system (finite resolution optical images) ‣ Nyquist-Shannon sampling theorem A continuous DSM can be reconstructed from the samples (if some conditions are satisfied) ‣ Spline interpolation , if need to interpolate: good compromise between accuracy and computational complexity ‣ The predicted error is with respect to a band-limited ideal DSM , not the infinite resolution “ground truth” Spline kernel Representation - sum of splines

Forward model 1. Underlying 2D reflected radiance map X Y 1 Y 2 ๏ Common reflected radiance map X ๏ Model common/change maps using Gaussian processes ‣ X : Spatially uncorrelated, zero-mean iid Gaussian noise N(0, σ 02 ) ‣ Additive change maps : same with mean μ , N( μ i , σ i2 ) ๏ Warping and resampling scheme ‣ Resampling via B-Spline interpolation ‣ Uniform vertical shift (small patch assumption)

Forward model 2a. Radiometric change model ๏ Spatially adaptive parametric model ‣ Multiplicative changes - include reflectance effects (non-Lambert, lighting variations), shadows, atmospheric attenuation, instrumental artifacts... ‣ Additive changes - include atmospheric haze, clouds, instrumental biases... ‣ Additive noise - approx. Gaussian, independent pixels changes are non-stationary, so should be the the model parameters σ 0 σ 1 σ 2 μ 1 μ 2 Y 2 Y 1

Forward model 2b. Parameter density and overlap ๏ Statistics : large number of samples ‣ Use large windows to estimate the Gaussian parameters ๏ Goal : high resolution DSM ‣ Don’t degrade the resolution too much! ‣ Solution: reduce the spacing between windows, create overlap ๏ Inference : independent data terms ‣ Some optimization algorithms require independent data terms ‣ Dependence if overlapping windows, unknown correlation! ‣ Independent estimation is easier! Optimal window shape/size/spacing max. resolution, min. correlation, reasonable statistics (20 DOF) Bilinear or Spline weighting, FWHM = spacing = 3 pixels

Forward model 3. Smoothness priors for natural surfaces ๏ Model the topography in the object (world) space ‣ Simplest space for modeling (vs. disparity space), no resampling required in the end ‣ Easier comparison of multi-date DSMs ‣ Natural surface modeling using self-similar process ‣ Approximation using Markov Random Fields : • spatial interaction only between nearest neighbors • energy term based on slope or curvature planetary surface modeling P(Z) ∝ e −Φ (Z, ω ) Markov Random Field: undirected graphical model, limited spatial interactions ω Σ (Z i+1,j -Z i,j ) 2 +(Z i,j+1 -Z i,j ) 2

The forward model (Bayesian network) P GCP Ref DSM smoothness parameters p ICP control orientation DSM � parameters control � � points model space Z Y 1 Y 2 � 0 observed � 1 � 1 � 2 � 2 images underlying change scene maps noise noise

Inference: invert the forward model ‣ Integrate out all the nuisance variables ‣ Use explicit values for known parameters ‣ Hierarchical inference scheme P GCP Ref DSM smoothness parameters p ICP control orientation DSM � parameters control � � points model space Z Y 1 Y 2 � 0 P( Z | data, Θ 1 , Θ 2 ) observed � 1 � 1 � 2 � 2 images underlying change scene maps noise noise P( Z | data )

Computing the data terms from object space to image spaces I 2 I 1 Z trajectory Image 1 Image 2 Projections Θ 1 Θ 2 Z k Z sampling grid pixels of X probability Object space density (cartesian ground frame)

Computing the data terms local elevation pdfs ๏ Marginalization of model variables of X+changes ‣ Requires the estimation of the parameters of a 2D Gaussian pdf I 1 = N(0, σ 02 ) + N( μ 1 , σ 12 ) I 2 = N(0, σ 02 ) + N( μ 2 , σ 22 ) ⇒ ( I 1 , I 2 ) = N(M, Σ ) ‣ Normalization term | Σ | -1/2 (one variable) ๏ Robustness issues ‣ The data being matched depend on the tested elevation Z k ! ‣ Use | Σ |/ Σ 11 Σ 22 instead of | Σ | ‣ Should not change the behavior near the optimum P( Z k ) ∝ ( Σ 11 Σ 22 /| Σ | ) D/2 = ( 1- ρ 2 ) -D/2 ρ = correlation coefficient D = effective degrees of freedom ≈ 20

Deterministic inference and approximations ๏ Compute the posterior marginal pdfs Iterative optimization of an energy functional (nonlinear search: conjugate gradient, LBP, graph cuts...) U(Z) = − log P(Z | Y 1 ,Y 2 ) = D(Z,Y 1 ,Y 2 ) + Φ (Z, ω ) data term smoothness penalty ๏ Gaussian approximation of the pdfs: optimum & covariance matrix vertical covar. optimal Z horiz. var. covar. optimal DSM + uncertainties

Optimization via Loopy Belief Propagation (LBP) ๏ Assumptions ‣ Precomputed, independent data terms (data term) ‣ First order interactions only (prior term) ๏ Exact inference on chains ‣ Forward/Backward message passing yields marginal pdfs ๏ Approximate inference on graphs with loops (MRF) ‣ Multiple sweeping for message passing (left-right, up-down) ‣ Fast convergence (less than 10 iterations) ‣ The optimum is accurate enough but the variance is not: use the data terms for uncertainty computation

Sub-pixel inference and pdf sampling ๏ Minimize information loss due to pdf sampling ‣ LBP optimization requires discrete probabilities: P(Z=k) ‣ The LBP complexity is O(n 2 ) (n=number of Z samples) ➥ work with arbitrary sampling Δ , e.g. ≈½ pixel ➥ band-limit (blur) before sampling to avoid information loss ‣ Computing correlations is expensive ➥ compute local Σ on ≈½ pixel grid then interpolate conserve the Original pdf peak location accuracy! Band-limited pdf (filtered)

From precision matrix to covariance matrix ๏ Data contribution ‣ Diagonal elements only ‣ Log data term pdfs ‣ Gaussian fitting or 2nd derivatives at the optimum ๏ Prior contribution ‣ If Φ = ω Z T QZ, then the prior precision matrix is 2 ω Q where Q=4 for diagonal elements, and -1 for near-diagonal ones ๏ Sparse, but very large matrix ‣ Restrict to the neighborhood of the current element Z k ‣ Perform a local inversion using a fast, conjugate gradient algo. ‣ In the end, we get an approximate covariance matrix

Probabilistic DSM and error structure Include the camera errors ∫ P( Z | data, Θ 1 , Θ 2 ) d Θ 1 d Θ 2 (second marginalization) ๏ Additive errors (independence assumption) ‣ Area matching errors: Markov Random Field structure ‣ Camera-related errors: fully correlated (functions of the same variables): will not contribute to spatial derivative errors (e.g. slopes) Σ C = S( Θ )S( Θ ) T where S k ( Θ )=f(P k , Σ Θ ) Σ total = Σ matching + Σ C1 + Σ C2 absolute relative errors errors