From f7bdc2acff3c13a6d632c28c4569690ab106eed7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Picca=20Fr=C3=A9d=C3=A9ric-Emmanuel?= Date: Fri, 18 Aug 2017 14:48:52 +0200 Subject: Import Upstream version 0.5.0+dfsg --- doc/source/Tutorials/Sift/sift.rst | 315 ++++++++++ doc/source/Tutorials/array_widget.rst | 246 ++++++++ doc/source/Tutorials/fit.rst | 635 +++++++++++++++++++++ doc/source/Tutorials/fitconfig.rst | 207 +++++++ doc/source/Tutorials/img/arraywidget3D_0.png | Bin 0 -> 35724 bytes doc/source/Tutorials/img/arraywidget3D_1.png | Bin 0 -> 23956 bytes doc/source/Tutorials/img/arraywidget5D_0.png | Bin 0 -> 34473 bytes doc/source/Tutorials/img/arraywidget5D_1.png | Bin 0 -> 24553 bytes .../Tutorials/img/custom_config_scale0.5.png | Bin 0 -> 186529 bytes .../Tutorials/img/custom_config_scale1.0.png | Bin 0 -> 181099 bytes .../Tutorials/img/custom_config_scale2.1.png | Bin 0 -> 183066 bytes doc/source/Tutorials/img/fitwidget1.png | Bin 0 -> 31627 bytes doc/source/Tutorials/img/fitwidget2.png | Bin 0 -> 106983 bytes doc/source/Tutorials/img/fitwidget3.png | Bin 0 -> 26141 bytes doc/source/Tutorials/img/fitwidget4.png | Bin 0 -> 121846 bytes doc/source/Tutorials/img/fitwidget5.png | Bin 0 -> 29784 bytes doc/source/Tutorials/img/stripbg_plot1.png | Bin 0 -> 41420 bytes doc/source/Tutorials/img/stripbg_plot2.png | Bin 0 -> 45619 bytes doc/source/Tutorials/specfile_to_hdf5.rst | 323 +++++++++++ 19 files changed, 1726 insertions(+) create mode 100644 doc/source/Tutorials/Sift/sift.rst create mode 100644 doc/source/Tutorials/array_widget.rst create mode 100644 doc/source/Tutorials/fit.rst create mode 100644 doc/source/Tutorials/fitconfig.rst create mode 100644 doc/source/Tutorials/img/arraywidget3D_0.png create mode 100644 doc/source/Tutorials/img/arraywidget3D_1.png create mode 100644 doc/source/Tutorials/img/arraywidget5D_0.png create mode 100644 doc/source/Tutorials/img/arraywidget5D_1.png create mode 100644 doc/source/Tutorials/img/custom_config_scale0.5.png create mode 100644 doc/source/Tutorials/img/custom_config_scale1.0.png create mode 100644 doc/source/Tutorials/img/custom_config_scale2.1.png create mode 100644 doc/source/Tutorials/img/fitwidget1.png create mode 100644 doc/source/Tutorials/img/fitwidget2.png create mode 100644 doc/source/Tutorials/img/fitwidget3.png create mode 100644 doc/source/Tutorials/img/fitwidget4.png create mode 100644 doc/source/Tutorials/img/fitwidget5.png create mode 100644 doc/source/Tutorials/img/stripbg_plot1.png create mode 100644 doc/source/Tutorials/img/stripbg_plot2.png create mode 100644 doc/source/Tutorials/specfile_to_hdf5.rst (limited to 'doc/source/Tutorials') diff --git a/doc/source/Tutorials/Sift/sift.rst b/doc/source/Tutorials/Sift/sift.rst new file mode 100644 index 0000000..7168cf3 --- /dev/null +++ b/doc/source/Tutorials/Sift/sift.rst @@ -0,0 +1,315 @@ + +SIFT image alignment tutorial +============================= + +SIFT (Scale-Invariant Feature Transform) is an algorithm developped by +David Lowe in 1999. It is a worldwide reference for image alignment and +object recognition. The robustness of this method enables to detect +features at different scales, angles and illumination of a scene. + +Silx provides an implementation of SIFT in OpenCL, meaning that it can +run on Graphics Processing Units and Central Processing Units as well. +Interest points are detected in the image, then data structures called +*descriptors* are built to be characteristic of the scene, so that two +different images of the same scene have similar descriptors. They are +robust to transformations like translation, rotation, rescaling and +illumination change, which make SIFT interesting for image stitching. + +In the fist stage, descriptors are computed from the input images. Then, +they are compared to determine the geometric transformation to apply in +order to align the images. This implementation can run on most graphic +cards and CPU, making it usable on many setups. OpenCL processes are +handled from Python with PyOpenCL, a module to access OpenCL parallel +computation API. + +This tutuorial explains the three subsequent steps: + +- keypoint extraction +- Keypoint matching +- image alignment + +All the tutorial has been made using the Jupyter notebook. + +.. code:: python + + import time + start_time = time.time() + %pylab nbagg + + +.. parsed-literal:: + + Populating the interactive namespace from numpy and matplotlib + + +.. code:: python + + # display test image + import silx + print("Silx version %s"%silx.version) + + from PIL import Image + from silx.test.utils import utilstest + path = utilstest.getfile("lena.png") + image = numpy.asarray(Image.open(path)) + fig, ax = subplots() + ax.imshow(image, cmap="gray") + + +.. parsed-literal:: + + Silx version 0.5.0-dev0 + + + +.. parsed-literal:: + + + + + +.. raw:: html + + + + + + +.. parsed-literal:: + + + + + +.. code:: python + + #Initialization of the sift object is time consuming: it compiles all the code. + import os + #set to 1 to see the compilation going on + os.environ["PYOPENCL_COMPILER_OUTPUT"] = "0" + from silx.image import sift + #switch to "GPU" to "CPU" to enable fail-save version. + %time sift_ocl = sift.SiftPlan(template=image, devicetype="GPU") + + +.. parsed-literal:: + + CPU times: user 112 ms, sys: 112 ms, total: 224 ms + Wall time: 240 ms + + +.. code:: python + + print("Time for calculating the keypoints on one image of size %sx%s"%image.shape[:2]) + %time keypoints = sift_ocl(image) + print("Number of keypoints: %s"%len(keypoints)) + print("Keypoint content:") + print(keypoints.dtype) + print("x: %.3f \t y: %.3f \t sigma: %.3f \t angle: %.3f" % + (keypoints[-1].x,keypoints[-1].y,keypoints[-1].scale,keypoints[-1].angle)) + print("descriptor:") + print(keypoints[-1].desc) + + +.. parsed-literal:: + + Time for calculating the keypoints on one image of size 512x512 + CPU times: user 652 ms, sys: 0 ns, total: 652 ms + Wall time: 649 ms + Number of keypoints: 411 + Keypoint content: + (numpy.record, [('x', ' + + + +.. raw:: html + + + + + + +.. parsed-literal:: + + [] + + + +.. code:: python + + #Diplaying keypoints by scale: + fig, ax = subplots() + ax.hist(keypoints[:].scale, 100) + ax.set_xlabel("scale") + + + +.. parsed-literal:: + + + + + +.. raw:: html + + + + + + +.. parsed-literal:: + + + + + +.. code:: python + + #One can see 3 groups of keypoints, boundaries at: 8 and 20. Let's display them using colors. + S = 8 + L = 20 + tiny = keypoints[keypoints[:].scale=S)] + bigger = keypoints[keypoints[:].scale>=L] + + fig, ax = subplots() + ax.imshow(image, cmap="gray") + ax.plot(tiny[:].x, tiny[:].y,",g", label="tiny") + ax.plot(small[:].x, small[:].y,".b", label="small") + ax.plot(bigger[:].x, bigger[:].y,"or", label="large") + ax.legend() + + + +.. parsed-literal:: + + + + + +.. raw:: html + + + + + + +.. parsed-literal:: + + + + + +Image matching and alignment +---------------------------- + +Matching can also be performed on the device (GPU) as every single +keypoint from an image needs to be compared with all keypoints from the +second image. + +In this simple example we will simple offset the first image by a few +pixels + +.. code:: python + + shifted = numpy.zeros_like(image) + shifted[5:,8:] = image[:-5, :-8] + shifted_points = sift_ocl(shifted) + +.. code:: python + + %time mp = sift.MatchPlan() + %time match = mp(keypoints, shifted_points) + print("Number of Keypoints with for image 1 : %i, For image 2 : %i, Matching keypoints: %i" % (keypoints.size, shifted_points.size, match.shape[0])) + from numpy import median + print("Measured offsets dx: %.3f, dy: %.3f"%(median(match[:,1].x-match[:,0].x),median(match[:,1].y-match[:,0].y))) + + +.. parsed-literal:: + + CPU times: user 20 ms, sys: 128 ms, total: 148 ms + Wall time: 167 ms + CPU times: user 8 ms, sys: 4 ms, total: 12 ms + Wall time: 14.9 ms + Number of Keypoints with for image 1 : 411, For image 2 : 399, Matching keypoints: 351 + Measured offsets dx: 8.000, dy: 5.000 + + +.. code:: python + + # Example of usage of the automatic alignment: + import scipy.ndimage + rotated = scipy.ndimage.rotate(image, 20, reshape=False) + sa = sift.LinearAlign(image) + fig,ax = subplots(1, 3, figsize=(12,4)) + ax[0].imshow(image) + ax[1].imshow(rotated) + ax[2].imshow(sa.align(rotated)) + + +.. parsed-literal:: + + /scisoft/users/jupyter/jupy34/lib/python3.4/site-packages/pyopencl/cffi_cl.py:1469: CompilerWarning: Non-empty compiler output encountered. Set the environment variable PYOPENCL_COMPILER_OUTPUT=1 to see more. + "to see more.", CompilerWarning) + + + +.. parsed-literal:: + + + + + +.. raw:: html + + + + + + +.. parsed-literal:: + + + + + +References +~~~~~~~~~~ + +- David G. Lowe, Distinctive image features from scale-invariant + keypoints, International Journal of Computer Vision, vol. 60, no 2, + 2004, p. 91–110 - "http://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf" + +.. code:: python + + print("Total execution time: %.3fs" % (time.time() - start_time)) + + +.. parsed-literal:: + + Total execution time: 6.190s + diff --git a/doc/source/Tutorials/array_widget.rst b/doc/source/Tutorials/array_widget.rst new file mode 100644 index 0000000..f2ccfdf --- /dev/null +++ b/doc/source/Tutorials/array_widget.rst @@ -0,0 +1,246 @@ + +.. currentmodule:: silx.gui.data.ArrayTableWidget + +ArrayTableWidget +================ + +:class:`ArrayTableWidget` is a widget designed to visualize numpy arrays or h5py datasets. + +3D example +---------- + +Let's look at a simple usage example: + +.. code-block:: python + + from silx.gui import qt + from silx.gui.data.ArrayTableWidget import ArrayTableWidget + import numpy + array = numpy.arange(1000) + array.shape = (5, 10, 20) + app = qt.QApplication([]) + w = ArrayTableWidget() + w.setArrayData(array, labels=True) + w.show() + app.exec_() + + +.. |imgArray0| image:: img/arraywidget3D_0.png + :height: 300px + :align: middle + +|imgArray0| + +We get a widget that allows us to see a *slice*, or a *frame*, of 2D data +with 10 lines and 20 columns, in a 3D array (5 x 10 x 20). +The column index corresponds to the last dimension of the array, and the +row index corresponds to the second to last dimension. The first index can be browsed +using a slider, icons or a text entry for access to any given slice among the 5 available. + +The parameter ``labels=True`` of :meth:`setArrayData` causes the browser to be labeled +*Dimension 0*. + +If we want to see slices in different perspective, we can use +:meth:`ArrayTableWidget.setPerspective`. The perspective is defined as the list +of dimensions that are not represented in the frame, orthogonal to it. +For a 3D array, there are 3 possible perspectives: *[0, ]* (the default perspective), +*[1, ]* and *[2, ]*. + +Lets change the perspective: + +.. code-block:: python + + w.setPerspective([1]) + +.. |imgArray1| image:: img/arraywidget3D_1.png + :height: 300px + :align: middle + +|imgArray1| + +What we see now is a frame of *5 x 20* values, and the browser now browses the second dimension +to select one of 10 available frames. The label is updated accordingly to show *Dimension 1*. + +To select a different frame programmatically, without using the browser, you can +use the :meth:`ArrayTableWidget.setIndex` method. To select the 9-th frame, use: + +.. code-block:: python + + w.setIndex([8]) + +More dimensions +--------------- + +This widget can be used for arrays with any numbers of dimensions. Let's create +a 5-dimensional array and display it: + +.. code-block:: python + + array = numpy.arange(10000) + array.shape = (5, 2, 10, 5, 20) + w.setArrayData(array, labels=True) + +.. |imgArray2| image:: img/arraywidget5D_0.png + :height: 300px + :align: middle + +|imgArray2| + +We now have 3 frames browsers, one for each one of the orthogonal dimensions. + +Let's look at a frame whose axes are along the second +and the fourth dimension, by setting the orthogonal axes to the first, +third and fifth dimensions: + +.. code-block:: python + + w.setPerspective([0, 2, 4]) + +.. |imgArray3| image:: img/arraywidget5D_1.png + :height: 300px + :align: middle + +|imgArray3| + + +Listing all the orthogonal dimensions might not feel very convenient for arrays +with more than 3 or 4 dimensions. +Fortunately, you can use the opposite approach of defining the two axes +parallel to the frame, using :meth:`ArrayTableWidget.setFrameAxes`: + +.. code-block:: python + + w.setFrameAxes(row_axis=1, col_axis=3) + +This achieves the exact same result as ``w.setPerspective([0, 2, 4])``. + +.. note:: + + Currently you cannot switch the row and column axes. The row axis + is always the lowest free dimension and the column axis is the + highest one with the current implementation. + So setting ``w.setFrameAxes(row_axis=3, col_axis=1)`` would not modify + the table axes. + + For the same reason, the order of the dimensions given as parameter to + :meth:`setPerspective` is not significant. + +To select a frame programmaticaly, you can again use :meth:`setFrameIndex`. +This time you must provide 3 unique indices: + +.. code-block:: python + + w.setIndex([2, 5, 14]) + +The 3 indices relate to the first, third and fifth dimensions. + +The frame index must always be defined as indices on the orthogonal axes/dimensions, +as defined by the *perspective*. + +Editing the data +---------------- + +By default, the data displayed in the table view can be edited. If you modify +a cell with a valid value, it will be modified in the internal data model. + +You can get the modified data with the following line: + +.. code-block:: python + + newdata = w.getData() + +This will give you a copy of the data, by default. + +If you want the data to be read-only, not editable, you must specify it when +you set the data: + +.. code-block:: python + + w.setDataArray(array, editable=False) + +More performances +----------------- + +By default, the method :meth:`setArrayData` creates a copy of the data array +for internal storage. This ensures that the original data object is not +modified when a cell of the table is changed interactively in the widget. + +This behavior has a negative impact on performances, especially for large data arrays. +To avoid this, you can explicitely disable the copy operation when setting the data: + +.. code-block:: python + + w.setArrayData(array, copy=False) + +The internal data array used by the widget is then a reference +to the same data object as the original *array*. The memory is shared and +is not duplicated. + +.. warning:: + + This can cause side-effects, if your original array is re-used elsewhere + in your program. + +Similarly, you can pass *copy=False* to the :meth:`getData` method, to avoid +doing a data copy operation: + +.. code-block:: python + + newdata = w.getData(copy=False) + +The variable *newdata* is then a reference to the internal widget data. + +.. warning:: + + Modifying the internal data used by the widget can have unpredictable + consequences. + +Background color +---------------- + +You can set the background color for each cell by passing a numpy array of +RGB colors to the :meth:`setArrayColors` method. + +The colors array must have one more dimension than the data array. This dimension +must be of length 3 for RGB colors or length 4 for RGBA colors. + +The colors array associates 3 (or 4) integers between 0 and 255 to each value +in the data array. The values represent the red, green, blue and alpha (opacity) +channels. + +In the following examples, we create a table displaying a complete palette +of RGB colors. + +.. code-block:: python + + import numpy + from silx.gui import qt + from silx.gui.widgets.ArrayTableWidget import ArrayTableWidget + + # data array + data = numpy.arange(256**3) + data.shape = 256, 256, 256 + + # RGB colors array + bcolors = numpy.empty((256, 256, 256, 3), dtype=numpy.uint8) + # fill red channel + bcolors[..., 0] = data[:] & 255 + # green + bcolors[..., 1] = (data[:] & (255 << 8)) >> 8 + # blue + bcolors[..., 2] = (data[:] & (255 << 16)) >> 16 + + # make text contrast with background (XOR) + fcolors = numpy.bitwise_xor(bcolors, 255) + + app = qt.QApplication([]) + + atw = ArrayTableWidget() + atw.setArrayData(data, copy=False) + atw.setArrayColors(bgcolors=bcolors, + fgcolors=fcolors) + atw.show() + + app.exec_() + + diff --git a/doc/source/Tutorials/fit.rst b/doc/source/Tutorials/fit.rst new file mode 100644 index 0000000..1305299 --- /dev/null +++ b/doc/source/Tutorials/fit.rst @@ -0,0 +1,635 @@ + +.. _fit-tutorial: + +Fit tools +--------- + +.. contents:: :local: + +.. _leastsq-tutorial: + +Using :func:`leastsq` ++++++++++++++++++++++ + +.. currentmodule:: silx.math.fit + +Running an iterative fit with :func:`leastsq` involves the following steps: + + - designing a fit model function that has the signature ``f(x, ...)``, + where ``x`` is an array of values of the independant variable and all + remaining parameters are the parameters to be fitted + - defining the sequence of initial values for all parameters to be fitted. + You can usually start with ``[1., 1., ...]`` if you don't know a better + estimate. The algorithm is robust enough to converge to a solution most + of the time. + - setting constraints (optional) + +Data required to perform a fit is: + + - an array of ``x`` values (abscissa, independant variable) + - an array of ``y`` data points + - the ``sigma`` array of uncertainties associated to each data point. + This is optional, by default each data point gets assigned a weight of 1. + +Standard fit +************ + +Let's demonstrate this process in a short example, using synthetic data. +We generate an array of synthetic data using a polynomial function of degree 4, +and try to use :func:`leastsq` to find back the functions parameters. + +.. code-block:: python + + import numpy + from silx.math.fit import leastsq + + # create some synthetic polynomial data + x = numpy.arange(1000) + y = 2.4 * x**4 - 10. * x**3 + 15.2 * x**2 - 24.6 * x + 150. + + # define our fit function: a generic polynomial of degree 4 + def poly4(x, a, b, c, d, e): + return a * x**4 + b * x**3 + c * x**2 + d * x + e + + # The fit is an iterative process that requires an initial + # estimation of the parameters. Let's just use 1s. + initial_parameters = numpy.array([1., 1., 1., 1., 1.]) + + # Run fit + fitresult = leastsq(model=poly4, + xdata=x, + ydata=y, + p0=initial_parameters, + full_output=True) + + # leastsq with full_output=True returns 3 objets + optimal_parameters, covariance, infodict = fitresult + # the first object is an array with the fitted parameters + a, b, c, d, e = optimal_parameters + + print("Fit took %d iterations" % infodict["niter"]) + print("Reduced chi-square: %f" % infodict["reduced_chisq"]) + print("Theoretical parameters:\n\t" + + "a=2.4, b=-10, c=15.2, d=-24.6, e=150") + print("Optimal parameters for y2 fitting:\n\t" + + "a=%f, b=%f, c=%f, d=%f, e=%f" % (a, b, c, d, e)) + +The output of this program is:: + + Fit took 35 iterations + Reduced chi-square: 682592.670690 + Theoretical parameters: + a=2.4, b=-10, c=15.2, d=-24.6, e=150 + Optimal parameters for y fitting: + a=2.400000, b=-9.999665, c=14.970422, d=31.683448, e=-3216.131136 + +We can see that this fit result is poor. In particular, parameters ``d`` and ``e`` +are very poorly fitted. +This is most likely due to numerical rounding errors. As we are dealing with +very large values in our ``y`` array, we are affected by the limits of how +floating point numbers are represented by computers. The larger a value, the +larger its rounding error. + +If you limit the ``x`` range to deal with +smaller ``y`` values, the fit result becomes perfect. In our example, replacing ``x`` +with:: + + x = numpy.arange(100) + +produces the following result:: + + Fit took 9 iterations + Reduced chi-square: 0.000000 + Theoretical parameters: + a=2.4, b=-10, c=15.2, d=-24.6, e=150 + Optimal parameters for y fitting: + a=2.400000, b=-10.000000, c=15.200000, d=-24.600000, e=150.000000 + + + +Constrained fit +*************** + +But let's revert back to our initial ``x = numpy.arange(1000)``, to experiment +with different approaches to improving the fit. + +The :func:`leastsq` functions provides +a way to set constraints on parameters. You can for instance assert that a given +parameter must remain equal to it's initial value, or define an acceptable range +for it to vary, or decide that a parameter must be equal to another parameter +multiplied by a certain factor. This is very useful in cases in which you have +enough knowledge to make reasonable assumptions on some parameters. + +In our case, we will set constraints on ``d`` and ``e``. We will quote ``d`` to +stay in the range between -25 and -24, and fix ``e`` to 150. + +Replace the call to :func:`leastsq` by following lines: + +.. code-block:: python + + # Define constraints + cons = [[0, 0, 0], # a: no constraint + [0, 0, 0], # b: no constraint + [0, 0, 0], # c: no constraint + [2, -25., -23.], # -25 < d < -24 + [3, 0, 0]] # e is fixed to initial value + fitresult = leastsq(poly4, x, y, + # initial values must be consistent with constraints + p0=[1., 1., 1., -24., 150.], + constraints=cons, + full_output=True) + +The output of this program is:: + + Constrained fit took 100 iterations + Reduced chi-square: 3.749280 + Theoretical parameters: + a=2.4, b=-10, c=15.2, d=-24.6, e=150 + Optimal parameters for y fitting: + a=2.400000, b=-9.999999, c=15.199648, d=-24.533014, e=150.000000 + +The chi-square value is much improved and the results are much better, at the +cost of more iterations. + +Weighted fit +************ +A third approach to improve our fit is to define uncertainties for the data. +The larger the uncertainty on a data sample, the smaller its weight will be +in the least-square problem. + +In our case, we do not know the uncertainties associated to our data. We could +determine the uncertainties due to numerical rounding errors, but let's just use +a common approach that requires less work: use the square-root of the data values +as their uncertainty value: + +.. code-block:: python + + sigma = numpy.sqrt(y) + + # Fit y + fitresult = leastsq(model=poly4, + xdata=x, + ydata=y, + sigma=sigma, + p0=initial_parameters, + full_output=True) + +This results in a great improvement:: + + Weighted fit took 6 iterations + Reduced chi-square: 0.000000 + Theoretical parameters: + a=2.4, b=-10, c=15.2, d=-24.6, e=150 + Optimal parameters for y fitting: + a=2.400000, b=-10.000000, c=15.200000, d=-24.600000, e=150.000000 + +The resulting fit is perfect. The very large ``y`` values with their very large +associated uncertainties have been practicaly rejected from the fit process. The fit +converged even faster than with the solution of limiting the ``x`` range to +0 -- 100. + +.. _fitmanager-tutorial: + +Using :class:`FitManager` ++++++++++++++++++++++++++ + +.. currentmodule:: silx.math.fit.fitmanager + +A :class:`FitManager` is a tool that provides a way of handling fit functions, +associating estimation functions to estimate the initial parameters, modify +the configuration parameters for the fit (enabling or disabling weights...) or +for the estimation function, and choosing a background model. + +It provides an abstraction layer on top of :func:`leastsq`. + +Weighted polynomial fit +*********************** + +The following program accomplishes the same weighted fit of a polynomial as in +the previous tutorial (`Weighted fit`_) + +.. code-block:: python + + import numpy + from silx.math.fit.fitmanager import FitManager + + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(1000).astype(numpy.float) + y = 2.4 * x**4 - 10. * x**3 + 15.2 * x**2 - 24.6 * x + 150. + + # define our fit function: a generic polynomial of degree 4 + def poly4(x, a, b, c, d, e): + return a * x**4 + b * x**3 + c * x**2 + d * x + e + + # define an estimation function to that returns initial parameters + # and constraints + def esti(x, y): + p0 = numpy.array([1., 1., 1., 1., 1.]) + cons = numpy.zeros(shape=(5, 3)) + return p0, cons + + # Fitting + fit = FitManager() + fit.setdata(x=x, y=y) + + fit.addtheory("polynomial", + function=poly4, + # any list of 5 parameter names would be OK + parameters=["A", "B", "C", "D", "E"], + estimate=esti) + fit.settheory('polynomial') + fit.configure(WeightFlag=True) + fit.estimate() + fit.runfit() + + print("\n\nFit took %d iterations" % fit.niter) + print("Reduced chi-square: %f" % fit.chisq) + print("Theoretical parameters:\n\t" + + "a=2.4, b=-10, c=15.2, d=-24.6, e=150") + a, b, c, d, e = (param['fitresult'] for param in fit.fit_results) + print("Optimal parameters for y2 fitting:\n\t" + + "a=%f, b=%f, c=%f, d=%f, e=%f" % (a, b, c, d, e)) + + +The result is the same as in our weighted :func:`leastsq` example, +as expected:: + + Fit took 6 iterations + Reduced chi-square: 0.000000 + Theoretical parameters: + a=2.4, b=-10, c=15.2, d=-24.6, e=150 + Optimal parameters for y2 fitting: + a=2.400000, b=-10.000000, c=15.200000, d=-24.600000, e=150.000000 + +Fitting gaussians +***************** + +The :class:`FitManager` object is especially useful for fitting multi-peak +gaussian-shaped spectra. The *silx* module :mod:`silx.math.fit.fittheories` +provides fit functions and their associated estimation functions that are +specifically designed for this purpose. + +These fit functions can handle variable number of parameters defining a +variable number of peaks, and the estimation functions use a peak detection +algorithm to determine how many initial parameters must be returned. + +For the sake of the example, let's test the multi-peak fitting on synthetic +data, generated using another *silx* module: :mod:`silx.math.fit.functions`. + +.. code-block:: python + + import numpy + from silx.math.fit.functions import sum_gauss + from silx.math.fit import fittheories + from silx.math.fit.fitmanager import FitManager + + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(1000).astype(numpy.float) + + # height, center x, fwhm + p = [1000, 100., 250, # 1st peak + 255, 690., 45, # 2nd peak + 1500, 800.5, 95] # 3rd peak + + y = sum_gauss(x, *p) + + # Fitting + fit = FitManager() + fit.setdata(x=x, y=y) + fit.loadtheories(fittheories) + fit.settheory('Gaussians') + fit.estimate() + fit.runfit() + + print("Searched parameters = %s" % p) + print("Obtained parameters : ") + dummy_list = [] + for param in fit.fit_results: + print(param['name'], ' = ', param['fitresult']) + dummy_list.append(param['fitresult']) + print("chisq = ", fit.chisq) + +And the result of this program is:: + + Searched parameters = [1000, 100.0, 250, 255, 690.0, 45, 1500, 800.5, 95] + Obtained parameters : + ('Height1', ' = ', 1000.0) + ('Position1', ' = ', 100.0) + ('FWHM1', ' = ', 250.0) + ('Height2', ' = ', 255.0) + ('Position2', ' = ', 690.0) + ('FWHM2', ' = ', 44.999999999999993) + ('Height3', ' = ', 1500.0) + ('Position3', ' = ', 800.5) + ('FWHM3', ' = ', 95.000000000000014) + ('chisq = ', 0.0) + +In addition to gaussians, we could have fitted several other similar type of +functions: asymetric gaussian functions, lorentzian functions, +Pseudo-Voigt functions or hypermet tailing functions. + +The :meth:`loadtheories` method can also be used to load user defined +functions. Instead of a module, a path to a Python source file can be given +as a parameter. This source file must adhere to certain conventions, as explained +in the documentation of :mod:`silx.math.fit.fittheories` and +:mod:`silx.math.fit.fittheory.FitTheory`. + +Subtracting a background +************************ + +:class:`FitManager` provides a few standard background theories, for cases when +a background signal is superimposed on the multi-peak spectrum. + +For example, let's add a linear background to our synthetic data, and see how +:class:`FitManager` handles the fitting. + +In our previous example, redefine ``y`` as follows: + +.. code-block:: python + + p = [1000, 100., 250, + 255, 690., 45, + 1500, 800.5, 95] + y = sum_gauss(x, *p) + # add a synthetic linear background + y += 0.13 * x + 100. + +Before the line ``fit.estimate()``, add the following line: + +.. code-block:: python + + fit.setbackground('Linear') + +The result becomes:: + + Searched parameters = [1000, 100.0, 250, 255, 690.0, 45, 1500, 800.5, 95] + Obtained parameters : + ('Constant', ' = ', 100.00000000000001) + ('Slope', ' = ', 0.12999999999999998) + ('Height1', ' = ', 1000.0) + ('Position1', ' = ', 100.0) + ('FWHM1', ' = ', 249.99999999999997) + ('Height2', ' = ', 255.00000000000003) + ('Position2', ' = ', 690.0) + ('FWHM2', ' = ', 44.999999999999993) + ('Height3', ' = ', 1500.0) + ('Position3', ' = ', 800.5) + ('FWHM3', ' = ', 95.0) + ('chisq = ', 3.1789004676997597e-27) + +The available background theories are: *Linear*, *Constant* and *Strip*. + +The strip background is a popular background model that can compute and +subtract any background shape as long as its curvature is significantly +lower than the peaks' curvature. In other words, as long as the background +signal is significantly smoother than the actual signal, it can be easily +computed. + +The main parameters required by the strip function are the strip width *w* +and the number of iterations. At each iteration, if the contents of channel *i*, +``y(i)``, is above the average of the contents of the channels at *w* channels of +distance, ``y(i-w)`` and ``y(i+w)``, ``y(i)`` is replaced by the average. +At the end of the process we are left with something that resembles a spectrum +in which the peaks have been "stripped". + +The following example illustrates the strip background removal process: + +.. code-block:: python + + from silx.sx import plot + from silx.gui import qt + import numpy + from silx.math.fit.filters import strip + from silx.math.fit.functions import sum_gauss + + x = numpy.arange(5000) + # (height1, center1, fwhm1, ...) 5 peaks + params1 = (50, 500, 100, + 20, 2000, 200, + 50, 2250, 100, + 40, 3000, 75, + 23, 4000, 150) + y0 = sum_gauss(x, *params1) + + # random values between [-1;1] + noise = 2 * numpy.random.random(5000) - 1 + # make it +- 5% + noise *= 0.05 + + # 2 gaussians with very large fwhm, as background signal + actual_bg = sum_gauss(x, 15, 3500, 3000, 5, 1000, 1500) + + # Add 5% random noise to gaussians and add background + y = y0 * (1 + noise) + actual_bg + + # compute strip background model + strip_bg = strip(y, w=5, niterations=5000) + + # plot results + app = qt.QApplication([]) + plot(x, y, x, actual_bg, x, strip_bg) + plot(x, y, x, (y - strip_bg)) + app.exec_() + +.. |imgStrip1| image:: img/stripbg_plot1.png + :height: 300px + :align: middle + +.. |imgStrip2| image:: img/stripbg_plot2.png + :height: 300px + :align: middle + +.. list-table:: + :widths: 1 2 + + * - |imgStrip1| + - Data with background in black (``y``), actual background in red, computed strip + background in green + * - |imgStrip2| + - Data with background in blue, data after subtracting strip background in black + +The strip also removes the statistical noise, so the computed strip background +will be slightly lower than the actual background. This can be solved by +performing a smoothing prior to the strip computation. + +See the `PyMca documentation `_ +for more information on the strip background. + +To configure the strip background model of :class:`FitManager`, use its :meth:`configure` +method to modify the following parameters: + + - *StripWidth*: strip width parameter *w*, mentionned earlier + - *StripNIterations*: number of iterations + - *StripThresholdFactor*: if this parameter is left to its default value 1, + the algorithm behaves as explained earlier: ``y(i)`` is compared to the average of + ``y(i-1)`` and ``y(i+1)``. + If this factor is set to another value *f*, ``y(i)`` is compared to the + average multiplied by ``f``. + - *SmoothStrip*: if this parameter is set to ``True``, a smoothing is applied + prior to the strip. + + +These parameters can be modified like this: + +.. code-block:: python + + # ... + fit.settheory('Strip') + fit.configure(StripWidth=5, + StripNIterations=5000, + StripThresholdFactor=1.1, + SmoothStrip=True) + # ... + +Using a strip background has performance implications. You should try to keep +the number of iterations as low as possible if you need to run batch fitting +using this model. Increasing the strip width can help reducing the number of +iterations, with the risk of underestimating the background signal. + +.. _fitwidget-tutorial: + +Using :class:`FitWidget` +++++++++++++++++++++++++ + +.. currentmodule:: silx.gui.fit.FitWidget + +Simple usage +************ + +The :class:`FitWidget` is a graphical interface for :class:`FitManager`. + +.. code-block:: python + + import numpy + from silx.gui import qt + from silx.gui.fit import FitWidget + from silx.math.fit.functions import sum_gauss + + x = numpy.arange(2000).astype(numpy.float) + constant_bg = 3.14 + + # gaussian parameters: height, position, fwhm + p = numpy.array([1000, 100., 30.0, + 500, 300., 25., + 1700, 500., 35., + 750, 700., 30.0, + 1234, 900., 29.5, + 302, 1100., 30.5, + 75, 1300., 75.]) + y = sum_gauss(x, *p) + constant_bg + + a = qt.QApplication([]) + + w = FitWidget() + w.setData(x=x, y=y) + w.show() + + a.exec_() + +.. |imgFitWidget1| image:: img/fitwidget1.png + :width: 300px + :align: middle + +.. |imgFitWidget2| image:: img/fitwidget2.png + :width: 300px + :align: middle + +.. |imgFitWidget3| image:: img/fitwidget3.png + :width: 300px + :align: middle + +.. |imgFitWidget4| image:: img/fitwidget4.png + :width: 300px + :align: middle + +Executing this code opens the following widget. + + |imgFitWidget1| + +The functions you can choose from are the standard gaussian-shaped functions +from :mod:`silx.math.fit.fittheories`. At the top of the list, you will find +the *Add Function(s)* option, that allows you to load your user defined fit +theories from a *.py* source file. + +After selecting the *Constant* background model and clicking the *Estimate* +button, the widget displays this: + + |imgFitWidget2| + +The 7 peaks have been detected, and their parameters estimated. +Also, the estimation function defined some constraints (positive height and +positive full-width at half-maximum). + +You can modify the values in the estimation column of the table, to use different +initial parameters for the fit. + +The individual constraints can be modified prior to fitting. It is also possible to +modify the constraints globally by clicking the *Configure* button' to open a +configuration dialog. To get help on the meaning of the various parameters, +hover the mouse on the corresponding check box or entry widget, to display a +tooltip help message. + + |imgFitWidget3| + +The other configuration tabs can be modified to change the peak search parameters +and the strip background parameters prior to the estimation. +After closing the configuration dialog, you must re-run the estimation +by clicking the *Estimate* button. + +After all configuration parameters and all constrants are set according to your +preferences, you can click the *Start Fit* button. This runs the fit and displays +the results in the *Fit Value* column of the table. + + |imgFitWidget4| + +Customizing the functions +************************* + +.. |imgFitWidget5| image:: img/fitwidget5.png + :width: 300px + :align: middle + +The :class:`FitWidget` can be initialized with a non-standard +:class:`FitManager`, to customize the available functions. + +.. code-block:: python + + from silx.gui import qt + from silx.math.fit import FitManager + from silx.gui.fit import FitWidget + + def linearfun(x, a, b): + return a * x + b + + # create synthetic data for the example + x = list(range(0, 100)) + y = [linearfun(x_, 2.0, 3.0) for x_ in x] + + # we need to create a custom fit manager and add our theory + myfitmngr = FitManager() + myfitmngr.setdata(x, y) + myfitmngr.addtheory("my linear function", + function=linearfun, + parameters=["a", "b"]) + + a = qt.QApplication([]) + + # our fit widget can now use our custom fit manager + fw = FitWidget(fitmngr=myfitmngr) + fw.show() + + a.exec_() + +In our previous example, we didn't load a custom :class:`FitManager`. +Therefore, the fit widget automatically initialized a fit manager and +loaded the custom gaussian functions. + +This time, we initialized our own :class:`FitManager` and loaded our +own function, so only this function is presented as an option in the GUI. + +Our custom function does not provide an associated estimation function, so +the default estimation function of :class:`FitManager` was used. This +default estimation function returns an array of ones the same length as the +list of *parameter* names, and set all constraints to *FREE*. + + |imgFitWidget5| diff --git a/doc/source/Tutorials/fitconfig.rst b/doc/source/Tutorials/fitconfig.rst new file mode 100644 index 0000000..1b5e8ae --- /dev/null +++ b/doc/source/Tutorials/fitconfig.rst @@ -0,0 +1,207 @@ +Custom fit configuration widgets +================================ + +The *silx* fit widget allows users to add custom fit theories. +A fit theory consists of several components, such as the model function +to be fitted, an estimation function... + +One of these component is the optional custom configuration widget. +This is the dialog widget that is opened when a user clicks the *Configure* +button next to the drop-down list used to select the fit theory or the +background theory. + +This tutorial explains how to define your own fit configuration widget +for your custom fit theories. + +Prerequisites +-------------- + +This tutorial assumes that you are already familiar with +the standard features of :class:`silx.gui.fit.FitWidget`. +See the :ref:`fitwidget-tutorial` tutorial. + +You should also be familiar with adding custom fit theories. + +You will find documentation about these subjects by clicking the following links: + + - :class:`silx.math.fit.fitmanager` (methods `addtheory` and `addbgtheories`) + - :class:`silx.math.fit.fittheory` + - :class:`silx.math.fit.fittheories` + - :class:`silx.math.fit.bgtheories` + +The widgets we will create in this tutorial are based on the PyQt library. +Some knowledge of *PyQt* is desirable. + + +Basic concepts +-------------- + +Modal window +++++++++++++ + +A fit configuration widget should be a modal dialog window, so that +when a user opens the dialog to modify the configuration, the rest of +the program is frozen until all the configuration parameters are properly +defined. The program usually resumes when the user clicks the *Ok* or the +*Cancel* button in the dialog. + +The widget must implement a number of methods and attributes to be used as a +dialog by FitWidget: + + - standard *QDialog* methods: + + - :meth:`show`: should cause the widget to become visible to the + user) + - :meth:`exec_`: should run while the user is interacting with the + widget, interrupting the rest of the program. It should + typically end (*return*) when the user clicks an *OK* + or a *Cancel* button. + - :meth:`result`: must return ``True`` if the new configuration + is to be accepted (*OK* clicked) or ``False`` if it should be + rejected (*Cancel* clicked) + + - an additional *output* attribute, a dictionary storing configuration parameters + to be interpreted by the corresponding fit theory. + + - an optional *setDefault* method to initialize the + widget values with values in a dictionary passed as a parameter. + This will be executed first. + +The first 3 methods can be automatically defined by inheriting :class:`QDialog`. + +Associate a dialog to a theory +++++++++++++++++++++++++++++++ + +After defining a custom dialog widget, it must be initialized and associated +with a theory. + +A fit theory in :class:`FitWidget` is defined by a name. For example, +one of the default theories is named *"Gaussians"*. +So if you define a configuration dialog :class:`MyGaussianConfigWidget` to define +configuration parameters understood by this theory, you can associate it the following +way. + +.. code-block:: python + + fw = FitWidget() + my_config_widget = MyGaussianConfigWidget(parent=fw) + fw.associateConfigDialog(theory_name="Gaussians", + config_widget=my_config_widget) + + +Example +------- + +The following example defines a very basic configuration dialog widget +with a simple text entry in which the user can type in a floating point value. + +The value is simply saved in a dictionary attribute +:attr:`CustomConfigWidget.output`. *FitWidget* will look-up this dictionary +and pass it to the theory's custom configuration function, :func:`fitconfig`. +The configuration function essentially updates the :const:`CONFIG` dictionary +used by our fit function to scale the *y* values. + +.. code-block:: python + + from silx.gui import qt + from silx.gui.fit import FitWidget + from silx.math.fit.fittheory import FitTheory + from silx.math.fit.fitmanager import FitManager + + app = qt.QApplication([]) + + # default fit configuration + CONFIG = {"scale": 1.0} + + # define custom fit config dialog + class CustomConfigWidget(qt.QDialog): + def __init__(self): + qt.QDialog.__init__(self) + self.setModal(True) + self.scalingFactorEdit = qt.QLineEdit(self) + self.scalingFactorEdit.setToolTip( + "Enter the scaling factor" + ) + self.scalingFactorEdit.setValidator(qt.QDoubleValidator()) + + self.ok = qt.QPushButton("ok", self) + self.ok.clicked.connect(self.accept) + cancel = qt.QPushButton("cancel", self) + cancel.clicked.connect(self.reject) + + layout = qt.QVBoxLayout(self) + layout.addWidget(self.scalingFactorEdit) + layout.addWidget(self.ok) + layout.addWidget(cancel) + + self.old_scale = CONFIG["scale"] + self.output = {} + + def accept(self): + self.output["scale"] = float(self.scalingFactorEdit.text()) + qt.QDialog.accept(self) + + def reject(self): + self.output["scale"] = self.old_scale + qt.QDialog.reject(self) + + # our actual fit model function + def fitfun(x, a, b): + return CONFIG["scale"] * (a * x + b) + + # fit configuration + def fitconfig(scale=None, **kw): + """Update global config dict CONFIG""" + if scale is not None: + CONFIG["scale"] = scale + return CONFIG + + # synthetic test data a=2, b=3 + x = list(range(0, 100)) + y = [fitfun(x_, 2, 3) for x_ in x] + + # register our custom fit theory + fitmngr = FitManager() + fitmngr.setdata(x, y) + fitmngr.addtheory("scaled linear", + FitTheory( + function=fitfun, + parameters=["a", "b"], + configure=fitconfig)) + + # open a fitwidget and associate an instance of our custom + # configuration dialog to our custom theory + fw = FitWidget(fitmngr=fitmngr) + fw.associateConfigDialog("scaled linear", CustomConfigWidget()) + fw.show() + + app.exec_() + +.. |img0| image:: img/custom_config_scale1.0.png + :height: 300px + :align: middle + +.. |img1| image:: img/custom_config_scale2.1.png + :height: 300px + :align: middle + +.. |img2| image:: img/custom_config_scale0.5.png + :height: 300px + :align: middle + + +.. list-table:: + :widths: 1 4 + :header-rows: 1 + + * - Screenshot + - Description + * - |img0| + - If the default value of 1.0 is used, the fit finds *a=2* and *b=3* + as expected. + * - |img1| + - Setting a scaling factor of 2.1 causes the fit to find results that are + less than about half of the normal expected result. + * - |img2| + - A scaling factor of 0.5 causes the fit to find the values to be double + of the ones used for generating the synthetic data. diff --git a/doc/source/Tutorials/img/arraywidget3D_0.png b/doc/source/Tutorials/img/arraywidget3D_0.png new file mode 100644 index 0000000..0b654d7 Binary files /dev/null and b/doc/source/Tutorials/img/arraywidget3D_0.png differ diff --git a/doc/source/Tutorials/img/arraywidget3D_1.png b/doc/source/Tutorials/img/arraywidget3D_1.png new file mode 100644 index 0000000..fb55827 Binary files /dev/null and b/doc/source/Tutorials/img/arraywidget3D_1.png differ diff --git a/doc/source/Tutorials/img/arraywidget5D_0.png b/doc/source/Tutorials/img/arraywidget5D_0.png new file mode 100644 index 0000000..11a56de Binary files /dev/null and b/doc/source/Tutorials/img/arraywidget5D_0.png differ diff --git a/doc/source/Tutorials/img/arraywidget5D_1.png b/doc/source/Tutorials/img/arraywidget5D_1.png new file mode 100644 index 0000000..817b507 Binary files /dev/null and b/doc/source/Tutorials/img/arraywidget5D_1.png differ diff --git a/doc/source/Tutorials/img/custom_config_scale0.5.png b/doc/source/Tutorials/img/custom_config_scale0.5.png new file mode 100644 index 0000000..58d4665 Binary files /dev/null and b/doc/source/Tutorials/img/custom_config_scale0.5.png differ diff --git a/doc/source/Tutorials/img/custom_config_scale1.0.png b/doc/source/Tutorials/img/custom_config_scale1.0.png new file mode 100644 index 0000000..9dc9462 Binary files /dev/null and b/doc/source/Tutorials/img/custom_config_scale1.0.png differ diff --git a/doc/source/Tutorials/img/custom_config_scale2.1.png b/doc/source/Tutorials/img/custom_config_scale2.1.png new file mode 100644 index 0000000..d1d3663 Binary files /dev/null and b/doc/source/Tutorials/img/custom_config_scale2.1.png differ diff --git a/doc/source/Tutorials/img/fitwidget1.png b/doc/source/Tutorials/img/fitwidget1.png new file mode 100644 index 0000000..be371ac Binary files /dev/null and b/doc/source/Tutorials/img/fitwidget1.png differ diff --git a/doc/source/Tutorials/img/fitwidget2.png b/doc/source/Tutorials/img/fitwidget2.png new file mode 100644 index 0000000..5ee43bb Binary files /dev/null and b/doc/source/Tutorials/img/fitwidget2.png differ diff --git a/doc/source/Tutorials/img/fitwidget3.png b/doc/source/Tutorials/img/fitwidget3.png new file mode 100644 index 0000000..5da223c Binary files /dev/null and b/doc/source/Tutorials/img/fitwidget3.png differ diff --git a/doc/source/Tutorials/img/fitwidget4.png b/doc/source/Tutorials/img/fitwidget4.png new file mode 100644 index 0000000..feff3a2 Binary files /dev/null and b/doc/source/Tutorials/img/fitwidget4.png differ diff --git a/doc/source/Tutorials/img/fitwidget5.png b/doc/source/Tutorials/img/fitwidget5.png new file mode 100644 index 0000000..1e1278f Binary files /dev/null and b/doc/source/Tutorials/img/fitwidget5.png differ diff --git a/doc/source/Tutorials/img/stripbg_plot1.png b/doc/source/Tutorials/img/stripbg_plot1.png new file mode 100644 index 0000000..260890e Binary files /dev/null and b/doc/source/Tutorials/img/stripbg_plot1.png differ diff --git a/doc/source/Tutorials/img/stripbg_plot2.png b/doc/source/Tutorials/img/stripbg_plot2.png new file mode 100644 index 0000000..70e2ca2 Binary files /dev/null and b/doc/source/Tutorials/img/stripbg_plot2.png differ diff --git a/doc/source/Tutorials/specfile_to_hdf5.rst b/doc/source/Tutorials/specfile_to_hdf5.rst new file mode 100644 index 0000000..31f8383 --- /dev/null +++ b/doc/source/Tutorials/specfile_to_hdf5.rst @@ -0,0 +1,323 @@ + +SpecFile as HDF5 +================ + +Introduction to SPEC data files +------------------------------- + +SPEC data files are ASCII files. +They contain two general types of block of lines: + + - header lines starting with a ``#`` immediately followed by one or more characters + identifying the information that follows + - data lines + +Header lines +++++++++++++ + +There are two types of headers. The first type is the *file header*. File headers always start +with a ``#F`` line. +The metadata stored in a file header applies to all the content of the data file, until a +new file header is encountered. There can be more than one file header, but a file with +multiple headers can be treated as multiple SPEC files concatenated into a single one. +File headers are sometimes missing. + +A file header contains general information: + + - ``#F`` - file name + - ``#E`` - epoch + - ``#D`` - file time and date + - ``#C`` - First comment (SPEC title, SPEC user) + - ``#O`` - Motor names (separated by at least two blank spaces) + +The second type of header is the *scan header*. A scan header must start with a ``#S`` line +and must be preceded by an empty line. This also applies to files without file headers: in +such a case, the file must start with an empty line. +The metadata stored in scan headers applies to a single block of data lines. + +A scan header contains following information: + + - ``#S`` - Mandatory first line showing the scan number and the + command that was used to record the scan + - ``#D`` - scan time and date + - ``#Q`` - *H, K, L* values + - ``#P`` - Motor positions (corresponding motor names are in file header ``#O``) + - ``#N`` - Number of data columns in the following data block + - ``#L`` - Column labels (``#N`` labels separated by two blank spaces) + +Users can also define their own type of header lines in their macros. + +There can sometimes be a block of scan header lines after a data block, but before the ``#S`` of the next +scan. + +Data lines +++++++++++ + +Data blocks are structured as 2D arrays. Each line contains ``#N`` values, each value +corresponding to the label with the same position in the ``#L`` scan header line. +This implies that each column corresponds to one series of measurements. + +A column typically contains motor positions for a given positioner, a timestamp or the measurement +of a sensor. + +MCA data +++++++++ + +Newer SPEC files can also contain multi-channel analyser data, in between each *normal* data line. +A multichannel analyser records multiple values per single measurement. +This is typically a histogram of number of counts against channels (*MCA spectrum*), to analyze energy distribution +of a process. + +SPEC data files containing MCA data have additional scan header lines: + + - ``#@MCA %16C`` - a spectrum will usually extend for more than one line. + This indicates a number of 16 values per line. + - ``#@CHANN`` - contains 4 values: + + - the number of channels per spectrum + - the first channel number + - the last channel number + - the increment between two channel numbers (usually 1) + - ``#@CALIB`` - 3 polynomial calibration values a, b, c. ( i.e. energy = a + b * channel + c * channel ^ 2) + - ``#@CTIME`` - 3 values: preset time, live time, elapsed time + +The actual MCA data for a single spectrum usually spans over multiple lines. +A spectrum starts on a new line with a ``@A``, and when it span over multiple lines, all +lines except the last one end with a continuation character ``\``. + +Example of SPEC file +++++++++++++++++++++ + +Example of file header:: + + #F ./data/binary_mixtures_mca1.100211 + #E 1295362398 + #D Thu Feb 10 22:43:43 2011 + #C id10b User = opid10 + #O0 delta gamma omega theta mu sigma sigmat xt + #O1 zt zt1 thd chid rhod xd yd zd + #O2 att0 arcf zf PhiD phigH chigH ygH + #O3 zgH phigV chigV xgV ygV zgV gslithg gslitho + #O4 gslitvo gslitvg slit1T slit1B slit1F slit1R slit1hg slit1ho + #O5 slit1vg slit1vo s0T s0B s0R s0F + #O6 s0hg s0ho s0vg s0vo TRT + #O7 pi trough hv1 mpxthl apdwin apdthl apdhv xcrl2 + #O8 thcrl2 zcrl2 picou picod vdrift vmulti vglo vghi + #O9 rien + +Example of scan and data block, without MCA:: + + #S 30 ascan tz3 29.35 29.75 100 0.5 + #D Sat Oct 31 15:43:21 1998 + #T 0.5 (Seconds) + #G0 0 + #G1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + #G2 0 + #Q + #P0 40.135381 40.262001 65.6 70 35 -1.83 0 -36.1 + #P1 0 0 -1.98 0 0 35.6 86.2 -29.5 + #P2 3.0688882 24.893749 295.98749 28 -27.249938 + #N 22 + #L TZ3 Epoch Seconds If2 If3 If5 If6 If7 If8 I0 It ItdI0 If1dI0 If2dI0 If3dI0 If4dI0 If5dI0 If6dI0 If7dI0 If8dI0 If1 If4 + 29.35 45246 0.000264 478 302 206 201 209 264 177860 646 0.00363207 0.00468346 0.00268751 0.00169796 0.00146745 0.00115821 0.0011301 0.00117508 0.00148431 833 261 + 29.353976 45249 0.000295 549 330 219 208 227 295 178021 684 0.00384224 0.00537577 0.00308391 0.00185371 0.00158408 0.00123019 0.0011684 0.00127513 0.00165711 957 282 + 29.357952 45251 0.000313 604 368 231 215 229 313 178166 686 0.00385034 0.00603931 0.0033901 0.00206549 0.00166698 0.00129654 0.00120674 0.00128532 0.00175679 1076 297 + 29.362028 45253 0.000333 671 390 237 226 236 333 178387 672 0.00376709 0.00683346 0.00376148 0.00218626 0.00176582 0.00132857 0.00126691 0.00132297 0.00186673 1219 315 + 29.366004 45256 0.000343 734 419 248 229 236 343 178082 664 0.00372862 0.00765939 0.0041217 0.00235285 0.00185308 0.00139262 0.00128592 0.00132523 0.00192608 1364 330 + 29.36998 45258 0.00036 847 448 254 229 248 360 178342 668 0.00374561 0.00857342 0.0047493 0.00251203 0.00194009 0.00142423 0.00128405 0.00139059 0.00201859 1529 346 + +Synthetic example of file with 3 scans. The last scan includes data of 3 multichannel analysers, sharing the +same MCA header. + +:: + + #F /tmp/sf.dat + #E 1455180875 + #D Thu Feb 11 09:54:35 2016 + #C imaging User = opid17 + #O0 Pslit HGap MRTSlit UP MRTSlit DOWN + #O1 Sslit1 VOff Sslit1 HOff Sslit1 VGap + #o0 pshg mrtu mrtd + #o2 ss1vo ss1ho ss1vg + + #S 1 ascan ss1vo -4.55687 -0.556875 40 0.2 + #D Thu Feb 11 09:55:20 2016 + #T 0.2 (Seconds) + #P0 180.005 -0.66875 0.87125 + #P1 14.74255 16.197579 12.238283 + #N 3 + #L MRTSlit UP second column 3rd_col + -1.23 5.89 8 + 8.478100E+01 5 1.56 + 3.14 2.73 -3.14 + 1.2 2.3 3.4 + + #S 25 ascan c3th 1.33245 1.52245 40 0.15 + #D Sat 2015/03/14 03:53:50 + #P0 80.005 -1.66875 1.87125 + #P1 4.74255 6.197579 2.238283 + #N 4 + #L column0 column1 col2 col3 + 0.0 0.1 0.2 0.3 + 1.0 1.1 1.2 1.3 + 2.0 2.1 2.2 2.3 + 3.0 3.1 3.2 3.3 + + #S 1 aaaaaa + #D Thu Feb 11 10:00:32 2016 + #@MCA %16C + #@CHANN 20 0 19 1 + #@CALIB 1.2 2.3 3.4 + #@CTIME 123.4 234.5 345.6 + #N 2 + #L uno duo + 1 2 + @A 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15\ + 16 17 18 19 + 3 4 + @A 0 0 2 4 15 10 5 1 0 0 0 0 1 0 0 0\ + 0 0 0 0 + 5 6 + @A 0 0 0 0 5 7 2 0 0 0 0 0 1 0 0 0\ + 0 0 0 1 + +Reading a SpecFile as an HDF5 file +---------------------------------- + +Introduction to the spech5 module ++++++++++++++++++++++++++++++++++ + +The *silx* module :mod:`silx.io.spech5` can be used to expose SPEC files in a hierarchical tree structure +and access them through an API that mimics the *h5py* Python library used to read HDF5 files. + +The structure exposed is as follows:: + + / + 1.1/ + title = "…" + start_time = "…" + instrument/ + specfile/ + file_header = ["…", "…", …] + scan_header = ["…", "…", …] + positioners/ + motor_name = value + … + mca_0/ + data = … + calibration = … + channels = … + preset_time = … + elapsed_time = … + live_time = … + + mca_1/ + … + … + measurement/ + colname0 = … + colname1 = … + … + mca_0/ + data -> /1.1/instrument/mca_0/data + info -> /1.1/instrument/mca_0/ + … + 2.1/ + … + +Scans appear as *Groups* at the root level. The name of a scan group is +made of two numbers, the first one being the *scan number* from the ``#S`` +header line, and the second one being the *scan order*. +If a scan number appears multiple times in a SPEC file, the scan order is incremented. +For examples, the scan *3.2* designates the second occurence of scan number 3 in a given file. + +Data is stored in the ``measurement`` subgroup, one dataset per column. The dataset name +is the column label as it appears on the ``#L`` header line. + +The ``instrument`` subgroup contains following subgroups: + + - ``specfile`` - contains two datasets, ``file_header`` and ``scan_header``, + containing all header lines as a long string. Lines are delimited by the ``\n`` character. + - ``positioners`` - contains one dataset per motor (positioner), containing + either the single motor position from the ``#P`` header line, or a complete 1D array + of positions if the motor names corresponds to a data column (i.e. if the motor name + from the ``#O`` header line is identical to a label on the ``#L`` header line) + - one subgroup per MCA analyser/device containing a 2D ``data`` array with all spectra + recorded by this analyser, as well as datasets for the various MCA metadata + (``#@`` header lines). The first dimension of the ``data`` array corresponds to the number + of points and the second one to the spectrum length. + + +In addition the the data columns, this group contains one subgroup per MCA analyser/device +with links to the data already contained in ``instrument/mca_...`` + +spech5 examples ++++++++++++++++ + +Accessing groups and datasets: + +.. code-block:: python + + from silx.io.spech5 import SpecH5 + + # Open a SpecFile + sfh5 = SpecH5("test.dat") + + # using SpecH5 as a regular group to access scans + scan1group = sfh5["1.1"] # This retrieves scan 1.1 + scan1group = sfh5[0] # This retrieves the first scan irrespectively of its number. + instrument_group = scan1group["instrument"] + + # alternative: full path access + measurement_group = sfh5["/1.1/measurement"] + + # accessing a scan data column by name as a 1D numpy array + data_array = measurement_group["Pslit HGap"] + + # accessing all mca-spectra for one MCA device as a 2D array + mca_0_spectra = measurement_group["mca_0/data"] + + +Files and groups can be treated as iterators, which allows looping through them. + +.. code-block:: python + + # get all column names (labels) in all scans in a file + for scan_group in SpecH5("test.dat"): + dataset_names = [item.name in scan_group["measurement"] if not + item.name.startswith("mca")] + print("Found labels in scan " + scan_group.name + " :") + print(", ".join(dataset_names)) + +Converting SPEC data to HDF5 +++++++++++++++++++++++++++++ + +The *silx* module :mod:`silx.io.spectoh5` can be used to convert a SPEC file into a +HDF5 file with the same structure as the one exposed by the :mod:`spech5` module. + +.. code-block:: python + + from silx.io.spectoh5 import convert + + convert("/home/pierre/myspecfile.dat", "myfile.h5") + + +You can then read the file with any HDF5 reader. + + +In addition to the function :func:`silx.io.spectoh5.convert`, which is simplified +on purpose, you can use the more flexible :func:`silx.io.spectoh5.write_spec_to_h5`. + +This way, you can choose to write scans into a specific HDF5 group in the output directory. +You can also decide whether you want to overwrite an existing file, or append data to it. +You can specify whether existing data with the same name as input data should be overwritten +or ignored. + +This allows you to repeatedly transfer new content of a SPEC file to an existing HDF5 file, in between +two scans. + +The following script is an example of a command line interface to :func:`write_spec_to_h5`. + +.. literalinclude:: ../../../examples/spectoh5.py + :lines: 42- + -- cgit v1.2.3