Richardson-Lucy Python documentation
Introduction
The Richardson-Lucy method is an iterative deconvolution method for spectral measurements.
Given knowledge about the line-spread function (the reversed bandpass function) and the
measured spectrum, this method estimates the actual input spectrum.
The Python Code can be used in several ways
- Using the GUI by invoking RichLucyGUI
- Calling the module RichLucy from a system command line
- Calling the methods of RichLucy from a Python shell or a Python function
Disclaimer
Copyright S. Eichstaedt, F. Schmaehling (PTB) 2013
sascha.eichstaedt@ptb.de
This software is licensed under the BSD-like license:
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution.
This software was developed at Physikalisch-Technische Bundesanstalt.
The software is made available ‘as is’ free of cost. PTB assumes
no responsibility whatsoever for its use by other parties, and makes no
guarantees, expressed or implied, about its quality, reliability, safety,
suitability or any other characteristic. In no event will PTB be liable
for any direct, indirect or consequential damage arising in connection
with the use of this software.
Citation when using this software is:
Eichstaedt, Schmaehling, Wuebbeler, Anhalt, Buenger, Krueger & Elster
‘Comparison of the Richardson-Lucy method and a classical approach for
spectrometer bandwidth correction’, Metrologia vol. 50, 2013
DOI: 10.1088/0026-1394/50/2/107
In case of questions or suggestions contact Sascha.Eichstaedt@ptb.de
This module contains the following functions:
- calcOptCurv
- caclOptDiff
- RichLucy
- EstMeas
- Interpolate2Equidist
- Interpolate2SameStep
- InterpolateBandpass
- LoadFromFile
- RichLucyUncertainty
- RichLucyDeconv
The different use cases
Graphical User Interface (GUI)
The GUI is called simply by running RichLucyGUI. It allows to load data files
containing the bandpass and the measured spectrum values, respectively. The Richardson-Lucy
deconvolution method is then called simply by clicking on “Start”. Several method
parameters can be adjusted and the results be saved to txt-files. The plotted figures
can also be saved to disk.
Running from system command
When the module RichLucy is called from a system command line:
python RichLucy.py [bandpass file] [measured spectrum file] and [output file]
it runs the full Richardson-Lucy deconvolution method. The [xxx file] entries must be the
full path to the corresponding file.
In addition the maximum number of iterations (maxiter) can be specified - the default value
is 500. If in addition to that a fifth argument is specified it is interpreted as the number
of Monte Carlo trials and uncertainty evaluation is carried out.
For example:
python RichLucy.py mybandpass.txt myspectrum.txt result.txt 800 1000
will run Richardson-Lucy deconvolution for the spectrum loaded from myspectrum.txt using the
bandpass from mybandpass.txt with 800 iterations maximum. Uncertainties will be calculated using
1000 Monte Carlo trials. The results will be stored in result.txt.
In case that interpolation of the bandpass and/or the measured spectrum is necessary, this
information will be printed to the command line and the intermediate results stored to txt-files
for inspection by the user.
Calling particular methods
When the module RichLucy is imported from a Python command or another Python module the whole
set of methods are available to carry out the Richardson-Lucy deconvolution.
This includes in particular the below described functions
In addition the following functions may be helpful to simulate measurements using a triangular
bandpass of prescribed width and skewness.
-
RichLucy.sim_bandpass(bandwidth, sampledist, skew=1, noise=0, noise_low=0)[source]
Generate a triangular bandpass function with given bandwidth
and skewnes (skew) on the wavelength axis [-bandwidth,bandwidth]
Relative measurement noise can be simulated using noise and noise_low
as maximum and minimum relative noise levels.
Parameters: |
- bandwidth – float - the bandwidth of the bandpass function
- sampledist – float - the wavelength step size
- skew – float - skewness of the triangle (left-skewed for <1.0, right-skewed for >1.0)
- noise – float - upper limit for the relative bandpass uncertainty
- noise_low – float - lower limit for the relative bandpass uncertainty
|
Returns: | hat, lamb - the bandpass function and the corresponding wavelengths
|
-
RichLucy.sim_measurement(Sfun, lamb, b, bandwidth, delta_b, noise_low=0, noise_up=0)[source]
Generate a Gaussian-shaped spectrum by simulating a measurement of a
clean Gaussian-shaped function Sfun defined on the wavelength axis lamb
with a triangular bandpass function b of given bandwidth and wavelength
step size delta_b.
Note that Sfun needs to be a function object.
Measurement noise can be simulated using the relative noise levels
noise_low and noise_up
Parameters: |
- Sfun – a function S(lambda) returning the spectrum at wavelengths lambda
- lamb – ndarray of shape (N,) with wavelengths for measurement
- b – ndarray of shape (L,) - the bandpass function values
- bandwidth – float - the bandwidth of the bandpass function
- delta_b – float - the wavelength step size of the bandpass function
- noise_low – float - lower limit for the relative measurement uncertainty
- noise_up – float - upper limit for the relative measurement uncertainty
|
Returns: | M - ndarray of shape (N,) with simulated measured spectrum values
|
The usage of these functions is outlined in the function demo()
-
RichLucy.demo(bandwidth=40, skew=1, FWHM=10, maxiter=2000)[source]
Demonstrate usage of the Richardson-Lucy method for bandpass correction.
A triangular bandpass of given bandwidth and skewness is used to simulate
a measurement of a Gaussian-shaped spectrum of given FWHM value.
The maximum number of Richardson-Lucy iterations is given by maxiter.
Two different automatic stopping criteria are employed and the results
are compared.
Main functions
-
RichLucy.RichLucy(b, delta, M, maxiter=500, autostop=True, display=True, optFun=<function calcOptCurv at 0x03638570>, regFun=None, calcCrit=True, returnAll=False, initialShat=None)[source]
Richardson Lucy iteration under the assumption that the sampling of b and M is equidistant and
that the wavelength step size of b and M is equal.
Parameters: |
- b – ndarray of shape (N,) with the bandpass function values
- delta – float, wavelength step size for b and M
- M – ndarray of shape (L,) with the measured spectrum values
|
optional parameters
Parameters: |
- maxiter – (optional) maximum number of iterations, default is 500
- autostop – (optional) boolean whether to automatically find optimal iteration number (default: True)
- display – (optional) boolean whether to print information to the command line (default: True)
- optFun – (optional) function handle to a function that calculates the optimal iteration number (default: calcOptCurv)
:param regFun: (optional) function handle to a function that regularizes each RichLucy iteration (e.g., Tikhonov)
- calcCrit – (optional) boolean whether to calculate the stopping criterion irrespective of autostop (default:True)
- returnAll – (optional) boolean whether to return all intermediate estimates (default: False)
- initialShat – (optional) ndarray of shape (L,) with an initial estimate of the input spectrum (default: None)
|
Output:
Returns: | Shat, ndarray of shape (L,) - the estimated spectrum |
Returns: | critVals1, ndarray of shape (maxiter,) containing values calculated from optFun (rms of progress for calcOptCurv) |
Returns: | critVals2, ndarray of shape (maxiter,) containing other values calculated from optFun (curvature of rms of progress for calcOptCurv) |
-
RichLucy.RichLucyUncertainty(bhat, b_unc, lambda_b, Mhat, M_unc, lambda_M, runs=100, maxiter=500, autostop=1, display=1, maxncpu=10, optFun=<function calcOptCurv at 0x03638570>, calcCrit=True, returnAll=False, plotRes=False)[source]
This function carries out uncertainty evaluation for the Richardson-Lucy
deconvolution by using the GUM Supplement 2 Monte Carlo method. It assumes
that knowledge about the input quantities (i.e., bandpass and measured spectrum)
is available in terms of an estimate and a vector of associated standard uncertainties.
TODO: Extend the uncertainty evaluation to arbitrary type of input knowledge.
Parameters: |
- bhat – ndarray of shape (N,) with estimate of bandpass values
- b_unc – ndarray of shape (N,) with associated standard uncertainties
- lambda_b – ndarray of shape (N,) with corresponding wavelengths
- Mhat – ndarray of shape (L,) with measured spectrum values
- M_unc – ndarray of shape (L,) with associated standard uncertainties
- lambda_M – ndarray of shape (L,) with corresponding wavelengths
|
optional arguments
Parameters: |
- runs – int - the number of Monte Carlo runs for uncertainty propagation
- maxiter – int - the maximum number of RichLucy iterations per run
- autostop – boolean - whether to determine optimal stopping point
- display – boolean - whether to print information to the command line
- maxncpu – int - the maximum number of parallel workers for the Monte Carlo
- optFun – function for the determination of the optimal stopping point
- calcCrit – boolean - whether to calculate the stopping criterion using optFun irrespective of autostop
- returnAll – boolean - whether to return all intermediate results
- plotRes – boolean - whether to plot result
|
Returns: | Shat, UShat - ndarrays of shape (L,) with the mean estimated spectrum and its associated uncertainty, respectively
|
The function RichLucyUncertainty() for uncertainty evaluation uses the
Python multiprocessing module to evaluate the Monte Carlo trials in parallel
if the system has more than 2 threads available.
Each Monte Carlo trial consists of the following steps
- Realizations of the bandpass and the measured spectrum are drawn from normal distributions corresponding to the provided mean and standard uncertainties
- If necessary, the obtained bandpass function is interpolated to equidistant wavelength axis
- If necessary, the obtained spectrum values are interpolated to the same step size as the bandpass
- If necessary, the bandpass is normalized to integrate to 1.0
- The Richardson-Lucy method RichLucy() is called.
-
RichLucy.RichLucyDeconv(lambda_bandpass, bandpass, lambda_spectrum, spectrum, display=True, returnAll=False, **kwargs)[source]
Application of Richardson-Lucy deconvolution for provided bandpass (line spread function)
and measured spectrum. This method carries out interpolation if the sampling of bandpass
and spectrum does not meet the RichLucy assumptions. If necessary, the estimation result
is returned interpolated to the original wavelength values.
Parameters: |
- lambda_bandpass – ndarray of shape (N,) with wavelengths for the bandpass function
- bandpass – ndarray of shape (N,) with corresponding bandpass values
- lambda_spectrum – ndarray of shape (L,) with wavelengths for the measured spectrum
- spectrum – ndarray of shape (L,) with corresponding measured spectrum values
|
optional arguments
Parameters: | maxiter – int - maximum number of Richardson-Lucy iterations |
Returns: | Shat - ndarray of shape (L,) with the estimated spectrum values |
If returnAll = True:
Returns: | RMS - ndarray of shape (maxiter,) with root-mean-squared values of the estimation progress |
Returns: | curv - ndarray of shape (maxiter,) - the curvature of RMS |
EXCEL
Some of the code can also be run from Microsoft EXCEL.
Note
Using RichLucy code in EXCEL requires to install the pyXLL EXCEL Add-In
Installing pyXLL
Installation and configuration of the pyXLL Add-In is carried out as follows
Download pyxll from www.pyxll.com
Unzip the downloaded file and copy it somewhere
- Open the pyxll.cfg and adjust the following entries:
- pythonpath: replace the default entry with the your Python code is located (incl RichLucy)
- modules: replace the default entries with the modules you would like to use from within EXCEL (in particular: RichLucy)
Open EXCEL and install the pyxll Add-In (see: EXCEL documentation)
Note
If the RichLucy.py is in a subdirectory of ‘pythonpath’ then write in modules folder.RichLucy, where ‘folder’ is the name of the folder where you installed RichLucy
The following functions can be called from Microsoft EXCEL:
Stopping rules
Collection of criteria for automatic stopping of RichLucy iterations.
The RichLucy.RichLucy function requires to call such a function as
Crit1,Crit2,indopt = calcCrit(SH,M,bdtilde,display,autostop)
with inputs
- SH ....... a (maxiter,len(M)) dimensional matrix of RichLucy intermediate results
- M ....... the measured spectrum
- bdtilde .. the (discrete) line-spread function
- display .. boolean, whether to print information to screen
- autostop.. boolean, whether to actually calculate the automatic stopping
and outputs
Crit1 .... a len(maxiter) vector of criterion values (e.g. rms of progress)
Crit2 .... a len(maxiter) vector of other criterion values (e.g. curvature)
indopt ... the determined optimal stopping
If anything went wrong, indopt should be equal to -1.
The stopping rules can be tested using the function Test_CalcOptFun in the
module test_modules.py.
implemented stopping rules
@author: Sascha Eichstaedt
-
stopping_rules.MaxCurv(SH, M, btilde, display=True, autostop=True)[source]
Calculate the root-mean-squared progress in the RichLucy iterations.
Determine the optimal number of iterations from the curvature of the rms values.
Parameters: |
- SH – ndarray of shape (N,K) with N the number of iterations and K the number of wavelengths
- M – ndarray of shape (K,) containing the measured spectrum values (possibly after interpolation)
- btilde – ndarray of shape (L,) containg the values of the reversed bandpass function (i.e., the line spread function)
- autostop – (optional) int, if equal to 1 then the point of maximum curvature is calculated
- display – (optional) boolean, if True then comments and intermediate results are printed to the command line
- orig_mode – (optional) boolean, if True then the optimal point is the global maximum curvature, if False then the optimal
point is taken to be that local maximum with smallest error in the estimated measurement
|
Returns: | RMS, ndarray of shape (N,) containing the root-mean-squared progress in the RL estimate per iteration
|
Returns: | curv ndarray of shape (N,) containg the curvature values of the RMS as a function of iteration
|
Returns: | indopt the optimal stopping point as calculated by this criterion
|
This criterion is published in:
Eichstaedt, Schmaehling, Wuebbeler, Anhalt, Buenger, Krueger & Elster
Comparison of the Richardson-Lucy method and a classical approach for
spectrometer bandwidth correction, Metrologia vol. 50, 2013
-
stopping_rules.EEM(SH, M, bdtilde, display=True, autostop=True)[source]
Determine the optimal stopping by calculating the estimated measurement
for each intermediate estimate and the weighted rms difference from the actual
measurement. The optimal stopping point is that with the first local minimal
weighted rms difference.
Parameters: |
- SH – ndarray of shape (N,K) with N the number of iterations and K the number of wavelengths
- M – ndarray of shape (K,) containing the measured spectrum values (possibly after interpolation)
- btilde – ndarray of shape (L,) containg the values of the reversed bandpass function (i.e., the line spread function)
- autostop – (optional) int, if equal to 1 then the point of maximum curvature is calculated
- display – (optional) boolean, if True then comments and intermediate results are printed to the command line
|
Returns: | wDiff, ndarray of shape (N,) containing the weighted rms errors of the estimated measurement for each iteration
|
Returns: | wDiff (twice, to have the same input-output structure as calcOptCurv)
|
Returns: | indopt, int – the index of the optimal stopping iteration
|
This criterion is published in:
Jan Audenaert et al. Bayesian deconvolution method applied to experimental bidirectional
transmittance distribution functions, Meas. Sci. Techol. 24, 2013
Regularization methods
Regularization of RichLucy is carried out by damping the update in each iteration.
The function RichLucy.RichLucy therefore uses
Shat_new = regFun(Shat,delta)*(RL_update)
where RL_update is conv(M/(conv(Shat,b)),b(-lambda)).
Example: RichLucy with Tikhonov regularization
Shat = RichLucy.RichLuy(b,delta,M,regFun=lambda Sh,delta: Tikhonov(Sh,delta,reglambda))
with reglambda the chosen regularization parameter (e.g. reglambda=1e-5)
Implemented regularization methods
- maxEntropy
- Tikhonov
- TotalVar
@author: Sascha Eichstaedt
-
regularization_methods.maxEntropy(Shat, deltaS, reglambda=1e-05)[source]
maximum entropy regularization
- Criterion introduced for RichLucy in
- J Boutet de Monvel, S Le Calvez, and M Ulfendahl
Image Restoration for Confocal Microscopy: Improving the Limits of
Deconvolution, with Application to the Visualization of the Mammalian Hearing Organ
Biophysical Journal Volume 80 May 2001 2455–2470
-
regularization_methods.Tikhonov(Shat, deltaS, reglambda=1e-05)[source]
Regularization using Tikhonov regularization
Therefore, each estimated spectrum in SH is weighted with the regularizer:

i.e., the smoothness of the previous iteration.
Criterion introduced for RichLucy in
G van Kempen, Image Restoration in Flourescence Microscopy, PhD Thesis,
Technische Universiteit Delft, Holland, 1999
-
regularization_methods.TotalVar(Shat, deltaS, reglambda=1e-05)[source]
Determination of optimal number of iterations by using total variation regularization
Therefore, each estimated spectrum in SH is weighted with the regularizer:
Criterion introduced for RichLucy in
N Dey, L Blanc-Feraud, C Zimmer, P Roux, Z Kam, J C Olivo-Marin and J Zerubia
3D Microscopy Deconvolution using Richardson-Lucy Algorithm with Total Variation Regularization
Rapport de Recherce, No 8272, July 2004, INRIA