Thread Subject: Contour spline

Is there anyway to get a contour plot to use spline interpolation rather than linear?

"Travis" <sinusoid2@hotmail.com> wrote in message <gnngj6$n8f$1@fred.mathworks.com>...
> Is there anyway to get a contour plot to use spline interpolation rather than linear?

Not simply, I think not.

Even if you took the contours and then fit a
spline through them, splines here could do
some nasty things, crossing in places where
they should never do so.

John

"Travis" <sinusoid2@hotmail.com> wrote in message <gnngj6$n8f$1@fred.mathworks.com>...
> Is there anyway to get a contour plot to use spline interpolation rather than linear?

At the low level, a 2D plot in Matlab contains basic objects (line/patch) that composed by a a linear connexions between two points, no matter what.

If the contour "looks" too coarse, an alternative is may be approximate you function by spline in a finer grid, then call the contour on the interpolated spline function.

Bruno

An illustration:

img=peaks(16);

xi=1:0.125:size(img,2);
yi=1:0.125:size(img,1);
[XI YI]=meshgrid(xi,yi);
imgs=interp2(img,XI,YI,'spline');

clf
subplot(1,2,1);
contour(img,(-7:7));
subplot(1,2,2);
contour(imgs,(-7:7));

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gno9ci$1k$1@fred.mathworks.com>...
> An illustration:
>
> img=peaks(16);
>
> xi=1:0.125:size(img,2);
> yi=1:0.125:size(img,1);
> [XI YI]=meshgrid(xi,yi);
> imgs=interp2(img,XI,YI,'spline');
>
> clf
> subplot(1,2,1);
> contour(img,(-7:7));
> subplot(1,2,2);
> contour(imgs,(-7:7));
>
> % Bruno

This works well, but my data has NaN's in the bottom corners often, is there anyway to get it to ignore those?

"Travis" <sinusoid2@hotmail.com> wrote in message <gnpcmt$kad$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gno9ci$1k$1@fred.mathworks.com>...
> > An illustration:
> >
> > img=peaks(16);
> >
> > xi=1:0.125:size(img,2);
> > yi=1:0.125:size(img,1);
> > [XI YI]=meshgrid(xi,yi);
> > imgs=interp2(img,XI,YI,'spline');
> >
> > clf
> > subplot(1,2,1);
> > contour(img,(-7:7));
> > subplot(1,2,2);
> > contour(imgs,(-7:7));
> >
> > % Bruno
>
> This works well, but my data has NaN's in the bottom corners often, is there anyway to get it to ignore those?

Anyone? I have tried putting 0's in there, but that of course doesn't work. It is more representative than the NaN's, but it still shows data where none was collected. Is there a way to get MATLab to ignore those areas?

"Travis" <sinusoid2@hotmail.com> wrote in message <gnvn35$fd1$1@fred.mathworks.com>...
> "Travis" <sinusoid2@hotmail.com> wrote in message <gnpcmt$kad$1@fred.mathworks.com>...

>
> Anyone? I have tried putting 0's in there, but that of course doesn't work. It is more representative than the NaN's, but it still shows data where none was collected. Is there a way to get MATLab to ignore those areas?

I have no idea what you meant by "ignore NaN". Please elaborate.

When interpolate using spline method, the destination array get fully filled, even if there are few NaNs in the source array (Those NaN points are discarded when the spline engine works - it removes correspondent rows in the matrix before solving for basis coefficients). Is that behavior you don't want?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <go02o9$10$1@fred.mathworks.com>...
> "Travis" <sinusoid2@hotmail.com> wrote in message <gnvn35$fd1$1@fred.mathworks.com>...
> > "Travis" <sinusoid2@hotmail.com> wrote in message <gnpcmt$kad$1@fred.mathworks.com>...
>
> >
> > Anyone? I have tried putting 0's in there, but that of course doesn't work. It is more representative than the NaN's, but it still shows data where none was collected. Is there a way to get MATLab to ignore those areas?
>
> I have no idea what you meant by "ignore NaN". Please elaborate.
>
> When interpolate using spline method, the destination array get fully filled, even if there are few NaNs in the source array (Those NaN points are discarded when the spline engine works - it removes correspondent rows in the matrix before solving for basis coefficients). Is that behavior you don't want?
>
> Bruno

The NaN's (or at least the corners where they are), I would like to be more or less ignored. There is no data there, and I would like it to stay that way. Do I have to just setup something that will pick them out after the spline?

"Travis" <sinusoid2@hotmail.com> wrote in message <go03q1$a57$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <go02o9$10$1@fred.mathworks.com>...
> > "Travis" <sinusoid2@hotmail.com> wrote in message <gnvn35$fd1$1@fred.mathworks.com>...
> > > "Travis" <sinusoid2@hotmail.com> wrote in message <gnpcmt$kad$1@fred.mathworks.com>...
> >
> > >
> > > Anyone? I have tried putting 0's in there, but that of course doesn't work. It is more representative than the NaN's, but it still shows data where none was collected. Is there a way to get MATLab to ignore those areas?
> >
> > I have no idea what you meant by "ignore NaN". Please elaborate.
> >
> > When interpolate using spline method, the destination array get fully filled, even if there are few NaNs in the source array (Those NaN points are discarded when the spline engine works - it removes correspondent rows in the matrix before solving for basis coefficients). Is that behavior you don't want?
> >
> > Bruno
>
> The NaN's (or at least the corners where they are), I would like to be more or less ignored. There is no data there, and I would like it to stay that way. Do I have to just setup something that will pick them out after the spline?

Yes, you have to post process the output : detect the fine grid points in the patch affected by NaN, then set them to NaN.

Bruno

OK, here is what my data looks like:

NaN 6 12 18 24 30
0 29.6 29.4 28.5 28.1 28.5
1 25.3 23.7 23.5 24.7 24.4
2 24.2 21.9 21.8 23.8 23.4
3 23.1 20.8 20.9 23 22.4
4 22.3 20.2 20.1 22.2 21.8
4.5 22 19.85 20.05 21.7 21.6
5 NaN 19.5 20 21.2 NaN
6 NaN 19.9 NaN NaN NaN

With the X across the top, and the Y along the side. When I use interp2 with the spline it not only fills the NaN's, but the values at the points are actually changed. I just want to make the graph of this data a little finer amd smoother

"Travis" <sinusoid2@hotmail.com> wrote in message <go80u1$s46$1@fred.mathworks.com>...
> OK, here is what my data looks like:
>

I'm disappointed that you have showed little any work on your own.

img=peaks(15);

% Select randomly 10 NaN points
idxnan=randperm(numel(img));
img(idxnan(1:10))=NaN;

% Interpolation
m=8;
step=1/m;
xi=1:step:size(img,2);
yi=1:step:size(img,1);
[XI YI]=meshgrid(xi,yi);
imgs=interp2(img,XI,YI,'spline');

% Post processing, put back NaN value in
% fine surface
temp=nan(size(imgs,1)+2*(m-1),...
         size(imgs,2)+2*(m-1));
temp(m:end-m+1,m:end-m+1)=imgs;

[i j]=find(isnan(img));
I = sub2ind(size(temp), i(:)*m, j(:)*m);
[di dj]=ndgrid(-m+1:m-1,-m+1:m-1);
DI = sub2ind(size(temp), m+di, m+dj)-sub2ind(size(temp), m, m);
temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
imgs=temp(m:end-m+1,m:end-m+1);
clear temp

% Contour plot
clf
levels=-7:7;
subplot(1,2,1);
contour(img,levels);
subplot(1,2,2);
contour(imgs,levels);

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <go87jr$jtf$1@fred.mathworks.com>...
> "Travis" <sinusoid2@hotmail.com> wrote in message <go80u1$s46$1@fred.mathworks.com>...
> > OK, here is what my data looks like:
> >
>
> I'm disappointed that you have showed little any work on your own.
>
> img=peaks(15);
>
> % Select randomly 10 NaN points
> idxnan=randperm(numel(img));
> img(idxnan(1:10))=NaN;
>
> % Interpolation
> m=8;
> step=1/m;
> xi=1:step:size(img,2);
> yi=1:step:size(img,1);
> [XI YI]=meshgrid(xi,yi);
> imgs=interp2(img,XI,YI,'spline');
>
> % Post processing, put back NaN value in
> % fine surface
> temp=nan(size(imgs,1)+2*(m-1),...
> size(imgs,2)+2*(m-1));
> temp(m:end-m+1,m:end-m+1)=imgs;
>
> [i j]=find(isnan(img));
> I = sub2ind(size(temp), i(:)*m, j(:)*m);
> [di dj]=ndgrid(-m+1:m-1,-m+1:m-1);
> DI = sub2ind(size(temp), m+di, m+dj)-sub2ind(size(temp), m, m);
> temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
> imgs=temp(m:end-m+1,m:end-m+1);
> clear temp
>
> % Contour plot
> clf
> levels=-7:7;
> subplot(1,2,1);
> contour(img,levels);
> subplot(1,2,2);
> contour(imgs,levels);
>
> % Bruno

Thank you for the code bruno a few questions. in this line;
temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
What does the bsxfun command do, and the @plus? I am not on a computer right now that has a compatible version of MATLab that understands that.

Secondly, I have been putting a lot of time into this, and getting very frustrated by messing with it. I do not have an extensive knowledge of MATLab and all its associated code, so I mainly work with the simple stuff, and obviously that does not work in this case.

"Travis" <sinusoid2@hotmail.com> wrote in message <go9p60$dp7$1@fred.mathworks.com>...

>
> Thank you for the code bruno a few questions. in this line;
> temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
> What does the bsxfun command do, and the @plus? I am not on a computer right now that has a compatible version of MATLab that understands that.
>

To learn about bsxfun, here is a great page:
http://blogs.mathworks.com/loren/2008/08/04/comparing-repmat-and-bsxfun-performance/
 
% Code with comments

img=peaks(20);

% Select randomly 10 NaN points in peaks image
idxnan=randperm(numel(img));
img(idxnan(1:10))=NaN;

% 2D Interpolation by spline
m=8; % how much finer
step=1/m;
xi=1:step:size(img,2);
yi=1:step:size(img,1);
[XI YI]=meshgrid(xi,yi);
imgs=interp2(img,XI,YI,'spline');
clear XI YI

% Post processing, put back NaN value in
% fine surface
% temp is imsg, enlarged by (m-1) pixels on each sides
temp=nan(size(imgs,1)+2*(m-1),...
         size(imgs,2)+2*(m-1));
temp(m:end-m+1,m:end-m+1)=imgs;

% look for where NaN are located in the coarse image
[i j]=find(isnan(img));

% (i*m,j*m) is corresponding indices in temp
I = sub2ind(size(temp), i(:)*m, j(:)*m);

% We will enlarge NaN-locations by m-1 in each sides
% left/right/bottom/top so as to find the pixels
% falling in coarse patch around a given NaN corner
[di dj]=deal(-m+1:m-1,-m+1:m-1);

% Corresponding linear index-shift (in temp)
% same as: DI(k,l) = di(k)+dj(l)*size(temp,1)
% "@plus" because we *add* the expanded column
% vector di(:) to expanded row vector dj(:).'*...
DI = bsxfun(@plus,di(:),dj(:).'*size(temp,1));

% Lock the array temp with NaN
% Corresponding linear index shift (in temp)
% K(p,q) = I(p) + DI(q); p=1,...,; q=1,...
% temp(K)=NaN;
temp(bsxfun(@plus,I(:),DI(:).'))=NaN;

% strink back temp to imgs
imgs=temp(m:end-m+1,m:end-m+1);
clear temp % go to waste basket

% Contour plot
figure(1); clf
levels=-7:7;
subplot(2,2,1);
contour(img,levels);
subplot(2,2,2);
contour(xi,yi,imgs,levels);
subplot(2,2,3);
surf(img,'edgecolor','none');
subplot(2,2,4);
surf(xi,yi,imgs,'edgecolor','none');

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <go9rkq$pcb$1@fred.mathworks.com>...
> "Travis" <sinusoid2@hotmail.com> wrote in message <go9p60$dp7$1@fred.mathworks.com>...
>
> >
> > Thank you for the code bruno a few questions. in this line;
> > temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
> > What does the bsxfun command do, and the @plus? I am not on a computer right now that has a compatible version of MATLab that understands that.
> >
>
> To learn about bsxfun, here is a great page:
> http://blogs.mathworks.com/loren/2008/08/04/comparing-repmat-and-bsxfun-performance/
>
> % Code with comments
>
> img=peaks(20);
>
> % Select randomly 10 NaN points in peaks image
> idxnan=randperm(numel(img));
> img(idxnan(1:10))=NaN;
>
> % 2D Interpolation by spline
> m=8; % how much finer
> step=1/m;
> xi=1:step:size(img,2);
> yi=1:step:size(img,1);
> [XI YI]=meshgrid(xi,yi);
> imgs=interp2(img,XI,YI,'spline');
> clear XI YI
>
> % Post processing, put back NaN value in
> % fine surface
> % temp is imsg, enlarged by (m-1) pixels on each sides
> temp=nan(size(imgs,1)+2*(m-1),...
> size(imgs,2)+2*(m-1));
> temp(m:end-m+1,m:end-m+1)=imgs;
>
> % look for where NaN are located in the coarse image
> [i j]=find(isnan(img));
>
> % (i*m,j*m) is corresponding indices in temp
> I = sub2ind(size(temp), i(:)*m, j(:)*m);
>
> % We will enlarge NaN-locations by m-1 in each sides
> % left/right/bottom/top so as to find the pixels
> % falling in coarse patch around a given NaN corner
> [di dj]=deal(-m+1:m-1,-m+1:m-1);
>
> % Corresponding linear index-shift (in temp)
> % same as: DI(k,l) = di(k)+dj(l)*size(temp,1)
> % "@plus" because we *add* the expanded column
> % vector di(:) to expanded row vector dj(:).'*...
> DI = bsxfun(@plus,di(:),dj(:).'*size(temp,1));
>
> % Lock the array temp with NaN
> % Corresponding linear index shift (in temp)
> % K(p,q) = I(p) + DI(q); p=1,...,; q=1,...
> % temp(K)=NaN;
> temp(bsxfun(@plus,I(:),DI(:).'))=NaN;
>
> % strink back temp to imgs
> imgs=temp(m:end-m+1,m:end-m+1);
> clear temp % go to waste basket
>
> % Contour plot
> figure(1); clf
> levels=-7:7;
> subplot(2,2,1);
> contour(img,levels);
> subplot(2,2,2);
> contour(xi,yi,imgs,levels);
> subplot(2,2,3);
> surf(img,'edgecolor','none');
> subplot(2,2,4);
> surf(xi,yi,imgs,'edgecolor','none');
>
> % Bruno

Again, thank you Bruno. It works great until the:
I = sub2ind(size(temp), i(:)*m, j(:)*m);
it is trying to access something outside the matrix dimensions. I have messed around with it a fair amount, but to no real avail. Any suggestions?

>
> Again, thank you Bruno. It works great until the:
> I = sub2ind(size(temp), i(:)*m, j(:)*m);
> it is trying to access something outside the matrix dimensions. I have messed around with it a fair amount, but to no real avail. Any suggestions?

I guess because the m value is not 2,4, 8, 16... (2 power something); thus the finer grid does not have the expected length.

Try to replace the variable xi and yi with these:

xi=linspace(1,size(img,2),(size(img,2)-1)*m+1);
yi=linspace(1,size(img,1),(size(img,1)-1)*m+1);

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <goaqbp$cq2$1@fred.mathworks.com>...
>
> >
> > Again, thank you Bruno. It works great until the:
> > I = sub2ind(size(temp), i(:)*m, j(:)*m);
> > it is trying to access something outside the matrix dimensions. I have messed around with it a fair amount, but to no real avail. Any suggestions?
>
> I guess because the m value is not 2,4, 8, 16... (2 power something); thus the finer grid does not have the expected length.
>
> Try to replace the variable xi and yi with these:
>
> xi=linspace(1,size(img,2),(size(img,2)-1)*m+1);
> yi=linspace(1,size(img,1),(size(img,1)-1)*m+1);
>
> Bruno

Thank you very much again Bruno, works perfectly, except for one thing. When I change the xi and yi to represent the distance and depth, 6 to 30, and 0 to -6 respectively, i get values from 50+ to over 1000. Should I simply leave the xi and yi as indicated above, and change the values on the graph, or is there a better way?

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <goaqbp$cq2$1@fred.mathworks.com>...
>
> >
> > Again, thank you Bruno. It works great until the:
> > I = sub2ind(size(temp), i(:)*m, j(:)*m);
> > it is trying to access something outside the matrix dimensions. I have messed around with it a fair amount, but to no real avail. Any suggestions?
>
> I guess because the m value is not 2,4, 8, 16... (2 power something); thus the finer grid does not have the expected length.
>
> Try to replace the variable xi and yi with these:
>
> xi=linspace(1,size(img,2),(size(img,2)-1)*m+1);
> yi=linspace(1,size(img,1),(size(img,1)-1)*m+1);
>
> Bruno

OK, nevermind, I found a way around it. Again, THANK YOU!!!

The code has been working great on most of the data thus far, but some of the times it takes areas that should be 28 and after the code is run they are 36+. This doesn't happen every time, so I don't know what to do about it. The max on the dataset I mentioned is only 28.8. Why would this happen?

"Travis" <sinusoid2@hotmail.com> wrote in message <gp4l99$jb0$1@fred.mathworks.com>...
> The code has been working great on most of the data thus far, but some of the times it takes areas that should be 28 and after the code is run they are 36+. This doesn't happen every time, so I don't know what to do about it. The max on the dataset I mentioned is only 28.8. Why would this happen?

I guess is because Matlab interp2 uses a spline with the so-called NOT-A-KNOT condition, which is IMHO not a great choice, since it is not very stable with NaN data. A better choice is using spline with NATURAL condition (second derivative on the boundary). In my understanding you need spline toolbox for this.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gp5140$t4g$1@fred.mathworks.com>...
> "Travis" <sinusoid2@hotmail.com> wrote in message <gp4l99$jb0$1@fred.mathworks.com>...
> > The code has been working great on most of the data thus far, but some of the times it takes areas that should be 28 and after the code is run they are 36+. This doesn't happen every time, so I don't know what to do about it. The max on the dataset I mentioned is only 28.8. Why would this happen?
>
> I guess is because Matlab interp2 uses a spline with the so-called NOT-A-KNOT condition, which is IMHO not a great choice, since it is not very stable with NaN data. A better choice is using spline with NATURAL condition (second derivative on the boundary). In my understanding you need spline toolbox for this.
>
> Bruno

Thank you for the information. I will look into the spline toolbox.

"Travis" <sinusoid2@hotmail.com> wrote in message <gp5v5q$okc$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gp5140$t4g$1@fred.mathworks.com>...
> > "Travis" <sinusoid2@hotmail.com> wrote in message <gp4l99$jb0$1@fred.mathworks.com>...
> > > The code has been working great on most of the data thus far, but some of the times it takes areas that should be 28 and after the code is run they are 36+. This doesn't happen every time, so I don't know what to do about it. The max on the dataset I mentioned is only 28.8. Why would this happen?
> >
> > I guess is because Matlab interp2 uses a spline with the so-called NOT-A-KNOT condition, which is IMHO not a great choice, since it is not very stable with NaN data. A better choice is using spline with NATURAL condition (second derivative on the boundary). In my understanding you need spline toolbox for this.
> >
> > Bruno

OK, so I now have the spline toolbox, but after about 30 minytes of looking through the new functions and examples, I am still unsure of how change the 'not-a-knot' to a nautral interpolation.

Have been struggling with this for a while now. Does anyone have any suggestions?