% est_fftpchirp2m.m  ...                                  Steve Mann, 1992,1993
%
% uses FFT abelian guess if no starting paramter, p0, given
%
% multi-resolution (and recursive: calls est_pchirp2r) 2-D image group est.
%
% estimates parameters of the plane projection group
%
% Use: p = est_pchirp2m(E1,E2);                      % units of pchirp2
% or : est_pchirp2m(E1,E2);                          % nicely formatted screen
%                                                    % output
%
%      est_pchirp2m(E1,E2,3);                        % does 3 itterations at
%                                                    % each level; default is
%                                                    % 2 itteration at each
%                                                    % level (min. number
%                                                    % needed to see result)
%
%      est_pchirp2m(E1,E2,3,4,p0);                   % 4 levels and starting p0
%      compared to est_pchirp2r: NOTE NEW PARAMETER IN MIDDLE OF PARAMETER LIST
%      not at end, because I still want to be able to leave off p0 if unknown
%      e.g. est_pchirp2m(E1,E2,2,4); does 2 itt/level, 4 levels, unknown p0.
%      est_pchirp2m(E1,E2,3,4,p0,numleveltoskip);    % skip some top levels
%      est_pchirp2m(E1,E2,3,4,p0,1);                 % skip top level (default)
%      est_pchirp2m(E1,E2,3,4,p0,0);                 % all levels (as before)

function argument_p1 = multi_parameter_function(E1,E2,K,num_levels,p0,numlevelstoskip);

more off  % in case it was on, so we don't stop the process for want of screen

if nargin==1
  error('est_pchirp2m: you must supply 2 input images... exiting...')
end%if

[M  N ] = size(E1);
[M2 N2] = size(E2);
if (M~=M2) ! (N~=N2)
  error('est_pchirp2m: The two input images are not same size... exiting...')
end%if

if nargin <3
  disp1='est_pchirp2m: number of iterations per level not specified; ';
  disp2='defaulting to 2';
  disp([disp1 disp2])
  K=2;
end%if

if nargin < 4
  num_levels=floor(min(log2(M/16),log2(N/16))); % at least 16 pixel image size
  disp(sprintf('est_pchipr2m: defaulting to %g levels',num_levels))
end%if

boxsizes=[1 3 5 9 17 33 65 129 257 513 1025 2049 4097 8193]; % kernel sizes
							     % not all used
level=num_levels;

if nargin < 5
  disp('est_pchirp2m: estimating init parameters from the two smoothed images')
  %p0 = est_pchirp2r(smooth(E1,boxsizes(level),'nan'),...
  p0 = est_fftpchirp2r(smooth(E1,boxsizes(level)),...  % smooth no longer uses
                    smooth(E2,boxsizes(level)),...  % third input argument
                    K         ...
                   );
  set(gcf,'Name',sprintf(' pchirp2r done (called by pchirp2m init.est.)'))
  num_levels=num_levels-1;  % if you pass initial guess, it _is_ the first level
end%if

if nargin < 6
  numlevelstoskip = 1;   % 0 would give all levels
  disp('est_pchirp2m: numlevels to skip defaulting to 1 (skipping top level)')
end%if

p_comp=p0;  % | initialize composite group parameters, either first guess or
            % | else it would have been specified by user input

% for level=num_levels:-1:1
for level=num_levels:-1:1+numlevelstoskip  % default down to level 2 only
  disp(sprintf('est_pchirp2m: now doing level %g',level));
  p_comp=est_pchirp2r(smooth(E1,boxsizes(level)),...
                      smooth(E2,boxsizes(level)),...
                      K,...
                      p_comp ...
                     ); % update composed p
  set(gcf,'Name',sprintf(' pchirp2r done (called by pchirp2m level %g)',level))
end%for level

if nargout == 1
  argument_p1 = p_comp;   % separate argument to avoid ans=... if nargout =0
end%if

if nargout == 0
  x=[0; 0; 1; 1];
  y=[0; 1; 0; 1];
  one=ones(4,1);
  xp=(p_comp(1)*x+p_comp(2)*y+p_comp(3))./(p_comp(7)*x+p_comp(8)+one);
  yp=(p_comp(4)*x+p_comp(5)*y+p_comp(6))./(p_comp(7)*x+p_comp(8)+one);
  disp('est_pchirp2m: 0 output args given so I''ll tell you a little more,')
  disp('              like where E1 corners correspond to E2')
  disp1=sprintf('       %+9.1f/%g %+9.1f/%g %+9.1f/%g %+9.1f/%g',...
                M*xp(1),M,M*xp(2),M,M*xp(3),M,M*xp(4),M);
  disp([disp1])
  disp1=sprintf('       %+9.1f/%g %+9.1f/%g %+9.1f/%g %+9.1f/%g',...
                N*yp(1),N,N*yp(2),N,N*yp(3),N,N*yp(4),N);
  disp(disp1)
  % not that we have left argument_p1 undefined so that if est_pchirp2m called
  % with no semicolon, ``ans'' will not appear
end%if

more on  % it usually is on (return it to what I think it was)