splin3d

// 例1
// =============================================================================
func =  "v=cos(2*%pi*x).*sin(2*%pi*y).*cos(2*%pi*z)";
deff("v=f(x,y,z)",func);
n = 10;  // n x n x n  補間点
x = linspace(0,1,n); y=x; z=x; // 補間グリッド
[X,Y,Z] = ndgrid(x,y,z);
V = f(X,Y,Z);
tl = splin3d(x,y,z,V,[5 5 5]);
m = 10000;
// 近似誤差を計算
xp = grand(m,1,"def"); yp = grand(m,1,"def"); zp = grand(m,1,"def");
vp_exact = f(xp,yp,zp);
vp_interp = interp3d(xp,yp,zp, tl);
er = max(abs(vp_exact - vp_interp))
// n=20で再試行し,誤差を見る

// 例 2 (linear_interpn のヘルプを参照ください, 
//       trilinear補間に関する同じ例があります)
// =============================================================================
exec("SCI/modules/interpolation/demos/interp_demo.sci")
func =  "v=(x-0.5).^2 + (y-0.5).^3 + (z-0.5).^2";
deff("v=f(x,y,z)",func);
n = 5; 
x = linspace(0,1,n); y=x; z=x;
[X,Y,Z] = ndgrid(x,y,z);
V = f(X,Y,Z);
tl = splin3d(x,y,z,V);
// いくつかの断面で3次元スプラインを計算(および表示)
m = 41;
dir = ["z="  "z="  "z="  "x="  "y="];
val = [ 0.1   0.5   0.9   0.5   0.5];
ebox = [0 1 0 1 0 1];
XF=[]; YF=[]; ZF=[]; VF=[];
for i = 1:length(val)
  [Xm,Xp,Ym,Yp,Zm,Zp] = slice_parallelepiped(dir(i), val(i), ebox, m, m, m);
  Vm = interp3d(Xm,Ym,Zm, tl);
  [xf,yf,zf,vf] = nf3dq(Xm,Ym,Zm,Vm,1);
  XF = [XF xf]; YF = [YF yf]; ZF = [ZF zf]; VF = [VF vf]; 
  Vp = interp3d(Xp,Yp,Zp, tl);
  [xf,yf,zf,vf] = nf3dq(Xp,Yp,Zp,Vp,1);
  XF = [XF xf]; YF = [YF yf]; ZF = [ZF zf]; VF = [VF vf]; 
end
nb_col = 128;
vmin = min(VF); vmax = max(VF);
color = dsearch(VF,linspace(vmin,vmax,nb_col+1));
xset("colormap",jetcolormap(nb_col));
clf(); xset("hidden3d",xget("background"));
colorbar(vmin,vmax)
plot3d(XF, YF, list(ZF,color), flag=[-1 6 4])
xtitle("3d spline interpolation of "+func)
show_window()