The problem with finding boundary chroma in CIELAB

Previously I had mentioned the use of CIELCHab and HuSL for image blending and adjustment.  Whether using LCH or a HuSL method, my approach to these operations requires the knowledge of the maximum chroma of the sRGB gamut as a function of L and H.  This isn't a simple task in CIELAB.

The extent of sRGB in CIELUV (left) and CIELAB (right)

To see why it's problematic, let's start with an existing method.  The reference implementation of HuSL uses LCHuv (CIELUV) as the parent space instead of LCHab.  The idea is the same.  For the L and H of a given pixel, we can calculate the maximum chroma for the purpose of normalization and clamping.  The approach used here is to simply calculate six boundary curves representing the intersection of the RGB gamut faces with the plane of constant L containing the given pixel.  The rest is the calculations of line-line intersections and distance minimization.

function Cout=luvboundcalc(L,H)
    % this function uses the method used by the other implementations
    % of HuSL in LUV
    Axyz=[3.240454162114103 -1.537138512797715 -0.49853140955601; ... 
        -0.96926603050518 1.876010845446694 0.041556017530349; ...
        0.055643430959114 -0.20402591351675 1.057225188223179];
    ka=903.2962962962963;
    ep=0.0088564516790356308;
    ref=[0 1];
    Cout=ones(size(L))*realmax;

    Hradians=H*pi/180;
    sinH=sin(Hradians);
    cosH=cos(Hradians);
    sub1=(L+16).^3/1560896;
    mask=(sub1>ep);
    sub2=sub1;
    sub2(~mask)=L(~mask)/ka;
    
    for r=1:3;
        row=Axyz(r,:);
        a1=row(1);
        a2=row(2);
        a3=row(3);
        top=(11120499*a1 + 11700000*a2 + 12739311*a3)*sub2;
        rbot=9608480*a3 - 1921696*a2;
        lbot=1441272*a3 - 4323816*a1;
        bot=(rbot.*sinH + lbot.*cosH).*sub2;
        
        for k=1:2;
            C=L.*(top - 11700000*ref(k))./(bot + 1921696*sinH*ref(k));
            mask=(C>0 & C<Cout);
            Cout(mask)=C(mask);
        end
    end
    Cout=min(Cout,175.2);
end

Chroma map for CIELUV at L=75 showing linear boundary curves

Although one might not suspect it at first glance, the level curves of this odd shape are indeed straight lines, and the described method works well.  We can calculate the maximal chroma for every pixel and use it to normalize or truncate an image in whole, or we can precalculate a lookup table if it makes things quicker.

The lookup table data for CIELCHuv

Discussion of HuSL elsewhere, and my own idle curiosity had to question why there was no implementation in LCHab.  Maybe I don't know enough about why it would be a dumb idea, but I like the symmetry of having both options.  Besides, if I want to do LCHab>RGB conversions, I'd want such a bounding method for the purposes of truncation anyway; again, maybe I just don't know enough.

Very quickly after trying to produce my own HuSL method in LAB, I had to wonder if the choice of LUV was in part a matter of convenience.  Unlike LUV, the level curves and meridians of the projected RGB cube are not linear.  Could we break this problem up by angle intervals and create a piecewise solution?  That's a bit problematic, as the face intersections don't lie neatly on a curve of constant H.  It would be a game of finding interval endpoints which straddle the face intersections, calculating an excess of data and then minimizing.  I start to question whether this will be costly. 

Chroma map for CIELAB at L=75, showing nonlinear boundary curves

In the process of trying to shove the LAB>RGB conversions through a solver to find a boundary curve equation, I noticed something.  At points near the yellow corner, the boundary function can't be solved as a simple function of H and L.  In simple terms, at these points there exist multiple solutions for the R=1 curve, and contrary to the required solver behavior when handling solutions from multiple curves, we want the maximum solution.

Chroma map for CIELAB at L=97.14, corresponding to the yellow corner

Detail of the above map, for H=90-120 degrees

Radials cast from the neutral axis to the primary-secondary edges intersect the R=1 face

I have yet to come up with even a conditional piecewise analytic solution. Looking around, I have not found any other approach involving an analytic solution of this bounding problem.  I ended up using a bisection solver to calculate the boundary chroma, though this approach is still very slow and a naive approach has issues converging to the apex.  Due to the undercut, convergence at the very edge would require a infinitesimally small initial step size.

I ended up performing additional calculations in the ROI near the corner, solving the R=1, G=1, and B=0 faces independently and forming a logical combination to locate the exact edge.  By this method, the numeric solver can converge to the very edge even with a large initial step size.

The lookup table data for CIELCHab, note vertical face near the problem area

function Cout=labboundcalc(L,H)
    % LUV approach won't be simple in LAB, since level curves and meridians are all nonlinear
    % furthermore, yellow corner is actually undercut and convergence at the edge 
    broadroi=H>85 & H<114 & L>85;
    pastcorner=(H-5)<L;
    ROI=(broadroi & ~pastcorner);
    
    % find boundary for entire area
    % then calculate logical combination of faces near ROI
    tic
    Cout=labsolver(L,H,1);    
    if any(ROI)
        Cbg=labsolver(L(ROI),H(ROI),2);
        Cr=labsolver(L(ROI),H(ROI),3);

        edge=(Cbg>=Cr);
        Croi=Cout(ROI);
        Croi(edge)=Cbg(edge);
        Cout(ROI)=Croi;
    end
    toc
end

function Cout=labsolver(L,H,mode)
    % adapted bisection solver for LAB
    
    % initial boundary generation for LAB
    Lc0=[0 5 33 61 67 88 98 100];
    Cc0=[30 60 135 115 95 119 97 15]+5;
    Lc=linspace(0,100);
    Cc=interp1(Lc0,Cc0,Lc);
    ind=L/100*(length(Lc)-1);
    Lp=round(ind)+1;
    
    s=size(H);
    C=Cc(Lp);
    C=reshape(C,s);
        
    % initial step sizes
    cstep=10;
    stepsize=-ones(s)*cstep;
    
    limitdelta=1E-7*prod(s);
    lastCsum=abs(sum(sum(C)));
    unmatched=true;
    out=true(s);
    first=true;
    while unmatched
        % CONVERSION MUST PASS OOG VALUES
        % bypass gamma correction for speed (results converge at the faces)
        rgb=lch2rgb(cat(3,L,C,H),'lab','nogc');

        % is point in-gamut?
        wasout=out;
        switch mode
            case 1
                out=(any(rgb<0,3) | any(rgb>1,3));
            case 2
                out=rgb(:,:,3)<0 | rgb(:,:,2)>1;
            case 3
                out=rgb(:,:,1)>1 | C<cstep;
                if first
                    fout=out;
                    first=false;
                end
        end
        neg=C<0;
        big=C>140;
        out=out & ~neg;

        change=xor(wasout,out);
        stepsize(change)=-stepsize(change)/2;
        stepsize(big)=-abs(stepsize(big));
        stepsize(neg)=abs(stepsize(neg));

        C=C+stepsize;

        Csum=abs(sum(sum(C)));
        dC=abs(Csum-lastCsum)
        lastCsum=Csum;

        if dC<limitdelta 
            unmatched=false;
        end
    end
    Cout=max(C,0);
    
    if mode==3
        Cout(fout)=150;
    end
end

This is the approach used by the function maxchroma() in my Matlab image mangling toolbox.  As mentioned, the HuSL and LCH conversion tools are based on this function.

Create a uniform colormap for plots in Matlab

This is just a quick one, and I know that there are plenty of other methods for doing this.  I figured that if I'm going to demonstrate that the symmetric normalized HuSLp and HSYp color models can be used for image editing, I might as well show other uses as well.

The symmetric models HuSLp (CIELUV and CIELAB) and HSYp (YPbPr)

Recall the nonuniformity addressed by the rotationally-symmetric normalization of HuSLp and HSYp.  While HuSL and HSY are attempts to bring the convenience of models like HSL and HSV to other, more uniform color spaces (e.g. CIELAB, CIELUV), their normalization method distorts the relationship between S and C.  A path of constant normalized saturation now has a radius dependent on hue.  To regain the uniformity of the parent space, HuSLp and HSYp are normalized to the maximal rotationally symmetric subset of the projected RGB space.  This means that HuSLp and HSYp retain visual uniformity, but at the expense of maximum chroma range.

Surfaces of maximum S: HSL (top); HuSLab, HSY (middle); HuSLpab, HSYp (bottom)

Comparing HuSL and HSY to HSL, we see the difficulty of naive color selection by simply varying hue alone.  The low resolution of these gradients makes the brightness variation much easier to notice.  The HuSLp and HSYp maps appear as a series of uniform rows.  Color selection in these maps is limited, but greatly simplified.  While care can be taken to make colormaps of higher chroma from discontinuous paths in CIELAB or YPbPr, these generic image conversion methods can be used quite easily with adequate results.

In this example, we create three colormaps corresponding to the hue angles of the primary and secondary colors.  Note the brightness variation in HSL; anyone who has spent any time making plots in Matlab has probably found themselves avoiding certain named line colors due to their poor visibility ('b' on black, 'y' on white).  While notably less saturated, the HSYp and HuSLp colormaps yield lines which are all equally visible.  Here, the HuSLp method uses CIELUV simply because the saturated corners of the RGB space are not 60 degrees apart in CIELAB. Also, take note that I'm inverting colors in the plot() call and in the imwrite() call because I'm operating with an inverted X display.

steps=6;

H=0:360/steps:(360-360/steps);
K=ones(size(H));

hslset=permute(colorspace('<hsl',cat(3,H,K,K*0.5)),[2 3 1]);
hsypset=permute(hsy2rgb(cat(3,H,K,K*0.6),'pastel'),[2 3 1]);
huslpset=permute(husl2rgb(cat(3,H,K*100,K*65),'luvp','aligned'),[2 3 1]);

lw=2;
x=0:0.01:1.4;
os=2;
d=9;
sl=1.2;
subplot_tight(1,1,1);
for n=1:1:steps;
    plot(x,     (atan(x/(n/d))-sl*x),'color',1-hslset(n,:),'linewidth',lw); hold on;
    plot(x,    -(atan(x/(n/d))-sl*x),'color',1-hslset(n,:),'linewidth',lw); hold on;
    plot(x,1*os+(atan(x/(n/d))-sl*x),'color',1-hsypset(n,:),'linewidth',lw); hold on;
    plot(x,1*os-(atan(x/(n/d))-sl*x),'color',1-hsypset(n,:),'linewidth',lw); hold on;
    plot(x,2*os+(atan(x/(n/d))-sl*x),'color',1-huslpset(n,:),'linewidth',lw); hold on;
    plot(x,2*os-(atan(x/(n/d))-sl*x),'color',1-huslpset(n,:),'linewidth',lw); hold on;
end
view(-90, 90);
set(gca,'ydir','reverse');

to=-0.45;
text(1,-to  ,'HSL','fontsize',16,'fontweight','bold');
text(1,-to+1*os,'HSYp','fontsize',16,'fontweight','bold');
text(1,-to+2*os,'HuSLpuv','fontsize',16,'fontweight','bold');

frame=getframe;
imwrite(255-frame.cdata,'blogshit/colormaps1.png','png');

Plots using colormaps generated in HSL, HSYp, HuSLp

Certainly one has to weigh uniformity against the ability to discriminate plot colors.  Strategic selection of colors based on adjacency of plot objects is often unavoidable, and at some point, one has to ask how many things should even be in a plot before a different approach to visualization is more appropriate.  Perhaps these HuSLp/HSYp palettes would be even more appropriate for creating color themes for your favorite terminal emulator.

Maybe I'll do more to extend these tools toward making discontinuous color sets at higher chroma, but I think there are plenty of people working on that already.  Then again, it's not like I'm unaccustomed to reinventing wheels.  It's not like I have anything better to do, either.

More thoughts on the 'color' blend mode

In the previous post, I discussed the 'color' image blend mode as used by GIMP and my efforts in Matlab to make a more useful variation.  To recap, methods such as HSV, HSI, and HSL in general offer poor separation of color and brightness information.  Of these methods, HSL tends to be the most useful for color adjustment and blending due to the symmetry of its saturation space, despite the fact that L is a poorer representation of image brightness than I.  Unconstrained operation in other color spaces (LAB, LUV, YPbPr) leaves us with the problem of how to deal with out-of-gamut conditions on conversion back to RGB.  Calculating boundary chroma and performing data truncation prior to RGB conversion can solve the problems caused by out-of-gamut points, resulting in good retention of background brightness.

When I had originally approached the use of LCH for image blending operations, I began by using a constrained and normalized version of LCH known originally as HuSL (now called HSLuv).  For a speed increase I also used a normalized polar version of YPbPr I refer to as HSY.  In these variations, saturation is normalized to the chroma range of the projected RGB space.  This is convenient and solves most of the same issues which a constrained LCH method solves.  It offers a few conveniences and features, but it has some weaknesses as well.  The normalization of C causes colors of constant S to no longer have a uniformity of appearance, but repeated hue adjustments do not cause chroma compression as with LCH truncation.

The projection of RGB in CIELAB and YPbPr as the extent of HuSLab and HSY

One of the drawbacks of normalizing or bounding chroma to the extent of the projected RGB space in LCH is that only a fraction of the projected RGB space lies on continuous locus of constant chroma. There exists only a small subset of the RGB space wherein color points can be subject to arbitrary hue rotation without leaving the gamut and having their chroma changed by truncation.  This maximal rotationally-symmetric boundary defines the normalizing chroma used for HuSLp, a variation of HuSL wherein the uniformity of the parent space (LAB or LUV) is retained.  As HuSLp and HSYp are subsets of the projected RGB space, they obviously cannot specify the full range of RGB colors.  Specifically, they are restricted to colors near the neutral axis, hence the suffix denoting pastel.  Compared to HuSLp, HSYp is not as uniform, but it has access to a greater fraction of the RGB space than HuSLp in either CIELAB or CIELUV modes.

Maximal biconic subsets of RGB in CIELUV, CIELAB, YPbPr

Uniformity in HuSLp (LAB): surfaces of 100% S and 50% L

Uniformity in HSYp: surfaces of 100% S and 50% Y

While this makes HuSLp and HSYp useful for color selection tasks, they are also useful for image editing so long as the chroma limiting is acceptable. Consider a color blending task with highly saturated foreground.  With its larger color range, let us compare a color blend in HSYp and LCHab.

Extreme test foreground and background

Color blends in CIELCHab and HSYp

While the limitations of HSYp are noticeable, the benefits are far more striking in this case.  The blue, red, and green stripes are moderated and the image regains its depth. 

To be realistic, this use-case is impractical but for testing.  A foreground with moderated colors would be less influenced by the limited range of HSYp.  For a sufficiently desaturated foreground, HSYp blending would produce results similar to a LCH blend, but it may still be quite a bit faster.  The major convenience of HSYp color blending is simply that the chroma is always compressed to the maximal uniform boundary.  There is no need to tweak foreground saturation to prevent uneven results.

Obviously, using stripes makes it difficult to produce a pleasing blend in any case, but that's kind of the point when testing.  In general, image editing operations need to consider content as much as method.  Blending a broad range of saturated hues with a photographic scene tends to create results which defy a viewer's expectations of reality.  Seeing a person with a purple face or a wood table lit with rainbow light doesn't make any sense to a pattern-recognition machine built on earth.

The image manipulation toolbox I've posted contains HSYp and HuSLp conversion tools, as well as support for these methods in the related image blending function (imblend) and the color adjustment function (imtweak). 

Exploring the 'color' blend mode

In my ongoing adventures in creating image manipulation tools for Matlab, I have spent quite a bit of time getting desirable behavior from blending and adjustment tools.  Although most of these problems share a common set of core difficulties, perhaps one of the more intuitive ways to explore the issue is through the common "color" layer blending mode.

In image editing tools such as GIMP and Photoshop, image layers can be combined or blended using different pixel blending modes.  While modes such as the contrast enhancing modes don't suggest their operating methods so directly, certainly the component and arithmetic modes do.  Addition and division modes suggest the underlying math.  Hue and saturation suggest the channels to be transferred between foreground and background.  The detail-oriented observer might begin to raise questions at this point.  Hue? Saturation?  In what model are we working?

What about the 'color' mode?  Mere observation suggests that it's a combination of hue and saturation information being transferred, but again, the details are critically important. Various sources online suggest that this mode performs a HS transfer in HSV... or HSB... or HSL... or HSY?  One begins to wonder if anybody knows what they're talking about.

As my experience lies with GIMP, I based my expectations on its behavior.  Certainly, I can't test my blending algo against software I don't have; however, it might be worth keeping in mind that other suites may use different approaches.  The passing impression of online resources suggests that GIMP uses HSV or HSL, Photoshop uses CIELAB, and according to simplefilter.de, PaintShopPro uses some bastard HSL/HSY swap that has people so baffled that they ask questions of forums only to get more nonsensical guidance.   Let's sort this out as best we can by simple experiment.

I assume at this point a few things about what is desired.  Forgive my loose language hereon, but in this pixel blending mode, what should be transferred is the chroma information independent of the luminosity-related information.  In other words, the hue and colorfulness should change, but not the perceived brightness.  White should remain white, black should remain black.  Acronyms make it sound like a simple task, don't they?

Let's start with HSV, HSL, and HSI.  These are polar models used to express the content of the sRGB color space.  The convenience of these familiar models is simply the manner in which their polar format allows intuitive selection and adjustment of colors.

HSV is commonly visualized as a hexagonal cone or prism defined by piecewise math.  HSI is expressed as a cylinder or cone defined by simple trigonometry.  HSL is often described as a hexagonal bicone.  What's with all the different geometry?  All we want to do is separate HS from L/V/I, right?  Feel free to look at the math, but I think it's more valuable to the novice (me) to see how the isosurfaces associated with these models are projected into the sRGB space.

Isosurfaces of Value and Intensity as projected into the RGB cube

Value, the V in HSV, is defined as the the maximum of the RGB channels.  This means that the isosurfaces of V are cube shells in RGB.  Intensity (I in HSI) is the simple average of the channels.  This translates to planar isosurfaces.  At this point, it might start to look like HSI would provide more meaningful separation of color information from brightness information.  Does it make sense that [1 1 1] and [1 0 0] have the same brightness metric as they do in HSV?  Let's try to do an HS swap in these models and see what happens.

Foreground image

Background Image

Color blend in HSV

Color blend in HSI

Well, those look terrible.  White background content is not retained, and the image brightness is severely altered.  What happens with HSL?

Color blend in HSL

Well that looks better.  In fact, that looks exactly like GIMP behaves.  Checking the source code for GIMP 2.8.14, this is exactly how GIMP behaves. Let's look for a moment at why the white content is preserved.  While I previously compared the projected brightness isosurfaces to mislead the reader as I had misled myself, the answer lies in the projection of the saturation isosurfaces.  Simply put, it is the conditional nature of the calculation of S in HSL which creates its bicone geometry and allows comparable S resolution near both the black and white corners of the RGB cube.

Saturation isosurfaces as projected into the RGB cube from HSV, HSI, HSL

So we know how GIMP does it.  What are better ways?  While HSL helps prevent the loss of white image regions, the overall brightness still suffers.  Let's start with the vague suggestion that Photoshop uses CIELAB.  At this point, I can only make educated guesses at the internals, but intuition would suggest that we might try using CIELCHab, the polar expression of CIELAB.  In naivety, I perform a CH swap to obtain the following.

Color blend in CIELCHab

Now things are actually worse.  Both black and white regions are lost.  To understand why, I put together a tool for visualizing the projection of sRGB into other color spaces such as CIELAB.

The trajectory of an out-of-gamut point in unconstrained LCHab and HuSLab

Here, we see the volume containing all RGB triplets whose values lie within standard data ranges (e.g. 0-255 for uint8).  Unlike the case with HSV, HSL, and HSI, the sRGB space is a subset of the CIELAB space.  A color point which lies outside the depicted volume would have pixel values that are either negative or larger than allowed by the datatype.  Upon conversion back to RGB, these values would be truncated to the standard range. In this manner, color points outside the cube will be mapped to points on its surface. Since the truncation occurs after conversion to RGB, the trajectory of out-of-gamut color points will be parallel to the axes of the cube.  It then stands to reason that the orientation of the neutral axis means that out-of-gamut points tend to be mapped away from the white and black corners of the cube and toward the primary-secondary edges.

If the goal remains to isolate the effects of color manipulation and image brightness, the geometry of this projection poses an issue.  To transfer the radial position of a primary or secondary corner to a pixel which is near the extent of the L axis would be to push the pixel outside the cube.  If we were to rotate the hue of a pixel which lies on a primary or secondary corner, it would again move outside the cube.  How can we manage the existence and graceful conversion of these points?   To be honest, I'm not sure I know the best methods.  All I can offer are my current approaches.

One approach would be to do data truncation prior to RGB conversion.  Of course, the maximum chroma for a given L and H is not simple to calculate, but this would keep the trajectories of out-of-gamut points in a plane of constant L.  This produces good results, but reveals an issue with working methods.

Color blend in CIELCHab

Ideally, we would try to do all our image operations in the uniform space (LAB/LUV) and only perform the conversion to RGB for final output or screen display.  At the very least, we could convert to RGB with no data truncation, allowing for negative and supramaximal pixel values to retain the information stored by OOG points.  For my efforts in creating standalone tools with RGB input and output, this best approach seems terribly impractical.  When constrained by compatibility with Matlab's image viewer and other image handling functions, we're forced to decide whether it's appropriate to lose chroma information by truncating in LCH.  Even though it causes no ill effects for a single operation, it means that successive operations will tend to compress the chroma content of the image to a small subset of the projected RGB space.

This is more of a problem for image and color adjustment than for a single channel swapping blend operation, but it's something to keep in mind.  I'll probably come back to the topic of normalized versions of LCH at a later date.

If implementing a constrained color model like this is beyond the realm of practical effort-costs, consider this simple dirty method to retain brightness in a color blend.  Simply copy the background luma, perform a HS channel swap in HSL, then re-establish the original luma using your favorite luma-chroma transformation.  This might seem like it would waste a lot of time with converting twice, but it's often the fastest method among those I've described.  The results are similar to the normalized methods, and existing conversion tools can typically be used.

Color blend in HSL with Y preservation

For an example of the methods, feel free to check out the image blending function I've posted on the Mathworks file exchange. The conversion tools for LCHab, LCHuv, HuSLab, HuSLuv, and HSY are included in the parent submission.