function [bmask, b, threshM] = bigstarget(a, varargin)
% [bmask, b, threshM] = bigstarget(a, varargin)
% inputs :
%  a      : amplitude or sound dB matrix:
%               (nbins x nprofiles)  or (nbins x nbeams x nprofs)
% vararegin: name, value pairs:
% name                default                      what
% thresh                 15               size of amplitude bump for a
%                                         shallow bottom, (like old waterfall
%                                         plots).
%                      scalar             - constant for all bins, all beams
%                      vector             - same for all beams, linearly
%                                              interpolate the vector in
%                                              depth for a vertical profile
%                     structure           - field names are "beam1", "beam2",
%                                             (etc), one per beam.  Each field
%                                             is treated as scalar or vector
%                                             as above.  Use this option
%                                             if one beam is noisy and has
%                                             low threshold of detection
%
% nside               1                run a box filter in vertical with
%                                              this halfwidth
% offset                  3                how far away to look (in the
%                                              vertical) to define the bump
% ampfilt               tol              "median" filter, 3 pts (hardwired)
%                                             default [] (disabled)
%                                              tol: eg 10 for rssi or amp,
%                                                      4 for sound (dB)
% nmask               (none)            mask of bins to ignore  (noisy)
%                                             when finding max amp bin
%                                             (mask is ones where bad)
%
% outputs:
% bmask   : 1 at max amp bin, 0 otherwise
% b       : detrended amplitude, smoothed away from highest amplitudes
%           (values mostly near 0, but high at the 'bottom return')
% threshM : matrix of thresholds used

% NOTE: no NaNs in field; use 'nmask'

cfg.thresh = 25;
cfg.ampfilt = [];
cfg.nside = 1;
cfg.offset = 3;
cfg.nmask = [];


cfg = fillstruct(cfg, varargin, 0);

if length(size(a)) == 3
   [nd, nbeams, nprofs]=size(a);
else
   [nd, nprofs]=size(a);
   nbeams = 1;
end


% standardize threshold as a structure with one field per beam
thresh_err = 0;
if (isstruct(cfg.thresh))
   fnames = fieldnames(cfg.thresh);
   if length(fnames) ~= nbeams
      thresh_err = 1;
   end
   for beami = 1:nbeams
      if ~ismember(sprintf('beam%d', beami), fnames)
         thresh_err = 1;
         break
      end
   end
else
   for beami = 1:nbeams
      eval(sprintf('thresh.beam%d = cfg.thresh;', beami))
   end
   cfg.thresh = thresh;
end
if (thresh_err == 1)
   fprintf('error in ''thresh'' specification.  Structure was wrong\n')
   return
end
% now expand to fit the profile length
tt1 = cfg.thresh.beam1;
Nt = length(tt1);
ttn = table1( [ [0:Nt]'/(Nt)  [tt1(:); tt1(end)]], [0:nd]'/(nd+1));
ttn(end)=[];
threshM = repmat(ttn,1,nprofs);

for beami = 2:nbeams
   eval(sprintf('tt1 = cfg.thresh.beam%d;', beami));
   Nt = length(tt1);
   ttn = table1( [ [0:Nt]'/(Nt)  [tt1(:); tt1(end)]], [0:nd]'/(nd+1));
   ttn(end)=[];
   threshM = cat(3, threshM, repmat(ttn,1,nprofs));
end

if (nbeams > 1)
   threshM = permute(threshM, [1,3,2]);
end



if (~isempty(cfg.ampfilt))
   if length(cfg.ampfilt) == 2
      fprintf('ERROR in bigstarget: \n')
      fprintf('specify filter tolerance only\n')

      bmask = [];
      b = [];
      return
   end
   if (nprofs >= 3)
      % Make a square version if it was 3-D.
      if nbeams > 1,
         aa = reshape(a, nd * nbeams, nprofs);
      else
         aa=a;
      end
      
      % Average of profiles i-i and i+1.
      aamn = 0.5 * ([aa(:, 2:end), aa(:,end)] + ...
                    [aa(:, 1), aa(:,1:(end-1))]);
      % Find large positive outliers.
      badmask = (aa - aamn) > cfg.ampfilt;
      aa(badmask) = aamn(badmask);
      % Make it 3-D again if necessary.
      if nbeams > 1,
         a = reshape(aa, nd, nbeams, nprofs);
      else
         a = aa;
      end
   end
end

%run a boxcar to smooth the amplitude
if (nbeams == 1)
   switch cfg.nside
    case 1
     aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :)) / 3;
    case 2
     aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :) + ...
           a([1, 1, 1:(nd-2)], :) + a([3:nd, nd, nd], :)) / 5;
    case 3
     aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :) + ...
           a([1, 1, 1:(nd-2)], :) + a([3:nd, nd, nd], :) + ...
           a([1, 1, 1, 1:(nd-3)], :) + a([4:nd, nd, nd, nd], :)) / 7;
    case 0
     aa = a;
    otherwise
     error('nside must be in range 0-3');
   end
else
   switch cfg.nside
    case 1
     aa = (a + a([1, 1:(nd-1)], : ,:) + a([2:nd, nd], : ,:)) / 3;
    case 2
     aa = (a + a([1, 1:(nd-1)], : ,:) + a([2:nd, nd], : ,:) + ...
           a([1, 1, 1:(nd-2)], : ,:) + a([3:nd, nd, nd], : ,:)) / 5;
    case 3
     aa = (a + a([1, 1:(nd-1)], : ,:) + a([2:nd, nd], : ,:) + ...
           a([1, 1, 1:(nd-2)], : ,:) + a([3:nd, nd, nd], : ,:) + ...
           a([1, 1, 1, 1:(nd-3)], : ,:) + a([4:nd, nd, nd, nd], : ,:)) / 7;
    case 0
     aa = a;
    otherwise
     error('nside must be in range 0-3');
   end
end

%taper for endpoints
wu = ones(nd, 1);
wu(1:cfg.offset) = (0:(cfg.offset-1)).' / cfg.offset;
wu(nd:-1:(nd - cfg.offset + 1)) = 2 - wu(1:cfg.offset);
wd = flipud(wu);
%propagate to the right for nprofiles
wu = wu(:, ones(nprofs, 1));
wd = wd(:, ones(nprofs, 1));


iu = [ones(1, cfg.offset), 1:(nd - cfg.offset)];
id = [(1 + cfg.offset):nd, nd * ones(1, cfg.offset)];

%  unsmoothed version at the apex and smoothed versions
% at the shallower and deeper ends
if nbeams == 1
   b = 0.5 * a - .25 * (wu .* aa(iu, :) + wd .* aa(id, :));

   % Ignore the top (cfg.offset + cfg.nside) bins
   topmask = ones(nd, nprofs);
   topmask(1:(cfg.offset+ cfg.nside), :) = 0;
else
   wuM = wu;
   wdM = wd;
   for beami = 2:nbeams
      wuM = cat(3, wuM, wu);
      wdM = cat(3, wdM, wd);
   end
   wu = permute(wuM, [1,3,2]);
   wd = permute(wdM, [1,3,2]);

   b = 0.5 * a - .25 * (wu .* aa(iu, :, :) + wd .* aa(id, :, :));

   % Ignore the top (cfg.offset + cfg.nside) bins
   topmask = ones(nd, nbeams, nprofs);
   topmask(1:(cfg.offset+ cfg.nside), :, :) = 0;
end

if isempty(cfg.nmask)
   cfg.nmask = zeros(size(b));
end

bmask = (b > threshM) & topmask & ~cfg.nmask;



