Transcript Automated analysis of Spatial Grids

Automated analysis of Spatial
Grids
[email protected]
Part of the short course on Artificial Intelligence
at the American Meteorological Society Annual
Meeting, Jan. 2006, Atlanta.
[email protected]
1
Information from grids

Observed data are often on grids
Satellite imagery
 Ground-based radar
 Other data can be gridded


You may need to extract information
from these grids
Is there a boundary?
 Where are the troughs?

[email protected]
2
Spatial grids as images
Our spatial grids are often two
dimensional
 Automated extraction of information from
images is a long-established field in
machine intelligence

We can treat spatial grids as images
 And apply image processing techniques to
our data.

[email protected]
3
Image processing terminology
What you think
of it as
What they call it
Grid
Image
Grid point
Pixel
Finding
signatures
Pattern recognition
Identifying machine
parts on an assembly
line
Finding storms
Segmentation
(ditto)
Tracking storms
Optical flow
MPEG compression
Smoothing,
finding changes
Filtering
Photo-editing software
Combining
gridded data
Stereography
Robotics
[email protected]
Where is it used?
4
Other terms for this field

There are subtle differences between
image processing and

Data mining


Knowledge discovery


Relationships between identified signatures
A human draws the inferences
Computer Vision

Pattern recognition when it doesn’t work
[email protected]
5
Limitations of image processing

Scenarios where pattern recognition
succeeds:

Highly controlled environments


Photography



Red, green, blue channels can be assumed independent
IR channels can not!
Objects have sharp boundaries


Assembly lines, robotics
Storms, etc. do not
The most mature sub-fields


Filters (smoothing, edge finding)
Segmentation
[email protected]
6
Workflow of AI applications
Edge finding
Gridding
Feature
extraction
Filtering
Classification
(NN, etc.)
Segmentation

Artificial Intelligence applications that operate on spatial data usually follow these
steps:






Gridding: take spatial data and make it a locally uniform resolution grid
Filtering: smoothing, pattern matching, etc.
Edge finding: find areas of drastic spatial change
Segmentation: find contiguous areas of constant value
Feature extraction: find features based on edges/regions
Classification: classify features using a rule base, neural network, support vector
machine, etc.


Normally done with all the features from all your data fields
This will be covered by the other speakers
[email protected]
7
A note on convention

Image processing uses matrix notation
(0,0) is the top-left corner of the image
 (0,1) is the first row, second column.


As long as you don’t think of the first
index as “x” and the second index as “y”,
you will be okay

In what direction does the first index
increase?
Answer: South or downwards
[email protected]
8
Image processing assumptions

Traditionally image processing assumes a
Markov random process

Two pixels are correlated if they are adjacent to
each other



If you have pixels: a b c
a and b are correlated as are b and c
But the correlation between a and c is captured solely by
the correlation between a-b and b-c


No “second-order” correlation term.
A simplifying assumption that usually causes no
problems.

If your field is random (no correlation at all) or if the
correlation is second-order image processing techniques
may not work
[email protected]
9
Image processing assumptions

An assumption that could cause problems:



The field has uniform resolution
Projected spatial grids are not uniform resolution.
Locally, most projected grids can be assumed
uniform.


Pixels 10 km apart don’t vary that much in size.
Be careful with global operators on projected grids


Segmentation is a global operation
Filtering, edge finding, etc. are fine.

Filters are local operations.
[email protected]
10
Radar data

Radar data are NOT uniform resolution



Size of a “pixel” varies with range
Dramatic changes
Two options:

Use varying-size windows for local operators


Hard to keep track!
Convert to a Cartesian grid before applying image
processing techniques


Will cause problems for local operators at far ranges
where data are “repeated”
This is an issue with most projected data as well
[email protected]
11
Missing Data

Some times, you will have spatial data
where some pixels are missing

Can you treat this data as some value?
Say zero?
 Will simplify your task considerably
 Can dramatically speed up operations


If not, you can vary the formulae to use only
those pixels for which you have data
Filtering results may have higher quality
 Operations will be much slower

[email protected]
12
Summary: Spatial grids

To treat spatial grids as images,

Operate on the raw, un-projected data




Change “missing” data to background value if
available.
Use local operators (e.g: filters) freely


Less chance of repeated values
For “gridded” data, operate on grids at the resolution
closest to half the spacing between observations.
Be careful with global operations (e.g: segmentation)
Convert computed features (such as area) to
projected systems
[email protected]
13
Gridding
From insitu observations to grids
[email protected]
14
In-situ observations

Some environmental
data are measured by
in-situ instruments



Not remotely sensed
These measurements
are at geographically
dispersed points
Need to be converted
into grids
[email protected]
15
Gridding observations

In-situ observations have to be
interpolated on to grids to process as
images
What is the pixel resolution?
 What interpolation techniques are
available?

[email protected]
16
Pixel resolution
If the chosen pixel resolution is too
coarse, lose some of the observations
 If the chosen pixel resolution is too fine,
strong gradients may not be apparent
 A good rule of thumb:


Choose pixel resolution to be half the mean
distance between observation points
[email protected]
17
Interpolation methods

Interpolation methods:
Cressman analysis
 Barnes analysis
 Kriging

I_xy
I_i
[email protected]
18
Cressman analysis

The value at each pixel of grid:
R2 ( x xi )2 ( y  yi )2
Ixy 
 I i R2 ( x xi )2 ( y  yi )2
i

R2 ( x xi )2 ( y  yi )2
2
2
2
i R ( x xi ) ( y  yi )

Every observation gets weighted based on its
distance from grid point


R is the “radius of influence”
Higher the R, the more distant points are
considered
[email protected]
19
Barnes analysis

The problem with a Cressman analysis:

Even if a pixel is collocated with an
observation, the final value of pixel will be
different


Due to (small) effect of pixels further away
Barnes analysis tries to account for this
[email protected]
20
Barnes analysis equation

Barnes also changed
the Cressman
interpolation function
to be exponential

Smoother than
simple polynomials
[email protected]
( x  xi )2 ( y  yi )2
Ixy 
 Iie
2  2
i
( x  xi )2 ( y  yi )2
e
2  2
i
21
Barnes analysis: the technique

Perform Cressman analysis on obs


Compute errors at observation points
Perform Cressman analysis on errors




Give the error field a weight 
Add the two fields
One pass of Barnes analysis
Can repeat Barnes analysis


N-pass Barnes analysis
Closer and closer to the value of the observations
at the observation points.
[email protected]
22
Kriging

Kriging is a principled method of interpolation




Compute correlation between observations
Use variogram to determine scaling parameter
Works only if the observations are very close to
each other
M. A. Oliver and R. Webster "Kriging: a method of
interpolation for geographical information system",
Int’l. J. Geographical Information Systems, 1990,
VOL. 4, No. 3, 313-332
[email protected]
23
Convolution Filters
A surprisingly powerful local
operation
[email protected]
24
Smoothing: an intuitive look

How would you
smooth the image on
the top right?


So that the edges
aren’t quite so
ragged?
To get an image like
bottom right?
[email protected]
25
Averaging is a smoothing operation

The most obvious way to
smooth is to replace each
pixel by a local average:
Ixy  ( 2 k 1).(1 2 k 1)

i xk j  y k
 I
ij
i  xk j  y k
I used k=3 in the WDSS-II
filter “average”
 http://www.wdssii.org/
 The dataset I’m using for
demonstration is available
online. See:
 http://cimms.ou.edu/~laks
hman/aitut/
[email protected]
26
How it works

What does this equation mean?
Ixy  ( 2 k 1).(1 2 k 1)




i xk j  y k
 I
ij
i  xk j  y k
For every pixel x,y, look in a 2D
neighborhood, add up all the
values and divide by the number
of pixels.
(2k+1) is the neighborhood size
k is the half-size
Because the neighborhood size
needs to be odd, most software
will ask you for the half-size, k
[email protected]
27
Convolution

Another way to think
about smoothing:



At each pixel, you
assign weights to the
neighboring pixels.
And compute a
weighted average.
W is a (2k+1)x(2k+1)
matrix shown on right.

Called the convolution
matrix or kernel

Technically, this is cross-correlation, but if you
are using symmetric kernels (as we almost
always will be), cross-correlation and
convolution are the same thing.
Ixy 
i xk j  yk
 W I
ij ij
i  xk j  y k
1 1 1
1

W  1 1 1
9
1 1 1
[email protected]
28
Smoothing kernels

The “boxcar” kernel




A better convolution kernel is
where the kernel function drops
off with distance.




Poor characteristics in
transition areas.
More spatial relocation of
maxima
Smudges more
Gauss:9
A Gaussian convolution kernel
offers the best
space/frequency trade off
Frequency = raggedness
Space = location
You can smooth a lot more
effectively because the features
don’t move or smudge as much.
Average:3
[email protected]
29
Gaussian Kernel
( x )2 ( y )2
Wxy 


2
e
 2 x 2
x and y range from –k to k.



1
Since the Gaussian has infinite range, you need to choose k
appropriately.
Choosing k to be 3 * largest sigma is a reasonable choice.
The sigmas are the scale

Divide by total weight (The W above doesn’t add up to 1)
[email protected]
30
Convolution kernels

By changing the
weights in the
convolution kernel,
you can extract
many features.
1 0  1
1

W  1 0  1
3
1 0  1
[email protected]
31
Matching filters


In general, make your convolution kernel be
the feature you are trying to extract.
1 0  1
1
W  1 0  1
3
1 0  1
The above kernel returns large magnitudes for
thin vertical lines



But not just thin vertical lines
Need to normalize by smoothed value also
Can use the matched filter idea to find objects
[email protected]
32
Scale and rotation

Convolution filters are:

Scale dependent



Orientation dependent



What if the gradient is 5 pixels apart?
Need a matching filter with columns 5 apart!
What if the boundary is horizontal?
Need a matching filter with that orientation.
Remember typical usage?


Do you know the size of a machine part before
hand?
How about a squall line?
[email protected]
33
Scale and orientation
( x cos  )2 ( y sin  )2
Wxy

1

e
2
2
 x  y
2
 ( x 2  y
)
If you specify your convolution filter as Gaussian functions, it is
easy to specify rotation and scale.
 x and y range from –k to k.





Since the Gaussian has infinite range, you need to choose k
appropriately.
Choosing k to be 3 * largest sigma is a reasonable choice.
Theta is the orientation from x-axis (anti-clockwise from south)
The sigmas are the scale (use same sigma for circular
kernels; different sigmas for elliptical ones)
Divide by total weight (The W above doesn’t add up to 1)
[email protected]
34
Speed

The larger the scale of interest, the larger the
half-size of the kernel.



Will take a long time to compute.
Since we need to compute weighted average at
every pixel.
Two speed ups possible


Convolution can be performed in Fourier Transform
space.
Convolution can be performed using separable
filters.
[email protected]
35
Convolution with Fourier Transforms

To convolve with Fourier Transforms:

Compute 2D FFT of original image



Compute 2D FFT of convolution kernel




The convolution kernel has to be same size as image (pad with zeros)
Multiply the FFT values pixel-by-pixel
Take inverse 2D FFT of result.
The larger your kernel, the more dramatic the speed up.



Can not have missing data
Fill with background value somehow
Can be 100-200 times faster.
Seconds instead of hours!
Resources


fftw is a free, nice Fourier Transform library.
Lakshmanan, V: 2000: Speeding up a large scale filter. J. of Oc. and
Atm. Tech., 17, 468-473
[email protected]
36
Convolution with separable filters

What if you have missing data?



Can not use FFT method
Problem with radar data in particular.
If the convolution filter is “separable”:



Compute weighted average row-wise on original
image.
Compute weighted average on the row-wise result
column-by-column.
If your filter was 129x129, 60 times faster than
doing the standard way.
[email protected]
37
What is a separable filter

A filter is separable if

For example, if the
convolution filter is a
simple Gaussian:

Wxy  WxWy
x2  y 2
Wxy 
1

2
e
 2
x2

1

e
 2
y2
1

e
 2
If you want to vary
orientation, see:

Lakshmanan, V 2004: A
separable filter for
directional smoothing.
IEEE Geosc. and
Remote Sensing
Letters, 1, 192-195
[email protected]
38
Edge finding

The “line-finding” filter is not a
robust way to find edges


Similar problems to using
boxcar for smoothing
1 0  1
1
W  1 0  1
3
1 0  1
Use Laplacian of a Gaussian
(Canny)




Differentiate the smoothing
kernel twice.
Convolve image with that
kernel
Look for zero crossings in the
convolution result.
J. Canny, ``A computational
approach to edge detection",
IEEE Trans. on Pattern
Analysis and Machine
Intelligence, 8(6),pp679-698
(1986)
1 x y 
Wxy  ( 
)e
2
2

[email protected]
2
2
x2  y 2
2
39
Multi scale or orientation analysis

Use filter banks:




Drawbacks:




To perform multi scale (multi resolution) or orientation analysis
A bunch of filters at different scales and orientations
Combine their results by a weighted average or maximum
Repeated convolutions on original image.
Results of the filters are not related.
Can get very expensive
Simplifications possible:


We can steer the results based on a few filters
The resulting features are related



The squall line contains storm cells
So far, we have only considered image-to-image
Filter 1
The resulting images are related.

Use wavelets
Image
Filter 2
Pixel
operation
Result
Filter N
[email protected]
40
Wavelets

Wavelet analysis



breaks up images
smoother and smoother results on the right
(s0,s1)
The lost “detail” information on the left






D1 (N/2 x N/2)
3
So that s1 is half the size of s0
Convolution filter has to be a wavelet function


3 images of the size of the corresponding s
Can not use any old filter

S0 (NxN)
Haar function (boxcar)
Daubechies p-functions
Many others
D2 (N/4 x N/4)
3


S2 (N/4 x N/4)
Fast Wavelet transformations exist
Use wavelets only if you need this
decomposition

S1 (N/2 x N/2)
D3 (N/8 x N/8)
3
S3 (N/8 x N/8)
Otherwise, use simple multi scale analysis
with lots of Gaussian filters.
Gaussian filters have lots of nice properties.
Plus everybody knows what they are!
[email protected]
41
Texture
In convolution, we compute weighted
average of the original image’s pixels
 Can also do other operations inside
neighborhood

Texture operators
 Operators based on relative ordering

[email protected]
42
Texture operators

Texture is a good way to distinguish different objects



2nd and 3rd order statistics



Not just on data value, but on distribution of data values
Lakshmanan, V., V. DeBrunner, and R. Rabin, 2000: Texturebased segmentation of satellite weather imagery. Int'l
Conference on Image Processing, Vancouver, 732-735.
Variance, standard deviation
Homogeneity
Entropy



Shannon’s measure of information content
Bin the data values and form a histogram
P is the fraction of pixels in each bin
[email protected]
 p log( p )
i
i
i
43
Percent filters

Operators based on sorting

Take the median value



Take the maximum value



Excellent smoothing filter
Optimal when you have speck noise
Effect of dilating the field spatially
Minimum erodes the field spatially
Use percentiles




50% is the median (maximum noise reduction)
0% is the min (heavily noise dependent)
100% is the max
Intermediate values trade off on noise reduction
[email protected]
44
Dilation and Erosion filters

Using the second-highest or
second-lowest value in the
neighborhood is an excellent
dilation/erosion filter.
 Implemented this way in
WDSSII.
erode
[email protected]
dilate
orig
45
Filter Banks and Neural Networks

Classification using NN






Use filter banks
Convolution filters Image
Texture operators
Percentile operators
Each pixel provides a
pattern
May have to select pixels –
not use all of them
 Since a 1000x1000
image has 1 million
patterns.
 But is only data case
[email protected]
Filter 1
Filter 2
Neural
Network
Result
Filter N
46
QCNN


Features computed in local
neighborhoods of radar gates Effect of
Pre-processing
(don’t care)
Classified by NN
 Returns [0,1]
 Is it perfect?


If we use this directly as
a mask, there will be
holes in the result!
Lakshmanan, V., A. Fritz,
T. Smith, K. Hondl, and
G. J. Stumpf, 2005: An
automated technique to
quality control radar
reflectivity data. J. Applied
Meteorology, subm.
Original
NN random
error
QC’ed
Effect of post
processing
[email protected]
47
Actual QC technique
Original
Effect of
Pre-processing
(don’t care)
NN random
error
QC’ed
Effect of post
processing
[email protected]
48
Images to Objects

So far, we have looked at ways of
processing images

Result is still image


Bunch of pixels
Need to group pixels together into objects

Classify, track, etc. these objects if possible.
[email protected]
49
Segmentation
From images to objects
[email protected]
50
What is segmentation?
Finding groups of pixels that together
comprise an object.
 This is how we (humans) make sense of
most images.
 Very hard for a computer algorithm.

Simple region-growing approach
 Complex vector segmentation approach

[email protected]
51
Binary images

Region growing can be applied only to
binary images.
Apply to images that can be thresholded
 Pixels below (above) a threshold are not of
interest.
 Your image now has only 0’s (not
interested) and 1’s (interested)


Threshold dependent on application
[email protected]
52
Hysterisis

Alternative is to employ “hysterisis”




You are interested in areas > 30 dBZ (T_i –
threshold of interest)
But if you see a 30 dBZ pixel, you lose interest only
if pixel value < 20 dBZ (T_d – threshold of drop off)
Less chance of too many disconnected regions.
Employed in QCNN


Lakshmanan, V., A. Fritz, T. Smith, K. Hondl, and G. J.
Stumpf, 2005: An automated technique to quality control
radar reflectivity data. J. Applied Meteorology, subm.
Choose T_i and T_d by experiment.
[email protected]
53
Region growing algorithm

Initialize:



label, an image with all pixels = 0.
curr_label = 0
Walk through image pixel-by-pixel.

If label(x,y) = 0 and image(x,y) > T_i



Region growing steps:


Set label(x,y) = curr_label
Check all 8 neighbors (i,j) of (x,y)


Increment curr_label
Perform region growing at x,y with curr_label
If image(i,j) > T_d, perform region growing at i,j
Two pixels with the same label are part of the same region.

Can compute region properties based on all the pixels within a
region.
[email protected]
54
Vector segmentation

Region growing:

Can lead to disconnected regions


Not hierarchical


Even with hysterisis
Can not always use hierarchy of thresholds
Need a segmentation approach that:
Can avoid disconnected regions
 Is hierarchical

[email protected]
55
Watershed Segmentation

Watershed segmentation
 Sort pixels in ascending order
of magnitude; then, flood the
plain:



Vincent, L., and Soille, P.,
Watersheds in Digital
Spaces: An Efficient
Algorithm Based on
Immersion Simulations,
PAMI(13), No. 6, June 1991,
pp. 583-598.
Use saliency to form hierarchy
Watershed segmentation
works very poorly in images
with statistical noise (most
meteorological images).
[email protected]
56
Contiguity-enhanced segmentation

Explicitly trade-off contiguity of region with the
idea that a region’s pixels ought to be similar




Similarity of pixels based on a number of texture or
wavelet features
Pixels are dissimilar if vector distance between
features is large
Lakshmanan V. 2001: A Heirarchical, Multiscale
Texture Segmentation Algorithm for Real-World
Scenes. Ph.D. thesis, U. Oklahoma, Norman, OK
The WDSS-II algorithm w2segmotionll
implements this segmentation approach

http://www.wdssii.org/
[email protected]
57
K-Means Segmentation
Lakshmanan, V., R. Rabin, and V. DeBrunner 2003: Multiscale
storm identification and forecast. J. Atm. Res., 367-380
[email protected]
58
Hierarchical Segmentation
[email protected]
59
Incorporating time
[email protected]
60
Frames and sequences
A spatial grid at a point time is called a
“frame”
 A bunch of frames form a sequence
 Terminology from video

Focus of most image processing research
 Used for security or monitoring applications

[email protected]
61
Detecting change

To detect change between frames

Simply difference the images pixel-by-pixel
Captures edges of moving objects
 Works well on video



Truck enters scene
Doesn’t work well when the objects themselves
are changing

Such as storms
[email protected]
62
Object association

One alternative is to associate objects


Segment the frames
Try to associate objects in one frame to next based
on:




object similarity (use texture features)
object overlap (fraction of pixels in common)
Closeness of centroids
A linear programming problem

Dixon, M. and Wiener, G., 1993: TITAN: Thunderstorm
Identification, Tracking, Analysis and Nowcasting: A Radar-Based
Methodology. J. Atmos. Oceanic. Tech, 10, 6
[email protected]
63
Object association (contd.)

Advantages of object association:
Can extract underlying trends.
 Can advect each storm individually.


Disadvantages of object association:
Segmentation results can change from
frame to frame.
 Splits/merges can not be handled easily.
 Accurate only for small storms
 Works only at one scale (not hierarchical)

[email protected]
64
Optical flow

Another way is to treat the images are 2D fluid
fields

Find movement that minimizes some cost measure
to move from one frame to another.

Cross-correlation


Lagrangian motion




Implemented as w2imgmotion in WDSS-II
Turner B., Zawadzki, I. and Germann, U. 2004: Predictability
of Precipitation from Continental Radar Images. Part III:
Operational Nowcasting Implementation (MAPLE), J.
Applied Meteorology, 43, 2, pp.231-248
Typically more accurate than object association.
Can extrapolate motion vector to get forecasts.
No way to extract trends for individual storms
[email protected]
65
Hybrid technique



Object association gives us tracks/trends
Optical flow gives us prediction accuracy
A hybrid technique




Identify objects at current frame
Match to image (not objects) in previous frame using cost
measure
Extract trends from objects at current frame and images from
previous frames based on movement of object
Lakshmanan, V., R. Rabin, and V. DeBrunner 2003:
Multiscale storm identification and forecast. J. Atm.
Res., 367-380

Implemented as w2segmotion in WDSS-II
[email protected]
66
Hybrid technique (contd.)

Advantage of hybrid technique

Can use shape of identified storms to track
other fields
Identify storms on radar images
 Track satellite cloud-top temperature within
storm


Possible because the exact shape of the
storm is known
Though it changes from frame to frame.
 So item of interest should not be at edges.

[email protected]
67
GIS systems

Trends can include geographic
properties
Just grid the GIS data into the same
resolution at your spatial grids
 Then, you can find the total population
affected by a storm, for example.
 Treat the GIS data as just another spatial
grid.

[email protected]
68
Polygons

Gridding the GIS data is not always
practical.
County boundaries are polygons.
 Want to find things that happen within that
polygon
 Maximum hail within that polygon


How can we extract spatial properties
within a earth-based polygon?

Implemented as w2alarm in WDSS-II
[email protected]
69
Polygon in/out tests

Put the polygon in the gridded
reference



The polygon corners have to
be converted to (x,y) locations.
We’ll assume the polygon is
completely inside the grid.
Walk row-wise



3
2
1
from the grid boundary to the
pixel of interest.
count the number of times you
cross the sides of the polygon


1
Odd: the pixel is inside the
polygon
Even: the pixel is outside the
polygon
2
0
Repeat column-wise

Has to be inside in both
directions.
[email protected]
70
Summary
Edge finding
Gridding
Feature
extraction
Filtering
Classification
(NN, etc.)
Segmentation

We looked at how image processing plays an important role in most AI
applications on spatial grids.




Filtering
Finding edges
Finding objects
Also at purely image-driven applications


Tracking
Extraction of spatial properties.
[email protected]
71
Thank you!
http://www.cimms.ou.edu/~lakshman/
[email protected]
http://www.wdssii.org/
[email protected]
72