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 |

The lookup table data for CIELCHuv |

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 |

*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 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.