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
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