diff --git a/EBSDAnalysis/@EBSD/plotUnitCells.m b/EBSDAnalysis/@EBSD/plotUnitCells.m index 3fccc6681..b4b360000 100644 --- a/EBSDAnalysis/@EBSD/plotUnitCells.m +++ b/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'); diff --git a/EBSDAnalysis/@grain2d/plot.m b/EBSDAnalysis/@grain2d/plot.m index 105331240..b7a713c48 100644 --- a/EBSDAnalysis/@grain2d/plot.m +++ b/EBSDAnalysis/@grain2d/plot.m @@ -157,7 +157,9 @@ 'parent', mP.ax,'DisplayName',grains.mineralList{k},varargin{:}); %#ok % 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 diff --git a/EBSDAnalysis/@parentGrainReconstructor/private/getC2CPairs.m b/EBSDAnalysis/@parentGrainReconstructor/private/getC2CPairs.m index 2422dfcce..1e4e3e344 100644 --- a/EBSDAnalysis/@parentGrainReconstructor/private/getC2CPairs.m +++ b/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); diff --git a/EBSDAnalysis/@parentGrainReconstructor/private/getP2PPairs.m b/EBSDAnalysis/@parentGrainReconstructor/private/getP2PPairs.m index 2cf3be703..3c58996f5 100644 --- a/EBSDAnalysis/@parentGrainReconstructor/private/getP2PPairs.m +++ b/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 \ No newline at end of file diff --git a/ODFAnalysis/@ODF/calcComponents.m b/ODFAnalysis/@ODF/calcComponents.m index 062c550c0..381f1d27a 100644 --- a/ODFAnalysis/@ODF/calcComponents.m +++ b/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(omega0; -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)) 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); diff --git a/S2Fun/@S2AxisFieldHarmonic/approximation.m b/S2Fun/@S2AxisFieldHarmonic/approximation.m index 3b7d52f27..de1940ef7 100644 --- a/S2Fun/@S2AxisFieldHarmonic/approximation.m +++ b/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); diff --git a/S2Fun/@S2AxisFieldHarmonic/quadrature.m b/S2Fun/@S2AxisFieldHarmonic/quadrature.m index b69fd9761..bbae14069 100644 --- a/S2Fun/@S2AxisFieldHarmonic/quadrature.m +++ b/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') diff --git a/S2Fun/@S2FunHarmonic/approximation.m b/S2Fun/@S2FunHarmonic/approximation.m index d9cba55d0..c23f7403a 100644 --- a/S2Fun/@S2FunHarmonic/approximation.m +++ b/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(:)); diff --git a/S2Fun/@S2FunHarmonic/entropy.m b/S2Fun/@S2FunHarmonic/entropy.m index 9985c7aa8..0d38832e1 100644 --- a/S2Fun/@S2FunHarmonic/entropy.m +++ b/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); diff --git a/S2Fun/@S2FunHarmonic/eval.m b/S2Fun/@S2FunHarmonic/eval.m index 61880e018..7092cca9a 100644 --- a/S2Fun/@S2FunHarmonic/eval.m +++ b/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,7 +10,6 @@ % f - double % -v = v(:); if sF.bandwidth == 0 f = ones(size(v)).*sF.fhat/sqrt(4*pi); return; @@ -18,9 +17,15 @@ % 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 diff --git a/S2Fun/@S2FunHarmonic/quadrature.m b/S2Fun/@S2FunHarmonic/quadrature.m index 1f35cd34f..26f12c88f 100644 --- a/S2Fun/@S2FunHarmonic/quadrature.m +++ b/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); diff --git a/S2Fun/@S2FunHarmonic/times.m b/S2Fun/@S2FunHarmonic/times.m index 345f80ae8..6266eb243 100644 --- a/S2Fun/@S2FunHarmonic/times.m +++ b/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 diff --git a/S2Fun/@S2VectorFieldHarmonic/approximation.m b/S2Fun/@S2VectorFieldHarmonic/approximation.m index 546be2555..476ded926 100644 --- a/S2Fun/@S2VectorFieldHarmonic/approximation.m +++ b/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; diff --git a/S2Fun/@S2VectorFieldHarmonic/quadrature.m b/S2Fun/@S2VectorFieldHarmonic/quadrature.m index 501a2bee3..85bc5c9fd 100644 --- a/S2Fun/@S2VectorFieldHarmonic/quadrature.m +++ b/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') diff --git a/S2Fun/S2BumpKernel.m b/S2Fun/S2BumpKernel.m index d6788c8f9..533f0ccc6 100644 --- a/S2Fun/S2BumpKernel.m +++ b/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); diff --git a/TensorAnalysis/@stiffnessTensor/energyVector.m b/TensorAnalysis/@stiffnessTensor/energyVector.m index 1aefe74ee..2fac68a05 100644 --- a/TensorAnalysis/@stiffnessTensor/energyVector.m +++ b/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 diff --git a/extern/nfft_openMP/nfsftmex.mexa64 b/extern/nfft_openMP/nfsftmex.mexa64 index ab96e0c4b..3bf500538 100755 Binary files a/extern/nfft_openMP/nfsftmex.mexa64 and b/extern/nfft_openMP/nfsftmex.mexa64 differ diff --git a/geometry/@orientation/byMiller.m b/geometry/@orientation/byMiller.m index eee04105b..6b74fac24 100644 --- a/geometry/@orientation/byMiller.m +++ b/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 diff --git a/geometry/@orientation/unique.m b/geometry/@orientation/unique.m index dd27810df..f4eeaea87 100644 --- a/geometry/@orientation/unique.m +++ b/geometry/@orientation/unique.m @@ -27,13 +27,18 @@ [~,ndx,pos] = unique@rotation(ori,varargin{:}); -else +elseif length(ori)<2000 rot = rotation(symmetrise(ori,varargin{:})); [~,~,pos] = unique(rot,varargin{:}); [~,ndx,pos] = unique(min(reshape(pos,size(rot)),[],1)); - + +else + + ori = ori.project2FundamentalRegion; + [~,ndx,pos] = unique@rotation(ori,varargin{:}); + end ori = ori.subSet(ndx); \ No newline at end of file diff --git a/geometry/@rotation/unique.m b/geometry/@rotation/unique.m index 4242014b8..38c4858c4 100644 --- a/geometry/@rotation/unique.m +++ b/geometry/@rotation/unique.m @@ -22,21 +22,35 @@ % unique % -a = r.a(:); -b = r.b(:); -c = r.c(:); -d = r.d(:); -i = r.i(:); -abcd = [a.^2,b.^2,c.^2,d.^2,a.*b,a.*c,a.*d,b.*c,b.*d,c.*d,i]; +a = r.a(:); b = r.b(:); c = r.c(:); d = r.d(:); -tol = get_option(varargin,'tolerance',1e-3); +if length(r) < 1000 && ~check_option(varargin,'tolerance') + + abcd = [a.^2,b.^2,c.^2,d.^2,a.*b,a.*c,a.*d,b.*c,b.*d,c.*d, r.i(:)]; + + tol = get_option(varargin,'tolerance',1e-3); + + % in case it should not be sorted + if check_option(varargin,'stable') + [~,ir,iu] = unique(round(abcd./tol),'rows','stable'); + else + [~,ir,iu] = uniquetol(abcd,tol,'ByRows',true,'DataScale',1); + end + +else % faster but less accurate + + + tol = get_option(varargin,'tolerance',5*degree); + + % ensure upper hemisphere + ind = d < 0; + a(ind) = -a(ind); b(ind) = -b(ind); + c(ind) = -c(ind); d(ind) = -d(ind); + + [~,ir,iu] = uniquetol([a,b,c,d],tol / 100 / degree, ... + 'ByRows',true,'DataScale',1); -% in case it should not be sorted -if check_option(varargin,'stable') - [~,ir,iu] = unique(round(abcd./tol),'rows','stable'); -else - [~,ir,iu] = uniquetol(abcd,tol,'ByRows',true,'DataScale',1); end % remove duplicated points @@ -44,4 +58,5 @@ r.b = b(ir); r.c = c(ir); r.d = d(ir); -r.i = i(ir); +r.i = r.i(ir); + diff --git a/geometry/geometry_tools/private/getPolarRange.m b/geometry/geometry_tools/private/getPolarRange.m index 6d91cc6f2..db27ce921 100644 --- a/geometry/geometry_tools/private/getPolarRange.m +++ b/geometry/geometry_tools/private/getPolarRange.m @@ -9,7 +9,11 @@ % % set up fundamental region -bounds.FR = {0,pi,0,2*pi}; +if check_option(varargin,'FSFT') + bounds.FR = {0,pi,-pi,pi}; +else + bounds.FR = {0,pi,0,2*pi}; +end bounds.VR{3} = get_option(varargin,'minRho',bounds.FR{3}); bounds.VR{4} = get_option(varargin,'maxRho',bounds.FR{4}); diff --git a/geometry/geometry_tools/quadratureS2Grid.m b/geometry/geometry_tools/quadratureS2Grid.m index 23634e627..8126a548e 100644 --- a/geometry/geometry_tools/quadratureS2Grid.m +++ b/geometry/geometry_tools/quadratureS2Grid.m @@ -1,55 +1,88 @@ -function [S2G, W, M2] = quadratureS2Grid(bandwidth, varargin) +function [S2G, W, bandwidth] = quadratureS2Grid(bandwidth, varargin) % % Syntax % [S2G, W, M2] = quadratureS2Grid(M) quadrature grid of type chebyshev % [S2G, W, M2] = quadratureS2Grid(M, 'gauss') quadrature grid of type gauss % - persistent S2G_p; persistent W_p; -persistent M2_p; +persistent bw_p; -if check_option(varargin, 'gauss') - load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),'gridIndex'); +if check_option(varargin, 'gauss') || check_option(varargin, 'chebyshev') + if check_option(varargin, 'gauss') + load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),'gridIndex'); + elseif check_option(varargin, 'chebyshev') + load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),'gridIndex'); + end + index = find(gridIndex.bandwidth >= bandwidth, 1); + if isempty(index) + index = size(gridIndex,1); + warning('M is too large, instead we are giving you the largest quadrature grid we got.'); + end + bandwidth = gridIndex.bandwidth(index); + N = gridIndex.N(index); +elseif check_option(varargin,'clenshawCurtis') + bandwidth = 2*ceil(bandwidth/2); % ensure bandwidth to be even + N = (bandwidth+1)*(bandwidth+2); else - load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),'gridIndex'); -end -index = find(gridIndex.bandwidth >= bandwidth, 1); -if isempty(index) - index = size(gridIndex,1); - warning('M is too large, instead we are giving you the largest quadrature grid we got.'); + bandwidth = 2*ceil(bandwidth/2); % ensure bandwidth to be even + N = (bandwidth+1)*(bandwidth/2+1); end -M2 = gridIndex.bandwidth(index); - -if ~isempty(M2_p) && M2_p == M2 +if ~isempty(bw_p) && bw_p == bandwidth && length(S2G_p) == N S2G = S2G_p; W = W_p; - M2 = M2_p; + bandwidth = bw_p; else - name = cell2mat(gridIndex.name(index)); - if check_option(varargin, 'gauss') - data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),name); - else - data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),name); - end + if check_option(varargin, 'gauss') || check_option(varargin, 'chebyshev') + name = cell2mat(gridIndex.name(index)); + if check_option(varargin, 'gauss') + data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),name); + else + data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),name); + end + + data = data.(name); + S2G = vector3d.byPolar(data(:, 1), data(:, 2)); - data = data.(name); - S2G = vector3d.byPolar(data(:, 1), data(:, 2)); - - if check_option(varargin, 'gauss') - W = data(:,3); + if check_option(varargin, 'gauss') + W = data(:,3); + else + W = 4*pi/size(data, 1) .* ones(size(S2G)); + end + elseif check_option(varargin,'clenshawCurtis') + S2G = regularS2Grid('clenshawCurtis','bandwidth',bandwidth); + + w = fclencurt2(size(S2G,1)); + W = 2*pi/size(S2G,2)*repmat(w,1,size(S2G,2)); + + %W = W(:); else - W = 4*pi/size(data, 1) .* ones(size(S2G)); - end - + S2G = regularS2Grid('FSFT','bandwidth',bandwidth); + + w = fclencurt2(size(S2G,1)); + W = 2*pi/size(S2G,2)*repmat(w,1,size(S2G,2)); + end % store the data S2G_p = S2G; W_p = W; - M2_p = M2; + bw_p = bandwidth; +end + end +function w = fclencurt2(n) + +c = zeros(n,2); +c(1:2:n,1) = (2./[1 1-(2:2:n-1).^2 ])'; +c(2,2)=1; +f = real(ifft([c(1:n,:);c(n-1:-1:2,:)])); +w = 2*([f(1,1); 2*f(2:n-1,1); f(n,1)])/2; +%x = 0.5*((b+a)+n*bma*f(1:n1,2)); + end + + diff --git a/geometry/geometry_tools/regularS2Grid.m b/geometry/geometry_tools/regularS2Grid.m index 04efeb703..33b235ab1 100644 --- a/geometry/geometry_tools/regularS2Grid.m +++ b/geometry/geometry_tools/regularS2Grid.m @@ -23,6 +23,41 @@ % See also % equispacedS2Grid plotS2Grid +if check_option(varargin,'FSFT') + % a regular grid that fits the FSFT + % required spherical harmonic bandwidth 2*n + % 2n+1 Clenshaw Curtis Quadrature rule + % -> (2n+1)x(2n+2) points + + N = ceil(get_option(varargin,'bandwidth',256)/2); + + theta = linspace(0,pi,N+2); + rho = (-N-1:N)/(2*N+2)*2*pi; + [rho,theta] = meshgrid(rho,theta); + S2G = vector3d.byPolar(theta,rho); + + S2G = S2G.addOption('using_fsft',N); + return + +elseif check_option(varargin,'ClenshawCurtis') + + % a regular grid for ClenshawCurtis quadrature + % required spherical harmonic bandwidth 2*n + % 2n+1 Clenshaw Curtis Quadrature rule + % -> (2n+1)x(2n+2) points + + N = ceil(get_option(varargin,'bandwidth',256)/2); + + theta = linspace(0,pi,2*N+1); + rho = (-N-1:N)/(2*N+2)*2*pi; + [rho,theta] = meshgrid(rho,theta); + S2G = vector3d.byPolar(theta,rho); + + S2G = S2G.addOption('using_fsft',N); + return + +end + % extract options bounds = getPolarRange(varargin{:}); diff --git a/mtex_settings.m b/mtex_settings.m index 785a7f382..4970c7816 100644 --- a/mtex_settings.m +++ b/mtex_settings.m @@ -204,6 +204,7 @@ setMTEXpref('FFTAccuracy',1E-2); setMTEXpref('maxBandwidth',512); +setMTEXpref('NFSFTBandwidth',256); %% degree character % MTEX sometimes experences problems when printing the degree character diff --git a/tools/math_tools/sphericalY.m b/tools/math_tools/sphericalY.m index 3bb560143..ff1b29c98 100644 --- a/tools/math_tools/sphericalY.m +++ b/tools/math_tools/sphericalY.m @@ -1,4 +1,4 @@ -function Y = sphericalY(l, v, varargin) +function Y = sphericalY(n, v, varargin) % spherical harmonics of degree l % % Description @@ -21,40 +21,45 @@ [theta,rho] = polar(v); % try NFSFT version -if check_option(varargin,'nfsft') && l > 0 +if check_option(varargin,'nfsft') && n > 0 % precomputation for nfsft - nfsft_precompute(l,1000); - plan = nfsft_init_advanced(l,1,NFSFT_NORMALIZED); + nfsftmex('precompute',n,1000,1,0); + plan = nfsftmex('init_advanced',n,length(v),1); - nfsft_set_x(plan,[rho;theta]); + nfsftmex('set_x',plan,[rho(:).';theta(:).']); % node-dependent precomputation - nfsft_precompute_x(plan); + nfsftmex('precompute_x',plan); - % Set Fourier coefficients. - nfsft_set_f(plan,1); + Y = zeros(length(v),2*n+1); + for k = -n:n - % transform - nfsft_adjoint(plan); - - % store results - Y = nfsft_get_f_hat_linear(plan); - Y = Y((end-2*l):end)'; + % set Fourier coefficients. + fhat = zeros((n+1)^2,1); + fhat(n^2+n+1+k) = 1; + nfsftmex('set_f_hat_linear',plan,fhat); + + % transform + nfsftmex('trafo',plan); - nfsft_finalize(plan); + % store results + Y(:,k+n+1) = nfsftmex('get_f',plan); + %Y = Y((end-2*n):end)'; + end + nfsftmex('finalize',plan); return end % calculate assoziated legendre functions -L = sqrt(2/(2*l+1))*reshape(legendre(l,cos(theta(:)),'norm').',numel(theta),l+1); % nodes x order +L = sqrt(2/(2*n+1))*reshape(legendre(n,cos(theta(:)),'norm').',numel(theta),n+1); % nodes x order % expand to negative order -if l>0 +if n>0 L = [fliplr(L(:,2:end)),L]; end % calcualte spherical harmonics -Y = sqrt((2*l+1)/4/pi) .* L .* exp(1i*rho(:) * (-l:l)); +Y = sqrt((2*n+1)/4/pi) .* L .* exp(1i*rho(:) * (-n:n)); diff --git a/tools/statistic_tools/maxVote.m b/tools/statistic_tools/maxVote.m new file mode 100644 index 000000000..4fd4d814c --- /dev/null +++ b/tools/statistic_tools/maxVote.m @@ -0,0 +1,33 @@ +function maxInd = maxVote(ind,data,varargin) +% returns +% +% Syntax +% +% maxInd = maxVote(ind,data) +% +% Input +% idList - list of ids +% voteList - list of votes +% +% Output +% vote - most frequent vote in voteList(idList == id) for each id +% numVotes - number occurences of vote in voteList(idList == id) for each id +% +% Options +% strict - only assign a vote if all votes coincide, otherwise set in to nan +% +% Example +% +% ind = [1,2,3,1,3,2,1]; +% data = [2,3,1,1,2,1,3]; +% maxInd = [7, 2, 5] +% maxInd = maxVote(ind,data) +% +% See also +% + +[~,a] = sort(data,varargin{:}); + +newInd = (1:length(data)).'; + +maxInd = accumarray(ind(a),newInd(a),[],@(x) x(end));