ita_sph_sampling_equalarea.m 2.13 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
function varargout = ita_sph_sampling_equalarea(varargin)
%ITA_SPH_SAMPLING_EQUALAREA - Equal area sampling on the sphere
%  This function generates a spherical sampling based on the center points
%  of faces with equal area on the sphere. Algorithm based on Leopardi's
%  equalarea algorithm.
%  In order to ensure a feasible SHT the number of sampling points is 
%  increased until the condition number of the SH basis matrix is sufficiently small.
%  For high SH orders it is recommended to start with a higher number of sampling points
%  than (Nmax+1)^2 in order to save computation time.
%  In order to ignore the feasibility criterion, use 'condSHT', inf
%
%  Syntax:
%   sampling = ita_sph_sampling_equalarea(Nmax,options)
%
%  Example:
%   sampling = ita_sph_sampling_equalarea(Nmax,'condSHT',2)
%
%  See also:
%   ita_toolbox_gui, ita_read, ita_write, ita_generate
%
%   Reference page in Help browser 
%        <a href="matlab:doc ita_sph_sampling_equalarea">doc ita_sph_sampling_equalarea</a>

% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved. 
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder. 
% </ITA-Toolbox>


% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de
% Created:  09-Nov-2016 

sArgs = struct('pos1_Nmax','integer',...
               'condSHT',2.5,...
			   'nPoints',[]);
[Nmax,sArgs] = ita_parse_arguments(sArgs,varargin);
           
if isempty(sArgs.nPoints)
	sArgs.nPoints = (Nmax+1)^2;
else
	if Nmax > 15
		ita_verbose_info('You may want to consider setting a number of points higher than (Nmax+1)^2 as starting point.',1);
	end
end


% find a sampling with a feasible SHT transform as this may not be the case
% for every set of sampling points resulting from the equal area
% partitioning
% use while true to exec the loop at least once
while true
    coordsCart = eq_point_set(2,sArgs.nPoints).';
    sampling = itaSamplingSph(coordsCart,'cart');
    Y = mbe_sph_base(sampling,Nmax);
    condNum = cond(Y);
    if condNum < sArgs.condSHT
        break;
    end
    sArgs.nPoints = sArgs.nPoints+1;
end

sampling.nmax = Nmax;

varargout{1} = sampling;

end