The goal of this log is to show the API of the spatial correlations and the options available. With this code, it is possible to run the spatial correlations on masked and unmasked data.
Also, it is possible to apply a correction, called symmetric averaging, which is a derivative of a method by Schatzel (1988):
Schätzel, Klaus, Martin Drewel, and Sven Stimac. "Photon correlation measurements at large lag times: improving statistical accuracy." Journal of Modern Optics 35.4 (1988): 711-718.
Technique adapted to arbitrary masks by Julien Lhermitte, Jan 2017
The correlation function in 1 dimension is: $$C = \frac{1}{N(k)} \sum \limits_j^{N_t} I_j I_{j+k} M_j M_{j+k}$$ We may normalize it by its average intensity in two different ways:
the normalized correlation function is just divided by the average squared: $$cc_{reg} = \frac{CC}{\bar{I}^2}$$ where: $$\bar{I} = \frac{1}{N(k)} \sum \limits_{j=1}^{N_t} I_j$$ is average intensity and where $$N(k) = \sum \limits_{j= 1}^{N_t}M_j M_{j+k}$$ (Note that in the limit of no mask, $N(k) = N_t$ as it should, mask has effect of inducing a $k$ dependence on the effective ''$N_t$'')
For symmetric averaging, we define two new averages, $I_p$ and $I_f$ (I 'past' and I 'future'): $$I_p = \frac{1}{N(k)} \sum \limits_j I_j M_j M_{j+k}$$ $$I_f = \frac{1}{N(k)} \sum \limits_l I_{l+k} M_l M_{l+k}$$ we define symmetric averaging as: $$cc_{sym} = \frac{CC}{\bar{I}_p \bar{I}_p}$$
Schatzel shows this averaging is superior for the case of a simple ''mask'' : a 1D time series (data outside of time range is ''masked'')
In [1]:
%matplotlib inline
import numpy as np
#from pyCXD.tools.CrossCorrelator import CrossCorrelator
from skbeam.core.correlation import CrossCorrelator
import matplotlib.pyplot as plt
from skbeam.core.roi import ring_edges, segmented_rings
# for some convolutions, used to smooth images (make spatially correlated images)
# avoid more dependencies for this example
def convol2d(a,b=None,axes=(-2,-1)):
''' convolve a and b along axes axes
if axes 1 element, then convolves along that dimension
only works with dimensions 1 or 2 (1 or 2 axes)
'''
from numpy.fft import fft2, ifft2
if(b is None):
b = a
return ifft2(fft2(a,axes=axes)*np.conj(fft2(b,axes=axes)),axes=axes).real
def pos2extent(pos):
# convenience routine to turn positions to extent
# left right bottom top. For 2D data
extent = [pos[1][0], pos[1][-1], pos[0][-1], pos[0][0]]
return extent
In [2]:
# test 1D data
sigma = .1
Npoints = 1000
x = np.linspace(-10, 10, Npoints)
y = convol2d(np.random.random(Npoints)*10, np.exp(-x**2/(2*sigma**2)),axes=(-1,))
mask_1D = np.ones_like(y)
mask_1D[10:20] = 0
mask_1D[60:90] = 0
mask_1D[111:137] = 0
mask_1D[211:237] = 0
mask_1D[411:537] = 0
mask_1D *= mask_1D[::-1]
y_masked = y*mask_1D
cc1D = CrossCorrelator(mask_1D.shape)
cc1D_symavg = CrossCorrelator(mask_1D.shape,normalization='symavg')
cc1D_masked = CrossCorrelator(mask_1D.shape,mask=mask_1D)
cc1D_masked_symavg = CrossCorrelator(mask_1D.shape, mask=mask_1D,normalization='symavg')
ycorr_1D = cc1D(y)
ycorr_1D_masked = cc1D_masked(y*mask_1D)
ycorr_1D_symavg = cc1D_symavg(y)
ycorr_1D_masked_symavg = cc1D_masked_symavg(y*mask_1D)
# the x axis
ycorr_1D_x = cc1D.positions
ycorr_1D_masked_x = cc1D_masked.positions
ycorr_1D_symavg_x = cc1D_symavg.positions
ycorr_1D_masked_symavg_x = cc1D_masked_symavg.positions
In [3]:
ycorr_1D[0].shape
Out[3]:
In [4]:
# plot 1D Data
plt.figure(0);plt.clf();
plt.plot(x,y)
plt.plot(x,y*mask_1D)
plt.xlabel("position")
plt.ylabel("intensity (arb. units)")
Out[4]:
In [50]:
plt.figure(1);plt.clf();
plt.plot(ycorr_1D_x, ycorr_1D,color='k',label='regular')
plt.plot(ycorr_1D_masked_x, ycorr_1D_masked,color='r',label='masked')
plt.plot(ycorr_1D_symavg_x, ycorr_1D_symavg,color='g',label='symavg')
plt.plot(ycorr_1D_masked_symavg_x, ycorr_1D_masked_symavg,color='b',label='masked + symavg')
plt.ylim(.9,1.2)
plt.xlabel("shift ($\Delta x$)")
plt.ylabel("Correlation")
plt.legend()
Out[50]:
In [23]:
# test 2D data
Npoints2 = 100
x2 = np.linspace(-10, 10, Npoints2)
X, Y = np.meshgrid(x2,x2)
Z = np.random.random((Npoints2,Npoints2))
Z = convol2d(Z, np.exp(-(X**2 + Y**2)/2./sigma**2))
mask_2D = np.ones_like(Z)
mask_2D[10:20, 10:20] = 0
mask_2D[73:91, 45:67] = 0
mask_2D[1:20, 90:] = 0
cc2D = CrossCorrelator(mask_2D.shape)
cc2D_symavg = CrossCorrelator(mask_2D.shape,normalization='symavg')
cc2D_masked = CrossCorrelator(mask_2D.shape,mask=mask_2D)
cc2D_masked_symavg = CrossCorrelator(mask_2D.shape, mask=mask_2D,normalization='symavg')
ycorr_2D = cc2D(Z)
ycorr_2D_masked = cc2D_masked(Z*mask_2D)
ycorr_2D_symavg = cc2D_symavg(Z)
ycorr_2D_masked_symavg = cc2D_masked_symavg(Z*mask_2D)
ycorr_2D_pos = cc2D.positions
ycorr_2D_masked_pos = cc2D_masked.positions
ycorr_2D_symavg_pos = cc2D_symavg.positions
ycorr_2D_masked_symavg_pos = cc2D_masked_symavg.positions
In [24]:
plt.figure(2);plt.clf();
plt.subplot(2,2,1)
plt.title("not masked")
plt.imshow(Z)
plt.subplot(2,2,2)
plt.title("masked")
plt.imshow(Z*mask_2D)
Out[24]:
In [38]:
vmin=.95; vmax=1.03
plt.figure(3);plt.clf();
plt.subplot(2,2,1)
plt.title("regular")
plt.imshow(ycorr_2D,extent = pos2extent(ycorr_2D_pos))
#plt.axhline(ycorr_2D_masked.shape[0]//2)
plt.clim(vmin,vmax)
plt.xlim(-30,30)
plt.ylim(-30,30)
plt.subplot(2,2,2)
plt.title("masked")
plt.imshow(ycorr_2D_masked, extent = pos2extent(ycorr_2D_masked_pos))
#plt.axhline(ycorr_2D_masked.shape[0]//2)
plt.clim(vmin,vmax)
plt.xlim(-30,30)
plt.ylim(-30,30)
plt.subplot(2,2,3)
plt.title("symavg")
plt.imshow(ycorr_2D_symavg, extent = pos2extent(ycorr_2D_symavg_pos))
#plt.axhline(ycorr_2D_masked.shape[0]//2)
plt.clim(vmin,vmax)
plt.xlim(-30,30)
plt.ylim(-30,30)
plt.subplot(2,2,4)
plt.title("mask + symavg")
plt.imshow(ycorr_2D_masked_symavg, extent = pos2extent(ycorr_2D_masked_symavg_pos))
#plt.axhline(ycorr_2D_masked.shape[0]//2)
plt.clim(vmin,vmax)
plt.xlim(-30,30)
plt.ylim(-30,30)
Out[38]:
In [49]:
plt.figure(4);plt.clf();
plt.plot(cc2D.positions[1], ycorr_2D[cc2D.centers[0]],label="reg")
plt.plot(cc2D_masked.positions[1], ycorr_2D_masked[cc2D_masked.centers[0]],label="masked")
plt.plot(cc2D_symavg.positions[1], ycorr_2D_symavg[cc2D_symavg.centers[0]],label="symavg")
plt.plot(cc2D_masked_symavg.positions[1], ycorr_2D_masked_symavg[cc2D_masked_symavg.centers[0]],label="masked+symavg")
plt.ylim(0.8, 1.2)
plt.xlabel("shift ($\Delta x$)")
plt.ylabel("Correlation")
plt.legend()
Out[49]:
In [10]:
# make id numbers
edges = ring_edges(1, 20, num_rings=2)
segments = 5
x0, y0 = np.array(mask_2D.shape)//2
maskids = segmented_rings(edges,segments,(y0,x0),mask_2D.shape)
cc2D_ids = CrossCorrelator(mask_2D.shape, mask=maskids)
cc2D_ids_symavg = CrossCorrelator(mask_2D.shape,mask=maskids,normalization='symavg')
ycorr_ids_2D = cc2D_ids(Z)
ycorr_ids_2D_symavg = cc2D_ids_symavg(Z)
In [11]:
plt.figure(2);plt.clf();
plt.imshow(maskids)
Out[11]:
Here, we see that without symmetric averaging, the correlations quickly come back at values higher than the point of initial correlation, whereas with symmetric averaging, the result looks more as what is expected, a nice Gaussian like curve centered in image. (Center of image is zero correlation)
In [52]:
vmin=.95; vmax=1.1
fig, axes = plt.subplots(2,4)
ax1 = axes[:len(axes)//2].ravel()
ax2 = axes[len(axes)//2:].ravel()
for i in range(len(ax1)):
plt.sca(ax1[i])
plt.title("regular")
plt.imshow(ycorr_ids_2D[i],extent=pos2extent(cc2D_ids.positions[i]))
plt.clim(vmin,vmax)
plt.sca(ax2[i])
plt.title("sym avg")
plt.imshow(ycorr_ids_2D_symavg[i],extent=pos2extent(cc2D_ids_symavg.positions[i]))
plt.clim(vmin,vmax)
In [13]:
## Cross correlate image with itself shifted
Z2 = np.roll(np.roll(Z, 4,axis=0),-5,axis=1)
In [14]:
ycorr_ids_2D_shift = cc2D_ids(Z, Z2)
centers_ids_2D_shift = cc2D_ids.centers
ycorr_ids_2D_symavg_shift = cc2D_ids_symavg(Z,Z2)
centers_ids_2D_symavg_shift = cc2D_ids_symavg.centers
In [56]:
vmin=.95; vmax=1.05
fig, axes = plt.subplots(2,4)
ax1 = axes[:len(axes)//2].ravel()
ax2 = axes[len(axes)//2:].ravel()
for i in range(len(ax1)):
plt.sca(ax1[i])
plt.title("regular")
plt.imshow(ycorr_ids_2D_shift[i], extent=pos2extent(cc2D_ids.positions[i]))
yc, xc = centers_ids_2D_shift[i]
plt.axvline(xc)
plt.axhline(yc)
plt.clim(vmin,vmax)
plt.sca(ax2[i])
plt.title("sym avg")
plt.imshow(ycorr_ids_2D_symavg_shift[i], extent=pos2extent(cc2D_ids_symavg.positions[i]))
yc, xc = centers_ids_2D_symavg_shift[i]
plt.axvline(xc)
plt.axhline(yc)
plt.clim(vmin,vmax)
In [54]:
mask_test = (maskids == 1).astype(float)
from scipy.signal import fftconvolve
from numpy.fft import fft2, ifft2, fftshift
cc = fftconvolve(mask_test, mask_test, mode='same')
cc2 = fftshift(ifft2((np.conj(fft2(mask_test)))*fft2(mask_test)).real)
In [ ]: