function doHybrid(infile, outfile)

% Usage: doHybrid(infile, outfile)
%
% A interesting algorithm which uses a linear combination
% of the scale factors obtained using the gray world
% assumption and those obtained using the perfect reflector
% assumption. The weights depend on R+G+B.
%
% KLY 98

load /local/class/psych221/tutorials/cMatch/monitorGam;
load /local/class/psych221/tutorials/cMatch/monitor;
load xyz
rgb2xyz = xyz'*phosphors;

[r g b] = tiffread(infile);

imgSize = size(r);

r = reshape(r, 1, imgSize(1)*imgSize(2))';
g = reshape(g, 1, imgSize(1)*imgSize(2))';
b = reshape(b, 1, imgSize(1)*imgSize(2))';

rgb = [r,g,b];
maxQ = size(monitorGam,1);
rgb = scale(rgb,1,maxQ);
linrgb = monitorGam(round(rgb));

imagexyz = rgb2xyz*linrgb';

% Locate specularities by looking for the pixel with the highest
% Y value.

maxY = max(imagexyz(2,:));
index = find(imagexyz(2,:)==maxY);
index = index(1);

% Display white point

disp('Located white point...')
whitePt = linrgb(index,:)
PRscaleFactors = 1./whitePt;

% Calculate averages (gray world assumption)

disp('Image averages...')
averagergb = sum(linrgb)
averageGray = sum(averagergb)./3
GWscaleFactors = averageGray./averagergb;

for i=1:imgSize(1)*imgSize(2)

  pixelsum = sum(linrgb(i,:));
    
  % Consider three regions
    
  if pixelsum > 2.5
    linrgb(i,:) = linrgb(i,:).*PRscaleFactors;

  elseif pixelsum < .5
    linrgb(i,:) = linrgb(i,:).*GWscaleFactors;

  else
    fractionS = (pixelsum - .5)/2;
    fractionG = 1 - fractionS;
    linrgb(i,:) = linrgb(i,:).*GWscaleFactors.*fractionG + ...
	linrgb(i,:).*PRscaleFactors.*fractionS;

  end
end

maxQ = size(monitorInvGam,1);
rgb = scale(linrgb,1,maxQ);
rgb = monitorInvGam(round(rgb))./255;

r = reshape(rgb(:,1), imgSize(1), imgSize(2));
g = reshape(rgb(:,2), imgSize(1), imgSize(2));
b = reshape(rgb(:,3), imgSize(1), imgSize(2));

tiffwrite(r, g, b, outfile);







