Ex11

Author:Jan Beeck

Correlation normalization in the frequency domain

For the correlation normalization in the frequency domain we are going to use the following equations (for image-processing applications in which the brightness of the image and template can vary due to lighting and exposure conditions):

where , and are the fourier transform of each signal (signal i is made of ones).

To calculate in the frequency domain we are going to use what the convolution theorem states:

where is the fourier transform.

 1 from ia636 import iapconv, iadft, ianormalize, iaidft, iaind2sub
 2 
 3 def correlation_norm(f,w):
 4     w1 = w[::-1,::-1]
 5     g = iapconv(f, w1)
 6     i = ones(shape(w1))
 7     fm2 = iapconv(f*f, i)
 8     g2 = g/sqrt(fm2)
 9 
10     return g2
11 
12 def correlation_norm_fd(f,w): #in the frecuency domain, w is the mask
13     w1 = w[::-1,::-1]
14     g = iapconv(f, w1)
15     i = ones(shape(w1))
16     FH = iadft(g)
17     ffi = iapconv(f*f,i)
18     FFI = iadft(ffi)
19     g1 = iaidft(FH)/sqrt(iaidft(FFI))
20 
21     return g1
22 
23 def correlation_norm_fd_fft(f,w): #in the frecuency domain, w is the mask
24     w1 = w[::-1,::-1]
25     g = iapconv(f, w1)
26     i = ones(shape(w1))
27     FH = fft.fft2(g)
28     ffi = iapconv(f*f,i)
29     FFI = fft.fft2(ffi)
30     g1 = fft.ifft2(FH)/sqrt(fft.ifft2(FFI))
31 
32     return g1
33 
34 
35 f = iaread('cameraman.pgm')
36 adshow(f, title='input image')
37 w = f[25:25+17,106:106+17]
38 adshow(w, title='mask')
39 
40 t1=time.time()
41 f1 = correlation_norm(f,w)
42 t2=time.time()
43 t = t2 - t1
44 print "Time correlation norm with convolution:", t
45 v, pos = max(ravel(f1)), argmax(ravel(f1))
46 (row, col) = iaind2sub(f1.shape, pos)
47 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
48 adshow(ianormalize(log(abs(f1)+1), [0,255]), title='correlation normalization - spatial domain')
49 
50 t1=time.time()
51 f2=correlation_norm_fd(f,w)
52 t2=time.time()
53 t = t2 - t1
54 print "Time correlation norm in the frecuency domain:", t
55 v, pos = max(ravel(f2)), argmax(ravel(f2))
56 (row, col) = iaind2sub(f2.shape, pos)
57 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
58 adshow(ianormalize(log(abs(f2)+1), [0,255]),title="correlation normalization in the frequency domain")
59 
60 t1=time.time()
61 f2=correlation_norm_fd_fft(f,w)
62 t2=time.time()
63 t = t2 - t1
64 print "Time correlation norm in the frecuency domain with fft:", t
65 v, pos = max(ravel(f2)), argmax(ravel(f2))
66 (row, col) = iaind2sub(f2.shape, pos)
67 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
68 adshow(ianormalize(log(abs(f2)+1), [0,255]),title="correlation normalization in the frequency domain with fft")
Time correlation norm with convolution: 0.426994085312
found best match at ( 37,122)

Time correlation norm in the frecuency domain: 1.38098883629
found best match at ( 37,122)

Time correlation norm in the frecuency domain with fft: 0.444134950638
found best match at ( 37,122)

input image

mask

correlation normalization - spatial domain

correlation normalization in the frequency domain

correlation normalization in the frequency domain with fft

Index correlation in the frequency domain

For the index correlation in the frequency domain we are going to use the following equations:

where , represents the fourier transform, the image, is the mask, is the mean where match with the mask and are the mean values of the mask.

To calculate in the frequency domain we are going to use what the convolution theorem states:

where is the fourier transform.

 1 from ia636 import iapconv, iaind2sub, ianormalize, iadft, iaidft
 2 
 3 def index_correlation(f,w):
 4     i = ones(shape(w))
 5     fc = iapconv(f,i)/product(w.shape)
 6     wc = mean(w)
 7     f1 = f - fc
 8     w1 = w - wc
 9     w1 = w1[::-1,::-1]
10     g = (iapconv(f1,w1))/(sqrt(iapconv(f1**2,w1**2)))
11 
12     return g
13 
14 def index_correlation_fd(f,w): #in the frequency domain
15     i = ones(shape(w))
16     fc = iapconv(f,i)/w.shape[0]
17     wc = mean(w)
18     f1 = f - fc
19     w1 = w - wc
20     F = iadft(f1)
21     W = iadft(w1)
22     x = iapconv(F,W)
23     w1 = w1[::-1,::-1]
24     F1 = iadft(f1)**2
25     W1 = iadft(w1)**2
26     y = iapconv(F1,W1)
27     g = iaidft(x)/sqrt(iaidft(y))
28 
29     return g
30 
31 def index_correlation_fd_fft(f,w): #in the frequency domain
32     i = ones(shape(w))
33     fc = iapconv(f,i)/w.shape[0]
34     wc = mean(w)
35     f1 = f - fc
36     w1 = w - wc
37     F = fft.fft2(f1)
38     W = fft.fft2(w1)
39     x = iapconv(F,W)
40     w1 = w1[::-1,::-1]
41     F1 = fft.fft2(f1)**2
42     W1 = fft.fft2(w1)**2
43     y = iapconv(F1,W1)
44     g = fft.ifft2(x)/sqrt(fft.ifft2(y))
45 
46     return g
47 
48 f = iaread('cameraman.pgm')
49 f = asarray(f).astype(float)
50 adshow(f, title='input image')
51 w = f[25:25+17,106:106+17]
52 adshow(w, title='mask')
53 t1=time.time()
54 f1 = index_correlation(f,w)
55 t2=time.time()
56 t = t2 - t1
57 print "Time index correlation with convolution:", t
58 v, pos = max(ravel(f1)), argmax(ravel(f1))
59 (row, col) = iaind2sub(f1.shape, pos)
60 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
61 adshow(ianormalize(f1,[0,255]), title='index correlation')
62 
63 t1=time.time()
64 f2 = index_correlation_fd(f,w)
65 t = t2 - t1
66 print "Time index correlation in the frecuency domain:", abs(t)
67 v, pos = max(ravel(f2)), argmax(ravel(f2))
68 (row, col) = iaind2sub(f2.shape, pos)
69 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
70 adshow(ianormalize(log(abs(f2)+1), [0,255]),title="index correlation in the frequency domain")
71 
72 t1=time.time()
73 f3 = index_correlation_fd_fft(f,w)
74 t = t2 - t1
75 print "Time index correlation in the frecuency domain:", abs(t)
76 v, pos = max(ravel(f3)), argmax(ravel(f3))
77 (row, col) = iaind2sub(f3.shape, pos)
78 print 'found best match at (%3.0f,%3.0f)\n' %(col,row)
79 adshow(ianormalize(log(abs(f3)+1), [0,255]),title="index correlation in the frequency domain with fft")
Time index correlation with convolution: 0.425409078598
found best match at (114, 33)

Time index correlation in the frecuency domain: 0.0440819263458
found best match at (129,249)

Time index correlation in the frecuency domain: 1.88814091682
found best match at (129,249)

input image

mask

index correlation

index correlation in the frequency domain

index correlation in the frequency domain with fft

Function iavarfilter

Synopse

Variance filter.

  • g = iavarfilter(f, h)
    • g: Image.
    • f: Image. input image.
    • h: Image. kernel.

Description

Computes the variance on the neighborhood of the pixel. The neighborhood is given by the set marked by the kernel elements.

 1 from numpy import *
 2 
 3 def iavarfilter(f, h):
 4     from iapconv import iapconv
 5 
 6     f = asarray(f).astype(float64)
 7     f = f + 1e-320*(f == 0) # change zero by a very small number (prevent 'math range error')
 8     n = sum(ravel(h))
 9     fm = iapconv(f, h) / n
10     f2m = iapconv(f*f, h) / n
11     g = sqrt(f2m - (fm*fm)) / fm
12     return g

Explanation

The variance is a measure of how far each value in the data set is from the mean. Here is how it is defined:

  1. Subtract the mean from each value in the data. This gives you a measure of the distance of each value from the mean.
  2. Square each of these distances (so that they are all positive values), and add all of the squares together.
  3. Divide the sum of the squares by the number of values in the data set.

The output image contains the calculated variance at each input pixel location. Variance calculations can be restricted by setting an image mask.

 1 from ia636 import *
 2 
 3 def iavarfilter(f, h):
 4     f = asarray(f).astype(float64)
 5     f = f + 1e-320*(f == 0) # change zero by a very small number (prevent 'math range error')
 6     n = sum(ravel(h))
 7     print "n=",n
 8     fm = iapconv(f, h) / n
 9     adshow(fm,title='fm')
10     f2m = iapconv(f*f, h) / n
11     adshow(f2m,title='f2m')
12     g = sqrt(f2m - (fm*fm)) / fm
13     return g
14 
15 f = adread('cameraman.pgm')
16 adshow(f,title='input')
17 g = iavarfilter(f, [[0,1,0],[1,1,1],[0,1,0]])
18 adshow(ianormalize(g),title='output')
n= 5

input

fm

f2m

output

Example 1

1 from ia636 import iavarfilter, ianormalize
2 
3 f = adread('cameraman.pgm')
4 adshow(f)
5 g = iavarfilter(f, [[0,1,0],[1,1,1],[0,1,0]])
6 adshow(ianormalize(g))

Example 2

1 from ia636 import iavarfilter, ianormalize
2 
3 f = adread('cameraman.pgm')
4 adshow(f)
5 g = iavarfilter(f,[[1,1,1],[1,1,1],[1,1,1]])
6 adshow(ianormalize(g))

Example 3

1 from ia636 import iavarfilter, ianormalize
2 
3 f = adread('cameraman.pgm')
4 adshow(f)
5 g = iavarfilter(f,[[0,1,0],[0,1,0],[0,1,0]])
6 adshow(ianormalize(g))

Example 4

1 from ia636 import iavarfilter, ianormalize
2 
3 f = adread('cameraman.pgm')
4 adshow(f)
5 g = iavarfilter(f,[[0,0,1],[0,0,1],[0,0,1]])
6 adshow(ianormalize(g))