{"id":5510,"date":"2022-05-03T16:15:43","date_gmt":"2022-05-03T20:15:43","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=5510"},"modified":"2022-05-05T07:51:09","modified_gmt":"2022-05-05T11:51:09","slug":"more-experiments-with-srgb-gamut-boundary-in-lab-space","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2022\/05\/03\/more-experiments-with-srgb-gamut-boundary-in-lab-space\/","title":{"rendered":"More experiments with sRGB gamut boundary in L*a*b* space"},"content":{"rendered":"
I'm still playing around with RGB gamut calculations in <\/span>$ L^* a^* b^* $<\/span> space. (See my last post on this topic, <\/span>\"Visualizing out-of-gamut colors in a Lab curve.\"<\/span><\/a>) Today's post features some new gamut-related visualizations, plus some computational tricks involving gamut boundaries and rays in <\/span>$ L^* a^* b^* $<\/span> space.<\/span><\/div>
First, here is another way to communicate the idea that the in-gamut area in the <\/span>$ (a^*,b^*) $<\/span> plane varies with <\/span>$ L^* $<\/span> (perceptual lightness). For 9 values of <\/span>$ L^* $<\/span>, (10, 20, ..., 90), I'll compute a 2-D <\/span>$ (a^*,b^*) $<\/span> gamut mask by brute-forcing it. Then, I'll use overlaid contour plots to show the variation in gamut boundaries.<\/span><\/div>
a = -110:.1:110;<\/span><\/span><\/div><\/div>
b = -110:.1:110;<\/span><\/span><\/div><\/div>
L = 10:20:90;<\/span><\/span><\/div><\/div>
[aa,bb,LL] = meshgrid(a,b,L);<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
for <\/span>k = 1:length(L)<\/span><\/span><\/div><\/div>
rgb = lab2rgb(cat(3,LL(:,:,k),aa(:,:,k),bb(:,:,k)));<\/span><\/span><\/div><\/div>
mask = all((0 <= rgb) & (rgb <= 1),3) * 2 - 1 + L(k);<\/span><\/span><\/div><\/div>
contour(a,b,mask,[L(k) L(k)],LineColor=[.8 .8 .8],LineWidth=1.5,ShowText=true,<\/span>...<\/span><\/span><\/div><\/div>
LabelSpacing=288)<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
axis <\/span>equal<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
box <\/span>on<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"a*\"<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>\"b*\"<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>\"Gamut boundary in the (a,b) plane for several values of L*\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Here's another visualization concept. People often show colors on the <\/span>$ (a^*,b^*) $<\/span> plane, to give an idea of the meaning of <\/span>$ a^* $<\/span> and <\/span>$ b^* $<\/span>, but that doesn't communicate very well the idea that there are usually multiple colors, corresponding to various <\/span>$ L^* $<\/span> values, at any one <\/span>$ (a^*,b^*) $<\/span> location. Below, I show both the brightest in-gamut color and the darkest in-gamut color at each <\/span>$ (a^*,b^*) $<\/span> location.<\/span><\/div>
a = -110:110;<\/span><\/span><\/div><\/div>
b = -110:110;<\/span><\/span><\/div><\/div>
[aa,bb] = meshgrid(a,b);<\/span><\/span><\/div><\/div>
L_max = zeros(size(aa));<\/span><\/span><\/div><\/div>
L_min = zeros(size(aa));<\/span><\/span><\/div><\/div>
for <\/span>p = 1:size(aa,1)<\/span><\/span><\/div><\/div>
<\/span>for <\/span>q = 1:size(bb,1)<\/span><\/span><\/div><\/div>
[L_min(p,q),L_max(p,q)] = Lrange(aa(p,q),bb(p,q));<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
 <\/div><\/div><\/div>
rgb = lab2rgb(cat(3,L_max,aa,bb));<\/span><\/span><\/div><\/div>
figure<\/span><\/span><\/div><\/div>
tiledlayout(1,2)<\/span><\/span><\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
imshow(rgb,XData=a([1 end]),YData=b([1 end]))<\/span><\/span><\/div><\/div>
axis <\/span>xy<\/span><\/span><\/div><\/div>
axis <\/span>on<\/span><\/span><\/div><\/div>
xlabel <\/span>a<\/span><\/span><\/div><\/div>
ylabel <\/span>b<\/span><\/span><\/div><\/div>
title(<\/span>\"Brightest in-gamut color\"<\/span>)<\/span><\/span><\/div><\/div>
rgb_min = lab2rgb(cat(3,L_min,aa,bb));<\/span><\/span><\/div><\/div>
nexttile<\/span><\/span><\/div><\/div>
imshow(rgb_min,XData=a([1 end]),YData=b([1 end]))<\/span><\/span><\/div><\/div>
axis <\/span>on<\/span><\/span><\/div><\/div>
axis <\/span>xy<\/span><\/span><\/div><\/div>
xlabel <\/span>a<\/span><\/span><\/div><\/div>
ylabel <\/span>b<\/span><\/span><\/div><\/div>
title(<\/span>\"Darkest in-gamut color\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Next, I find myself sometimes wanting to draw a ray in <\/span>$ L^* a^* b^* $<\/span> space and find the gamut boundary location along that ray. To that end, I wrote a simple utility function (<\/span>findNonzeroBoundary<\/span> below) that performs a binary search to find where a function goes from positive to 0. Then, I wrote some anonymous functions to find the desired gamut boundary point. Specifically, I was interested in this question: For a given <\/span>$ L^* $<\/span> value and a given <\/span>$ (a^*,b^*) $<\/span> plane angle, <\/span>h<\/span>, what is the in-gamut color with the maximum chroma, <\/span>c<\/span>, or distance from <\/span>$ (0,0) $<\/span> in the <\/span>$ (a^*,b^*) $<\/span> plane?<\/span><\/div>
Fair warning<\/span>: the code below gets tricky with anonymous functions. You might hate it. If so, I totally understand, and I hope you'll forgive me. :-)<\/span><\/div>
I'll start by creating an anonymous function that converts from <\/span>$ L^* c h $<\/span> to <\/span>$ L^* a^* b^* $<\/span>:<\/span><\/div>
lch2lab = @(lch) [lch(1) lch(2)*cosd(lch(3)) lch(2)*sind(lch(3))];<\/span><\/span><\/div><\/div><\/div>
Next, here is an anonymous function that returns whether or not a particular <\/span>$ L^* a^* b^* $<\/span> point is in gamut.<\/span><\/div>
inGamut = @(lab) all(0 <= lab2rgb(lab),2) & all(lab2rgb(lab) <= 1,2);<\/span><\/span><\/div><\/div><\/div>
Finally, a third anonymous function uses <\/span>findNonzeroBoundary<\/span> to find the gamut boundary point that I'm interested in.<\/span><\/div>
maxChromaAtLh = @(L,h) findNonzeroBoundary(@(c) inGamut(lch2lab([L c h])), 0, 200);<\/span><\/span><\/div><\/div><\/div>
Let's exercise this function to find a high chroma dark color at <\/span>$ h=0^{\\circ} $<\/span>.<\/span><\/div>
L = 35;<\/span><\/span><\/div><\/div>
h = 0;<\/span><\/span><\/div><\/div>
c = maxChromaAtLh(L,h)<\/span><\/span><\/div>
c = 61.0229<\/div><\/div><\/div><\/div>
And here's what that color looks like.<\/span><\/div>
rgb_out = lab2rgb(lch2lab([L c h]));<\/span><\/span><\/div>
rgb_out = <\/span>1\u00d73<\/span><\/div><\/div>
0.6414 0.0001 0.3316\r\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div>
figure<\/span><\/span><\/div><\/div>
colorSwatches(rgb_out)<\/span><\/span><\/div><\/div>
axis <\/span>equal<\/span><\/span><\/div><\/div>
axis <\/span>off<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
(The function <\/span>colorSwatches<\/span> is from <\/span>Digital Image Processing Using MATLAB<\/span><\/a>, and it is included in the <\/span>MATLAB code files for the book<\/span><\/a>.)<\/span><\/div>
What happens when we try to find a high chroma color, at the same hue angle, that is bright instead of dark?<\/span><\/div>
L = 90;<\/span><\/span><\/div><\/div>
h = 0;<\/span><\/span><\/div><\/div>
c = maxChromaAtLh(L,h)<\/span><\/span><\/div>
c = 14.8666<\/div><\/div><\/div><\/div>
You can can see that the maximum <\/span>c<\/span> value is much lower for the higher value of <\/span>$ L^* $<\/span>. What does it look like?<\/span><\/div>
rgb_out = lab2rgb(lch2lab([L c h]));<\/span><\/span><\/div><\/div>
figure<\/span><\/span><\/div><\/div>
colorSwatches(rgb_out)<\/span><\/span><\/div><\/div>
axis <\/span>equal<\/span><\/span><\/div><\/div>
axis <\/span>off<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
When I was doing these experiments to prepare for this blog post, my original intent was to just show examples for a couple of different values of <\/span>h<\/span> and <\/span>$ L^* $<\/span>. But I couldn't stop! It was too much fun, and I kept trying different values.<\/span><\/div>
After about 15 minutes or so, I decided it would be best to write some simple loops to generate a relatively large number of <\/span>$ (L^*,h) $<\/span> combinations to look at. Here's the code to generate the high <\/span>c<\/span> colors for a variety of <\/span>$ (L^*,h) $<\/span> combinations.<\/span><\/div>
dh = 30;<\/span><\/span><\/div><\/div>
h = -180:dh:150;<\/span><\/span><\/div><\/div>
L = 35:15:95;<\/span><\/span><\/div><\/div>
dL = 15;<\/span><\/span><\/div><\/div>
rgb = zeros(length(h),length(L),3);<\/span><\/span><\/div><\/div>
for <\/span>q = 1:length(L)<\/span><\/span><\/div><\/div>
<\/span>for <\/span>k = 1:length(h)<\/span><\/span><\/div><\/div>
c = maxChromaAtLh(L(q),h(k));<\/span><\/span><\/div><\/div>
rgb(k,q,:) = lab2rgb(lch2lab([L(q) c h(k)]));<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div><\/div>
And here's the code to view all those colors on a grid, with labeled <\/span>h<\/span> and <\/span>$ L^* $<\/span> axes.<\/span><\/div>
rgb2 = reshape(fliplr(rgb),[],3);<\/span><\/span><\/div><\/div>
p = colorSwatches(rgb2,[length(L) length(h)]);<\/span><\/span><\/div><\/div>
p.XData = (p.XData - 0.5) * (dh\/1.5) + h(1);<\/span><\/span><\/div><\/div>
p.YData = (p.YData - 0.5) * (dL\/1.5) + L(1);<\/span><\/span><\/div><\/div>
xticks(h);<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"h\"<\/span>)<\/span><\/span><\/div><\/div>
yticks(L)<\/span><\/span><\/div><\/div>
ylabel(<\/span>\"L*\"<\/span>)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
title(<\/span>\"Highest chroma (most saturated) colors for different L* and h values\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

Utility Functions<\/span><\/h3>
function <\/span>[L_min,L_max] = Lrange(a,b)<\/span><\/span><\/div><\/div>
<\/span>arguments<\/span><\/span><\/div><\/div>
a <\/span>(1,1) {mustBeFloat}<\/span><\/span><\/div><\/div>
b <\/span>(1,1) {mustBeFloat}<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
 <\/div><\/div>
L = 0:0.01:100;<\/span><\/span><\/div><\/div>
a = a*ones(size(L));<\/span><\/span><\/div><\/div>
b = b*ones(size(L));<\/span><\/span><\/div><\/div>
lab = [L ; a ; b]';<\/span><\/span><\/div><\/div>
rgb = lab2rgb(lab);<\/span><\/span><\/div><\/div>
gamut_mask = all((0 <= rgb) & (rgb <= 1),2);<\/span><\/span><\/div><\/div>
j = find(gamut_mask,1,<\/span>'first'<\/span>);<\/span><\/span><\/div><\/div>
k = find(gamut_mask,1,<\/span>'last'<\/span>);<\/span><\/span><\/div><\/div>
<\/span>if <\/span>isempty(j)<\/span><\/span><\/div><\/div>
L_min = NaN;<\/span><\/span><\/div><\/div>
L_max = NaN;<\/span><\/span><\/div><\/div>
<\/span>else<\/span><\/span><\/div><\/div>
L_min = L(j);<\/span><\/span><\/div><\/div>
L_max = L(k);<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
 <\/div><\/div>
function <\/span>x = findNonzeroBoundary(f,x1,x2,abstol)<\/span><\/span><\/div><\/div>
<\/span>arguments<\/span><\/span><\/div><\/div>
f <\/span>(1,1) function_handle<\/span><\/span><\/div><\/div>
x1 <\/span>(1,1) {mustBeFloat}<\/span><\/span><\/div><\/div>
x2 <\/span>(1,1) {mustBeFloat}<\/span><\/span><\/div><\/div>
abstol <\/span>(1,1) {mustBeFloat}<\/span> = 1e-4<\/span><\/span><\/div><\/div>
 <\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
 <\/div><\/div>
<\/span>if <\/span>(f(x1) == 0) || (f(x2) ~= 0)<\/span><\/span><\/div><\/div>
error(<\/span>\"Function must be nonzero at initial starting point and zero at initial ending point.\"<\/span>)<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
 <\/div><\/div>
xm = mean([x1 x2]);<\/span><\/span><\/div><\/div>
<\/span>if <\/span>abs(xm - x1) \/ max(abs(xm),abs(x1)) <= abstol<\/span><\/span><\/div><\/div>
x = x1;<\/span><\/span><\/div><\/div>
<\/span>elseif <\/span>(f(xm) == 0)<\/span><\/span><\/div><\/div>
x = findNonzeroBoundary(f,x1,xm);<\/span><\/span><\/div><\/div>
<\/span>else<\/span><\/span><\/div><\/div>
x = findNonzeroBoundary(f,xm,x2);<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div><\/div>\r\n<\/div>