Why surface fitting is not working in logarithmic domain?

Algorithm:

- Original Image without any bias field is introduced with the polynomial profile, which results in an image having non-uniform illumination.

- Biased image is converted in log domain and surface is approximated using polynomial fit

- approximated surface is subtracted from the Biased image which results in original image back with no bias fields.(Ideally).

- measure RMSE between approximated surface and introduced polynomial field in step 1. Measure RMSE between the Biased Image and the image we get back at the end after subtraction.

Code:

clear;clc;close all;

%read the image, on which profile is to be generated
I = ones(300);
I = padarray(I,[20,20],'both','symmetric'); % padding
%creating a bias profile using polynomial modeling
[x,y] = meshgrid(1:size(I,1),1:size(I,2));
profile = -2.5.*x.^3 - 2.5.* y.^3 + 0.25 .*(x.* y.^2) - 0.3*( x.^2 .* y ) - 0.5.* x .* y - x + y  - 2.5*( x.^2) - y.^2 + 0.5 .* x .*y + 1;

% come to scale [0 1] 
profile = profile - min(profile(:));
profile = profile / max(profile(:));

figure,imshow(profile,[]); %introduced bias profile

%%corrupt the image 
biasedImage = (I .* profile);

figure,imshow(biasedImage,[]); %biased Image

cImage = log(biasedImage+1);% conversion to log domain

%%forming the input for prediction of surface
colorChannel = cImage;
rows = size(colorChannel, 1);
columns = size(colorChannel, 2);
[X,Y] = meshgrid(1:columns, 1:rows);
z = colorChannel;
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
z1d = double(reshape(z, numel(z), 1)); %dependent variables
x = [x1d y1d]; % two regressors

%%surface fitting
m = 3; %define the order of polynomial
p = polyfitn(x,z1d,m); %prediction step

zg = polyvaln(p, x);

modeledColorChannel = reshape(zg, [rows columns]); % predicted surface
%modeledColorChannel = exp(modeledColorChannel)-1; if we use this step then the step below will be division instead of subtraction

%%correction

f  = cImage - modeledColorChannel; %getting the original image back. 

%grayImage = f(21:end-20,21:end-20);
%modeledColorChannel = modeledColorChannel(21:end-20,21:end-20); %to remove the padding

figure,imshow(f,[]);
figure,imshow(modeledColorChannel,[]);

%%measure the RMSE for image
y = (I - f);
RMSE = sqrt(mean(y(:).^2));
disp(RMSE);

% RMSE for profile
z = (modeledColorChannel - profile);
RMSE = sqrt(mean(z(:).^2));
disp(RMSE);

Now, the question is: I am getting lower RMSE value at the end if I do division in the log domain. How does that possible? If I keep all of my calculation in log domain image division will become image subtraction, which ultimately eliminates the need to apply exp(correctedImage/f here).It is obvious if you run the code and see the image f at the end of the correction.

In both the cases, RMSE between the introduced and perceived profile is same as I am doing my estimation in log domain in both the cases.Either image division or in image subtraction. Only difference is in image division we will have to come back to the original image domain from log domain.

Use polyfitn tool box for estimation.

Thanks.