乗算された2つの画像を分離するために使用できるアルゴリズムまたはC++/Matlabライブラリを探しています。この問題の視覚的な例を以下に示します。
画像1は何でもかまいません（比較的複雑なシーンなど）。画像2は非常に単純で、数学的に生成できます。画像2は常に同様の形態を持っています（つまり、下降傾向）。画像1に画像2を乗算することにより（ポイントごとの乗算を使用）、変換された画像を取得します。
変換された画像のみを前提として、画像1または画像2を推定したいと思います。これを実行できるアルゴリズムはありますか？
Matlabのコードと画像は次のとおりです。
load('trans.mat');
imageA = imread('room.jpg');
imageB = abs(response); % loaded from MAT file
[m,n] = size(imageA);
image1 = rgb2gray( imresize(im2double(imageA), [m n]) );
image2 = imresize(im2double(imageB), [m n]);
figure; imagesc(image1); colormap gray; title('Image 1 of Room')
colorbar
figure; imagesc(image2); colormap gray; title('Image 2 of Response')
colorbar
% This is image1 and image2 multiplied together (point-by-point)
trans = image1 .* image2;
figure; imagesc(trans); colormap gray; title('Transformed Image')
colorbar
アップデート
この問題に取り組む方法はいくつかあります。これが私の実験の結果です。私の質問に答えてくれたすべての人に感謝します！
1.画像のローパスフィルタリング
duskwuffが指摘しているように、変換された画像のローパスフィルターを使用すると、画像2の近似値が返されます。この場合、ローパスフィルターはガウス分布になっています。ローパスフィルターを使用すると、画像内の乗法性ノイズを識別できることがわかります。
2.準同型フィルタリング
EitenTによって提案されたように、私は準同型フィルタリングを調べました。このタイプの画像フィルタリングの名前を知っていたので、同様の問題を解決するのに役立つと思う多くの参照を見つけることができました。
SPバンク、信号処理、画像処理、およびパターン認識。ニューヨーク：プレンティスホール、1990年。
A. Oppenheim、R。Schafer、およびJ. Stockham、T。、「乗算および畳み込み信号の非線形フィルタリング」、IEEE Transactions on Audio and Electroacoustics、vol。16、いいえ。3、pp。437 – 466、1968年9月。
ブラインドイメージデコンボリューション：理論と応用。ボカラトン：CRC Press、2007年。
ブラインドイメージデコンボリューションブックの第5章は特に優れており、準同型フィルタリングへの多くの参照が含まれています。これはおそらく、多くの異なるアプリケーションでうまく機能する最も一般化されたアプローチです。
3.を使用した最適化fminsearch
Sergが提案したように、私はで目的関数を使用しましたfminsearch。ノイズの数学的モデルを知っているので、これを最適化アルゴリズムへの入力として使用することができました。 このアプローチは完全に問題固有であり、すべての状況で常に役立つとは限りません。
これが画像2の再構成です。
これは、画像2の再構成で除算することによって形成された画像1の再構成です。
ノイズを含む画像は次のとおりです。
ソースコード
これが私の問題のソースコードです。コードで示されているように、これは非常に特殊なアプリケーションであり、すべての状況でうまく機能するとは限りません。
N = 1001;
q = zeros(N, 1);
q(1:200) = 55;
q(201:300) = 120;
q(301:400) = 70;
q(401:600) = 40;
q(601:800) = 100;
q(801:1001) = 70;
dt = 0.0042;
fs = 1 / dt;
wSize = 101;
Glim = 20;
ginv = 0;
[R, ~, ~] = get_response(N, q, dt, wSize, Glim, ginv);
rows = wSize;
cols = N;
cut_val = 200;
figure; imagesc(abs(R)); title('Matrix output of algorithm')
colorbar
figure;
imagesc(abs(R)); title('abs(response)')
figure;
imagesc(imag(R)); title('imag(response)')
imageA = imread('room.jpg');
% images should be of the same size
[m,n] = size(R);
image1 = rgb2gray( imresize(im2double(imageA), [m n]) );
% here is the multiplication (with the image in complex space)
trans = ((image1.*1i)) .* (R(end:-1:1, :));
figure;
imagesc(abs(trans)); colormap(gray);
% take the imaginary part of the response
imagLogR = imag(log(trans));
% The beginning and end points are not usable
Mderiv = zeros(rows, cols-2);
for k = 1:rows
val = deriv_3pt(imagLogR(k,:), dt);
val(val > cut_val) = 0;
Mderiv(k,:) = val(1:end-1);
end
% This is the derivative of the imaginary part of R
% d/dtau(imag((log(R)))
% Do we need to remove spurious values from the matrix?
figure;
imagesc(abs(log(Mderiv)));
disp('Running iteration');
% Apply curve-fitting to get back the values
% by cycling over the cols
q0 = 10;
q1 = 500;
NN = cols - 2;
qout = zeros(NN, 1);
for k = 1:NN
data = Mderiv(:,k);
qout(k) = fminbnd(@(q) curve_fit_to_get_q(q, dt, rows, data),q0,q1);
end
figure; plot(q); title('q value input as vector');
ylim([0 200]); xlim([0 1001])
figure;
plot(qout); title('Reconstructed q')
ylim([0 200]); xlim([0 1001])
% make the vector the same size as the other
qout2 = [qout(1); qout; qout(end)];
% get the reconstructed response
[RR, ~, ~] = get_response(N, qout2, dt, wSize, Glim, ginv);
RR = RR(end:-1:1,:);
figure; imagesc(abs(RR)); colormap gray
title('Reconstructed Image 2')
colorbar;
% here is the reconstructed image of the room
% NOTE the division in the imagesc function
check0 = image1 .* abs(R(end:-1:1, :));
figure; imagesc(check0./abs(RR)); colormap gray
title('Reconstructed Image 1')
colorbar;
figure; imagesc(check0); colormap gray
title('Original image with noise pattern')
colorbar;
function [response, L, inte] = get_response(N, Q, dt, wSize, Glim, ginv)
fs = 1 / dt;
Npad = wSize - 1;
N1 = wSize + Npad;
N2 = floor(N1 / 2 + 1);
f = (fs/2)*linspace(0,1,N2);
omega = 2 * pi .* f';
omegah = 2 * pi * f(end);
sigma2 = exp(-(0.23*Glim + 1.63));
sign = 1;
if(ginv == 1)
sign = -1;
end
ratio = omega ./ omegah;
rs_r = zeros(N2, 1);
rs_i = zeros(N2, 1);
termr = zeros(N2, 1);
termi = zeros(N2, 1);
termr_sub1 = zeros(N2, 1);
termi_sub1 = zeros(N2, 1);
response = zeros(N2, N);
L = zeros(N2, N);
inte = zeros(N2, N);
% cycle over cols of matrix
for ti = 1:N
term0 = omega ./ (2 .* Q(ti));
gamma = 1 / (pi * Q(ti));
% calculate for the real part
if(ti == 1)
Lambda = ones(N2, 1);
termr_sub1(1) = 0;
termr_sub1(2:end) = term0(2:end) .* (ratio(2:end).^-gamma);
else
termr(1) = 0;
termr(2:end) = term0(2:end) .* (ratio(2:end).^-gamma);
rs_r = rs_r - dt.*(termr + termr_sub1);
termr_sub1 = termr;
Beta = exp( -1 .* -0.5 .* rs_r );
Lambda = (Beta + sigma2) ./ (Beta.^2 + sigma2); % vector
end
% calculate for the complex part
if(ginv == 1)
termi(1) = 0;
termi(2:end) = (ratio(2:end).^(sign .* gamma) - 1) .* omega(2:end);
else
termi = (ratio.^(sign .* gamma) - 1) .* omega;
end
rs_i = rs_i - dt.*(termi + termi_sub1);
termi_sub1 = termi;
integrand = exp( 1i .* -0.5 .* rs_i );
L(:,ti) = Lambda;
inte(:,ti) = integrand;
if(ginv == 1)
response(:,ti) = Lambda .* integrand;
else
response(:,ti) = (1 ./ Lambda) .* integrand;
end
end % ti loop
function sse = curve_fit_to_get_q(q, dt, rows, data)
% q = trial q value
% dt = timestep
% rows = number of rows
% data = actual dataset
fs = 1 / dt;
N2 = rows;
f = (fs/2)*linspace(0,1,N2); % vector for frequency along cols
omega = 2 * pi .* f';
omegah = 2 * pi * f(end);
ratio = omega ./ omegah;
gamma = 1 / (pi * q);
% calculate for the complex part
termi = ((ratio.^(gamma)) - 1) .* omega;
% for now, just reverse termi
termi = termi(end:-1:1);
%
% Do non-linear curve-fitting
% termi is a column-vector with the generated noise pattern
% data is the log-transformed image
% sse is the value that is returned to fminsearchbnd
Error_Vector = termi - data;
sse = sum(Error_Vector.^2);
function output = deriv_3pt(x, dt)
N = length(x);
N0 = N - 1;
output = zeros(N0, 1);
denom = 2 * dt;
for k = 2:N0
output(k - 1) = (x(k+1) - x(k-1)) / denom;
end