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

๏ 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!

Goals and objectives

Bayesian Inference

p(θ | observations) = p(observations | θ)× p(θ) p(observations)

evidence (useful for model comparison) likelihood image formation model prior model (a priori knowledge about the observed object) parameters of interest (unknown solution)

OBJECTIVE: posterior probability density function (pdf)

• Integrate out the unwanted parameters
• Obtain the a posteriori probability distribution function (pdf)
• f all the parameters of interest
• This pdf contains both the optimum and the predictive uncertainties

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

• 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

Y B

• A

unknown known pdf known (observed data) unknown, unwanted

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

๏ 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)

X Y1 Y2

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

Y1 Y2

changes are non-stationary, so should be the the model parameters

σ0 σ1 σ2 μ1 μ2

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

๏ 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

Markov Random Field: undirected graphical model, limited spatial interactions

Forward model

• 3. Smoothness priors for natural surfaces

P(Z) ∝ e−Φ(Z,ω)

ωΣ(Zi+1,j-Zi,j)2+(Zi,j+1-Zi,j)2

The forward model (Bayesian network)

Y2

• bserved

images

underlying scene

11 22

Y1

change maps

smoothness parameters

Z

• Ref

DSM control DSM

• rientation

parameters

noise

• noise

model space

P p

control points GCP ICP

Y2

• bserved

images

underlying scene

11 22

Y1

change maps

smoothness parameters

Z

• Ref

DSM control DSM

• rientation

parameters

noise

• noise

model space

P p

control points GCP ICP

Inference: invert the forward model

• Integrate out all the nuisance variables
• Use explicit values for known parameters
• Hierarchical inference scheme

P(Z | data) P(Z | data, Θ1,Θ2)

Computing the data terms

from object space to image spaces

Image 1 Image 2 Object space (cartesian ground frame)

Z sampling grid Zk

probability density Projections

Θ1 Θ2

pixels of X

Z trajectory

I1 I2

๏ Marginalization of model variables of X+changes

• Requires the estimation of the parameters of a 2D Gaussian pdf

I1 = N(0,σ02) + N(μ1,σ12) I2 = N(0,σ02) + N(μ2,σ22) ⇒ (I1,I2) = N(M,Σ)

• Normalization term |Σ|-1/2 (one variable)

๏ Robustness issues

• The data being matched depend on the tested elevation Zk!
• Use |Σ|/Σ11Σ22 instead of |Σ|
• Should not change the behavior near the optimum

P(Zk) ∝ (Σ11Σ22/|Σ|)D/2 = (1-ρ2)-D/2

D = effective degrees of freedom ≈ 20 ρ = correlation coefficient

Computing the data terms

local elevation pdfs

Deterministic inference and approximations

Iterative optimization of an energy functional (nonlinear search: conjugate gradient, LBP, graph cuts...)

U(Z) = −log P(Z | Y1,Y2) = D(Z,Y1,Y2) + Φ(Z,ω)

data term smoothness penalty

var. vertical covar. horiz. covar.

• ptimal

Z

• ptimal DSM + uncertainties

๏ Compute the posterior marginal pdfs ๏ Gaussian approximation of the pdfs:

• ptimum & covariance matrix

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(n2) (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

Original pdf Band-limited pdf (filtered)

conserve the peak location accuracy!

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 Φ = ωZTQZ, 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 Zk
• Perform a local inversion using a fast, conjugate gradient algo.
• In the end, we get an approximate covariance matrix

Probabilistic DSM and error structure

∫ P(Z | data, Θ1,Θ2) dΘ1dΘ2

๏ 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 Sk(Θ)=f(Pk,ΣΘ)

Include the camera errors (second marginalization)

Σtotal = Σmatching + ΣC1 + ΣC2

relative errors absolute errors

Probabilistic DSM comparison

๏ DSM differencing

• Take advantage of the object space reconstruction: no resampling
• Probabilistic surface difference model: DSM1-DSM2 with Σ1+Σ2
• Hypothesis testing: “the elevation has changed”

Local difference var. σd2=σ12+σ22 |Z1-Z2|2>4(σ12+σ22) (95% confidence interval)

๏ Joint DSM analysis

• Feature matching (cliff retreat, sand dune displacement)

Least-squares matching using Σ1, Σ2

• Flow computation (erosion and transport mechanisms)

Dense 2D deformation field inference using Σ1, Σ2

Y2 Y1 Y3 Y4

date t1 date t2

Camera attitude errors

๏ Inertial Navigation Systems are not perfect...

• The localization accuracy is usually sufficient (a few cm)
• Attitude angles drift, accuracy between 0.005°-0.01°

(a few pixels or 50cm on the ground)

๏ Relative camera attitude calibration

• Re-estimate the orientation angles of camera 2 (camera 1 fixed)

using a Bayesian framework

• In a multi-date framework (2 stereo pairs) re-estimate the angles
• f cameras 2, 3, 4 (camera 1 fixed)
• The absolute error only comes from camera 1,

and it will not influence DSM comparison

Θ2 Θ1 Θ3 Θ4

date t1 date t2

Work in progress, future work...

Methodology:

Automatic DSM-based orientation (via matching) Bayesian GCP-based orientation Visualization: probabilistic contour maps 3D surface recovery from n images simultaneously Uncertainty validation

(sparse GCP, LIDAR points, accurate ref. DEM)

Applications:

Aerial coastal imaging: large scale change detection and monitoring for the Troia-Sines coastline (60km @20cm GSD) DEM on Mars (DEM=DSM) (Mars Express/HRSC and MRO/HiRise)

Conclusions

Advantages

Predictive accuracy for the DSM (no validation datasets) Predictive accuracy for the multi-date DSM differences (change hypothesis testing, error bars on elevation variations) Object space method (multi-date comparisons without resampling) Fast computation (deterministic method)

Limitations

20 samples may not be sufficient to do statistics! ➥ use triplet stereo or more images to reject outliers Imaging cameras don’t see through vegetation ➥ use LiDAR data!