/[MITgcm]/MITgcm/utils/matlab/inpaint_nans.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/inpaint_nans.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Wed Oct 19 18:30:02 2005 UTC (15 years, 9 months ago) by baylor
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58e_post, checkpoint58u_post, checkpoint58r_post, checkpoint57y_post, checkpoint58n_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint58j_post, checkpoint58, checkpoint58f_post, checkpoint57x_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint58a_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint58b_post, checkpoint58m_post
Adapt interpickups to be much more reliable regarding topography.  Now uses the inpaint_nans.m interpolation (which is now added to repository with the permission of the author) to fill through topography.  Then, new topography is used to zero the fields.  This should allow seamounts or walls to appear and dissapear between pickup files without unexpected consequences.

1 function B=inpaint_nans(A,method)
2 % inpaint_nans: in-paints over nans in an array
3 % usage: B=inpaint_nans(A)
4 %
5 % solves approximation to one of several pdes to
6 % interpolate and extrapolate holes
7 %
8 % Written by John D'Errico (woodchips@rochester.rr.com)
9 %
10 % arguments (input):
11 % A - nxm array with some NaNs to be filled in
12 %
13 % method - (OPTIONAL) scalar numeric flag - specifies
14 % which approach (or physical metaphor to use
15 % for the interpolation.) All methods are capable
16 % of extrapolation, some are better than others.
17 % There are also speed differences, as well as
18 % accuracy differences for smooth surfaces.
19 %
20 % methods {0,1,2} use a simple plate metaphor
21 % methods {3} use a better plate equation,
22 % but will be slower
23 % methods 4 use a spring metaphor
24 %
25 % method == 0 --> (DEFAULT) see method 1, but
26 % this method does not build as large of a
27 % linear system in the case of only a few
28 % NaNs in a large array.
29 % Extrapolation behavior is linear.
30 %
31 % method == 1 --> simple approach, applies del^2
32 % over the entire array, then drops those parts
33 % of the array which do not have any contact with
34 % NaNs. Uses a least squares approach, but it
35 % does not touch existing points.
36 % In the case of small arrays, this method is
37 % quite fast as it does very little extra work.
38 % Extrapolation behavior is linear.
39 %
40 % method == 2 --> uses del^2, but solving a direct
41 % linear system of equations for nan elements.
42 % This method will be the fastest possible for
43 % large systems since it uses the sparsest
44 % possible system of equations. Not a least
45 % squares approach, so it may be least robust
46 % to noise on the boundaries of any holes.
47 % This method will also be least able to
48 % interpolate accurately for smooth surfaces.
49 % Extrapolation behavior is linear.
50 %
51 % method == 3 --+ See method 0, but uses del^4 for
52 % the interpolating operator. This may result
53 % in more accurate interpolations, at some cost
54 % in speed.
55 %
56 % method == 4 --+ Uses a spring metaphor. Assumes
57 % springs (with a nominal length of zero)
58 % connect each node with every neighbor
59 % (horizontally, vertically and diagonally)
60 % Since each node tries to be like its neighbors,
61 % extrapolation is as a constant function where
62 % this is consistent with the neighboring nodes.
63 %
64 %
65 % arguments (output):
66 % B - nxm array with NaNs replaced
67
68 % I always need to know which elements are NaN,
69 % and what size the array is for any method
70 [n,m]=size(A);
71 nm=n*m;
72 k=isnan(A(:));
73
74 % list the nodes which are known, and which will
75 % be interpolated
76 nan_list=find(k);
77 known_list=find(~k);
78
79 % how many nans overall
80 nan_count=length(nan_list);
81
82 % convert NaN indices to (r,c) form
83 % nan_list==find(k) are the unrolled (linear) indices
84 % (row,column) form
85 [nr,nc]=ind2sub([n,m],nan_list);
86
87 % both forms of index in one array:
88 % column 1 == unrolled index
89 % column 2 == row index
90 % column 3 == column index
91 nan_list=[nan_list,nr,nc];
92
93 % supply default method
94 if (nargin<2)|isempty(method)
95 method = 0;
96 elseif ~ismember(method,0:4)
97 error 'If supplied, method must be one of: {0,1,2,3,4}.'
98 end
99
100 % for different methods
101 switch method
102 case 0
103 % The same as method == 1, except only work on those
104 % elements which are NaN, or at least touch a NaN.
105
106 % horizontal and vertical neighbors only
107 talks_to = [-1 0;0 -1;1 0;0 1];
108 neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
109
110 % list of all nodes we have identified
111 all_list=[nan_list;neighbors_list];
112
113 % generate sparse array with second partials on row
114 % variable for each element in either list, but only
115 % for those nodes which have a row index > 1 or < n
116 L = find((all_list(:,2) > 1) & (all_list(:,2) < n));
117 nl=length(L);
118 if nl>0
119 fda=sparse(repmat(all_list(L,1),1,3), ...
120 repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
121 repmat([1 -2 1],nl,1),nm,nm);
122 else
123 fda=spalloc(n*m,n*m,size(all_list,1)*5);
124 end
125
126 % 2nd partials on column index
127 L = find((all_list(:,3) > 1) & (all_list(:,3) < m));
128 nl=length(L);
129 if nl>0
130 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
131 repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
132 repmat([1 -2 1],nl,1),nm,nm);
133 end
134
135 % eliminate knowns
136 rhs=-fda(:,known_list)*A(known_list);
137 k=find(any(fda(:,nan_list(:,1)),2));
138
139 % and solve...
140 B=A;
141 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
142
143 case 1
144 % least squares approach with del^2. Build system
145 % for every array element as an unknown, and then
146 % eliminate those which are knowns.
147
148 % Build sparse matrix approximating del^2 for
149 % every element in A.
150 % Compute finite difference for second partials
151 % on row variable first
152 [i,j]=ndgrid(2:(n-1),1:m);
153 ind=i(:)+(j(:)-1)*n;
154 np=(n-2)*m;
155 fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
156 repmat([1 -2 1],np,1),n*m,n*m);
157
158 % now second partials on column variable
159 [i,j]=ndgrid(1:n,2:(m-1));
160 ind=i(:)+(j(:)-1)*n;
161 np=n*(m-2);
162 fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
163 repmat([1 -2 1],np,1),nm,nm);
164
165 % eliminate knowns
166 rhs=-fda(:,known_list)*A(known_list);
167 k=find(any(fda(:,nan_list),2));
168
169 % and solve...
170 B=A;
171 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
172
173 case 2
174 % Direct solve for del^2 BVP across holes
175
176 % generate sparse array with second partials on row
177 % variable for each nan element, only for those nodes
178 % which have a row index > 1 or < n
179 L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n));
180 nl=length(L);
181 if nl>0
182 fda=sparse(repmat(nan_list(L,1),1,3), ...
183 repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
184 repmat([1 -2 1],nl,1),n*m,n*m);
185 else
186 fda=spalloc(n*m,n*m,size(nan_list,1)*5);
187 end
188
189 % 2nd partials on column index
190 L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m));
191 nl=length(L);
192 if nl>0
193 fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
194 repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
195 repmat([1 -2 1],nl,1),n*m,n*m);
196 end
197
198 % fix boundary conditions at extreme corners
199 % of the array in case there were nans there
200 if ismember(1,nan_list(:,1))
201 fda(1,[1 2 n+1])=[-2 1 1];
202 end
203 if ismember(n,nan_list(:,1))
204 fda(n,[n, n-1,n+n])=[-2 1 1];
205 end
206 if ismember(nm-n+1,nan_list(:,1))
207 fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
208 end
209 if ismember(nm,nan_list(:,1))
210 fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
211 end
212
213 % eliminate knowns
214 rhs=-fda(:,known_list)*A(known_list);
215
216 % and solve...
217 B=A;
218 k=nan_list(:,1);
219 B(k)=fda(k,k)\rhs(k);
220
221 case 3
222 % The same as method == 0, except uses del^4 as the
223 % interpolating operator.
224
225 % del^4 template of neighbors
226 talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
227 0 1;0 2;1 -1;1 0;1 1;2 0];
228 neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
229
230 % list of all nodes we have identified
231 all_list=[nan_list;neighbors_list];
232
233 % generate sparse array with del^4, but only
234 % for those nodes which have a row & column index
235 % >= 3 or <= n-2
236 L = find( (all_list(:,2) >= 3) & ...
237 (all_list(:,2) <= (n-2)) & ...
238 (all_list(:,3) >= 3) & ...
239 (all_list(:,3) <= (m-2)));
240 nl=length(L);
241 if nl>0
242 % do the entire template at once
243 fda=sparse(repmat(all_list(L,1),1,13), ...
244 repmat(all_list(L,1),1,13) + ...
245 repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ...
246 repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm);
247 else
248 fda=spalloc(n*m,n*m,size(all_list,1)*5);
249 end
250
251 % on the boundaries, reduce the order around the edges
252 L = find((((all_list(:,2) == 2) | ...
253 (all_list(:,2) == (n-1))) & ...
254 (all_list(:,3) >= 2) & ...
255 (all_list(:,3) <= (m-1))) | ...
256 (((all_list(:,3) == 2) | ...
257 (all_list(:,3) == (m-1))) & ...
258 (all_list(:,2) >= 2) & ...
259 (all_list(:,2) <= (n-1))));
260 nl=length(L);
261 if nl>0
262 fda=fda+sparse(repmat(all_list(L,1),1,5), ...
263 repmat(all_list(L,1),1,5) + ...
264 repmat([-n,-1,0,+1,n],nl,1), ...
265 repmat([1 1 -4 1 1],nl,1),nm,nm);
266 end
267
268 L = find( ((all_list(:,2) == 1) | ...
269 (all_list(:,2) == n)) & ...
270 (all_list(:,3) >= 2) & ...
271 (all_list(:,3) <= (m-1)));
272 nl=length(L);
273 if nl>0
274 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
275 repmat(all_list(L,1),1,3) + ...
276 repmat([-n,0,n],nl,1), ...
277 repmat([1 -2 1],nl,1),nm,nm);
278 end
279
280 L = find( ((all_list(:,3) == 1) | ...
281 (all_list(:,3) == m)) & ...
282 (all_list(:,2) >= 2) & ...
283 (all_list(:,2) <= (n-1)));
284 nl=length(L);
285 if nl>0
286 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
287 repmat(all_list(L,1),1,3) + ...
288 repmat([-1,0,1],nl,1), ...
289 repmat([1 -2 1],nl,1),nm,nm);
290 end
291
292 % eliminate knowns
293 rhs=-fda(:,known_list)*A(known_list);
294 k=find(any(fda(:,nan_list(:,1)),2));
295
296 % and solve...
297 B=A;
298 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
299
300 case 4
301 % Spring analogy
302 % interpolating operator.
303
304 % list of all springs between a node and a horizontal
305 % or vertical neighbor
306 hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1];
307 hv_springs=[];
308 for i=1:4
309 hvs=nan_list+repmat(hv_list(i,:),nan_count,1);
310 k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m);
311 hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]];
312 end
313
314 % delete replicate springs
315 hv_springs=unique(sort(hv_springs,2),'rows');
316
317 % build sparse matrix of connections, springs
318 % connecting diagonal neighbors are weaker than
319 % the horizontal and vertical springs
320 nhv=size(hv_springs,1);
321 springs=sparse(repmat((1:nhv)',1,2),hv_springs, ...
322 repmat([1 -1],nhv,1),nhv,nm);
323
324 % eliminate knowns
325 rhs=-springs(:,known_list)*A(known_list);
326
327 % and solve...
328 B=A;
329 B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs;
330
331 end
332
333 % ====================================================
334 % end of main function
335 % ====================================================
336 % ====================================================
337 % begin subfunctions
338 % ====================================================
339 function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)
340 % identify_neighbors: identifies all the neighbors of
341 % those nodes in nan_list, not including the nans
342 % themselves
343 %
344 % arguments (input):
345 % n,m - scalar - [n,m]=size(A), where A is the
346 % array to be interpolated
347 % nan_list - array - list of every nan element in A
348 % nan_list(i,1) == linear index of i'th nan element
349 % nan_list(i,2) == row index of i'th nan element
350 % nan_list(i,3) == column index of i'th nan element
351 % talks_to - px2 array - defines which nodes communicate
352 % with each other, i.e., which nodes are neighbors.
353 %
354 % talks_to(i,1) - defines the offset in the row
355 % dimension of a neighbor
356 % talks_to(i,2) - defines the offset in the column
357 % dimension of a neighbor
358 %
359 % For example, talks_to = [-1 0;0 -1;1 0;0 1]
360 % means that each node talks only to its immediate
361 % neighbors horizontally and vertically.
362 %
363 % arguments(output):
364 % neighbors_list - array - list of all neighbors of
365 % all the nodes in nan_list
366
367 if ~isempty(nan_list)
368 % use the definition of a neighbor in talks_to
369 nan_count=size(nan_list,1);
370 talk_count=size(talks_to,1);
371
372 nn=zeros(nan_count*talk_count,2);
373 j=[1,nan_count];
374 for i=1:talk_count
375 nn(j(1):j(2),:)=nan_list(:,2:3) + ...
376 repmat(talks_to(i,:),nan_count,1);
377 j=j+nan_count;
378 end
379
380 % drop those nodes which fall outside the bounds of the
381 % original array
382 L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m);
383 nn(L,:)=[];
384
385 % form the same format 3 column array as nan_list
386 neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn];
387
388 % delete replicates in the neighbors list
389 neighbors_list=unique(neighbors_list,'rows');
390
391 % and delete those which are also in the list of NaNs.
392 neighbors_list=setdiff(neighbors_list,nan_list,'rows');
393
394 else
395 neighbors_list=[];
396 end
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412

  ViewVC Help
Powered by ViewVC 1.1.22