Prev: restoration of signal
Next: writing liklihood function for estimation parameters in time series models
From: Ben on 13 Jul 2010 09:10
I have a medium to large set of 2D experimental data, which is spaced non-uniformly, and I need to interpolate it to a uniform grid. The size of each 2D field is usually around 2000x30 data points, and I can have as many as 5000 of these 2D fields. The other catch is every field is missing about 30% of it's data. So far, the missing data has always been at the edges, leading to concave and convex regions, but there is no reason why I might not get a large hole in the middle of my data.

At this point, I have a working solution for my data, but it's quite slow, and I don't think it will work well with a large hole in the data. Here is a stripped down version of my code that only has the essential parts in it:

%Cycle thru the fields
for f=1:5000
%Isolate data needed for this for loop
x_in=x_data(:,:,f);
y_in=y_data(:,:,f);
lgc_in=lgc_data(:,:,f);
%Create a Delaunay Triangulation only using the valid data points (lgc_in is a logical that specifies whether a data point is valid or missing)
x=x_in(lgc_in);
y=y_in(lgc_in);
dt = DelaunayTri(x,y);

%Find the circumcenters of the triangles
[cc,cr]=circumcenters(dt); %#ok<ASGLU>
%Find the triangles whose circumcenters are sufficiently small
%(large circumcenter radii indicate large missing areas of data,
%which should not be interpolated over)
tri_ndx=cr<5;
%Create a triangle representation with only the triangles that
%are sufficiently small
tri_edges=dt.Triangulation(tri_ndx,:);
tri=TriRep(tri_edges,x,y);
%Find the free edges of this triangle representation
fe=freeBoundary(tri);
%Add the free edges as a constraint to the Delauney
%Triangulation. (The TriScatteredInterp function below
%cannot use a generic triangle representation. It requires a
%Delauney Triangulation.)
%(A constraint can sometimes cut thru a triangle requiring the
%Delauney triangulation to be reconstructed, but here only
%entire triangles are eliminated by the constraint.)
dt.Constraints=fe;

%Find the limits of the field
AOI_lmt=round([min(x),max(x);min(y),max(y)]);

%Create interpolation list
[xi,yi]=meshgrid((AOI_lmt(1,1):AOI_lmt(1,2))',(AOI_lmt(2,1):AOI_lmt(2,2))');
[i,j]=size(xi);
n_pts=numel(xi);
xi=reshape(xi,n_pts,1);
yi=reshape(yi,n_pts,1);

%Find the triangle index for each xi-yi data point
tri_ndx=pointLocation(dt,xi,yi);
%Find the data points that are within the triangulation region
in_tri=~isnan(tri_ndx);
%Get the list of which triangles are inside constraint boundary
io=dt.inOutStatus();
%Of the points that are inside the triangulation region, some are inside
%triangles that are outside the constraint boundary. Find the points that
%are within the constraint boundary.
in_bndry=in_tri;
in_bndry(in_tri)=io(tri_ndx(in_tri));
%Set the points that are outside the constraint boundary to NaN. (Don't
%remove them since that would produce holes in the grid points, and we
%need a full grid to make a contour plot or image of the field data.)
xi(~in_bndry)=NaN;
yi(~in_bndry)=NaN;

%Find the valid field points
field_in=field_data(:,:,f);
z=field_in(lgc_in);
%Create the interpolation object
F=TriScatteredInterp(dt,z);
%Interpolate at the grid points (regions outside the
%constraint boundary have NaN, so interpolation will result
%in NaN at those data points).
zi=F(xi,yi);
end

I've run profiler on this code and the slowest parts are pointLocation(dt,xi,yi) and TriScatteredInterp(dt,z), which take about 40 seconds each. I guess I have to accept the computation time required for TriScatteredInterp, but I'm hoping that there may be a way to speed up the code associated with pointLocation. As you can see, I'm using a Delauney constraint boundary to specify where I want to interpolate (inside the boundary there is enough data so interpolate; outside the boundary there isn't enough data so do not interpolate). Then I use pointLocation to help me determine whether data points on the interpolation grid are inside or outside the Delauney constraint boundary.

Does anyone know a better (faster) way to do this? Furthermore I'm not sure if the freeBoundary(tri) function call will find a large hole in the data.

Thanks for any help you can lend.
From: us on 13 Jul 2010 09:48
"Ben " <breedlun(a)hotmail.com> wrote in message <i1hoju$p2$1(a)fred.mathworks.com>...
> I have a medium to large set of 2D experimental data, which is spaced non-uniformly, and I need to interpolate it to a uniform grid. The size of each 2D field is usually around 2000x30 data points, and I can have as many as 5000 of these 2D fields. The other catch is every field is missing about 30% of it's data. So far, the missing data has always been at the edges, leading to concave and convex regions, but there is no reason why I might not get a large hole in the middle of my data.
>
> At this point, I have a working solution for my data, but it's quite slow, and I don't think it will work well with a large hole in the data. Here is a stripped down version of my code that only has the essential parts in it:
>
> %Cycle thru the fields
> for f=1:5000
> %Isolate data needed for this for loop
> x_in=x_data(:,:,f);
> y_in=y_data(:,:,f);
> lgc_in=lgc_data(:,:,f);
> %Create a Delaunay Triangulation only using the valid data points (lgc_in is a logical that specifies whether a data point is valid or missing)
> x=x_in(lgc_in);
> y=y_in(lgc_in);
> dt = DelaunayTri(x,y);
>
> %Find the circumcenters of the triangles
> [cc,cr]=circumcenters(dt); %#ok<ASGLU>
> %Find the triangles whose circumcenters are sufficiently small
> %(large circumcenter radii indicate large missing areas of data,
> %which should not be interpolated over)
> tri_ndx=cr<5;
> %Create a triangle representation with only the triangles that
> %are sufficiently small
> tri_edges=dt.Triangulation(tri_ndx,:);
> tri=TriRep(tri_edges,x,y);
> %Find the free edges of this triangle representation
> fe=freeBoundary(tri);
> %Add the free edges as a constraint to the Delauney
> %Triangulation. (The TriScatteredInterp function below
> %cannot use a generic triangle representation. It requires a
> %Delauney Triangulation.)
> %(A constraint can sometimes cut thru a triangle requiring the
> %Delauney triangulation to be reconstructed, but here only
> %entire triangles are eliminated by the constraint.)
> dt.Constraints=fe;
>
> %Find the limits of the field
> AOI_lmt=round([min(x),max(x);min(y),max(y)]);
>
> %Create interpolation list
> [xi,yi]=meshgrid((AOI_lmt(1,1):AOI_lmt(1,2))',(AOI_lmt(2,1):AOI_lmt(2,2))');
> [i,j]=size(xi);
> n_pts=numel(xi);
> xi=reshape(xi,n_pts,1);
> yi=reshape(yi,n_pts,1);
>
> %Find the triangle index for each xi-yi data point
> tri_ndx=pointLocation(dt,xi,yi);
> %Find the data points that are within the triangulation region
> in_tri=~isnan(tri_ndx);
> %Get the list of which triangles are inside constraint boundary
> io=dt.inOutStatus();
> %Of the points that are inside the triangulation region, some are inside
> %triangles that are outside the constraint boundary. Find the points that
> %are within the constraint boundary.
> in_bndry=in_tri;
> in_bndry(in_tri)=io(tri_ndx(in_tri));
> %Set the points that are outside the constraint boundary to NaN. (Don't
> %remove them since that would produce holes in the grid points, and we
> %need a full grid to make a contour plot or image of the field data.)
> xi(~in_bndry)=NaN;
> yi(~in_bndry)=NaN;
>
> %Find the valid field points
> field_in=field_data(:,:,f);
> z=field_in(lgc_in);
> %Create the interpolation object
> F=TriScatteredInterp(dt,z);
> %Interpolate at the grid points (regions outside the
> %constraint boundary have NaN, so interpolation will result
> %in NaN at those data points).
> zi=F(xi,yi);
> end
>
> I've run profiler on this code and the slowest parts are pointLocation(dt,xi,yi) and TriScatteredInterp(dt,z), which take about 40 seconds each. I guess I have to accept the computation time required for TriScatteredInterp, but I'm hoping that there may be a way to speed up the code associated with pointLocation. As you can see, I'm using a Delauney constraint boundary to specify where I want to interpolate (inside the boundary there is enough data so interpolate; outside the boundary there isn't enough data so do not interpolate). Then I use pointLocation to help me determine whether data points on the interpolation grid are inside or outside the Delauney constraint boundary.
>
> Does anyone know a better (faster) way to do this? Furthermore I'm not sure if the freeBoundary(tri) function call will find a large hole in the data.
>
> Thanks for any help you can lend.

a hint:
- look at these two great FEX contributions by john d'errico

http://www.mathworks.com/matlabcentral/fileexchange/4551-inpaintnans

http://www.mathworks.com/matlabcentral/fileexchange/8998-surface-fitting-using-gridfit

us
From: Ben on 13 Jul 2010 10:05
> a hint:
> - look at these two great FEX contributions by john d'errico
>
> http://www.mathworks.com/matlabcentral/fileexchange/4551-inpaintnans
>
> http://www.mathworks.com/matlabcentral/fileexchange/8998-surface-fitting-using-gridfit
>
> us

Thanks for the reply, but I believe inpaintnans is for interpolating to fill in holes in your data and gridfit is for smoothing/extropolating your data. I do not want to interpolate/extrapolate into the areas where I have missing data. I only want to interpolate if I have enough data points nearby. Sorry if I was not clear above.
 | 
Pages: 1
Prev: restoration of signal
Next: writing liklihood function for estimation parameters in time series models