Re: bicubic spline interpolation in octave BIG OOPS!

[Robert A. Macy]

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
> if length(tout)>0, F(tout) = NaN; end
> 
> endfunction
> 
> 
> 
> by changeing the interp2.m code in matlab, I can do
> bicubic spline interpolation in some case, you can see
> a demo ploted by gnuplot.
> 
> But, you know, the work haven't finished, it not work
> with interp2.m, and it use the code of matlab, there
> may be some problem in copyright. I plan to rewrite
> this codes.   
> 
> I write this mail to ask for some suggestions.
> 
> thx.
> 
> 
>       
> 
>       
>               
>
___________________________________________________________
> 
> 雅虎免费G邮箱－中国第一绝无垃圾邮件骚扰超大邮箱
> http://cn.mail.yahoo.com/?id=77071



-------------------------------------------------------------
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
-------------------------------------------------------------