Matlab中用于转换3D图像的插值或重采样算法,最好是Sinc插值 [英] Interpolating or resampling algorithm in matlab for transformed 3D image, preferably sinc interpolation
问题描述
我有一个3D数据集和一个2D数据集,这是第一个体积的一部分.它们具有不同的比例,分辨率和坐标系,但是我都知道对世界坐标的仿射变换.然后,我想我知道如何应用这些图像,但是如何使用Sinc插值再次从这些转换后的坐标中获得图像呢?我想学习如何执行此操作/如何工作.下面的第一条评论已经向我介绍了matlab中执行线性插值的现有函数,但是我也想知道如何自己执行此操作,以便可以使用sinc插值(和其他函数).
我可以四舍五入并获得这些坐标的值,这将是最近邻插值.无论计算时间如何,我都希望丢失尽可能少的信息,然后我应该使用sinc插值.坐标变换后,如何制作(例如sinc)插值算法?
例如:
%%get data
A = rand(200,250,250); % I want to get the slice from this that corresponds to B
B = rand(1200,1200); % I want to get the data from A that corresponds to the same pixels
%%get coordinates for A
siza = size(A); clear xx;clear yy;clear zz;
[xx,yy,zz] = meshgrid(1:siza(1), 1:siza(2), 1:siza(3));
coora = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%%get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2)); zz = zeros(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%%define affine transformation matrices
T3d = [-0.02 0.02 1 -88 ;
-0.98 0 -0.02 130;
0 0.98 -0.02 -110;
0 0 0 1 ];
T2d = [-0.2 0 0 126;
0 0.2 -0.2 -131;
0 0 2 43 ;
0 0 0 1 ];
%%transform A coordinates to world coordinates and world coordinates to B coordinates
cooraInBref = T3d*inv(T2d)*coora;
aslice = zeros(size(B));
%% then nearest neighbor would go something like this (dont exactly know how to do this either):
cooraInBround = round(cooraInBref);
for idx = 1:length(coorb);
if cooraInBround(3,idx) == 0
aslice(cooraInBround(1,idx),cooraInBround(2,idx)) = ...;% dont know how to do this
end
end
%% how would I implement sinc interpolation instead of rounding the transformed coordinates
不能进一步帮助我的相关问题:
如何将仿射变换(4x4矩阵)应用于ndgrid/meshgrid结果?
如chappjc anonsubmitter85和cape代码的注释中所指出的那样,可以立即使用matlab中的函数:
http://mathworks.com/help/matlab/math/interpolating-scattered-data.html
http://mathworks.com/help/matlab/ref/griddata.html?refresh=true
使用注释中建议的散点插值法,所链接的问题花费了将近十分钟,并导致包含所有NaN的切片,现在正在尝试其他方法.
我现在以另一种方式开始:获取要采样的2d点,然后将其转换为3d体积,然后使用interp3计算值.它的速度很快,效果很好.该代码仅适用于获取单个切片,但是我认为您可以轻松地对其进行调整以获取整个转换后的体积.我仍然不知道如何进行正弦插值.
%% transformation from one image to the other
Affine = inv(T2d)*T3d
%% get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2));
zz = ones(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%% transformed coordinates
coorb_t = Affine*coorb;
idxX = reshape(coorb_t(:,1), sizb(1), sizb(2), 1);
idxY = reshape(coorb_t(:,2), sizb(1), sizb(2), 1);
idxZ = reshape(coorb_t(:,3), sizb(1), sizb(2), 1);
%% interpolate
Asliced = interp3(A, idxX, idxY, idxZ, 'cubic');
仍然不确定我是否应该对Z使用零或一.
I have one 3D dataset and one 2D dataset which is a slice through the first volume. They are at different scales, resolutions and in a different coordinate system, but of both I know the affine transformation to world coordinates. I then think I know how to apply those, but how do I get an image back from those transformed coordinates again using sinc interpolation? I would like to learn how to perform this/how this works. The first comments below already pointed me at existing functions inside matlab that perform linear interpolation, but I would also like to know how to do this myself so I can use sinc interpolation (and others).
I can round the coordinates and get values for those, which would be nearest neighbor interpolation. I would like to lose as little information as possible irrespective of computation time, I think I should use sinc interpolation then. When I have the transformed coordinates, how do I make (for example sinc) an interpolation algorithm?
for example:
%%get data
A = rand(200,250,250); % I want to get the slice from this that corresponds to B
B = rand(1200,1200); % I want to get the data from A that corresponds to the same pixels
%%get coordinates for A
siza = size(A); clear xx;clear yy;clear zz;
[xx,yy,zz] = meshgrid(1:siza(1), 1:siza(2), 1:siza(3));
coora = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%%get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2)); zz = zeros(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%%define affine transformation matrices
T3d = [-0.02 0.02 1 -88 ;
-0.98 0 -0.02 130;
0 0.98 -0.02 -110;
0 0 0 1 ];
T2d = [-0.2 0 0 126;
0 0.2 -0.2 -131;
0 0 2 43 ;
0 0 0 1 ];
%%transform A coordinates to world coordinates and world coordinates to B coordinates
cooraInBref = T3d*inv(T2d)*coora;
aslice = zeros(size(B));
%% then nearest neighbor would go something like this (dont exactly know how to do this either):
cooraInBround = round(cooraInBref);
for idx = 1:length(coorb);
if cooraInBround(3,idx) == 0
aslice(cooraInBround(1,idx),cooraInBround(2,idx)) = ...;% dont know how to do this
end
end
%% how would I implement sinc interpolation instead of rounding the transformed coordinates
Related questions that dont quite help me further:
Nearest-neighbor interpolation algorithm in MATLAB
"Nearest neighbor"-like interpolation in MATLAB
How to apply an affine transformation (4x4 matrix) to ndgrid/meshgrid results?
Python/PIL affine transformation
How do I rotate a 3D matrix by 90 degrees counterclockwise?
Resizing a 3D image (and resampling)
Ready to use functions from matlab as pointed out in comments by chappjc anonsubmitter85 and cape code:
http://mathworks.com/help/matlab/math/interpolating-scattered-data.html
http://mathworks.com/help/matlab/ref/griddata.html?refresh=true
My other SE question I use to get the affine matrix
Using scatteredInterpolant as suggested in the comments and the linked question took almost ten minutes and resulted in a slice with all NaN`s, am now trying other methods.
I now start the other way around: get the 2d points to sample and then transform those to the 3d volume and then use interp3 to calculate the values. Its pretty quick and works well. This code only works for getting a single slice, but i think you can easily adapt it to get a whole transformed volume. I still dont know how to do sinc interpolation though.
%% transformation from one image to the other
Affine = inv(T2d)*T3d
%% get coordinates for B
sizb = size(B); clear xx;clear yy;clear zz;
[xx,yy] = meshgrid(1:sizb(1),1:sizb(2));
zz = ones(size(xx));
coorb = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
%% transformed coordinates
coorb_t = Affine*coorb;
idxX = reshape(coorb_t(:,1), sizb(1), sizb(2), 1);
idxY = reshape(coorb_t(:,2), sizb(1), sizb(2), 1);
idxZ = reshape(coorb_t(:,3), sizb(1), sizb(2), 1);
%% interpolate
Asliced = interp3(A, idxX, idxY, idxZ, 'cubic');
Still not sure if I should have used zeros or ones for Z.
这篇关于Matlab中用于转换3D图像的插值或重采样算法,最好是Sinc插值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!