## -*- texinfo -*- 
## @deftypefn {} {@var{ah} =} PlotLego3 (@var{z}, ...)
## @deftypefnx {} {@var{ah} =} PlotLego3 (@var{x}, @var{y}, @var{z},...)
## @deftypefnx {} {@var{ah} =} PlotLego3 (...,@var{options})
##    Create a 3d lego plot. Each matrix element in z is ploted as a rectangular
##    patch or as a cuboid (default). The function can handle missing data points
##    marked by NaN values. 
##
## Parameters:
##
##    @var{z} - 2d array with data to visualize
##
##    @var{x}, @var{y} - Both 1d arrays with x and y coordinates for each column
##              and row or 2d arrays with x,y coordinates for each data point.
##              Irregular grids can also be displayed correctly.
##
##    @var{options} - option value pairs:
##
##                    'curtain' - boolean: set to true to draw a curtain from
##                                each horizontal facet down to minimal z value.
##                                Defaults to true.
##
##                    'facecolor' - Defaults to 'interp'. Same values as 'facecolor' in patch() properties. 
##
##                    'edgecolor' - Defaults to 'k'. Same values as 'edgecolor' in patch() properties. 
##
##                    'floorcolor' - Color specification for the graph's floor to close shape
##                                   from bottom. This does also close all holes
##                                   in the date at minimal z value. Defaults to 'none'.
##
##                    'zmin'       - Minimal z value to use as bottom. Defaults to min(z(:)).
##
##                    All other options are passed to patch() function.
## @example
##  >> [x,y] = meshgrid(-2:3,-1:2);
##  >> z=10-x.*x-y.*y;
##  >> PlotLego3 (x,y,z)
##  >> PlotLego3 (x,y,z,'curtain',false)
##     % create an irregular grid
##  >> x1 = x+rand(size(x))*0.12-0.06;
##  >> y1 = y+rand(size(y))*0.12-0.06;
##  >> PlotLego3 (x1,y1,z)
##     % Fixed facet color
##  >> PlotLego3 (x1,y1,z,'facecolor','g')
##     % Holes in dataset
##     z1=z;z1(2:3,3)=NaN;
##  >> PlotLego3 (x1,y1,z1,'facecolor','g')
##  >> PlotLego3 (x1,y1,z1,'zmin',3)
##  >> PlotLego3 (x1,y1,z1,'floorcolor',[0.7,0.7,0.7])
## @end example
## @seealso{patch}
## @end deftypefn

## Author: Dr. Georg Wiora <gwiora@courage-khazaka.de>
## Created: 2019-07-10

function ah = PlotLego3 (varargin)
  p = inputParser ();
  p.addRequired('x',@(x) isnumeric(x) & (ismatrix(x) | isvector(x)));
  p.addOptional('y',[],@(y) isnumeric(y) & (ismatrix(y) | isvector(y)));
  p.addOptional('z',[],@(z) isnumeric(z) & (ismatrix(z) | isvector(z)));
  p.addParameter('curtain',true,@(x) isbool(x) & isscalar(x));
  p.addParameter('facecolor','interp');
  p.addParameter('edgecolor','k');
  p.addParameter('floorcolor','none');
  p.addParameter('zmin',NaN,@(zmin) isscalar(zmin) & isnumeric(zmin));
  % Keep additional parameters to forward to patch.
  p.KeepUnmatched = true;
  
  p.parse(varargin{:});
  
  % Get curtain argument
  bDrawCurtain = p.Results.curtain;
  % Get facet color argument (forwared to patch command)
  faceColor = p.Results.facecolor;
  % Get edge color argument (forwared to patch command)
  edgeColor = p.Results.edgecolor;
  % minz value
  zmin = p.Results.zmin;
  % Bottom color
  floorColor = p.Results.floorcolor;
  
  % Check input
  if isempty(p.Results.y) && isempty(p.Results.z)
    % first argument is z
    z=p.Results.x;
    % create index for x and y
    [x,y]=meshgrid(1:size(z,2),1:size(z,1));
  elseif ~isempty(p.Results.y) && ~isempty(p.Results.z)
    % if x and y are vectors make meshgrid
    if isvector(p.Results.x) && numel(p.Results.x)==size(p.Results.z,2) && isvector(p.Results.y) && numel(p.Results.y)==size(p.Results.z,1)
      [x,y] = meshgrid(p.Results.x,p.Results.y);
      % Get z data
      z = p.Results.z;
    elseif all(size(p.Results.z)==size(p.Results.x)) && all(size(p.Results.z)==size(p.Results.y))
      % Three matrices same size or point lists.
      x = p.Results.x;
      y = p.Results.y;
      z = p.Results.z;
    else
      error('x and y must be vectors matching size of z or matrices of same size as z!');
    endif

  else
    error('Call with two arguments is not supported!');
  endif

  % Create a patch for each matrix element
  % Points list: each data point needs 4 points for one patch
  pts = zeros(numel(z)*4,3);

  % Size of data array  
  [sy,sx]=size(z);
  
  % Interpolation grid for coordinates
  [xi,yi] = meshgrid(linspace(0.5,sx+0.5,sx+1),linspace(0.5,sy+0.5,sy+1));
  % Inter- and extrapolated coordinates of patch corners
  % Spline interpolation is ok, but a linear extrapolation would be much better.
  xc = interp2(x,xi,yi,'spline');
  yc = interp2(y,xi,yi,'spline');

  % Fix corners by linear extrapolate after spline interpolate
  % 1,1 corner
  fc = ([2*xc(1,2) - xc(1,3),2*yc(1,2)-yc(1,3)] + [2*xc(2,1) - xc(3,1),2*yc(2,1)-yc(3,1)])./2;
  xc(1,1)=fc(1);
  yc(1,1)=fc(2);
  % 1,end corner
  fc = ([2*xc(1,end-1) - xc(1,end-2),2*yc(1,end-1)-yc(1,end-2)] + [2*xc(2,end) - xc(3,end),2*yc(2,end)-yc(3,end)])./2;
  xc(1,end)=fc(1);
  yc(1,end)=fc(2);

  % end,end corner
  fc = ([2*xc(end,end-1) - xc(end,end-2),2*yc(end,end-1)-yc(end,end-2)] + [2*xc(end-1,end) - xc(end-2,end),2*yc(end-1,end)-yc(end-2,end)])./2;
  xc(end,end)=fc(1);
  yc(end,end)=fc(2);

  % end,1 corner
  fc = ([2*xc(end,2) - xc(end,3),2*yc(end,2)-yc(end,3)] + [2*xc(end-1,1) - xc(end-2,1),2*yc(end-1,1)-yc(end-2,1)])./2;
  xc(end,1)=fc(1);
  yc(end,1)=fc(2);

  
  % create points 4 times
  for i=1:4
    % Z is 4 times the same
    pts(i:4:end,3) = z(:);
  endfor
  % x and y
  % top left
  pts(1:4:end,1) = xc(1:end-1,1:end-1)(:);
  pts(1:4:end,2) = yc(1:end-1,1:end-1)(:);
  % top right
  pts(2:4:end,1) = xc(1:end-1,2:end)(:);
  pts(2:4:end,2) = yc(1:end-1,2:end)(:);
  % bottom right
  pts(3:4:end,1) = xc(2:end,2:end)(:);
  pts(3:4:end,2) = yc(2:end,2:end)(:);
  % bottom left
  pts(4:4:end,1) = xc(2:end,1:end-1)(:);
  pts(4:4:end,2) = yc(2:end,1:end-1)(:);

  % facet list
  % Each data point needs 1 facets with 4 points
  fac = reshape(1:size(pts,1),4,numel(z))';
  
  if bDrawCurtain
    % Zero plane for z if not preset
    if isnan(zmin)
      zmin = min(z(:));
    endif
    
    % Add vertical side 1 facets (left) + one column for right side
    ptsleft = zeros(numel(z)*4+sy*4,3);
    ptsleft(1:4:end,3) = horzcat(repmat(zmin,sy,1), z(:,1:end-1), repmat(zmin,sy,1))(:);
    ptsleft(2:4:end,3) = horzcat(z,z(:,end))(:);
    ptsleft(3:4:end,3) = horzcat(z,z(:,end))(:);
    ptsleft(4:4:end,3) = ptsleft(1:4:end,3);
    % replace all nan values by zmin
    ptsleft(isnan(ptsleft(:,3)),3) = zmin;
    
    ptsleft(1:4:end,1) = xc(1:end-1,:)(:);
    ptsleft(2:4:end,1) = xc(1:end-1,:)(:);
    ptsleft(3:4:end,1) = xc(2:end,:)(:);
    ptsleft(4:4:end,1) = xc(2:end,:)(:);

    ptsleft(1:4:end,2) = yc(1:end-1,:)(:);
    ptsleft(2:4:end,2) = yc(1:end-1,:)(:);
    ptsleft(3:4:end,2) = yc(2:end,:)(:);
    ptsleft(4:4:end,2) = yc(2:end,:)(:);
    
    % Vertex indices corrected
    facleft = reshape(1:size(ptsleft,1),4,numel(z)+sy)' + size(pts,1);
    
    % add to patch list
    pts = vertcat(pts,ptsleft);
    fac = vertcat(fac,facleft);
    
    % Add vertical side 1 facets (top) + one column for bottom side
    ptstop = zeros(numel(z)*4+size(z,2)*4,3);
    % z coordinates of vertical facets
    ptstop(1:4:end,3) = vertcat(repmat(zmin,1,sx), z(1:end-1,:), repmat(zmin,1,sx))(:);
    ptstop(2:4:end,3) = vertcat(z,z(end,:))(:);
    ptstop(3:4:end,3) = vertcat(z,z(end,:))(:);
    ptstop(4:4:end,3) = ptstop(1:4:end,3);
    % replace all nan values by zmin
    ptstop(isnan(ptstop(:,3)),3) = zmin;

    % x coordinates of vertical facets
    ptstop(1:4:end,1) = xc(:,2:end)(:);
    ptstop(2:4:end,1) = xc(:,2:end)(:);
    ptstop(3:4:end,1) = xc(:,1:end-1)(:);
    ptstop(4:4:end,1) = xc(:,1:end-1)(:);

    % y coordinates of vertical facets
    ptstop(1:4:end,2) = yc(:,2:end)(:);
    ptstop(2:4:end,2) = yc(:,2:end)(:);
    ptstop(3:4:end,2) = yc(:,1:end-1)(:);
    ptstop(4:4:end,2) = yc(:,1:end-1)(:);
    
    % Vertex indices corrected
    factop = reshape(1:size(ptstop,1),4,numel(z)+sx)' + size(pts,1);
    pts = vertcat(pts,ptstop);
    fac = vertcat(fac,factop);
  endif

  % plot facets and hand over additional parameters
  ah=patch('Faces', fac, 'Vertices', pts,'cdata',pts(:,3),'cdatamapping','scaled','facecolor',faceColor,'edgecolor',edgeColor, ...
        horzcat(fieldnames(p.Unmatched),struct2cell(p.Unmatched)){:}); 

  % Add a floor if requested
  if ~isstring(floorColor) || strcmpi(floorColor,'none')!=1
    % Add a floor vertex and facet
    ptsbot = zeros(2*sx+2*sy,3);
    % Place bottom a tini bit below zmin to avoid flicker
    ptsbot(:,3)=zmin-range(z(:))*0.001;
    % Top border
    ptsbot(1:sx+1,1) = xc(1,:)(:);
    ptsbot(1:sx+1,2) = yc(1,:)(:);
    cnt=sx+2;
    % right border
    ptsbot(cnt:cnt+sy-2,1) = xc(2:end-1,end)(:);
    ptsbot(cnt:cnt+sy-2,2) = yc(2:end-1,end)(:);
    cnt+=sy-1;
    % bottom border
    ptsbot(cnt:cnt+sx,1) = xc(end,end:-1:1)(:);
    ptsbot(cnt:cnt+sx,2) = yc(end,end:-1:1)(:);
    cnt+=sx+1;
    % left border
    ptsbot(cnt:cnt+sy-2,1) = xc(end-1:-1:2,1)(:);
    ptsbot(cnt:cnt+sy-2,2) = yc(end-1:-1:2,1)(:);
    
    % Facets: All points in one polygone
    facbot = 1:size(ptsbot,1);

    patch('Faces', facbot, 'Vertices', ptsbot,'cdata',ptsbot(:,3),'cdatamapping','scaled','facecolor',floorColor,'edgecolor',edgeColor, ...
        horzcat(fieldnames(p.Unmatched),struct2cell(p.Unmatched)){:});     
  endif
endfunction
