function [dgr dr dgb db] = show_tca(file)
% experimental tool to display and estimate tca correction coefficients,
% based on a hugin project file that specifies control points between
% r and g, and g and b channel.
%
% It will read the specified pto file, extract the distortion
% correction values and control points. A graph, with the
% scaling factor over the distance from the image center will
% be plotted.
%
% Additionally, a correction polynomial is estimated useing these
% distances. This is better than the optimized correction polynomial,
% because PTOptimizer minimizes the distance between the two points,
% and this script minimizes the difference between the distances
% from the center.
%
% Licensed under the GPL version 2 by Pablo d'Angelo
%
% Changes:
% 2006-11-11  do not use pcregexp anymore, should work with standard octave
%             installations.
%

plot1 = 1;
plot2 = 1;
save_plots=0;

% read control points from file
f = fopen(file,'r');

red = [];
blue = [];

ptodist = [];
ptov = [];

while ~ feof(f)
  line = fgetl(f);
  if (length(line) == 0)
    continue;
  end
  % new code to parse i line without pcregexp, which seems to behave differently
  % with various octave versions
  if (line(1) == 'i')
    str = line(3:end);
    imsize = [0 0];
    distc = [0 0 0];
    v = 0;
    while (length(str) > 0)
      [t str] = strtok(str);
      if t(1) == 'w'; imsize(1) = str2num(t(2:end)); end
      if t(1) == 'h', imsize(2) = str2num(t(2:end)); end
      if t(1) == 'a', distc(1) = str2num(t(2:end)); end
      if t(1) == 'b', distc(2) = str2num(t(2:end)); end
      if t(1) == 'c', distc(3) = str2num(t(2:end)); end
      if t(1) == 'v', v = str2num(t(2:end)); end
    end
    ptodist = [ptodist ; distc ];
    ptov = [ptov ; v];
  end
  % extract image size
  % this code works from some versions of octave, with octave-forge
  %s = pcregexp('^i.* w(\\d+) .*h(\\d+) .* a([-\\w\\.]+) b([-\\w\\.]+) c([-\\w\\.]+) .* v([\\w\\.]+)', line);
  %if size(s,1) == 7
  %  imsize = [str2num(line(s(2,1):s(2,2))), str2num(line(s(3,1):s(3,2)))];
  %  distc = [str2num(line(s(4,1):s(4,2))), str2num(line(s(5,1):s(5,2))), str2num(line(s(6,1):s(6,2))) ]; 
  %  ptodist = [ptodist ; distc];
  %  ptov = [ ptov ; str2num(line(s(7,1):s(7,2))) ];
  %end

  % c n0 N1 x20 y2724 X20.6432 Y2723.36 t0
  if (line(1) == 'c')
    t = sscanf(line, 'c n%d N%d x%f y%f X%f Y%f t%d');
    if t(1) == 1
      a = 1;
      b = 2;
      ai = 3;
      bi = 5;
    elseif t(2) == 1
      a = 2;
      b = 1;
      ai = 5;
      bi = 3;
    else
      continue;
    end
    if t(b) == 0,
      % red channel
      v = [ t(ai:ai+1)' t(bi:bi+1)' ];
      red = [ red ; v];
    end
    if t(b) == 2,
      v = [ t(ai:ai+1)' t(bi:bi+1)' ];
      blue = [ blue ; v];
    end
  end
end

fclose(f);

% calculate distortion parameters estimated by PTOptimizer
scale = ptov(2) ./ ptov;
ptdist = ones(3,1) - sum(ptodist,2);
ptdist = [ptodist ptdist];
ptdist = ptdist .* [scale.^4 scale.^3 scale.^2 scale];

center = imsize / 2;

rfactor = min(imsize) / 2;
maxr = norm(imsize/2);

[dr, dgr] = dist(red, center);
[db, dgb] = dist(blue, center);

disp('correction parameters read from pto file:')
disp(sprintf(' -r %.7f:%.7f:%.7f:%.7f', ptdist(1,:)));
disp(sprintf(' -b %.7f:%.7f:%.7f:%.7f', ptdist(3,:)));


% fit own correction function

% fit polynom to the function here
coefr2 = polyfit(dgr/rfactor, dgr./dr,3);
coefb2 = polyfit(dgb/rfactor, dgb./db,3);

disp('new fit, based on distance to image center:')
disp(sprintf(' -r %.7f:%.7f:%.7f:%.7f', coefr2));
disp(sprintf(' -b %.7f:%.7f:%.7f:%.7f', coefb2));


if plot1,
figure
plot(dgr/rfactor, dgr./dr, 'r+');
hold on;

plot(dgb/rfactor, dgb./db, 'b+');
xlabel('distance from center');
ylabel('scale factor');

% plot correction read from pto file
x = (1:50:maxr)/rfactor;
plot(x, polyval(ptdist(1,:), x), 'r-');
plot(x, polyval(ptdist(2,:), x), 'g-');
plot(x, polyval(ptdist(3,:), x), 'b-');


plot(x, polyval(coefr2,x), 'r-@*');
plot(x, polyval(coefb2,x), 'b-@*');

%legend('red channel', 'blue channel', 'PT red', 'PT green', 'PT blue', 'new red', 'new blue', 3);

if save_plots
	print('corr_curves.eps','-color', '-F:18');
end

end

% calculate distance perpendicular to ray from image center
distr = red(:,1:2) - red(:,3:4);
distr = sqrt(sum(distr.^2,2));
dr_perp = sqrt(distr.^2 - (abs(dr-dgr)).^2);

if plot2,
figure
plot(dgr, dr_perp,'r*');
hold on;
plot(dgr, abs(dr-dgr),'bo');
xlabel('distance from center (pixels)');
ylabel('distance between green and red channel (pixels)');
legend('red: tangential distance', 'blue: sagittal distance');
if save_plots
  print('distance.eps','-color','-F:18');
end
end

function [dc, dg] = dist(a, center)
% calculate distance between the control points..
dchannel = a(:,1:2) - repmat(center, size(a,1),1);
dgreen =  a(:,3:4) - repmat(center, size(a,1),1);

% calculate difference of distance from center
dg =  sqrt(sum(dgreen.^2,2));
dc = sqrt(sum(dchannel.^2,2));
distance = dg - dc;