Учитывая наложенное изображение (A ^ B) и одно исходное изображение (A), как рассчитать другое исходное изображение (B)

Question

Предположим следующее:

• Два изображения были наложены друг на друга и объединены в одно изображение (назовем это: изображение Y)
• У меня есть одно из исходных изображений(давайте назовем это: изображение A). Это изображение является полупрозрачной маской, которая была наложена поверх другого изображения. Но я не знаю математической операции. Это было не просто добавлено, но, вероятно, умножено. Или, может быть, объединены в альфа-композит.

Как получить другое исходное изображение (изображение 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))

• Я хотел бы либо обратить реверсирование умножения (используя что-то вроде деления), либо реверсировать альфа-композитинг. Я предполагаю, что проблема может быть решена. Это должно быть что-то вроде линейной алгебры с 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 здесь .

Я проверил его на 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

Зависит от того, как слились изображения. Если это простое альфа-смешение Y = alpha * A + (1-alpha) * B, то вы можете использовать превосходный ответ Андрея Смородова , если у вас есть более одного наложенного изображения (которого у вас, вероятно, нет).
Не знаюдумаю, что вы можете вычислить 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)

Когда вы устанавливаете альфа-слайдер на 30 (альфа = 0,3), вы 'Увидим оригинальную оправку, для 100 смешанного входного изображения.

---

Но вы упомянули прозрачность: если у вас есть наложение A с прозрачными и непрозрачными областями, и вы комбинируете это наложение с фоном B, чтобы Y = B для прозрачных областей A и Y = A в противном случае (как в этом вопросе SO ), то нет шансов восстановить B для непрозрачных областей A, так как 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, поскольку нам известны и 1013 *, и 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.

. Перепишем (4), взяв e=1/d, предполагая, что d не равно нулю.

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

Предположим, мы находим все отношения, то есть эти e s и f s для всех остальных пикселей изображения B в терминах B(0, 0), и устанавливаем 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 пикселей до аффинного преобразования.

Это очень интересная проблема. Я попытаюсь написать код, чтобы попробовать это, когда у меня будет немного свободного времени.