FPD¶
A python package for fast pixelated detector data storage, analysis and visualisation.
Web: https://gitlab.com/fpdpy/fpd/
HTML documentation (in progress):
- https://gitlab.com/fpdpy/fpd/builds/artifacts/master/file/pages_development/index.html?job=pages_development (development)
- https://fpdpy.gitlab.io/fpd/ (release)
This notebook:
Notebook¶
This notebook outlines the main functionality of the fpd package, but does not cover every detail. All functions and classes are documented in the code, some with examples not included here.
Code documentation may be read in the notebook by running in a code cell the command followed by a ?, or by depressing shift + tab when in the parenthesis of a class or function. Pressing tab multiple times while shift is held brings up successively more detail.
# print python version
import sys
print(sys.version)
import fpd
import fpd.fpd_processing as fpdp
import fpd.fpd_file as fpdf
import scipy as sp
import numpy as np
import matplotlib.pylab as plt
plt.ion()
print('fpd package version:', fpd.__version__)
BLAS libraries¶
If importing fpd prints a message about BLAS libraries, then you might find improvements in multiprocessing performace by following the advice. However, it is not nessecary in order for everything in the package to work.
On linux with OpenBLAS, the number of threads can be set by running export OMP_NUM_THREADS=1 in the shell before running python (or more accurately, before importing numpy). 1 can be changed to 0 to return to the default behaviour.
The library check function is in the cell below and may be called at any time.
fpdp._check_libs()
A more thorough check of if BLAS is multithreaded can be performed with:
fpdp.cpu_thread_lib_check()
Interactive notebook plotting¶
The nbagg backend allows interactive plotting but is slow and can't be switched mid-notebook.
#%matplotlib nbagg
so, we will use qt5 (or qt4 for older systems) and instead run the following 'magic' to switch to interactive plotting
%matplotlib qt
and switch to static inline images for documenting results
%matplotlib inline
ransac_tools¶
RANSAC fits are implemented for 1-D data and for images.
https://en.wikipedia.org/wiki/Random_sample_consensus
These allow fitting to data with outliers. They are useful for large number of applications and particularly useful in processing DPC images where descan correction needs to be applied. In this context they can be used to correct for bright field disc shifts due do imperfections in the setup of microscope, while ignoring magnetic contrast.
from fpd.ransac_tools import ransac_1D_fit, ransac_im_fit
1-D data¶
Fitting to 1-D data with outliers using regular approaches produces poor fits.
x = np.linspace(0, np.pi*2, 50)
y = 2*x + 4
# add noise
y += (np.random.random(50)-1) * 0.5
# add outliers
bad_inds = np.arange(50)%3.0==0
y[bad_inds] = np.random.random(bad_inds.sum()) + 1
# do regular fit
p = np.polyfit(x, y, 1)
yfit = np.polyval(p, x)
# plot
plt.plot(x, y, label='Data')
plt.plot(x, yfit, label='Fit')
plt.legend()
The fit can be improved with RANSAC, which clasifies the data into inliers and outliers:
# Ransac fit
fit, inliers, model = ransac_1D_fit(x, y, mode=1, residual_threshold=0.3, plot=True)
The mode parameter controls the fit function and can be:
- linear
- quadratic
- cubic
- smoothing spline
- or any user supplied function.
For example, for a cubic polynomial:
x = np.linspace(0, np.pi*2, 50)
x -= x.mean()
y = 0.5*x**3 - x**2 + 0.2*x + 1
# add outliers
bad_inds = np.arange(50)%4.0==0
y[bad_inds] = np.random.random(bad_inds.sum()) + 1
# plot data
plt.plot(x, y)
fit, inliers, model = ransac_1D_fit(x, y, mode=3, residual_threshold=0.3, plot=True)
Using the same data, we can fit with our own function, defined below. Through the param_dict dictionary, we can set sigma values to weight the fit.
def f(p, x):
return p[0]*x**3 + p[1]*x**2 + p[2]*x + p[3]
p0 = (1,1,1,1)
param_dict = {'sigma' : x*0.1+1}
fit, inliers, model = ransac_1D_fit(x, y, mode=0, model_f=f, p0=p0,
param_dict=param_dict, plot=True)
Image data¶
These are 2-D regular arrays. The cell below generates a parabolic surface with a plane backgound. The background could be due to descan in a DPC signal, while the parabolic surface is magnetic contrast.
# generate data
yy, xx = np.indices((128,)*2)
x0 = xx.mean() + 32
y0 = yy.mean() + 32
im_plane = (xx*2 + yy*1 -9.1)
im_para = 0.16*((xx-x0)**2 + (yy-y0)**2) -200
# combination of the two
im_pp = np.minimum(im_plane, im_para)
plt.matshow(im_pp)
The entire image can be fitted to without RANSAC by setting max_trials=1, min_samples=1.0, andresidual_threshold=x, where x is a suitably large value.
This will produce a poor result:
fit, inliers, _ = ransac_im_fit(im_pp, mode=1, min_samples=1.0, max_trials=1, residual_threshold=10000, plot=True)
plt.matshow(im_pp - fit)
plt.colorbar()
Enabling RANSAC gives a much better fit:
fit, inliers, _ = ransac_im_fit(im_pp, mode=1, max_trials=2000, plot=True)
plt.matshow(im_pp - fit)
plt.colorbar()
The mode parameter controls the fit function and can be:
- plane
- quadratic
- concave paraboloid with offset.
- smoothing spline
- or a user supplied function.
In all RANSAC fittings, the residual_threshold parameter must be specified carefully and the likelyhood of the fit being optimum is improved be increasing max_trials.
Masks may be supplied with the mask parameter to select the region to which to fit, which can increase the likelyhood of a good fit.
Note that multiple images may be processed at once. The axes parameter defines the image axes and by default are the last two axes.
fft_tools¶
Images can be filtered in the frequency domain interactively
%matplotlib qt
The slider and buttons control the plot setting and the mask applied to the Fourier space data.
py, px = 32.0, 16.0
y, x = np.indices((256,)*2)
im = np.zeros(y.shape, dtype=float)
im += np.sin(y/py*2*np.pi)
im += np.sin(x/px*2*np.pi)
b = fpd.fft_tools.BandpassPlotter(im, fminmax=(0.045, 0.08), gy=1, fact=1, mode='vert')
The mode can be one of ['circ', 'horiz', 'vert'] and controls the mask shape.
%matplotlib inline
The results can be accessed from the class:
plt.matshow(b.im_pass)
Or images can be filtered without the GUI.
This also alows custom masks to be used through the mask parameter (not shown here)
rtn = fpd.fft_tools.bandpass(im, fminmax=(0.045, 0.08), gy=1, mask=None, full_out=True, mode='circ')
im_pass, (im_fft, mask, extent), (rrf, yf, xf, imf) = rtn
plt.matshow(im_fft, norm=plt.mpl.colors.LogNorm())
plt.colorbar()
plt.matshow(mask)
plt.matshow(im_pass)
The radial distribution can also be calculated from the results:
r, rdf = fpd.fft_tools.fft2rdf(rrf, im_fft)
# plot distribution of angles
plt.figure()
plt.plot(r, rdf)
# plot periods of sinusoids in source data
peak_f = 1.0/np.array([px, py])
peak_amp = np.interp(peak_f, r, rdf)
plt.plot(peak_f, peak_amp, 'or')
plt.xlim(None, 0.15)
The non-gui function accepts an array of images and will process them all together. The image axes are taken as the last 2 dimensions.
DPC_Explorer¶
2-D vector data such as that from DPC analyses can be plotted interactively with the DPC_Explorer class.
Normally, the data passed to the class will be the array or similar object containing the y and x data in the first dimension. However, if an integer is passed instrad of an array, synthetic radial data is generated. If positive (negative), it is noise free (has Poissonian noise).
%matplotlib qt
b = fpd.DPC_Explorer(-64)
Again, the sliders and buttons control settings. In addition, some settings are controlled via the plots and figures. These are described in the docstring and include:
On the histogram window:
- right / left arrows for sequence navigation.
- Ctrl+ mouse drag on circle is rotate angle.
- Mouse drag on circle sets r_max.
- Shift+ mouse drag on circle sets r_min.
- Mouse drag on square moves centre.
- Ctrl+arrows for Y descan.
- Alt+arrows for X descan.
The decan correction increment is fraction of r_max.
On yx plots:
- Click+drag in y for histogram region select.
- Open Y (or X) opens the processed data in Gwyddion.
- 'lim<-r' links the scale of the y and x component plots to the radius set in the histgram (the centre is always linked).
Small changes in slider positions can be made by scrolling with the mouse wheel.
Median line correction is sometimes useful for correcting applifier noise. These may be used separately from the utils module.
Gwyddion¶
If the data requires further processing or analysis, clicking open Y or open X in the yx-plots will open the file in Gwyddion, a multipurpose data analysis program principly designed for SPM analysis.
Files may be written using:
filename = fpd.gwy.writeGSF(filename, data, open=False)
where open controls whether the image is automatically opened.
and read with:
data, dsf = fpd.gwy.readGSF(filename)
Image sequences¶
When supplied with a 4-D dataset, the 1st dimension is taken as multiple yx images.
Navigation is by clicking the prev or next buttons, or by mouse wheel scroll over them, or by left and right arrow keys.
Below, we create a sequence of 4 images of random data, and plot them.
# create random data
yx = np.random.rand(4, 2, 64, 64)
yx -= yx.mean((2,3))[..., None, None]
# stager the displacements in y
yx[-2,0] += 0.5
yx[-1,0] += 1.0
# plot
b = fpd.DPC_Explorer(d=yx)
For image sequences, all processed images can be returned using:
tr_ims = b.get_image_sequence(im_types='tr_im', images=False)
print(tr_ims.shape)
plt.imshow(tr_ims[0])