Skip to content
GitLab
Commit bcf1f0e9 authored by Carsten Richter's avatar Carsten Richter
Browse files

Merge branch 'spherical-coords' into 'master' 

Add support of spherical coordinates for QSpace

Closes #59

See merge request !89
parents 1609865d e037f994
Pipeline #6533 passed with stages
in 5 minutes and 8 seconds
Expand all Hide whitespace changes
Inline Side-by-side
xsocs/examples/id01_peak.py
View file @ bcf1f0e9
import numpy as np
import os
#!/usr/bin/python
# coding: utf-8
"""
This scripts illustrates the API of the gaussian/center-of-mass processing
"""
from xsocs.process.fit import peak_fit
import sys
import numpy
from xsocs.process.fit import PeakFitter, FitTypes, BackgroundTypes
# output directory (some temporary files will also be written there)
workdir = '/path/to/workdir/'
# path to the hdf5 file written by the img_to_qspace function
# path to the hdf5 file written by the kmap_2_qspace function
qspace_f = '/path/to/qspace.h5'
# result file
result_file = os.path.join(workdir, 'results.txt')
# Name of the text file where to store the resuls
result_file = 'results.txt'
# positions (on the sample) to convert to qspace
# indices = array with indices (of the sample positions array)
# List of QSpace indices to process
# This selects sample positions for which QSpace are processed.
# If None, all QSpace are processed
indices = None
# number of processes to use
# If None, will use the number of availble core (see multiprocessing.cpu_count)
# Number of processors to use
# If None, all available cores are used.
n_proc = None
results, success = peak_fit.peak_fit(qspace_f,
indices=indice,
fit_type=peak_fit.FitTypes.GAUSSIAN,
n_proc=n_proc)
with open(result_file, 'w+') as res_f:
res_f.write('# X Y qx qy qz q I valid\n')
for i, s in enumerate(success):
x_pos = results[i, 0]
y_pos = results[i, 1]
xpeak = results[i, 3]
ypeak = results[i, 6]
zpeak = results[i, 9]
xpeak_max = results[i, 2]
q = np.sqrt(xpeak**2 + ypeak**2 + zpeak**2)
r = (x_pos, y_pos, xpeak, ypeak, zpeak, q, xpeak_max, s)
res_str = '{0} {1} {2} {3} {4} {5} {6} {7} ({8})\n'.format(i, *r)
res_f.write(res_str)
# Set-up the processing for a gaussian fit without background subtraction
fitter = PeakFitter(qspace_f,
fit_type=FitTypes.GAUSSIAN,
indices=indices,
n_proc=n_proc,
roi_indices=None, # QSpace ROI
background=BackgroundTypes.NONE)
# Run the processing and get the results
results = fitter.peak_fit()
# Check for errors
if fitter.status != fitter.DONE:
print("Fit process failed")
sys.exit()
# Prepare columns header and data
# Start with position (x, y) on the sample
headers = ['sample_x', 'sample_y']
values = [results.sample_x, results.sample_y]
# Add one column for each parameter of the fit/COM result
# For each axis of the QSpace
for dimension, axis_name in enumerate(results.qspace_dimension_names):
# For each fitted/computed parameter of the result
for parameter_name in results.available_result_names:
# Add its name to the header and its data to the values
headers.append(axis_name + '_' + parameter_name)
values.append(results.get_results(dimension, parameter_name))
# Transpose values from (parameters x points) to (points x parameters)
values = numpy.transpose(numpy.array(values))
# Write results to text file
with open(result_file, 'w') as res_f:
res_f.write('\t'.join(headers) + '\n')
for row_values in values:
res_f.write('\t'.join(str(v) for v in row_values) + '\n')
xsocs/fit/Centroid.py deleted 100644 → 0
View file @ 1609865d
#!/usr/bin/python
# coding: utf8
# /*##########################################################################
#
# Copyright (c) 2015-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.
#
# ###########################################################################*/
from __future__ import absolute_import
__authors__ = ["D. Naudet"]
__date__ = "01/01/2017"
__license__ = "MIT"
import numpy as np
from .Plotter import Plotter
from ..process.fit.Fitter import Fitter
from ..process.fit.fitresults import FitStatus
from ..process.fit.sharedresults import FitSharedResults
from ..process.fit.fitresults import FitResult
class CentroidFitter(Fitter):
def fit(self, i_fit, i_cube, qx_profile, qy_profile, qz_profile):
profiles = qx_profile, qy_profile, qz_profile
for iax, axis in enumerate(self._AXIS_NAMES):
y = profiles[iax]
Sum = y.sum()
if Sum != 0:
x = getattr(self, "_%s"%axis)
com = x.dot(y) / Sum
#idx = np.abs(x - com).argmin()
max_idx = y.argmax()
I_max = y[max_idx]
x_max = x[max_idx]
self._shared_results.set_results(axis, i_fit,
[com, Sum, I_max, x_max],
FitStatus.OK)
else:
self._shared_results.set_results(axis, i_fit,
[np.nan, np.nan, np.nan, np.nan],
FitStatus.FAILED)
class CentroidResults(FitSharedResults):
def __init__(self,
n_points=None,
shared_results=None,
shared_status=None,
**kwargs):
super(CentroidResults, self).__init__(n_points=n_points,
n_params=4,
n_peaks=1,
shared_results=shared_results,
shared_status=shared_status)
def fit_results(self, *args, **kwargs):
fit_name = 'Centroid'
fitresults = FitResult(fit_name, *args, **kwargs)
for axis in self._AXIS_NAMES:
results = getattr(self, "_npy_%s_results"%axis)
status = getattr(self, "_npy_%s_status"%axis)
for i_p, param in enumerate(("COM", "I_sum", "I_max", "Pos_max")):
fitresults.add_result(axis, 'centroid', param, results[:, 0+i_p].ravel())
fitresults.set_status(axis, status)
return fitresults
class CentroidPlotter(Plotter):
def plotFit(self, plot, x, peakParams):
plot.setGraphTitle('center of mass')
for peakName, peak in peakParams.items():
center = peak.get('COM')
xmax = peak.get('Pos_max')
if np.isfinite(center):
plot.addXMarker(center, legend='center of mass', text="com")
plot.addXMarker(xmax, legend='maximum position',
text="max",
color="gray")
def getPlotTitle(self):
return 'Center Of Mass'
# process = fitH5.processes(entry)[0]
#
# positions = fitH5.get_result(entry, process, 'COM')
#
# plots[0].addCurve(xAcqQX, yAcqQX, legend='measured')
# plots[0].addXMarker(positions.qx[index], legend='center of mass')
# plots[0].setGraphTitle('QX center of mass')
#
# plots[1].addCurve(xAcqQY, yAcqQY, legend='measured')
# plots[1].addXMarker(positions.qy[index], legend='center of mass')
# plots[1].setGraphTitle('QY center of mass')
#
# plots[2].addCurve(xAcqQZ, yAcqQZ, legend='measured')
# plots[2].addXMarker(positions.qz[index], legend='center of mass')
# plots[2].setGraphTitle('QZ center of mass')
#
xsocs/fit/Gaussian.py deleted 100644 → 0
View file @ 1609865d
#!/usr/bin/python
# coding: utf8
# /*##########################################################################
#
# Copyright (c) 2015-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.
#
# ###########################################################################*/
from __future__ import absolute_import
__authors__ = ["D. Naudet"]
__date__ = "01/01/2017"
__license__ = "MIT"
import numpy as np
from scipy.optimize import leastsq
from .Plotter import Plotter
from ..process.fit.Fitter import Fitter
from ..process.fit.fitresults import FitStatus
from ..process.fit.sharedresults import FitSharedResults
from ..process.fit.fitresults import FitResult
# Some constants
_const_inv_2_pi_ = np.sqrt(2 * np.pi)
class GaussianFitter(Fitter):
def __init__(self, *args, **kwargs):
super(GaussianFitter, self).__init__(*args, **kwargs)
self._n_peaks = self._shared_results._n_peaks
def fit(self, i_fit, i_cube, qx_profile, qy_profile, qz_profile):
profiles = dict(zip(self._AXIS_NAMES,
[qx_profile, qy_profile, qz_profile]))
for axis in self._AXIS_NAMES: # qx, qy, qz
x = getattr(self, "_%s"%axis)
y = profiles[axis]
# compute guess:
area = y.sum()/self._n_peaks * (x[-1]-x[0])/len(x)
cen = x[y.argmax()]
sigma = area / (y.max() * _const_inv_2_pi_)
guess = np.tile([area, cen, sigma],self._n_peaks)
if self._n_peaks>1: # we don't actually use many peaks yet
# a stupid way to distribute the starting values...
idx = np.where(y>(y.max()/20.))[0]
i_cen = idx[np.arange(0, len(idx), len(idx)//self._n_peaks)]
cens = x[i_cen]
guess[1::3] = cens
fit, success = gaussian_fit(x, y, guess)
self._shared_results.set_results(axis, i_fit, fit, success)
class GaussianResults(FitSharedResults):
def __init__(self,
n_points=None,
n_peaks=1,
shared_results=None,
shared_status=None):
super(GaussianResults, self).__init__(n_points=n_points,
n_params=3,
n_peaks=n_peaks,
shared_results=shared_results,
shared_status=shared_status)
def fit_results(self, *args, **kwargs):
fit_name = 'Gaussian'
fitresults = FitResult(fit_name, *args, **kwargs)
for axis in self._AXIS_NAMES:
results = getattr(self, "_npy_%s_results"%axis)
status = getattr(self, "_npy_%s_status"%axis)
for i_peak in range(self._n_peaks):
peak_name = 'gauss_{0}'.format(i_peak)
i_start = i_peak * 3 # 3 parameters
for i_p, param in enumerate(("Area", "Center", "Sigma")):
fitresults.add_result(axis, peak_name, param,
results[:, i_start+i_p].ravel())
fitresults.set_status(axis, status)
return fitresults
# 1d Gaussian func
# TODO : optimize
def gaussian(x, a, c, s):
"""
Returns (a / (sqrt(2 * pi) * s)) * exp(- 0.5 * ((x - c) / s)^2)
:param x: values for which the gaussian must be computed
:param a: area under curve ( amplitude * s * sqrt(2 * pi) )
:param c: center
:param s: sigma
:return: (a / (sqrt(2 * pi) * s)) * exp(- 0.5 * ((x - c) / s)^2)
"""
# s /= _two_sqrt_2_log_2
return (a * (1. / (_const_inv_2_pi_ * s)) *
np.exp(-0.5 * ((x - c) / s) ** 2))
# 1d Gaussian fit
# TODO : optimize
def gaussian_fit_err(p, x, y):
"""
:param p:
:param x:
:param y:
:return:
"""
n_p = len(p) // 3
gSum = 0.
for i_p in range(n_p):
gSum += gaussian(x, *p[i_p*3:i_p*3 + 3])
return gSum - y
# return gaussian(x, *p) - y
_two_sqrt_2_log_2 = 2 * np.sqrt(2 * np.log(2))
def gaussian_fit(x, y, p):
"""
Fits (leastsq) a gaussian on the provided data f(x) = y.
p = (a, c, s)
and f(x) = (a / (sqrt(2 * pi) * s)) * exp(- 0.5 * ((x - c) / s)^2)
:param x:
:param y:
:param p:
:return: area, position, fwhm
"""
result = leastsq(gaussian_fit_err,
p,
args=(x, y,),
maxfev=100000,
full_output=True)
if result[4] not in [1, 2, 3, 4]:
return [np.nan] * len(p), FitStatus.FAILED
return result[0], FitStatus.OK
def _gauss_first_guess(x, y):
i_max = y.argmax()
y_max = y[i_max]
p1 = x[i_max]
i_fwhm = np.where(y >= y_max / 2.)[0]
fwhm = (x[1] - x[0]) * len(i_fwhm)
p2 = fwhm / np.sqrt(2 * np.log(2)) # 2.35482
p0 = y_max * np.sqrt(2 * np.pi) * p2
return [p0, p1, p2]
class GaussianPlotter(Plotter):
def plotFit(self, plot, x, peakParams):
for peakName, peak in peakParams.items():
height = peak.get('Area')
position = peak.get('Center')
width = peak.get('Sigma')
params = [height, position, width]
if np.all(np.isfinite(params)):
fitted = gaussian(x, *params)
plot.addCurve(x,
fitted,
legend='{0}'.format(peakName),
color='red')
def getPlotTitle(self):
return 'Gaussian'
xsocs/fit/SilxFitter.py deleted 100644 → 0
View file @ 1609865d
#!/usr/bin/python
# coding: utf8
# /*##########################################################################
#
# Copyright (c) 2015-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.
#
# ###########################################################################*/
from __future__ import absolute_import
__authors__ = ["D. Naudet"]
__date__ = "01/01/2017"
__license__ = "MIT"
import numpy as np
from silx.math.fit import fittheories, sum_agauss
from silx.math.fit.fitmanager import FitManager
from .Plotter import Plotter
from ..process.fit.Fitter import Fitter
from ..process.fit.fitresults import FitStatus
from ..process.fit.sharedresults import FitSharedResults
from ..process.fit.fitresults import FitResult
class SilxFitter(Fitter):
p_types = ['A', 'P', 'F']
def __init__(self, *args, **kwargs):
super(SilxFitter, self).__init__(*args, **kwargs)
self._n_peaks = self._shared_results._n_peaks
self._fit = FitManager()
self._fit.loadtheories(fittheories)
self._fit.settheory('Area Gaussians')
self._results = np.zeros(3 * self._n_peaks)
def fit(self, i_fit, i_cube, qx_profile, qy_profile, qz_profile):
profiles = qx_profile, qy_profile, qz_profile
fit = self._fit
results = self._results
for iax, axis in enumerate(self._AXIS_NAMES):
failed = False
x = getattr(self, "_%s"%axis)
y = profiles[iax]
fit.setdata(x=x, y=y)
try:
fit.estimate()
fit.runfit()
except Exception as ex:
failed = True
results[:] = np.nan
if not failed:
for param in fit.fit_results:
p_name = param['name']
p_type = p_name[0]
peak_idx = int(p_name[-1]) - 1
if peak_idx >= self._n_peaks:
continue
# TODO : error management
param_idx = self.p_types.index(p_type)
results[peak_idx * 3 + param_idx] = param['fitresult']
self._shared_results.set_results(axis, i_fit, results,
FitStatus.OK)
else:
self._shared_results.set_results(axis, i_fit, results,
FitStatus.FAILED)
class SilxResults(FitSharedResults):
def __init__(self,
n_points=None,
n_peaks=1,
shared_results=None,
shared_status=None):
super(SilxResults, self).__init__(n_points=n_points,
n_params=3,
n_peaks=n_peaks,
shared_results=shared_results,
shared_status=shared_status)
def fit_results(self, *args, **kwargs):
fit_name = 'SilxFit'
fitresults = FitResult(fit_name, *args, **kwargs)
for axis in self._AXIS_NAMES:
results = getattr(self, "_npy_%s_results"%axis)
status = getattr(self, "_npy_%s_status"%axis)
for i_peak in range(self._n_peaks):
peak_name = 'gauss_{0}'.format(i_peak)
i_start = i_peak * 3
for i_p, param in enumerate(("Area", "Center", "FWHM")):
fitresults.add_result(axis, peak_name, param,
results[:, i_start+i_p].ravel())
fitresults.set_status(axis, status)
return fitresults
class SilxPlotter(Plotter):
def plotFit(self, plot, x, peakParams):
for peakName, peak in peakParams.items():
area = peak.get('Area')
position = peak.get('Center')
width = peak.get('FWHM')
params = [area, position, width]
fitSum = None
if np.all(np.isfinite(params)):
fitted = sum_agauss(x, *params)
plot.addCurve(x,
fitted,
legend='{0}'.format(peakName),
color='red')
if fitSum is None:
fitSum = fitted
else:
fitSum += fitted
if fitSum is not None:
plot.addCurve(x, fitSum, legend='Sum')
def getPlotTitle(self):
return 'Silx Gaussian Fit'
xsocs/gui/XsocsGui.py
View file @ bcf1f0e9
......@@ -214,7 +214,7 @@ class XsocsGui(Qt.QMainWindow):
view = node.getView(self)
if view not in self.__openedViews:
self.__openedViews.append(view)
view.sigFitDone.connect(self.__slotFitDone)
view.sigFitDone.connect(self.__slotFitDone)
view.show()
view.setAttribute(Qt.Qt.WA_DeleteOnClose, True)
if bringToFront:
......
xsocs/gui/model/Model.py
View file @ bcf1f0e9
......@@ -41,10 +41,8 @@ class RootNode(Node):
def __init__(self, *args, **kwargs):
super(RootNode, self).__init__(*args, **kwargs)
for colIdx in range(len(self.ColumnNames)):
self._setDataInternal(colIdx,
self.ColumnNames[colIdx],
Qt.Qt.DisplayRole)
for index, name in enumerate(self.ColumnNames):
self._setDataInternal(index, name, Qt.Qt.DisplayRole)
class Model(Qt.QAbstractItemModel):
......@@ -56,9 +54,9 @@ class Model(Qt.QAbstractItemModel):
sigRowsRemoved = Qt.Signal(object, int, int)