Warning, /utils/matlab/inpaint_nans.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
36eeb8dc0c Bayl*0001 function B=inpaint_nans(A,method)
0002 % inpaint_nans: in-paints over nans in an array
0003 % usage: B=inpaint_nans(A)
0004 %
0005 % solves approximation to one of several pdes to
0006 % interpolate and extrapolate holes
0007 %
0008 % Written by John D'Errico (woodchips@rochester.rr.com)
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
