Re: Re: bicubic spline interpolation in octave BIG OOPS!

[Hoxide Ma]

Thanks for your reply and your programs, that's great!

bicubic spline is the one of the three traditional
method sof the 2d interpolation problem, the other
methods is 'nearst' and 'linear.

you can get more information about the bicubic spline
interpolation from :

http://members.bellatlantic.net/%7Evze2vrva/design.html

.

I coding will based on this artile, too.

BTW: I have found that the code about bicubic in
Matlab    maybe is not nice.

--- "Robert A. Macy" <address@hidden>写道:

> Sent my reply before reading your email.  Standard
> protocol, right?
> 
> You're trying to create data where none existed
> before?  
> 
> For that, I use a program that linearizes and/or
> changes
> the function to a new ordinate base.  For example, I
> have a
> log function, I can then make it linear and add as
> many
> points as I'd like.  That was the only way I could
> use
> gsplot 3D and change linearize x ordinate, while
> maintaining quick plotting.  
> 
> But that program uses strict linear interpolation
> between
> *any* points which would cause "ramping" for you. 
> But then
> I run my smoothing programs and that removes all the
> ramping and I get very nice curves.  
> 
> However, it takes more than 2 passes *if* there are
> more
> than two "created" points between known data.  And
> so on.  
> 
> 
> 
> 
> On Sun, 17 Jul 2005 00:40:12 +0800 (CST)
>  Hoxide Ma <address@hidden> wrote:
> > When I try to use interp2.m in octave-forge, i
> found
> > that in the source code, there is no bicubic
> > interpolation method. 
> > "At the moment only 'linear' and 'nearest'
> methods" !
> > 
> > there is the code:
> > 
> >
>
%------------------------------------------------------
> > function F = cubic(arg1,arg2,arg3,arg4,arg5)
> > %CUBIC 2-D bicubic data interpolation.
> > %   CUBIC(...) is the same as LINEAR(....) except
> that
> > it uses
> > %   bicubic interpolation.
> > %   
> > %   This function needs about 7-8 times SIZE(XI)
> > memory to be available.
> > %
> > %   See also LINEAR.
> > 
> > %   Clay M. Thompson 4-26-91, revised 7-3-91,
> 3-22-93
> > by CMT.
> > 
> > %   Based on "Cubic Convolution Interpolation for
> > Digital Image
> > %   Processing", Robert G. Keys, IEEE Trans. on
> > Acoustics, Speech, and
> > %   Signal Processing, Vol. 29, No. 6, Dec. 1981,
> pp.
> > 1153-1160.
> > do_fortran_indexing  =1
> > 
> > if nargin==1, % cubic(z), Expand Z
> >   [nrows,ncols] = size(arg1);
> >   s = 1:.5:ncols; sizs = size(s);
> >   t = (1:.5:nrows)'; sizt = size(t);
> >   s = s(ones(sizt),:);
> >   t = t(:,ones(sizs));
> > 
> > elseif nargin==2, % cubic(z,n), Expand Z n times
> >   [nrows,ncols] = size(arg1);
> >   ntimes = floor(arg2);
> >   s = 1:1/(2^ntimes):ncols; sizs = size(s);
> >   t = (1:1/(2^ntimes):nrows)'; sizt = size(t);
> >   s = s(ones(sizt),:);
> >   t = t(:,ones(sizs));
> > 
> > elseif nargin==3, % cubic(z,s,t), No X or Y
> specified.
> >   [nrows,ncols] = size(arg1);
> >   s = arg2; t = arg3;
> > 
> > elseif nargin==4,
> >   error('Wrong number of input arguments.');
> > 
> > elseif nargin==5, % cubic(x,y,z,s,t), X and Y
> > specified.
> >   [nrows,ncols] = size(arg3);
> >   mx = prod(size(arg1)); my = prod(size(arg2));
> >   if any([mx my] ~= [ncols nrows]) & ...
> >      ~isequal(size(arg1),size(arg2),size(arg3))
> >     error('The lengths of the X and Y vectors must
> > match Z.');
> >   end
> >   if any([nrows ncols]<[3 3]), error('Z must be at
> > least 3-by-3.'); end
> >   s = 1 +
> (arg4-arg1(1))/(arg1(mx)-arg1(1))*(ncols-1);
> >   t = 1 +
> (arg5-arg2(1))/(arg2(my)-arg2(1))*(nrows-1);
> >   
> > end
> > 
> > if any([nrows ncols]<[3 3]), error('Z must be at
> least
> > 3-by-3.'); end
> > if ~isequal(size(s),size(t)), 
> >   error('XI and YI must be the same size.');
> > end
> > 
> > % Check for out of range values of s and set to 1
> > sout = find((s<1)|(s>ncols));
> > if length(sout)>0, s(sout) = ones(size(sout)); end
> > 
> > % Check for out of range values of t and set to 1
> > tout = find((t<1)|(t>nrows));
> > if length(tout)>0, t(tout) = ones(size(tout)); end
> > 
> > % Matrix element indexing
> > ndx = floor(t)+floor(s-1)*(nrows+2);
> > 
> > % Compute intepolation parameters, check for
> boundary
> > value.
> > if isempty(s), d = s; else d = find(s==ncols); end
> > s(:) = (s - floor(s));
> > if length(d)>0, s(d) = s(d)+1; ndx(d) =
> > ndx(d)-nrows-2; end
> > 
> > % Compute intepolation parameters, check for
> boundary
> > value.
> > if isempty(t), d = t; else d = find(t==nrows); end
> > t(:) = (t - floor(t));
> > if length(d)>0, t(d) = t(d)+1; ndx(d) = ndx(d)-1;
> end
> > d = [];
> > 
> > if nargin==5,
> >   % Expand z so interpolation is valid at the
> > boundaries.
> >   zz = zeros(size(arg3)+2);
> >   zz(1,2:ncols+1) =
> 3*arg3(1,:)-3*arg3(2,:)+arg3(3,:);
> >   zz(2:nrows+1,2:ncols+1) = arg3;
> >   zz(nrows+2,2:ncols+1) =
> > 3*arg3(nrows,:)-3*arg3(nrows-1,:)+arg3(nrows-2,:);
> >   zz(:,1) = 3*zz(:,2)-3*zz(:,3)+zz(:,4);
> >   zz(:,ncols+2) =
> > 3*zz(:,ncols+1)-3*zz(:,ncols)+zz(:,ncols-1);
> >   nrows = nrows+2; ncols = ncols+2;
> > else
> >   % Expand z so interpolation is valid at the
> > boundaries.
> >   zz = zeros(size(arg1)+2);
> >   zz(1,2:ncols+1) =
> 3*arg1(1,:)-3*arg1(2,:)+arg1(3,:);
> >   zz(2:nrows+1,2:ncols+1) = arg1;
> >   zz(nrows+2,2:ncols+1) =
> > 3*arg1(nrows,:)-3*arg1(nrows-1,:)+arg1(nrows-2,:);
> >   zz(:,1) = 3*zz(:,2)-3*zz(:,3)+zz(:,4);
> >   zz(:,ncols+2) =
> > 3*zz(:,ncols+1)-3*zz(:,ncols)+zz(:,ncols-1);
> >   nrows = nrows+2; ncols = ncols+2;
> > end
> > 
> > ndx
> > % Now interpolate using computationally efficient
> > algorithm.
> > t0 = ((2-t).*t-1).*t;
> > t1 = (3*t-5).*t.*t+2;
> > t2 = ((4-3*t).*t+1).*t;
> > t(:) = (t-1).*t.*t;
> > F     = ( zz(ndx).*t0 + zz(ndx+1).*t1 +
> zz(ndx+2).*t2
> > + zz(ndx+3).*t ) ...
> >         .* (((2-s).*s-1).*s);
> > ndx(:) = ndx + nrows;
> > F(:)  = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 +
> > zz(ndx+2).*t2 + zz(ndx+3).*t ) ...
> >         .* ((3*s-5).*s.*s+2);
> > ndx(:) = ndx + nrows;
> > F(:)  = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 +
> > zz(ndx+2).*t2 + zz(ndx+3).*t ) ...
> >         .* (((4-3*s).*s+1).*s);
> > ndx(:) = ndx + nrows;
> > F(:)  = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 +
> > zz(ndx+2).*t2 + zz(ndx+3).*t ) ...
> >        .* ((s-1).*s.*s);
> > F(:) = F/4;
> > 
> > % Now set out of range values to NaN.
> > if length(sout)>0, F(sout) = NaN; end
> 
=== message truncated ===


A student of Math in Soochow University in China.
Intrested in Math, Python, Twisted, Prolog and .NET/MONO.

__________________________________________________
赶快注册雅虎超大容量免费邮箱?
http://cn.mail.yahoo.com



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------