Alternative function to fit ellipses to grains

[kilir]

Hi Ralf, I was trying out some ellispe fit methods, mainly because principalComponents is sometimes a bit exaggerating when it comes to aspect ratios (see attached figures). Below I used the method described in Mulchrone& Choudhury,2004 (https://doi.org/10.1016/S0191-8141(03)00093-2) - although it looks quite familiar to what happens in principalComponents, - the result is for my "feeling" for aspect ratio more satisfying.

That's what principalComponents does (colorcbar gives a/b, aspect ratio):

and this is what I can come up with below:

function [omega,a,b] = fitEllipse(grains,varargin)
% fit an ellipse to a grains using moments
% 
% Input:
%  grains   - @grain2d
%
% Output:
%  omega    - angle  giving trend of ellipse long axis
%  a        - long axis radius
%  b        - short axis radius
% 
% Options:
%  hull     - use convex hull tu fit ellipse
%

Vin=grains.V; % sort of carries on ALL! V

omega=zeros(length(grains),1);
a=omega;
b=omega;

% todo: get rid of that loop
%       scale to boundary length

for i = 1:length(grains)

V= Vin(grains(i).idV,:); % vertices of ith grains

if check_option(varargin,'hull')
    %use hull
    ind = convhulln(V);
    V = V(ind(:,1),:);
end

n=size(V,1);

% get centers and 2nd moments of area
xb= sum(V(:,1))/n;
yb= sum(V(:,2))/n;

u20=sum((V(:,1)-xb).^2)/n;
u02=sum((V(:,2)-yb).^2)/n;
u11=abs(sum( (V(:,1)-xb).*(V(:,2)-yb) )/n); % since this can be negative
                                            % but we like real angles

% directions
omega(i)= 0.5*atan(2*u11/(u20-u02));
% axes; - eigenvalues of [u20 u11; u11 u02]
a(i)=sqrt(2*( u20+u02+sqrt(4*u11^2+(u20-u02).^2 ))/u11);
b(i)=sqrt(2*( u20+u02-sqrt(4*u11^2+(u20-u02).^2 ))/u11);

end
% scale to area
areaElli= pi.*a.*b;
scalef= sqrt(grains.area./areaElli);

a = a.*scalef;
b = b.*scalef;

end

Currently it's slower than principalComponents, since I did not manage to get rid of the loop, so any hints on how to do that are welcome.

Cheers, Rüdiger

[marcoalopez]

Hello Rüdiger,

as an alternative to the calculation of grain aspect ratios based on fitting ellipses, perhaps you can use the ratio between the maximum and minimum feret (calliper) diameters. In principle, feret diameters algorithms do not suffer from these undesirable effects that you referred above. If you are interested, you can find an example of the implementation of the algorithm for the calculation of feret diameters in the following link: https://blogs.mathworks.com/steve/2017/09/29/feret-diameter-introduction/

hope this helps.

All the best, Marco

[kilir]

Hi Marco, thanks for the suggestion. My problem with calliper diameters is, that not necessarily the longest and the shortest may not be orthogonal and hence one could expect those to flucuate quite a bit on grain shape. One could think of using the longets and the length orthogonal to it - or reverse, taking the shortest caliper and the corresponding orthognal length. I think one might be seeing the issue when e.g. grain are preferentially rhomb shaped vs. grains which are rather elongated rectangles. Using the longest calliper on the latter case would preferentially find one of the diagonals which are systematically off the direction of the long axis of the grains. The nice thing about ellipses is that - if we have a good reason to believe so - they can realte to strain- However, having the directions of the min or max callipers or arbitrary callipare length into a given direction is something nice to have in general, so that's a good idea. Cheers, Rüdiger

[zmichels]

Hi all.

Just chiming in. First of all... Rüdiger, I think your attempt you showed is a notable improvement already over the existing Principal Components.

Maybe the following doesn’t make much sense and I am not sure how to do it offhand but: What if you used the calipers — although not necessarily orthogonal — to define the long axis... and use the short caliper to constrain the boundary of the ellipse? So that a fit ellipse boundary passes through the large caliper and small caliper points, and uses the large caliper to define the long axis...(?)

Probably that suggestion is not useful, but I am following this discussion. As I said before, I think what you’ve done is already a step towards something more usable. I look forward to whatever else you do to improve it, and thank you for tackling this.

Cheers, Zach

On Mar 16, 2018, at 12:38 PM, Rüdiger Kilian notifications@github.com<mailto:notifications@github.com> wrote:

Hi Marco, thanks for the suggestion. My problem with calliper diameters is, that not necessarily the longest and the shortest may not be orthogonal and hence one could expect those to flucuate quite a bit on grain shape. One could think of using the longets and the length orthogonal to it - or reverse, taking the shortest caliper and the corresponding orthognal length. I think one might be seeing the issue when e.g. grain are preferentially rhomb shaped vs. grains which are rather elongated rectangles. Using the longest calliper on the latter case would preferentially find one of the diagonals which are systematically off the direction of the long axis of the grains. The nice thing about ellipses is that - if we have a good reason to believe so - they can realte to strain- However, having the directions of the min or max callipers or arbitrary callipare length into a given direction is something nice to have in general, so that's a good idea. Cheers, Rüdiger

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/mtex-toolbox/mtex/issues/332#issuecomment-373789969, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AIXQUzCbjG3TFAfh_qF2NzgsuudCirVHks5te_iTgaJpZM4SsUdK.

[kilir]

The thing with the longest calliper is that it'll give you the diagonal in case you have sort of brick shaped grains, e.g. just think of idiomorphic feldspar laths (I agree this isn't somethign very common at the scale of ebsd maps, but still, possible to have these types of grain shapes).

Hoverever, having some sort of calliper functionality might comne handy as well - here's a bit of a result with some artifical grains:

principalComponents:

fitEllipse:

calliper (here the ellipse is paralle the longest calliper and the short axis is of length perpendicular to longest - and just not scaled to area):

Cheers, Rüdiger

[zmichels]

I understand what you mean and the examples show it well. The “fitEllipse” approach seems the best fit to me, and also the most conservative. Looks good!

Cheers, Zach

On Mar 18, 2018, at 3:57 PM, Rüdiger Kilian notifications@github.com<mailto:notifications@github.com> wrote:

The thing with the longest calliper is that it'll give you the diagonal in case you have sort of brick shaped grains, e.g. just think of idiomorphic feldspar laths (I agree this isn't somethign very common at the scale of ebsd maps, but still, possible to have these types of grain shapes).

Hoverever, having some sort of calliper functionality might comne handy as well - here's a bit of a result with some artifical grains:

principalComponents: [principal_components]https://user-images.githubusercontent.com/11697557/37570960-e1b55ff8-2af6-11e8-9c71-bf6d457d0929.png

fitEllipse: [moment_ellipses]https://user-images.githubusercontent.com/11697557/37570965-f3721326-2af6-11e8-9886-fb8ca2167fca.png

calliper (here the ellipse is paralle the longest calliper and the short axis is of length perpendicular to longest - and just not scaled to area):

[calliper_longest]https://user-images.githubusercontent.com/11697557/37570974-222c4182-2af7-11e8-97b6-7536b08393d0.png

Cheers, Rüdiger

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/mtex-toolbox/mtex/issues/332#issuecomment-374046434, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AIXQU0c3tdDN85RWVzx0UwZ0nsed_lGuks5tfsopgaJpZM4SsUdK.

[ralfHielscher]

Hi Rüdiger,

I have included the code into MTEX, so that we can have a base for further discussions:

1) For me it does not work in all cases. E.g.

mtexdata csl
grains = calcGrains(ebsd); grains = smooth(grains,5);

plot(grains)
[omega,a,b] = fitEllipse(grains);
plotEllipse(grains.centroid,a,b,omega,'lineColor','r')

results in some ellipse wrongly oriented. Especially the grain with id=676.

2. I very much like your artificial grains. May you include this script into MTEX? We may use it for testing or documenting the differences.

3. It would be nice to have some documentation page on ellipse fitting with focus what could be done further with this fit. I remember that the function export_VPSC should include the ellipse parameters.

Ralf.

[kilir]

Wrt 1: Ok, I think the wrong long axis should be fixed by baea4d846efafaf481abba616d22c7055e90fbdf I also added the calliper function used above since it might be usefull in some cases.

Wrt 2: It's sort of a very hackish approach I use to "import" manually derived grain maps and/or phase maps as a "fake" ebsd to use them with the grain functionalities. For now I've put it here (https://gist.github.com/kilir/3f731daf0cd42fea95e3b01e0fb71a68), please have a look (also if one could speed up the segmentation - which is obviously the required but non-sensical part ). I also attached some test images and results to the gist.

Wrt 3: I can try and work on that.

Cheers, Rüdiger

[ralfHielscher]

Hi Rüdiger,

somehow the alignment is still wrong for about 50% of the grains in the CSL example.

Ralf.

[kilir]

Now (a054992f2a49dd83dab350a40fa21990aebed73c) it should work, at least here's what I get:

mtexdata csl
grains = calcGrains(ebsd); grains = smooth(grains,5);
[om,a,b]=fitEllipse(grains);
plot(grains)
hold on
plotEllipse(grains.centroid,a,b,om,'linewidth',2,'lineColor','r')
hold off

[ralfHielscher]

This what I get: No clue why

[kilir]

Still with the latest version? I can't get it to misbehave in such a way now, but I needed to clear the toolbox path cache or just make sure somehow that the recent file is used.

[ralfHielscher]

Hi Rüdiger,

sorry, your latest version was correct. I took the chance and went a bit deeper into the code and recognized the following:

1) grain2d/principleComponents had a bug, replacing there the long and short axes of the ellipse with its square root results in almost the same result as fitEllipse

2) actually, both algorithm do compute almost the same, have a look at the rewritten code. The only difference is that the fitEllipse code does not take into account that the length of the segments may be not the same everywhere. I would expect that this code works worse if vertices are not uniformaly distributed around the grain

3) I included a function grain2d/hull which replaces all grains by its convex hull. Allowing to remove this option from the fitEllipse function

4) fitEllipse included code for vectorised computation of the eigen values and vectors of a symmetric 2x2 matrix. I moved this into a separate function.

Ralf

[kilir]

Hi Ralf, thanks for looking into that. Yes, I already thought that it might be "the same" somewhere...

For now I'd say that for pixelated data (almost constant segment length) it's idential (with minor exceptions - see below) and I thought that the consideration of segment length is what is taken into account in principalComponets in the weighting by segments length, I just didn't know how to implement it.

If the previous behaviour in principalComponents was actually a bug, we might wonder if we need two functions which produce essentially the same.

There's just a little issue right now with principalComponents: Some shapes won't be fitted (NaN axes length), I guess those with a high symmetry with resepct to the coordiante axes. It works again as soon as grains are smoothed; also then minor differences between both functions become visisble.

thick black stippled line: prinicpalComponents thin red line: fitEllipse

pixel-grains (doesn't work for some grains with principalComponents)

%smooth grains

Cheers, Rüdiger

[ralfHielscher]

Hi Rüdiger,

I think I have fixed the issue with the small grains and I think the old code is a bit better. I will keep the function name fitEllipse (I like it better) but will redirect it to principleComponent.

Ralf.

[kilir]

Ok, now (as of 3d10f6256a73a7f400f9ff4d985a4b741f9c7f8e) it seems to be fine for both, issues are gone for pC

again, black stippled line : principalComponents red line: fitEllipse stepped

smoothed:

[kilir]

Yes, for sure, espceially since it is able to take account for unequal segment sizes and since they are doing the same thing, I also think there's no need for any redundancy. Cheers, Rüdiger

[kilir]

I guess that's settled and works all fine. One moment the grain2d/calliper might have to be optimized, but so far works for me. Closing