cross-peak-DH

corr-peak-DH.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

import abel
from scipy.ndimage import zoom
from scipy import signal

import peakutils

fn = "O-10N21024.txt"   # raw VMI  binned to 1024x1024 pixel image
im = np.loadtxt(fn)
im = zoom(im, 2)  # improves angular resolution for small radius
im /= im.max()  # normalize, consistent between images

# centre
imc = abel.tools.center.center_image(im, center='convolution')

# map to polar coordinates
polarim, r, t = abel.tools.polar.reproject_image_into_polar(imc, Jacobian=False)

R = r[:, 0]  # slice radial grid

nslices = 32  # divide image into 32 slices

pie = np.array_split(polarim, nslices, axis=1) 

angle = np.array(np.hsplit(t[0], nslices)).mean(axis=1)
n2 = len(angle)//2   #  mid point is zero angle

# restrict radial range to peak of interest
subr = np.logical_and(R >= 0, R <= 1000)

previous = pie[0][subr].sum(axis=1) # -pi/2 slice intensity profile

xcorrpk = []
profpk = []

for ang, aslice in zip(angle[1:], pie[1:]):
    # slice profile
    profile = aslice[subr].sum(axis=1)

    # cross correlation
    corr = signal.correlate(previous, profile, "same")
    plt.plot(corr)

    # find xcorr position
    ind = peakutils.indexes(corr, thres=0.3, min_dist=40)
    
    # fine the peak location, setting the center of the R-array to zero
    # I'm not quite sure why the "-0.5" is needed. It might be something to do with the pixels?
    peaks = peakutils.interpolate( R[subr]-np.max(R[subr])*0.5-0.5, corr, ind=ind)

    # sort the peaks by smallest absolute value, to find the one in the center
    peaks = peaks[np.argsort(np.abs(peaks))]
    print peaks
    peak = peaks[0]
    xcorrpk.append(peak)

    # peak position using raw profile
    indraw = peakutils.indexes(profile, thres=0.3, min_dist=40)
    peaks = peakutils.interpolate(R[subr], profile, ind=indraw)
    profpk.append(peaks[0]) 

    previous = profile



plt.figure()

xcorrpk = np.array(xcorrpk)
# --- plots --------------
gs = gridspec.GridSpec(1, 1)
ax0 = plt.subplot(gs[0, 0])

ax0.plot(angle[1:], np.cumsum(-xcorrpk)+np.max(R[subr])*0.5, 'o-', label='xcorr')
ax0.plot(angle[1:],  profpk, 'o-', label='peak')

ax0.axis(xmin=-np.pi, xmax=np.pi)
ax0.set_xticks([-np.pi, 0, np.pi])
ax0.set_xticklabels([r"-$\pi$", "0", r"$\pi$"])
ax0.set_xlabel("angle (radians)")
ax0.set_ylabel("relative position")

plt.legend(frameon=False)
plt.savefig("correlation.png", dpi=75)
plt.show()