MATLAB 中的自适应椭圆结构元素
adaptive elliptical structuring element in MATLAB
我正在尝试为图像创建一个自适应椭圆结构元素以扩大或侵蚀它。我写了这段代码,但不幸的是所有的结构元素都是 ones(2*M+1)
.
I = input('Enter the input image: ');
M = input('Enter the maximum allowed semi-major axes length: ');
% determining ellipse parameteres from eigen value decomposition of LST
row = size(I,1);
col = size(I,2);
SE = cell(row,col);
padI = padarray(I,[M M],'replicate','both');
padrow = size(padI,1);
padcol = size(padI,2);
for m = M+1:padrow-M
for n = M+1:padcol-M
a = (l2(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
b = (l1(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
if e1(m-M,n-M,1)==0
phi = pi/2;
else
phi = atan(e1(m-M,n-M,2)/e1(m-M,n-M,1));
end
% defining structuring element for each pixel of image
x0 = m;
y0 = n;
se = zeros(2*M+1);
row_se = 0;
for i = x0-M:x0+M
row_se = row_se+1;
col_se = 0;
for j = y0-M:y0+M
col_se = col_se+1;
x = j-y0;
y = x0-i;
if ((x*cos(phi)+y*sin(phi))^2)/a^2+((x*sin(phi)-y*cos(phi))^2)/b^2 <= 1
se(row_se,col_se) = 1;
end
end
end
SE{m-M,n-M} = se;
end
end
a
、b
和 phi
是半长轴和半短轴长度,phi 是 a
和 x 轴之间的角度。
我使用了 2 个 MATLAB 函数来计算图像的局部结构张量,然后是每个像素的特征值和特征向量。这些是矩阵 l1
、l2
、e1
和 e2
.
这是您的代码中我不理解的部分:
a = (l2(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
b = (l1(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
我将 b
的表达式简化为(仅删除索引):
b = (l1+eps/l1+l2+2*eps)*M;
对于正常范围内的 l1
和 l2
,我们得到:
b =(approx)= (l1+0/l1+l2+2*0)*M = (l1+l2)*M;
因此,b
可以很容易地大于 M
,我认为这不是您的意图。 eps
在这种情况下也不能防止被零除,这通常是添加 eps
的目的:如果 l1
为零,则 eps/l1
为 Inf
.
看这个表达式,在我看来你的本意是这样的:
b = (l1+eps)/(l1+l2+2*eps)*M;
在这里,您将 eps
添加到每个特征值,使它们保证非零(结构张量是对称的,半正定的)。然后将 l1
除以特征值之和,然后乘以 M
,这会导致每个轴的值介于 0
和 M
之间。
所以,这似乎是括号放错地方的情况。
仅作记录,这是您在代码中需要的内容:
a = (l2(m-M,n-M)+eps ) / ( l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
b = (l1(m-M,n-M)+eps ) / ( l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
^ ^
added parentheses
请注意,您可以通过在循环外定义来简化代码:
[se_x,se_y] = meshgrid(-M:M,-M:M);
i
和 j
的内部两个循环构造 se
可以简单地写成:
se = ((se_x.*cos(phi)+se_y.*sin(phi)).^2)./a.^2 + ...
((se_x.*sin(phi)-se_y.*cos(phi)).^2)./b.^2 <= 1;
(请注意 .*
和 .^
运算符,它们执行逐元素乘法和幂运算。)
进一步的细微改进来自于认识到 phi
首先从 e1(m,n,1)
和 e1(m,n,2)
计算,然后用于对 cos
和 [=43= 的调用].如果我们假设特征向量被正确归一化,那么
cos(phi) == e1(m,n,1)
sin(phi) == e1(m,n,2)
但您始终可以确保它们已标准化:
cos_phi = e1(m-M,n-M,1);
sin_phi = e1(m-M,n-M,2);
len = hypot(cos_phi,sin_phi);
cos_phi = cos_phi / len;
sin_phi = sin_phi / len;
se = ((se_x.*cos_phi+se_y.*sin_phi).^2)./a.^2 + ...
((se_x.*sin_phi-se_y.*cos_phi).^2)./b.^2 <= 1;
考虑到三角函数运算相当昂贵,这应该会稍微加快您的代码速度。
我正在尝试为图像创建一个自适应椭圆结构元素以扩大或侵蚀它。我写了这段代码,但不幸的是所有的结构元素都是 ones(2*M+1)
.
I = input('Enter the input image: ');
M = input('Enter the maximum allowed semi-major axes length: ');
% determining ellipse parameteres from eigen value decomposition of LST
row = size(I,1);
col = size(I,2);
SE = cell(row,col);
padI = padarray(I,[M M],'replicate','both');
padrow = size(padI,1);
padcol = size(padI,2);
for m = M+1:padrow-M
for n = M+1:padcol-M
a = (l2(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
b = (l1(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
if e1(m-M,n-M,1)==0
phi = pi/2;
else
phi = atan(e1(m-M,n-M,2)/e1(m-M,n-M,1));
end
% defining structuring element for each pixel of image
x0 = m;
y0 = n;
se = zeros(2*M+1);
row_se = 0;
for i = x0-M:x0+M
row_se = row_se+1;
col_se = 0;
for j = y0-M:y0+M
col_se = col_se+1;
x = j-y0;
y = x0-i;
if ((x*cos(phi)+y*sin(phi))^2)/a^2+((x*sin(phi)-y*cos(phi))^2)/b^2 <= 1
se(row_se,col_se) = 1;
end
end
end
SE{m-M,n-M} = se;
end
end
a
、b
和 phi
是半长轴和半短轴长度,phi 是 a
和 x 轴之间的角度。
我使用了 2 个 MATLAB 函数来计算图像的局部结构张量,然后是每个像素的特征值和特征向量。这些是矩阵 l1
、l2
、e1
和 e2
.
这是您的代码中我不理解的部分:
a = (l2(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M; b = (l1(m-M,n-M)+eps/l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
我将 b
的表达式简化为(仅删除索引):
b = (l1+eps/l1+l2+2*eps)*M;
对于正常范围内的 l1
和 l2
,我们得到:
b =(approx)= (l1+0/l1+l2+2*0)*M = (l1+l2)*M;
因此,b
可以很容易地大于 M
,我认为这不是您的意图。 eps
在这种情况下也不能防止被零除,这通常是添加 eps
的目的:如果 l1
为零,则 eps/l1
为 Inf
.
看这个表达式,在我看来你的本意是这样的:
b = (l1+eps)/(l1+l2+2*eps)*M;
在这里,您将 eps
添加到每个特征值,使它们保证非零(结构张量是对称的,半正定的)。然后将 l1
除以特征值之和,然后乘以 M
,这会导致每个轴的值介于 0
和 M
之间。
所以,这似乎是括号放错地方的情况。
仅作记录,这是您在代码中需要的内容:
a = (l2(m-M,n-M)+eps ) / ( l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
b = (l1(m-M,n-M)+eps ) / ( l1(m-M,n-M)+l2(m-M,n-M)+2*eps)*M;
^ ^
added parentheses
请注意,您可以通过在循环外定义来简化代码:
[se_x,se_y] = meshgrid(-M:M,-M:M);
i
和 j
的内部两个循环构造 se
可以简单地写成:
se = ((se_x.*cos(phi)+se_y.*sin(phi)).^2)./a.^2 + ...
((se_x.*sin(phi)-se_y.*cos(phi)).^2)./b.^2 <= 1;
(请注意 .*
和 .^
运算符,它们执行逐元素乘法和幂运算。)
进一步的细微改进来自于认识到 phi
首先从 e1(m,n,1)
和 e1(m,n,2)
计算,然后用于对 cos
和 [=43= 的调用].如果我们假设特征向量被正确归一化,那么
cos(phi) == e1(m,n,1)
sin(phi) == e1(m,n,2)
但您始终可以确保它们已标准化:
cos_phi = e1(m-M,n-M,1);
sin_phi = e1(m-M,n-M,2);
len = hypot(cos_phi,sin_phi);
cos_phi = cos_phi / len;
sin_phi = sin_phi / len;
se = ((se_x.*cos_phi+se_y.*sin_phi).^2)./a.^2 + ...
((se_x.*sin_phi-se_y.*cos_phi).^2)./b.^2 <= 1;
考虑到三角函数运算相当昂贵,这应该会稍微加快您的代码速度。