# nan_interp

## PURPOSE

Modifed Parsons Aug 03

## SYNOPSIS

function B=nan_interp(A,method)

## DESCRIPTION

     Modifed Parsons Aug 03

     nan_interp: in-paints over nans in an array
     usage: B=inpaint_nans(A)

     solves approximation to one of several pdes to
     interpolate and extrapolate holes

     arguments (input):
       A - nxm array with some NaNs to be filled in

       method - (OPTIONAL) scalar numeric flag - specifies
           which approach (or physical metaphor to use
           for the interpolation.) All methods are capable
           of extrapolation, some are better than others.
           There are also speed differences, as well as
           accuracy differences for smooth surfaces.

           methods {0,1,2} use a simple plate metaphor
           methods {3} use a better plate equation,
                       but will be slower
           methods 4 use a spring metaphor

           method == 0 --> (DEFAULT) see method 1, but
             this method does not build as large of a
             linear system in the case of only a few
             NaNs in a large array.
             Extrapolation behavior is linear.

           method == 1 --> simple approach, applies del^2
             over the entire array, then drops those parts
             of the array which do not have any contact with
             NaNs. Uses a least squares approach, but it
             does not touch existing points.
             In the case of small arrays, this method is
             quite fast as it does very little extra work.
             Extrapolation behavior is linear.

           method == 2 --> uses del^2, but solving a direct
             linear system of equations for nan elements.
             This method will be the fastest possible for
             large systems since it uses the sparsest
             possible system of equations. Not a least
             squares approach, so it may be least robust
             to noise on the boundaries of any holes.
             This method will also be least able to
             interpolate accurately for smooth surfaces.
             Extrapolation behavior is linear.

           method == 3 --+ See method 0, but uses del^4 for
             the interpolating operator. This may result
             in more accurate interpolations, at some cost
             in speed.

           method == 4 --+ Uses a spring metaphor. Assumes
             springs (with a nominal length of zero)
             connect each node with every neighbor
             (horizontally, vertically and diagonally)
             Since each node tries to be like its neighbors,
             extrapolation is as a constant function where
             this is consistent with the neighboring nodes.

     arguments (output):
       B - nxm array with NaNs replaced

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SUBFUNCTIONS

- function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)

## SOURCE CODE

    0001 function B=nan_interp(A,method)
    0002 % Modifed Parsons Aug 03
    0003 %
    0004 % nan_interp: in-paints over nans in an array
    0005 % usage: B=inpaint_nans(A)
    0006 %
    0007 % solves approximation to one of several pdes to
    0008 % interpolate and extrapolate holes
    0009 %
    0010 % arguments (input):
    0011 %   A - nxm array with some NaNs to be filled in
    0012 %
    0013 %   method - (OPTIONAL) scalar numeric flag - specifies
    0014 %       which approach (or physical metaphor to use
    0015 %       for the interpolation.) All methods are capable
    0016 %       of extrapolation, some are better than others.
    0017 %       There are also speed differences, as well as
    0018 %       accuracy differences for smooth surfaces.
    0019 %
    0020 %       methods {0,1,2} use a simple plate metaphor
    0021 %       methods {3} use a better plate equation,
    0022 %                   but will be slower
    0023 %       methods 4 use a spring metaphor
    0024 %
    0025 %       method == 0 --> (DEFAULT) see method 1, but
    0026 %         this method does not build as large of a
    0027 %         linear system in the case of only a few
    0028 %         NaNs in a large array.
    0029 %         Extrapolation behavior is linear.
    0030 %
    0031 %       method == 1 --> simple approach, applies del^2
    0032 %         over the entire array, then drops those parts
    0033 %         of the array which do not have any contact with
    0034 %         NaNs. Uses a least squares approach, but it
    0035 %         does not touch existing points.
    0036 %         In the case of small arrays, this method is
    0037 %         quite fast as it does very little extra work.
    0038 %         Extrapolation behavior is linear.
    0039 %
    0040 %       method == 2 --> uses del^2, but solving a direct
    0041 %         linear system of equations for nan elements.
    0042 %         This method will be the fastest possible for
    0043 %         large systems since it uses the sparsest
    0044 %         possible system of equations. Not a least
    0045 %         squares approach, so it may be least robust
    0046 %         to noise on the boundaries of any holes.
    0047 %         This method will also be least able to
    0048 %         interpolate accurately for smooth surfaces.
    0049 %         Extrapolation behavior is linear.
    0050 %
    0051 %       method == 3 --+ See method 0, but uses del^4 for
    0052 %         the interpolating operator. This may result
    0053 %         in more accurate interpolations, at some cost
    0054 %         in speed.
    0055 %
    0056 %       method == 4 --+ Uses a spring metaphor. Assumes
    0057 %         springs (with a nominal length of zero)
    0058 %         connect each node with every neighbor
    0059 %         (horizontally, vertically and diagonally)
    0060 %         Since each node tries to be like its neighbors,
    0061 %         extrapolation is as a constant function where
    0062 %         this is consistent with the neighboring nodes.
    0063 %
    0064 %
    0065 % arguments (output):
    0066 %   B - nxm array with NaNs replaced
    0067
    0068 % I always need to know which elements are NaN,
    0069 % and what size the array is for any method
    0070 [n,m]=size(A);
    0071 nm=n*m;
    0072 k=isnan(A(:));
    0073
    0074 % list the nodes which are known, and which will
    0075 % be interpolated
    0076 nan_list=find(k);
    0077 known_list=find(~k);
    0078
    0079 % how many nans overall
    0080 nan_count=length(nan_list);
    0081
    0082 % convert NaN indices to (r,c) form
    0083 % nan_list==find(k) are the unrolled (linear) indices
    0084 % (row,column) form
    0085 [nr,nc]=ind2sub([n,m],nan_list);
    0086
    0087 % both forms of index in one array:
    0088 % column 1 == unrolled index
    0089 % column 2 == row index
    0090 % column 3 == column index
    0091 nan_list=[nan_list,nr,nc];
    0092
    0093 % supply default method
    0094 if (nargin<2)|isempty(method)
    0095   method = 0;
    0096 elseif ~ismember(method,0:4)
    0097   error 'If supplied, method must be one of: {0,1,2,3,4}.'
    0098 end
    0099
    0100 % for different methods
    0101 switch method
    0102  case 0
    0103   % The same as method == 1, except only work on those
    0104   % elements which are NaN, or at least touch a NaN.
    0105
    0106   % horizontal and vertical neighbors only
    0107   talks_to = [-1 0;0 -1;1 0;0 1];
    0108   neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
    0109
    0110   % list of all nodes we have identified
    0111   all_list=[nan_list;neighbors_list];
    0112
    0113   % generate sparse array with second partials on row
    0114   % variable for each element in either list, but only
    0115   % for those nodes which have a row index > 1 or < n
    0116   L = find((all_list(:,2) > 1) & (all_list(:,2) < n));
    0117   nl=length(L);
    0118   if nl>0
    0119     fda=sparse(repmat(all_list(L,1),1,3), ...
    0120       repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
    0121       repmat([1 -2 1],nl,1),nm,nm);
    0122   else
    0123     fda=spalloc(n*m,n*m,size(all_list,1)*5);
    0124   end
    0125
    0126   % 2nd partials on column index
    0127   L = find((all_list(:,3) > 1) & (all_list(:,3) < m));
    0128   nl=length(L);
    0129   if nl>0
    0130     fda=fda+sparse(repmat(all_list(L,1),1,3), ...
    0131       repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
    0132       repmat([1 -2 1],nl,1),nm,nm);
    0133   end
    0134
    0135   % eliminate knowns
    0136   rhs=-fda(:,known_list)*A(known_list);
    0137   k=find(any(fda(:,nan_list(:,1)),2));
    0138
    0139   % and solve...
    0140   B=A;
    0141   B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
    0142
    0143  case 1
    0144   % least squares approach with del^2. Build system
    0145   % for every array element as an unknown, and then
    0146   % eliminate those which are knowns.
    0147
    0148   % Build sparse matrix approximating del^2 for
    0149   % every element in A.
    0150   % Compute finite difference for second partials
    0151   % on row variable first
    0152   [i,j]=ndgrid(2:(n-1),1:m);
    0153   ind=i(:)+(j(:)-1)*n;
    0154   np=(n-2)*m;
    0155   fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
    0156       repmat([1 -2 1],np,1),n*m,n*m);
    0157
    0158   % now second partials on column variable
    0159   [i,j]=ndgrid(1:n,2:(m-1));
    0160   ind=i(:)+(j(:)-1)*n;
    0161   np=n*(m-2);
    0162   fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
    0163       repmat([1 -2 1],np,1),nm,nm);
    0164
    0165   % eliminate knowns
    0166   rhs=-fda(:,known_list)*A(known_list);
    0167   k=find(any(fda(:,nan_list),2));
    0168
    0169   % and solve...
    0170   B=A;
    0171   B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
    0172
    0173  case 2
    0174   % Direct solve for del^2 BVP across holes
    0175
    0176   % generate sparse array with second partials on row
    0177   % variable for each nan element, only for those nodes
    0178   % which have a row index > 1 or < n
    0179   L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n));
    0180   nl=length(L);
    0181   if nl>0
    0182     fda=sparse(repmat(nan_list(L,1),1,3), ...
    0183       repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
    0184       repmat([1 -2 1],nl,1),n*m,n*m);
    0185   else
    0186     fda=spalloc(n*m,n*m,size(nan_list,1)*5);
    0187   end
    0188
    0189   % 2nd partials on column index
    0190   L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m));
    0191   nl=length(L);
    0192   if nl>0
    0193     fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
    0194       repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
    0195       repmat([1 -2 1],nl,1),n*m,n*m);
    0196   end
    0197
    0198   % fix boundary conditions at extreme corners
    0199   % of the array in case there were nans there
    0200   if ismember(1,nan_list(:,1))
    0201     fda(1,[1 2 n+1])=[-2 1 1];
    0202   end
    0203   if ismember(n,nan_list(:,1))
    0204     fda(n,[n, n-1,n+n])=[-2 1 1];
    0205   end
    0206   if ismember(nm-n+1,nan_list(:,1))
    0207     fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
    0208   end
    0209   if ismember(nm,nan_list(:,1))
    0210     fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
    0211   end
    0212
    0213   % eliminate knowns
    0214   rhs=-fda(:,known_list)*A(known_list);
    0215
    0216   % and solve...
    0217   B=A;
    0218   k=nan_list(:,1);
    0219   B(k)=fda(k,k)\rhs(k);
    0220
    0221  case 3
    0222   % The same as method == 0, except uses del^4 as the
    0223   % interpolating operator.
    0224
    0225   % del^4 template of neighbors
    0226   talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
    0227       0 1;0 2;1 -1;1 0;1 1;2 0];
    0228   neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
    0229
    0230   % list of all nodes we have identified
    0231   all_list=[nan_list;neighbors_list];
    0232
    0233   % generate sparse array with del^4, but only
    0234   % for those nodes which have a row & column index
    0235   % >= 3 or <= n-2
    0236   L = find( (all_list(:,2) >= 3) & ...
    0237             (all_list(:,2) <= (n-2)) & ...
    0238             (all_list(:,3) >= 3) & ...
    0239             (all_list(:,3) <= (m-2)));
    0240   nl=length(L);
    0241   if nl>0
    0242     % do the entire template at once
    0243     fda=sparse(repmat(all_list(L,1),1,13), ...
    0244         repmat(all_list(L,1),1,13) + ...
    0245         repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ...
    0246         repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm);
    0247   else
    0248     fda=spalloc(n*m,n*m,size(all_list,1)*5);
    0249   end
    0250
    0251   % on the boundaries, reduce the order around the edges
    0252   L = find((((all_list(:,2) == 2) | ...
    0253              (all_list(:,2) == (n-1))) & ...
    0254             (all_list(:,3) >= 2) & ...
    0255             (all_list(:,3) <= (m-1))) | ...
    0256            (((all_list(:,3) == 2) | ...
    0257              (all_list(:,3) == (m-1))) & ...
    0258             (all_list(:,2) >= 2) & ...
    0259             (all_list(:,2) <= (n-1))));
    0260   nl=length(L);
    0261   if nl>0
    0262     fda=fda+sparse(repmat(all_list(L,1),1,5), ...
    0263       repmat(all_list(L,1),1,5) + ...
    0264         repmat([-n,-1,0,+1,n],nl,1), ...
    0265       repmat([1 1 -4 1 1],nl,1),nm,nm);
    0266   end
    0267
    0268   L = find( ((all_list(:,2) == 1) | ...
    0269              (all_list(:,2) == n)) & ...
    0270             (all_list(:,3) >= 2) & ...
    0271             (all_list(:,3) <= (m-1)));
    0272   nl=length(L);
    0273   if nl>0
    0274     fda=fda+sparse(repmat(all_list(L,1),1,3), ...
    0275       repmat(all_list(L,1),1,3) + ...
    0276         repmat([-n,0,n],nl,1), ...
    0277       repmat([1 -2 1],nl,1),nm,nm);
    0278   end
    0279
    0280   L = find( ((all_list(:,3) == 1) | ...
    0281              (all_list(:,3) == m)) & ...
    0282             (all_list(:,2) >= 2) & ...
    0283             (all_list(:,2) <= (n-1)));
    0284   nl=length(L);
    0285   if nl>0
    0286     fda=fda+sparse(repmat(all_list(L,1),1,3), ...
    0287       repmat(all_list(L,1),1,3) + ...
    0288         repmat([-1,0,1],nl,1), ...
    0289       repmat([1 -2 1],nl,1),nm,nm);
    0290   end
    0291
    0292   % eliminate knowns
    0293   rhs=-fda(:,known_list)*A(known_list);
    0294   k=find(any(fda(:,nan_list(:,1)),2));
    0295
    0296   % and solve...
    0297   B=A;
    0298   B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
    0299
    0300  case 4
    0301   % Spring analogy
    0302   % interpolating operator.
    0303
    0304   % list of all springs between a node and a horizontal
    0305   % or vertical neighbor
    0306   hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1];
    0307   hv_springs=[];
    0308   for i=1:4
    0309     hvs=nan_list+repmat(hv_list(i,:),nan_count,1);
    0310     k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m);
    0311     hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]];
    0312   end
    0313
    0314   % delete replicate springs
    0315   hv_springs=unique(sort(hv_springs,2),'rows');
    0316
    0317   % build sparse matrix of connections, springs
    0318   % connecting diagonal neighbors are weaker than
    0319   % the horizontal and vertical springs
    0320   nhv=size(hv_springs,1);
    0321   springs=sparse(repmat((1:nhv)',1,2),hv_springs, ...
    0322      repmat([1 -1],nhv,1),nhv,nm);
    0323
    0324   % eliminate knowns
    0325   rhs=-springs(:,known_list)*A(known_list);
    0326
    0327   % and solve...
    0328   B=A;
    0329   B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs;
    0330
    0331 end
    0332
    0333 % ====================================================
    0334 %      end of main function
    0335 % ====================================================
    0336 % ====================================================
    0337 %      begin subfunctions
    0338 % ====================================================
    0339 function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)
    0340 % identify_neighbors: identifies all the neighbors of
    0341 %   those nodes in nan_list, not including the nans
    0342 %   themselves
    0343 %
    0344 % arguments (input):
    0345 %  n,m - scalar - [n,m]=size(A), where A is the
    0346 %      array to be interpolated
    0347 %  nan_list - array - list of every nan element in A
    0348 %      nan_list(i,1) == linear index of i'th nan element
    0349 %      nan_list(i,2) == row index of i'th nan element
    0350 %      nan_list(i,3) == column index of i'th nan element
    0351 %  talks_to - px2 array - defines which nodes communicate
    0352 %      with each other, i.e., which nodes are neighbors.
    0353 %
    0354 %      talks_to(i,1) - defines the offset in the row
    0355 %                      dimension of a neighbor
    0356 %      talks_to(i,2) - defines the offset in the column
    0357 %                      dimension of a neighbor
    0358 %
    0359 %      For example, talks_to = [-1 0;0 -1;1 0;0 1]
    0360 %      means that each node talks only to its immediate
    0361 %      neighbors horizontally and vertically.
    0362 %
    0363 % arguments(output):
    0364 %  neighbors_list - array - list of all neighbors of
    0365 %      all the nodes in nan_list
    0366
    0367 if ~isempty(nan_list)
    0368   % use the definition of a neighbor in talks_to
    0369   nan_count=size(nan_list,1);
    0370   talk_count=size(talks_to,1);
    0371
    0372   nn=zeros(nan_count*talk_count,2);
    0373   j=[1,nan_count];
    0374   for i=1:talk_count
    0375     nn(j(1):j(2),:)=nan_list(:,2:3) + ...
    0376         repmat(talks_to(i,:),nan_count,1);
    0377     j=j+nan_count;
    0378   end
    0379
    0380   % drop those nodes which fall outside the bounds of the
    0381   % original array
    0382   L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m);
    0383   nn(L,:)=[];
    0384
    0385   % form the same format 3 column array as nan_list
    0386   neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn];
    0387
    0388   % delete replicates in the neighbors list
    0389   neighbors_list=unique(neighbors_list,'rows');
    0390
    0391   % and delete those which are also in the list of NaNs.
    0392   neighbors_list=setdiff(neighbors_list,nan_list,'rows');
    0393
    0394 else
    0395   neighbors_list=[];
    0396 end
    0397
    0398
    0399
    0400
    0401
    0402
    0403
    0404
    0405
    0406
    0407
    0408
    0409
    0410
    0411
    0412