function [C0,C] = modecontours(x,y,z,varargin);

  % fix args
  (nargin > 3) || (OPT="");
  OPT      = structmerge(...
    "{zero:1e-4, inc:0.1, slvl: 0.2, decade_lvl:1, Z0: 0, debug: false}",varargin{});
  
  OPT.zero = abs(OPT.zero);
  OPT.slvl = abs(OPT.slvl);
  OPT.inc  = min(OPT.slvl, abs(OPT.inc));
  
  [N,M,K] = size(z);  
  z       = mat2cell(z,N,M,ones(1,K));

  if iscolumn(x) && isrow(y) || length(x) != M;
    z  = foreach(@transpose,z);
  endif
  
  edge = _detect_edge_(z{1}); 
  
  % align x/y to symm edge and cellularize
  [x,OPT.Nc,edge.x,cedge.i] = _align_xy_(x,edge.x);
  [y,OPT.Nr,edge.y,cedge.j] = _align_xy_(y,edge.y);
 
  % contour each z
  C0 = foreach(@(z) _contour_by_cell_(x,y,z,cedge,OPT), z);
  
  if K == 1;
    C0 = C0{1};
  endif
  
  if nargout > 1;
    C = parameterizeC(C0,OPT);
  endif
  
endfunction

###############################    CONTOURS   ##############################

function C = _contour_by_cell_(X,Y,Z,cedge,OPT);

  % clip FP precsion to prevent redundant contours levels
%  z = _prec_(z,OPT.zero/100);  ## no need, happens in solver
    
  if !OPT.Z0;
    abz = abs(Z);
    OPT.Z0 = max(abz(:));
  endif

  OPT.slvl *= OPT.Z0;
  OPT.inc  *= OPT.Z0;
  OPT.zero *= 10 * OPT.inc;
  OPT.debug && (count=1);

  Z = mat2cell(Z,OPT.Nr,OPT.Nc);  % convert z to cells according to discontinuities
  
  [N,M] = size(Z);

  for i = 1:M;
    x = X{i};
    for j = 1:N;                  OPT.debug && _show_(X,Y,Z,i,j,count++);
      
      y = Y{j};
      z = Z{j,i};
      ctype     = _celltype_([cedge.i(i), cedge.j(j)], z);
      
      if ctype.polar == 0;  
        z = _crease_(z);
      endif     
    
      [c,peaks] = mc_contour(x,y,z,ctype,OPT);
      
      % clean up contours
      c  = foreach(@(c) _fix_(c,peaks), c);
   
      if any(ctype.toInf);
        c = c(![c.closed]);
      endif
      
      c(foreach(@(c) c.length==0,c,true)) = [];
      C{j,i}  = _export_(c);
      
    endfor
  endfor

  function C = _export_(c);

    if any(ctype.toInf);
      data = [c.data];
      bbox = {_bounds_([data.x]),_bounds_([data.y])};
    else
      bbox = {x([1,end])(:),y([1,end])(:)}; 
    end
    
  % struct cell contours
    C.bbox     = xyz(bbox{});
    C.cell     = xyz((i==M)-(i==1),(j==N)-(j==1));
    C.clad     = ctype.toInf;
    C.peaks    = peaks;
    C.contours = c;
    C.id       = id('region');
  endfunction  
            
endfunction

################### EDGE SYMMETRY and POLAR DETECT #######################

function [edge] = _detect_edge_(z);   % look for symmetric edge

  abz   = abs(z);
  Z0    = max(abz(:));
  DZ    = Z0 ./ max(size(z));

  xedge = abz(:,[1,3,end-3,end]);    % sample x-edges
  yedge = abz([1,3,end-2,end],:);    %    and y-edges

  xymax = [max(xedge,[],1); max(yedge,[],2).'] ./ Z0;       % stack and normalize x and y edge 

  issymm = xymax(:,[1,4]) > DZ;                           % test for symmetric edges
  isasym = xymax(:,[1,4]) < 1e-14 & xymax(:,[2,3]) > DZ;  % test for antisymmetric edges
  
  % compute edge types: -1 for antisym, +1 for sym, Inf for evanescent to Inf
  edge_type = (issymm - isasym).\1;
  
  edge.x = edge_type(1,:);
  edge.y = edge_type(2,:);

endfunction

function ctype = _celltype_(cedge, z);
  
  ctype.toInf = (cedge == Inf);
  ctype.asym  = (cedge == -1);
  
  if any(ctype.toInf);      % get polarity of region
    [a b] = _bounds_(z([4,end-3],[2,end-1]));  % test near corners but avoid sampling the diagonal
    ctype.polar = (b>0) - (0>a);
    ctype.asym  = ctype.polar==0;              % in cladding regions, only 1 antisymm line is possible
    
  else  
    [a,b] = _bounds_(z(3:4,:));        % ***always*** look along x for polarity change 
    ctype.polar    = (b>0)-(0>a);      % even if asymmetry along x is known (need z polarity)
    ctype.asym(1) |= ctype.polar == 0;
    
    if !ctype.asym(2);  % conditionally do the same along y
      [c,d] = _bounds_(z(:,3:4)); 
      polar =  (d>0) - (0>c);
      ctype.polar == 0 || (ctype.polar = polar);  % do NOT overide ctype.polar == 0
      ctype.asym(2) = polar==0;
    endif 

    ctype.asym  = sum(ctype.asym);   % count the asymm edges
    
  endif
endfunction

#############################    CREASE   ##############################
 %%%%  TODO:  CREASE HORIZONTAL OR VERTICAL LINES ONLY
function z = _crease_(z,Z0);
% where z transistions from pos to neg, set to zero if close enough
%  IMPORTANT: do not assume z is symmetric exclusive along x or y

    az  = abs(z); 

    [N M] = size(z);
    N = 2:N-1;
    M = 2:M-1;
    
  % look for crease along row and col seperately
    cv = az(:,M-1) > az(:,M) & az(:,M) < az(:,M+1); 
    ch = az(N-1,:) > az(N,:) & az(N,:) < az(N+1,:); 

  % fill in gaps due to orthogonal crease
    cv = cv(N,:) | (cv(N-1,:) & cv(N+1,:));
    ch = ch(:,M);    % no need to look for gaps, but need to match size
    
    c = (cv | ch)  &  (az(N,M) < max(az(:)) * 1e-2);
    
  % restore size of c to that of z by cloning ends
    c = [c(:,1) c c(:,end)]; 
    c = [c(1,:);c;c(end,:)];
     
    z(c) = 0;
  
endfunction
    
###############################    FIX   ##############################

function C = _fix_(C,peaks);

    c = complex(C.x,C.y);
    C = rmfield(rmfield(C,'x'),'y');
    
    _fpequals_(c) && (c=c(1));

    N = length(c);
    closed = N > 2 && _fpequals_(c([1,end]));

    if closed;
      c = _ccw_circle_(c);
      CUTOFF = N/200;  % CUTOFF criteria for duplicate sample filtering
    else
      c = _ccw_curve_(c);
      CUTOFF = N/100;
    endif
    
    c = _dump_dups_(c,min(0.5, CUTOFF));   % IMPORTANT: do this after _ccw_circle_
    
    x = real(c);
    y = imag(c);
    
    C.bbox       = xyz(_bounds_(x),_bounds_(y));  
    C.localpeak  = _nearest_peak_(c,peaks);
    C.length     = sum(abs(diff(c))); 
    C.closed     = closed;
    C.id         = id('contour');
    C.type       = 'list';
    C.data       = struct('t',[0 length(c)-1],'x',x,'y',t);
    
endfunction

function peak = _nearest_peak_(c,peak);

    if length(peaks) < 2;  return; end;

    dist = abs(c' - complex([peak.x], [peak.y])); 
    [dist,idx] = max(dist);
    
    dist( any(idx == [1;numel(c)]) ) = Inf;  % ignore if max dist is end of c (curves wrong way)
    peak = peak([~,i] = min(dist));          % my peak has the shortest max distance
    
endfunction

function c = _ccw_circle_(c);
  
    [~,idx0] = max(real(c));                       % idx of right most point on circle (to be 0deg)
    (idx0 == 1) || (c = [c(idx0:end) c(2:idx0)]);  % remeber: c(1) == c(end)
    imag(c(2)) > imag(c(1)) || (c = flip(c));

endfunction

function [c] = _ccw_curve_(c);
  
%TODO: for speed up, change this to use slopes instead of arg
    ang  = arg( c([ceil(end/2),end]) - c(1) ) * 180/pi;
    abs(ang(1)) < 90  || (ang = mod(ang,360));
    (ang(1) < ang(2)) || (c = flip(c));       % ccw || (c=flip(c))
    
endfunction

function [c] = _dump_dups_(c, CUTOFF);
  
    dt   = abs(diff(c(:)'));
    dtn  = [dt(1:end-2);dt(3:end)];  % neigboring dt
    dt   = dt(2:end-1);
    
    toosmall = dt < max(dtn)*CUTOFF;
    smallest = dt < min(dtn);
    
    dup  = find(smallest & toosmall); % this point to offending interval
    
    if dup;
      dtn  = dtn(:,dup);                % get neighboring intervals
      dup += 1 + (1==([~,i]=max(dtn))); % get index of offending **sample**
      c(dup) = [];
    endif
    
endfunction  

function [x,N,edge,i_edge] = _align_xy_(x,edge);
  % this is going to align x with the symmetry edge
  %  BUT, it is sufficient and necesarry to change the numerical 
  % values WIHTOUT sorting (otherwise z has to be flipped too)
  
  _edge = edge;
  
  if edge(1) != Inf;
    x -= x(1);
  elseif edge(2) != Inf;
    x = x(end) - x;       % do NOT flip x
    edge = [edge(2) Inf]; % but DO flip edge
  end     
  
  % divide x into cells
  N  = diff( [0; find(x(1:end-1) == x(2:end))(:); length(x)] );

  if isrow(x);
    x = mat2cell(x,1,N);
    N = N.';
  else
    x = mat2cell(x,N,1);
  endif
  
  % record presence of symmetry edge for each cell
  % NOTE: if domain is single cell (_N == 1), then BOTH edges must 
  % be symmetric, prioritize recording the presence of anti-symmetric
  _N = length(N);  % number of cells
  
  if _N == 1;
    i_edge = min(_edge);
  else
    i_edge([1,_N]) = _edge;    % [_edge(1), 0,0..., _edge(end)]
  endif
  
endfunction

function [a,b] = _bounds_(x);
  
  x = x(:);
  a = min(x);
  b = max(x);
  nargout == 1 && (a = [a,b]);
  
endfunction

function istrue = _fpequals_(x);
  istrue = abs(x-x(1)) < abs(x(1)) * 1e-15;
endfunction


function _show_(x,y,z,i,j,count);
    f=gcf;figure(100); clf;
    imagesc(x{i},y{j},z{j,i});  axis xy;
    title(sprintf("[ %2d, %2d ] %2d",j,i,count)); 
    figure(f);
endfunction
   
   

################################ SCRAPS ###########################  




##      
##      function _crop_();
##
##          ok = absz >= OPT.zero/100;  % capture 1 row/column beyond limit
##
##          ix = any(ok,1); % test along y (dim=1)
##          iy = any(ok,2); % test along x (dim=2)
##          
##          % prevent cropping if it is only the edge
##          ix([1,end]) = ix([2,end-1]);
##          iy([1,end]) = iy([2,end-1]);
##          
##          y = y(iy);
##          x = x(ix);
##          z = z(iy,ix);
##          
##          absz  = abs(z);
##          [N M] = size(z);
##
##      endfunction

##
##function lvl = _append_(lvl,nlvl,INC);
##  % discard any level near other more critical ones
##  if lvl;
##    del = abs(nlvl - lvl.');
##    lvl = [lvl nlvl(all(del >= INC,1))];
##  else
##    lvl = nlvl;
##  end
##  
##endfunction    
## 
##function mlvl = _make_mlvl_(pv,N);
##  
##  N \= 1;
##  mlvl = N:N:(1-N);
##
##  ab = _bounds_([0 pv.z]);
##
##  if all(ab);
##    mlvl = (ab .* mlvl.')(:).';
##  else
##    mlvl = ab(1) + diff(ab).*mlvl;
##  endif
##  
##endfunction
##
##function elvl = _make_elvl_(epv,INC);
##%  elvl = [epv.z] - double([epv.type]) .* INC;  % level by edge peaks
##  elvl = [epv.z];
##endfunction
##

