光纤对准的傅里叶变换

Question

光纤对准的傅里叶变换

3

我正在开发一个应用程序，用于从图像中确定纤维网络的对齐程度。我已经阅读了几篇关于这个问题的论文，它们基本上是这样做的：

1. 找到图像（灰度，范围为0-255）的2D离散傅里叶变换（DFT = F（u，v））。 2. 找到傅里叶频谱（FS = abs（F（u，v）））和功率谱（PS = FS ^ 2）。 3. 将频谱转换为极坐标并将其分成1º间隔。 4. 计算每个间隔（theta）的数平均线强度（FI），即所有与水平轴形成“theta”度的强度（像素）的平均值。 5. 将FI（theta）转换为笛卡尔坐标 Cxy（theta）= [FI * cos（theta），FI * sin（theta)] 6. 找到矩阵Cxy'* Cxy的特征值（lambda1和lambda2）。 7. 将对齐指数作为alpha = 1-lamda2 / lambda1找到。

我已经在MATLAB中实现了这个过程（下面是代码），但是我不确定是否正确，因为第3点和第4点对我来说不是很清楚（虽然我得到了类似于论文的结果，但并非在所有情况下都是如此）。例如，在第3点中，“spectrum”是指FS还是PS？而在第4点中，应该如何进行平均？是否考虑所有像素？（即使对角线上有更多像素）。

rgb = imread('network.tif');%513x513 pixels
im = rgb2gray(rgb);
im = imrotate(im,-90);%since FFT space is rotated 90º
FT = fft2(im) ;
FS = abs(FT); %Fourier spectrum
PS = FS.^2; % Power spectrum
FS = fftshift(FS);
PS = fftshift(PS);

xoffset = (513-1)/2;
yoffset = (513-1)/2;

% Avoid low frequency points
x1 = 5;
y1 = 0;

% Maximum high frequency pixels
x2 = 255;
y2 = 0;

for theta = 0:pi/180:pi
    % Transposed rotation matrix
    Rt = [cos(theta) sin(theta); 
         -sin(theta) cos(theta)]; 

    % Find radial lines necessary for improfile
    xy1_rot = Rt * [x1; y1] + [xoffset; yoffset]; 
    xy2_rot = Rt * [x2; y2] + [xoffset; yoffset];

    plot([xy1_rot(1) xy2_rot(1)], ...
         [xy1_rot(2) xy2_rot(2)], ...
         'linestyle','none', ...
         'marker','o', ...
         'color','k');

     prof = improfile(F,[xy1_rot(1) xy2_rot(1)],[xy1_rot(2) xy2_rot(2)]);
     i = i + 1;
     FI(i) = sum(prof(:))/length(prof);
     Cxy(i,:) = [FI(i)*cos(theta), FI(i)*sin(theta)];
end

C = Cxy'*Cxy;
[V,D] = eig(C)
lambda2 = D(1,1);
lambda1 = D(2,2);

alpha = 1 - lambda2/lambda1

图示：A）原始图像，B）log（P + 1）绘图，C）FI极坐标图。

我的主要问题是，当我选择完美对齐的人工图像（如附图），我得到alpha = 0.91，而实际上应该是准确的1。非常感谢任何帮助。

注：中间图中的黑点只是被improfile使用的点。

- Carlos Borau

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Suever · Accepted Answer

我认为这里可能存在几个潜在的误差来源，导致您无法得到完美的α值。

离散傅里叶变换

您有离散成像数据，这强制您进行离散傅里叶变换，这不可避免地（取决于输入数据的分辨率）会存在一些精度问题。

分箱与沿线采样

您所做的分箱方式是，实际上画了一条直线（以特定角度旋转），并使用improfile沿着该线对图像进行采样。使用improfile会在该线上对数据进行插值，引入了另一个潜在的误差来源。默认情况下是最近邻插值，如下面示例所示，可能会导致多个“轮廓”都选择相同的点。

这是在垂直线偏离1度的情况下进行的旋转，而从技术上讲，您希望这些峰值仅出现在完全垂直的线上。很明显，这种傅里叶谱的插值方式会导致“正确”答案周围的扩散。 数据欠采样 与傅里叶域中的奈奎斯特采样类似，空间域中的采样也有一些要求。

想象一下，如果你想使用45度的条宽而不是1度，你的方法仍然会沿着一条细线采样，并使用该样本来代表45度的数据。显然，这是对数据的严重欠采样，可以想象结果不会非常准确。

随着你离图像中心越远，问题就越来越严重，因为这个“bin”中的数据实际上是馅饼楔形的，而你正在用一条线来近似它。 一个潜在的解决方案 一种不同的分箱方法是确定图像中所有像素中心的极坐标（r、theta）。然后将theta组件分成1度的箱子。然后将落入该箱子中的所有值相加。

这有几个优点：

它消除了我们谈论过的欠采样，并从整个“饼图楔形”中绘制样本，而不考虑采样角度。
它确保每个像素属于一个且仅属于一个角度区间。

我已经在下面的代码中实现了这种替代方法，并使用一些错误的水平线数据实现了0.988的alpha值，考虑到数据的离散性，我认为这是相当不错的。

% Draw a bunch of horizontal lines
data = zeros(101);
data([5:5:end],:) = 1;

fourier = fftshift(fft2(data));

FS = abs(fourier);
PS = FS.^2;

center = fliplr(size(FS)) / 2;

[xx,yy] = meshgrid(1:size(FS,2), 1:size(FS, 1));

coords = [xx(:), yy(:)];

% De-mean coordinates to center at the middle of the image
coords = bsxfun(@minus, coords, center);

[theta, R] = cart2pol(coords(:,1), coords(:,2));

% Convert to degrees and round them to the nearest degree
degrees = mod(round(rad2deg(theta)), 360);

degreeRange = 0:359;

% Band pass to ignore high and low frequency components;
lowfreq = 5;
highfreq = size(FS,1)/2;

% Now average everything with the same degrees (sum over PS and average by the number of pixels)
for k = degreeRange
    ps_integral(k+1) = mean(PS(degrees == k & R > lowfreq & R < highfreq));
    fs_integral(k+1) = mean(FS(degrees == k & R > lowfreq & R < highfreq));
end

thetas = deg2rad(degreeRange);

Cxy = [ps_integral.*cos(thetas);
       ps_integral.*sin(thetas)]';

C = Cxy' * Cxy;
[V,D] = eig(C);

lambda2 = D(1,1);
lambda1 = D(2,2);

alpha = 1 - lambda2/lambda1;