Understanding Color Space Transform Using The Moving Least Squares Method

By utilizing simple mathematics and parallel computing, we can efficiently overcome performance constraints when manipulating large images. This article will explore how the geometric distortion features of ImageMetalling can be used to transform color spaces, as demonstrated in the following example:

To comprehend the code discussed in this article, it would be beneficial to have knowledge of the Color Lookup Table, geometric distortion, and how to work with pixels using the GPU.

Main idea

Up until now, to alter the color ratio in an image's pixels, we've employed a variety of intuitive and familiar tools such as curves, levels, and Hue/Saturation, among others. Essentially, these tools manipulate the image's color space.

What if we could replace all these tools with a single, definitive method? This method is unique because it can replace all other tools simultaneously. A method that's accurate due to its perceptual impact, aligning with the nature of human perception, rather than a literal or simplistic approach, like how Photoshop operates.

This is the method we will explore in this article. Let's see how we can manipulate not just the image itself but also its color representation. This will be done in a simple, mechanical way that is accessible to all.

Deformation of objects using the MLS method

There are various methods of image deformation, such as affine transformation or nonlinear transformation on the B-Spline. These are indeed fundamental and excellent tools. However, they often suffer from limitations such as the restricted number of control points that can be manipulated, or the excessively smooth extrapolation of the target distortions.

When dealing with color spaces, we may want to tightly associate a control point and its surroundings with their deformation. Furthermore, the set of points can be arbitrary. This issue has been effectively addressed for animation tasks in computer graphics, but it is still not widely used as a color correction tool in major photo editing software.

MLS is a method for reconstructing (approximating) an objective function on a set of known values using the values of the weighted mean square sum. This remarkable technique has many applications in regression analysis, statistics, Big Data, and deep learning.

Solving the problem

So, imagine we have N control points on a plane. To visualize this, think of a rubber surface scattered with balls that lay flat on a hard table. Each of these balls can be moved and fixed to a target point. Now, envision this surface colored with a slice of color space, for instance, RGB, HSV, or, for the sake of clarity and smoothness, Lab.

L=80

Then, with respect to affine transformation, we need to find a solution for the equation:

- root mean square weight,

- desired transformation function of point X,

- initial position of control points, 

- target position of control points

Let's minimize this: in our case, we need to find such a function:

in which the sum of all differences between points (where we arrived)

and the value of the function from the starting point (where we started)

— would be minimal, or in mathematical terms, tend to zero.

On the other hand, any affine transformation can be written through the equation:

where M is the transformation matrix (like CCM), and T is the translation matrix.

We can assume that for each variable in equation (1), we can find the minimum value for equation (2).

Because the derivatives with respect to each of the variables in are zero for equation (2), we can directly solve for T through the matrix M in equation (2). The partial derivatives with respect to the free variables in T form a linear system of equations.

Now, as you can see, we are able to express T as a function of M:

Where  and  — are weighted centroids:

This enables us to express equation (2) as a function of M:

 As you can see, it makes no difference to us which points we look for in the minimum of equation (1). Whether to search for points

or distances between weighted centroids and control points. Therefore, we will not be solving equation (1) for the original control points but instead for new ones: and  Considering equation (3), the solution to the problem is simplified to finding M from:

Next, the problem is simplified by differentiating equation (4) and determining M as a solution to a system of linear equations using inverse matrices. To find the inverse matrix, we must reduce the system to a solution that includes square matrices on both sides of the equation. And so from (4), we first obtain:

Boom!

That’s it, the required function has been found, use (5) to substitute into (3):

The only remaining task is to code the solution (6), as exemplified in the ImageMetalling-16 example. Given that there is no necessity for a system of equations solver, and that the matrices involved are merely 2x2, the entire mechanism should, in theory, operate fairly rapidly.

Expansion of optimization conditions

Upon closer examination of M, we can observe certain characteristics associated with it, as we have solved the problem by considering it as an affine transformation task, i.e., the task of scaling and shifting. It's clear that in reality, deformations are not confined to just these two parameters. Therefore, M can be expressed in a more complex form. For instance, consider similarity transformation which permits uneven translations of the original field. Therefore, let's envision M as: At the same time

where M’ and M’’ are connected through the similarity coefficient

with other vectors  and Since all  are similar to each other as plane vectors, we can write the similarity equality as an equation for the two-dimensional case:

Or rewrite it as: From which it follows that (1) we can now write as:

Since we expressed M through  and and both of them were found, then (5) we can rewrite by solving (7), as we did in the previous step:

Where

In conclusion, the original publication presents a method to link similarity transformation with the covariance matrix and to identify M for what's termed the rigid transformation. This approach allows for maximum deformation of spatial properties. At first glance, this method appears quite intuitive, and we anticipate it will be effective in deforming certain color spaces for our intended purposes.

However, this hypothesis still needs to be empirically tested, which you can do yourself. This can be accomplished by downloading the example and experimenting with the deformations in various spaces, and with different settings for finding M.

Implementation

We will write the main code of the solver in C++, taking into account the fact that the code will be executed both on the GPU in Metal cores and on the CPU; for this, it will have to be wrapped in an Objective-C++ extension. Obviously, running the solver code on the GPU will be beneficial if the size of the array of control points and the space that we will be modifying are large.

Code of kernels using the solver
kernel void kernel_mlsSolver(
                             constant float2  *input_points  [[buffer(0)]],
                             device float2    *output_points [[buffer(1)]],
                             constant float2  *p      [[buffer(2)]],
                             constant float2  *q      [[buffer(3)]],
                             constant int    &count   [[buffer(4)]],
                             constant MLSSolverKind  &kind [[buffer(5)]],
                             constant float    &alpha [[buffer(6)]],
                             uint gid [[thread_position_in_grid]]
                             ){
    ///
    /// The core of the solver calculates new coordinates for each point 
    ///
    float2 point = input_points[gid];
    IMPMLSSolver solver = IMPMLSSolver(point,p,q,count,kind,alpha);
    output_points[gid] = solver.value(point);
}

General MLS solver code

class IMPMLSSolver{
     
public:
         
    /**
     Create a solver for MNC
      
     @param point the current point relative to which we calculate the deformation 
     @param p array of initial control points
     @param q array of control target points
     @param count number of control points
     @param kind solver type: affine transformation, similarity transformation, transformation with maximum possible rigidity
     @param alpha - degree of deformation
     */
    IMPMLSSolver(const float2 point, 
#ifndef __METAL_VERSION__              
              const float2 *p, 
              const float2 *q,
#else
//
// It's important to understand that this is not just a different data type, but a different layout mechanism
// data in memory. In a specific solution for Metal Java, we pass arrays to
// GPU constant memory.                 
//                 
              constant float2 *p, 
              constant float2 *q,
#endif
              const int count, 
              const MLSSolverKind kind = MLSSolverKindAffine, 
              const float alpha = 1.0): 
    point_(point),
    kind_(kind), 
    alpha_(alpha), 
    count_(count),
    weight_(0),
    p_(p), q_(q)
    {
         
#ifndef __METAL_VERSION__
        //
        // For the assembly CPU, allocate temporary buffers for calculating intermediate values
        // solver
        //
        w_ = new float[count_];   
        pHat_ = new float2[count_];   
        qHat_ = new float2[count_];   
#endif
         
        // Finding the weights  
        solveW();
         
        // Finding a vector with a "header"
        solveHat();
         
        // Solving the main optimization equation
        solveM();  
         
#ifndef __METAL_VERSION__
        delete pHat_;
        delete qHat_;
        pHat_=qHat_=0;
#endif
         
    }    
         
    /**
     Returns the deformation point by coordinate
 
     @param point starting position
     @return new position 
     */
    float2 value(float2 point) {
        if (count_ <= 0) return point;   
        return (point - pStar_) * M + qStar_;
    }
     
    ~IMPMLSSolver() {
#ifndef __METAL_VERSION__
        delete w_;
        if (pHat_) delete pHat_;
        if (qHat_) delete qHat_;
#endif
    }
     
private:
     
    MLSSolverKind kind_;
    float   alpha_;
    int     count_;
    float2  point_;
#ifndef __METAL_VERSION__              
    const float2 *p_;
    const float2 *q_;
    float *w_;
    float2 *pHat_;
    float2 *qHat_;    
#else
    constant float2 *p_;
    constant float2 *q_;
    thread float w_[MLS_MAXIMUM_POINTS];
    thread float2 pHat_[MLS_MAXIMUM_POINTS];
    thread float2 qHat_[MLS_MAXIMUM_POINTS];
#endif
     
    float weight_;
    float2 pStar_;
    float2 qStar_;
    float mu_;
     
    float2x2 M;
     
    void solveW() {
         
        weight_ = 0;
        pStar_ = float2(0);
        qStar_ = float2(0);
 
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
             
            float d =  pow(distance(p_[i], point_), 2*alpha_);
             
            if (d < FLT_EPSILON)  d = FLT_EPSILON; 
             
            w_[i] = 1.0 / d;
            weight_ = weight_ + w_[i];
             
            pStar_ += w_[i] * p_[i];
            qStar_ += w_[i] * q_[i];
 
        }
         
        if (weight_ < FLT_EPSILON)  weight_ = FLT_EPSILON;
         
        pStar_ = pStar_ / weight_;                
        qStar_ = qStar_ / weight_;
    }  
     
    void solveHat(){        
         
        mu_ = 0;
         
        float _rmu1 = 0;
        float _rmu2 = 0;
         
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
             
            pHat_[i] = p_[i] - pStar_;                        
            qHat_[i] = q_[i] - qStar_;
             
            switch (kind_) {            
                case MLSSolverKindSimilarity:
                    mu_ += similarityMu(i);
                    break;                
                case MLSSolverKindRigid:
                    _rmu1 += rigidMu1(i); 
                    _rmu2 += rigidMu2(i);
                    break;
                default:
                    break;
            }
        }
         
        switch (kind_) {            
            case MLSSolverKindRigid:
                mu_ = sqrt(_rmu1*_rmu1 + _rmu2*_rmu2);
                break;
            default:
                break;
        }
         
        if (mu_ < FLT_EPSILON)  mu_ = FLT_EPSILON; 
         
        mu_ = 1/mu_;
    }
     
    void solveM() {
        switch (kind_) {
            case MLSSolverKindAffine:
                M = affineM();
                break;
            case MLSSolverKindSimilarity:
            case MLSSolverKindRigid:
                M = similarityM(point_);
                break;
        }
    }  
                    
    float2x2 affineM() {
        float2x2 mi = (float2x2){(float2){0,0},(float2){0,0}};
        float2x2 mj = (float2x2){(float2){0,0},(float2){0,0}};
 
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
             
            float2x2 pti({
                pHat_[i], 
                float2(0)                
            });
             
            float2x2 ppi({
                (float2){w_[i] * pHat_[i].x, 0.0}, 
                (float2){w_[i] * pHat_[i].y, 0.0}                
            });
             
            mi = mi + (float2x2)(pti * ppi);
             
             
            float2x2 ptj({
                w_[i] * pHat_[i], 
                float2(0)                
            });
             
            float2x2 qpj({
                (float2){qHat_[i].x, 0.0}, 
                (float2){qHat_[i].y, 0.0}                
            });
             
            mj = mj + (float2x2)(ptj * qpj);
        }        
        return inverse(mi) * mj;
    }
     
     
    float similarityMu(int i)  {
        return w_[i]*dot(pHat_[i], pHat_[i]);
    }
     
    float2x2 similarityM(float2 value) {
         
        float2x2 m = (float2x2){(float2){0,0},(float2){0,0}};
         
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
             
            float2 sp = -1 * __slashReflect(pHat_[i]);
            float2 sq = -1 * __slashReflect(qHat_[i]);
             
            float2x2 _p({
                (float2){w_[i] * pHat_[i].x, w_[i] * sp.x}, 
                (float2){w_[i] * pHat_[i].y, w_[i] * sp.y}                
            });
             
            float2x2 _q({qHat_[i], sq});
                         
            m = m + (float2x2)(_p * _q);
        }
        return  mu_ * m; 
    }
     
    float rigidMu1(int i) {
        return w_[i]*dot(qHat_[i], pHat_[i]);        
    }
     
    float rigidMu2(int i) {
        return w_[i]*dot(qHat_[i],  __slashReflect(pHat_[i]));        
    }
};

Examples

In the application, you can interact with the colors of the picture; the colors will be marked on the Lab color space plane: a-b (in the code, you can replace it with any other).

Colors that you want to fix and not change are also worth noting. Then, pull the color dot that you want to change to another. In this example, I highlighted the color of the bottle and chaotically marked a few more. As you can see, by significantly changing one color, the rest are practically unaffected, which can be very convenient for, for example, commercial photographers engaged in digital subject photography or for enthusiasts who dream of matching the color of a digital image to a film one.

And this is what the RGB cube looks like after our Lab space settings. For reference, the application also draws a hald-LUT and an HV plane from HSV space.

In this example, light color correction has been made: the cold tint of the original slide has been removed.

This is how the LUT of the original RGB cube is approximately transformed.

It is obvious that the presented framework for working with color spaces can become a universal color correction tool for future generations of digital photographers, when film finally leaves the scene, and color and sound design become a common practice.

Also published here.