octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

More text and figures for interpolation chapter


From: David Bateman
Subject: More text and figures for interpolation chapter
Date: Thu, 14 Jun 2007 03:51:43 +0200
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

Here is a patch that adds a little more text to the interpolation
chapter. It also adds additional figures to this chapter. It uses a
script to create these in the same manner as the sparse figures.

D.
Index: doc/interpreter/Makefile.in
===================================================================
RCS file: /usr/local/cvsroot/octave/doc/interpreter/Makefile.in,v
retrieving revision 1.147
diff -u -r1.147 Makefile.in
--- doc/interpreter/Makefile.in 13 Jun 2007 05:57:02 -0000      1.147
+++ doc/interpreter/Makefile.in 14 Jun 2007 01:39:36 -0000
@@ -18,7 +18,7 @@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_DATA = @INSTALL_DATA@
 
-SCRIPT_SOURCES = sparseimages.m
+SCRIPT_SOURCES = sparseimages.m interpimages.m
 
 EXAMPLE_FILES_NODIR = \
   addtwomatrices.cc \
@@ -43,16 +43,20 @@
 
 EXAMPLE_FILES = $(addprefix $(top_srcdir)/examples/, $(EXAMPLE_FILES_NODIR))
 
-SPARSEIMAGES_1 = gplot grid spmatrix spchol spcholperm
+INTERPIMAGES = interpft interpn interpderiv
+INTERPIMAGES_EPS = $(addsuffix .eps, $(INTERPIMAGES))
+INTERPIMAGES_PDF = $(addsuffix .pdf, $(INTERPIMAGES))
+INTERPIMAGES_PNG = $(addsuffix .png, $(INTERPIMAGES))
 
+SPARSEIMAGES_1 = gplot grid spmatrix spchol spcholperm
 SPARSEIMAGES_EPS = $(addsuffix .eps, $(SPARSEIMAGES_1))
 SPARSEIMAGES_PDF = $(addsuffix .pdf, $(SPARSEIMAGES_1))
 SPARSEIMAGES_PNG = $(addsuffix .png, $(SPARSEIMAGES_1))
 SPARSEIMAGES_TXT = $(addsuffix .txt, $(SPARSEIMAGES_1))
 
-IMAGES_EPS = $(SPARSEIMAGES_EPS)
-IMAGES_PDF = $(SPARSEIMAGES_PDF)
-IMAGES_PNG = $(SPARSEIMAGES_PNG)
+IMAGES_EPS = $(SPARSEIMAGES_EPS) $(INTERPIMAGES_EPS)
+IMAGES_PDF = $(SPARSEIMAGES_PDF) $(INTERPIMAGES_PDF)
+IMAGES_PNG = $(SPARSEIMAGES_PNG) $(INTERPIMAGES_PNG)
 IMAGES_TXT = $(SPARSEIMAGES_TXT)
 
 HTML_IMAGES_PNG = $(addprefix HTML/, $(IMAGES_PNG))
@@ -207,6 +211,9 @@
     --eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst 
.%,%, $(suffix $@))'); sleep (1);"
 endef
 
+$(INTERPIMAGES_EPS) $(INTERPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
+       $(run-octave)
+
 $(SPARSEIMAGES_EPS) $(SPARSEIMAGES_PNG) $(SPARSEIMAGES_TXT): sparseimages.m
        $(run-octave)
 
Index: doc/interpreter/interp.txi
===================================================================
RCS file: /usr/local/cvsroot/octave/doc/interpreter/interp.txi,v
retrieving revision 1.3
diff -u -r1.3 interp.txi
--- doc/interpreter/interp.txi  12 Jun 2007 21:39:27 -0000      1.3
+++ doc/interpreter/interp.txi  14 Jun 2007 01:39:36 -0000
@@ -15,6 +15,43 @@
 
 @DOCSTRING(interp1)
 
+There are some important differences between the various interpolation
+methods. The 'spline' method enforces that both the first and second
+derivatives of the interpolated values have a continuous derivative,
+whereas the other methods do not. This can be demonstrated by the code
+
address@hidden
address@hidden
+t = 0 : 0.3 : pi; dt = t(2)-t(1);
+n = length (t); k = 100; dti = dt*n/k;
+ti = t(1) + [0 : k-1]*dti;
+y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+      ti(2:end-1),ddyp,'c^');
+legend('cubic','spline','pchip');
address@hidden group
address@hidden example
+
address@hidden
address@hidden
+The result of which can be seen in @ref{fig:interpderiv}.
+
address@hidden Figure,fig:interpderiv
address@hidden,8cm}
address@hidden of second derivative of interpolated values for
+various interpolation methods}
address@hidden float
address@hidden ifnotinfo
+
+This means that in general the 'spline' method results in smooth
+results. If the function to be interpolated is in fact smooth, then
+'spline' will give excellent results. However, if the function to be
+evaluated is in some manner discontinuous, then 'cubic' or 'pchip'
+interpolation might give better results.
+
 Fourier interpolation, is a resampling technique where a signal is
 converted to the frequency domain, padded with zeros and then
 reconverted to the time domain.
@@ -40,8 +77,20 @@
 @end group
 @end example
 
address@hidden
 which demonstrates the poor behavior of Fourier interpolation for non
 periodic functions.
address@hidden ifinfo
address@hidden
+which demonstrates the poor behavior of Fourier interpolation for non
+periodic functions, as can be seen in @ref{fig:interpft}.
+
address@hidden Figure,fig:interpft
address@hidden,8cm}
address@hidden of @code{interp1} and @code{interpft} for non
+periodic data}
address@hidden float
address@hidden ifnotinfo
 
 In additional the support function @code{spline} and @code{lookup} that
 underlie the @code{interp1} function can be called directly.
@@ -81,11 +130,12 @@
 f = @@(x,y,z) x.^2 - y - z.^2;
 [xx, yy, zz] = meshgrid (x, y, z);
 v = f (xx,yy,zz);
-xi = yi = zi = -1:0.5:1;
+xi = yi = zi = -1:0.1:1;
 [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
-vi = interp3(x, y, z, v, xxi, yyi, zzi);
+vi = interp3(x, y, z, v, xxi, yyi, zzi, 'spline');
 [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
-vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+vi2 = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+mesh (yi, zi, squeeze (vi2(1,:,:)));
 @end group
 @end example
 
@@ -93,6 +143,14 @@
 where @code{vi} and @code{vi2} are identical. The reversal of the
 dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions
 respectively.
address@hidden
+The result of this code can be seen in @ref{fig:interpn}.
+
address@hidden Figure,fig:interpn
address@hidden,8cm}
address@hidden of the use of @code{interpn}}
address@hidden float
address@hidden ifnotinfo
 
 In additional the support function @code{bicubic} that underlies the
 cubic interpolation of @code{interp2} function can be called directly.
Index: scripts/general/__splinen__.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/__splinen__.m,v
retrieving revision 1.1
diff -u -r1.1 __splinen__.m
--- scripts/general/__splinen__.m       12 Jun 2007 21:39:27 -0000      1.1
+++ scripts/general/__splinen__.m       14 Jun 2007 01:39:37 -0000
@@ -28,14 +28,14 @@
   if (nargin != 5)
     error ("Incorrect number of arguments");
   endif
-
-  if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (@isvector, x)) ||
-      !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (@isvector, xi)))
+  isvec = @(x) prod(size(x)) == length(x);   # ND isvector function
+  if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (isvec, x)) ||
+      !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (isvec, xi)))
     error ("%s: non gridded data or dimensions inconsistent", f);
   endif
   yi = y;
   for i = length(x):-1:1
-    yi = spline (x{i}, yi, xi{i}).';
+    yi = permute (spline (x{i}, yi, xi{i}), [length(x),1:length(x)-1]);
   endfor
 
   [xi{:}] = ndgrid (xi{:});
Index: scripts/general/interp1.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/interp1.m,v
retrieving revision 1.5
diff -u -r1.5 interp1.m
--- scripts/general/interp1.m   12 Jun 2007 21:39:27 -0000      1.5
+++ scripts/general/interp1.m   14 Jun 2007 01:39:37 -0000
@@ -345,6 +345,19 @@
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
+%!demo
+%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
+%! n = length (t); k = 100; dti = dt*n/k;
+%! ti = t(1) + [0 : k-1]*dti;
+%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+%! ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+%! ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+%! ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+%! plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+%!       ti(2:end-1),ddyp,'c^');
+%! legend('cubic','spline','pchip');
+%! title("Second derivative of interpolated 'sin (4*t + 0.3) .* cos (3*t - 
0.1)'");
+
 ## For each type of interpolated test, confirm that the interpolated
 ## value at the knots match the values at the knots.  Points away
 ## from the knots are requested, but only 'nearest' and 'linear'
Index: scripts/general/interpn.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/interpn.m,v
retrieving revision 1.1
diff -u -r1.1 interpn.m
--- scripts/general/interpn.m   12 Jun 2007 21:39:27 -0000      1.1
+++ scripts/general/interpn.m   14 Jun 2007 01:39:37 -0000
@@ -128,22 +128,22 @@
       endif
       idx (1 : nd) = {1};
       idx (i) = ":";
-      x{i} = x{i}(idx{:});
+      x{i} = x{i}(idx{:})(:);
     endfor
     idx (1 : nd) = {1};
     idx (1) = ":";
-    x{1} = x{1}(idx{:});
+    x{1} = x{1}(idx{:})(:);
   endif
 
   if (strcmp (method, "linear") || strcmp (method, "nearest"))
     if (all (cellfun (@isvector, y)))
       [y{:}] = ndgrid (y{:});
     endif
-  elseif (any (! cellfun (@isvector, x)))
+  elseif (any (! cellfun (@isvector, y)))
     for i = 1 : nd
       idx (1 : nd) = {1};
       idx (i) = ":";
-      y{i} = y{i}(idx{:});
+      y{i} = y{i}(idx{:})(:).';
     endfor
   endif
 
@@ -170,7 +170,7 @@
     endfor
     vi(idx) = extrapval;
     vi = reshape (vi, yshape); 
-  elseif (strcmp (method, "spline")) 
+  elseif (strcmp (method, "spline"))
     vi = __splinen__ (x, v, y, extrapval, "interpn");
   elseif (strcmp (method, "cubic")) 
     error ("cubic interpolation not yet implemented");
@@ -216,3 +216,14 @@
 %! [x,y] = meshgrid(x,y); 
 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 
+
+%!demo
+%! x = y = z = -1:1;
+%! f = @(x,y,z) x.^2 - y - z.^2;
+%! [xx, yy, zz] = meshgrid (x, y, z);
+%! v = f (xx,yy,zz);
+%! xi = yi = zi = -1:0.1:1;
+%! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+%! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+%! mesh (yi, zi, squeeze (vi(1,:,:)));
+
Index: scripts/polynomial/mkpp.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/polynomial/mkpp.m,v
retrieving revision 1.3
diff -u -r1.3 mkpp.m
--- scripts/polynomial/mkpp.m   22 Jan 2007 17:29:53 -0000      1.3
+++ scripts/polynomial/mkpp.m   14 Jun 2007 01:39:37 -0000
@@ -52,7 +52,7 @@
     d = round (rows (P) / pp.n); 
   endif
   pp.d = d;
-  if (pp.n*d != rows (P))
+  if (pp.n * prod (d) != rows (P))
     error ("mkpp: num intervals in x doesn't match num polynomials in P");
   endif
 endfunction
2007-06-14  David Bateman  <address@hidden>

        * interpreter/Makefile.in (SCRIPT_SORCES): add interimages.m
        (INTERPIMAGES*): New variables. Add targets for them
        (IMAGES_EPS, IMAGES_PDF, IMAGES_PNG): Add the INTERPIMAGES.
        * interpreter/interpimages.m: New function
        * interpreter/interp.txi: Add text about second derivation of
        splines and add figures.

2007-06-14  David Bateman  <address@hidden>

        * general/__splinen__.m: Check also for ND vectors. Fix for N > 2,
        as permutation of results was incorrect.
        * general/interp1.m: Add demo on second derivative
        * general/interpn.m: Convert "y" to vectors for __splinen__
        call. Add 3D demo.

        * polynomial/mkpp.m: Correction for matrices of 3 or more dimensions.

reply via email to

[Prev in Thread] Current Thread [Next in Thread]