NanSurf.zip

Make surface plot when X,Y,Z contain NaNs, taking care of NaN-discontinuities

% NANSURF(X,Y,Z,varargin) - makes a trisurf plot of X,Y,Z choosing a triangulation
% that minimizes zigzag edges due to NaN's in Z. No data interpolation takes place.
%
% inputs : X,Y,Z 2D matrixes representing a surface
% outputs: ph plot handle
% trgList the chosen triangulation
%
% NANSURF(X,Y,Z) plots the colored parametric surface defined by matrix arguments.
% it works the same as trisurf but choosing triangulation according to NaN's in Z.
%
% NANSURF(X,Y,Z,'PropertyName',PropertyValue,...) plots the data
% stored in the vectors X and Y, according to Property/Value pairs:
%
% QuadEdgesOnly
% 1 : only show quad edges as in surf, hides the triangle edges
% 0 : (Default) plots all triangle edges (can be changed outside with 'edgecolor')
%
% trgExtend
% 1 : (Default) enable extending triangles to better reduce NaN-edge-aliasing
% in close to horizonal/vertical lines.
% 0 : maximizes nr of rendered triangles without creating/extending them
% should only be used to compare/debug triangulation
% DBGTrg
% 1 : debug plot showing which triangles are affected by NaNs
% 0 : (Default) no debug plots
%
% IDEA BEHIND:
% When a face, e.g. quadruple of vertexes in a meshgrid, contain 1 Nan, surf does not render it.
% Choose triangles to maximize # rendered ones, maintaining correct face orientation.
% Correct face orientation (based on gradient) improves color rendering.
%
% METHOD:
% (1) NaNsurf makes a triangulation which maximizes the number of valid triangles
% a valid triangle contains no NaNs. output: TRG, a list of L/R triangles based on NaNs in Z
% (2) TRG list from step (1) is used as input to trisurf / patch: trisurf(TRG,X,Y,Z)
%
% EXAMPLE:
% plotSpec = {'edgecolor','k','facecolor','interp'};
% NaNsurf(X,Y,Z,plotSpec{:});
%
% EXAMPLE:
% NaNsurf(X,Y,Z,plotSpec{:},'trgExtend',1,'DBGTrg',1);
%
% EXAMPLE:
% M = 22; N = 34;
% Mh = (M+1)/2; Nh = (N+1)/2;
% [X,Y] = meshgrid((1:N)-Nh,(1:M)-Mh);
% Rn = sqrt((X/Nh).^2+(Y/Mh).^2);
% Z = sin(3*pi*Rn)./(3*pi*Rn);
% Z(find(X.^2+Y.^2 > Mh^2)) = NaN;
% NaNsurf(X,Y,Z,'edgecolor','k','facecolor','interp')
% axis([-Mh Mh -Nh Nh -.5 1])
%
% rev.1.8, Massimo Ciacci, September 20, 2019