mtex-toolbox / mtex

MTEX is a free Matlab toolbox for quantitative texture analysis. Homepage:
http://mtex-toolbox.github.io/
GNU General Public License v2.0
281 stars 185 forks source link

neighbors problem #865

Open neilmanck opened 3 years ago

neilmanck commented 3 years ago

pairs = grains('Quartz-hex').neighbors; qtz_qtz_mean = inv(grains(pairs(:,1)).meanOrientation) .* grains(pairs(:,2)).meanOrientation;

With these commands I get the following error message. From my understanding the first command should only return neighbors for the phase 'Quartz-hex' so I should not get the error message. A bug?

regards, Neil

Error using phaseList/checkSinglePhase (line 307)


Your variable contains the phases: Quartz-hex, Albite

However, your are executing a command that is only permitted for a single phase!

Please read the chaper "select EBSD data" for how to restrict EBSD data or grains to a single phase.

Error in phaseList/get.CS (line 154) id = checkSinglePhase(pL);

Error in grain2d/get.meanOrientation (line 266) ori = orientation(grains.prop.meanRotation,grains.CS);

Error in grain2d/subsref (line 44) [varargout{1:nargout}] = builtin('subsref',grains,s);

neilmanck commented 3 years ago

Sorry, I should have also noted that this is with Version 5.6.0

kilir commented 3 years ago

Hi Neil, I guess the index of the grains does not coincide with the grainId anymore (e.g. some grains from the list were deleted) after computing grains.

Does calling them by id solve the problem?

qtz_qtz_mean = inv(grains('id',pairs(:,1)).meanOrientation) .* grains('id',pairs(:,2)).meanOrientation;

Cheers, Rüdiger

neilmanck commented 3 years ago

Hi Rüdiger, Thanks for the suggestion. No, just changing that line did not solve the problem but you are certainly on to root of the problem. Indeed some poorly constrained grains were deleted following the approach of Cross et al. . If I comment these lines out and do not remove the "poorly constrained grains" then everything works fine with the original code: pairs = grains('Quartz-hex').neighbors; qtz_qtz_mean = inv(grains(pairs(:,1)).meanOrientation) .* grains(pairs(:,2)).meanOrientation;

So can I please use your knowledge on how I should update the grainId correctly after removing grains? Otherwise I will just skip the lines removing the "poorly constrained grains" - it is probably overkill for the misorientation analysis anyway, but more important when considering grain size.

The code I lifted from Cross et al. is below.

cheers, Neil

%% Remove poorly constrained grains

stepsize = ebsd(2).x - ebsd(1).x; % EBSD step-size (i.e. resolution) in microns

grains_q = grains('Quartz-hex');

% Calculate the fraction of each grain's area that is made up of indexed
% pixels (0-1)
% - for example, fraction = 0.1 means that only 10% of that grain is made 
% up of indexed pixels (in other words, it's not very well constrained)
fraction = (full(grains_q.grainSize).*(stepsize^2))./grains_q.area; 

% Use trade-off curve to find cutoff between well-constrained and 
% poorly-constrained grains
knee = tradeOff(fraction);
    xlabel('Number of grains (cumulative)')
    ylabel('Indexed fraction')

% Keep only the well-constrained quartz grains, and those made up of 4 or more pixels
condition = (full(grains_q.grainSize) <= 4) & (fraction < knee);
grains(grains_q.id(condition)) = [];
neilmanck commented 3 years ago

Hi Rüdiger, In the meantime I have found a work around by not deleting "poorly constrained grains" but simply setting their phase to zero. There may be a "more correct" way but this achieves the desired result.

thanks for your help, Neil

% Keep only the well-constrained quartz grains, and those made up of 4 or more pixels condition = (full(grains_q.grainSize) <= 4) & (fraction < knee); % grains(grains_q.id(condition)) = []; grains(grains_q.id(condition)).phase = 0;

kilir commented 3 years ago

Hi Neil, does this snipplet actually produce satisfying results in taking care of poorly indexed grains?

plot(ebsd)
hold on
plot(grains_q(grains_q.id(condition)).boundary)
hold off

since this line

grains(grains_q.id(condition)) = [];

addresses grains by their index (position in the list) using the id's of grains_q .

I'd think, this one would be more robust.

grains('id',grains_q.id(condition)) = [];

One could reset the grain ids to the indices. However this might end up being more complicated since it would also require to take care of other places e.g. in the ebsd.grainId and potentially using grains('id',...) or grains{} might be easier.

I guess you want to compute the misorientations between neighboring grains and it's a bit puzzling why they won't like to be addressed by their ids.

Here's a little demo maybe this helps in some way:

mtexdata forsterite
grains = ebsd.calcGrains

%%
% delete something
grains(grains.grainSize < 200) = [];

% also delete all not Indexed grains
grains = grains('indexed')

%%

plot(grains)

%%
% get IDs of forsterite neighbors
pairs = grains('f').neighbors;

%% example plot a grain and it's neighbors
% some grains with id 6471
idtest = 6471

plot(grains,'FaceAlpha',0.3)
hold on
% plot the neighbors
nIds = unique(pairs(any(ismember(pairs,idtest),2),:));
% this is identical to
% nIds = unique(grains.neighbors(grains{idtest}))

plot(grains('id',nIds),'faceColor','DarkGoldenrod','linecolor','k')
hold on% plot the central grains we asked about it's neighbors
plot(grains('id',idtest),'faceColor','Fuchsia','linecolor','k')
text(grains('id',idtest),  {num2str(idtest)})
hold off

%% get paris of neighboring orientations
o = grains('id',pairs).meanOrientation;
% misorientations between those 
mori = inv(o(:,1)).*o(:,2);

Cheers, Rüdiger