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.