Commit 6ac7d907 authored by payno's avatar payno

[core mapping] rework Intensity mapping too fit several dimension

parent d5e1cbbc
...@@ -325,21 +325,21 @@ class Experiment(object): ...@@ -325,21 +325,21 @@ class Experiment(object):
self.__dims.clear() self.__dims.clear()
_fail = False, None # register fail status and latest reason _fail = False, None # register fail status and latest reason
for axis, dim in dims.items(): for axis, dim in dims.items():
if dim.size is None: try:
assert isinstance(dim, _Dim) unique_dim_value = metadatautils.getUnique(self,
# try to deduce dimension size kind=dim.kind,
try: key=dim.name,
unique_dim_value = metadatautils.getUnique(self, relative_prev_val=dim.relative_prev_val,
kind=dim.kind, cycle_length=dim.size,
key=dim.name, axis=axis)
relative_prev_val=dim.relative_prev_val, except Exception as e:
cycle_length=dim.size, if _fail[0] is False:
axis=axis) _fail = True, e
except Exception as e: else:
if _fail[0] is False: if dim.size is None:
_fail = True, e
else:
dim._setSize(size=len(unique_dim_value)) dim._setSize(size=len(unique_dim_value))
dim._setUniqueValues(unique_dim_value)
self.__dims.add_dim(axis=axis, dim=dim) self.__dims.add_dim(axis=axis, dim=dim)
if _fail[0] is True: if _fail[0] is True:
...@@ -390,6 +390,10 @@ class AcquisitionDims(object): ...@@ -390,6 +390,10 @@ class AcquisitionDims(object):
kind=self.__dims[axis].kind, kind=self.__dims[axis].kind,
size=size) size=size)
def __iter__(self):
for iAxis, dim in self.__dims.items():
yield (iAxis, dim)
class _Dim(object): class _Dim(object):
def __init__(self, kind, name, size=None, relative_prev_val=False): def __init__(self, kind, name, size=None, relative_prev_val=False):
...@@ -418,6 +422,8 @@ class _Dim(object): ...@@ -418,6 +422,8 @@ class _Dim(object):
self.__name = name self.__name = name
self._size = size self._size = size
self.__relative_prev_val = relative_prev_val self.__relative_prev_val = relative_prev_val
self.__unique_values = []
"""Ordered values through the dimension"""
@property @property
def kind(self): def kind(self):
...@@ -450,3 +456,16 @@ class _Dim(object): ...@@ -450,3 +456,16 @@ class _Dim(object):
(see :class:`_DimensionItem`) (see :class:`_DimensionItem`)
""" """
self._size = size self._size = size
@property
def unique_values(self):
return self.__unique_values
def _setUniqueValues(self, values):
if len(values) != self.size:
_logger.warning('Unique values set for %s, is incoherent with size' % self)
raise ValueError('Unique values set for %s, is incoherent with size' % self)
self.__unique_values = values
def __str__(self):
return " ".join((str(self.kind), str(self.name), 'size:', str(self.size)))
...@@ -62,7 +62,7 @@ class COM(AdditiveOperation): ...@@ -62,7 +62,7 @@ class COM(AdditiveOperation):
assert self._map.ndim is 1 assert self._map.ndim is 1
self.__maps = [] self.__maps = []
for dim in self._map.dim: for dim in self._map.dims:
_map1 = dim[:, :, 0] _map1 = dim[:, :, 0]
_map1 = _map1 - numpy.nan_to_num(_map1.mean()) _map1 = _map1 - numpy.nan_to_num(_map1.mean())
......
...@@ -29,12 +29,17 @@ __license__ = "MIT" ...@@ -29,12 +29,17 @@ __license__ = "MIT"
__date__ = "15/10/2018" __date__ = "15/10/2018"
from . import AdditiveOperation from . import AdditiveOperation
from collections import namedtuple
import numpy import numpy
import logging import logging
_logger = logging.getLogger(__file__) _logger = logging.getLogger(__file__)
_IMap = namedtuple('_IMap', ['name', 'kind', 'mean', 'variance', 'skewness',
'kurtosis'])
class _MappingBase(AdditiveOperation): class _MappingBase(AdditiveOperation):
""" """
Base class used for mapping Base class used for mapping
...@@ -73,7 +78,8 @@ class IntensityMapping(_MappingBase): ...@@ -73,7 +78,8 @@ class IntensityMapping(_MappingBase):
self._threshold = threshold self._threshold = threshold
self.__intensity_map = None self.__intensity_map = None
self.__dim = [] self.__dim = {}
"""Associate to each dimension the _IMap"""
def key(self): def key(self):
return self.KEY return self.KEY
...@@ -85,18 +91,15 @@ class IntensityMapping(_MappingBase): ...@@ -85,18 +91,15 @@ class IntensityMapping(_MappingBase):
_logger.warning('Experiment as no data defined, unable to apply' _logger.warning('Experiment as no data defined, unable to apply'
'mapping') 'mapping')
return None return None
elif self._experiment.geometry is None:
_logger.warning('No geometry defined to apply mapping')
return None
else: else:
assert self.data.ndim is 3 assert self.data.ndim > 2
# TODO: avoid recomputing several time the intensity map # TODO: avoid recomputing several time the intensity map
self._resetIntensityMap() self._resetIntensityMap()
self._resetDim() self._resetDim()
for x in range(self.data.shape[2]): for x in range(self.data.shape[-1]):
self.updateProgress(int(x / self.data.shape[2] * 100.0)) self.updateProgress(int(x / self.data.shape[2] * 100.0))
for y in range(self.data.shape[1]): for y in range(self.data.shape[-2]):
# TODO: do this operation using a grid or something close # TODO: do this operation using a grid or something close
# but would have to find a way to write at the same position... # but would have to find a way to write at the same position...
reciprocalVol = numpy.squeeze(self.data[:, y, x]) reciprocalVol = numpy.squeeze(self.data[:, y, x])
...@@ -106,25 +109,21 @@ class IntensityMapping(_MappingBase): ...@@ -106,25 +109,21 @@ class IntensityMapping(_MappingBase):
self.__intensity_map[y, x] = maxIntensity self.__intensity_map[y, x] = maxIntensity
# Here: TODO, later when can have several dimension (motors scanned/rockers) # Here: TODO, later when can have several dimension (motors scanned/rockers)
# do it for all dims # do it for all dims
for iDim in range(self._experiment.ndim): for axis, dim in self._experiment.dims:
intensity = numpy.squeeze(reciprocalVol.sum()) intensity = numpy.squeeze(reciprocalVol.sum())
# TODO: move calculation of angle_mean, angle_variance, # TODO: move calculation of angle_mean, angle_variance,
# angle_skewness, angle_kurtosis in :class:`Experiment` dim_unique_values = self.expriment.get_unique_values(axis)
angles = self._experiment.angles
# TODO: why mean is obtained by the sum ? # TODO: why mean is obtained by the sum ?
# angles_mean = (intensity * angles).mean() # angles_mean = (intensity * angles).mean()
angles_mean = (intensity * angles).sum() _mean = (intensity * dim_unique_values).sum()
angles_variance = ( diff = (dim_unique_values - _mean)
intensity * (angles - angles_mean) ** 2).sum() _variance = (intensity * diff ** 2).sum()
angles_skewness = ( _skewness = (intensity * diff ** 3).sum()
intensity * (angles - angles_mean) ** 3).sum() _kurtosis = (intensity * diff ** 4).sum()
angles_kurtosis = ( self.dim[axis].mean[y, x] = _mean
intensity * (angles - angles_mean) ** 4).sum() self.dim[axis].variance[y, x] = _variance
# TODO: pb, this is overwritten at each image if maxIntensity > threshold ??? self.dim[axis].skewness[y, x] = _skewness
self.dim[iDim][y, x, 0] = angles_mean self.dim[axis].kurtosis[y, x] = _kurtosis
self.dim[iDim][y, x, 1] = angles_variance
self.dim[iDim][y, x, 2] = angles_skewness
self.dim[iDim][y, x, 3] = angles_kurtosis
self.registerOperation() self.registerOperation()
return self.dim return self.dim
...@@ -133,9 +132,16 @@ class IntensityMapping(_MappingBase): ...@@ -133,9 +132,16 @@ class IntensityMapping(_MappingBase):
self.__intensity_map = numpy.zeros(self.data.shape[1:]) self.__intensity_map = numpy.zeros(self.data.shape[1:])
def _resetDim(self): def _resetDim(self):
self.__dim = [] self.__dim = {}
for dim in range(self._experiment.ndim): for axis, dim in self._experiment.dims:
self.__dim.append(numpy.zeros((*(self.data.shape[1:]), 4))) _map = _IMap(
name=dim.name,
kind=dim.kind,
mean=numpy.zeros((self.data.shape[-2:])),
variance=numpy.zeros((self.data.shape[-2:])),
skewness=numpy.zeros((self.data.shape[-2:])),
kurtosis=numpy.zeros((self.data.shape[-2:])))
self.__dim[axis] = _map
@property @property
def intensity_map(self): def intensity_map(self):
...@@ -192,7 +198,7 @@ class GradientRemoval(_MappingBase): ...@@ -192,7 +198,7 @@ class GradientRemoval(_MappingBase):
# TODO: if necessary compute mean... # TODO: if necessary compute mean...
for iDim in range(mapping.ndim): for iDim in range(mapping.ndim):
gradient = mapping.dim[iDim][...] gradient = mapping.dims[iDim][...]
# TODO: for now gradient corrcetion is only apply on mean, why ? # TODO: for now gradient corrcetion is only apply on mean, why ?
# TODO: this should be in enumerate(Mode) but intensity is stored # TODO: this should be in enumerate(Mode) but intensity is stored
......
...@@ -32,10 +32,10 @@ MEAN = 'mean' ...@@ -32,10 +32,10 @@ MEAN = 'mean'
VARIANCE = 'variance' VARIANCE = 'variance'
SKENESS = 'skeness' SKEWNESS = 'skewness'
KURTOSIS = 'kurtosis' KURTOSIS = 'kurtosis'
INTENSITY = 'intensity' INTENSITY = 'intensity'
MODES = (MEAN, VARIANCE, SKENESS, KURTOSIS, INTENSITY) MODES = (MEAN, VARIANCE, SKEWNESS, KURTOSIS, INTENSITY)
...@@ -34,13 +34,13 @@ from id06workflow.core.geometry.TwoThetaGeometry import TwoThetaGeometry ...@@ -34,13 +34,13 @@ from id06workflow.core.geometry.TwoThetaGeometry import TwoThetaGeometry
class TestTwoThetaExpSetup(unittest.TestCase): class TestTwoThetaExpSetup(unittest.TestCase):
""" """
Test that RefCopy process is correct Test the :class:`TwoThetaGeometry`
""" """
def testCreation(self): def testCreation(self):
""" """
Some stupid test to make sur TwoThetaExpSetup is correctly Some stupid test to make sur TwoThetaExpSetup is correctly
instanciated instantiated
""" """
exp_setup = TwoThetaGeometry(orientation=TwoThetaGeometry.VERTICAL) exp_setup = TwoThetaGeometry(orientation=TwoThetaGeometry.VERTICAL)
self.assertTrue(exp_setup.orientation == TwoThetaGeometry.VERTICAL) self.assertTrue(exp_setup.orientation == TwoThetaGeometry.VERTICAL)
......
# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/
__authors__ = ["H. Payno"]
__license__ = "MIT"
__date__ = "03/10/2018"
import unittest
import os
from id06workflow.core.experiment.operation.mapping import IntensityMapping
from id06workflow.core.experiment import Experiment, Dataset, POSITIONER_METADATA, _Dim
from id06workflow.core.geometry.TwoThetaGeometry import TwoThetaGeometry
@unittest.skipIf(os.path.exists('/nobackup/linazimov/payno/datasets/id06/strain_scan') is False, reason='data files not available')
class TestMappingOperation(unittest.TestCase):
"""
Test Mapping operation.
"""
def setUp(self):
unittest.TestCase.setUp(self)
# define experiment data
root_folder = '/nobackup/linazimov/payno/datasets/id06/strain_scan'
data_file_pattern = os.path.join(root_folder, 'reduced_strain/strain_0000.edf')
assert os.path.exists(data_file_pattern)
data_bb_files = []
bb_folder = os.path.join(root_folder, 'bg_ff_5s_1x1')
for _file in os.listdir(bb_folder):
data_bb_files.append(os.path.join(bb_folder, _file))
self.dataset = Dataset(data_files_pattern=data_file_pattern,
ff_files=data_bb_files)
# define geometry
geometry = TwoThetaGeometry(
twotheta=0.03, # for now those are defined in degree but should
# be defined in radians
xmag=1.0, # what is the unit of this ?
xpixelsize=1.0,
ypixelsize=1.0,
orientation=TwoThetaGeometry.VERTICAL)
# define experiment
self.experiment = Experiment(dataset=self.dataset, geometry=geometry)
# define dimensions
dim1 = _Dim(kind=POSITIONER_METADATA, name='diffry',
relative_prev_val=True, size=31)
dim2 = _Dim(kind=POSITIONER_METADATA, name='obpitch')
self.experiment.set_dims(dims={0: dim1, 1: dim2})
def testIntensityMapping(self):
"""
Test IntensityMapping operation
"""
operation = IntensityMapping(experiment=self.experiment)
res = operation.compute()
self.assertTrue(len(res) is 2)
self.assertTrue(res[0].name == 'diffry')
self.assertTrue(res[0].kind == POSITIONER_METADATA)
self.assertTrue(res[1].name == 'obpitch')
self.assertTrue(res[1].kind == POSITIONER_METADATA)
self.assertTrue(res[0].mean.shape == (2048, 2048))
self.assertTrue(res[0].variance.shape == (2048, 2048))
self.assertTrue(res[0].skewness.shape == (2048, 2048))
self.assertTrue(res[0].kurtosis.shape == (2048, 2048))
self.assertTrue(res[1].mean.shape == (2048, 2048))
def suite():
test_suite = unittest.TestSuite()
for ui in (TestMappingOperation, ):
test_suite.addTest(unittest.defaultTestLoader.loadTestsFromTestCase(ui))
return test_suite
if __name__ == '__main__':
unittest.main(defaultTest="suite")
...@@ -65,7 +65,7 @@ class ComWidget(qt.QWidget): ...@@ -65,7 +65,7 @@ class ComWidget(qt.QWidget):
dim = self._control.getDim() dim = self._control.getDim()
_iMap = self._control.getMapType() _iMap = self._control.getMapType()
if dim >= 0: if dim >= 0:
# TODO: maps[dim] should be also a list with key values: 'value', # TODO: maps[dims] should be also a list with key values: 'value',
# 'mean', 'gradients' ... # 'mean', 'gradients' ...
if len(self._operation.maps[dim]) <= _iMap: if len(self._operation.maps[dim]) <= _iMap:
_logger.error('map of index %s has not been created' % _iMap) _logger.error('map of index %s has not been created' % _iMap)
...@@ -94,7 +94,7 @@ class _ComControlWidget(qt.QWidget): ...@@ -94,7 +94,7 @@ class _ComControlWidget(qt.QWidget):
def __init__(self, parent=None): def __init__(self, parent=None):
qt.QWidget.__init__(self, parent) qt.QWidget.__init__(self, parent)
self.setLayout(qt.QHBoxLayout()) self.setLayout(qt.QHBoxLayout())
self.layout().addWidget(qt.QLabel('dim:', parent=self)) self.layout().addWidget(qt.QLabel('dims:', parent=self))
self._dim = qt.QComboBox(parent=self) self._dim = qt.QComboBox(parent=self)
self.layout().addWidget(self._dim) self.layout().addWidget(self._dim)
self.layout().addWidget(qt.QLabel('map type:', parent=self)) self.layout().addWidget(qt.QLabel('map type:', parent=self))
......
...@@ -31,7 +31,7 @@ __date__ = "10/10/2018" ...@@ -31,7 +31,7 @@ __date__ = "10/10/2018"
from silx.gui import qt from silx.gui import qt
from silx.gui.plot import Plot2D from silx.gui.plot import Plot2D
from id06workflow.gui import icons from id06workflow.gui import icons
from id06workflow.core.mapping import MEAN, VARIANCE, SKENESS, KURTOSIS, MODES, INTENSITY from id06workflow.core.mapping import MEAN, VARIANCE, SKEWNESS, KURTOSIS, MODES, INTENSITY
from id06workflow.core.experiment.operation.mapping import _MappingBase from id06workflow.core.experiment.operation.mapping import _MappingBase
from id06workflow.gui.settings import DEFAULT_COLORMAP from id06workflow.gui.settings import DEFAULT_COLORMAP
from collections import OrderedDict from collections import OrderedDict
...@@ -48,6 +48,7 @@ class MappingPlot(qt.QMainWindow): ...@@ -48,6 +48,7 @@ class MappingPlot(qt.QMainWindow):
""" """
def __init__(self, parent=None): def __init__(self, parent=None):
qt.QMainWindow.__init__(self, parent) qt.QMainWindow.__init__(self, parent)
self._operation = None
self.setWindowFlags(qt.Qt.Widget) self.setWindowFlags(qt.Qt.Widget)
self._plot = Plot2D(parent=self) self._plot = Plot2D(parent=self)
...@@ -57,7 +58,7 @@ class MappingPlot(qt.QMainWindow): ...@@ -57,7 +58,7 @@ class MappingPlot(qt.QMainWindow):
self._data = { self._data = {
MEAN: None, MEAN: None,
VARIANCE: None, VARIANCE: None,
SKENESS: None, SKEWNESS: None,
KURTOSIS: None, KURTOSIS: None,
INTENSITY: None INTENSITY: None
} }
...@@ -65,6 +66,12 @@ class MappingPlot(qt.QMainWindow): ...@@ -65,6 +66,12 @@ class MappingPlot(qt.QMainWindow):
# deal with toolbar # deal with toolbar
self._modesToolbar = qt.QToolBar(parent=self) self._modesToolbar = qt.QToolBar(parent=self)
_axis_widget = qt.QWidget(parent=self)
_axis_widget.setLayout(qt.QHBoxLayout())
self._dimCB = qt.QComboBox(parent=_axis_widget)
_axis_widget.layout().addWidget(qt.QLabel('dims:', self))
_axis_widget.layout().addWidget(self._dimCB)
# self._modesToolbar.addWidget(_axis_widget)
self._intensityAction = None self._intensityAction = None
self._meanAction = None self._meanAction = None
self._varianceAction = None self._varianceAction = None
...@@ -74,7 +81,7 @@ class MappingPlot(qt.QMainWindow): ...@@ -74,7 +81,7 @@ class MappingPlot(qt.QMainWindow):
(INTENSITY, self._intensityAction), (INTENSITY, self._intensityAction),
(MEAN, self._meanAction), (MEAN, self._meanAction),
(VARIANCE, self._varianceAction), (VARIANCE, self._varianceAction),
(SKENESS, self._skenessAction), (SKEWNESS, self._skenessAction),
(KURTOSIS, self._kurtosisAction), (KURTOSIS, self._kurtosisAction),
]) ])
...@@ -86,10 +93,16 @@ class MappingPlot(qt.QMainWindow): ...@@ -86,10 +93,16 @@ class MappingPlot(qt.QMainWindow):
action.changed.connect(callback) action.changed.connect(callback)
self._modesToolbar.addAction(action) self._modesToolbar.addAction(action)
self._actions[mode] = action self._actions[mode] = action
self._modesToolbar.setIconSize(qt.QSize(50, 50)) if mode == INTENSITY:
self._modesToolbar.addWidget(_axis_widget)
self._modesToolbar.setIconSize(qt.QSize(70, 70))
self._plot.addToolBar(qt.Qt.LeftToolBarArea, self._modesToolbar) self._plot.addToolBar(qt.Qt.LeftToolBarArea, self._modesToolbar)
self.setMode(INTENSITY) self.setMode(INTENSITY)
# connect signal / SLOT
self._dimCB.currentIndexChanged.connect(self._updatePlot)
def setMode(self, mode): def setMode(self, mode):
""" """
Change the display mode Change the display mode
...@@ -107,11 +120,18 @@ class MappingPlot(qt.QMainWindow): ...@@ -107,11 +120,18 @@ class MappingPlot(qt.QMainWindow):
def _updatePlot(self): def _updatePlot(self):
self._plot.clear() self._plot.clear()
if self._operation is None:
return
mode = self._getMode() mode = self._getMode()
if self._data[mode] is not None: dim = self._getDim()
self._plot.addImage(self._data[mode])
assert dim in self._operation.dims
if mode == INTENSITY:
data = self._operation.intensity_map
else: else:
self._plot.clear() assert hasattr(self._operation.dims[dim], mode)
data = getattr(self._operation.dims[dim], mode)
self._plot.addImage(data)
def _getMode(self): def _getMode(self):
for action_mode, action in self._actions.items(): for action_mode, action in self._actions.items():
...@@ -119,18 +139,38 @@ class MappingPlot(qt.QMainWindow): ...@@ -119,18 +139,38 @@ class MappingPlot(qt.QMainWindow):
return action_mode return action_mode
return None return None
def _getDim(self):
"""
:return: dimension the user want to display
:rtype: int or None
"""
if self._dimCB.size() is 0 or self._dimCB.currentText() == "":
return None
else:
return int(self._dimCB.currentText())
def setOperation(self, operation): def setOperation(self, operation):
self._data[MEAN] = self._data[VARIANCE] = None self._operation = operation
self._data[SKENESS] = self._data[KURTOSIS] = None
if operation is not None: if operation is not None and operation.ndim > 0:
assert isinstance(operation, _MappingBase) assert isinstance(operation, _MappingBase)
self._data[INTENSITY] = operation.intensity_map self._updateAxis()
self._data[MEAN] = operation.dim[0][:, :, 0] self._updatePlot()
self._data[VARIANCE] = operation.dim[0][:, :, 1]
self._data[SKENESS] = operation.dim[0][:, :, 2] def _updateAxis(self):
self._data[KURTOSIS] = operation.dim[0][:, :, 3] self._dimCB.blockSignals(True)
self._dimCB.clear()
if self._operation is None:
return
current_dim = self._getDim()
for dim in self._operation.dims:
self._dimCB.addItem(str(dim))
idx = self._dimCB.findText(current_dim)
if idx >= 0:
self._dimCB.setCurrentIndex(idx)
self._dimCB.blockSignals(False)
self._updatePlot() self._updatePlot()
......
...@@ -133,7 +133,7 @@ class ShiftCorrectionWidget(qt.QWidget): ...@@ -133,7 +133,7 @@ class ShiftCorrectionWidget(qt.QWidget):
self._rawData = self._operation.data self._rawData = self._operation.data
self._lastShift = (0, 0, 0) self._lastShift = (0, 0, 0)
self._stack.setStack(self._operation.data_flatten) self._stack.setStack(self._operation.data_flatten)
self._sumPlot.addImage(self._operation.data.sum(axis=0)) self._sumPlot.addImage(self._operation.data_flatten.sum(axis=0))
self._updateShift() self._updateShift()