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 ##