On this page we consider the problem of determining a smooth orientation dependent function \(f(\mathtt{ori})\) given a list of orientations \(\mathtt{ori}_n\) and a list of corresponding values \(v_n\). These values may be the volume of crystals with a specific orientation, as in the case of an ODF, or any other orientation dependent physical property.
Orientation dependent data may be stored in ASCII files with lines of Euler angles, representing the orientations, and values. Such data files can be imported by the command orientation.load
, where we have to specify the position of the columns of the Euler angles as well as of the additional properties.
fname = fullfile(mtexDataPath, 'orientation', 'dubna.csv');
[ori, S] = orientation.load(fname,'columnNames',{'phi1','Phi','phi2','values'});
As a result the command returns a list of orientations ori
and a struct S
. The struct contains one field for each additional column in our data file. In our toy example it is the field S.values
. Lets generate a discrete plot of the given orientations ori
together with the values S.values
.
plotSection(ori, S.values,'all','sigma');

The process of finding a function which coincides with the given function values in the nodes reasonably well is called approximation. MTEX support different approximation schemes: approximation by harmonic expansion, approximation by radial functions and approximation by a Bingham distribution.
Approximation by Harmonic Expansion
An approximation by harmonic expansion is computed by the command SO3FunHarmonic.approximation
% SO3F = SO3Fun.interpolate(ori,S.values,'harmonic')
SO3F = SO3FunHarmonic.approximation(ori,S.values)
plot(SO3F,'sigma')
SO3F = SO3FunHarmonic (1 → xyz)
bandwidth: 30
weight: 0.99

Note that SO3FunHarmonic.approximation
does not aim at replicating the values exactly. In fact the relative error between given data and the function approximation is
norm(SO3F.eval(ori) - S.values) / norm(S.values)
ans =
0.7743
The reason for this difference is that MTEX by default applies regularization. The default regularization parameter is \(\lambda = 0.0001\). We can switch off regularization by setting this value to \(0\).
SO3F = SO3FunHarmonic.approximation(ori,S.values,'regularization',0)
% the relative error
norm(SO3F.eval(ori) - S.values) / norm(S.values)
plot(SO3F,'sigma')
SO3F = SO3FunHarmonic (1 → xyz)
bandwidth: 30
weight: 0.67
ans =
8.8649e-07

We observe that the relative error is much smaller, however the oscillatory behavior of the approximated function indicates overfitting. A more detailed discussion about choosing a good regularization parameter can be found in the section harmonic approximation theory.
An alternative way of regularization is to reduce the harmonic bandwidth
SO3F = SO3FunHarmonic.approximation(ori,S.values,'bandwidth',16)
% the relative error
norm(SO3F.eval(ori) - S.values) / norm(S.values)
plot(SO3F,'sigma')
SO3F = SO3FunHarmonic (1 → xyz)
bandwidth: 16
weight: 1
ans =
0.7752

One big disadvantage of harmonic approximation is that the resulting function is not guarantied to be non negative, even if all given function values have been non negative.
min(SO3F)
ans =
-0.5800
This possibility to guaranty non negativity is the central advantage of kernel based approximation.
Approximation by Radial Functions
The command for approximating orientation dependent data by a superposition of radial functions is SO3FunRBF.approximation
.
% SO3F = SO3Fun.interpolate(ori,val,'odf');
SO3F = SO3FunRBF.approximation(ori,S.values,'odf');
% the relative error
norm(SO3F.eval(ori) - S.values) / norm(S.values)
plot(SO3F,'sigma')
Warning: Maximum number of iterations reached, result may not have converged to
the optimum yet.
ans =
0.0125

The option 'ODF'
ensures that the resulting function is nonnegative and is normalized to \(1\)
min(SO3F)
mean(SO3F)
ans =
0
ans =
1.0000
The key parameter when approximation by radial functions is the halfwidth of the kernel function. This can be set by the option 'halfwidth'
. A large halfwidth results in a very smooth approximating function whereas a very small halfwidth may result in overfitting
SO3F = SO3FunRBF.approximation(ori,S.values,'halfwidth',2.5*degree,'odf');
plot(SO3F,'sigma')
Warning: Maximum number of iterations reached, result may not have converged to
the optimum yet.

A more detailed discussion about the correct choice of the halfwidth and other options can be found in the section theory of RBF approximation.
If we omit the option 'odf'
the resulting function may have negative values similar to the harmonic setting
SO3F = SO3FunRBF.approximation(ori,S.values);
% the relative error
norm(SO3F.eval(ori) - S.values) / norm(S.values)
plot(SO3F,'sigma')
ans =
0.0313

Approximating using the Bingham distribution
Approximation with the Bingham distribution currently works only with no symmetry. TODO!
cs = crystalSymmetry("1");
ori = equispacedSO3Grid(cs)
odf = fibreODF(fibre.rand(ori.CS))
figure(1)
plot(odf)
ori = SO3Grid (1 → xyz)
grid: 119088 orientations, resolution: 5°
odf = SO3FunCBF (1 → xyz)
kernel: de la Vallee Poussin, halfwidth 10°
fibre : (71-1) || -9,4,1
weight: 1

SO3F = calcBinghamODF(ori,'weights',odf.eval(ori))
figure(2)
plot(SO3F)
SO3F = SO3FunBingham (1 → xyz)
kappa: 0 0.066 92 93
weight: 1
