INPAINT_NANS: in-paints over nans in an array
usage: B=INPAINT_NANS(A) % default method
usage: B=INPAINT_NANS(A,method) % specify method used
Solves approximation to one of several pdes to
interpolate and extrapolate holes in an array
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.
method 3 uses a better plate equation,
but may be much slower and uses
more memory.
method 4 uses a spring metaphor.
method 5 is an 8 neighbor average, with no
rationale behind it compared to the
other methods. I do not recommend
its use.
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 modify known values.
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.
Note: method 2 has problems in 1-d, so this
method is disabled for vector inputs.
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.
method == 5 --+ See method 2, but use an average
of the 8 nearest neighbors to any element.
This method is NOT recommended for use.
arguments (output):
B - nxm array with NaNs replaced
Example:
[x,y] = meshgrid(0:.01:1);
z0 = exp(x+y);
znan = z0;
znan(20:50,40:70) = NaN;
znan(30:90,5:10) = NaN;
znan(70:75,40:90) = NaN;
z = inpaint_nans(znan);
See also: griddata, interp1
Author: John D'Errico
e-mail address: woodchips@rochester.rr.com
Release: 2
Release date: 4/15/060001 function B=inpaint_nans(A,method) 0002 % INPAINT_NANS: in-paints over nans in an array 0003 % usage: B=INPAINT_NANS(A) % default method 0004 % usage: B=INPAINT_NANS(A,method) % specify method used 0005 % 0006 % Solves approximation to one of several pdes to 0007 % interpolate and extrapolate holes in an array 0008 % 0009 % arguments (input): 0010 % A - nxm array with some NaNs to be filled in 0011 % 0012 % method - (OPTIONAL) scalar numeric flag - specifies 0013 % which approach (or physical metaphor to use 0014 % for the interpolation.) All methods are capable 0015 % of extrapolation, some are better than others. 0016 % There are also speed differences, as well as 0017 % accuracy differences for smooth surfaces. 0018 % 0019 % methods {0,1,2} use a simple plate metaphor. 0020 % method 3 uses a better plate equation, 0021 % but may be much slower and uses 0022 % more memory. 0023 % method 4 uses a spring metaphor. 0024 % method 5 is an 8 neighbor average, with no 0025 % rationale behind it compared to the 0026 % other methods. I do not recommend 0027 % its use. 0028 % 0029 % method == 0 --> (DEFAULT) see method 1, but 0030 % this method does not build as large of a 0031 % linear system in the case of only a few 0032 % NaNs in a large array. 0033 % Extrapolation behavior is linear. 0034 % 0035 % method == 1 --> simple approach, applies del^2 0036 % over the entire array, then drops those parts 0037 % of the array which do not have any contact with 0038 % NaNs. Uses a least squares approach, but it 0039 % does not modify known values. 0040 % In the case of small arrays, this method is 0041 % quite fast as it does very little extra work. 0042 % Extrapolation behavior is linear. 0043 % 0044 % method == 2 --> uses del^2, but solving a direct 0045 % linear system of equations for nan elements. 0046 % This method will be the fastest possible for 0047 % large systems since it uses the sparsest 0048 % possible system of equations. Not a least 0049 % squares approach, so it may be least robust 0050 % to noise on the boundaries of any holes. 0051 % This method will also be least able to 0052 % interpolate accurately for smooth surfaces. 0053 % Extrapolation behavior is linear. 0054 % 0055 % Note: method 2 has problems in 1-d, so this 0056 % method is disabled for vector inputs. 0057 % 0058 % method == 3 --+ See method 0, but uses del^4 for 0059 % the interpolating operator. This may result 0060 % in more accurate interpolations, at some cost 0061 % in speed. 0062 % 0063 % method == 4 --+ Uses a spring metaphor. Assumes 0064 % springs (with a nominal length of zero) 0065 % connect each node with every neighbor 0066 % (horizontally, vertically and diagonally) 0067 % Since each node tries to be like its neighbors, 0068 % extrapolation is as a constant function where 0069 % this is consistent with the neighboring nodes. 0070 % 0071 % method == 5 --+ See method 2, but use an average 0072 % of the 8 nearest neighbors to any element. 0073 % This method is NOT recommended for use. 0074 % 0075 % 0076 % arguments (output): 0077 % B - nxm array with NaNs replaced 0078 % 0079 % 0080 % Example: 0081 % [x,y] = meshgrid(0:.01:1); 0082 % z0 = exp(x+y); 0083 % znan = z0; 0084 % znan(20:50,40:70) = NaN; 0085 % znan(30:90,5:10) = NaN; 0086 % znan(70:75,40:90) = NaN; 0087 % 0088 % z = inpaint_nans(znan); 0089 % 0090 % 0091 % See also: griddata, interp1 0092 % 0093 % Author: John D'Errico 0094 % e-mail address: woodchips@rochester.rr.com 0095 % Release: 2 0096 % Release date: 4/15/06 0097 0098 0099 % I always need to know which elements are NaN, 0100 % and what size the array is for any method 0101 [n,m]=size(A); 0102 A=A(:); 0103 nm=n*m; 0104 k=isnan(A(:)); 0105 0106 % list the nodes which are known, and which will 0107 % be interpolated 0108 nan_list=find(k); 0109 known_list=find(~k); 0110 0111 % how many nans overall 0112 nan_count=length(nan_list); 0113 0114 % convert NaN indices to (r,c) form 0115 % nan_list==find(k) are the unrolled (linear) indices 0116 % (row,column) form 0117 [nr,nc]=ind2sub([n,m],nan_list); 0118 0119 % both forms of index in one array: 0120 % column 1 == unrolled index 0121 % column 2 == row index 0122 % column 3 == column index 0123 nan_list=[nan_list,nr,nc]; 0124 0125 % supply default method 0126 if (nargin<2) || isempty(method) 0127 method = 0; 0128 elseif ~ismember(method,0:5) 0129 error 'If supplied, method must be one of: {0,1,2,3,4,5}.' 0130 end 0131 0132 % for different methods 0133 switch method 0134 case 0 0135 % The same as method == 1, except only work on those 0136 % elements which are NaN, or at least touch a NaN. 0137 0138 % is it 1-d or 2-d? 0139 if (m == 1) || (n == 1) 0140 % really a 1-d case 0141 work_list = nan_list(:,1); 0142 work_list = unique([work_list;work_list - 1;work_list + 1]); 0143 work_list(work_list <= 1) = []; 0144 work_list(work_list >= nm) = []; 0145 nw = numel(work_list); 0146 0147 u = (1:nw)'; 0148 fda = sparse(repmat(u,1,3),bsxfun(@plus,work_list,-1:1), ... 0149 repmat([1 -2 1],nw,1),nw,nm); 0150 else 0151 % a 2-d case 0152 0153 % horizontal and vertical neighbors only 0154 talks_to = [-1 0;0 -1;1 0;0 1]; 0155 neighbors_list=identify_neighbors(n,m,nan_list,talks_to); 0156 0157 % list of all nodes we have identified 0158 all_list=[nan_list;neighbors_list]; 0159 0160 % generate sparse array with second partials on row 0161 % variable for each element in either list, but only 0162 % for those nodes which have a row index > 1 or < n 0163 L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); 0164 nl=length(L); 0165 if nl>0 0166 fda=sparse(repmat(all_list(L,1),1,3), ... 0167 repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... 0168 repmat([1 -2 1],nl,1),nm,nm); 0169 else 0170 fda=spalloc(n*m,n*m,size(all_list,1)*5); 0171 end 0172 0173 % 2nd partials on column index 0174 L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); 0175 nl=length(L); 0176 if nl>0 0177 fda=fda+sparse(repmat(all_list(L,1),1,3), ... 0178 repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... 0179 repmat([1 -2 1],nl,1),nm,nm); 0180 end 0181 end 0182 0183 % eliminate knowns 0184 rhs=-fda(:,known_list)*A(known_list); 0185 k=find(any(fda(:,nan_list(:,1)),2)); 0186 0187 % and solve... 0188 B=A; 0189 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); 0190 0191 case 1 0192 % least squares approach with del^2. Build system 0193 % for every array element as an unknown, and then 0194 % eliminate those which are knowns. 0195 0196 % Build sparse matrix approximating del^2 for 0197 % every element in A. 0198 0199 % is it 1-d or 2-d? 0200 if (m == 1) || (n == 1) 0201 % a 1-d case 0202 u = (1:(nm-2))'; 0203 fda = sparse(repmat(u,1,3),bsxfun(@plus,u,0:2), ... 0204 repmat([1 -2 1],nm-2,1),nm-2,nm); 0205 else 0206 % a 2-d case 0207 0208 % Compute finite difference for second partials 0209 % on row variable first 0210 [i,j]=ndgrid(2:(n-1),1:m); 0211 ind=i(:)+(j(:)-1)*n; 0212 np=(n-2)*m; 0213 fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... 0214 repmat([1 -2 1],np,1),n*m,n*m); 0215 0216 % now second partials on column variable 0217 [i,j]=ndgrid(1:n,2:(m-1)); 0218 ind=i(:)+(j(:)-1)*n; 0219 np=n*(m-2); 0220 fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ... 0221 repmat([1 -2 1],np,1),nm,nm); 0222 end 0223 0224 % eliminate knowns 0225 rhs=-fda(:,known_list)*A(known_list); 0226 k=find(any(fda(:,nan_list),2)); 0227 0228 % and solve... 0229 B=A; 0230 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); 0231 0232 case 2 0233 % Direct solve for del^2 BVP across holes 0234 0235 % generate sparse array with second partials on row 0236 % variable for each nan element, only for those nodes 0237 % which have a row index > 1 or < n 0238 0239 % is it 1-d or 2-d? 0240 if (m == 1) || (n == 1) 0241 % really just a 1-d case 0242 error('Method 2 has problems for vector input. Please use another method.') 0243 0244 else 0245 % a 2-d case 0246 L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n)); 0247 nl=length(L); 0248 if nl>0 0249 fda=sparse(repmat(nan_list(L,1),1,3), ... 0250 repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... 0251 repmat([1 -2 1],nl,1),n*m,n*m); 0252 else 0253 fda=spalloc(n*m,n*m,size(nan_list,1)*5); 0254 end 0255 0256 % 2nd partials on column index 0257 L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m)); 0258 nl=length(L); 0259 if nl>0 0260 fda=fda+sparse(repmat(nan_list(L,1),1,3), ... 0261 repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... 0262 repmat([1 -2 1],nl,1),n*m,n*m); 0263 end 0264 0265 % fix boundary conditions at extreme corners 0266 % of the array in case there were nans there 0267 if ismember(1,nan_list(:,1)) 0268 fda(1,[1 2 n+1])=[-2 1 1]; 0269 end 0270 if ismember(n,nan_list(:,1)) 0271 fda(n,[n, n-1,n+n])=[-2 1 1]; 0272 end 0273 if ismember(nm-n+1,nan_list(:,1)) 0274 fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1]; 0275 end 0276 if ismember(nm,nan_list(:,1)) 0277 fda(nm,[nm,nm-1,nm-n])=[-2 1 1]; 0278 end 0279 0280 % eliminate knowns 0281 rhs=-fda(:,known_list)*A(known_list); 0282 0283 % and solve... 0284 B=A; 0285 k=nan_list(:,1); 0286 B(k)=fda(k,k)\rhs(k); 0287 0288 end 0289 0290 case 3 0291 % The same as method == 0, except uses del^4 as the 0292 % interpolating operator. 0293 0294 % del^4 template of neighbors 0295 talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ... 0296 0 1;0 2;1 -1;1 0;1 1;2 0]; 0297 neighbors_list=identify_neighbors(n,m,nan_list,talks_to); 0298 0299 % list of all nodes we have identified 0300 all_list=[nan_list;neighbors_list]; 0301 0302 % generate sparse array with del^4, but only 0303 % for those nodes which have a row & column index 0304 % >= 3 or <= n-2 0305 L = find( (all_list(:,2) >= 3) & ... 0306 (all_list(:,2) <= (n-2)) & ... 0307 (all_list(:,3) >= 3) & ... 0308 (all_list(:,3) <= (m-2))); 0309 nl=length(L); 0310 if nl>0 0311 % do the entire template at once 0312 fda=sparse(repmat(all_list(L,1),1,13), ... 0313 repmat(all_list(L,1),1,13) + ... 0314 repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ... 0315 repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm); 0316 else 0317 fda=spalloc(n*m,n*m,size(all_list,1)*5); 0318 end 0319 0320 % on the boundaries, reduce the order around the edges 0321 L = find((((all_list(:,2) == 2) | ... 0322 (all_list(:,2) == (n-1))) & ... 0323 (all_list(:,3) >= 2) & ... 0324 (all_list(:,3) <= (m-1))) | ... 0325 (((all_list(:,3) == 2) | ... 0326 (all_list(:,3) == (m-1))) & ... 0327 (all_list(:,2) >= 2) & ... 0328 (all_list(:,2) <= (n-1)))); 0329 nl=length(L); 0330 if nl>0 0331 fda=fda+sparse(repmat(all_list(L,1),1,5), ... 0332 repmat(all_list(L,1),1,5) + ... 0333 repmat([-n,-1,0,+1,n],nl,1), ... 0334 repmat([1 1 -4 1 1],nl,1),nm,nm); 0335 end 0336 0337 L = find( ((all_list(:,2) == 1) | ... 0338 (all_list(:,2) == n)) & ... 0339 (all_list(:,3) >= 2) & ... 0340 (all_list(:,3) <= (m-1))); 0341 nl=length(L); 0342 if nl>0 0343 fda=fda+sparse(repmat(all_list(L,1),1,3), ... 0344 repmat(all_list(L,1),1,3) + ... 0345 repmat([-n,0,n],nl,1), ... 0346 repmat([1 -2 1],nl,1),nm,nm); 0347 end 0348 0349 L = find( ((all_list(:,3) == 1) | ... 0350 (all_list(:,3) == m)) & ... 0351 (all_list(:,2) >= 2) & ... 0352 (all_list(:,2) <= (n-1))); 0353 nl=length(L); 0354 if nl>0 0355 fda=fda+sparse(repmat(all_list(L,1),1,3), ... 0356 repmat(all_list(L,1),1,3) + ... 0357 repmat([-1,0,1],nl,1), ... 0358 repmat([1 -2 1],nl,1),nm,nm); 0359 end 0360 0361 % eliminate knowns 0362 rhs=-fda(:,known_list)*A(known_list); 0363 k=find(any(fda(:,nan_list(:,1)),2)); 0364 0365 % and solve... 0366 B=A; 0367 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); 0368 0369 case 4 0370 % Spring analogy 0371 % interpolating operator. 0372 0373 % list of all springs between a node and a horizontal 0374 % or vertical neighbor 0375 hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1]; 0376 hv_springs=[]; 0377 for i=1:4 0378 hvs=nan_list+repmat(hv_list(i,:),nan_count,1); 0379 k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m); 0380 hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]]; 0381 end 0382 0383 % delete replicate springs 0384 hv_springs=unique(sort(hv_springs,2),'rows'); 0385 0386 % build sparse matrix of connections, springs 0387 % connecting diagonal neighbors are weaker than 0388 % the horizontal and vertical springs 0389 nhv=size(hv_springs,1); 0390 springs=sparse(repmat((1:nhv)',1,2),hv_springs, ... 0391 repmat([1 -1],nhv,1),nhv,nm); 0392 0393 % eliminate knowns 0394 rhs=-springs(:,known_list)*A(known_list); 0395 0396 % and solve... 0397 B=A; 0398 B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs; 0399 0400 case 5 0401 % Average of 8 nearest neighbors 0402 0403 % generate sparse array to average 8 nearest neighbors 0404 % for each nan element, be careful around edges 0405 fda=spalloc(n*m,n*m,size(nan_list,1)*9); 0406 0407 % -1,-1 0408 L = find((nan_list(:,2) > 1) & (nan_list(:,3) > 1)); 0409 nl=length(L); 0410 if nl>0 0411 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0412 repmat(nan_list(L,1),1,2)+repmat([-n-1, 0],nl,1), ... 0413 repmat([1 -1],nl,1),n*m,n*m); 0414 end 0415 0416 % 0,-1 0417 L = find(nan_list(:,3) > 1); 0418 nl=length(L); 0419 if nl>0 0420 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0421 repmat(nan_list(L,1),1,2)+repmat([-n, 0],nl,1), ... 0422 repmat([1 -1],nl,1),n*m,n*m); 0423 end 0424 0425 % +1,-1 0426 L = find((nan_list(:,2) < n) & (nan_list(:,3) > 1)); 0427 nl=length(L); 0428 if nl>0 0429 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0430 repmat(nan_list(L,1),1,2)+repmat([-n+1, 0],nl,1), ... 0431 repmat([1 -1],nl,1),n*m,n*m); 0432 end 0433 0434 % -1,0 0435 L = find(nan_list(:,2) > 1); 0436 nl=length(L); 0437 if nl>0 0438 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0439 repmat(nan_list(L,1),1,2)+repmat([-1, 0],nl,1), ... 0440 repmat([1 -1],nl,1),n*m,n*m); 0441 end 0442 0443 % +1,0 0444 L = find(nan_list(:,2) < n); 0445 nl=length(L); 0446 if nl>0 0447 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0448 repmat(nan_list(L,1),1,2)+repmat([1, 0],nl,1), ... 0449 repmat([1 -1],nl,1),n*m,n*m); 0450 end 0451 0452 % -1,+1 0453 L = find((nan_list(:,2) > 1) & (nan_list(:,3) < m)); 0454 nl=length(L); 0455 if nl>0 0456 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0457 repmat(nan_list(L,1),1,2)+repmat([n-1, 0],nl,1), ... 0458 repmat([1 -1],nl,1),n*m,n*m); 0459 end 0460 0461 % 0,+1 0462 L = find(nan_list(:,3) < m); 0463 nl=length(L); 0464 if nl>0 0465 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0466 repmat(nan_list(L,1),1,2)+repmat([n, 0],nl,1), ... 0467 repmat([1 -1],nl,1),n*m,n*m); 0468 end 0469 0470 % +1,+1 0471 L = find((nan_list(:,2) < n) & (nan_list(:,3) < m)); 0472 nl=length(L); 0473 if nl>0 0474 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... 0475 repmat(nan_list(L,1),1,2)+repmat([n+1, 0],nl,1), ... 0476 repmat([1 -1],nl,1),n*m,n*m); 0477 end 0478 0479 % eliminate knowns 0480 rhs=-fda(:,known_list)*A(known_list); 0481 0482 % and solve... 0483 B=A; 0484 k=nan_list(:,1); 0485 B(k)=fda(k,k)\rhs(k); 0486 0487 end 0488 0489 % all done, make sure that B is the same shape as 0490 % A was when we came in. 0491 B=reshape(B,n,m); 0492 0493 0494 % ==================================================== 0495 % end of main function 0496 % ==================================================== 0497 % ==================================================== 0498 % begin subfunctions 0499 % ==================================================== 0500 function neighbors_list=identify_neighbors(n,m,nan_list,talks_to) 0501 % identify_neighbors: identifies all the neighbors of 0502 % those nodes in nan_list, not including the nans 0503 % themselves 0504 % 0505 % arguments (input): 0506 % n,m - scalar - [n,m]=size(A), where A is the 0507 % array to be interpolated 0508 % nan_list - array - list of every nan element in A 0509 % nan_list(i,1) == linear index of i'th nan element 0510 % nan_list(i,2) == row index of i'th nan element 0511 % nan_list(i,3) == column index of i'th nan element 0512 % talks_to - px2 array - defines which nodes communicate 0513 % with each other, i.e., which nodes are neighbors. 0514 % 0515 % talks_to(i,1) - defines the offset in the row 0516 % dimension of a neighbor 0517 % talks_to(i,2) - defines the offset in the column 0518 % dimension of a neighbor 0519 % 0520 % For example, talks_to = [-1 0;0 -1;1 0;0 1] 0521 % means that each node talks only to its immediate 0522 % neighbors horizontally and vertically. 0523 % 0524 % arguments(output): 0525 % neighbors_list - array - list of all neighbors of 0526 % all the nodes in nan_list 0527 0528 if ~isempty(nan_list) 0529 % use the definition of a neighbor in talks_to 0530 nan_count=size(nan_list,1); 0531 talk_count=size(talks_to,1); 0532 0533 nn=zeros(nan_count*talk_count,2); 0534 j=[1,nan_count]; 0535 for i=1:talk_count 0536 nn(j(1):j(2),:)=nan_list(:,2:3) + ... 0537 repmat(talks_to(i,:),nan_count,1); 0538 j=j+nan_count; 0539 end 0540 0541 % drop those nodes which fall outside the bounds of the 0542 % original array 0543 L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m); 0544 nn(L,:)=[]; 0545 0546 % form the same format 3 column array as nan_list 0547 neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn]; 0548 0549 % delete replicates in the neighbors list 0550 neighbors_list=unique(neighbors_list,'rows'); 0551 0552 % and delete those which are also in the list of NaNs. 0553 neighbors_list=setdiff(neighbors_list,nan_list,'rows'); 0554 0555 else 0556 neighbors_list=[]; 0557 end 0558 0559 0560 0561 0562 0563 0564 0565 0566 0567