In [ ]:
# let's make a test for subpixel localization
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from openpiv.pyprocess import find_first_peak
In [3]:
N = 64

corr = np.zeros((N,N))

corr[2:5,2:5] = 1
corr[3,3] = 2
corr[3,4] = 3
corr[3,5] = 1
corr
Out[3]:
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 1., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
In [4]:
pos,height = find_first_peak(corr)
In [5]:
pos,height
Out[5]:
((3, 4), 3.0)
In [6]:
from openpiv.pyprocess import find_subpixel_peak_position
In [7]:
find_subpixel_peak_position(corr)
Out[7]:
(3.0, 3.769577293545741)

let's find some corner cases

In [8]:
# peak on the border 
corr = np.zeros((N,N))

corr[:3,:3] = 1
corr[0,0] = 2
corr[0,2] = 3
corr[0,3] = 1
corr
Out[8]:
array([[2., 1., 3., ..., 0., 0., 0.],
       [1., 1., 1., ..., 0., 0., 0.],
       [1., 1., 1., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

Corner case 1: peak on the border

it is disregarded in our function because we cannot define well the subpixel position. Or do we?

In [9]:
find_subpixel_peak_position(corr)
Out[9]:
(nan, nan)
In [10]:
# peak on the border 
corr = np.flipud(corr)
corr
Out[10]:
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [1., 1., 1., ..., 0., 0., 0.],
       [1., 1., 1., ..., 0., 0., 0.],
       [2., 1., 3., ..., 0., 0., 0.]])
In [11]:
find_subpixel_peak_position(corr)
Out[11]:
(nan, nan)
In [12]:
corr = np.fliplr(corr)
corr[-2,-1]=5
corr
Out[12]:
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 1., 1., 1.],
       [0., 0., 0., ..., 1., 1., 5.],
       [0., 0., 0., ..., 3., 1., 2.]])
In [13]:
find_subpixel_peak_position(corr)
Out[13]:
(nan, nan)
In [14]:
## Corner case 2: zero next to the peak - the log(0) fails
In [15]:
corr = np.zeros((N,N))

corr[2:5,2:5] = 1
corr[3,3] = 2
corr[3,4] = 3
# corr[3,5] = 1
corr
Out[15]:
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 1., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
In [16]:
find_subpixel_peak_position(corr)
Out[16]:
(3.0, 3.5230088020336483)
In [17]:
fig, ax = plt.subplots(figsize=(12,12))
ax.pcolor(corr[:8,:8])
for eps in np.logspace(-15,5):
    i,j = find_subpixel_peak_position(corr+eps)
    # print(i,j)
    ax.plot(i,j,'rx')
In [18]:
import matplotlib.pyplot as plt
In [19]:
plt.pcolor(corr[:8,:8])
plt.plot(i,j,'ro')
Out[19]:
[<matplotlib.lines.Line2D at 0x7ff38e6d5e20>]
In [20]:
for method in ['gaussian','parabolic','centroid']:
    i,j = find_subpixel_peak_position(corr,method)
    print(i,j)
    i,j = find_subpixel_peak_position(corr+eps,method)
    print(i,j)
3.0 3.5239508804084627
3.0 3.7499981250225445
3.0 3.75
3.0 3.75
3.0 3.600000095999977
3.0 3.9999933334444426
In [ ]: