Embedding a shadow 1d invariant image

I implemented a shadow removal method based on the invariant color functions found in Minimizing Entropy to Remove Shadows. My implementation seems to sometimes produce similar computation results, but they are always off and the grayscale image is blocky, possibly due to a miscalculation of the geometric mean.

Here is an example of the dependence of the information potential on the image of the horse in the document, as well as my invariant image. Multiply the x-axis by 3 to get theta (which goes from 0 to 180):

And here is the grayscale image. My code outputs for correct maximum theta (mine is off by 10):

You can see that their image is not blocky:

Here is their informational potential:

When dividing by the geometric mean, I tried to use NaN and transform the image, so the smallest possible value is .01, but it doesn't seem to change my result.

Here is my code:

I = im2double(imread(strname));
[m,n,d] = size(I);
I = max(I, .01);
chrom = zeros(m, n, 3, 'double');
for i = 1:m
    for j = 1:n
       % if ((I(i,j,1)*I(i,j,2)*I(i,j,3))~= 0)
              chrom(i,j, 1) = I(i,j,1)/((I(i,j,1)*I(i,j,2)*I(i,j, 3))^(1/3));
              chrom(i,j, 2) = I(i,j,2)/((I(i,j,1)*I(i,j,2)*I(i,j, 3))^(1/3));
              chrom(i,j, 3) = I(i,j,3)/((I(i,j,1)*I(i,j,2)*I(i,j, 3))^(1/3));
          % else
          %    chrom(i,j, 1) = 1;
          %    chrom(i,j, 2) = 1;
          %    chrom(i,j, 3) = 1;
        % end
    end
end

p1 = mat2gray(log(chrom(:,:,1)));
p2 = mat2gray(log(chrom(:,:,2)));
p3 = mat2gray(log(chrom(:,:,3)));
X1 = mat2gray(p1*1/(sqrt(2)) - p2*1/(sqrt(2)));
X2 = mat2gray(p1*1/(sqrt(6)) + p2*1/(sqrt(6)) - p3*2/(sqrt(6)));
maxinf = 0;
maxtheta = 0;
data2 = zeros(1, 61);
for theta = 0:3:180
M = X1*cos(theta*pi/180) - X2*sin(theta*pi/180);
s = sqrt(std2(X1)^(2)*cos(theta*pi/180) + std2(X2)^(2)*sin(theta*pi/180));
s = abs(1.06*s*((m*n)^(-1/5)));
[m, n] = size(M);
length = m*n;
sources = zeros(1, length, 'double');
count = 1;
for x=1:m
    for y = 1:n
    sources(1, count) = M(x , y);
    count = count + 1;
end
end
weights = ones(1, length);
sigma = 2*s;
[xc , Ak] = fgt_model(sources , weights , sigma , 10, sqrt(length) , 6  );
sum1 = sum(fgt_predict(sources , xc , Ak , sigma , 10 ));
sum1 = sum1/sqrt(2*pi*2*s*s);
data2(theta/3 + 1) = sum1;
if (sum1 > maxinf)
   maxinf = sum1;
   maxtheta = theta;
end
end
InvariantImage2 = cos(maxtheta*pi/180)*X1 + sin(maxtheta*pi/180)*X2;

      

        
        
        
      

    

Let's assume the Fast Gauss transform is correct.