Automatic separation of two images that have been multiplied together

I am searching 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 (such as a relatively complicated scene). Image 2 is very simple, and can be mathematically generated. Image 2 always has similar morphology (i.e. downward trend). By multiplying Image 1 by Image 2 (using point-by-point multiplication), we get a transformed image. Given only the transformed image, I would like to estimate Image 1 or Image 2. Is there an algorithm that can do this? Here are 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), [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 ``` UPDATE There are a number of ways to approach this problem. Here are the results of my experiments. Thank you to all who responded to my question! 1. Low-pass filtering of image As noted by duskwuff, taking the low-pass filter of the transformed image returns an approximation of Image 2. In this case, the low-pass filter has been Gaussian. You can see that it is possible to identify multiplicative noise in the image using the low-pass filter. 2. Homomorphic Filtering As suggested by EitenT I examined homomorphic filtering. Knowing the name of this type of image filtering, I managed to find a number of references that I think would be useful in solving similar problems. S. P. Banks, Signal processing, image processing, and pattern recognition. New York: Prentice Hall, 1990. A. Oppenheim, R. Schafer, and J. Stockham, T., “Nonlinear filtering of multiplied and convolved signals,” IEEE Transactions on Audio and Electroacoustics, vol. 16, no. 3, pp. 437 – 466, Sep. 1968. Blind image Deconvolution: theory and applications. Boca Raton: CRC Press, 2007. Chapter 5 of the Blind image deconvolution book is particularly 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 objective function with `fminsearch`. Since I know the mathematical model of the noise, I was able to use this as input to an optimization algorithm. This approach is entirely problem-specific, and may not be always useful in all situations. Here is a reconstruction of Image 2: Here is a reconstruction of Image 1, formed by dividing by the reconstruction of Image 2: Here is the image containing the noise: Source code Here is the source code for my problem. As shown by 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), [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 ```