给定一张叠加图（A^B）和一张原图（A），如何计算另一张原图（B）

Question

假设如下：

• 两张图片相互叠加并合并为一张图片（我们称之为：图片 Y）
• 我有一张原始图像（我们称之为：图像 A）。此图像是覆盖在另一张图像上的半透明蒙版。但我不知道数学运算。它不仅被添加，而且可能成倍增加。或者可以组合成一个 alpha 复合材料。

在给定其他两个图像的情况下，如何获得另一个原始图像（图像 B）？所有图像均为 RGBA 图像。

我曾尝试使用 PIL 的 PIL.ImageChops.difference() 但这并不能解决问题。数学运算不仅仅是简单的加法。

我还尝试使用我找到的一些代码。但它肯定不会产生正确的结果。

# 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))

• 我想反转乘法（通过使用除法之类的方法）或反转 alpha 合成。我认为问题可以解决。它一定是线性代数，f(A,B) = Y。我有 A 和 Y，但是 B 是什么？
• 我的问题的正确名称是什么？

Answer 1

如果您有多个线性组合的分量，如果您有相同数量或更多的混合（以不同的系数混合），则可以恢复这些分量。

例如：

你有：

Mix1=I1*c1+I2*c2 和 Mix2=I1*c3+I2*c4，其中 Mix1 和 Mix2 是您拥有的混合物。 I1和I2是源图像，c1-c4是混合系数。

使用盲源分离可以恢复 I1 和 I2。

要快速开始使用 pyton 版本，请查看 here。

我对它进行了 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();

}

输入 2 个混合物：

因此我们有 2 张分离的图像。

Answer 2

这取决于图像的合并方式。如果它是一个简单的 alpha 混合 Y = alpha * A + (1-alpha) * B，那么您可以使用 Andrey Smorodov 的出色 answer，因为您有多个叠加图像（您可能没有）。
我不认为你可以从上面的等式计算B，因为你有一个带有两个未知数的等式（B 和alpha）。您可以做的是将上述两个图像反混合并根据视觉印象任意调整alpha：

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)

当您将 alpha 滑块调整为 30 (alpha = 0.3) 时，您将看到原始山魈，100 为混合输入图像。

---

但是您提到了透明度：如果您有一个具有透明和不透明区域的叠加层 A，并将此叠加层与背景 B 组合，使得 A 的透明区域的 Y = B 和 Y = A 否则（如this SO question），那么就有没有机会为 A 的不透明区域重建 B，因为 B 在此处 复制 到 Y 中，没有 A 的份额 (Y = 1.0 * B + 0.0 * A)。

Answer 3

如果图像是线性混合的，并且如果我们考虑一个理想的情况，即由于像素可以表示的值范围有限（比如单通道图像中的 8 位），我们不会丢失精度，那么，我认为我们可以将问题表述如下（我认为这应该可行，至少在理论上，但如果有任何缺陷，请指出）：

如果

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 可供我们使用。

然后

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

现在，我们知道Y-A，因为我们都知道Y 和A。让Y-A=C。那么

C = c.(B-A)

现在让我们考虑(x, y)位置的图像像素

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

让我们考虑位于(x1,y1) 和(x2,y2) 位置的两个图像像素

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

如果 C(x2,y2) 不为零，我们可以将 (1) 除以 (2) 并抵消 c。

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

让C(x1,y1)/C(x2,y2) = d。我们知道d，它的值取决于我们选择的像素对。 重新排列 (3) 我们得到

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

在这里，我们知道d、A(x1,y1) 和A(x2,y2)。因此我们推导出了图像B中两个像素值之间的线性关系。

假设 d 不为零，让我们用 e=1/d 重写 (4)。

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

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

在这里，我们知道e 和f，它们的值取决于我们选择的像素对。 使用(5)我们可以建立B的像素值之间的递归关系，如：

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
:

e2.e1 = e3 和 e2.f1+f2=f3 的位置。

例如，

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

假设我们根据B(0, 0)找到图像B的所有其他像素的所有关系，即这些es和fs，并设置B(0, 0) = 1。

那么，

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

因此，我们将所需的图像B 像素推导出来进行仿射变换。

这是一个非常有趣的问题。有空我会试着写一些代码来试试这个。