Title: Grayscale Morphology in 1D Author: Thomas Breuel Institution: UniKL
from pylab import *
from scipy import stats
from scipy.ndimage import morphology,filters
Let's construct a detection problem. As usual, we have three parts:
Construction of background signal.
background = filters.gaussian_filter(randn(1000),10.0)
background -= amin(background)
background /= amax(background)
plot(background)
[<matplotlib.lines.Line2D at 0x7a48bd0>]
Construction of template.
template = concatenate([linspace(0,0.2,100),linspace(0.2,0,100)])
template += 0.3
plot(template)
[<matplotlib.lines.Line2D at 0x7cbfad0>]
Construction of combined signal (template + background + noise).
combined = background.copy()
combined[200:400] = template
combined += randn(len(combined))*0.02
combined = filters.gaussian_filter(combined,2.0)
plot(combined)
[<matplotlib.lines.Line2D at 0x7ced150>]
The general idea behind detection is to...
def window(a,r,i):
return roll(a,r//2-i)[:r]
A common way of determining similarity is by looking at the difference of the window $w_i$ from the template $t$ as vectors.
$$ d_i = ||w^{(i)}-t||^2 $$result = [sum((window(combined,r,i)-template)**2)**.5 for i in range(n)]
plot(result)
print argmin(result)
301
We can also look at the sum of the absolute differences; it's a slightly different measure that's more robust to outliers.
$$ d_i = \sum_j |w^{(i)}_j-t_j| $$sad = [sum(abs(window(combined,r,i)-template)) for i in range(n)]
plot(sad)
print argmin(sad)
300
The maximum difference is very sensitive to outliers, but still outputs the correct location. This forms the basis of mathematical morphology.
$$ d_i = \max_j |w^{(i)}_j-t_j| $$mad = [amax(abs(window(combined,r,i)-template)) for i in range(n)]
plot(mad)
print argmin(mad)
301
We can also use a percentile filter on the distribution of differences, which is a little more robust.
$$ d_i = \hbox{perc}_{90}\\{ |w^{(i)}_j-t_j| \\}$$result = [stats.scoreatpercentile(abs(window(combined,r,i)-template),90) for i in range(n)]
plot(result)
print argmin(result)
301
We might also be interested in matching regardless of a constant offset of the template. We can achieve this by subtracting out, for example, the difference between the average of the window and template.
result = [sum(abs(window(combined,r,i)-template-(mean(window(combined,r,i))-mean(template)))) for i in range(n)]
plot(result)
print argmin(result)
300
We used the maximum of the difference between the window and the signal as a quality of match measure:
$$ d_i = \max_j |w^{(i)}_j-t_j| $$We can actually break this up into two pieces to get a closely related measure. First, we measure how far above the template the signal is:
$$ d_i = \max_j (w^{(i)}_j-t_j) $$We can do the equivalent for the minimum:
$$ d_i = \min_j (w^{(i)}_j-t_j) $$You can think of this as taking the template and moving it up/down as far as possible until it touches the graph of the signal somewhere.
# upper difference
upper = array([amax(window(combined,r,i)-template) for i in range(n)])
plot(result)
print argmin(result)
163
# lower difference
lower = array([amin(window(combined,r,i)-template) for i in range(n)])
plot(result)
print argmin(result)
163
We can now find a quality of match as a difference between the upper and lower differences.
Note that this is not exactly the same as the maximum of the differences over the window, but it also serves to localize the template.
diff = abs(upper)+abs(lower)
plot(diff)
plot(mad)
print argmin(diff)
301
If the template is offset (in value) relative to the signal, then the upper and lower matches are also going to be offset the same way, so we can match signals by looking at the difference between the upper and lower matches.
diff = upper-lower
plot(diff)
plot(mad)
print argmin(diff)
301
# matching with offset template
upper = array([amax(window(combined,r,i)-(template+0.7)) for i in range(n)])
lower = array([amin(window(combined,r,i)-(template+0.7)) for i in range(n)])
plot(abs(upper)+abs(lower),color="red")
plot(upper-lower,color="blue")
[<matplotlib.lines.Line2D at 0xf325d50>]