傅里叶变换用于光纤对准
我正在研究一个应用程序，以根据图像确定光纤网络的对齐程度。 我已经阅读了关于这个问题的几篇论文，他们基本上这样做：
DFT = F(u,v) ）， FS = abs(F(u,v)) ）和功率谱（ PS = FS^2 ） theta ）的数均线强度（ FI ），即相对于水平轴形成“θ”度的所有强度（像素）的平均值。 将FI（theta）转换为笛卡尔坐标
Cxy(theta) = [FI*cos(theta), FI*sin(theta)]
查找特征值（ lambda1和lambda2矩阵） Cxy'*Cxy
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：中间图中的那些黑点只是由简档所使用的点。
我相信这里有一些潜在的错误来源会导致你没有获得完美的alpha值。
离散傅里叶变换
你有离散的成像数据，迫使你采取离散傅里叶变换，这不可避免地（取决于输入数据的分辨率）有一些准确性问题。
分档与沿线取样
你完成分箱的方式是你从字面上画了一条线（旋转一个特定的角度），并使用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;