Mean grain orientation calculation

[gmcastroguiza]

Hello everybody,

I was wondering, can anybody explain to me how the mean orientation of a grain is calculated? I am wondering since this is not a scalar value, and I am not sure how the orientations in each point inside the grain are averaged.

Kind regards,

[zmichels]

Well... this is the function as I find it in MTEX 5.11.2 for orientations... and it looks like it references a function written for quaternion versions of the orientations...

function [m, o, lambda, eigv]  = mean(o,varargin)
% mean of a list of orientations, principle axes and moments of inertia
%
% Syntax
%   [m, q, lambda, V] = mean(ori)
%   [m, q, lambda, V] = mean(ori,'robust')
%   [m, q, lambda, V] = mean(ori,'weights',weights)
%
% Input
%  ori      - list of @orientation
%
% Options
%  weights  - list of weights
%
% Output
%  m      - mean @orientation
%  o      - crystallographic equivalent @orientation projected to fundamental region
%  lambda - principle moments of inertia
%  V      - principle axes of inertia (@orientation)
%
% See also
% BinghamODF

if isempty(o)
  m = o;
  m.a = NaN; m.b = NaN; m.c = NaN; m.d = NaN; m.i = false;
  if nargout > 2, lambda = zeros(1,4); end
  if nargout > 3, eigv = eye(4); end
  return  
elseif isscalar(o) 
  m = o;
  if nargout > 1
    eigv = eye(4);
    lambda = [1,0,0,0];
  end
  return;
end

if check_option(varargin,'noSymmetry')
  
  [m, lambda, eigv] = mean@quaternion(o,varargin{:});
    
else
  
  s = size(o);

  % first approximation
  m = get_option(varargin,'q0');
  if isempty(m), m = o.subSet(find(~isnan(o.a),1)); end
  
  % project around q_mean
  o = project2FundamentalRegion(o,m);
  
  % compute mean without symmetry
  [m, lambda, eigv] = mean@quaternion(o,varargin{:});

  d = abs(quat_dot(o,m));
  if min(d(:)) < cos(10*degree)
    o = project2FundamentalRegion(o,m);
    [m, lambda, eigv] = mean@quaternion(o,varargin{:});
  end

  if nargout > 1, o = reshape(project2FundamentalRegion(o,m),s); end
 
end
m.i = false;

the quaternion "mean" function it references is this one, I believe:

function [qm, lambda, V] = mean(q,varargin)
% mean of a list of quaternions, principle axes and moments of inertia
%
% Syntax
%
%   [m, lambda, V] = mean(q)
%   [m, lambda, V] = mean(q,'robust')
%   [m, lambda, V] = mean(q,'weights',weights)
%
% Input
%  q        - list of @quaternion
%  weights  - list of weights
%
% Output
%  m      - mean orientation
%  lambda - principle moments of inertia
%  V      - principle axes of inertia (@quaternion)
%
% See also
% orientation/mean


qm = q;

if isempty(q) || all(isnan(q.a(:)))
  
  [qm.a,qm.b,qm.c,qm.d] = deal(nan,nan,nan,nan);
  lambda = [0 0 0 1];
  if nargout == 3, V = nan(4); end
  return
elseif isscalar(q)
  lambda = [0 0 0 1];
  if nargout == 3
    T = qq(q,varargin{:});
    [V, lambda] = eig(T);
    lambda = diag(lambda).';
  end
  return
end

% shall we apply the robust algorithm?
isRobust = check_option(varargin,'robust');
if isRobust, varargin = delete_option(varargin,'robust'); end

T = qq(q,varargin{:});
[V, lambda] = eig(T);
lambda = diag(lambda).';
[~,pos] = max(lambda);

VV = V(:,pos);
qm.a = VV(1); qm.b = VV(2); qm.c = VV(3); qm.d = VV(4);
if isa(qm,'rotation'), qm.i = false; end

if isRobust && length(q)>4
  omega = angle(qm,q,'noSymmetry');
  id = omega <= quantile(omega,0.8)*2.5;
  if all(id), return; end
  if nargout == 3
    [qm,lambda, V] = mean(q.subSet(id),varargin{:});
  else
    [qm,lambda] = mean(q.subSet(id),varargin{:});
  end
end