Commit 7a5f82fd authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

3 colour flag analysis

parent 7fca8d34
...@@ -62,11 +62,11 @@ For a discrete number of $N$ wavelengths, we can now write for the absorbance of ...@@ -62,11 +62,11 @@ For a discrete number of $N$ wavelengths, we can now write for the absorbance of
\end{align} \end{align}
However, this column for $\col{k}$ has a (definite) length in this discrete formulation, and this length affects the (strength of) absorption of a certain `colour' as much as the factor $a$. For our purposes, we would like to separate `dye amount' and 'dye colour effect' completely, i.e.\ only $a$ should affect the strength of absorption. However, this column for $\col{k}$ has a (definite) length in this discrete formulation, and this length affects the (strength of) absorption of a certain `colour' as much as the factor $a$. For our purposes, we would like to separate `dye amount' and 'dye colour effect' completely, i.e.\ only $a$ should affect the strength of absorption.
So we divide $\col{k}$ by its length $l_k$ (`normalize') to get $\col{\hat{k}}$ and add the original effect of this length to $a$ to get $\hat{a}$: So we divide $\col{k}$ by its length $\abs{\col{k}}$ (`normalize') to get $\col{\hat{k}}$ and add the original effect of this length to $a$ to get $\hat{a}$:
\begin{align} \begin{align}
l_k & = \sqrt{k_1^2 + k_2^2 + \dots +k_N^3} = \sqrt{\sum_{n=1}^N k_n^2}\\ \abs{\col{k}} & = \sqrt{k_1^2 + k_2^2 + \dots +k_N^3} = \sqrt{\sum_{n=1}^N k_n^2} \label{klength}\\
\col{\hat{k}} & = \frac{\col{k}}{l_k} = \begin{bmatrix}\frac{k_1}{l_k}\\[.5em] \frac{k_2}{l_k}\\[.5em] \vdots\\[.5em] \frac{k_N}{l_k}\end{bmatrix} \\ \col{\hat{k}} & = \frac{\col{k}}{\abs{\col{k}}} = \begin{bmatrix}\frac{k_1}{\abs{\col{k}}}\\[.5em] \frac{k_2}{\abs{\col{k}}}\\[.5em] \vdots\\[.5em] \frac{k_N}{\abs{\col{k}}}\end{bmatrix} \\
\hat{a} & = al_k\\ \hat{a} & = a\abs{\col{k}}\\
\col{A} & = \hat{a}\hat{\col{k}} \col{A} & = \hat{a}\hat{\col{k}}
\end{align} \end{align}
Thus, the column $\hat{\col{k}}$ contains the \textsl{relative} absorption coefficients for our (three) discrete wavelengths such that the combined effect of the (three) absorption coefficients adds up to 1. Thus, the column $\hat{\col{k}}$ contains the \textsl{relative} absorption coefficients for our (three) discrete wavelengths such that the combined effect of the (three) absorption coefficients adds up to 1.
...@@ -164,8 +164,8 @@ You cannot just derive the (normalized) absorption vector $\hat{k}$ for your dy ...@@ -164,8 +164,8 @@ You cannot just derive the (normalized) absorption vector $\hat{k}$ for your dy
\end{equation} \end{equation}
The values in this column are the product of the amount of dye (a single number) and the (normalised) absorbance of the dye: $\hat{a}\hat{k}$. And therefore: The values in this column are the product of the amount of dye (a single number) and the (normalised) absorbance of the dye: $\hat{a}\hat{k}$. And therefore:
\begin{align} \begin{align}
l_d & = \sqrt{\hat{A}_R^2 + \hat{A}_G^2 + \hat{A}_B^2}\\ \abs{\col{\hat{A}}_d} & = \sqrt{\hat{A}_R^2 + \hat{A}_G^2 + \hat{A}_B^2} \label{Alength}\\
\col{\hat{k}}_d & = \frac{\col{\hat{A}}_d}{l_d} \label{kpixeld} \col{\hat{k}}_d & = \frac{\col{\hat{A}}_d}{\abs{\col{\hat{A}}_d}} \label{kpixeld}
\end{align} \end{align}
That is: the normalized column $\col{\hat{k}}$ equals the normalized column $\col{\hat{A}}$. That is: the normalized column $\col{\hat{k}}$ equals the normalized column $\col{\hat{A}}$.
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
\begin{savequote} \begin{savequote}
'Out Of Many, One People'\qauthor{Chairman of the Independence Celebrations Committee \cite{Sealy1962}} 'Out Of Many, One People'\qauthor{Chairman of the Independence Celebrations Committee \cite{Sealy1962}}
\end{savequote} \end{savequote}
\chapter{A flag and some discoveries}\label{testJA} \chapter{A flag and some lessons learned}\label{testJA}
In this chapter we will first benchmark the colour deconvolution algorithm (chapter \ref{algebra}) with a simple flag with three pure colours: yellow, green and black (no blends or `intensity' differences, figure \ref{flagJA}). In this chapter we will first benchmark the colour deconvolution algorithm (chapter \ref{algebra}) with a simple flag with three pure colours: yellow, green and black (no blends or `intensity' differences, figure \ref{flagJA}).
...@@ -28,17 +28,25 @@ I_\mathrm{max} = 2^8 - 1 = 255 ...@@ -28,17 +28,25 @@ I_\mathrm{max} = 2^8 - 1 = 255
We will first see whether we can deconvolve and reconstruct this simple three colour flag with the algorithm in chapter \ref{algebra}. We will first see whether we can deconvolve and reconstruct this simple three colour flag with the algorithm in chapter \ref{algebra}.
We will also write out (many) numerical examples in order to help the reader get some feeling and insight in what is going on.
\subsection{Deconvolution}
So, we calculate the absorbance columns of the three dyes (equation \ref{Apixeld}): So, we calculate the absorbance columns of the three dyes (equation \ref{Apixeld}):
\begin{align} \begin{align}
\col{\hat{A}}_y & = -\ln\left( \frac{ \col{R}_y + 1}{256} \right) = \begin{bmatrix} 0.0039\\ 0.1981 \\ 5.5452 \end{bmatrix} \label{JAAy}\\ \col{\hat{A}}_y & = -\ln\left( \frac{ \col{R}_y + 1}{256} \right) = \begin{bmatrix} 0.0039\\ 0.1981 \\ 5.5452 \end{bmatrix} \label{JAAy}\\
\col{\hat{A}}_g & = -\ln\left( \frac{ \col{R}_g + 1}{256} \right) = \begin{bmatrix} 5.5452 \\ 0.4953 \\ 1.4676 \end{bmatrix} \label{JAAg}\\ \col{\hat{A}}_g & = -\ln\left( \frac{ \col{R}_g + 1}{256} \right) = \begin{bmatrix} 5.5452 \\ 0.4953 \\ 1.4676 \end{bmatrix} \label{JAAg}\\
\col{\hat{A}}_b & = -\ln\left( \frac{ \col{R}_y + 1}{256} \right) = \begin{bmatrix} 5.5452 \\ 5.5452 \\ 5.5452 \end{bmatrix} \label{JAAb} \col{\hat{A}}_b & = -\ln\left( \frac{ \col{R}_y + 1}{256} \right) = \begin{bmatrix} 5.5452 \\ 5.5452 \\ 5.5452 \end{bmatrix} \label{JAAb}
\end{align} \end{align}
Normalise these columns (equation \ref{kpixeld}): Normalise these columns (equations \ref{Alength}--\ref{kpixeld}):
\begin{align} \begin{align}
\col{\hat{k}}_y & = \frac{\col{\hat{A}}_y}{\abs{\col{\hat{A}}_y}} = \begin{bmatrix} 0.0007 \\ 0.0357\\ 0.9994 \end{bmatrix} \label{JAky}\\ \abs{\col{\hat{A}}_y} & = 5.5487 \label{flagJAabsy} \\
\col{\hat{k}}_g & = \frac{\col{\hat{A}}_g}{\abs{\col{\hat{A}}_g}} = \begin{bmatrix} 0.9631 \\ 0.0860 \\ 0.2549 \end{bmatrix} \label{JAkg}\\ \col{\hat{k}}_y = \frac{\col{\hat{A}}_y}{\abs{\col{\hat{A}}_y}} & = \begin{bmatrix} 0.0007 \\ 0.0357\\ 0.9994 \end{bmatrix} \label{JAky}\\
\col{\hat{k}}_b & = \frac{\col{\hat{A}}_g}{\abs{\col{\hat{A}}_g}} = \begin{bmatrix} 0.5774 \\ 0.5774 \\ 0.5774 \end{bmatrix} \label{JAkb} \abs{\col{\hat{A}}_g} & = 5.7575 \label{flagJAabsg}\\
\col{\hat{k}}_g = \frac{\col{\hat{A}}_g}{\abs{\col{\hat{A}}_g}} & = \begin{bmatrix} 0.9631 \\ 0.0860 \\ 0.2549 \end{bmatrix} \label{JAkg}\\
\abs{\col{\hat{A}}_b} & = 9.6045\label{flagJAabsb}\\
\col{\hat{k}}_b = \frac{\col{\hat{A}}_g}{\abs{\col{\hat{A}}_g}} & = \begin{bmatrix} 0.5774 \\ 0.5774 \\ 0.5774 \end{bmatrix} \label{JAkb}
\end{align} \end{align}
And collect those three columns in a matrix (equation \ref{Kna}) and calculate the inverse: And collect those three columns in a matrix (equation \ref{Kna}) and calculate the inverse:
\begin{align} \begin{align}
...@@ -58,17 +66,97 @@ and next calculate the column with (estimated) dye contributions for this pixel ...@@ -58,17 +66,97 @@ and next calculate the column with (estimated) dye contributions for this pixel
\begin{equation} \begin{equation}
\col{\hat{a}}_p = \mat{K}^{-1} \col{A}_p = \begin{bmatrix} \hat{a}_y\\ \hat{a}_g \\ \hat{a}_b\end{bmatrix} \col{\hat{a}}_p = \mat{K}^{-1} \col{A}_p = \begin{bmatrix} \hat{a}_y\\ \hat{a}_g \\ \hat{a}_b\end{bmatrix}
\end{equation} \end{equation}
And for this simple pure colour flag, this works like a charm: we find either zero for the amount of dye or a constant value for the amount of dye, in the expected and appropriate places (figure \ref{flagJAas}).
\begin{figure} The actual (non-zero) value for `the amount of dye' that we find is precisely the length of that dye's absorption column $\abs{\col{\hat{A}}}$ before normalisation to $\col{\hat{k}}$. This is a direct result of that normalisation step: if we had \textsl{not} normalised the absorption columns before collecting them in $\mat{K}$, we would have found $a_y = a_g = a_b = 1$ for the non-zero dye amounts.\\
\begin{figure}[h]
\tiny \tiny
\subfloat[\label{flagJAay}]{% \subfloat[\label{flagJAay}]{%
\def\svgwidth{0.31\linewidth}\includesvg{pics/flagJAay}}\hfill \def\svgwidth{0.47\linewidth}\includesvg{pics/flagJAay}}\hfill
\subfloat[\label{flagJAag}]{% \subfloat[\label{flagJAag}]{%
\def\svgwidth{0.31\linewidth}\includesvg{pics/flagJAag}}\hfill \def\svgwidth{0.47\linewidth}\includesvg{pics/flagJAag}}\\
\subfloat[\label{flagJAab}]{% \subfloat[\label{flagJAab}]{%
\def\svgwidth{0.31\linewidth}\includesvg{pics/flagJAab}}\\ \def\svgwidth{0.47\linewidth}\includesvg{pics/flagJAab}}\\
\caption{Deconvolved (pure) due amounts for the three colour flag. \label{flagJAas}} \caption{Deconvolved (pure) dye amounts for the three colour flag. For each dye we find a binary image with dye amounts of either zero (black) or the length of the dye's absorption column $\abs{{\hat{A}}}$ before normalisation to ${\hat{k}}$ (white). With \textbf{(a)} the amount of yellow is either zero or $\hat{a}_y = 5.5487$ (equation \ref{flagJAabsy}); \textbf{(b)} the amount of green is either zero or $\hat{a}_g = 5.7575$ (equation \ref{flagJAabsg}); and \textbf{(c)} the amount of yellow is either zero or $\hat{a}_b = 9.6045$ (equation \ref{flagJAabsb}). \label{flagJAas}}
\end{figure}
\subsection{Reconstruction}
With the estimated dye amounts and the (normalised) absorption matrix $\mat{K}$, we can reconstruct `single dye RGB-images' and the complete three-dye RGB image. And if all went well, the reconstructed three-dye RGB image matches the original flag (near) perfectly (figure \ref{flagJArecon}).
\begin{figure}[h]
\subfloat[RGB = 252, 207, 0\label{flagJAy}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAy.png}}\hfill
\subfloat[RGB = 0, 154, 58\label{flagJAg}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAg.png}}\\
\subfloat[RGB = 0, 0, 0\label{flagJAb}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAb.png}}\hfill
\subfloat[\label{flagJAygb}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAygb.png}}\\
\caption{Reconstructed RGB images from estimated dye amounts and normalised dye absorption columns. With: \textbf{(a)} RGB image with contributions of only the yellow dye (equation \ref{flagJARy}); \textbf{(b)} RGB image with contributions of only the green dye (equation \ref{flagJARg}); \textbf{(c)} RGB image with contributions of only the black dye (equation \ref{flagJARy}); and \textbf{(d)} RGB image with contributions of all three dyes together (equation \ref{flagJARygb}). \label{flagJArecon}}
\end{figure}
We can reconstruct an RGB image for each pixel $p$ and for each dye $n$ with only the (absorption) contribution of that dye $n$ with (equations \ref{Apixela}--\ref{Ipixela}):
\begin{equation}
\col{\hat{R}}_{p, n} = \left(I_\mathrm{max} + 1\right)\mathrm{e}^{-\col{\hat{k}}_n\hat{a}_{p, n}} - 1
\end{equation}
So the reconstructed yellow pixels will have RGB values\footnote{Note that these values are rounded to integers, RGB doesn't know fractions.} (figure \ref{flagJAy}):
\begin{equation}
\col{\hat{R}}_y = 256\cdot \mathrm{e}^{-\begin{bmatrix} 0.0007 \\ 0.0357\\ 0.9994 \end{bmatrix} \cdot 5.5487} - 1 = \begin{bmatrix} 252 \\ 207 \\ 0 \end{bmatrix} \label{flagJARy}
\end{equation}
The reconstructed green pixels will have RGB values (figure \ref{flagJAg}):
\begin{equation}
\col{\hat{R}}_g = 256\cdot \mathrm{e}^{-\begin{bmatrix} 0.9631 \\ 0.0860 \\ 0.2549 \end{bmatrix} \cdot 5.7575} - 1 = \begin{bmatrix} 0 \\ 154 \\ 58 \end{bmatrix} \label{flagJARg}
\end{equation}
And the reconstructed black pixels will have RGB values (figure \ref{flagJAb}):
\begin{equation}
\col{\hat{R}}_b = 256\cdot \mathrm{e}^{-\begin{bmatrix} 0.5774 \\ 0.5774 \\ 0.5774 \end{bmatrix} \cdot 9.6045} - 1 = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \label{flagJARb}
\end{equation}
If we wish to reconstruct the original flag with all three colours combined, we can \textsl{not} add the RGB images in figures \ref{flagJAy}--\ref{flagJAb}: that would be additive colour mixing for dyes that exhibit subtractive colour mixing. This is also obvious from adding e.g.\ the red values of the three dyes for the upper left (yellow) pixel: 252+255+255 (because that pixel in white in figures \ref{flagJAg}--\ref{flagJAb}) is not a `proper' RGB value.
Instead we have to combine the \textsl{absorption} contributions of the dyes and use:
\begin{equation}
\col{\hat{R}}_p = \left(I_\mathrm{max} + 1\right)\mathrm{e}^{-\mat{K}\col{\hat{a}}} - 1 \label{flagJARygb}
\end{equation}
By multiplying the complete dye absorption matrix $\mat{K}$ with the complete column of (estimated) dye amounts $\col{\hat{a}}$ we ensure that all dye absorption effects are properly combined and accounted for (figure \ref{flagJAygb}). Note, by the way, that all RGB values in the fully reconstructed image (figure \ref{flagJAygb}) are a perfect match with the RGB values in the original flag (figure \ref{flagJA}).
\begin{figure}[h!]
\subfloat[RGB = 253, 244, 93\label{flagJAyn}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAyn.png}}\hfill
\subfloat[RGB = 96, 232, 196\label{flagJAgn}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAgn.png}}\\
\subfloat[RGB = 142, 142, 142\label{flagJAbn}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAbn.png}}\hfill
\subfloat[\label{flagJAygbn}]{%
\includegraphics[width=0.47\linewidth]{pics/flagJAygbn.png}}\\
\caption{Reconstructed RGB images with normalised dye absorption columns and $\hat{a} = 1$ for all dye amounts. With: \textbf{(a)} RGB image with contributions of only the normalised yellow dye; \textbf{(b)} RGB image with contributions of only the normalised green dye; \textbf{(c)} RGB image with contributions of only the normalised black dye; and \textbf{(d)} RGB image with contributions of all three normalised dyes together. \label{flagJAreconn}}
\end{figure} \end{figure}
\subsection{Normalisation}
The dye absorption columns that were calculated from the dyes RGB values all had a certain length, a certain `total absorption' over the three RGB channels. We normalised these absorption columns so that they all have a `total absorption' over the three RGB channels of 1. \\
\noindent As a result, deconvolution of these pure colours with the actual pure colours as input gives the length of the
dyes absorption column before normalisation as `the amount of dye'. \\
We did not ask
`how much of this yellow with total absorption $\abs{\col{A}_y}$ is present';\\
we asked
`how much of this yellow with total absorption of 1 is present'.\\
And that answer is, of course, $\abs{\col{A}_y}$.\\
\noindent So, this normalisation makes sure that the amount of dye that you used to obtain your reference `dye colours' (RGB) or reference dye absorption columns, does not affect the outcome of the analysis (the estimated dye amounts).
We therefore did not deconvolve with the actual colours (absorptions) in the flag, but with `normalised colours' (absorptions) of the flag (figure \ref{flagJAreconn}).
\section{Two colour flag analysis} \section{Two colour flag analysis}
......
...@@ -26,70 +26,64 @@ for d = 1:3 ...@@ -26,70 +26,64 @@ for d = 1:3
end end
figure figure
imshow(uint8(P)) imshow(sum(amounts, 3)/norm(A(:, 3)))
colormap(gray)
colorbar
figure figure
imshow(uint8(Q)) imshow(uint8(round(P)))
% imwrite(uint8(P), 'pics/flagJAy.png', 'PNG')
figure figure
imshow(uint8(R)) imshow(uint8(round(Q)))
% imwrite(uint8(Q), 'pics/flagJAg.png', 'PNG')
figure figure
imshow(uint8(RGB)) imshow(uint8(round(R)))
% imwrite(uint8(R), 'pics/flagJAb.png', 'PNG')
% test image green/yellow figure
m = 3; imshow(uint8(round(RGB)))
[ay, ag] = meshgrid(linspace(0, norm(A(:, 1)), 256), linspace(0, norm(A(:, 2)), 256)); % imwrite(uint8(RGB), 'pics/flagJAygb.png', 'PNG')
tim = zeros([size(ag) 3]);
for r = 1:size(ag, 1)
for c = 1:size(ag, 2) % normalised deconvolution colours
tim(r, c, :) = cld_OD2RGB(K, [ay(r, c), ag(r, c), 0]'); p = 0; done = -1; Imax = 255;
for r = 1:size(im, 1)
for c = 1:size(im, 2)
p = p + 1; % update pixel count
if round(100*p/(numel(im)/3)) > done
% update feedback
done = round(100*p/(numel(im)/3));
fprintf('\b\b\b\b\b\b\b\b\b%3i%% done', done)
end
pam = squeeze(amounts(r, c, :));
pam = pam/max(pam);
% reconstruct pixels RGB contributions
[RGB(r, c, :), PQR] = cld_od2rgb(K, pam, Imax);
P(r, c, :) = PQR(:, 1);
Q(r, c, :) = PQR(:, 2);
R(r, c, :) = PQR(:, 3);
end end
end end
figure
imshow(uint8(tim))
[tamounts, tP, tQ, tR, tRGB, tA, tK, tiOD] = cld_decon(dyes(:, 1:2), tim, 'rgb');
figure figure
imagesc(ay) imshow(uint8(round(P)))
colorbar % imwrite(uint8(P), 'pics/flagJAyn.png', 'PNG')
figure figure
imagesc(squeeze(tamounts(:, :, 1))) imshow(uint8(round(Q)))
colorbar % imwrite(uint8(Q), 'pics/flagJAgn.png', 'PNG')
p = ay(:);
q = squeeze(tamounts(:, :, 1)); q = q(:);
figure figure
plot(p, q, '.', [0 6], [0 6]) imshow(uint8(round(R)))
% imwrite(uint8(R), 'pics/flagJAbn.png', 'PNG')
figure figure
imagesc(ag) imshow(uint8(round(RGB)))
colorbar % imwrite(uint8(RGB), 'pics/flagJAygbn.png', 'PNG')
figure
imagesc(squeeze(tamounts(:, :, 2)))
colorbar
p = ag(:);
q = squeeze(tamounts(:, :, 2)); q = q(:);
figure
plot(p, q, '.', [0 6], [0 6])
figure
imagesc(zeros(size(ay)))
colorbar
figure
imagesc(squeeze(tamounts(:, :, 3)))
colorbar
p = ay(:);
q = squeeze(tamounts(:, :, 3)); q = q(:);
figure
plot(p, q, '.', [0 6], [0 6])
clear all; close all;
% read image
im = imread('pics/Flag_of_Jamaica.png');
% image colours
dyes = [254 209 0; % yellow
0 155 58 % green
0 0 0]'; % black
% do the deconvolve
[amounts, P, Q, R, RGB, A, K, iOD] = cld_decon(dyes, im, 'rgb');
% show image
figure
imshow(im)
% test image green/yellow
m = 3;
[ay, ag] = meshgrid(linspace(0, norm(A(:, 1)), 256), linspace(0, norm(A(:, 2)), 256));
tim = zeros([size(ag) 3]);
for r = 1:size(ag, 1)
for c = 1:size(ag, 2)
tim(r, c, :) = cld_OD2RGB(K, [ay(r, c), ag(r, c), 0]');
end
end
figure
imshow(uint8(tim))
[tamounts, tP, tQ, tR, tRGB, tA, tK, tiOD] = cld_decon(dyes(:, 1:2), tim, 'rgb');
figure
imagesc(ay)
colorbar
figure
imagesc(squeeze(tamounts(:, :, 1)))
colorbar
p = ay(:);
q = squeeze(tamounts(:, :, 1)); q = q(:);
figure
plot(p, q, '.', [0 6], [0 6])
figure
imagesc(ag)
colorbar
figure
imagesc(squeeze(tamounts(:, :, 2)))
colorbar
p = ag(:);
q = squeeze(tamounts(:, :, 2)); q = q(:);
figure
plot(p, q, '.', [0 6], [0 6])
figure
imagesc(zeros(size(ay)))
colorbar
figure
imagesc(squeeze(tamounts(:, :, 3)))
colorbar
p = ay(:);
q = squeeze(tamounts(:, :, 3)); q = q(:);
figure
plot(p, q, '.', [0 6], [0 6])
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment