% comparahist: compares histograms; S. Mann
% see also Contrast Limited Adaptive Histogram Specification (local stats...?)

frame1=13
frame2=16
%frame2=34 % frame 13 and 34 are the same exposure, straight line where certain
%frame2=37;

eval(sprintf('image1=loadpnm(''s%03d.pgm'');',frame1));
eval(sprintf('image2=loadpnm(''s%03d.pgm'');',frame2));

[M,N]=size(image1);
%image1=image1(:,1:N-25);
%image2=image2(:,1+25:N);
%image1=image1(:,1:N-50);
%image2=image2(:,1+50:N);

J=comparagram(image1,image2);
%Jl=J;
%Jl=sqrt(J); % sqrt() or log(()+1) usually improves results
myeps=.001;
Jl=log(J+myeps) - log(myeps); % make it be nonnegative
%h1=hist256(image1);
%h2=hist256(image2);
h1=sum(Jl.');
h2=sum(Jl);
x=0:255;
axisscale=max(max(h1),max(h2));
subplot(221); plot(x,h1,x,h2); axis([0,255,0,axisscale]); title('HISTOGRAMS');
ylabel('h_1 and h_2');

H1=cumsum(h1);
H2=cumsum(h2);
axisscale=max(max(H1),max(H2)); % max(H1) should equal max(H2)
p=1.05; % looks nicer if we can see clearly where they meet at the upper right
subplot(222); plot(x,H1,x,H2); axis([0,255*p,0,axisscale*p]);
title('CUMULATIVES')
ylabel('H_1 and H_2');

g=NaN*ones(length(H1),1);
for gind=1:256
 gh=min(find(H1(gind)<=H2)); % max(find(x<4)) % gh=8
 gl=max(find(H1(gind)>=H2));                  % gl=7
 %H2(gh) ans = 12655
 %H2(gl) ans = 8055
 %H1(gind) ans = 9504
 if gl>gh % in certain rare cases, g low is greater than g high
   ghigh=gl;
   gl=gh;
   gh=ghigh;
 end%if
 if gh==gl
   g(gind) = gl; % both the same
 elseif isempty(gl)
   g(gind) = gh;
 elseif isempty(gh)
   g(gind) = gl;
 else%if
   if H2(gh)==H2(gl) % interpolating between identical points; zero denominator
    g(gind) = gl+H2(gl); % doesn't matter which
   else
    g(gind)= gl+(H2(gh)-H1(gind))/(H2(gh)-H2(gl)); % linear interp
   end%if
 end%if
 % (12655-9504)/(12655-8055) = 0.68500
end%for
g = g - 1; % indices were 1:256, but g is a real value on iterval [0..255]
g = monotonize(g); % force it to be monotonic; that's the best we can do
                   % because forcing derivative to be monotonic is unrealistic
subplot(223); plot(g); axis([0,255,0,255]); axis('square'); title('CUMULAGRAM')
xlabel('f'); ylabel('g(f)=H_2^{-1}(H_1(f))');

%J=comparagram(image1,image2);
subplot(224);
%tvs(-Jl.');
axis('xy');
axis('square');
title('COMPARAGRAM'); xlabel('f'); ylabel('g(f)');
hold on
plot(g); axis([0,255,0,255]); axis('square');
hold off