scikit-beam/scikit-beam

Data analysis tools for X-Ray, Neutron and Electron sciences

MAILING LIST

scikit-beam

Documentation

Examples

scikit-beam-examples repository

• Powder calibration (still needs tilt correction)
• 1-time correlation
  • dir
  • Jupyter notebook
• Differential Phase Contrast
  • dir
  • Jupyter notebook
• Fast conversion to reciprocal space
• X-Ray Speckle Visibility Spectroscopy
• Basic Plotting of X-Ray Fluorescence Elemental Lines

Quick start

install with conda

conda install scikit-beam -c nsls2forge

install development version with setuptools

git clone [email protected]:scikit-beam/scikit-beam.git
cd scikit-beam
python setup.py install

set up for development

git clone [email protected]:scikit-beam/scikit-beam.git
cd scikit-beam
python setup.py develop
pip install pytest coverage setuptools

make sure all the tests pass!

python run_tests.py

and you can check the code coverage with

coverage run run_tests.py
coverage report -m

scikit-beam's Issues

Fast conversion to Q for other diffractometers

Description

Stuart's advice is to write a c function (or set of functions) for each class of diffractometer. The classes are:

• 6 circle (4 circle is a subclass of this) (CHX/CSX)
• 4 + 2 circle (4 on the sample, 2 on the detector)
• 5 circle

3-D Reconstruction via binning

Description

Problems to be solved with this issue

• CSX/CHX: Reconstruct reciprocal space from a 4 x N list of (H, K, L, i)

Existing code related to this issue

pyspec: https://github.com/stuwilkins/pyspec

Look at pyspec/src/ctrans.c

Line 310: static PyObject* gridder_3D(PyObject *self, PyObject *args, PyObject *kwargs){

and Line 375:

unsigned long c_grid3d(double *dout, unsigned long *nout, double *standarderror, 
               double *data, double *grid_start, double *grid_stop, int max_data, 
               int *n_grid, int norm_data){

sort out relative vs absolute imports

determine if best practice is to use

from nsls2.foo import bar

or

from .foo import bar

Image Construction

Description

This use case is primarily related to SRX. IIRC, one of the data collection modes on SRX is to raster the sample through the beam, collecting an MCA spectrum for each sample position. This MCA spectrum at each position is then used to reconstruct a 2-D image of the sample. @tacaswell can probably make this description significantly more verbose

Problems that need to be solved for this Issue

• SRX Approval

Pool of physics functions

It is convenient if we can collect more physics functions for fitting. This will speed up our work to cover more beamlines.

I wrote some functions for fluorescence fitting, as shown in https://github.com/licode/NSLS2/blob/xrf_fit/nsls2/fitting/model/physics_peak.py (PR #40)

Also some can be found in lmfit
https://github.com/licode/lmfit-py/blob/master/lmfit/utilfuncs.py

To be done:

• collect functions or models from beamline scientists or other resources
• modularize all the physics functions in an efficient way

[TAC edited to add githubfoo]

Convert MCA Channel to Energy

Description

X-ray fluorescence is measured by exciting a sample with some incident energy and then measuring the energy emitted by the sample as it's excited state decays. This emitted energy is measured with a Multi-Channel Analyzer (MCA). SRX has 3 MCA detectors with 4096 channels. To obtain useful information from these MCA measurements, the channel needs to be converted to an energy, presumably through some sort of known standard.

+----------------+              +------+
|Fluorescence MCA|  ~~~~~~~~~~~ |sample|
+----------------+              +------+
                                   ^
                                   |
                                   |
                                   X
                                   R
                                   A
                                   Y
                                   S
                                   |
                                   |
                                   |

Problems to be solved

• Successfully calibrate HXN MCAs
• Successfully calibrate SRX MCAs

And either:

• Incorporate existing calibration code into nsls2.calib

or

• Measure a bunch of known samples and compare the peak position with the expected peak position. FIt a curve to be able to convert bins to energy

Automated Detector Calibration

We need to have a way to automatically calibrate the detector distance, wavelength, and detector rotations. Currently this is done by a user in either pyFAI or Fit2D. The user selects rings on the diffraction pattern of the standard. The locations of these rings are compared to a file which contains the d-spacings for the calibrant, which is then used to refine the detector distance and other parameters. pyFAI has a fairly good GUI version of this, but it would be nice if we could use an algorithm to pick out the appropriate rings. This may require knowledge of the beam point of direct incidence.

CPU-based Fourier Transforms

Description

fftw is a high performance fourier transform library

fftw performs heuristics on your machine to determine the fastest way to compute 1-D (and maybe 2D?) fourier transforms specific to the hardware on your machine. It will run slowly the first time but be very snappy after that.

Problem

• Use pyFFTW

Find and fit peaks , diffraction peaks (1 D)

• Using a reference sample find and fit the peaks (1D)
• Fit peak to find the maximum

Calibrate Strip Detector

Description

Inelastic x-ray scattering is measured with a position sensitive strip detector (PSD) to detect the position of the scattered photons. These 6 PSDs need to be converted to energy, presumably through a known standard or some set of standards?

Problems that need to be solved

• Figure out if code exists to calibrate strip detectors

Either:

• Incorporate existing calibration code into nsls2.calib

or

• Write new calibration code for nsls2.calib
• Compute tilt angle for 1-D detector
• Fitting a scatter plot of energy versus position

Determining ideal bin size for Powder Diffraction

We will need a way to calculate the size of the bins for any integration algorithm we implement. If the bins are too large, then we loose resolution in our diffraction pattern. If the bins are too small, then some may contain no pixels, causing the pattern to become jagged. This will also determine our maximum Q resolution, an important part of propagating diffraction pattern error. Most likely this should be calculated from the detector setup/calibration information.

Energy <-> Wavelength converter

Description

It seems that there is some bimodal distribution of scientists who prefer wavelength (myself included) or energy when dealing with X-rays. Thus it would appear that a converter between the two would be a useful addition to this library.

Write conda recipe

pick intensity information out of `process_to_q` and `process_grid`

Mixing positional and intensity data in the same array is bad practice. It will make our lives more difficult if we decide to work with quantities later.

It is better to carry the information around in a pair of arrays, one hkl and one intensity.

Determine detector tilt angles

Description

2-D detectors are never perfectly perpendicular to the incident beam (call it z). In order to compute pixel positions, the detector rotation about the x- and y- directions need to be taken into account.

Problems to be solved

• Determine if code exists to compute detector tilt angles from a calibration standard
• Wrap that code if possible

Otherwise:

• Figure out what calibration standards are required
• Write code to compute detector tilt angles from the required calibration standards

See Also

#41

Integrate continuous benchmarking to the NSLS2 repository

Description

airspeed velocity is a tool for benchmarking Python packages over their lifetime. It is primarily is designed to benchmark a single project over its lifetime using a given suite of benchmarks. The results are displayed in an interactive web frontend that requires only a basic static webserver to host.

scipy talk on airspeed velocity

This would be a useful integration into this analysis library

Image Corrections: Dark Current Subtraction

Description

The dark current subtraction at the lowest level is actually very simple, as it is merely the measured detector image with the shutter open minus some reference detector image with the shutter closed. The trouble with this requirement is that there is no pre-determined method by which dark current images are collected during the experiment. In the worst case, the dark current images (shutter closed) will be randomly interspersed with the shutter open images.

Someone, somewhere is going to have to correlate the shutter_open images with the shutter_closed images so that they can be subtracted from each other. We can probably make a function that works by taking in a run_header nested dictionary from the data broker and spits out a series of dark current subtracted images. In fact, that might be the best approach, for now.

Convert the gridder to use multi-threading

Description

In order to make the gridder faster, convert it to use ipython parallel, OpenMP, etc...

Issues that need to be solved

• Prototype with ipython parallel
• Prototype with OpenMP

Nonlinear 1-D Fitting

Description

Provides interface to non-linear which will wrap scipy.optimize.

http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html

Dan brought up 'lmfit' which is a much improved API on scipy.curvefit

• http://cars9.uchicago.edu/software/python/lmfit/
• lmfit documentation notebook
  • http://nbviewer.ipython.org/github/danielballan/lmfit-py/blob/master/examples/lmfit-model.ipynb

Problems that need to be solved for this Issue

• CSX/CHX: Fit 1-time correlation
• IXS: Fit 6 elastic peaks
• SRX/HXN: Reduce MCA spectrum via peak fitting
• XPD: Fit diffraction peaks in 1-D

Owner search key

In the search front end GUI for userapi.commands.search(), turn it into a combobox with the possible user names based on your current authentication

Array Algebra

Description

Array algebra is mathematical operations on arrays that can be exclusively handled by numpy whenever a user is directly writing python and thus should not be included in the NSLS2 library. For doing array algebra in a GUI, we will need to wrap these numpy operations and so these wrappers belong in VisTrails layer. Thus, see NSLS-II/userpackages Issue 2 https://github.com/NSLS-II/userpackages/issues/2

Linear 1-D Fitting

Provides simple interface to least squares fitting which will wrap scipy.optimize.

http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html

Decide on data structure for image stacks

Description

It is unclear how we are going to represent stacks of images in python data structures.

Option 1: list of 2-D images (x, y)
Option 2: 3-D numpy array (time, x, y)

Atomic Constants

Description

Need easy access to reasonably standard atomic constants (mass, etc...)

Need access to atomic constants that are specific to elastic scattering (scattering factors, Cromer-Mann coefficients, etc...), inelastic scattering (has some constants), X-ray fluorescence (also has some contants, etc...

Existing Options

Existing options may not solve the entire problem, but they serve as a starting point, either to beg/borrow/incorporate or simply as education on how someone else did this

• pyspec
• xraylib

Installation Procedures

To get xraylib set up 'the easy way'

from @tacaswell:

• checkout my fork of conda-recipes and the xraylib branch (https://github.com/tacaswell/conda-recipes/tree/xraylib)
• conda build ~/other_source/conda-recipes/xraylib/

(modify this to point at the actual location your system built)
conda install -c file:///home/tcaswell/miniconda/conda-bld/ xraylib

At that point you should be able to import xraylib in python.

You will need autoreconf installed for this to work

logging

python has very good built-in logging tools, we should use them.

Messages have an importance level associated with them so at runtime you can select the verbosity level.

There is also a flexible framework for piping the messages to storage so it should not be too bad to pipe (a sub set of) them back into the data broker/exlog/olog.

2-D Feature Finding

Description

2-D images frequently have many interesting features which are of scientific interest. In reciprocal space, these are frequently diffraction peaks (2-D Voigt profile), planes of diffuse scattering (can be seen as cross sections through 2-D sheets) or rods of diffuse scattering (can be seen as cross sections through 1-D cigar). For examples of diffuse scattering, see

• High Energy X-ray Diffuse Scattering
• Neutron Diffuse Scattering

Problems to be solved

Reciprocal Space Images

Typically these are symmetric about the image center

• Find Bragg Peaks
• Find Diffuse Sheets
• Find Diffuse Rods
• Find Bragg Rings

Geometric shapes

Five sided geometries are exceedingly rare in crystallography, since pentagons cannot be tessellated

• Triangles
• Squares
• Rectangles
• Hexagon
• Rhombus
• Octagon

Real Space Images

Not necessarily symmetric about the "image" center

• Circles
• Ellipses
• Curved Lines
• Straight Lines
• Connected Regions: e.g.: Threshold the image to produce a binary image. Determine the boundaries of the connected regions that are all True or all False. Related to segmentation and (water shed)ing

Input Validation

Description

To harden our analysis library, we will need to validate input parameters for externally facing APIs.

• Type
• Range
• Size (lists/dicts/etc)
• Validate that lists (arrays, etc.) are the same length

I think that @tacaswell is advocating the use of decorators to do this

Masking

Description

Measured diffraction images contain a variety of artifacts (dead pixels, saturated pixels, beam stop arm, etc...) that need to be removed in order for data analysis to proceed. Masking covers all of those tasks.

Problems that need to be solved

• Mask dead pixels
• Mask saturated pixels
• Mask beam stop arm
• Mask "zingers" (cosmic rays)
• Interactive masking

Data structure for XRF fit

Description

Data structure needs to be well designed for x-ray fluorescence fitting. There are so many parameters, such as the ones related to all the physics peaks, instruments, and all the elements information (with xraylib). Moreover, this data structure should also have easy interface to match lmfit, because we need to pass these data into the fitting engine eventually.

Problems to be solved

• Design data structure, such as using dictionary or class object
• Develop interface to wrap xraylib or other xray database to obtain information of all the elements

GUI for fluorescence fitting

Description

Once the fitting engine for fluorescence is finished, we need to create a user-friendly GUI for beamline scientists to use. The GUI will be designed to fit into the current framework VisTrails. We will implement the code according to vistools.

Problems to be solved

• Obtain user requirements from beamline scientists from HXN, SRX and so on
• Code design and implementation based on vistools
• Code test and maintain

gracefully degrade if missing external dependencies

Find an easy to use (from the dev side) way of marking functions to not be available due to missing external dependencies.

Cake

Description

Problems to be solved for this issue

• XPD
• CHX

XPD

XPD needs to be able to select a "cake-slice" of a 2-D image relative to the image center (xc, yc) by specifying either:

• r_min, r_max

or

• r0, delta_r

and

• theta_min, theta_max

or

• theta0, delta_theta

Output
XPD's output should be a two 1-D numpy arrays. The x-axis should be the radius from the image center in terms of:

• pixels, if no additional information is provided
• two-theta, if the required information is provided (sample distance, wavelength)
• q, with the same required information as two-theta.
• Q or 2-theta return value should be an input variable flag

The y-axis should be the average/sum/stdev/median/etc. of all pixels which share that x-value

CHX

CHX needs to be able to select a series of rings with ring centers (r0) and ring thicknesses (delta_r) relative to the image center (xc, yc)

Output
CHX's output should be one number for each ring. This single number represents the sum or average of all pixels bounded by the rings

References

Image from http://oldwww.iucr.org/iucr-top/cif/cbf/PAPER/olivine.html

1-D Peak Finding

Problems that need to be solved for this Issue:

• IXS: Find elastic peak in the spectrum (they will be the biggest peak in the spectrum)
• XPD: Find the diffraction peaks in 1-D

Existing code related to this issue

Link to code (https://gist.github.com/tacaswell/3048730) that will find peaks via zero crossings:

But does not deal with noise very well or a curvy background:

@jammcc has been working with this code already for the automated masking stuff

make c-extensions python 3 compatible

This should be done to ctrans in conjunction with #23

Extract Region of Interest (ROI) from data stack

Description

There are numerous cases where region(s) of interest in an image need to be extracted or masked. The two scientific use cases that come to mind for extracting a region of interest from an image stack are:

1. Measuring the change in some statistical property of Bragg peaks as the temperature is changed.
2. Measuring crystallization kinetics with diffraction (a large portion of my graduate work). This is an in-situ experiment where the sample is melted then rapidly quenched to a temperature where crystallization is observed. As the sample crystallizes, Bragg peaks will appear and then increase in intensity. An integral of each set of pixels that correspond to a Bragg peak will result in a curve that shows the time-dependence of crystal growth which can be normalized and fit to crystallization models. Long story short, multiple (can be up to a few hundred) diffraction spots might be of interest per data stack.

This issue is related to #30, at least in regards to applying the mask to the image.

1. One approach to ROI extraction that @tacaswell mentioned is to simply turn the ROI into a binary mask and in-place multiply, i.e., numpy.multiply(arr1, arr2), the binary mask with the image data.
2. Another approach to ROI extraction is to use the pixel coordinates of the ROI to directly access and extract the elements from a data stack.

The binary mask is a constant time approach for ROI extraction and is faster than when the ROI is greater than about 9000 pixels in a 2048x2048 image. The other method is faster when the ROI is smaller than 9000 pixels. This benchmark script is online here and the output for 10 2048x2048 images (stored in a list) is shown below.

Open Questions

• Are there faster/smarter/better ways to do this ROI extraction?
• Also, this issue seems primarily directed at obtaining ROI's from 2-D image data. Should we work on a consistent ROI API for N-D data? Or is it simpler to treat 1-D, 2-D, 3-D, N-D ROI's separately?

Image Corrections: Grid Correction

Description

Basically, this: http://toothwalker.org/optics/distortion.html

IIRC, in order to determine this type of distortion you:

1. put a grid in front of the beam and look at its projection onto the detector.
2. Then magic happens
3. Then you have some sort of transformation (matrix?) which converts the raw image to an 'unwarped' image

I talked with Sanjit and Milinda to try to get more information about this from them. They were confused about what I was talking about, so I must have heard this from Eric Dooryhee. I will speak with him when he gets back next week.

project logo

Draft logos from @doraf

For use as the organization avatar:

For use in the header on the documentation website/other branding (the text can obviously be tuned):

Avizo file converter

Description

Avizo is a popular tool for quantification in tomography. @giltis has (correctly) made the point that this NSLS-II toolbox will get broader usage if a file converter exists to translate the avizo file format into python and numpy primitives.

Problems to be solved

• Correctly separate header from binary data
• Parse header into python primitives
• Parse binary data into numpy arrays
• Construct metadata dictionary whose keys match those needed by the NSLS-II toolbox
• Create appropriate test functions for the conversion functions and add basic test data to test_data folder

Image Registration

@tacaswell should most likely fill this one out.

Description

Problems to be solved

• FXI: [description]
• SRX: [description]
• HXN?

Image Corrections: Flat field correction

Description

This is my best guess as to the requirements of flat field correction:

1. Expose the detector to as uniform of a field as possible (called the 'flat field')
2. Use the detector response to the uniform field as a correction, presumably by normalizing the flat field to 1 (removing the zero/negative pixels).
3. Divide the measured images by the normalized flat field images.

Open Questions about this Issue

1. Is the flat field something that is collected relatively infrequently and used for many experiments?
2. Presumably, anyone who is using a 2-D detector will require this functionality?

Use PIMS for image series/stack file I/O

Description

We will need to read in multi-image sequences from disk. Instead of doing this ad-hoc, extend PIMS (https://github.com/soft-matter/pims) which already has code to provide a uniform interface to different file-format backends

There are already several implementations of tiff readers as well as png, jpeg, etc

There is a prototype of a sensible dispatch function as well (soft-matter/pims#74)

getting proper docs is on the TODO list over there (soft-matter/pims#26)

To add a new file type you just need to sub-class pims.FramesSequence and override 3 functions and 2 properties. You then get lazy loading + numpy style slicing.

Problems to Solve

These can either be done in our code base or pushed back up to PIMS

• determine if fabio deals with any file formats PIMS can't already do, and if so wrap for PIMS (should go to PIMS)
• if avizo can be multi-frame, wrap with PIMS (should probably go to PIMS)
• wrap pySpec with PIMS (should probably stay here)

Convert to reciprocal space

Description

Converting to reciprocal space from the lab frame is commonly used in two contexts: control and analysis.

In the control context, the desired position in reciprocal space is computed with the purpose of rotating the sample, the detector or both. This is a calculation whose speed is not necessarily that important. Also, this conversion is being handled by @arkilic for the controls side.

However, in the analysis context, every pixel in the measured 2-d image needs to be converted into reciprocal space. Assuming the detector is 2048x2048, there are 4 million pixels that need to be converted to reciprocal space. For a small data set, there might be ~100 images, so the reciprocal space conversion needs to be performed 400 million times. In this analysis context, speed becomes absolutely critical.

@stuwilkins wrote a fast implementation in C and wrapped by python to perform 6-circle conversions, called pyspec: https://github.com/stuwilkins/pyspec

Python function call: def processToQ(self) https://github.com/stuwilkins/pyspec/blob/master/pyspec/ccd/transformations.py#L1277

Python-C interface: static PyObject* ccdToQ(PyObject *self, PyObject *args, PyObject *kwargs)
https://github.com/stuwilkins/pyspec/blob/master/src/ctrans.c#L33

fast threaded C code: int calcQTheta(_float* diffAngles, _float theta, _float mu, _float *qTheta, _int n, _float lambda)
https://github.com/stuwilkins/pyspec/blob/master/src/ctrans.c#L209

Work in progress

Pull Request: #59

Problems to be solved for this issue

• Wrap the fast C code that @stuwilkins wrote into the NSLS2 library

pip requirements.txt

Description

We need to start keeping track of all the required packages for our libraries. I believe that if these requirements are contained in a file called requirements.txt in the root of the repository then pip will ensure that those packages are installed. We can target git, mercurial and svn repositories with this requirements.txt file and even specific branches, versions or commits. All of this information is located here.

Azimuthal Integration

Description

Azimuthal Integration is a process by which a 2-D image (x,y) is converted to a 1-D array based on rings around a known image center (xc, yc). The resultant 1-D array has x-values of radius from the image center and y-values of the average intensity around the rings at the given radius.

Problems that need to be solved for this issue

• Given an image center (xc, yc) compute a 1-D histogram
• Test the code

Register 1-D curves: Fit, Scale & Shift

Description

The Inelastic scattering beamline at NSLS-II, IXS (10-ID), has 6 strip detectors which nominally should produce the same datasets. However, due to measurement error, detector issues, crystal issues, etc., these six curves might not be identical. Given that the elastic scattering peak should be identical for each of the six measured curves, the elastic peak can be fit and used as a handle with which to scale and shift all curves. The elastic peak should be centered on the zero of the x-axis, since the elastic peak, by definition, has an energy transfer of 0 (ΔE=0). It is unclear how the width of the elastic peak should be determined.

Problems to be solved

• Determine elastic peak function that should be used to fit the elastic peak
• Fit elastic peak
• Shift peaks to be centered at ΔE=0
• Determine how to scale peaks to the appropriate width
• Scale peaks to the appropriate width

Automated Masking

Description

Automatic masking is a complicated problem. The default behavior should be:

• find saturated pixels
• find dead pixels

Optional behavior should be:

• find the beam stop
• find dead or malfunctioning scintillators

@jammc should be able to add quite a bit of verbosity to this description and some actual code that does this

Problems that need to be solved for this issue

• Coherent: Remove bad pixels from 2-D images
• Coherent: Remove beam stop and holder from 2-D images
• XPD: Automated masking for powder samples

improve test script

see how matplotlib + scikit-image + numpy do this

• add ability to use KnownFail tests
• add pep8 test

Quantities/units

It is possible to propagate dimensional units through computations. The code is aware that a variable has both a value and a unit so that if you try to add a wave length to an energy, it will raise an exception and if you try to add a value in mm to a value in m do auto-unit conversion.

If we decide we want to use this, we should use one of the existing packages rather than rolling our own.

• https://pypi.python.org/pypi/quantities
• http://astropy.readthedocs.org/en/latest/units/

Image Corrections: General

The more basic image corrections (flat field, dark current, background, etc...) are very simple from the analysis perspective.

The flat field correction requires a normalization, removal of zeros and then it becomes the divisor for the experimental images (see #26)

The dark current is simply the experimental image minus the dark current. The tricky part is in correctly aligning the dark current images with the experimental images. (see #24)

The background requires a normalization for relative beam time exposure and then it is subtracted from the experimental image.

In all three of these cases, the tricky part is to correctly align the experimental image and the correction image. For this reason, I think there should simply be helper functions to parse the nested dictionary we get back from the arkilic/metadataStore that will align the experimental images with the appropriate correction images. The real question is where to put these functions? In nsls2/broker.py? As another module/folder in arkilic/metadataStore? Thoughts @tacaswell?

Calibrate 2-D Detector

Description

Corrections Required

Simple array math:

• gain (flat field)
• dark current
• normalize by beam intensity

Masking:

• Saturated pixels
• Dead pixels
• beam stop
• beam stop holder
• dead scintillators
• zingers (cosmic rays)
• other bad stuff

Warping:

• Pitch, roll and yaw: tilt angles.
  • if y is the beam direction, x is parallel to the ground and z is perpendicular to the ground, roll is the angle around y, yaw is the angle around z and pitch is the angle around x
• Grid/point distortion
  • example

Problems to be solved

Calibrate

• Fit ellipses to diffraction patterns of calibration standards (things like Si, Lab6, CeO2, etc...) to obtain pitch and yaw, refine detector distance, beam center
• Algorithmically find dead pixels (for CCD)
• Algorithmically find dead pixels (for insert_other_thing_here)
• Algorithmcially find saturated pixels (for CCD) (has implications for peak fitting, e.g., the peak is truncated on a Bragg peak, or something)
• Algorithmically find saturated pixels (for insert_other_thing_here)

Correct

• Apply mask with known pixel issues to image
• Algorithmically find beam stop holder
• Compute pixel-by-pixel positions given pitch, yaw and roll

See #10:

• Gain ()
• Dark current (array subtraction)
• Norm by beam intensity (array division by scalar)

See Also

#30
#11
#32
#37