matlab - 用于光纤对准的傅里叶变换

Question

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

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 点中，“频谱”指的是 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。 任何帮助将不胜感激。

PD:中间图中的那些黑点只是improfile使用的点。

Answer 1

我相信这里有几个潜在的错误来源会导致您无法获得完美的 alpha 值。

离散傅立叶变换

您的离散成像数据迫使您进行离散傅里叶变换，这不可避免地(取决于输入数据的分辨率)会出现一些准确性问题。

分箱与沿线采样

您完成分箱的方式是您实际上画了一条线(旋转了特定角度)并使用 improfile 沿该线对图像进行采样。使用 improfile 沿该行执行数据插值，引入另一个潜在的错误源。默认值是最近邻插值，在下面显示的示例中，它会导致多个“配置文件”全部拾取相同的点。

这是在偏离垂直方向旋转 1 度时，从技术上讲，您希望这些峰值仅出现在一条完全垂直的线上。很明显，这种傅立叶频谱插值如何导致围绕“正确”答案展开。

数据欠采样

与傅里叶域的奈奎斯特采样类似，空间域的采样也有一些要求。

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

离图像中心越远，问题就越严重，因为这个“bin”中的数据实际上是饼楔形的，你用一条线来近似它。

一个潜在的解决方案

合并的另一种方法是确定图像中所有像素中心的极坐标 (r, theta)。然后将 theta 分量分箱到 1 度箱中。然后将所有落入该 bin 的值相加。

这有几个优点:

• 它消除了我们谈到的欠采样，并从整个“饼楔”中抽取样本，而不管采样角度如何。
• 确保每个像素属于一个且仅属于一个 angular bin

我已经在下面的代码中使用一些错误的水平线数据实现了这种替代方法，并且能够实现 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;