matlab - matlab中结构张量的特征值分解
问题描述
我有一个合成图像。为了某些边缘检测目的,我想对其局部结构张量(LST)进行特征值分解。我使用LST的eigenvaluesl1
和l2
eigenvectors为图像的每个像素生成一个自适应椭圆。不幸的是,我得到不相等的 eigenvalues ,因此对于我的图形的均匀区域,椭圆的半轴长度不相等:
e1
e2
l1
l2
我不知道我的代码有什么问题:
function [H,e1,e2,l1,l2] = LST_eig(I,sigma1,rw)
% LST_eig - compute the structure tensor and its eigen
% value decomposition
%
% H = LST_eig(I,sigma1,rw);
%
% sigma1 is pre smoothing width (in pixels).
% rw is filter bandwidth radius for tensor smoothing (in pixels).
%
n = size(I,1);
m = size(I,2);
if nargin<2
sigma1 = 0.5;
end
if nargin<3
rw = 0.001;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pre smoothing
J = imgaussfilt(I,sigma1);
% compute gradient using Sobel operator
Sch = [-3 0 3;-10 0 10;-3 0 3];
%h = fspecial('sobel');
gx = imfilter(J,Sch,'replicate');
gy = imfilter(J,Sch','replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute tensors
gx2 = gx.^2;
gy2 = gy.^2;
gxy = gx.*gy;
% smooth
gx2_sm = imgaussfilt(gx2,rw); %rw/sqrt(2*log(2))
gy2_sm = imgaussfilt(gy2,rw);
gxy_sm = imgaussfilt(gxy,rw);
H = zeros(n,m,2,2);
H(:,:,1,1) = gx2_sm;
H(:,:,2,2) = gy2_sm;
H(:,:,1,2) = gxy_sm;
H(:,:,2,1) = gxy_sm;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% eigen decomposition
l1 = zeros(n,m);
l2 = zeros(n,m);
e1 = zeros(n,m,2);
e2 = zeros(n,m,2);
for i = 1:n
for j = 1:m
Hmat = zeros(2);
Hmat(:,:) = H(i,j,:,:);
[V,D] = eigs(Hmat);
D = abs(D);
l1(i,j) = D(1,1); % eigen values
l2(i,j) = D(2,2);
e1(i,j,:) = V(:,1); % eigen vectors
e2(i,j,:) = V(:,2);
end
end
任何帮助表示赞赏。
这是我的椭圆绘图代码:
% determining ellipse parameteres from eigen value decomposition of LST
M = input('Enter the maximum allowed semi-major axes length: ');
I = input('Enter the input data: ');
row = size(I,1);
col = size(I,2);
a = zeros(row,col);
b = zeros(row,col);
cos_phi = zeros(row,col);
sin_phi = zeros(row,col);
for m = 1:row
for n = 1:col
a(m,n) = (l2(m,n)+eps)/(l1(m,n)+l2(m,n)+2*eps)*M;
b(m,n) = (l1(m,n)+eps)/(l1(m,n)+l2(m,n)+2*eps)*M;
cos_phi1 = e1(m,n,1);
sin_phi1 = e1(m,n,2);
len = hypot(cos_phi1,sin_phi1);
cos_phi(m,n) = cos_phi1/len;
sin_phi(m,n) = sin_phi1/len;
end
end
%% plot elliptic structuring elements using parametric equation and superimpose on the image
figure; imagesc(I); colorbar; hold on
t = linspace(0,2*pi,50);
for i = 10:10:row-10
for j = 10:10:col-10
x0 = j;
y0 = i;
x = a(i,j)/2*cos(t)*cos_phi(i,j)-b(i,j)/2*sin(t)*sin_phi(i,j)+x0;
y = a(i,j)/2*cos(t)*sin_phi(i,j)+b(i,j)/2*sin(t)*cos_phi(i,j)+y0;
plot(x,y,'r','linewidth',1);
hold on
end
end
这是新的情节axis equal
:
解决方案
我创建了一个类似于您的测试图像(可能不太复杂),如下所示:
pos = yy([400,500]) + 100 * sin(xx(400)/400*2*pi);
img = gaussianlineclip(pos+50,7) + gaussianlineclip(pos-50,7);
I = double(stretch(img));
(这需要DIPimage运行)
然后运行你LST_eig
的它(sigma1=1
和rw=3
)和你的代码来绘制椭圆(没有改变,除了添加axis equal
),并得到这个结果:
我怀疑图像的某些蓝色区域存在一些不均匀性,这会导致出现非常小的渐变。使用椭圆时定义椭圆的问题在于,对于足够定向的模式,即使该模式难以察觉,您也会得到一条线。您可以通过如下定义椭圆轴长度来解决此问题:
a = repmat(M,size(l2)); % longest axis is always the same
b = M ./ (l2+1); % shortest axis is shorter the more important the largest eigenvalue is
最小的特征值l1
在具有强梯度但没有明确方向的区域中很高。上面没有考虑到这一点。一种选择可能是a
依赖于能量和各向异性测量,并且b
仅依赖于能量:
T = 1000; % some threshold
r = M ./ max(l1+l2-T,1); % circle radius, smaller for higher energy
d = (l2-l1) ./ (l1+l2+eps); % anisotropy measure in range [0,1]
a = M*d + r.*(1-d); % use `M` length for high anisotropy, use `r` length for high isotropy (circle)
b = r; % use `r` width always
这样,如果有强梯度但没有明确的方向,则整个椭圆会缩小,而当只有弱梯度或没有梯度时,它会保持大而圆形。阈值T
取决于图像强度,根据需要进行调整。
您可能还应该考虑取特征值的平方根,因为它们对应于方差。
一些建议:
你可以写
a = (l2+eps)./(l1+l2+2*eps) * M; b = (l1+eps)./(l1+l2+2*eps) * M; cos_phi = e1(:,:,1); sin_phi = e1(:,:,2);
没有循环。注意
e1
是按定义归一化的,不需要再次归一化。使用高斯梯度而不是高斯平滑,然后使用 Sobel 或 Schaar 滤波器。有关一些 MATLAB 实现细节,请参见此处。
当您需要所有特征值时,请使用
eig
,而不是。eigs
尤其是对于这么小的矩阵,使用 . 没有任何优势eigs
。eig
似乎会产生更一致的结果。没有必要取特征值 (D = abs(D)
) 的绝对值,因为它们根据定义是非负的。您的默认值
rw = 0.001
太小,该大小的 sigma 对图像没有影响。这种平滑的目标是平均局部邻域中的梯度。我用rw=3
了,效果很好。使用 DIP 图像。有一个
structuretensor
函数、高斯梯度和更多有用的东西。3.0 版本(仍在开发中)是一个重大改写,在处理矢量和矩阵值图像方面有显着改进。我可以将您的所有内容LST_eig
写成如下:I = dip_image(I); g = gradient(I, sigma1); H = gaussf(g*g.', rw); [e,l] = eig(H); % Equivalences with your outputs: l1 = l{2}; l2 = l{1}; e1 = e{2,:}; e2 = e{1,:};
推荐阅读
- python - 有没有一种方法可以限制 scipy 生成随机变量集?
- azure-devops - 如何为测试和发布环境构建 Azure Devops Pipelines?
- javascript - 如何从字符串中拆分某个数字但在该字符串中保留另一个数字仍在javascript中
- c# - 是什么导致视频播放期间“尚未为此流或其依赖的流设置有效类型”异常
- python-3.x - 计算会影响 python-3.x 中其他地方的堆栈中相同(但不同)的变量?
- react-native - TypeError:无法读取 /node_modules/galio-framework/src/Toast.js 中未定义的属性“样式”
- java - Lombok 条件 ToString.Include
- three.js - threejs - 重新定位后使用控制器在 AR 中进行光线投射
- java - 使用 java 和 jsoup 进行网页抓取
- node.js - Socket.io 在线用户数