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 can be used to transform color spaces, as demonstrated in the following example: ImageMetalling 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 , among others. Essentially, these tools manipulate the image's color space. Hue/Saturation 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 or . 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. affine transformation nonlinear transformation on the B-Spline 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. 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. MLS 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 is the transformation matrix (like ), and is the . M CCM T translation matrix We can assume that for each variable in equation , we can find the minimum value for equation . (1) (2) Because the derivatives with respect to each of the variables in are zero for equation , we can directly solve for through the matrix in equation . The partial derivatives with respect to the free variables in form a linear system of equations. (2) T M (2) T Now, as you can see, we are able to express as a function of : T M Where and — are : weighted centroids This enables us to express equation as a function of : (2) M As you can see, it makes no difference to us which points we look for in the minimum of equation . Whether to search for points (1) or distances between weighted centroids and control points. Therefore, we will not be solving equation for the original control points but instead for new ones: and Considering equation , the solution to the problem is simplified to finding M from: (1) (3) Next, the problem is simplified by differentiating equation and determining 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 we first obtain: (4) M (4), Boom! That’s it, the required function has been found, use to substitute into : (5) (3) The only remaining task is to code the solution , as exemplified in the . 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. (6) ImageMetalling-16 example Expansion of optimization conditions Upon closer examination of , 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, can be expressed in a more complex form. For instance, consider transformation which permits uneven translations of the original field. Therefore, let's envision M as: At the same time M M similarity where and are connected through the similarity coefficient M’ M’’ 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 we can now write as: (1) Since we expressed M through and and both of them were found, then we can rewrite by solving , as we did in the previous step: (5) (7) Where In conclusion, the presents a method to link similarity transformation with the covariance matrix and to identify 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. original publication M However, this hypothesis still needs to be empirically tested, which you can do yourself. This can be accomplished by downloading and experimenting with the deformations in various spaces, and with different settings for finding . the example M Implementation We will write the main code of the solver in , taking into account the fact that the code will be executed both on the in Metal cores and on the ; for this, it will have to be wrapped in an 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. C++ GPU CPU Objective-C++ 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

Too Long; Didn't Read

Making security fun one episode at a time

Understanding Color Space Transform Using The Moving Least Squares Method

Too Long; Didn't Read

Denis Svinarchuk

About Author

TOPICS

THIS ARTICLE WAS FEATURED IN...

Understanding Color Space Transform Using The Moving Least Squares Method

Too Long; Didn't Read

Denis Svinarchuk

About Author

TOPICS

THIS ARTICLE WAS FEATURED IN...

RELATED STORIES