Properly Clamping Y'IQ Color Space?

So, I've been taking a look at the Y'IQ color space, and I've come across a rather obvious problem. The Y'IQ color space has a larger gamut than sRGB, so it's quite easy to have a color go out of range.

If I just clamp the resulting sRGB value to between 1 and 0 it seems to give somewhat strange colors. Y'IQ colors for Y' = 0 or Y' = 1 come out in a range of colors instead of just black and white. I've been trying to figure out a way of clamping just the IQ values to the nearest acceptable color but haven't had much luck.

I've been poking around the Internet a bit and seen some other people mention this, but I haven't seen any proposed solutions. I imagine it would probably give fairly acceptable results to clamp the IQ coordinates to a simpler shape (that surrounds the gamut) and then clamp the resulting sRGB to between 0 and 1, but I'm wondering if anyone has any clue what a more mathematically accurate method might be.

If anyone's wondering why I'm even bothering with this... I've been doing a lot of research on color spaces and got a bit curious about this. So it's no big deal if I can't actually figure out a way of doing this... if anything it might be more useful to have a color space where I can just make all the colors go weird just by changing the Luma....

 

I've come up with a rather useful solution (and no I haven't been working on this nonstop since the last time posted).

Code (csharp):

1.  
2. inline float3 YIQTolinRGB(float3 YIQ)
3. {
4.    //The luma needs to be clamped between 0 and 1 here to avoid artifacts when clamping the chrominance
5.    float luma = saturate(YIQ.x);
6.    float3 sRGB = mul(YIQTosRGB, float3(luma, YIQ.yz));
7.    //luma.xxx is equivalent to multiplying YIQTosRGB by float3(luma, 0.0f, 0.0f)
8.    return GammaToLinearSpace(ClampChrominanceToCube(sRGB, luma.xxx));
9. }
10.  

Code (csharp):

1.  
2. //ClampChrominanceToCube clamps the chrominance of a cubical colour space (like sRGB) to the nearest point on the cube in the direction of that color's achromatic point.
3. //It's useful for clamping chrominance in non-cubicle color spaces after being transformed into a cubicle color space.
4. inline float3 ClampChrominanceToCube(float3 ColorPos, float3 AchromaticPos)
5. {
6.    AABB3D box;//Three-dimensional axis aligned bounding box
7.    box.InitAABB3D(float3(0.0f, 0.0f, 0.0f), float3(1.0f, 1.0f, 1.0f));//Size of the cube
8.  
9.    Ray3D ray;//The ray is initialized with the position of the color and the direction of the achromatic point
10.    ray.InitRay3D(ColorPos, AchromaticPos - ColorPos);
11.  
12.    float3 hitPoint;
13.    ray.HitAABB3D(box, hitPoint);//In this situation the ray will always hit (except occasionally when it misses a corner, but the result still looks fine after using the saturate function)
14.  
15.    return saturate(lerp(ColorPos, hitPoint, any(ColorPos != saturate(ColorPos))));//Decides which color to use depending on whether ColorPos was outside or inside the cube to begin with
16. }
17.  

Code (csharp):

1.  
2. //Three-dimensional axis aligned bounding box
3. struct AABB3D
4. {
5.    float3 Bounds[2];
6.  
7.    inline void InitAABB3D(float3 aMinBound, float3 aMaxBound)
8.    {
9.        Bounds[0] = aMinBound;
10.        Bounds[1] = aMaxBound;
11.    }
12. };
13.  

Code (csharp):

1.  
2. struct Ray3D
3. {
4.    float3 Pos;//Origin of the ray
5.    float3 Dir;//Direction the ray pointing
6.    float3 InvDir;//Reciprocal of the direction
7.    bool3 Sign;//Sign of the reciprocal (with zero and negative numbers assigned false)
8.  
9.    inline void InitRay3D(float3 aPos, float3 aDir)
10.    {
11.        float3 normalDir = normalize(aDir + (aDir == 0.0f) * 0.00000000000001f);//Avoid normalizing float3(0.0f, 0.0f, 0.0f) and gets rid of division by zero errors
12.  
13.        Pos = aPos;
14.        Dir = normalDir;
15.        InvDir = 1.0f / normalDir;
16.        Sign = (InvDir < 0.0f);
17.    }
18.  
19.    //See where (if at all) the ray hits the AABB3D
20.    inline bool HitAABB3D(AABB3D aBox, out float3 HitPoint)
21.    {
22.        bool didHit;
23.        float3 tMin;//Time values for the closest plains hit
24.        float3 tMax;//Time values for the farthest plains hit
25.  
26.        //Generates time values for both X and Y planes
27.        tMin.x = (aBox.Bounds[Sign.x].x - Pos.x) * InvDir.x;
28.        tMax.x = (aBox.Bounds[!Sign.x].x - Pos.x) * InvDir.x;
29.        tMin.y = (aBox.Bounds[Sign.y].y - Pos.y) * InvDir.y;
30.        tMax.y = (aBox.Bounds[!Sign.y].y - Pos.y) * InvDir.y;
31.  
32.        didHit = !((tMin.x > tMax.y) + (tMin.y > tMax.x));//Check if they hit box
33.        tMin.x = lerp(tMin.x, tMin.y, (tMin.y > tMin.x));//The min X value should always be the closest
34.        tMax.x = lerp(tMax.x, tMax.y, (tMax.y < tMax.x));//The max X value should always be the farthest
35.  
36.        //Generates time values for both Z planes
37.        tMin.z = (aBox.Bounds[Sign.z].z - Pos.z) * InvDir.z;
38.        tMax.z = (aBox.Bounds[!Sign.z].z - Pos.z) * InvDir.z;
39.  
40.        didHit = didHit * !((tMin.x > tMax.z) + (tMin.z > tMax.x));//Check if they hit box
41.        tMin.x = lerp(tMin.x, tMin.z, (tMin.z > tMin.x));//The min X value should always be the closest
42.        //Making sure the max X value is still the farthest is unnecessary at this point
43.        HitPoint = Pos + Dir * tMin.x;//Generate the closest hit point
44.        return didHit;
45.    }
46. };
47.  

Here's a picture of the color solid without clamping:

And with clamping: