MATLAB中的自适应椭圆结构元素 [英] adaptive elliptical structuring element in MATLAB

问题描述

我正在尝试为图像创建自适应的椭圆结构元素以使其扩张或腐蚀.我编写了这段代码，但不幸的是，所有结构元素都是ones(2*M+1).

I'm trying to create an adaptive elliptical structuring element for an image to dilate or erode it. I write this code but unfortunately all of the structuring elements are 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轴之间的角度.

a, b and phi are semi-major and semi-minor axes length and phi is angle between a and x axis.

我使用了2个MATLAB函数来计算图像的局部结构张量，然后计算每个像素的特征值和特征向量.这些是矩阵l1，l2，e1和e2.

I used 2 MATLAB functions to compute the Local Structure Tensor of the image, and then its eigenvalues and eigenvectors for each pixel. These are the matrices l1, l2, e1 and 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的表达式简化为(只是删除了索引):

I simplified the expression for b to (just removing the indexing):

b = (l1+eps/l1+l2+2*eps)*M;

对于l1和l2在正常范围内，我们得到:

For l1 and l2 in the normal range we get:

b =(approx)= (l1+0/l1+l2+2*0)*M = (l1+l2)*M;

因此，b可以容易地大于M，我认为这不是您的意图.在这种情况下，eps也不能防止被零除，这通常是添加eps的目的:如果l1为零，则eps/l1为Inf.

Thus, b can easily be larger than M, which I don't think is your intention. The eps in this case also doesn't protect against division by zero, which is typically the purpose of adding eps: if l1 is zero, eps/l1 is Inf.

看着这个表达式，在我看来您原本打算这样做的:

Looking at this expression, it seems to me that you intended this instead:

b = (l1+eps)/(l1+l2+2*eps)*M;

在这里，您要向每个特征值添加eps，以确保它们不为零(结构张量是对称的，正半定数).然后，将l1除以特征值之和，再乘以M，则得出每个轴的0和M之间的值.

Here, you're adding eps to each of the eigenvalues, making them guaranteed non-zero (the structure tensor is symmetric, positive semi-definite). Then you're dividing l1 by the sum of eigenvalues, and multiplying by M, which leads to a value between 0 and M for each of the axes.

因此，这似乎是括号放置错误的情况.

So, this seems to be a case of misplaced parenthesis.

仅作记录，这就是您在代码中所需要的:

Just for the record, this is what you need in your code:

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

---

请注意，您可以通过在循环之外进行定义来简化代码:

---

Note that you can simplify your code by defining, outside of the loops:

[se_x,se_y] = meshgrid(-M:M,-M:M);

构造se的内部两个循环(在i和j上)可以简单地写为:

The inner two loops, over i and j, to construct se can then be written simply as:

se = ((se_x.*cos(phi)+se_y.*sin(phi)).^2)./a.^2 + ...
     ((se_x.*sin(phi)-se_y.*cos(phi)).^2)./b.^2 <= 1;

(请注意.*和.^运算符，它们执行逐元素乘法和幂运算.)

(Note the .* and .^ operators, these do element-wise multiplication and power.)

认识到首先从e1(m,n,1)和e1(m,n,2)计算出phi，然后在对cos和sin的调用中使用phi，这带来了进一步的细微改进.如果我们假设特征向量已正确归一化，那么

A further slight improvement comes from realizing that phi is first computed from e1(m,n,1) and e1(m,n,2), and then used in calls to cos and sin. If we assume that the eigenvector is properly normalized, then

cos(phi) == e1(m,n,1)
sin(phi) == e1(m,n,2)

但是您始终可以确保将它们标准化:

But you can always make sure they are normalized:

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;

考虑三角运算是相当昂贵的，这应该会加快代码的速度.

Considering trigonometric operations are fairly expensive, this should speed up your code a bit.

这篇关于MATLAB中的自适应椭圆结构元素的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！