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

  1. Using the GUI by invoking RichLucyGUI
  2. Calling the module RichLucy from a system command line
  3. 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:

  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
  2. 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

  1. Realizations of the bandpass and the measured spectrum are drawn from normal distributions corresponding to the provided mean and standard uncertainties
  2. If necessary, the obtained bandpass function is interpolated to equidistant wavelength axis
  3. If necessary, the obtained spectrum values are interpolated to the same step size as the bandpass
  4. If necessary, the bandpass is normalized to integrate to 1.0
  5. 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

Toolbox

The following functions are also part of the RichLucy module:

Reading and processing files

RichLucy.LoadFromFile(filename)[source]

Load wavelength scale axis, measured values and uncertainties from file. Output of this method is the matrix “data”.

Parameters:filename – string, name of the file to load (including full path)
Returns:data matrix
Raises:ValueError if data file could not be loaded.
RichLucy.Interpolate2Equidist(lamb_axis, values, display=True, interp_kind='cubic', tol_equi=1e-05)[source]

Interpolation of measured values if wavelength scale axis is not equidistant. Input of this function are the original measured values and the corresponding wavelength scale. The output is the interpolated version with the new scale, a flag whether interpolated or not and a tuple (numNaN,numNeg) with the number of NaN and negative values after interpolation.

Parameters:
  • lamb_axis – ndarray of shape (N,) with wavelength values
  • values – ndarray of shape (N,) with measured values

optional parameters

Parameters:
  • display – boolean whether to print information to the command line (default:True)
  • interp_kind – string which is handed over to the scipy interp1d function (default:’cubic’)
  • tol_equi – float, tolerance value to judge whether wavelength axis is equidistant (default:1e-5)

Outputs

Returns:newAx, ndarray of shape (L,) with equidistant wavelength scale
Returns:newVals, ndarray of shape (L,) with corresponding interpolated spectrum values
Returns:did_interp, boolean whether interpolation was actually necessary
Returns:numNaN, numNeg, ints - the number of NaN and negative values after interpolation
RichLucy.Interpolate2SameStep(M, lambda_M, delta_b, display=True, interp_kind='cubic', tol_equi=1e-05)[source]

Interpolation of measured spectrum if its wavelength step size is not equal to that of the bandpass function or if it is not equidistant. The output is the interpolated version with the new scale, a flag whether interpolated or not and a tuple (numNaN,numNeg) with the number of NaN and negative values after interpolation.

Parameters:
  • M – ndarray of shape (N,) with the measured spectrum values
  • lambda_M – ndarray of shape (N,) with the corresponding wavelengths
  • delta_b – float, the wavelength step size for the bandpass function

optional parameters

Parameters:
  • display – boolean, whether to print information to the command line
  • interp_kind – string, which is passed to the scipy interp1d function
  • tol_equi – float, tolerance for judgin whether the wavelength step sizes are equal

Outputs

Returns:lm_interp - ndarray of shape (K,) of new wavelength values for measured spectrum
Returns:M_interp - ndarray of shape (K,) of corresponding interpolated spectrum values
Returns:did_interp - boolean whether interpolation was actually necessary
Returns:numNaN,numNeg - ints, the number of NaN and negative values after interpolation
RichLucy.checkBandpass(lambda_b, b, print_output=False)[source]

For given bandpass values b and corresponding wavelengths lambda_b, this function determines whether:: * the lambda axis is strictly monotone increasing * the lambda axis is equidistant * the bandpass is normalized to integrate to 1

Parameters:
  • lambda_b – ndarray of shape (N,) with wavelengths for the bandpass function
  • b – ndarray of shape (N,) with corresponding bandpass values
Returns:

result - a list: [lambda monotonically increasing?, is equidistant?, is normalized?, delta_b]

RichLucy.checkSpectrum(lambda_M, M, print_output=False)[source]

For given vector of measured spectrum values and corresponding wavelengths this function determines whether:: * the lambda axis is strictly monotonically increasing * the lambda axis is equidistant

Parameters:
  • lambda_M – ndarray of shape (N,) with wavelengths for the measured spectrum
  • M – ndarray of shape (N,) with corresponding spectrum values
Returns:

result - the list: [lambda_M monotonically increasing?, is equidistant?, delta_M]

RichLucy.InterpolateBandpass(lamb_axis, values)[source]

Interpolation of bandpass values (line spread function). This is just a wrapper to allow for a convenient function call from EXCEL.

RichLucy.fix_inputs(lambda_bandpass, bandpass, lambda_spectrum, spectrum, display=True)

For given wavelength scales (bandpass and measurement each), this functions verifies whether interpolation of the data is necessary in order to carry out the Richardson-Lucy method and returns the interpolated results.

If no interpolation is necessary, then the original data is returned.

Parameters:
  • lambda_bandpass – ndarray of wavelength scale values for bandpass
  • bandpass – ndarray of bandpass values
  • lambda_spectrum – ndarray of wavelength scale values for measured spectrum
  • spectrum – ndarray of measured spectrum amplitudes
Returns lb,b,lM,M:
 

each corrected (interpolated) if necessary

Returns changes:
 

an array of 2 boolean (bandpass interpolated?, measured spectrum interpolated?)

RichLucy.fix_raw_inputs(bdata, Mdata, display=True)

Same functionality as fix_inputs, but taking as inputs the matrices which result from reading the data text files

Parameters:
  • bdata – ndarray of shape (N,2) or (N,3) with N the number of bandpass values
  • Mdata – ndarray of shape (K,2) or (K,3) with K the number of spectrum values
Returns lb,b,lM,M:
 

wavelength scales and values interpolated if necessary

Return changes:

an array of 2 boolean (bandpass interpolated?, measured spectrum interpolated?)

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

  1. Download pyxll from www.pyxll.com

  2. Unzip the downloaded file and copy it somewhere

  3. 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)
  4. 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

  • MaxCurv
  • EME
  • RMSthreshold

@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: 1+lambda*Laplace(SH(k-1,:))

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: 1-lambda*divergence(gradient(SH(k-1,:)/abs(gradient(SH(k-1,:))) 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

Contents

Indices and tables