function [pntr, trir] = refine(pnt, tri, method, varargin)

% REFINE a 3D surface that is described by a triangulation
%
% Use as
%   [pnt, tri] = refine(pnt, tri)
%   [pnt, tri] = refine(pnt, tri, 'banks', selection)
%   [pnt, tri] = refine(pnt, tri, 'updown', numtri)
%
% If no method is specified, the default is to refine the mesh globally by bisecting
% each edge according to the algorithm described in Banks, 1983.
%
% The Banks method allows the specification of a subset of triangles to be refined
% according to Banks' algorithm. Adjacent triangles will be gracefully dealt with.
%
% The alternative 'updown' method refines the mesh a couple of times
% using Banks' algorithm, followed by a downsampling using the REDUCEPATCH
% function.
%
% The Banks method is a memory efficient implementation which remembers the
% previously inserted vertices. The refinement algorithm executes in linear
% time with the number of triangles. It is mentioned in
% http://www.cs.rpi.edu/~flaherje/pdf/fea8.pdf, which also contains the original
% reference.

% Copyright (C) 2002-2014, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: refine.m 10197 2015-02-11 09:35:58Z roboos $

if nargin<3
  method = 'banks';
end

if strcmp(method, 'banks') && nargin==4
  % FIXME this is work in progress
  keyboard
end

switch lower(method)
  case 'banks'
    npnt   = size(pnt,1);
    ntri   = size(tri,1);
    insert = spalloc(3*npnt,3*npnt,3*ntri);
    
    trir  = zeros(4*ntri,3);      % allocate memory for the new triangles
    pntr  = zeros(npnt+3*ntri,3); % allocate memory for the maximum number of new vertices
    pntr(1:npnt,:) = pnt;         % insert the original vertices
    current = npnt;
    
    for i=1:ntri
      
      if ~insert(tri(i,1),tri(i,2))
        current = current + 1;
        pntr(current,:) = (pnt(tri(i,1),:) + pnt(tri(i,2),:))/2;
        insert(tri(i,1),tri(i,2)) = current;
        insert(tri(i,2),tri(i,1)) = current;
        v12 = current;
      else
        v12 = insert(tri(i,1),tri(i,2));
      end
      
      if ~insert(tri(i,2),tri(i,3))
        current = current + 1;
        pntr(current,:) = (pnt(tri(i,2),:) + pnt(tri(i,3),:))/2;
        insert(tri(i,2),tri(i,3)) = current;
        insert(tri(i,3),tri(i,2)) = current;
        v23 = current;
      else
        v23 = insert(tri(i,2),tri(i,3));
      end
      
      if ~insert(tri(i,3),tri(i,1))
        current = current + 1;
        pntr(current,:) = (pnt(tri(i,3),:) + pnt(tri(i,1),:))/2;
        insert(tri(i,3),tri(i,1)) = current;
        insert(tri(i,1),tri(i,3)) = current;
        v31 = current;
      else
        v31 = insert(tri(i,3),tri(i,1));
      end
      
      % add the 4 new triangles with the correct indices
      trir(4*(i-1)+1, :) = [tri(i,1) v12 v31];
      trir(4*(i-1)+2, :) = [tri(i,2) v23 v12];
      trir(4*(i-1)+3, :) = [tri(i,3) v31 v23];
      trir(4*(i-1)+4, :) = [v12 v23 v31];
      
    end
    
    % remove the space for the vertices that was not used
    pntr = pntr(1:current, :);
    
  case 'updown'
    ntri = size(tri,1);
    while ntri<varargin{1}
      % increase the number of triangles by a factor of 4
      [pnt, tri] = refine(pnt, tri, 'banks');
      ntri = size(tri,1);
    end
    % reduce number of triangles using MATLAB function
    [trir, pntr] = reducepatch(tri, pnt, varargin{1});
    
  otherwise
    error(['unsupported method: ' method]);
end

