Feature/regular s2quadrature

EBSDAnalysis/@EBSD/plotUnitCells.m

@@ -13,16 +13,14 @@
 
   xy = xy(ind,:);
 
-  d = d(ind,:);
+  if numel(d) == numel(ind) || numel(ind) == size(d,1)
+    d = d(ind,:);
+  end
 
 end
 
-
 ax = get_option(varargin,'parent',gca);
 
-
-
-
 if ~isempty(unitCell)
 
   type = get_flag(varargin,{'unitcell','points','measurements'},'unitcell');

EBSDAnalysis/@grain2d/plot.m

@@ -157,7 +157,9 @@
       'parent', mP.ax,'DisplayName',grains.mineralList{k},varargin{:}); %#ok<AGROW>
 
     % reactivate legend information
-    set(get(get(h{k}(end),'Annotation'),'LegendInformation'),'IconDisplayStyle','on');
+    if ~isempty(h{k})
+      set(get(get(h{k}(end),'Annotation'),'LegendInformation'),'IconDisplayStyle','on');
+    end
 
   end
 
@@ -330,6 +332,7 @@ function spatialSelection(src,eventdata)
 obj.FaceColor = 'flat';
 obj.EdgeColor = 'None';
 obj.hitTest = 'off';
+h = [];
 
 for p=numel(Parts):-1:1
   zOrder = Parts{p}(end:-1:1); % reverse

EBSDAnalysis/@parentGrainReconstructor/private/getC2CPairs.m

@@ -1,5 +1,20 @@
-function [pairs,ori] = getC2CPairs(job,varargin)
+function [pairs,ori,weights] = getC2CPairs(job,varargin)
+% extract child to child grain pairs with their orientations 
+%
 % 
+% Options
+%  index    - return as index to job.grains instread of grainId
+%  minDelta - ignore grain pairs with similar orientations as they can not 
+%             vote for a unique parent orientation
+%  curvatureFactor - compute boundary weight from the curvature (default 0)
+%  job.useBoundaryOrientations - instead of the grain mean orientations
+%  average the orientations along the grain boundaries
+%
+% Output
+%  pairs - N x 2 list of grainIds
+%  ori   - N x 2 list of grain orientations
+%  weights - N x 1 list of boundary weights (1 if curvature is not considered)
+%
 
 % all child to child pairs
 pairs = neighbors(job.grains);

EBSDAnalysis/@parentGrainReconstructor/private/getP2PPairs.m

@@ -26,8 +26,14 @@
 else 
 
   % simply the mean orientations of the grains
-  ori = job.grains('id',pairs).meanOrientation;
+  %ori = job.grains('id',pairs).meanOrientation;
+  ori = orientation(job.grains.prop.meanRotation(job.grains.id2ind(pairs)),...
+    job.csParent);
+  ori = reshape(ori,size(pairs));
 
 end
 
+% translate to index if required
+if check_option(varargin,'index'), pairs = job.grains.id2ind(pairs); end
+
 end

ODFAnalysis/@ODF/calcComponents.m

@@ -16,6 +16,7 @@
 %
 % Options
 %  resolution - search-grid resolution
+%  angle      - maximum component width used for volume computation
 %  exact      - do not dismiss very small modes at the end
 %
 % Example
@@ -28,41 +29,44 @@
 % See also
 % ODF/max
 
+% extract options
 maxIter = get_option(varargin,'maxIter',100);
 res = get_option(varargin,'resolution',0.05*degree);
 omega = 1.5.^(-7:1:4) * degree;
 omega(omega<res) = [];
 omega = [0,omega];
-tol = get_option(varargin,'tolerance',0.05);
+tol = get_option(varargin,'tolerance',0.5*degree);
+maxAngle = get_option(varargin,{'radius','angle'},inf);
 
 % initial seed
 if check_option(varargin,'seed')
-  modes = reshape(get_option(varargin,'seed'),[],1);
+  seed = reshape(get_option(varargin,'seed'),[],1);
   weights = get_option(varargin,'weights',...
-    odf.eval(modes));
+    odf.eval(seed));
 elseif isa(odf.components{end},'unimodalComponent')
-  modes = odf.components{end}.center;
+  seed = odf.components{end}.center;
   weights = odf.components{end}.weights;
 else
-  modes = equispacedSO3Grid(odf.CS,odf.SS,'resolution',5*degree);
-  weights = ones(length(modes),1) ./ length(modes);
+  seed = equispacedSO3Grid(odf.CS,odf.SS,'resolution',2.5*degree);
+  weights = ones(length(seed),1) ./ length(seed);
 end
 id = weights>0;
-modes = reshape(modes(id),[],1);
+seed = reshape(seed(id),[],1);
 weights = weights(id);
 weights = weights ./ sum(weights);
 
-centerId = 1:length(modes);
+centerId = 1:length(seed);
+modes = seed;
 
 % join orientations if possible
-[modes,~,id2] = unique(modes,'tolerance',tol);
-centerId = id2(centerId);
-weights = accumarray(id2,weights);
+%[modes,~,id2] = unique(modes,'tolerance',tol);
+%centerId = id2(centerId);
+%weights = accumarray(id2,weights);
 
 finished = false(size(modes));
 
 for k = 1:maxIter
-  progress(k,maxIter,' finding ODF components: ');
+  %progress(k,maxIter,' finding ODF components: ');
 
   % gradient
   g = normalize(odf.grad(modes(~finished)),1);
@@ -74,21 +78,30 @@
   line_v = odf.eval(line_ori);
 
   % take the maximum
-  [~,id] = max(line_v,[],2);
+  [v_max(~finished),id] = max(line_v,[],2);
 
   % update orientions
   modes(~finished) = line_ori(sub2ind(size(line_ori),(1:length(g)).',id));
 
-  finished(~finished) = id == 1;
-  if all(finished), break; end
+  nnz(id>1)
+  if all(id == 1), break; end
 
   % join orientations if possible
-  [modes,~,id2] = unique(modes,'tolerance',tol);
+  [~,~,id2] = unique(modes,'tolerance',tol);
+
+  modes = modes(maxVote(id2,v_max));
+
   centerId = id2(centerId);
-
-  weights = accumarray(id2,weights);
-  finished = accumarray(id2,finished,[],@all);
-  %[length(modes), k, sum(finished)]
+  finished = accumarray(id2,finished,[],@any);
+  if maxAngle == inf, weights = accumarray(id2,weights); end
+  length(modes);
+
+end
+
+% accumulate only weights that are sufficiently close to the centers
+if maxAngle < inf
+  inRadius = angle(seed,modes(centerId))<maxAngle;
+  weights = accumarray(centerId(inRadius), weights(inRadius));
 end
 
 % sort components according to volume
@@ -98,7 +111,7 @@
 centerId = iid(centerId);
 
 if ~check_option(varargin,'exact')
-  id = weights > 0.01;
+  id = weights > min([0.01, numSym(odf.CS) * maxAngle^3 ./ 8*pi^2,0.5 * max(weights)]);  
   weights = weights(id);
   modes = modes(id);
   ids = 1:length(id);

S2Fun/@S2AxisFieldHarmonic/approximation.m

@@ -12,7 +12,7 @@
 %   sAF - @S2AxisFieldHarmonic
 %
 % Options
-%   bw - degree of the spherical harmonic (default: 128)
+%   bw - degree of the spherical harmonic (default: 256)
 %
 
 [x,y,z] = double(y);

S2Fun/@S2AxisFieldHarmonic/quadrature.m

@@ -11,7 +11,7 @@
 %  f - function handle in @vector3d
 %
 % Options
-%  M - degree of the spherical harmonic (default: 128)
+%  M - degree of the spherical harmonic (default: 256)
 %
 
 if isa(f,'vector3d')

S2Fun/@S2FunHarmonic/approximation.m

@@ -57,7 +57,7 @@
 if isempty(W) 
   W = nodes.calcVoronoiArea;
 else
-  W = W(ind);
+  W = accumarray(ind,W);
 end
 W = sqrt(W(:));
 

S2Fun/@S2FunHarmonic/entropy.m

@@ -16,7 +16,7 @@
 %
 
 % get some quadrature nodes
-[nodes, W] = quadratureS2Grid(max(128,sF.bandwidth));
+[nodes, W] = quadratureS2Grid(max(getMTEXpref('NFSFTBandwidth'),sF.bandwidth));
 
 f = sF.eval(nodes);
 

S2Fun/@S2FunHarmonic/eval.m

@@ -1,4 +1,4 @@
-function f =  eval(sF,v)
+function f = eval(sF,v)
 % evaluates the spherical harmonic on a given set of points
 % Syntax
 %   f = eval(sF,v)
@@ -10,17 +10,22 @@
 %   f - double
 %
 
-v = v(:);
 if sF.bandwidth == 0
   f = ones(size(v)).*sF.fhat/sqrt(4*pi);
   return;
 end
 
 % initialize nfsft
 nfsftmex('precompute', sF.bandwidth, 1000, 1, 0);
-plan = nfsftmex('init_advanced', sF.bandwidth, length(v), 1);
-nfsftmex('set_x', plan, [v.rho'; v.theta']); % set vertices
-nfsftmex('precompute_x', plan);
+
+if v.isOption('using_fsft')
+  sF.bandwidth = v.opt.using_fsft;
+  plan = nfsftmex('init_advanced', sF.bandwidth, length(v), 1+2^17);
+else
+  plan = nfsftmex('init_advanced', sF.bandwidth, length(v), 1);
+  nfsftmex('set_x', plan, [v.rho(:).'; v.theta(:).']); % set vertices
+  nfsftmex('precompute_x', plan);
+end
 
 f = zeros([length(v) size(sF)]);
 % nfsft

S2Fun/@S2FunHarmonic/quadrature.m

@@ -14,29 +14,33 @@
 %  sF - @S2FunHarmonic
 %
 % Options
-%  bandwidth - minimal degree of the spherical harmonic (default: 128)
+%  bandwidth - minimal degree of the spherical harmonic (default: 256)
 %
 
 persistent keepPlan;
 
 % kill plan
 if check_option(varargin,'killPlan')
-  nfsftmex('finalize',keepPlan);
+  if ~isempty(keepPlan), nfsftmex('finalize',keepPlan); end
   keepPlan = [];
   return
 end
 
-bw = get_option(varargin, 'bandwidth', 128);
+bw = get_option(varargin, 'bandwidth', getMTEXpref('NFSFTBandwidth'));
 
 if isa(f,'S2Fun'), f = @(v) f.eval(v); end
 
 if isa(f,'function_handle')
+  % we need a qudrature grid for twice the bandwidth to compute the Fourier
+  % coefficients
   if check_option(varargin, 'gauss')
-    [nodes, W] = quadratureS2Grid(2*bw, 'gauss');
+    [nodes, W] = quadratureS2Grid(2*bw, 'gauss');   
+  elseif check_option(varargin, 'chebyshev')
+    [nodes, W] = quadratureS2Grid(2*bw, 'chebyshev');
   else
     [nodes, W] = quadratureS2Grid(2*bw);
   end
-  values = f(nodes(:));
+  values = f(nodes);
 else
   nodes = f(:);
   values = varargin{1};
@@ -70,9 +74,13 @@
 % initialize nfsft
 if isempty(plan)
   nfsftmex('precompute', bw, 1000, 1, 0);
-  plan = nfsftmex('init_advanced', bw, length(nodes), 1);
-  nfsftmex('set_x', plan, [nodes.rho'; nodes.theta']); % set vertices
-  nfsftmex('precompute_x', plan);
+  if nodes.isOption('using_fsft')
+    plan = nfsftmex('init_advanced', bw, length(nodes), 1+2^17);
+  else
+    plan = nfsftmex('init_advanced', bw, length(nodes), 1);
+    nfsftmex('set_x', plan, [nodes.rho'; nodes.theta']); % set vertices
+    nfsftmex('precompute_x', plan);
+  end
 end
 
 s = size(values);

S2Fun/@S2FunHarmonic/times.m

@@ -23,7 +23,7 @@
   f = @(v) sF1.eval(v) .* sF2.eval(v);
   sF = S2FunHarmonic.quadrature(f, 'bandwidth', min(getMTEXpref('maxBandwidth'),sF1.bandwidth + sF2.bandwidth));
 else
-  sF = sF2.*sF1
+  sF = sF2.*sF1;
 end
 
 end

S2Fun/@S2VectorFieldHarmonic/approximation.m

@@ -12,7 +12,7 @@
 %   sVF - @S2VectorFieldHarmonic
 %
 % Options
-%   bw - degree of the spherical harmonic (default: 128)
+%   bw - degree of the spherical harmonic (default: 256)
 %
 
 y = y.xyz;

S2Fun/@S2VectorFieldHarmonic/quadrature.m

@@ -14,7 +14,7 @@
 %   sVF - @S2VectorFieldHarmonic
 %
 % Options
-%   bw - degree of the spherical harmonic (default: 128)
+%   bw - degree of the spherical harmonic (default: 256)
 %
 
 if isa(f,'vector3d')

S2Fun/S2BumpKernel.m

@@ -30,7 +30,7 @@
       end
 
       % extract bandwidth
-      L = get_option(varargin,'bandwidth',128);
+      L = get_option(varargin,'bandwidth',getMTEXpref('NFSFTBandwidth'));
 
       % compute Legendre coefficients
       psi.A = nan(1,L+1);

TensorAnalysis/@stiffnessTensor/energyVector.m

@@ -38,7 +38,8 @@
   if isa(V,'S2FunTri')
     E = S2VectorFieldTri(V.tri,energyVector(C,V.vertices,V.values,P.values,varargin{:}));
   else
-    E = S2VectorFieldHarmonic.quadrature(@(x) energyVector(C,x,V,P,varargin{:}),'bandwidth',128,C.CS);
+    E = S2VectorFieldHarmonic.quadrature(@(x) energyVector(C,x,V,P,varargin{:}),...
+      'bandwidth',getMTEXpref('NFSFTBandwidth'),C.CS);
   end
   return
 end

geometry/@orientation/byMiller.m

@@ -8,17 +8,17 @@
 %
 % Syntax
 %   
-%   hkl = Miller(h,k,l,CS,'hkl')
-%   uvw =  Miller(u,v,w,CS,'uvw')
+%   hkl = Miller(h,k,l,cs,'hkl')
+%   uvw =  Miller(u,v,w,cs,'uvw')
 %   ori = orientation.byMiller(hkl,uvw)
 %
-%   ori = orientation.byMiller([h k l],[u v w],CS)
+%   ori = orientation.byMiller([h k l],[u v w],cs)
 %
 % Input
 %  hkl, uvw - @Miller
 %  h,k,l  - Miller indece (double)
 %  u,v,w  - Miller indece (double)
-%  CS     - @crystalSymmetry
+%  cs     - @crystalSymmetry
 %
 % Output
 %  ori - @orientation