Given a overlayed image (A^B) and one original image (A), how to calculate the other original image (B)

Question

Assume the following:

• Two images were overlayed over each other and merged into a single image (let's call it: image Y)
• I have one of the original images (let's call it: image A). This image is a half-transparent mask that was overlayed over the other image. But I don't know the mathematical operation. It was not just added but probably multiplied. Or maybe combined as an alpha composite.

How do I obtain the other original image (image B) given the other two images? All images are RGBA images.

I have tried to use PIL's PIL.ImageChops.difference() but this does not do the trick. The mathematical operation was not just a simple addition.

I also tried to play around with some code I found. But it certainly does not produce the right results.

# split RGBA images into 4 channels
rA, gA, bA, aA = pil_im.split()
rB, gB, bB, aB = mask_img.split()

# divide each channel (image1/image2)
rTmp = ImageMath.eval("int(a/((float(b)+1)/256))", a=rA, b=rB).convert('L')
gTmp = ImageMath.eval("int(a/((float(b)+1)/256))", a=gA, b=gB).convert('L')
bTmp = ImageMath.eval("int(a/((float(b)+1)/256))", a=bA, b=bB).convert('L')
aTmp = ImageMath.eval("int(a/((float(b)+1)/256))", a=aA, b=aB).convert('L')

# merge channels into RGBA image
imgOut = Image.merge("RGBA", (rTmp, gTmp, bTmp, aTmp))

• I would like to either reverse a multiplication (by using something like a division) or reverse an alpha compositing. I assume that the problem can be solved. It must just be something like linear algebra with f(A,B) = Y. I have A and Y, but what is B?
• What is the correct name of my problem?

Answer 1

If you have several linearly combined components, you can restore these components if you have the same number or more of mixtures (mixed with different coefficients).

For instance:

you have:

Mix1=I1*c1+I2*c2 and Mix2=I1*c3+I2*c4, where Mix1 and Mix2 are mixtures you have. I1 and I2 are source images, c1-c4 are mixing coefficients.

Using blind sourse separation you can restore I1 and I2.

For fast start with pyton version check here.

I tested it for CPP:

#include <cmath>
#define EIGEN_RUNTIME_NO_MALLOC // Define this symbol to enable runtime tests for allocations
#include <Eigen/Dense>
#include <vector>
#include "opencv2/opencv.hpp"
#include "opencv2/core/eigen.hpp"

using namespace Eigen;
using namespace cv;
using namespace std;
// --------------------------------------------------------
// 
// --------------------------------------------------------
class CFastICA
{
public:
    // --- МЕТОДЫ ----
    CFastICA();
    ~CFastICA();
    // Главный метод, раскладывающий сигнал по компонентам
    void apply(Mat& src,Mat& dst,Mat& W);
private:
    // --- Свойства ---
    // Матрица масштабных коэффициентов
    MatrixXd m_mixing_matrix;
    // Верхнее значение предела по кол-ву итераций алгоритма.
    int max_iter;
    // Малое значение, для задания допустимой погрешности алгоритма.
    double tol;

    static MatrixXd CFastICA::sym_decorrelation(MatrixXd m_mixing_matrix);
    static double CFastICA::gx(double x);
    static double CFastICA::g_x(double x);
    // Главный метод, раскладывающий сигнал по компонентам
    MatrixXd apply(MatrixXd& features);
}; 
// --------------------------------------------------------
// 
// --------------------------------------------------------
MatrixXd CFastICA::sym_decorrelation(MatrixXd m_mixing_matrix)
{
    MatrixXd K = m_mixing_matrix * m_mixing_matrix.transpose();
    SelfAdjointEigenSolver<MatrixXd> eig;
    eig.compute(K);
    return ((eig.eigenvectors() * eig.eigenvalues().cwiseSqrt().asDiagonal().inverse()) * eig.eigenvectors().transpose()) * m_mixing_matrix;
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
double alpha = 1.0; // alpha must be in range [1.0 - 2.0]

double CFastICA::gx(double x)
{
    return tanh(x*alpha);
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
double CFastICA::g_x(double x)
{
    return alpha * (1.0 - gx(x)*gx(x));
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
CFastICA::CFastICA() 
{
    max_iter = 2000;
    tol = 1e-10;
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
CFastICA::~CFastICA()
{

}
// --------------------------------------------------------
// 
// --------------------------------------------------------
void CFastICA::apply(Mat& src,Mat& dst,Mat& W)
{
    MatrixXd X=MatrixXd(src.rows,src.cols);
    // Переведем изображение из openCV в Eigen 
    cv2eigen(src,X);
    apply(X);
    eigen2cv(X,dst);
    eigen2cv(m_mixing_matrix,W);
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
Eigen::MatrixXd CFastICA::apply(MatrixXd& X)
{
    int n = X.rows();
    int p = X.cols();
    int m = n;

    // Whiten
    MatrixXd K;
    MatrixXd WX;

    VectorXd mean = (X.rowwise().sum() / (double)p);
    MatrixXd SPX = X.colwise() - mean;

    Eigen::JacobiSVD<MatrixXd> svd;
    svd.compute(SPX, Eigen::ComputeThinU);

    MatrixXd u = svd.matrixU();
    MatrixXd d = svd.singularValues();

    // see Hyvarinen (6.33) p.140
    K = u.transpose();
    for (int r = 0; r < K.rows(); r++)
    {
        K.row(r) /= d(r);
    }
    // see Hyvarinen (13.6) p.267 Here WX is white and data
    // in X has been projected onto a subspace by PCA
    WX = K * SPX;
    WX *= sqrt((double)p);

    cv::RNG rng;
    // Initial mixing matrix estimate
    if (m_mixing_matrix.rows() != m || m_mixing_matrix.cols() != m)
    {
        m_mixing_matrix = MatrixXd(m,m);
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < m; j++)
            {
                m_mixing_matrix(i,j) = rng.gaussian(1);
            }
        }
    }

    m_mixing_matrix = sym_decorrelation(m_mixing_matrix);

    int iter = 0;
    double lim = tol+1;
    while (lim > tol && iter < max_iter)
    {
        MatrixXd wtx = m_mixing_matrix * WX;
        MatrixXd gwtx  = wtx.unaryExpr(std::ptr_fun(&gx));
        MatrixXd g_wtx = wtx.unaryExpr(std::ptr_fun(&g_x));
        MatrixXd W1 = (gwtx * WX.transpose()) / (double)p - (g_wtx.rowwise().sum()/(double)p).asDiagonal() * m_mixing_matrix;
        W1 = sym_decorrelation(W1);
        lim = ((W1 * m_mixing_matrix.transpose()).diagonal().cwiseAbs().array() - 1.0).abs().maxCoeff();
        m_mixing_matrix = W1;
        iter++;
    }

    // Unmix
    m_mixing_matrix = (m_mixing_matrix*K);
    X = m_mixing_matrix * X;
    // set m_mixing_matrix
    m_mixing_matrix = m_mixing_matrix.inverse();
    return X;
}
// --------------------------------------------------------
// 
// --------------------------------------------------------
void main()
{
    int N_pts=1000;
    /*
    Mat X(2,N_pts,CV_64FC1);

    for(int i=0;i<N_pts;i++)
    {
        X.at<double>(0,i)=1.0*sin((double)i*0.1)+0.4*cos((double)i*0.01);
        X.at<double>(1,i)=0.4*sin((double)i*0.1)+1.0*cos((double)i*0.01);
    }
    */

    Mat I1=imread("D:\\ImagesForTest\\mandril.jpg",0);
    Mat I2=imread("D:\\ImagesForTest\\peppers.jpg",0);
    I1.convertTo(I1,CV_64FC1,1.0/255.0);
    I2.convertTo(I2,CV_64FC1,1.0/255.0);
    Mat Mix1=I1*0.2+I2*0.8;
    Mat Mix2=I1*0.3+I2*0.7;
    namedWindow("Mix1");
    namedWindow("Mix2");
    imshow("Mix1",Mix1);
    imshow("Mix2",Mix2);
    cv::waitKey(10);

    N_pts=I1.rows*I1.cols;

    Mat X(2,N_pts,CV_64FC1);
    Mix1=Mix1.reshape(1,1);
    Mix2=Mix2.reshape(1,1);
    Mix1.copyTo(X.row(0));
    Mix2.copyTo(X.row(1));

    CFastICA fica;

    Mat Y,W;
    fica.apply(X,Y,W);
    cout << W << endl;

    double M,m;
    cv::minMaxLoc(Y,&M,&m);

     cout << M << endl;
     cout << m << endl;
     cout << endl;

    // double scale=M-m;
    // Y=Y/scale;

    // cv::Mat img(400,N_pts,CV_8UC3);
    // img=0;
    // cv::line(img,cv::Point(0,200),cv::Point(N_pts,200),Scalar(255,255,255),1);
    // for(int i=0;i<N_pts-1;i++)
    // {
    //  cv::line(img,cv::Point(i,-Y.at<double>(0,i)*50+300),cv::Point(i,-Y.at<double>(0,i+1)*50+300),Scalar(0,255,0),1);
    //  cv::line(img,cv::Point(i,-Y.at<double>(1,i)*50+300),cv::Point(i,-Y.at<double>(1,i+1)*50+300),Scalar(255,0,0),1);    
    // 
    //  cv::line(img,cv::Point(i,-X.at<double>(0,i)*50+100),cv::Point(i,-X.at<double>(0,i+1)*50+100),Scalar(0,255,0),1);
    //  cv::line(img,cv::Point(i,-X.at<double>(1,i)*50+100),cv::Point(i,-X.at<double>(1,i+1)*50+100),Scalar(255,0,0),1);
    // }

    //namedWindow("img");

namedWindow("img1");
namedWindow("img2");
Mat res1(1,N_pts,CV_64FC1);
Mat res2(1,N_pts,CV_64FC1);
Y.row(0).copyTo(res1);
Y.row(1).copyTo(res2);
res1=res1.reshape(1,I1.rows);
res2=res2.reshape(1,I1.rows);
cv::normalize(res1,res1,0,1,cv::NORM_MINMAX);
cv::normalize(1.0-res2,res2,0,1,cv::NORM_MINMAX);
imshow("img1",res1);
imshow("img2",res2);
    cv::waitKey();

}

It takes as input 2 mixtures:

And as result we have 2 separated images.

Answer 2

It depends on how the images where merged. If it is a simple alpha blending Y = alpha * A + (1-alpha) * B, then you can use Andrey Smorodov's excellent answer given you have more than one overlayed image (which you probably don't have).
I don't think you can calculate B from the above equation as you have one equation with two unknowns (B and alpha). What you could do is de-blend the two images reversing the above formula and arbitrarily adjusting alpha according to visual impression:

import cv2 as cv

def on_trackbar(val):
    alpha = max(val / alpha_slider_max, 0.00001)
    beta = ( 1.0 - alpha )
    dst = cv.addWeighted(comb, 1/alpha, src2, -beta/alpha, 0.0)
    cv.imshow(title_window, dst)

# sample images from http://sipi.usc.edu/database/database.php?volume=misc
src1 = cv.imread('mandrill.tiff')
src2 = cv.imread('peppers.png')
# combination with alpha = 0.3
comb = cv.addWeighted(src1, 0.3, src2, 0.7, 0.0)

alpha_slider_max = 100
title_window = 'Linear De-Blend'

cv.namedWindow(title_window)
cv.createTrackbar('Alpha', title_window , 0, alpha_slider_max, on_trackbar)

on_trackbar(0)

When you adjust the alpha slider to 30 (alpha = 0.3) you'll see the original mandrill, for 100 the blended input image.

---

But you mentioned transparency: if you have an overlay A with transparent and opaque regions and combine this overlay with a background B so that Y = B for the transparent regions of A and Y = A otherwise (as in this SO question), then there's no chance reconstructing B for the opaque regions of A as B is copied into Y here with no share of A (Y = 1.0 * B + 0.0 * A).

Answer 3

If the images are blended linearly, and if we consider an ideal case where we don't lose precision due to the limited range of values a pixel can represent (say 8-bits in a single channel image), then, I think we can formulate the problem as follows (I think this should work, at least in theory, but please point out if there are any flaws):

If

A : image1 that is available to us.
B : image2 which we have to estimate.
Y = (1-c)A + c.B : blended image, where c is the blending coefficient, 0 <= c <= 1.

Y is available to us.

Then

Y   = A + c.(B-A)
Y-A = c.(B-A)

Now, we know Y-A as both Y and A are known to us. Let Y-A=C. Then

C = c.(B-A)

Let's now consider image pixels at location (x, y)

C(x,y) = c.(B(x,y)-A(x,y))

Let's consider two image pixels at locations (x1,y1) and (x2,y2)

C(x1,y1) = c.(B(x1,y1)-A(x1,y1))  - (1)
C(x2,y2) = c.(B(x2,y2)-A(x2,y2))  - (2)

if C(x2,y2) is non-zero, we can divide (1) by (2) and cancel out c.

C(x1,y1)/C(x2,y2) = (B(x1,y1)-A(x1,y1))/(B(x2,y2)-A(x2,y2)) - (3)

Let C(x1,y1)/C(x2,y2) = d. We know d, its value varies depending on the pixel pair we choose. Rearranging (3) we get

d.B(x2,y2) = B(x1,y1)+[d.A(x2,y2) - A(x1,y1)] - (4)

here, we know d, A(x1,y1) and A(x2,y2). Therefore we have derived a linear relationship between two pixel values in image B.

Let's rewrite (4) by taking e=1/d, assuming d is non-zero.

B(x2,y2) = e.B(x1,y1) + [A(x2,y2)-e.A(x1,y1)]
B(x2,y2) = e.B(x1,y1) + f - (5)

where f=[A(x2,y2)-e.A(x1,y1)].

Here, we know both e and f, and their values are dependent on the pixel pair we choose. Using (5) we can establish a recurrence relationship between pixel values of B, like:

B(x1, y1) = e1.B(x0, y0) + f1
B(x2, y2) = e2.B(x1, y1) + f2 = e2.(e1.B(x0, y0) + f1) + f2 = e3.B(x0, y0) + f3
:

where e2.e1 = e3 and e2.f1+f2=f3.

For example,

B(1, 0) = e1.B(0, 0) + f1
B(2, 0) = e2.B(0, 0) + f2
B(3, 0) = e3.B(0, 0) + f3
:
B(i, j) = ei.B(0, 0) + fi

Suppose we find all the relations, that is, these es and fs for all the other pixels of image B in terms of B(0, 0), and set B(0, 0) = 1.

Then,

B(0, 0) = 1
B(1, 0) = e1 + f1
B(2, 0) = e2 + f2
B(3, 0) = e3 + f3
:
B(i, j) = ei + fi

Therefore, we have derived the required image B pixels up to an affine transformation.

This is a very interesting problem. I'll try to write some code to try this out when I get some free time.