Kubelka-Munk Color Mixing in VEX

Color mixing in computer graphics often relies on simple linear blending in RGB space. Different approaches, like HSV blending, attempt to create more intuitive color transitions. However, none of the RGB mixing techniques accurately replicate how real-world materials mix pigments.

In the 1930s, Paul Kubelka and Franz Munk, two Austrian scientists, published a paper titled “Ein Beitrag zur Optik der Farbanstriche” (A Contribution to the Optics of Paint Layers).
The Kubelka-Munk (KM) model provides a physically based method for color mixing, particularly useful for rendering paints, and translucent materials. This article explores implementing the Kubelka-Munk model in Houdini’s VEX language, and some use cases. Take every explanation with a grain of salt though since I’m by no means an expert on Color spaces or mixing.

Why Use the Kubelka-Munk Model?

Consider what happens when you mix yellow and blue. Using standard RGB interpolation, you get a muddy grayish result. But mix yellow and blue paint in real life, and you get a vibrant green. In real-world pigment mixing (e.g., paint, ink, or dyes), colors interact based on absorption and scattering rather than simple averaging. The Kubelka-Munk model provides a more accurate representation by considering:

Absorption (K) – The amount of light absorbed by the material.
Scattering (S) – The amount of light scattered or reflected within the material.

This leads to more realistic pigment mixing, where combining yellow and blue, for example, produces a natural green instead of the desaturated result seen in RGB blending.

Simple example

SPH gif Note: this example is still skewed by houdinis linear interpolation between points

When you mix blue and yellow paint, you get green because of the way pigments absorb and reflect light. This follows subtractive color mixing, which applies to paints, inks, and dyes.

Blue pigment absorbs red and reflects blue and some green.
Yellow pigment absorbs blue and reflects red and green.
When you mix blue and yellow paint, the only common reflected color is green, since both pigments absorb different parts of the spectrum. This is why mixing blue and yellow pigments results in green.

Implementation

Our main focus is the absorption of different pigments, we will ignore the scattering component reducing the complexity while still achieving much more realistic results than standard RGB blending. We’ll convert RGB colors into their K absorption coefficients, mix them in the “spectral space” where all the magic happens, and then convert them back to RGB.

This implementation is based on nearpoints, it works best on particle systems, like fluids but you can basically plug it into anything.

Here is the function to get K absorption coefficients from RGB values

// Convert RGB to K coefficients
rgb_to_ks(vector rgb) {     
vector k;    
k.x = pow(1-rgb.x, 2)/(2*rgb.x);    
k.y = pow(1-rgb.y, 2)/(2*rgb.y);    
k.z = pow(1-rgb.z, 2)/(2*rgb.z);    

return k;    // S is assumed to be 1 for all channels 
}

and back to RGB

// Convert K coefficients to RGB
vector ks_to_rgb(vector k) {     
// Assuming S=1 for all channels    
vector rgb;    
rgb.x = 1 + k.x - sqrt(k.x * (k.x + 2));    
rgb.y = 1 + k.y - sqrt(k.y * (k.y + 2));    
rgb.z = 1 + k.z - sqrt(k.z * (k.z + 2));    

return rgb; 
}

These functions are then used to to mix them with a certain ratio. This ratio will later be determined by the distance of nearby points

// Mix two colors according to Kubelka-Munk theory
// c1, c2: RGB colors to mix
// ratio: mixing ratio (0-1), where 0 is 100% c1 and 1 is 100% c2
vector mix_km(vector c1; vector c2; float ratio) {
    // Convert RGB to K coefficients
    vector k1 = rgb_to_k(c1);
    vector k2 = rgb_to_k(c2);
    
    // Mix K coefficients linearly based on concentration
    vector k_mix = k1 * (1-ratio) + k2 * ratio;
    
    // Convert back to RGB
    return k_to_rgb(k_mix);
}

since we probably want to mix particles with multiple colors we should implement a function that gets all of the neighbor colors and their concentrations and calculates them at once.

// More advanced version with custom concentrations (for actual fluid simulation)
vector mix_km_concentrations(vector colors[]; float concentrations[]; int num_colors) {
    // Initialize absorption coefficients
    vector k_total = {0,0,0};
    float total_concentration = 0;
    
    // Sum K values weighted by concentrations
    for (int i = 0; i < num_colors; i++) {
        vector k = rgb_to_k(colors[i]);
        k_total += k * concentrations[i];
        total_concentration += concentrations[i];
    }
    
    // Normalize if needed
    if (total_concentration > 0) {
        k_total /= total_concentration;
    }
    
    // Convert back to RGB
    return k_to_rgb(k_total);
}

This is how I usually use the function, retrieving nearby points and their colors to compute a blended result based on the Inverse distance.

// get nearpoints and the amount
int npts[] = nearpoints(0,@P,chf("distance"),50);
i@neighbours = len(npts);

// setup array for color and concentration values
vector colors[];
float concentrations[];

// Add current point's color first
append(colors, @Cd);
// Give the current point's color a base concentration (adjust as needed)
append(concentrations, 1.0); 

// Exit early if no neighbors
if(i@neighbours == 0) {
    return;
}

// Get particle colors and densities
foreach(int point; npts) {
    vector color = vector(point(0, "Cd", point));
    append(colors, color);
    
    // Inverse distance weighting - closer points have higher concentrations
    float dist = distance(@P, point(0, "P", point));
    float weight = 1.0 - min(dist / chf("distance"), 1.0);
    append(concentrations, weight);
}

// Calculate mixed color using Kubelka-Munk
vector mixed_color = mix_km_concentrations(colors, concentrations, i@neighbours);

@Cd = mixed_color;

Complete Code

// Kubelka-Munk color mixing implementation for Houdini VEX

// Convert RGB color to K (absorption) coefficients
// Assumes S (scattering) = 1 for simplicity
vector rgb_to_k(vector rgb) {
    vector k;
    // Avoid division by zero with small epsilon
    float epsilon = 0.001;
    
    // Apply K-M formula to each channel
    k.x = pow(1-rgb.x, 2) / (2 * max(rgb.x, epsilon));
    k.y = pow(1-rgb.y, 2) / (2 * max(rgb.y, epsilon));
    k.z = pow(1-rgb.z, 2) / (2 * max(rgb.z, epsilon));
    
    return k;
}

// Convert K coefficients back to RGB
// Assumes S (scattering) = 1 for simplicity
vector k_to_rgb(vector k) {
    vector rgb;
    
    // Apply inverse K-M formula to each channel
    rgb.x = 1 + k.x - sqrt(k.x * (k.x + 2));
    rgb.y = 1 + k.y - sqrt(k.y * (k.y + 2));
    rgb.z = 1 + k.z - sqrt(k.z * (k.z + 2));
    
    // Clamp values to valid RGB range
    rgb = clamp(rgb, 0, 1);
    
    return rgb;
}

// Mix two colors according to Kubelka-Munk theory
// c1, c2: RGB colors to mix
// ratio: mixing ratio (0-1), where 0 is 100% c1 and 1 is 100% c2
vector mix_km(vector c1; vector c2; float ratio) {
    // Convert RGB to K coefficients
    vector k1 = rgb_to_k(c1);
    vector k2 = rgb_to_k(c2);
    
    // Mix K coefficients linearly based on concentration
    vector k_mix = k1 * (1-ratio) + k2 * ratio;
    
    // Convert back to RGB
    return k_to_rgb(k_mix);
}

// More advanced version with custom concentrations (for actual fluid simulation)
vector mix_km_concentrations(vector colors[]; float concentrations[]; int num_colors) {
    // Initialize absorption coefficients
    vector k_total = {0,0,0};
    float total_concentration = 0;
    
    // Sum K values weighted by concentrations
    for (int i = 0; i < num_colors; i++) {
        vector k = rgb_to_k(colors[i]);
        k_total += k * concentrations[i];
        total_concentration += concentrations[i];
    }
    
    // Normalize if needed
    if (total_concentration > 0) {
        k_total /= total_concentration;
    }
    
    // Convert back to RGB
    return k_to_rgb(k_total);
}



int npts[] = nearpoints(0,@P,chf("distance"),50);
i@neighbours = len(npts);

vector colors[];
float concentrations[];

append(colors, @Cd);
append(concentrations, 1);

// Exit early if no neighbors
if(i@neighbours == 0) {
    return;
}

// Get particle colors and densities
foreach(int point; npts) {
    vector color = vector(point(0, "Cd", point));
    append(colors, color);
    
    // Inverse distance weighting - closer points have higher concentrations
    float dist = distance(@P, point(0, "P", point));
    float weight = 1.0 - min(dist / chf("distance"), 1.0); // Linear falloff
    append(concentrations, weight);
}

// Calculate mixed color using Kubelka-Munk
vector mixed_color = mix_km_concentrations(colors, concentrations, i@neighbours);

// Assign the resulting color
@Cd = mixed_color;

Use Cases

Once you apply this to a grid with some nice colors and run it in a SOP solver, you’ll see the mixing in action. The interaction points between colors generate new shades based on their absorption values, creating natural looking blends.

Like I mentiond before, this method is especially useful for point-based fluid simulations, with FLIP, POP, or Vellum fluids being the most obvious applications. The video above showcases this technique in action using my custom SPH fluid solver, which I’ve detailed in this blog.

Here, you can see a comparison of a FLIP simulation: one without the color mixing solver, one where the solver is limited to high-velocity regions, and one where it only runs based on nearby points.

This implementation is highly simplified and could be improved by using real pigment data for more accurate color mixing. Tools like Mixbox, Rebelle, or spectral.js leverage spectral data and more advanced algorithms to better simulate how real-world pigments interact. However considering performance and complexity i feel like this approach strikes a nice balance between cool looking and easy to implement.