I am looking for an algorithm or C ++ / Matlab library that can be used to separate two images multiplied together. A visual example of this problem is given below.
Image 1 can be anything (for example, a relatively complex scene). Image 2 is very simple and can be mathematically generated. Image 2 always has a similar morphology (i.e., a downward trend). Multiplying image 1 by image 2 (using point multiplication), we get the converted image.
Given only the converted image, I would like to evaluate image 1 or image 2. Is there an algorithm that can do this?
Here is the Matlab code and images:
load('trans.mat'); imageA = imread('room.jpg'); imageB = abs(response); % loaded from MAT file [m,n] = size(imageA); image1 = rgb2gray( imresize(im2double(imageA), [mn]) ); image2 = imresize(im2double(imageB), [mn]); 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
UPDATE
There are several ways to solve this problem. Here are the results of my experiments. Thanks to everyone who answered my question!
1. Low-pass image filtering
As duskwuff noted, receiving a low-pass filter of the converted image returns an approximation of image 2. In this case, the low-pass filter was Gaussian. You can see that it is possible to identify the multiplicative noise in the image using a low-pass filter.
2. Homomorphic filtration
As suggested by EitenT, I looked at homomorphic filtering. Knowing the name of this type of image filtering, I managed to find a number of links that, I think, would be useful in solving such problems.
- S
. P. Banks, signal processing, image processing and pattern recognition. New York: Prentice Hall, 1990.
but. Oppenheim, R. Schafer, J. Stockham, T., Nonlinear Filtering of Multiplied and Coiled Signals, IEEE Transactions on Audio and Electroacoustics, vol. 16, no. 3, p. 437-446, September 1968
Blind image Deconvolution: theory and applications. Boca Raton: CRC Press, 2007.
Chapter 5 of the book of deconvolution of blind images is especially good and contains many references to homomorphic filtering. This is perhaps the most generalized approach that will work well in many different applications.
3. Optimization using fminsearch
As suggested by Serg, I used an object function with fminsearch . Since I know the mathematical model of noise, I was able to use this as an input to an optimization algorithm. This approach is completely specific for a specific task and may not always be useful in all situations.
The following is a reconstruction of image 2:
Here is a reconstruction of image 1 formed by dividing by image restoration 2:
Here is an image containing noise:
Source
Here is the source code for my problem. As shown in the code, this is a very specific application and will not work well in all situations.
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), [mn]) ); % 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