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 ^

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

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005