Readme

deconvolution

A Python module providing Deconvolution class that implements and generalises Ruifrok-Johnston color deconvolution algorithm [RJ], [IJ]. It allows one to split an image into distinct color layers in just a few lines of code:

from deconvolution import Deconvolution
from PIL import Image

img = Image.open("image.jpg")

# Declare an instance of Deconvolution, with image loaded and with color basis defining what layers are interesting
decimg = Deconvolution(image=img, basis=[[1, 0.1, 0.2], [0, 0.1, 0.8]])

# Constructs new PIL Images, with different color layers
layer1, layer2 = decimg.out_images(mode=[1, 2])

Installation

You can install the package using pip:

pip install deconvolution

Alternatively, you can clone the repository and run:

make install

Since then you can import use the module in your scripts:

from deconvolution import Deconvolution
d = Deconvolution()

Testing

# For Python 3 users
make test

# For Python 2 users
make comp

# Check the code coverage
make coverage

# Check the coverage interactively, using a web browser
make html

Deconvolve

For better usage experience we created a script allowing one to deconvolve images from the shell. Copy deconvolve.py file into /usr/local/bin or, if you want to use it locally:

mkdir ~/bin
cp deconvolve.py ~/bin
export PATH=~/bin:$PATH

Since then you can deconvolve images using:

deconvolve.py image1.png image2.png ...
# For help
deconvolve.py -h

Documentation

Check out our documentation at Read The Docs.

Contributors

Method developed by Frederic Grabowski generalising Ruifrok-Johnston algorithm [RJ]. and implemented by Frederic Grabowski [FG] and Paweł Czyż [PC]. We are very grateful to prof. Daniel Wójcik and dr Piotr Majka [N1], [N2] who supervised the project. We also would like to thank prof. Gabriel Landini [GL], who implemented the colour deconvolution in ImageJ [IJ] and allowed us to test the algorithm on his data.

Ruifrok-Johnston deconvolution

Assume that we have an image and we know three stains. We can deconvolve image and get density layers:

"""
 In this example, we have an image and a known basis with three vectors. We want to get the density layers.

 In this example we use a cropped image [1], released under CC licence.

 [1] https://en.wikipedia.org/wiki/Chronic_lymphocytic_leukemia#/media/File:Chronic_lymphocytic_leukemia_-_high_mag.jpg
"""
from deconvolution import Deconvolution
from PIL import Image


def join_horizontally(*args):
    """Joins many PIL images of the same dimensions horizontally"""
    w, h = args[0].size
    n = len(args)
    joined = Image.new("RGB", (n*w, h))
    for x_off, img in zip(range(0, n*w, w), args):
        joined.paste(img, (x_off, 0))
    return joined


if __name__ == "__main__":
    # Load an image
    original = Image.open("cropped.jpg")

    # Create a deconvolution object with the image
    dec = Deconvolution(image=original, basis=[[0.91, 0.38, 0.71], [0.39, 0.47, 0.85], [1, 0, 0]])

    # Produce density matrices for both colors. Be aware, as Beer's law do not always hold.
    first_density, second_density, third_density = dec.out_scalars()
    print(first_density.shape, second_density.shape, third_density.shape)

    # Produce reconstructed image, first layer, second layer, third layer and rest
    out_images = dec.out_images(mode=[0, 1, 2, 3, -1])
    # Original image, reconstructed, layers and rest
    join_horizontally(original, *out_images).show()

Two stain deconvolution

Alternatively, we can use two stains. If we know them, we can just mimic the procedure for two stains:

"""
 In this example, we have an image and a known basis with two vectors. We want to get the density layers.

 In this example we use a cropped image [1], released under CC licence.

 [1] https://en.wikipedia.org/wiki/Chronic_lymphocytic_leukemia#/media/File:Chronic_lymphocytic_leukemia_-_high_mag.jpg
"""
from deconvolution import Deconvolution
from PIL import Image


def join_horizontally(*args):
    """Joins many PIL images of the same dimensions horizontally"""
    w, h = args[0].size
    n = len(args)
    joined = Image.new("RGB", (n*w, h))
    for x_off, img in zip(range(0, n*w, w), args):
        joined.paste(img, (x_off, 0))
    return joined


if __name__ == "__main__":
    # Load an image
    original = Image.open("cropped.jpg")

    # Create a deconvolution object with the image
    dec = Deconvolution(image=original, basis=[[0.91, 0.38, 0.71], [0.39, 0.47, 0.85]])

    # Produce density matrices for both colors. Be aware, as Beer's law do not always hold.
    first_density, second_density = dec.out_scalars()
    print(first_density.shape, second_density.shape)

    # Produce reconstructed image, first layer, second layer and rest
    out_images = dec.out_images(mode=[0, 1, 2, -1])
    # Original image, reconstructed, layers and rest
    join_horizontally(original, *out_images).show()

But if we do not know at least one of them, we should first find appropriate vectors:

"""
 In this example, we have an image and an insufficient number of color vectors - we will use Deconvolution
 complete_basis method to find other vectors.
 Then we try to make them more independent applying resolve_dependencies with different belligerency parameter.
 We show all the bases found.

 In this example we use a cropped image [1], released under CC licence.

 [1] https://en.wikipedia.org/wiki/Chronic_lymphocytic_leukemia#/media/File:Chronic_lymphocytic_leukemia_-_high_mag.jpg
"""
from deconvolution import Deconvolution
from PIL import Image


def join_horizontally(*args):
    """Joins many PIL images of the same dimensions horizontally"""
    w, h = args[0].size
    n = len(args)
    joined = Image.new("RGB", (n*w, h))
    for x_off, img in zip(range(0, n*w, w), args):
        joined.paste(img, (x_off, 0))
    return joined


def join_vertically(*args):
    """Joins many PIL images of the same dimensions vertically"""
    w, h = args[0].size
    n = len(args)
    joined = Image.new("RGB", (w, n*h))
    for y_off, img in zip(range(0, n*h, h), args):
        joined.paste(img, (0, y_off))
    return joined

if __name__ == "__main__":
    # Load an image
    original = Image.open("cropped.jpg")

    # Create a deconvolution object with the image
    dec = Deconvolution(image=original, sample_density=6)
    # Complete basis - as we did not provide 2 or 3 vectors, it needs to be found
    dec.complete_basis()
    # We can get the basis found
    print("Basis before resolve:\n{}\n".format(dec.pixel_operations.get_basis()))

    # Produce reconstructed image, first layer, second layer and rest
    out_images1 = dec.out_images(mode=[0, 1, 2, -1])
    # Original image, reconstructed, layers and rest
    before_resolve = join_horizontally(original, *out_images1)

    # Resolve dependencies - make the vectors more independent
    dec.resolve_dependencies(belligerency=0.1)
    # We can get the basis found
    print("Basis after resolve:\n{}\n".format(dec.pixel_operations.get_basis()))
    # And produce reconstructed image, first layer, second layer and rest
    out_images2 = dec.out_images(mode=[0, 1, 2, -1])
    # Construct a joined image from original image, reconstructed, layers and rest
    after_resolve = join_horizontally(original, *out_images2)

    # Resolve dependencies with huge belligerency - make the vectors very independent
    dec.resolve_dependencies(belligerency=1)
    print("Basis after aggressive resolve:\n{}\n".format(dec.pixel_operations.get_basis()))
    # Produce reconstructed image, first layer, second layer and rest
    out_images3 = dec.out_images(mode=[0, 1, 2, -1])
    # Show original image, reconstructed, layers and rest
    after_huge_resolve = join_horizontally(original, *out_images3)

    # Show all images for visual comparision
    join_vertically(before_resolve, after_resolve, after_huge_resolve).show()

Mathematical description of color deconvolution

Introduction

Let us imagine a light source that sends light through layers of different substances - each of them absorbing some of the light passing through. Given a digital recording of the light that passed through this setup (and information of the input light), we would like to reconstruct the layout of the substances.

Each pixel is represented by a three-component vector (it’s RGB components). I will make frequent use of the Hadamard (or entry-wise) product:

\[\begin{split}\begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} \bullet \begin{bmatrix} y_1\\ y_2\\ y_3 \end{bmatrix}= \begin{bmatrix} x_1y_1\\ x_2y_2\\ x_3y_3 \end{bmatrix}\end{split}\]

The Hadamard product is similar to the standard multiplication of real numbers and I will use the same notation for inverses and powers.

Consider light represented by a RGB-vector \(\vec i\) passing through a unit-wide layer of some substance. I will assume that Beer’s law [LB] hold’s, that is the RGB-vector of the outgoing light will be:

\[\begin{split}\vec i \bullet \vec s= \begin{bmatrix} i_1s_1\\ i_2s_2\\ i_3s_3 \end{bmatrix},\end{split}\]

where \(\vec s\) is specific for this substance (the choice of “unitary width” does change this value; this issue will be addressed later on). If, for example, \(\vec s=(1,1,1)\) the layer does not absorb any light, for \(\vec s=(0,1,1)\) the red channel is completely absorbed, but the remaining channels are untouched. Should the layer be \(a\) times wider, the RGB components of the outgoing light would be:

\[\begin{split}\vec i \bullet \vec s^{\ a}=\begin{bmatrix} i_1\cdot s_1^a\\ i_2\cdot s_2^a\\ i_3\cdot s_3^a \end{bmatrix},\end{split}\]

for multiple substances,

\[\begin{split}\vec i \bullet \vec p^{\ a} \bullet \vec q^{\ b} \bullet ...=\begin{bmatrix} i_1\cdot p_1^a \cdot q_1^b \cdot ...\\ i_2\cdot p_2^a \cdot q_2^b \cdot ...\\ i_3\cdot p_3^a \cdot q_3^b \cdot ... \end{bmatrix}.\end{split}\]

This equation implies that changing the order of layers or splitting some of them does not change the outgoing light. I will work under the assumption that this also holds for mixed substances.

Ruifrok-Johnston deconvolution

In the case studied by Ruifrok and Johnston [RJ], the light \(\vec i\) passes through three substances. Vectors \(($\vec v$, $ \vec u$, $\vec w$)\) describing absorption rates for all substances (that is, absorption coefficients for unit-wide substance layers) are assumed to be known. The width of each substance layer (which may change from point to point) has to be calculated given the output light. Suppose the camera registers a single vector \(\vec r\) at some arbitrary pixel. We wish to express this vector according to the equation:

\[\label{eq:decon3} \vec i \bullet \vec v^{\ a} \bullet \vec u^{\ b} \bullet \vec w^{\ c} = \vec r.\]

Solving this equation for \(a, b\) and \(c\) gives the layers widths in unit lengths. We may then compute how much light would have passed through each layer separately: the first deconvolution is just \(\vec i \bullet \vec v^{\ a}\), the second \(\vec i \bullet \vec u^{\ b}\) and the third \(\vec i \bullet \vec w^{\ c}\). It is possible that no real non-negative \(a, b, c\) solve this equation (due to data noise, imperfect digitization, traces of other substances, etc.). In that case, the reconstructed image will differ from the original: this difference can be visualized by considering “rest picture”:

\[\vec r \bullet \vec v^{\ -a} \bullet \vec u^{\ -b} \bullet \vec w^{\ -c}.\]

Let us briefly return to the issue of picking a unit width. Notice that changing the unit width of a substance by a factor of \(\lambda\), changes the constant $a$ at each pixel to \(a/ \lambda\), and thus the density distribution \(a\) times the unit width does not change. Similarly, the deconvolved images don’t change:

\[\vec i \bullet \vec v^{\ a} = \vec i \bullet (\vec v^{\ \lambda})^{a/ \lambda}.\]

Equation (ref{eq:decon3}) is in fact a set of three equations - one for each component. Lets rewrite it for the \(k\)-th component and transform equivalently:

\[\begin{split}v_k^a\cdot u_k^b \cdot w_k^c &= \frac{r_k}{i_k}\\ a \log v_k + b \log u_k + c \log w_k &= \log(r_k/i_k),\end{split}\]

going back to vector representation:

\[\begin{split}\label{eq:3vec} \begin{bmatrix} \log v_1 & \log u_1 & \log w_1\\ \log v_2 & \log u_2 & \log w_2\\ \log v_3 & \log u_3 & \log w_3 \end{bmatrix} \begin{bmatrix} a\\ b\\ c \end{bmatrix} - \begin{bmatrix} \log(r_1/i_1)\\\log(r_2/i_2)\\\log(r_3/i_3) \end{bmatrix}=0.\end{split}\]

This equation is solvable if and only if the left matrix is invertible, and \(i_k \neq 0\). In physical terms, the first assumption states that no substance can be “faked” by mixing the remaining two, and the second that the input light is nonzero for all channels.

Given an input image, at each pixel \(\vec r\) we may solve for \(a,b,c\). Should we get any negative results, it means that this particular pixels color can not be obtained by mixing the given substances. In this case the assumption that the image was obtained by mixing the given substances is violated - hence it’s reasonable to disregard such results data noise and drop all negative parts. Having that done, we are able to construct five pixels:

  • Reconstructed pixel: \(\vec i \bullet \vec v^{\ a} \bullet \vec u^{\ b} \bullet \vec w^{\ c}\)
  • Difference from the original image, due to negative cut-off: \(\vec r \bullet \vec v^{\ -a} \bullet \vec u^{\ -b} \bullet \vec w^{\ -c}\)
  • Three single-substance pixels: + \(\vec i \bullet \vec v^{\ a}\) + \(\vec i \bullet \vec u^{\ b}\) + \(\vec i \bullet \vec w^{\ c}\)

After processing every pixel in this manner, the reconstructed image, three single substance images, and one remainder image (showing the error) are obtained. An example is shown in Figure ref{fig:ruifrok}.

Until now the approach was based on Ruifrok and Johnston - however, the choice of formalism makes it easier to look for further development. Handle any number of channels straightforward: in fact, for the general case, the notation does not change. Secondly, if only two substances are of interest, Ruifrok and Johnston suggest measuring the absorbances of those two substances, and then choosing the third so that it minimizes the negative cut-off. The third single-substance image is then used as a measure of error. This arbitrarity seems a bit artificial - now I introduce a method developed by Frederic Grabowski.

Two-substance deconvolution

Choosing the third substance by hand so that it minimizes data loss due to negative cut-offs introduces ambiguity into the measurement, and seems artificial. In order to fix this problem, drop the third substance entirely, and look for \(a, b\) that minimize the squared error of the approximation:

\[\begin{split}\begin{bmatrix} \log v_1 & \log u_1 \\ \log v_2 & \log u_2 \\ \log v_3 & \log u_3 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} - \begin{bmatrix} \log(r_1/i_1)\\\log(r_2/i_2)\\\log(r_3/i_3) \end{bmatrix} \approx 0.\end{split}\]

Where both the matrix and the 3-vector are given. Clean up the notation:

\[\inf_{\vec x\ \in\ \mathbb{R}^2} f(x) = \inf_{\vec x\ \in\ \mathbb{R}^2}||A\vec x-\vec y||^2\]

We want to know the \(\vec x\) for which this infimum is obtained. This \(\vec x\) always exists, because it is the orthogonal projection of \(\vec y\) onto \(A(\mathbb{R}^2)\).

For each minimum \(\text{grad}\, f=0\). That is, for the first component:

\[\begin{split}0=\frac12 \frac{\partial f}{\partial x_1} (x_1,x_2)=\sum_{k=1}^3 A_{k1}\left(A_{k1}x_1+A_{k2}x_2-r_k\right) \\ x_1\sum_{k=1}^3 A_{k1}^2 + x_2\sum_{k=1}^3 A_{k1}A_{k2} = \sum_{k=1}^3 A_{k1}r_k\end{split}\]

Combining (ref{eq:solAB}) with the equation for the second component we get a set of two linear equations:

\[\begin{split}\sum_k \begin{bmatrix} A_{k1}^2 & A_{k1}A_{k2} \\ A_{k1}A_{k2} & A_{k2}^2 \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix} = \begin{bmatrix} \sum_k A_{k1}r_k\\ \sum_k A_{k2}r_k \end{bmatrix}.\end{split}\]

The equation above can be easily solved if and only if it’s determinant is not 0:

\[\begin{split}\det \begin{bmatrix} \sum_k A_{k1}^2 & \sum_k A_{k1}A_{k2} \\ \sum_k A_{k1}A_{k2} & \sum_k A_{k2}^2 \end{bmatrix} = \\ \sum_k A_{k1}^2 \cdot \sum_k A_{k2}^2 - \left(\sum_k A_{k1}A_{k2} \right)^2 \neq 0\end{split}\]

The Cauchy-Schwarz inequality states that the considered determinant is 0 if and only if there is a number \(t\) for which \(A_{k1}=t\cdot A_{k2}\). This is again the mixing independence of the basis. If \(\dim A(\mathbb{R}^2) <2\), then there does not exist a unique \(\vec x\) for which the projected vector is obtained.

We now have a method for finding “the best” \(a,\ b\) solving equation (ref{eq:2vec}). This means that for each pixel and a basis of two given substances, we are able to calculate four pixels: - Best reconstructed pixel: \(\vec i \bullet \vec v^{\ a} \bullet \vec u^{\ b}\) - Difference from the original image: \(\vec r \bullet \vec v^{\ -a} \bullet \vec u^{\ -b}\) - Two single-substance pixels: \(\vec i \bullet \vec v^{\ a}\) and \(\vec i \bullet \vec u^{\ b}\)

After processing every pixel in this manner the reconstructed image, two single substance images, and one remainder image (showing the error) are obtained. An example is shown in Figure ref{fig:auto}.

Formulation of the optimization problem

Considering deconvolutions with two substances has another advantage - it gives a criterium for comparing bases. Taking

\[\sum_{p\ \in\ pixels}\ \inf_{a,b,c\ \in\ \mathbb{R}} ||\text{LHS\ of eq.}(\ref{eq:3vec})||^2\]

does not work, because the equation is always soluble and the expression is identically zero (at least for all independent bases). Decreasing the number of degrees of freedom (that is, the number of substances to match) solves this difficulty:

\[\begin{split}\begin{bmatrix} \log v_1 & \log u_1 \\ \log v_2 & \log u_2 \\ \log v_3 & \log u_3 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} - \begin{bmatrix} \log(r_1/i_1)\\\log(r_2/i_2)\\\log(r_3/i_3) \end{bmatrix} \approx 0.\end{split}\]
\[\begin{split}\sum_{p\ \in\ pixels}\ \inf_{a,b\ \in\ \mathbb{R}} \lvert\lvert \begin{bmatrix} \log v_1 & \log u_1 \\ \log v_2 & \log u_2 \\ \log v_3 & \log u_3 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} - \begin{bmatrix} \log(r_1/i_1)\\\log(r_2/i_2)\\\log(r_3/i_3) \end{bmatrix} \rvert\rvert^2\end{split}\]

To solve this optimization problem, first clean up the notation. Let:

\[\begin{split}A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32} \end{bmatrix},\end{split}\]
\[\begin{split}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}= \vec x = \begin{bmatrix} a\\ b \end{bmatrix}\end{split}\]
\[\begin{split}\begin{bmatrix} y_1(p)\\ y_2(p)\\ y_3(p) \end{bmatrix}= \vec y(p) = \begin{bmatrix} \log(r_1/i_1)\\\log(r_2/i_2)\\\log(r_3/i_3) \end{bmatrix}\end{split}\]

Given some \(\vec y(p)\), the problem is to find a \(2\times3\) matrix \(A\), that minimizes the expression:

\[\label{eq:optOrig} f(A) = \sum_{p\ \in\ pixels}\ \inf_{\vec x\ \in\ \mathbb{R}^2}||A\vec x-\vec y(p)||^2\]

Solving the optimization problem

For any \(A\) in equation ref{eq:optOrig}:

\[\inf_{\vec x\ \in\ \mathbb{R}^2}||A\vec x-\vec y(p)||^2 = d(\vec y(p),\ A(\mathbb{R}^2))^2,\]

hence we want to minimize the mean squared distance of the points $y$ from the image space of \(A\). There are two cases: either \(\dim A(\mathbb{R}^2) = 2\) or is strictly less than 2. In the second case, we are always able to choose such a matrix \(A\), that the previous image is a subspace of the new image, but then the distances can only be smaller. It thus suffices to find the best two dimensional space. Every such space has a normal vector, which we choose so that the third component in non-negative (this convention is arbitrary, but does not matter). We can now rewrite (ref{eq:optOrig}) as:

\[\label{eq:optG} f(A) = g(\vec n) = \sum_{p\ \in\ pixels}\ \lvert\vec n \cdot \vec y(p)\rvert^2\]

Any \(\vec n\) minimizing (ref{eq:optG}), determines a class of matrices minimizing (ref{eq:optOrig}) - precisely those, whose image is perpendicular to \(\vec n\). We only need to consider \(\vec n\) such that \(||\vec n||=1\). Because the set of such \(\vec n\) is compact in \(\mathbb{R}^3\), we can apply the method of Lagrange multipliers:

\[\nabla ( g(\vec n) -\lambda ||\vec n||^2 ) = 0,\]

after expanding and rearranging

\[\label{eq:eigenOrig} \sum_{p\ \in\ pixels}\ (\ \vec y(p) \cdot \vec n\ )\ \vec y(p) = \lambda \vec n.\]

The left hand side is a linear operator from \(\mathbb{R}^3\) to \(\mathbb{R}^3\) applied to \(\vec n\). Equation (ref{eq:eigenOrig}) is just the eigenvalue equation for this operator. Moreover, multiplying both sides by \(\vec n\) we notice \(\lambda = g(\vec n)\). The smallest-eigenvalue eigenvector is the \(\vec n\) minimizing \(g(\vec n)\), and the corresponding eigenvalue the value \(g(\vec n)\).

Computing the deconvolution

Rewrite the left hand side of equation (ref{eq:eigenOrig}) using index notation:

\[\sum_{p\ \in\ pixels}\ \sum_{i=1}^3 y_i(p) y_j(p)\cdot n_i,\]

is the \(j\)-th component of the resulting vector. Hence the matrix of the linear transformation is:

\[(Y)_{ij} = \sum_{p\ \in\ pixels}\ \sum_{i=1}^3 y_i(p) y_j(p)\]

Given an input image, first calculate \(Y\), and find it’s eigenvalue decomposition. Pick the eigenvector \(\vec n\) with smallest eigenvalue. Choosing an \(A\) such that \(\vec n\) is perpendicular to \(A\)’s image is equivalent to choosing a basis with both elements perpendicular to \(\vec n\). To have any preference, let’s return to the physical interpretation. Naively, all these bases allow us to mix the same set of colors: but not for all bases will this mixing will be physically meaningful. Consider the following example: the basis consisting of two substances, one absorbing only red the other only blue, will be equivalent to a basis of one substance absorbing both red and blue, and the other only blue. However, only in the first base is it possible to construct a color with only the red channel absorbed. This happens because we cannot physically have negative widths of substances. It seems advantageous to choose the basis that allows us to mix the widest range of colors physically. It turns out that this choice is not always optimal. For now, we stick with the maximal physical color range. The basis of our choosing is the one for which both vectors are non-negative (so that the resulting substances absorb light, and not amplify it), and have the biggest angle between them. This determines them uniquely up to a rearrangement.