Commit 692a7713 authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

create dedicated functions and generalise

parent dc63abd0
function [RGB, PQR] = cld_OD2RGB(K, pam, Imax)
% [RGB, P, Q, R] = cld_OD2RGB(K, pam, Imax)
%
% Calculate RGB value from amounts and optical density with
% RGB = (Imax + 1) * exp(-K*pam) - 1
%
% Input:
% K 3x3 matrix with a column per dye or `colour' with normalised
% RGB absoption values.
% pam a column with 3 (deconvolved) dye or `colour' amounts
% Imax maximum intensity for the RGB values (default for a missing
% argument is 8-bits RGB: Imax = 2^8-1 = 255.
%
% Output:
% RGB RGB values with (absorption) contributions of all three dyes
% PQR RGB values with (absorption) contribution of dye or `colour' c
% only in column c
%
% (note that ClDlab uses a slightly different definition than in the
% <a href="https://imagej.net/Colour_Deconvolution">ImageJ plugin</a>)
%
% See also cld_RGB2OD, cld_decon
if nargin < 3, Imax = 255; end
if size(K, 2) < 3, K(3, 3) = 0; end
if numel(pam) < 3, pam(3) = 0; end
% reconstruct RGB with all dyes, set negative values to zero with subplus
RGB = subplus( (Imax+1)*exp(-K*pam) - 1);
if nargout > 1
% calculate individual dye RGBs
PQR = zeros(3);
for c = 1:3 % loop (backwards) over dyes
% RGB of this dye only in column c
PQR(:, c) = subplus( (Imax-1)*exp(-pam(c)*K(:, c)) - 1);
end
end
% part of ClDlab: https://gitlab.tue.nl/STEM/ClDlab
end
function OD = cld_RGB2OD(RGB, Imax)
% OD = cld_RGB2OD(RGB, Imax)
%
% Calculate optical density OD from RGB with
% OD = -log( (RGB + 1)/(Imax + 1) )
%
% (note that ClDlab uses a slightly different definition than in the <a
% href="https://imagej.net/Colour_Deconvolution">ImageJ plugin</a>)
%
% Imax = 2^8-1 = 255 by default for a missing Imax argument (8-bits RGB).
% Calculations are performed on double(RGB) when RGB is not a float
%
% See also cld_OD2RGB, cld_decon
if nargin < 2, Imax = 255; end
if ~isfloat(RGB), RGB = double(RGB); end
% add one to avoid taking log(0) (which would be bad)
OD = -log( (RGB + 1)/(Imax + 1) );
% part of ClDlab: https://gitlab.tue.nl/STEM/ClDlab
end
function [amounts, P, Q, R, RGB, A, K, iOD] = cld_decon(dyes, im, rod)
% [amounts, P, Q, R, RGB, A, K, iOD] = cld_decon(dyes, im, 'rod')
%
% Input:
% dyes matrix with dye RGB or OD values, a column [R G B]^T for each
% of either two or three dyes (3x2 or 3x3 matrix)
% im RGB image to deconvolve with dyes
% rod (optional) a string to identify whether dyes is in RGB or OD
% rod = 'RGB' dyes contains dye RGB values (default)
% rod = 'OD' dyes contains dye OD (optical density) values
%
% Output:
% amounts 3D matrix of size(im) with dye amounts for dye n in
% amounts(:, :, n)
% P RGB image with the RGB contribution to the final image of
% the first dye (note: double!)
% Q RGB image with the RGB contribution to the final image of
% the second dye (note: double!)
% R RGB image with the RGB contribution to the final image of
% the third dye (note: double!)
% RGB complete reconstructed RGB image from deconvolved dye
% amounts and input colours in dyes (P+Q+R, double!)
% A OD matrix for dyes (3x3)
% K normalised OD matrix for dyes (3x3)
% iOD OD matrix of size(im) for im (includes added third dye)
%
% When dyes contains only two columns (for two dyes), a third dye will be
% created with an OD vector perpendicluar to those of the two input dyes.
%
% Note that ClDlab uses OD = -log( (RGB + 1)/(Imax + 1) ); which differs
% slightly from OD = -log( (RGB + 1)/Imax ); as it is used in the <a
% href="https://imagej.net/Colour_Deconvolution">ImageJ plugin</a>
%
% See also cld_RGB2OD, cld_OD2RGB
% parse input
if nargin < 3
rod = 'RGB';
end
% initialise
amounts = zeros(size(im));
P = amounts; Q = P; R = P; RGB = P;
Imax = 255;
% RGB or OD?
if strcmpi(rod(1), 'r')
% dyes contains RGB values, convert to OD
A = cld_RGB2OD(dyes, Imax);
else
A = dyes;
end
% (only) two colours?
if size(A, 2) < 3
% add third `perpendicular' colour
A(:, 3) = cross(A(:, 1), A(:, 2));
end
% normalise: matrix K contains `k-hat'
K = A;
for c = 1:3
K(:, c) = K(:, c)/norm(K(:, c));
end
% convert image RGB to image OD
iOD = cld_RGB2OD(double(im), Imax);
for r = 1:size(im, 1)
for c = 1:size(im, 2)
% get pixels optical density
pOD = squeeze(iOD(r, c, :));
% pixels dye contributions: 3 by 1 column
% pixel amounts = inv(K)*pOD
pam = K\pOD; % is faster in Matlab
% add dye amounts for this pixel to 3D matrix
amounts(r, c, :) = 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
% part of ClDlab: https://gitlab.tue.nl/STEM/ClDlab
end
\ No newline at end of file
function [amounts, P, Q, R, RGB, A, K, iOD] = clrdecon(dyes, im)
% dyes: 3by2 or 3by3 matrix with a column with RGB values for each dye
if size(dyes, 2) < 3
dyes(3, 3) = 0;
adddye = 1;
else
adddye = 0;
end
amounts = zeros(size(im));
P = amounts; Q = P; R = P; RGB = P; iOD = RGB;
% dye OD columns: a*k
A = -log( (dyes + 1)/256 ); % add one to avoid taking log(0) (which would be bad)
if adddye > 0
% replace third colour by a perpendicluar one
A(:, 3) = cross(A(:, 1), A(:, 2));
end
% normalise: matrix K contains `k-hat'
K = A;
for c = 1:3
K(:, c) = K(:, c)/norm(K(:, c));
end
for r = 1:size(im, 1)
for c = 1:size(im, 2)
% pixel RGB values: 3 by 1 column
Ip = double(squeeze(im(r, c, :)));
% pixel OD
Ap = -log( (Ip + 1)/256 );
% add to image optical density
iOD(r, c, :) = Ap;
% dye contributions: 3 by 1 column
% amounts = inv(K)*Ap
ap = K\Ap; % is faster in Matlab
% ap = subplus(ap) % set negative amounts to zero
amounts(r, c, :) = ap;
% reconstruct pixels RGB contributions
P(r, c, :) = 256*exp(-ap(1)*K(:, 1))-1;
Q(r, c, :) = 256*exp(-ap(2)*K(:, 2))-1;
R(r, c, :) = 256*exp(-ap(3)*K(:, 3))-1;
RGB(r, c, :) = subplus(256*exp(-K*ap)-1); % set negative amounts to zero
end
end
end
\ No newline at end of file
......@@ -5,15 +5,9 @@ dyes = [254 209 0; % yellow
0 155 58 % green
0 0 0]'; % black
im = imread('pics/Flag_of_Jamaica.png');
[amounts, P, Q, R, RGB, A, K, iOD] = clrdecon(dyes, im);
[amounts, P, Q, R, RGB, A, K, iOD] = cld_decon(dyes, im);
figure
quiver3([0 0 0], [0 0 0], [0 0 0], K(1, :), K(2, :), K(3, :))
hold on
[vec, val] = eig(K)
quiver3([0 0 0], [0 0 0], [0 0 0], vec(1, :), vec(2, :), vec(3, :), 'r')
figure
imshow(im)
......
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