adaptive elliptical structuring element in MATLAB

This is the bit of your code I didn't understand:

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;

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

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

For l1 and l2 in the normal range we get:

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

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;

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);

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.)

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.