Richardson–Lucy deconvolution
The Richardson–Lucy algorithm, also known as Lucy–Richardson deconvolution, is an iterative procedure for recovering an underlying image that has been blurred by a known point spread function. It was named after William Richardson and Leon Lucy, who described it independently. Description. When an image is produced using an optical system and detected using photographic film or a charge-coupled device, for instance, it is inevitably blurred, with an ideal point source not appearing as a point but being spread out into what is known as the point spread function. Extended sources can be decomposed into  the sum of many individual point sources, thus the observed image can be represented in terms of a transition matrix p operating on an underlying image: where formula_2 is the intensity of the underlying image at pixel formula_3 and formula_4 is the detected intensity at pixel formula_5. In general, a matrix whose elements are formula_6 describes the portion of light from source pixel j that is detected in pixel i. In most good optical systems (or in general, linear systems that are described as shift invariant) the transfer function p can be expressed simply in terms of the spatial "offset  between the source pixel j and the observation pixel i: where formula_8 is called a point spread function. In that case the above equation becomes a convolution. This has been written for one spatial dimension, but of course most imaging systems are two dimensional, with the source, detected image, and point spread function all having two indices. So a two dimensional detected image is a convolution of the underlying image with a two dimensional point spread function formula_9 plus added detection noise. In order to estimate formula_10 given the observed formula_11 and a known formula_12 we employ the following iterative 	Richardson–Lucy deconvolution
9662620	procedure in which the estimate" of formula_10 which we call formula_14 for iteration number "t" is updated as follows: where It has been shown empirically that if this iteration converges, it converges to the maximum likelihood solution for formula_10. Writing this more generally for two (or more) dimensions in terms of convolution with a point spread function P: where the division and multiplication are element wise, formula_19 indicates a 2D convolution, and formula_20 is the flipped point spread function. In problems where the point spread function formula_21 is not known "a priori", a modification of the Richardson–Lucy algorithm has  been proposed, in order to accomplish blind deconvolution. Derivation. In the context of fluorescence microscopy, the probability of measuring a set of number of photons (or digitalization counts proportional to detected light) formula_22 for expected values formula_23 for a detector with formula_24 pixels is given by it is convenient to work with formula_26 since in the context of maximum likelihood estimation we want to find the position of the "maximum" of the likelihood function and we are not interested in its absolute value. Again since formula_28 is a constant, it will not give any additional information regarding the position of  the maximum, so let's consider where formula_30 is something that shares the same maximum position as formula_31. Now let's consider that formula_32 comes from a "ground truth" formula_33 and a measurement formula_34 which we assume to be linear. Then where a matrix multiplication is implied. We can also write this in the form where we can see how formula_37, mixes or blurs the ground truth. It can also be shown that the derivative of an element of formula_38, formula_39 with respect to some other element of formula_33 can be written as: Tip: it's easy to see this  by writing a matrix formula_34 of say (5 x 5) and two arrays formula_38 and formula_43 of 5 elements and check it. This last equation can interpreted as how much "one" element of formula_43, say element formula_5 influences the "other" elements formula_46 (and of course the case formula_47 is also taken into account). For example in a typical case an element of the ground truth formula_43 will influence nearby elements in formula_38 but not the very distant ones (a value of formula_50 is expected on those matrix elements). Now, the key and arbitrary step: we don't know formula_33 but  we want to estimate it with formula_52, let's call formula_53 and formula_54 the estimated ground truths while we are using the RL algorithm, where the "hat" symbol is used to distinguish ground truth from estimator of the ground truth Where formula_55 stands for a formula_56-dimensional gradient. If we work on the derivative of formula_57 we get and if we now use () we get But we can also note that formula_60 by definition of transpose matrix. And hence Then if we consider formula_3 spanning all the elements from formula_62 to formula_56 this equation can be rewritten in its vectorial  form where formula_65 is a matrix and formula_66, formula_32 and formula_68 are vectors. Let's now propose the following arbitrary and key step where formula_68 is a vector of ones of size formula_56 (same as formula_66, formula_32 and formula_33) and the division is element-wise. Using () and () we can rewrite () as which yields Where division refers to element-wise matrix division and formula_65 operates as a matrix but the division and the product (implicit after formula_76) are element-wise. Also, formula_77 can be calculated because we assume - The initial guess formula_78 is known (and is typically set to be the  experimental data) - The "measurement" function formula_34 is known On the other hand formula_80 is the experimental data. Therefore, equation () applied successively, provides an algorithm to estimate our ground truth formula_81 by ascending (since it moves in the direction of the gradient of the likelihood) in the likelihood "landscape". It has not been demonstrated in this derivation that it converges and no dependence on the initial choice is shown. Note that equation () provides a way of following the direction that increases the likelihood but the choice of the log-derivative is arbitrary. On the other hand equation () introduces a way of  weighting the movement from the previous step in the iteration. Note that if this term was not present in () then the algorithm would output a movement in the estimation even if formula_82. It's worth noting that the only strategy used here is to maximize the likelihood at all cost, so artifacts on the image can be introduced. It is worth noting that no prior knowledge on the shape of the ground truth formula_43 is used in this derivation.