% -- External Function: mk_2 = remove_ponds( mask );
%
% Take as input a two-dimensional array `mask' containing zeros or ones, with
%   values of zero representing land and values of one representing water. The
%   function discards all isolated bodies of water except for the largest one.
% Assumptions: The `largest body of water' in area is assumed to be the one with
%   the longest perimeter.
% Implementation: Individual bodies of water are identified with the intrinsic
%   function `contour'. Then, the small bodies of water are discarded one by one
%   using the intrinsic function `inpolygon'. (This is considerably faster than
%   running `inpolygon' over the main body of water.)
% Example of usage:
%   wdmk = ncread( 'roms_qck_0001.nc', 'wetdry_mask_rho' );
%   for irec = 1 : size( wdmk, 3 )
%     mk_2 = remove_ponds( wdmk(:, :, irec) );
%     figure(1); imagesc( mk_2' ); axis xy; axis equal; drawnow;
%   end; clear irec;
%
% Function arguments:
% `mask' (input ): Two-dimensional array of type single, double, or logical.
% `mk_2' (output): Two-dimensional array of type double (same size as mask)
%                    where ponds have been discarded.
% Tested R2021b Update 2 (9.11.0.1837725) 64-bit (glnxa64), GNU Octave v.6.2.0.
% psl@fastmail.net Pierre St-Laurent 2023-08-16.
function outp = remove_ponds( inpu )
  if ( islogical(  inpu ) )
    inpu = double( inpu );
  end
  nx   = size( inpu, 1 );
  ny   = size( inpu, 2 );

% nfig = length( findobj( 'type', 'figure' ) );
  figure(           'visible', 'off' )
% figure( nfig + 1, 'visible', 'off' )
    ccon = contour( inpu, 0.5 * [1, 1] );
  close;
% close(  nfig + 1 );
  ccon = ccon';

  ibeg = find( ccon(:, 1) == 0.5 );              % Beginning of  segment.
  iend = [ ibeg(2 : end); size( ccon, 1 ) + 1 ]; % End       of  segment.
  leng = ccon( ibeg, 2 );                        % Length    of  segment.
  main = find( leng == max( leng ) );            % Index of main segment.

  outp = inpu;
  [jj, ii] = meshgrid( 1 : ny, 1 : nx );
  for ibas = 1 : length( leng )
    if ( ibas == main )
      continue;
    end
    xpol = ccon( ibeg(ibas) + 1 : iend(ibas) - 1, 2 );
    ypol = ccon( ibeg(ibas) + 1 : iend(ibas) - 1, 1 );
    imin = max([ 1; floor(min(xpol)) - 1]);
    jmin = max([ 1; floor(min(ypol)) - 1]);
    imax = min([nx; ceil( max(xpol)) + 1]);
    jmax = min([ny; ceil( max(ypol)) + 1]);
    wrk1 = outp( imin : imax, jmin : jmax );
    inpo = inpolygon( ii(imin : imax, jmin : jmax), ...
                      jj(imin : imax, jmin : jmax), ...
                      xpol, ypol );
    wrk1( find( inpo ) ) = 0;
    outp( imin : imax, jmin : jmax ) = wrk1;
    clear xpol ypol imin imax jmin jmax wrk1 inpo;
  end; clear ibas;

  outp( find( outp > 0.5 & inpu < 0.5 ) ) = 0;
end