Find a 8 in a forest of 9
The goal here is to find one non-repeting pattern in an image. First, find the pattern, then second, compute correlation.
Let's first initialize the notebook:
In [27]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import numpy as np
phi = (np.sqrt(5)+1)/2
fig_width = 10
figsize = (fig_width, fig_width/phi)
In [28]:
import imageio as io
im = io.imread('../files/2021-12-01-99999999999998.jpg')
im = im.sum(axis=-1)
In [29]:
N_X, N_Y = im.shape
N_X, N_Y
Out[29]:
Original image
In [30]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(im, cmap=plt.gray());
Zoom
In [31]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(im[:40, :40]);
In [32]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(im[-40:, -40:]);
Crop
In [33]:
im = im[:, 24:1194]
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(im, cmap=plt.gray());
Average on one axis
In [34]:
im_x = im.mean(axis=1)
im_x = np.roll(im_x, N_X//2) # avoid border effects
#im_x = np.correlate(im_x, [-1, 2, 1], 'same')*1. # contrast detection
im_x -= im_x.mean()
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(im_x[:100]);
Cross-correlation
In [35]:
xcorr = np.correlate(im_x, im_x, 'same')*1.
xcorr /= xcorr.max()
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(xcorr);
In [36]:
N_X//2, im_x.shape, xcorr.shape, im.shape
Out[36]:
In [37]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(xcorr[(N_X//2+1):(50+N_X//2)]);
In [38]:
period_X = np.argmax(xcorr[(1+N_X//2):]) + 1
In [39]:
print(f'{period_X=}')
The other dimension:
In [40]:
im_y = im.mean(axis=0)
im_y -= im_y.mean()
#im_y[-1] = im_y[0]
#im_y[-2] = im_y[1]
im_y = np.roll(im_y, N_Y//2) # avoid border effects
#im_y = np.correlate(im_y, [-1, 2, -1], 'same')*1. # contrast detection
In [41]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(im_y);
In [42]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(im_y[:40]);
In [43]:
xcorr = np.correlate(im_y, im_y, 'same')*1.
In [44]:
xcorr /= xcorr.max()
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(xcorr[(N_Y//2+3):]);
In [45]:
xcorr /= xcorr.max()
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**2))
ax.plot(xcorr[(N_Y//2+3):(40+N_Y//2)]);
In [46]:
period_Y = np.argmax(xcorr[(3+N_Y//2):]) + 3
In [47]:
print(f'{period_Y=}')
One sample (other random choices would most probably fit - do not go in the borders!)
In [48]:
idx, idy = 0, 26
kernel = im[idx:idx+period_X, idy:idy+period_Y]
In [49]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(kernel);
compute cross correlation using https://laurentperrinet.github.io//sciblog/posts/2017-09-20-the-fastest-2d-convolution-in-the-world.html
In [50]:
from numpy.fft import fft2, ifft2
xcorr = np.real(ifft2(fft2(im)*fft2(kernel, s=im.shape)))
Now you can spot where there is a difference (and then in the rest due to kerning...)
In [51]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.matshow(xcorr, cmap=plt.magma());
some book keeping for the notebook¶
In [52]:
%load_ext watermark
%watermark -i -h -m -v -p numpy,matplotlib,imageio -r -g -b