Source code for libertem.analysis.radialfourier

import logging
import inspect
from functools import partial

import numpy as np
import sparse

from libertem import masks
from libertem.common.math import prod
from .base import AnalysisResult, AnalysisResultSet
from .masks import BaseMasksAnalysis
from .helper import GeneratorHelper

log = logging.getLogger(__name__)

class RadialTemplate(GeneratorHelper):

    short_name = "radial"
    api = "create_radial_fourier_analysis"
    temp = GeneratorHelper.temp_analysis
    temp_analysis = temp + ["print(radial_result)"]

    def __init__(self, params):
        self.params = params

    def get_dependency(self):
        return [
            "import as cm",
            "from empyre.vis.colors import ColormapCubehelix, ColormapPerception"

    def get_docs(self):
        title = "Radial Fourier Analysis"
        from libertem.api import Context
        docs_rst = inspect.getdoc(Context.create_radial_fourier_analysis)
        docs = self.format_docs(title, docs_rst)
        return docs

    def convert_params(self):
        params = ['dataset=ds']
        for k in ['cx', 'cy', 'ri', 'ro', 'n_bins', 'max_order']:
        return ', '.join(params)

    def get_plot(self):
        cells = []
            "fig, axes = plt.subplots()",
            "plt.imshow(radial_result.dominant_0, cmap=cm.tab20, vmin=0, vmax=20)",
            "fig, axes = plt.subplots()",
            "imag = radial_result.complex_0_1.raw_data.imag",
            "real = radial_result.complex_0_1.raw_data.real",
            "ch = ColormapCubehelix(start=1, rot=1, minLight=0.5, maxLight=0.5, sat=2)",
            "fig, axes = plt.subplots()",
            "plt.imshow(ch.rgb_from_vector((real, imag, 0)))",
            "fig, axes = plt.subplots()",
            'plt.imshow(radial_result.phase_0_1.raw_data, cmap=ColormapPerception())'
        return ['\n'.join(cell) for cell in cells]

    def get_save(self):
        save = []
        channels = ["absolute_0_0", "absolute_0_1"]
        for channel in channels:
            result = f"radial_result['{channel}'].raw_data"
            save.append(f"'radial_result_{channel}.npy', {result})")
        return '\n'.join(save)

[docs]class RadialFourierResultSet(AnalysisResultSet): """ Result set of a :class:`RadialFourierAnalysis` Running a :class:`RadialFourierAnalysis` via :meth:`` on a dataset returns an instance of this class. It contains the Fourier coefficients for each bin. See :meth:`libertem.api.Context.create_radial_fourier_analysis` for available parameters and :ref:`radialfourier app` for a description of the application! .. versionadded:: 0.3.0 Attributes ---------- dominant_0, absolute_0_0, absolute_0_1, ..., absolute_0_<max_order>,\ phase_0_0, ..., phase_0_<max_order>,\ complex_0_0, ..., complex_0_<max_order>;\ dominant_1, absolute_1_0, ..., complex_1_<max_order>;\ dominant_<nbins-1>, ..., complex_<nbins-1>_<max_order> : libertem.analysis.base.AnalysisResult Results for each bin: dominant Fourier coefficient (indicates symmetry), absolute values of each Fourier coefficient, phase values of each Fourier coefficient, complex values of each Fourier coefficient. The results have the shape of the navigation dimension. raw_results : numpy.ndarray Complex numbers, shape (<n_bins>, <max_order + 1>, \\*(<ds.shape.nav>)) """ pass
def radial_mask_factory(detector_y, detector_x, cx, cy, ri, ro, n_bins, max_order, use_sparse): def stack(): rings = masks.radial_bins( centerX=cx, centerY=cy, imageSizeX=detector_x, imageSizeY=detector_y, radius=ro, radius_inner=ri, n_bins=n_bins, use_sparse=use_sparse, dtype=np.complex64 ) orders = np.arange(max_order + 1) r, phi = masks.polar_map( centerX=cx, centerY=cy, imageSizeX=detector_x, imageSizeY=detector_y ) modulator = np.exp(phi * orders[:, np.newaxis, np.newaxis] * 1j) if use_sparse: rings = rings.reshape((rings.shape[0], 1, *rings.shape[1:])) ring_stack = [rings] * len(orders) ring_stack = sparse.concatenate(ring_stack, axis=1) ring_stack *= modulator else: ring_stack = rings[:, np.newaxis, ...] * modulator return ring_stack.reshape((-1, detector_y, detector_x)) return stack class RadialFourierAnalysis(BaseMasksAnalysis, id_="RADIAL_FOURIER"): ''' The Radial Fourier Analysis can be used to characterize atomic ordering in materials, in particular for low intensities where Fluctualtion EM :cite:`Gibson1997` has a hard time to distinguish speckle from shot noise. Reference :cite:`6980942` describes a previous application of this method to characterize features in images. This analysis doesn't use fast Fourier transforms, but calculates the Fourier coefficients using sparse matrices in a dot product following the `definition of Fourier series <>`_. See :meth:`libertem.api.Context.create_radial_fourier_analysis` for available parameters and :ref:`radialfourier app` for a description of the application! ''' TYPE = 'UDF' def get_udf_results(self, udf_results, roi, damage): ''' The AnalysisResults are calculated lazily in this function to reduce overhead. ''' shape = tuple(self.dataset.shape.nav) # NOTE: transposed for historical reasons udf_results = udf_results['intensity'].data.reshape((prod(shape), -1)).T orders = self.parameters['max_order'] + 1 n_bins = self.parameters['n_bins'] udf_results = udf_results.reshape((n_bins, orders, *shape)) def resultlist(): from libertem.viz import CMAP_CIRCULAR_DEFAULT, visualize_simple, cmaps import as cm sets = [] absolute = np.absolute(udf_results) normal = np.maximum(1, absolute[:, 0]) # New local variable since this is a closure over damage dam = damage & np.all(np.isfinite(absolute), axis=(0, 1)) normalized = absolute[:, 1:, ...] / normal[:, np.newaxis, ...] min_absolute = np.min(normalized[..., dam]) max_absolute = np.max(normalized[..., dam]) angle = np.angle(udf_results) threshold_map = absolute[:, 1:, ...].reshape((n_bins, -1)).max(axis=1) * 0.2 below_threshold = np.all( absolute[:, 1:, ...] < threshold_map[:, np.newaxis, np.newaxis, np.newaxis], axis=1 ) dominant = np.argmax(absolute[:, 1:], axis=1) + 1 dominant[below_threshold] = 0 for b in range(n_bins): sets.append( AnalysisResult( raw_data=dominant[b], visualized=partial( visualize_simple, dominant[b], colormap=cm.tab20, vmin=0, vmax=20 ), key="dominant_%s" % b, title="dominant order of bin %s" % b, desc="Dominant Fourier component", ) ) sets.append( AnalysisResult( raw_data=absolute[b, 0], visualized=partial(visualize_simple, absolute[b, 0], damage=dam), key=f"absolute_{b}_{0}", title=f"absolute of bin {b} order {0}", desc="Absolute value of Fourier component", ) ) for o in range(1, orders): sets.append( AnalysisResult( raw_data=absolute[b, o], visualized=partial(visualize_simple, absolute[b, o] / normal[b], vmin=min_absolute, vmax=max_absolute, damage=dam ), key=f"absolute_{b}_{o}", title=f"absolute of bin {b} order {o}", desc="Absolute value of Fourier component", ) ) for b in range(n_bins): for o in range(orders): sets.append( AnalysisResult( raw_data=angle[b, o], visualized=partial(visualize_simple, angle[b, o], colormap=cmaps['perception_circular'], damage=dam ), key=f"phase_{b}_{o}", title=f"phase of bin {b} order {o}", desc="Phase of Fourier component", ) ) for b in range(n_bins): data = udf_results[b, 0] f = partial( CMAP_CIRCULAR_DEFAULT.rgb_from_vector, (data.real, data.imag, 0), vmax=np.max(np.abs(data[..., dam])) ) sets.append( AnalysisResult( raw_data=data, visualized=f, key=f"complex_{b}_{0}", title=f"bin {b} order {0}", desc="Fourier component", ) ) for o in range(1, orders): data = udf_results[b, o] / normal[b] f = partial( CMAP_CIRCULAR_DEFAULT.rgb_from_vector, (data.real, data.imag, 0), vmax=max_absolute ) sets.append( AnalysisResult( raw_data=data, visualized=f, key=f"complex_{b}_{o}", title=f"bin {b} order {o}", desc="Fourier component", ) ) return sets return RadialFourierResultSet(resultlist, raw_results=udf_results) def get_mask_factories(self): if self.dataset.shape.sig.dims != 2: raise ValueError("can only handle 2D signals currently") (detector_y, detector_x) = self.dataset.shape.sig p = self.parameters return radial_mask_factory( detector_y=detector_y, detector_x=detector_x, cx=p['cx'], cy=p['cy'], ri=p['ri'], ro=p['ro'], n_bins=p['n_bins'], max_order=p['max_order'], use_sparse=p['use_sparse'], ) def get_parameters(self, parameters): (detector_y, detector_x) = self.dataset.shape.sig cx = parameters.get('cx', detector_x / 2) cy = parameters.get('cy', detector_y / 2) ri = parameters.get('ri', 0) ro = parameters.get( 'ro', masks.bounding_radius(cx, cy, detector_x, detector_y) ) n_bins = parameters.get('n_bins', 1) max_order = parameters.get('max_order', 24) mask_count = n_bins * (max_order + 1) bin_width = (ro - ri) / n_bins bin_area = np.pi * ro**2 - np.pi * (ro - bin_width)**2 stack_size = mask_count * detector_y * detector_x * 8 default = 'scipy.sparse' # If the mask stack comfortably fits the L3 cache # FIXME more testing for optimum backend if stack_size < 2**18: default = False # Masks are actually dense elif bin_area / (detector_x * detector_y) > 0.05 and n_bins < 10: default = False use_sparse = parameters.get('use_sparse', default) return { 'cx': cx, 'cy': cy, 'ri': ri, 'ro': ro, 'n_bins': n_bins, 'max_order': max_order, 'use_sparse': use_sparse, 'mask_count': mask_count, 'mask_dtype': np.complex64, } @classmethod def get_template_helper(cls): return RadialTemplate