/[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.2 - (show annotations) (download)
Sat Feb 17 23:49:43 2007 UTC (14 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint58w_post, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint58x_post, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint58y_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, HEAD
Changes since 1.1: +3 -0 lines
add $Header:  $ and $Name:  & (for CVS)

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 % $Header: $
69 % $Name: $
70
71 % I always need to know which elements are NaN,
72 % and what size the array is for any method
73 [n,m]=size(A);
74 nm=n*m;
75 k=isnan(A(:));
76
77 % list the nodes which are known, and which will
78 % be interpolated
79 nan_list=find(k);
80 known_list=find(~k);
81
82 % how many nans overall
83 nan_count=length(nan_list);
84
85 % convert NaN indices to (r,c) form
86 % nan_list==find(k) are the unrolled (linear) indices
87 % (row,column) form
88 [nr,nc]=ind2sub([n,m],nan_list);
89
90 % both forms of index in one array:
91 % column 1 == unrolled index
92 % column 2 == row index
93 % column 3 == column index
94 nan_list=[nan_list,nr,nc];
95
96 % supply default method
97 if (nargin<2)|isempty(method)
98 method = 0;
99 elseif ~ismember(method,0:4)
100 error 'If supplied, method must be one of: {0,1,2,3,4}.'
101 end
102
103 % for different methods
104 switch method
105 case 0
106 % The same as method == 1, except only work on those
107 % elements which are NaN, or at least touch a NaN.
108
109 % horizontal and vertical neighbors only
110 talks_to = [-1 0;0 -1;1 0;0 1];
111 neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
112
113 % list of all nodes we have identified
114 all_list=[nan_list;neighbors_list];
115
116 % generate sparse array with second partials on row
117 % variable for each element in either list, but only
118 % for those nodes which have a row index > 1 or < n
119 L = find((all_list(:,2) > 1) & (all_list(:,2) < n));
120 nl=length(L);
121 if nl>0
122 fda=sparse(repmat(all_list(L,1),1,3), ...
123 repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
124 repmat([1 -2 1],nl,1),nm,nm);
125 else
126 fda=spalloc(n*m,n*m,size(all_list,1)*5);
127 end
128
129 % 2nd partials on column index
130 L = find((all_list(:,3) > 1) & (all_list(:,3) < m));
131 nl=length(L);
132 if nl>0
133 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
134 repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
135 repmat([1 -2 1],nl,1),nm,nm);
136 end
137
138 % eliminate knowns
139 rhs=-fda(:,known_list)*A(known_list);
140 k=find(any(fda(:,nan_list(:,1)),2));
141
142 % and solve...
143 B=A;
144 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
145
146 case 1
147 % least squares approach with del^2. Build system
148 % for every array element as an unknown, and then
149 % eliminate those which are knowns.
150
151 % Build sparse matrix approximating del^2 for
152 % every element in A.
153 % Compute finite difference for second partials
154 % on row variable first
155 [i,j]=ndgrid(2:(n-1),1:m);
156 ind=i(:)+(j(:)-1)*n;
157 np=(n-2)*m;
158 fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
159 repmat([1 -2 1],np,1),n*m,n*m);
160
161 % now second partials on column variable
162 [i,j]=ndgrid(1:n,2:(m-1));
163 ind=i(:)+(j(:)-1)*n;
164 np=n*(m-2);
165 fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
166 repmat([1 -2 1],np,1),nm,nm);
167
168 % eliminate knowns
169 rhs=-fda(:,known_list)*A(known_list);
170 k=find(any(fda(:,nan_list),2));
171
172 % and solve...
173 B=A;
174 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
175
176 case 2
177 % Direct solve for del^2 BVP across holes
178
179 % generate sparse array with second partials on row
180 % variable for each nan element, only for those nodes
181 % which have a row index > 1 or < n
182 L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n));
183 nl=length(L);
184 if nl>0
185 fda=sparse(repmat(nan_list(L,1),1,3), ...
186 repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
187 repmat([1 -2 1],nl,1),n*m,n*m);
188 else
189 fda=spalloc(n*m,n*m,size(nan_list,1)*5);
190 end
191
192 % 2nd partials on column index
193 L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m));
194 nl=length(L);
195 if nl>0
196 fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
197 repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
198 repmat([1 -2 1],nl,1),n*m,n*m);
199 end
200
201 % fix boundary conditions at extreme corners
202 % of the array in case there were nans there
203 if ismember(1,nan_list(:,1))
204 fda(1,[1 2 n+1])=[-2 1 1];
205 end
206 if ismember(n,nan_list(:,1))
207 fda(n,[n, n-1,n+n])=[-2 1 1];
208 end
209 if ismember(nm-n+1,nan_list(:,1))
210 fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
211 end
212 if ismember(nm,nan_list(:,1))
213 fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
214 end
215
216 % eliminate knowns
217 rhs=-fda(:,known_list)*A(known_list);
218
219 % and solve...
220 B=A;
221 k=nan_list(:,1);
222 B(k)=fda(k,k)\rhs(k);
223
224 case 3
225 % The same as method == 0, except uses del^4 as the
226 % interpolating operator.
227
228 % del^4 template of neighbors
229 talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
230 0 1;0 2;1 -1;1 0;1 1;2 0];
231 neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
232
233 % list of all nodes we have identified
234 all_list=[nan_list;neighbors_list];
235
236 % generate sparse array with del^4, but only
237 % for those nodes which have a row & column index
238 % >= 3 or <= n-2
239 L = find( (all_list(:,2) >= 3) & ...
240 (all_list(:,2) <= (n-2)) & ...
241 (all_list(:,3) >= 3) & ...
242 (all_list(:,3) <= (m-2)));
243 nl=length(L);
244 if nl>0
245 % do the entire template at once
246 fda=sparse(repmat(all_list(L,1),1,13), ...
247 repmat(all_list(L,1),1,13) + ...
248 repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ...
249 repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm);
250 else
251 fda=spalloc(n*m,n*m,size(all_list,1)*5);
252 end
253
254 % on the boundaries, reduce the order around the edges
255 L = find((((all_list(:,2) == 2) | ...
256 (all_list(:,2) == (n-1))) & ...
257 (all_list(:,3) >= 2) & ...
258 (all_list(:,3) <= (m-1))) | ...
259 (((all_list(:,3) == 2) | ...
260 (all_list(:,3) == (m-1))) & ...
261 (all_list(:,2) >= 2) & ...
262 (all_list(:,2) <= (n-1))));
263 nl=length(L);
264 if nl>0
265 fda=fda+sparse(repmat(all_list(L,1),1,5), ...
266 repmat(all_list(L,1),1,5) + ...
267 repmat([-n,-1,0,+1,n],nl,1), ...
268 repmat([1 1 -4 1 1],nl,1),nm,nm);
269 end
270
271 L = find( ((all_list(:,2) == 1) | ...
272 (all_list(:,2) == n)) & ...
273 (all_list(:,3) >= 2) & ...
274 (all_list(:,3) <= (m-1)));
275 nl=length(L);
276 if nl>0
277 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
278 repmat(all_list(L,1),1,3) + ...
279 repmat([-n,0,n],nl,1), ...
280 repmat([1 -2 1],nl,1),nm,nm);
281 end
282
283 L = find( ((all_list(:,3) == 1) | ...
284 (all_list(:,3) == m)) & ...
285 (all_list(:,2) >= 2) & ...
286 (all_list(:,2) <= (n-1)));
287 nl=length(L);
288 if nl>0
289 fda=fda+sparse(repmat(all_list(L,1),1,3), ...
290 repmat(all_list(L,1),1,3) + ...
291 repmat([-1,0,1],nl,1), ...
292 repmat([1 -2 1],nl,1),nm,nm);
293 end
294
295 % eliminate knowns
296 rhs=-fda(:,known_list)*A(known_list);
297 k=find(any(fda(:,nan_list(:,1)),2));
298
299 % and solve...
300 B=A;
301 B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
302
303 case 4
304 % Spring analogy
305 % interpolating operator.
306
307 % list of all springs between a node and a horizontal
308 % or vertical neighbor
309 hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1];
310 hv_springs=[];
311 for i=1:4
312 hvs=nan_list+repmat(hv_list(i,:),nan_count,1);
313 k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m);
314 hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]];
315 end
316
317 % delete replicate springs
318 hv_springs=unique(sort(hv_springs,2),'rows');
319
320 % build sparse matrix of connections, springs
321 % connecting diagonal neighbors are weaker than
322 % the horizontal and vertical springs
323 nhv=size(hv_springs,1);
324 springs=sparse(repmat((1:nhv)',1,2),hv_springs, ...
325 repmat([1 -1],nhv,1),nhv,nm);
326
327 % eliminate knowns
328 rhs=-springs(:,known_list)*A(known_list);
329
330 % and solve...
331 B=A;
332 B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs;
333
334 end
335
336 % ====================================================
337 % end of main function
338 % ====================================================
339 % ====================================================
340 % begin subfunctions
341 % ====================================================
342 function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)
343 % identify_neighbors: identifies all the neighbors of
344 % those nodes in nan_list, not including the nans
345 % themselves
346 %
347 % arguments (input):
348 % n,m - scalar - [n,m]=size(A), where A is the
349 % array to be interpolated
350 % nan_list - array - list of every nan element in A
351 % nan_list(i,1) == linear index of i'th nan element
352 % nan_list(i,2) == row index of i'th nan element
353 % nan_list(i,3) == column index of i'th nan element
354 % talks_to - px2 array - defines which nodes communicate
355 % with each other, i.e., which nodes are neighbors.
356 %
357 % talks_to(i,1) - defines the offset in the row
358 % dimension of a neighbor
359 % talks_to(i,2) - defines the offset in the column
360 % dimension of a neighbor
361 %
362 % For example, talks_to = [-1 0;0 -1;1 0;0 1]
363 % means that each node talks only to its immediate
364 % neighbors horizontally and vertically.
365 %
366 % arguments(output):
367 % neighbors_list - array - list of all neighbors of
368 % all the nodes in nan_list
369
370 if ~isempty(nan_list)
371 % use the definition of a neighbor in talks_to
372 nan_count=size(nan_list,1);
373 talk_count=size(talks_to,1);
374
375 nn=zeros(nan_count*talk_count,2);
376 j=[1,nan_count];
377 for i=1:talk_count
378 nn(j(1):j(2),:)=nan_list(:,2:3) + ...
379 repmat(talks_to(i,:),nan_count,1);
380 j=j+nan_count;
381 end
382
383 % drop those nodes which fall outside the bounds of the
384 % original array
385 L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m);
386 nn(L,:)=[];
387
388 % form the same format 3 column array as nan_list
389 neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn];
390
391 % delete replicates in the neighbors list
392 neighbors_list=unique(neighbors_list,'rows');
393
394 % and delete those which are also in the list of NaNs.
395 neighbors_list=setdiff(neighbors_list,nan_list,'rows');
396
397 else
398 neighbors_list=[];
399 end
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415

  ViewVC Help
Powered by ViewVC 1.1.22