function B=inpaint_nans(A,method)
% 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.
%
% 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/06
% I always need to know which elements are NaN,
% and what size the array is for any method
[n,m]=size(A);
A=A(:);
nm=n*m;
k=isnan(A(:));
% list the nodes which are known, and which will
% be interpolated
nan_list=find(k);
known_list=find(~k);
% how many nans overall
nan_count=length(nan_list);
% convert NaN indices to (r,c) form
% nan_list==find(k) are the unrolled (linear) indices
% (row,column) form
[nr,nc]=ind2sub([n,m],nan_list);
% both forms of index in one array:
% column 1 == unrolled index
% column 2 == row index
% column 3 == column index
nan_list=[nan_list,nr,nc];
% supply default method
if (nargin<2) || isempty(method)
method = 0;
elseif ~ismember(method,0:5)
error 'If supplied, method must be one of: {0,1,2,3,4,5}.'
end
% for different methods
switch method
case 0
% The same as method == 1, except only work on those
% elements which are NaN, or at least touch a NaN.
% horizontal and vertical neighbors only
talks_to = [-1 0;0 -1;1 0;0 1];
neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
% list of all nodes we have identified
all_list=[nan_list;neighbors_list];
% generate sparse array with second partials on row
% variable for each element in either list, but only
% for those nodes which have a row index > 1 or < n
L = find((all_list(:,2) > 1) & (all_list(:,2) < n));
nl=length(L);
if nl>0
fda=sparse(repmat(all_list(L,1),1,3), ...
repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
repmat([1 -2 1],nl,1),nm,nm);
else
fda=spalloc(n*m,n*m,size(all_list,1)*5);
end
% 2nd partials on column index
L = find((all_list(:,3) > 1) & (all_list(:,3) < m));
nl=length(L);
if nl>0
fda=fda+sparse(repmat(all_list(L,1),1,3), ...
repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
repmat([1 -2 1],nl,1),nm,nm);
end
% eliminate knowns
rhs=-fda(:,known_list)*A(known_list);
k=find(any(fda(:,nan_list(:,1)),2));
% and solve...
B=A;
B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
case 1
% least squares approach with del^2. Build system
% for every array element as an unknown, and then
% eliminate those which are knowns.
% Build sparse matrix approximating del^2 for
% every element in A.
% Compute finite difference for second partials
% on row variable first
[i,j]=ndgrid(2:(n-1),1:m);
ind=i(:)+(j(:)-1)*n;
np=(n-2)*m;
fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
repmat([1 -2 1],np,1),n*m,n*m);
% now second partials on column variable
[i,j]=ndgrid(1:n,2:(m-1));
ind=i(:)+(j(:)-1)*n;
np=n*(m-2);
fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
repmat([1 -2 1],np,1),nm,nm);
% eliminate knowns
rhs=-fda(:,known_list)*A(known_list);
k=find(any(fda(:,nan_list),2));
% and solve...
B=A;
B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
case 2
% Direct solve for del^2 BVP across holes
% generate sparse array with second partials on row
% variable for each nan element, only for those nodes
% which have a row index > 1 or < n
L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n));
nl=length(L);
if nl>0
fda=sparse(repmat(nan_list(L,1),1,3), ...
repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
repmat([1 -2 1],nl,1),n*m,n*m);
else
fda=spalloc(n*m,n*m,size(nan_list,1)*5);
end
% 2nd partials on column index
L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m));
nl=length(L);
if nl>0
fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
repmat([1 -2 1],nl,1),n*m,n*m);
end
% fix boundary conditions at extreme corners
% of the array in case there were nans there
if ismember(1,nan_list(:,1))
fda(1,[1 2 n+1])=[-2 1 1];
end
if ismember(n,nan_list(:,1))
fda(n,[n, n-1,n+n])=[-2 1 1];
end
if ismember(nm-n+1,nan_list(:,1))
fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
end
if ismember(nm,nan_list(:,1))
fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
end
% eliminate knowns
rhs=-fda(:,known_list)*A(known_list);
% and solve...
B=A;
k=nan_list(:,1);
B(k)=fda(k,k)\rhs(k);
case 3
% The same as method == 0, except uses del^4 as the
% interpolating operator.
% del^4 template of neighbors
talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
0 1;0 2;1 -1;1 0;1 1;2 0];
neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
% list of all nodes we have identified
all_list=[nan_list;neighbors_list];
% generate sparse array with del^4, but only
% for thos
机载LiDAR点云滤波-SMRF简单形态学滤波(MATLAB代码).zip
版权申诉
110 浏览量
2024-05-18
14:28:18
上传
评论
收藏 19KB ZIP 举报
![avatar](https://profile-avatar.csdnimg.cn/864ffdc5a26342a6add0026479aef1e5_matlab_dingdang.jpg!1)
![avatar-vip](https://csdnimg.cn/release/downloadcmsfe/public/img/user-vip.1c89f3c5.png)
matlab科研助手
- 粉丝: 2w+
- 资源: 2506
最新资源
- 12位双通道高速ADC芯片AD9238评估板开发模块ALTIUM设计硬件(原理图+PCB)工程文件.zip
- 基于74LS160的20到70的置数仿真节线图
- 解决mac上qt链接mysql方案加使用软件,经测试连接成功
- 锅炉引风机控制;变频调速技术;PLC;组态软件
- Java项目-电影院售票管理系统(java+Servlet+JSP+JDBC+Mysql)
- SSM整合开发-图书管理系统
- 计二202301020210蒋怡.zip
- 基于74LS160的70进制计数器仿真节线图
- 基于74LS160的30进制与70进制转化仿真节线图
- 【用360解压工具解压】springboot+vue实验室(预约)管理系统【www.java1234.com】.zip
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
![feedback](https://img-home.csdnimg.cn/images/20220527035711.png)
![feedback](https://img-home.csdnimg.cn/images/20220527035711.png)
![feedback-tip](https://img-home.csdnimg.cn/images/20220527035111.png)