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

frame1=13
frame2=16
%frame2=36 %13 and 35 with Jl=J; try: gbad=g;gbad(123:178)=NaN;plot(gbad)
%frame2=34 % frame 13 and 34 are the same exposure, straight line where certain
%frame2=37;
%frame1=46;
%frame2=55;

eval(sprintf('image1=loadpnm(''s%03d.pgm'');',frame1));
eval(sprintf('image2=loadpnm(''s%03d.pgm'');',frame2));
%eval(sprintf('image1=loadpnm(''pork%02d.pgm'');',frame1));
%eval(sprintf('image2=loadpnm(''pork%02d.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.^1.01;
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_f and h_g');

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('CUMULAGRAM')
ylabel('H_f and H_g');

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 isempty(gl)
   g(gind) = gh;
 elseif isempty(gh)
   g(gind) = gl;
 elseif gh==gl
   g(gind) = gl; % both the same
 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('COMPARAGRAPH')
xlabel('f'); ylabel('g(f)=monotonize(H_g^{-1}(H_f(f)))');

%J=comparagram(image1,image2);
%subplot(224); tvs(-sqrt(sqrt(J)).'); axis('xy'); axis('square');
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

print('to use this figure, save as eps, then run pstoedit, for idraw')
print('print -deps pork.eps')
print('pstoedit -f idraw pork.eps pork.idraw')

disp('now assuming youre doing s13 versus s35, some data needs interp')
% use the interp1, as in this example:
%>> x = 0:10; y = sin(x); xi = 0:.01:10;
%>> yi = interp1(x,y,xi,'cubic'); plot(x,y,'o',xi,yi)

gbad=[g(1:122);g(179:256)];
x=[1:122 179:256];
x=x-1;
xi=0:255;

%gi=interp1(x,gbad,xi,'cubic'); plot(xi,gi) %%% plot(x,gbad,'o',xi,gi)
% 'cubic' and 'spline' fail with noisy data
gi=interp1(x,gbad,xi,'linear'); plot(xi,gi) %%% plot(x,gbad,'o',xi,gi)