Miyazaki's Specular-free method MATLAB implementation

miyazaki.m

%%=== Miyazaki Specular-free MATLAB Implementation ==%%

%%
% Based on Daisuke Miyazaki et al. publication : "Polarization-based Inverse Rendering from a Single View ",
% and awesome explanation at http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche40.html.
% Sorry I don't know the author who post the explanation on that page.
% Illumination Compenstion doesn't implemented yet on this script
%%

%% Miyazaki Function %%
% @params im => image in floating number
% @params arb => arbitrary degree, if arbitrary value is large, image will saturate.
function imnew = miyazaki(im, arb)

%%=== Transform RGM space to M Space ===%

 [h, w, ~] = size(im); % get image size

 % Extract RGB Value %
 imr = im(:,:,1); % Get Red Value
 img = im(:,:,2); % Get Green Value
 imb = im(:,:,3); % Get Blue Value

 % Create Zeros for Store Image in M Space %
 m1 = zeros(h,w);
 m2 = zeros(h,w);
 m3 = zeros(h,w);

 % RGB to M Space Transformation Value %
 a = [1 -1/2 -1/2; 0 sqrt(3)/2 -sqrt(3)/2; 1/3 1/3 1/3];

 % Do RGB to M Space Transformation %
 for i=1:h
  for j=1:w
   rgbmat = [imr(i,j);img(i,j);imb(i,j)];
   imr(i,j) = a(1,:)*rgbmat;
   img(i,j) = a(2,:)*rgbmat;
   imb(i,j) = a(3,:)*rgbmat;
  end
 end

 % Store RGM Transformation to m1, m2, m3 %
 m1 = m1 + imr;
 m2 = m2 + img;
 m3 = m3 + imb;

%%=== Do Miyazaki's Specular-free method ==%%
 m1a = m1;
 m2a = m2;
 m1ak = m1a.^2;
 m2ak = m2a.^2;

%% square root m1ak + m2ak %%
 sm12ak = sqrt(m1ak + m2ak );

 arbitrary = arb;

 m3a = arbitrary*sm12ak;

%%=== Inverse from M to RGB Space ===%%

 % Invers Transformation Value %
 specular =  [2/3 0 1; -1/3 1/sqrt(3) 1; -1/3 -1/sqrt(3) 1];

 r = zeros(h,w);
 g = zeros(h,w);
 b = zeros(h,w);

 % Do Transformation from M to RGB Space %
 for i=1:h
  for j=1:w
   mmat = [m1a(i,j); m2a(i,j); m3a(i,j)];
   ra(i,j) = specular(1,:)*mmat;
   ga(i,j) = specular(2,:)*mmat;
   ba(i,j) = specular(3,:)*mmat;
  end
 end

 % Store M to RGB Transformation space to r, g, b %
 r = r+ra;
 g = g+ga;
 b = b+ba;

 % Rebuild RGB Image as output %
 imnew = cat(3, r, g, b);

end

I'm really new at image processing, please correct me if what I did is wrong