三次元線形補間(MEX)

spm_sample_vol.dllはなかなか秀逸なのですが、ndgridを使って得られたx,y,zの三次元配列も必要となるため、メモリを圧迫します。 そこでhns_trilinear.dllなるMEXファイルを作成しました。
使い方は
new_3D=hns_trilinaer(old_3D,x_vector,y_vector,z_vector);
old_3Dはtargetのdoubleの三次元配列、以下x,y,z軸上の位置を意味します。
spm_sample_vol.dllと異なり、reshapeで三次元配列にしなおす必要はありません。
old_3Dが、double以外のときはSegmentation violation・・・とエラーが出ます。
hns_trilinear.dll (binary)
hns_trilinear.cpp (source)
確かめてみます。
A=rand(100,100,100);
x=linspace(1,100,200);
tic;B=hns_trilinear(A,x,x,x);toc
elapsed_time=3.0780
[x,y,z]=ndgrid(x,x,x);
tic;C=spm_sample_vol(A,x,y,z,1);toc
elapsed_time=3.5310;
max(abs(B(:)-C(:)))
ans =4.4409e-016
実際のデータで試して見ます。
clear all
[filename,pathname]=uigetfile('*.img'); % *.img fileを選択
VF=spm_vol(strcat(pathname,filename)); % 被験者脳のheaderを読む
VF.dat=spm_read_vols(VF); % 被験者脳の三次元配列を読む
x=linspace(1,VF.dim(1),abs(VF.mat(1,1)*VF.dim(1))); % =linspace(1,256,240);
y=linspace(1,VF.dim(2),abs(VF.mat(2,2)*VF.dim(2))); % =linspace(1,256,240);
z=linspace(1,VF.dim(3),abs(VF.mat(3,3)*VF.dim(3))); % =linspace(1,172,86);
tic;W=hns_trilinear(VF.dat,x,y,z);toc
figure;imagesc(rot90(squeeze(W(120,:,:))),[0,384]);
axis equal;axis tight;colormap(gray);
elapsed_time=3.8120

メモリ消費量が1/4になり経済的になりました。 また、ちょっとだけ高速化できました。
どうせなら、浮動小数点8byteのdoubleでなく、整数2byteのint16で補間したいものです。
hns_trilinear16.dll (binary)
hns_trilinear16.cpp (source)
x=linspace(1,256,480);y=x;z=linspace(1,183,86*4);
W=hns_trilinear16(int16(VF.dat),x,y,z);
figure;imagesc(rot90(squeeze(W(240,:,:))),[0,384]);
axis equal;axis tight;colormap(gray);

メモリ不足でspm_sample_vol.dllではできなかったことができました。めでたしめでたし。