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

Question

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

找到图像（灰色，范围 0-255）的 2D 离散傅立叶变换 (DFT = F(u,v))
求傅立叶谱 (FS = abs(F(u,v))) 和功率谱 (PS = FS^2)
将光谱转换为极坐标并将其划分为 1º 间隔。
计算每个间隔（theta）的数均线强度（FI），即所有强度（像素）形成"theta"度的平均值到水平轴。
将 FI(theta) 变换为笛卡尔坐标

Cxy(theta) = [FI*cos(theta), FI*sin(theta)]
求矩阵Cxy'*Cxy

lambda1

lambda2

查找比对索引为 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。任何帮助将不胜感激。

PD：中间图中的那些黑点就是improfile使用的点。

Answer 1

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

离散傅立叶变换

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

合并与沿线采样

您完成分箱的方式是您从字面上画一条线（旋转特定角度）并使用 improfile 沿该线对图像进行采样。使用 improfile 沿该行执行数据插值，引入另一个潜在的错误源。默认值是最近邻插值，在下面显示的示例中可以导致多个 "profiles" 全部拾取相同的点。

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

数据欠采样

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

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

离图像中心越远，问题就越严重，因为 "bin" 中的数据实际上是扇形楔形，而您正在用一条线对其进行近似。

一个潜在的解决方案

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

这有几个优点：

它消除了我们谈到的欠采样，并从整个 "pie wedge" 中抽取样本，而不考虑采样角度。
确保每个像素属于一个且仅属于一个 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;

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

Fourier transform for fiber alignment

matlab

fft

image-processing