How can I calculate the blending factor from a linear combination of two pictures?

Question

I have a sequence of images that are very similar and can be approximately represented as a combination of two 'ground truth' images with added noise, in the form:

Image_i ≈ x_i × Image(1) + (1−x_i) × Image(2) + Noise

where x_i is the mixing coefficient.

I'm looking for the best method in Python to determine the value of xx for each image in the sequence. The images are noisy and I expect x_i to range between 40% and 70%.

The images are relatively small, but I have thousands of them, so the solution should be computationally efficient.

What would be the best approach to solve this problem?

Edit: I added an image 3 which could be a potential image_i

Answer 1

there are many ways how you can do this.

Usually in case of normal errors you use a loss function that minimize root mean squared error (similar to how Linear regression is estimated).

Solution 1

After creating a loss function you can consider using fsolve to get an optimal x parameter.

import numpy as np
from scipy.optimize import fsolve

np.random.seed(88)

# --> Generating random input data

width, height = 640, 480
image_1 = np.random.normal(loc = 4, scale = 3, size = [ 640, 470])
image_2 = np.random.normal(loc = 8, scale = 3, size = [640, 470])
noise = np.random.normal(loc = 0, scale = 1, size = [640, 470])

# <-- End of random data generation


# sample unknown parameter
x_i = np.random.rand()
print(f"x_i: {x_i}") 
# x_i: 0.9129151984340369

# create a target image
image_new = x_i * image_1 + (1 - x_i) * image_2 + noise


def rmse_loss(x):
    """ 
        Root mean squared error loss
        loss = sqrt(sum_i[(y_pred[i] - y_true[i])**2])
    """

    estimation =  (x * image_1 + (1 - x) * image_2)
    return np.sqrt(np.sum(np.power(image_new  - estimation, 2)))

# for attribute `x0` you can use any value, it will converge to the solution eventually
res = fsolve(rmse_loss, x0 = 0)
print(res)
# [0.91293059]

Solution 2

Another option is to note that your equation is equivalent to the following:

image_new - image_2 = x_i * (image_1 - image_2) + noise

And this is exactly the definition of linear regression (b = Ax + noise), which can be solved using linalg.lstsq method:

from scipy import linalg

b = (image_new - image_2).reshape([-1,1])
a = (image_1 - image_2).reshape([-1,1])

res = linalg.lstsq(a,b)

print(res[0][0][0]) 
# 0.91302692

# how good is the fit: R2 value
r2 = 1 - res[1] / np.power(b,2).sum()

print(r2) # 0.965

Several caveats:

• I'm considering the situation where Noise have a zero mean;
• Methods like fsolve could be hard to scale if you have millions of pictures;