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 B. Lucy, who described it independently. == Description == When an image is produced using an optical system and detected using photographic film, a charge-coupled device or a CMOS sensor, for example, 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: d i = ∑ j p i , j u j , {\displaystyle d_{i}=\sum _{j}p_{i,j}u_{j},} where u j {\displaystyle u_{j}} is the intensity of the underlying image at pixel j {\displaystyle j} , and d i {\displaystyle d_{i}} is the detected intensity at pixel i {\displaystyle i} . In general, a matrix whose elements are p i , j {\displaystyle p_{i,j}} 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: p i , j = P ( i − j ) , {\displaystyle p_{i,j}=P(i-j),} where P ( Δ i ) {\displaystyle P(\Delta i)} is called a point spread function. In that case the above equation becomes a convolution. This has been written for one spatial dimension, but 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 P ( Δ x , Δ y ) {\displaystyle P(\Delta x,\Delta y)} plus added detection noise. In order to estimate u j {\displaystyle u_{j}} given the observed d i {\displaystyle d_{i}} and a known P ( Δ i x , Δ j y ) {\displaystyle P(\Delta i_{x},\Delta j_{y})} , the following iterative procedure is employed in which the estimate of u j {\displaystyle u_{j}} (called u ^ j ( t ) {\displaystyle {\hat {u}}_{j}^{(t)}} ) for iteration number t is updated as follows: u ^ j ( t + 1 ) = u ^ j ( t ) ∑ i d i c i p i j , {\displaystyle {\hat {u}}_{j}^{(t+1)}={\hat {u}}_{j}^{(t)}\sum _{i}{\frac {d_{i}}{c_{i}}}p_{ij},} where c i = ∑ j p i j u ^ j ( t ) , {\displaystyle c_{i}=\sum _{j}p_{ij}{\hat {u}}_{j}^{(t)},} and ∑ j p i j = 1 {\displaystyle \sum _{j}p_{ij}=1} is assumed. It has been shown empirically that if this iteration converges, it converges to the maximum likelihood solution for u j {\displaystyle u_{j}} . Writing this more generally for two (or more) dimensions in terms of convolution with a point spread function P: u ^ ( t + 1 ) = u ^ ( t ) ⋅ ( d u ^ ( t ) ⊗ P ⊗ P ∗ ) , {\displaystyle {\hat {u}}^{(t+1)}={\hat {u}}^{(t)}\cdot \left({\frac {d}{{\hat {u}}^{(t)}\otimes P}}\otimes P^{}\right),} where the division and multiplication are element-wise, ⊗ {\displaystyle \otimes } indicates a 2D convolution, and P ∗ {\displaystyle P^{}} is the mirrored point spread function, or the inverse Fourier transform of the Hermitian transpose of the optical transfer function. In problems where the point spread function p i j {\displaystyle p_{ij}} 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) m = [ m 0 , … , m K ] {\displaystyle \mathbf {m} =[m_{0},\dots ,m_{K}]} for expected values E = [ E 0 , … , E K ] {\displaystyle \mathbf {E} =[E_{0},\dots ,E_{K}]} for a detector with K + 1 {\displaystyle K+1} pixels is given by P ( m ∣ E ) = ∏ i K Poisson ( E i ) = ∏ i K E i m i e − E i m i ! . {\displaystyle P(\mathbf {m} \mid \mathbf {E} )=\prod _{i}^{K}\operatorname {Poisson} (E_{i})=\prod _{i}^{K}{\frac {E_{i}^{m_{i}}e^{-E_{i}}}{m_{i}!}}.} Since in the context of maximum-likelihood estimation the aim is to locate the maximum of the likelihood function without concern for its absolute value, it is convenient to work with ln ( P ) {\displaystyle \ln(P)} : ln P ( m ∣ E ) = ∑ i K [ ( m i ln E i − E i ) − ln ( m i ! ) ] . {\displaystyle \ln P(\mathbf {m} \mid \mathbf {E} )=\sum _{i}^{K}[(m_{i}\ln E_{i}-E_{i})-\ln(m_{i}!)].} Moreover, since ln ( m i ! ) {\displaystyle \ln(m_{i}!)} is a constant, it does not give any additional information regarding the position of the maximum, so consider α ( m ∣ E ) = ∑ i K [ m i ln E i − E i ] , {\displaystyle \alpha (\mathbf {m} \mid \mathbf {E} )=\sum _{i}^{K}[m_{i}\ln E_{i}-E_{i}],} where α {\displaystyle \alpha } is something that shares the same maximum position as P ( m ∣ E ) {\displaystyle P(\mathbf {m} \mid \mathbf {E} )} . Now consider that E {\displaystyle \mathbf {E} } comes from a ground truth x {\displaystyle \mathbf {x} } and a measurement H {\displaystyle \mathbf {H} } which is assumed to be linear. Then E = H x , {\displaystyle \mathbf {E} =\mathbf {H} \mathbf {x} ,} where a matrix multiplication is implied. This can also be written in the form E m = ∑ n K H m n x n , {\displaystyle E_{m}=\sum _{n}^{K}H_{mn}x_{n},} where it can be seen how H {\displaystyle H} mixes or blurs the ground truth. It can also be shown that the derivative of an element of E {\displaystyle \mathbf {E} } , ( E i ) {\displaystyle (E_{i})} with respect to some other element of x j {\displaystyle x_{j}} can be written as It is easy to see this by writing a matrix H {\displaystyle \mathbf {H} } of, say, 5 × 5 and two arrays E {\displaystyle \mathbf {E} } and x {\displaystyle \mathbf {x} } of 5 elements and check it. This last equation can be interpreted as how much one element of x {\displaystyle \mathbf {x} } , say element i {\displaystyle i} , influences the other elements j ≠ i {\displaystyle j\neq i} (and of course the case i = j {\displaystyle i=j} is also taken into account). For example, in a typical case an element of the ground truth x {\displaystyle \mathbf {x} } will influence nearby elements in E {\displaystyle \mathbf {E} } but not the very distant ones (a value of 0 {\displaystyle 0} is expected on those matrix elements). Now, the key and arbitrary step: x {\displaystyle \mathbf {x} } is not known but may be estimated by x ^ {\displaystyle {\hat {\mathbf {x} }}} . Let's call x ^ old {\displaystyle {\hat {\mathbf {x} }}_{\text{old}}} and x ^ new {\displaystyle {\hat {\mathbf {x} }}_{\text{new}}} the estimated ground truths while using the RL algorithm, where the hat symbol is used to distinguish ground truth from estimator of the ground truth where ∂ ∂ x {\displaystyle {\frac {\partial }{\partial \mathbf {x} }}} stands for a K {\displaystyle K} -dimensional gradient. Performing the partial derivative of α ( m ∣ E ( x ) ) {\displaystyle \alpha (\mathbf {m} \mid \mathbf {E} (\mathbf {x} ))} yields the following expression: ∂ α ( m ∣ E ( x ) ) ∂ x j = ∂ ∂ x j ∑ i K [ m i ln E i − E i ] = ∑ i K [ m i E i ∂ ∂ x j E i − ∂ ∂ x j E i ] = ∑ i K ∂ E i ∂ x j [ m i E i − 1 ] . {\displaystyle {\frac {\partial \alpha (\mathbf {m} \mid \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\sum _{i}^{K}[m_{i}\ln E_{i}-E_{i}]=\sum _{i}^{K}\left[{\frac {m_{i}}{E_{i}}}{\frac {\partial }{\partial x_{j}}}E_{i}-{\frac {\partial }{\partial x_{j}}}E_{i}\right]=\sum _{i}^{K}{\frac {\partial E_{i}}{\partial x_{j}}}\left[{\frac {m_{i}}{E_{i}}}-1\right].} By substituting (1), it follows that ∂ α ( m ∣ E ( x ) ) ∂ x j = ∑ i K H i j [ m i E i − 1 ] . {\displaystyle {\frac {\partial \alpha (\mathbf {m} \mid \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}=\sum _{i}^{K}H_{ij}\left[{\frac {m_{i}}{E_{i}}}-1\right].} Note that H j i T = H i j {\displaystyle H_{ji}^{T}=H_{ij}} by the definition of a matrix transpose. And hence Since this equation is true for all j {\displaystyle j} spanning all the elements from 1 {\displaystyle 1} to K {\displaystyle K} , these K {\displaystyle K} equations may be compactly rewritten as a single vectorial equation ∂ α ( m ∣ E ( x ) ) ∂ x = H T [ m E − 1 ] , {\displaystyle {\frac {\partial \alpha (\mathbf {m} \mid \mathbf {E} (\mathbf {x} ))}{\partial \mathbf {x} }}=\mathbf {H} ^{T}\left[{\frac {\mathbf {m} }{\mathbf {E} }}-\mathbf {1} \right],} where H T {\displaystyle \mathbf {H} ^{T}} is a matrix, and m {\displaystyle \mathbf {m} } , E {\displaystyle \mathbf {E} } and 1 {\displaystyle \mathbf {1} } are vectors. Now, as a seemingly arbitrary but key step, let where 1 {\displaystyle \mathbf {1} } is a vector of ones of size K {\displaystyle K} (same as m {\displaystyle \mathbf {m} } , E {\displaystyle \mathbf {E} } and x {\displaystyle \mathbf {x} } ), and the d
Read more →
