%  [P,nextOctave] = process_octave(Octave,octSize)
%  takes Octave and first downsamples it by 2 to render it half its
%  original size.  then generates images at gaussians of successive
%  multiples of k=2.^(1/S) and computes Difference-of-Gaussians,
%  finding and collecting extrema along the way. Calls
%  find_dog_peaks and scales resulting points back to their
%  positions on the original image.

function [P,nextOctave] = process_octave(Octave,octSize)

   function sd = sigmadelta(k,sigma_cur)
   %% returns the sigma_delta (what we need to convolve with to
   %% achieve an blur of k more than sigma_cur
      sd = sigma_cur*sqrt(k.^2-1);
   end;

   function img = downsampleImageByTwo(Image) 
       % downsamples Image by 2 in each direction
       img = downsample(downsample(Image,2)',2)';
   end;

   P = [];

   %% get parameters
   [InitSigma,Scales,PeakThresh,R] = get_dog_params;
   
   k = 2.^(1/Scales);
   
   % compute our current "sigma" based upon passed in octSize
   % ie, 
   % octsize = 0.5, then sigma = 1/(2*0.5) = 1.0 * InitSigma; 
   % octsize = 0.25 then sigma = 1/(2*0.25) = 2*InitSigma   

%% Confusion
%% I was quite confused here. I was not sure whether we wanted
%% to re-compute the sigma relative to the previous octave, or whether we
%% essentially "started over" with InitSigma; OR, if we wanted to
%% set our starting sigma to 1.0 for all but the first octave...

%   sigma_cur = 1/(octSize*2) * InitSigma; 
%   sigma_cur = InitSigma; 
% if (octSize == 0.5)
       sigma_cur = InitSigma;
%     else
%       sigma_cur = 1.0;
%     end
   
   %% to create base sublevel, 
   %% downsample incoming Octave image to 1/2 size
   sublevel = downsampleImageByTwo(Octave);
   sublevels(:,:,1) = sublevel;
   
   % debug ------------------------------------------
   %%  sizesublevelstack = size(sublevels)
   % end of debug -----------------------------------

   %% first generate image levels
   for sublevel_i = 2:Scales+3,
     
     % debug ------------------------------------------
     blurgaussoutsize = ...
	 size(blurGauss(sublevel,sigmadelta(k,sigma_cur)))
     % end of debug -----------------------------------
     
     
     %  using incremental sigma 
     sublevels(:,:,sublevel_i) = ...
       blurGauss(sublevel,sigmadelta(k,sigma_cur));
     sublevel = sublevels(:,:,sublevel_i);
     sigma_cur = k*sigma_cur;

     if false,
     % debugging, to see what happens here
     sublevels(:,:,sublevel_i) = ...
             blurGauss(sublevel,k^(sublevel_i-1)*sigma_cur);
     sprintf('sigma_cur is %d',sigma_cur)
     end;
     %% end of debug
   end;

   %% compute difference of gaussians (DOG)s for this octave
   for dog_i = 1:Scales+2,
     dogs(:,:,dog_i) = sublevels(:,:,dog_i + 1)  ...
	 - sublevels(:,:,dog_i);
   end;

   size(dogs)
   
   P = find_dog_peaks(dogs,octSize);
   nextOctave = sublevels(:,:,Scales+1);
end