逆離散cosine変換

離散cosine変換(discrete cosine transformation :DCT)はspm_normalise.m内のサブ関数snbasisで行います。 snbasis内で予め三次元離散cosine変換のx,y,z方向の基底ベクトルとその微分を計算した後、 本体であるspm_brainwarp.dllに上記演算結果と、標準脳、被験者脳の三次元配列、Affine変換行列、半値幅などを渡し、 最終目的であるTrという4次元配列を得ます。
Matlabのsignal processing toolboxの関数idctを使い、Trの逆離散cosine変換を行ってみます。
以下のような 関数Mファイル を作成しました。
function R=DrawTr(Tr,cutoff);
scale=1;
nx=1:size(Tr,1); % x方向のvoxel数
ny=1:size(Tr,2); % y方向のvoxel数
nz=1:size(Tr,3); % z方向のvoxel数
[nx,ny,nz]=ndgrid(nx,ny,nz);
nx=nx*cutoff;
ny=ny*cutoff;
nz=nz*cutoff;
vx=squeeze(Tr(:,:,:,1)); % 各voxelのx方向の移動量
vy=squeeze(Tr(:,:,:,2)); % 各voxelのy方向の移動量
vz=squeeze(Tr(:,:,:,3)); % 各voxelのz方向の移動量
for x=1:size(Tr,1);for y=1:size(Tr,2);...
vx(x,y,:)=idct(squeeze(vx(x,y,:)));...
vy(x,y,:)=idct(squeeze(vy(x,y,:)));...
vz(x,y,:)=idct(squeeze(vz(x,y,:)));...
end;end;
for y=1:size(Tr,2);for z=1:size(Tr,3);...
vx(:,y,z)=idct(squeeze(vx(:,y,z)));...
vy(:,y,z)=idct(squeeze(vy(:,y,z)));...
vz(:,y,z)=idct(squeeze(vz(:,y,z)));...
end;end;
for x=1:size(Tr,1);for z=1:size(Tr,3);...
vx(x,:,z)=idct(squeeze(vx(x,:,z)));...
vy(x,:,z)=idct(squeeze(vy(x,:,z)));...
vz(x,:,z)=idct(squeeze(vz(x,:,z)));...
end;end;
R.vx=vx;
R.vy=vy;
R.vz=vz;
figure;set(gcf,'color',[1,1,1],'position',[100,100,600,250]);
width=0.28;width2=[0.05,width,0.9];
subplot('position',[0.05,width2]);
DrawQuiver(nx,ny,nz,vx,vy,vz,scale);view([0,0]);
subplot('position',[0.10+width,width2]);
DrawQuiver(nx,ny,nz,vx,vy,vz,scale);view([0,90]);
subplot('position',[0.15+width*2,width2]);
DrawQuiver(nx,ny,nz,vx,vy,vz,scale);view([90,0]);

function DrawQuiver(nx,ny,nz,vx,vy,vz,scale);
h=quiver3(nx,ny,nz,vx,vy,vz,scale);
set(findobj(h,'type','line'),'LineWidth',2);
hold on;
plot3(nx(:),ny(:),nz(:),'b.');
daspect([1,1,1]);
axis tight;
voxelの大きさと思われるcutoff値を20,25,30,35,40と変えてみます。 因みにcutoff値15ではメモリが足りませんと表示されました。
V{1}=DrawTr(Tr,20); % cutoff=20mm
V{2}=DrawTr(Tr,25); % cutoff=25mm
V{3}=DrawTr(Tr,30); % cutoff=30mm
V{4}=DrawTr(Tr,35); % cutoff=35mm
V{5}=DrawTr(Tr,40); % cutoff=40mm
20mm
25mm
30mm
35mm
40mm
図は正面、側面、上から見たところの順になっています。
同じような傾向ですが、cutoff値が大きいと、矢印が大きくなっているようです。
x,y,z方向の絶対値の平均をとりました。
X=[mean(abs(V{1}.vx(:))),mean(abs(V{1}.vy(:))),mean(abs(V{1}.vz(:)));...% 20mm
mean(abs(V{2}.vx(:))),mean(abs(V{2}.vy(:))),mean(abs(V{2}.vz(:)));...% 25mm
mean(abs(V{3}.vx(:))),mean(abs(V{3}.vy(:))),mean(abs(V{3}.vz(:)));...% 30mm
mean(abs(V{4}.vx(:))),mean(abs(V{4}.vy(:))),mean(abs(V{4}.vz(:)));...% 35mm
mean(abs(V{5}.vx(:))),mean(abs(V{5}.vy(:))),mean(abs(V{5}.vz(:)))] % 40mm
cutoffvxvyvz
2021.021231.967429.6468
2529.923045.425442.2949
3040.727360.534355.6385
3552.129278.410972.9350
4057.402184.963379.5914

cutoff=25の値で補正したものを図示します。
x=20:5:40;
figure;set(gcf,'color',[1,1,1]);
plot(x,X./(ones(5,1)*X(2,:)),'Marker','o');grid on;
b=[ones(4,1),x(1:4)']\mean(X(1:4,:)./(ones(4,1)*X(2,:)),2)
b =ans
-0.6896
0.0684
hold on;plot(x,x*b(2)+b(1),'k','LineWidth',2);