在 Matlab 中计算 3d 结构的梯度

Question

我尝试使用 Matlab 函数gradient 来计算体积的梯度。我使用quiver 来显示切片的渐变。

我使用一个立方体状的体积，它关于 x、y 和 z 轴对称。令我惊讶的是，所有切片的结果都不相同。实际上只有 xy 平面（Z-slice，最后一个图像）的结果是预期的结果。

我知道在计算图像边界处的渐变时会出现问题。但对我来说，边界的结果并不重要，所以我不在乎边界旁边的结果是否正确。对我来说，重要的是所有三张图片看起来都像最后一张。

谁能告诉我我的代码有什么问题？谢谢！

f=zeros(20,20,20);
space = 5;
f(:,:,space) = 1; f(:,:,end-space) = 1;
f(:,space,:) = 1; f(:,end-space,:) = 1;
f(space,:,:) = 1; f(end-space,:,:) = 1;
space = 4;
f(:,:,space) = 1; f(:,:,end-space) = 1;
f(:,space,:) = 1; f(:,end-space,:) = 1;
f(space,:,:) = 1; f(end-space,:,:) = 1;

size_iso = size(f);
x_slice = round(size_iso(1)/2);
y_slice = round(size_iso(2)/2);
z_slice = round(size_iso(3)/2);

% display the gradient of the edge map
[fx,fy,fz] = gradient(f,0.1); 
figure;
image(squeeze(f(x_slice,:,:))*50); colormap(gray(64)); hold on;
quiver(squeeze(fy(x_slice,:,:)),squeeze(fz(x_slice,:,:))); 
axis equal;
title(['edge map gradient of X-slice ', num2str(x_slice)]);

figure;
image(squeeze(f(:,y_slice,:))*50); colormap(gray(64)); hold on;
quiver(squeeze(fx(:,y_slice,:)),squeeze(fz(:,y_slice,:))); 
axis equal;
title(['edge map gradient of Y-slice ', num2str(y_slice)]);

figure;
image(squeeze(f(:,:,z_slice))*50); colormap(gray(64)); hold on;
quiver(squeeze(fx(:,:,z_slice)),squeeze(fy(:,:,z_slice))); 
axis equal;
title(['edge map gradient of Z-slice ', num2str(z_slice)]);

Answer 1

使用 3D 矩阵和坐标会稍微复杂一些。

例如

img = rand(10,30);
imagesc(img);
axis equal;

将显示 30 像素宽和 10 像素高的图像。 在 MatLab 中，当您显示图像时，它的第一个维度（行）实际上是绘图上的 Y 轴。第二维（列）是图上的 X 轴。 例如，请参阅http://www.mathworks.com/help/matlab/math/multidimensional-arrays.html

为了说明代码中的错误，请考虑简化示例：

% we need a 3D matrix with
% 10 points along the X-axis
% 20 points along the Y-axis
% 30 points along the Z-axis
f = rand(20,10,30); % note the order of numbers

size_iso = size(f),               % gives [20 10 30]
x_slice = round(size_iso(2)/2)    % gives 5
y_slice = round(size_iso(1)/2)    % gives 10
z_slice = round(size_iso(3)/2)    % gives 15

figure;
image(squeeze(f(:,x_slice,:))*50); colormap(gray(64)); hold on;
axis equal;
title(['X-slice ', num2str(x_slice)]);
% this code produces image 30 pixels wide and 20 pixels high
% Thus 1st dimension (vertical   axis) is actually the Y-axis
% Thus 2nd dimension (horizontal axis) is actually the Z-axis

figure;
image(squeeze(f(y_slice,:,:))*50); colormap(gray(64)); hold on;
axis equal;
title(['Y-slice ', num2str(y_slice)]);
% this code produces image 30 pixels wide and 10 pixels high
% Thus 1st dimension (vertical   axis) is actually the X-axis
% Thus 2nd dimension (horizontal axis) is actually the Z-axis

figure;
image(squeeze(f(:,:,z_slice))*50); colormap(gray(64)); hold on;
axis equal;
title(['Z-slice ', num2str(z_slice)]);
% this code produces 10 pixels wide and 20 pixels high
% Thus 1st dimension (vertical   axis) is actually the Y-axis
% Thus 2nd dimension (horizontal axis) is actually the X-axis

为了使您的代码正常工作，您不仅应注意切片图像中的尺寸顺序，还应注意它们被squeeze 函数移动的方式。 因此，您应该为后续的quiver 函数调用提供适当的坐标组合。

我修改了您的代码，以使用唯一值填充垂直于给定轴的平板，因此您应该能够更容易地区分它们。此外，为了相同的目的，我沿每个轴使用不同的尺寸。

xvalue=0.33;
yvalue=0.66;
zvalue=1.00;

% we need a 3D matrix with
% 10 points along the X-axis
% 20 points along the Y-axis
% 30 points along the Z-axis
f = zeros(20,10,30); % note the order of numbers

space = 3;
f(:,space,:) = xvalue; f(:,end-space,:) = xvalue;
f(space,:,:) = yvalue; f(end-space,:,:) = yvalue;
f(:,:,space) = zvalue; f(:,:,end-space) = zvalue;

size_iso = size(f);
x_slice = round(size_iso(2)/2); % note dimension number here for x_slice
y_slice = round(size_iso(1)/2); % note dimension number here for y_slice
z_slice = round(size_iso(3)/2);

% display the gradient of the edge map
[fx,fy,fz] = gradient(f,0.1);

figure;
image(squeeze(f(:,x_slice,:))*50); colormap(gray(64)); hold on;
quiver(squeeze(fz(:,x_slice,:)),squeeze(fy(:,x_slice,:)));
axis equal;
title(['edge map gradient of X-slice ', num2str(x_slice)]);
xlabel('Z')
ylabel('Y')

figure;
image(squeeze(f(y_slice,:,:))*50); colormap(gray(64)); hold on;
quiver(squeeze(fz(y_slice,:,:)),squeeze(fx(y_slice,:,:)));
axis equal;
title(['edge map gradient of Y-slice ', num2str(y_slice)]);
xlabel('Z')
ylabel('X')

figure;
image(squeeze(f(:,:,z_slice))*50); colormap(gray(64)); hold on;
quiver(squeeze(fx(:,:,z_slice)),squeeze(fy(:,:,z_slice)));
axis equal;
title(['edge map gradient of Z-slice ', num2str(z_slice)]);
xlabel('X')
ylabel('Y')

是的，这很棘手，一开始很难理解，但通过练习你会习惯的。