function grayWorld(filename,outFile,catType,maxIter,plots)
%grayWorld(filename,outfile,catType,maxIter,plot)
% Performs color balancing via the gray world assumption and a chromatic
% adpatation transform.  Can be run interatively to possibly improve the
% result.  Set plot = 0 or 1 to turn diagnostic plots on or off.

tutorialinit

%% import image
im_orig = cbimread(filename);
% figure
% imshow(im_orig)
% title('Original Image')

%% various constants
xyz_D65 = [95.04; 100; 108.88]; %http://en.wikipedia.org/wiki/D65, normalized Y = 100

sRGBtoXYZ =   [0.4124564  0.3575761  0.1804375; ...
               0.2126729  0.7151522  0.0721750; ...
               0.0193339  0.1191920  0.9503041];

%% grayworld
% catType = 'vonKries';
% maxIter = 1;
b = 0.001; %convergence threshold

imRGB_orig = cbreshape(im_orig)*255;
imRGB = imRGB_orig;

grayDiff =[];
for iter = 1:maxIter
    rgbEst = mean(imRGB,2); %grayworld, average everything
    grayDiff = [grayDiff norm([rgbEst(1)-rgbEst(2),rgbEst(1)-rgbEst(3),rgbEst(2)-rgbEst(3)])];
    
    if grayDiff(end) < b
        disp(['Converged. RGB difference vector < ' num2str(b) ' in magnitude.'])
        break
    elseif iter >= 2 && abs(grayDiff(end-1)-grayDiff(end)) < 10^-6
        disp(['RGB difference vector no longer improving.'])
        break
    end
    
    xyEst = XYZ2xy(sRGBtoXYZ*rgbEst); %calculate xy chromaticity
    xyzEst = xy2XYZ(xyEst,100); %normalize Y to 100 so D65 luminance comparable
    imRGB = cbCAT(xyzEst,xyz_D65,catType)*imRGB;
end

imwrite(cbunshape(imRGB,size(im_orig))/255,outFile,'png');

if plots
    length(grayDiff) %number of iterations done
%     figure
%     imshow(cbunshape(imRGB,size(im_orig))/255)
%     title('Gray World Corrected')
    
    figure
    plot(grayDiff)
    title('GW: Norm of RGB Difference Vector vs Iterations')
    print(gcf,'-dpng',[outFile '-fig1'])
end
