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

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

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


Revision 1.2 - (hide annotations) (download)
Sat Feb 17 23:49:43 2007 UTC (14 years, 5 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 baylor 1.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 jmc 1.2 % $Header: $
69     % $Name: $
70    
71 baylor 1.1 % 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