Radon变换在Matlab中的实现,输出大小
Implementation of Radon transform in Matlab, output size
由于我的问题的性质,我想评估 Matlab 中 Radon 变换的数值实现(即不同的插值方法给出不同的数值)。
在尝试编写我自己的 Radon 代码并将其与 Matlab 的输出进行比较时,我发现我的 Radon 投影大小与 Matlab 的不同。
如果需要氡样本,我如何计算数量的一些直觉。让我们做二维案例吧。
想法是当对角线(至少是矩形)部分在 radon 变换中投影时最大尺寸,所以 diago=sqrt(size(I,1),size(I,2))
。因为我们什么都不想,n_r=ceil(diago)
。 n_r
应该是 radon 变换的离散样本量,以确保没有数据被遗漏。
我注意到 Matlab 的 radon
输出总是偶数,这是有道理的,因为您希望 "ray" 始终通过旋转中心。我注意到在所有情况下数组的端点都有 2 个零。
所以在那种情况下,n_r=ceil(diago)+mod(ceil(diago)+1,2)+2;
不过,似乎我与 Matlab 之间存在一些小差异。
MWE:
% Try: 255,256
pixels=256;
I=phantom('Modified Shepp-Logan',pixels);
rd=radon(I,pi/4);
size(rd,1)
s=size(I);
diagsize=sqrt(sum(s.^2));
n_r=ceil(diagsize)+mod(ceil(diagsize)+1,2)+2
rd=
367
n_r =
365
由于Matlab的Radon变换是一个我无法研究的函数,我想知道为什么会出现这种差异。
这是一个相当专业的问题,所以我会在不完全确定它是否是您特定问题的答案的情况下提供一个想法(通常我会通过并让其他人回答,但我不确定如何许多Whosebug的读者都研究过radon
)。我认为您可能忽略了 radon
函数调用文档中的 floor
函数。来自文档:
The radial coordinates returned in xp are the values along the x'-axis, which is
oriented at theta degrees counterclockwise from the x-axis. The origin of both
axes is the center pixel of the image, which is defined as
floor((size(I)+1)/2)
For example, in a 20-by-30 image, the center pixel is (10,15).
对于您传入的奇数或偶数大小的问题,这会给出不同的行为。因此,在您的示例 ("Try: 255, 256") 中,您需要针对奇数与偶数的不同情况,这可能涉及(实际上)用零行和零列填充。
我重新审视了这个问题,我相信这实际上是正确的答案。来自 radon.m 的 "hidden documentation"(输入 edit radon.m
并滚动到底部)
Grandfathered syntax
R = RADON(I,THETA,N) returns a Radon transform with the
projection computed at N points. R has N rows. If you do not
specify N, the number of points the projection is computed at
is:
2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
This number is sufficient to compute the projection at unit
intervals, even along the diagonal.
我没有尝试重新推导这个公式,但我想这就是你要找的。
由于我的问题的性质,我想评估 Matlab 中 Radon 变换的数值实现(即不同的插值方法给出不同的数值)。
在尝试编写我自己的 Radon 代码并将其与 Matlab 的输出进行比较时,我发现我的 Radon 投影大小与 Matlab 的不同。
如果需要氡样本,我如何计算数量的一些直觉。让我们做二维案例吧。
想法是当对角线(至少是矩形)部分在 radon 变换中投影时最大尺寸,所以 diago=sqrt(size(I,1),size(I,2))
。因为我们什么都不想,n_r=ceil(diago)
。 n_r
应该是 radon 变换的离散样本量,以确保没有数据被遗漏。
我注意到 Matlab 的 radon
输出总是偶数,这是有道理的,因为您希望 "ray" 始终通过旋转中心。我注意到在所有情况下数组的端点都有 2 个零。
所以在那种情况下,n_r=ceil(diago)+mod(ceil(diago)+1,2)+2;
不过,似乎我与 Matlab 之间存在一些小差异。
MWE:
% Try: 255,256
pixels=256;
I=phantom('Modified Shepp-Logan',pixels);
rd=radon(I,pi/4);
size(rd,1)
s=size(I);
diagsize=sqrt(sum(s.^2));
n_r=ceil(diagsize)+mod(ceil(diagsize)+1,2)+2
rd=
367
n_r =
365
由于Matlab的Radon变换是一个我无法研究的函数,我想知道为什么会出现这种差异。
这是一个相当专业的问题,所以我会在不完全确定它是否是您特定问题的答案的情况下提供一个想法(通常我会通过并让其他人回答,但我不确定如何许多Whosebug的读者都研究过radon
)。我认为您可能忽略了 radon
函数调用文档中的 floor
函数。来自文档:
The radial coordinates returned in xp are the values along the x'-axis, which is oriented at theta degrees counterclockwise from the x-axis. The origin of both axes is the center pixel of the image, which is defined as
floor((size(I)+1)/2)
For example, in a 20-by-30 image, the center pixel is (10,15).
对于您传入的奇数或偶数大小的问题,这会给出不同的行为。因此,在您的示例 ("Try: 255, 256") 中,您需要针对奇数与偶数的不同情况,这可能涉及(实际上)用零行和零列填充。
我重新审视了这个问题,我相信这实际上是正确的答案。来自 radon.m 的 "hidden documentation"(输入 edit radon.m
并滚动到底部)
Grandfathered syntax
R = RADON(I,THETA,N) returns a Radon transform with the projection computed at N points. R has N rows. If you do not specify N, the number of points the projection is computed at is:
2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
This number is sufficient to compute the projection at unit intervals, even along the diagonal.
我没有尝试重新推导这个公式,但我想这就是你要找的。