Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Open sidebar
OpenWalnut
OpenWalnut Core
Commits
874e625f
Commit
874e625f
authored
Apr 14, 2011
by
Sebastian Eichelbaum
Browse files
[MERGE]
parents
8281c831
e1645ebf
Changes
74
Hide whitespace changes
Inline
Side-by-side
Showing
20 changed files
with
666 additions
and
291 deletions
+666
-291
src/common/WPropertyTypes.h
src/common/WPropertyTypes.h
+82
-5
src/common/WSharedObject.h
src/common/WSharedObject.h
+2
-2
src/common/math/WLinearAlgebraFunctions.cpp
src/common/math/WLinearAlgebraFunctions.cpp
+22
-23
src/common/math/WLinearAlgebraFunctions.h
src/common/math/WLinearAlgebraFunctions.h
+3
-3
src/common/math/WMatrix.h
src/common/math/WMatrix.h
+72
-1
src/common/math/WMatrix4x4.h
src/common/math/WMatrix4x4.h
+7
-5
src/common/math/WSymmetricSphericalHarmonic.cpp
src/common/math/WSymmetricSphericalHarmonic.cpp
+78
-65
src/common/math/WSymmetricSphericalHarmonic.h
src/common/math/WSymmetricSphericalHarmonic.h
+19
-18
src/common/math/WValue.h
src/common/math/WValue.h
+27
-0
src/common/math/WVector3D.cpp
src/common/math/WVector3D.cpp
+124
-0
src/common/math/WVector3D.h
src/common/math/WVector3D.h
+116
-112
src/common/math/test/WLinearAlgebraFunctions_test.h
src/common/math/test/WLinearAlgebraFunctions_test.h
+47
-8
src/common/math/test/WSymmetricSphericalHarmonic_test.h
src/common/math/test/WSymmetricSphericalHarmonic_test.h
+24
-12
src/dataHandler/WDataSetSphericalHarmonics.cpp
src/dataHandler/WDataSetSphericalHarmonics.cpp
+14
-12
src/dataHandler/WDataSetSphericalHarmonics.h
src/dataHandler/WDataSetSphericalHarmonics.h
+1
-0
src/dataHandler/WDataTexture3D_2.cpp
src/dataHandler/WDataTexture3D_2.cpp
+11
-5
src/dataHandler/WDataTexture3D_2.h
src/dataHandler/WDataTexture3D_2.h
+13
-0
src/dataHandler/WGridRegular3D.cpp
src/dataHandler/WGridRegular3D.cpp
+1
-1
src/dataHandler/WGridRegular3D.h
src/dataHandler/WGridRegular3D.h
+1
-15
src/dataHandler/WGridTransformOrtho.cpp
src/dataHandler/WGridTransformOrtho.cpp
+2
-4
No files found.
src/common/WPropertyTypes.h
View file @
874e625f
...
...
@@ -27,6 +27,11 @@
#include <stdint.h>
#include <string>
#include <vector>
#include <boost/lexical_cast.hpp>
#include <string>
#include <list>
#include <utility>
...
...
@@ -35,7 +40,7 @@
#include <boost/lexical_cast.hpp>
#include "math/WPosition.h"
#include "math/WMatrix
4x4
.h"
#include "math/WMatrix.h"
#include "WItemSelector.h"
#include "WColor.h"
...
...
@@ -98,9 +103,9 @@ namespace WPVBaseTypes
typedef
std
::
string
PV_STRING
;
//!< base type used for every WPVString
typedef
boost
::
filesystem
::
path
PV_PATH
;
//!< base type used for every WPVFilename
typedef
WItemSelector
PV_SELECTION
;
//!< base type used for every WPVSelection
typedef
WPosition
PV_POSITION
;
//!< base type used for every WPVPosition
typedef
WPosition
PV_POSITION
;
//!< base type used for every WPVPosition
typedef
WColor
PV_COLOR
;
//!< base type used for every WPVColor
typedef
WMatrix4x4
PV_MATRIX4X4
;
//!< base type used for every WPVMatrix4X4
typedef
WMatrix4x4
_2
PV_MATRIX4X4
;
//!< base type used for every WPVMatrix4X4
/**
* Enum denoting the possible trigger states. It is used for trigger properties.
...
...
@@ -543,7 +548,22 @@ namespace PROPERTY_TYPE_HELPER
*/
WPVBaseTypes
::
PV_MATRIX4X4
create
(
const
WPVBaseTypes
::
PV_MATRIX4X4
&
/*old*/
,
const
std
::
string
str
)
{
return
fromString
(
str
);
WMatrix4x4_2
c
;
std
::
vector
<
std
::
string
>
tokens
;
tokens
=
string_utils
::
tokenize
(
str
,
";"
);
WAssert
(
tokens
.
size
()
>=
16
,
"There weren't 16 values for a 4x4 Matrix"
);
size_t
idx
=
0
;
for
(
size_t
row
=
0
;
row
<
4
;
++
row
)
{
for
(
size_t
col
=
0
;
col
<
4
;
++
col
)
{
c
(
row
,
col
)
=
boost
::
lexical_cast
<
double
>
(
tokens
[
idx
]
);
idx
++
;
}
}
return
c
;
}
/**
...
...
@@ -555,7 +575,64 @@ namespace PROPERTY_TYPE_HELPER
*/
std
::
string
asString
(
const
WPVBaseTypes
::
PV_MATRIX4X4
&
v
)
{
return
toString
(
v
);
std
::
ostringstream
out
;
for
(
size_t
row
=
0
;
row
<
4
;
++
row
)
{
for
(
size_t
col
=
0
;
col
<
4
;
++
col
)
{
out
<<
v
(
row
,
col
)
<<
";"
;
}
}
return
out
.
str
();
}
};
/**
* Class helping to create a new instance of the property content from an old one. Selections need this special care since they contain not
* serializable content which needs to be acquired from its predecessor instance.
*/
template
<
>
class
WStringConversion
<
WPVBaseTypes
::
PV_POSITION
>
{
public:
/**
* Creates a new instance of the type from a given string. Some classes need a predecessor which is also specified here.
*
* \param str the new value as string
*
* \return the new instance
*/
WPVBaseTypes
::
PV_POSITION
create
(
const
WPVBaseTypes
::
PV_POSITION
&
/*old*/
,
const
std
::
string
str
)
{
WPVBaseTypes
::
PV_POSITION
c
;
std
::
vector
<
std
::
string
>
tokens
;
tokens
=
string_utils
::
tokenize
(
str
,
";"
);
WAssert
(
tokens
.
size
()
>=
3
,
"There weren't 3 values for a 3D vector"
);
size_t
idx
=
0
;
for
(
size_t
col
=
0
;
col
<
3
;
++
col
)
{
c
[
col
]
=
boost
::
lexical_cast
<
double
>
(
tokens
[
idx
]
);
idx
++
;
}
return
c
;
}
/**
* Creates a string from the specified value.
*
* \param v the value to convert
*
* \return the string representation
*/
std
::
string
asString
(
const
WPVBaseTypes
::
PV_POSITION
&
v
)
{
std
::
ostringstream
out
;
for
(
size_t
col
=
0
;
col
<
3
;
++
col
)
{
out
<<
v
[
col
]
<<
";"
;
}
return
out
.
str
();
}
};
}
...
...
src/common/WSharedObject.h
View file @
874e625f
...
...
@@ -83,7 +83,7 @@ public:
*
* \return the condition
*/
boost
::
shared_ptr
<
WCondition
>
getChangeCondition
();
boost
::
shared_ptr
<
WCondition
>
getChangeCondition
()
const
;
protected:
...
...
@@ -122,7 +122,7 @@ WSharedObject< T >::~WSharedObject()
}
template
<
typename
T
>
boost
::
shared_ptr
<
WCondition
>
WSharedObject
<
T
>::
getChangeCondition
()
boost
::
shared_ptr
<
WCondition
>
WSharedObject
<
T
>::
getChangeCondition
()
const
{
return
m_changeCondition
;
}
...
...
src/common/math/WLinearAlgebraFunctions.cpp
View file @
874e625f
...
...
@@ -24,17 +24,16 @@
#include <vector>
#include "../../ext/Eigen/SVD"
#include "../WAssert.h"
#include "WLinearAlgebraFunctions.h"
#include "WMatrix.h"
#include "WMatrix4x4.h"
#include "WValue.h"
#include "WVector3D.h"
#ifdef OW_USE_OSSIM
#include "WOSSIMHelper.h"
#endif
WVector3D
multMatrixWithVector3D
(
WMatrix
<
double
>
mat
,
WVector3D
vec
)
{
WVector3D
result
;
...
...
@@ -301,32 +300,32 @@ bool linearIndependent( const WVector3D& u, const WVector3D& v )
return
true
;
}
void
computeSVD
(
const
WMatrix
<
double
>
&
A
,
WMatrix
<
double
>
&
U
,
WMatrix
<
double
>
&
V
,
WV
alue
<
double
>
&
S
)
void
computeSVD
(
const
WMatrix
_2
&
A
,
WMatrix
_2
&
U
,
WMatrix
_2
&
V
,
WV
ector_2
&
S
)
{
#ifdef OW_USE_OSSIM
WOSSIMHelper
::
computeSVD
(
A
,
U
,
V
,
S
);
#else
(
void
)
A
;
(
void
)
U
;
(
void
)
V
;
(
void
)
S
;
// NOLINT to prevent "unused variable" warnings
WAssert
(
false
,
"OpenWalnut must be compiled with OSSIM to support SVD."
);
#endif
Eigen
::
JacobiSVD
<
WMatrix_2
>
svd
(
A
,
Eigen
::
ComputeFullU
|
Eigen
::
ComputeFullV
);
U
=
svd
.
matrixU
();
V
=
svd
.
matrixV
();
S
=
svd
.
singularValues
();
}
WMatrix
<
double
>
pseudoInverse
(
const
WMatrix
<
double
>
&
input
)
WMatrix
_2
pseudoInverse
(
const
WMatrix
_2
&
input
)
{
// calc pseudo inverse
WMatrix
<
double
>
U
(
input
.
getNbR
ows
(),
input
.
getNbC
ols
()
);
WMatrix
<
double
>
V
(
input
.
getNbC
ols
(),
input
.
getNbC
ols
()
);
WV
alue
<
double
>
Svec
(
input
.
getNbC
ols
()
);
WMatrix
_2
U
(
input
.
r
ows
(),
input
.
c
ols
()
);
WMatrix
_2
V
(
input
.
c
ols
(),
input
.
c
ols
()
);
WV
ector_2
Svec
(
input
.
c
ols
()
);
computeSVD
(
input
,
U
,
V
,
Svec
);
// create diagonal matrix S
WMatrix
<
double
>
S
(
input
.
getNbCols
(),
input
.
getNbCols
()
);
for
(
size_t
i
=
0
;
i
<
Svec
.
size
()
&&
i
<
S
.
getNbRows
()
&&
i
<
S
.
getNbCols
();
i
++
)
S
(
i
,
i
)
=
(
Svec
[
i
]
==
0.0
)
?
0.0
:
1.0
/
Svec
[
i
];
WMatrix_2
S
(
input
.
cols
(),
input
.
cols
()
);
S
.
setZero
();
for
(
int
i
=
0
;
i
<
Svec
.
size
()
&&
i
<
S
.
rows
()
&&
i
<
S
.
cols
();
i
++
)
{
S
(
i
,
i
)
=
(
Svec
[
i
]
==
0.0
)
?
0.0
:
1.0
/
Svec
[
i
];
}
return
WMatrix
<
double
>
(
V
*
S
*
U
.
transpose
d
()
);
return
WMatrix
_2
(
V
*
S
*
U
.
transpose
()
);
}
src/common/math/WLinearAlgebraFunctions.h
View file @
874e625f
...
...
@@ -26,6 +26,7 @@
#define WLINEARALGEBRAFUNCTIONS_H
#include "../WExportCommon.h"
#include "WMatrix.h"
namespace
osg
{
...
...
@@ -115,8 +116,7 @@ bool OWCOMMON_EXPORT linearIndependent( const WVector3D& u, const WVector3D& v )
*
* \note This function needs the OSSIM library.
*/
void
OWCOMMON_EXPORT
computeSVD
(
const
WMatrix
<
double
>&
A
,
WMatrix
<
double
>&
U
,
WMatrix
<
double
>&
V
,
WValue
<
double
>&
S
);
void
OWCOMMON_EXPORT
computeSVD
(
const
WMatrix_2
&
A
,
WMatrix_2
&
U
,
WMatrix_2
&
V
,
WVector_2
&
S
);
/**
* Calculates for a matrix the pseudo inverse.
...
...
@@ -126,6 +126,6 @@ void OWCOMMON_EXPORT computeSVD( const WMatrix< double >& A, WMatrix< double >&
* \return Inverted Matrix
*
*/
WMatrix
<
double
>
OWCOMMON_EXPORT
pseudoInverse
(
const
WMatrix
<
double
>
&
input
);
WMatrix
_2
OWCOMMON_EXPORT
pseudoInverse
(
const
WMatrix
_2
&
input
);
#endif // WLINEARALGEBRAFUNCTIONS_H
src/common/math/WMatrix.h
View file @
874e625f
...
...
@@ -31,12 +31,83 @@
#include "WValue.h"
#include "WVector3D.h"
#include "../WDefines.h"
#include "../../ext/Eigen/Core"
#include "../../ext/Eigen/LU"
/**
* A double 3 times 3 matrix. Stack-allocated. Column Major!
*
* Column Major indexes:
* [0 3 6
* 1 4 7
* 2 5 8]
* If you want to access coefficients using the operator( size_t, size_t ), the first parameter is still the row index, starting with 0.
*
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html
*/
typedef
Eigen
::
Matrix
<
double
,
3
,
3
>
WMatrix3x3_2
;
/**
* A double 4 times 4 matrix. Stack-allocated. Column Major!
*
* Column Major indexes:
* [0 4 8 12
* 1 5 9 13
* 2 6 10 14
* 3 7 11 15]
* If you want to access coefficients using the operator( size_t, size_t ), the first parameter is still the row index, starting with 0.
*
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html
*/
typedef
Eigen
::
Matrix
<
double
,
4
,
4
>
WMatrix4x4_2
;
/**
* A double matrix of dynamic size. Heap-allocated. Column Major!
* If you want to access coefficients using the operator( size_t, size_t ), the first parameter is still the row index, starting with 0.
*
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html
*/
typedef
Eigen
::
MatrixXd
WMatrix_2
;
/**
* A complex double matrix of dynamic size. Heap-allocated.
* If you want to access coefficients using the operator( size_t, size_t ), the first parameter is still the row index, starting with 0.
*
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html
*/
typedef
Eigen
::
MatrixXcd
WMatrixComplex_2
;
/**
* Converts a given WMatrix4x4_2 to an osg matrix.
*
* \param m the matrix to convert
*
* \return the converted matrix
*/
inline
osg
::
Matrixd
toOsgMatrixd
(
WMatrix4x4_2
m
)
{
osg
::
Matrixd
m2
;
for
(
size_t
row
=
0
;
row
<
4
;
++
row
)
{
for
(
size_t
col
=
0
;
col
<
4
;
++
col
)
{
m2
(
row
,
col
)
=
m
(
row
,
col
);
}
}
return
m2
;
}
/**
* Matrix template class with variable number of rows and columns.
* The access function are row-major, which means that the rows
* are the first parameter or index.
*/
template
<
typename
T
>
class
WMatrix
:
public
WValue
<
T
>
template
<
typename
T
>
class
OW_API_DEPRECATED
WMatrix
:
public
WValue
<
T
>
{
public:
/**
...
...
src/common/math/WMatrix4x4.h
View file @
874e625f
...
...
@@ -36,10 +36,12 @@
#include "../WAssert.h"
#include "../WStringUtils.h"
#include "../WDefines.h"
/**
* Use osg 4x4 matrices as WMatrix4x4
*/
typedef
osg
::
Matrixd
WMatrix4x4
;
OW_API_DEPRECATED
typedef
osg
::
Matrixd
WMatrix4x4
;
/**
* Write a 4x4 matrix in string representation.
...
...
@@ -48,7 +50,7 @@ typedef osg::Matrixd WMatrix4x4;
*
* \return the matrix as string
*/
inline
std
::
string
toString
(
const
WMatrix4x4
&
c
)
/*
inline std::string toString( const WMatrix4x4& c )
{
std::ostringstream out;
for ( size_t row = 0; row < 4; ++row )
...
...
@@ -60,7 +62,7 @@ inline std::string toString( const WMatrix4x4& c )
}
return out.str();
}
*/
/**
* Read a 4x4 matrix in string representation from the given string.
*
...
...
@@ -68,7 +70,7 @@ inline std::string toString( const WMatrix4x4& c )
*
* \return the matrix
*/
inline
WMatrix4x4
fromString
(
std
::
string
str
)
/*
inline WMatrix4x4 fromString( std::string str )
{
WMatrix4x4 c;
std::vector< std::string > tokens;
...
...
@@ -86,6 +88,6 @@ inline WMatrix4x4 fromString( std::string str )
}
return c;
}
}
*/
#endif // WMATRIX4X4_H
src/common/math/WSymmetricSphericalHarmonic.cpp
View file @
874e625f
...
...
@@ -34,15 +34,16 @@
#include "WSymmetricSphericalHarmonic.h"
#include "WTensorSym.h"
#include "WUnitSphereCoordinates.h"
#include "WV
alue
.h"
#include "WV
ector3D
.h"
WSymmetricSphericalHarmonic
::
WSymmetricSphericalHarmonic
()
:
m_order
(
0
),
m_SHCoefficients
(
0
)
m_order
(
0
)
// ,
// m_SHCoefficients( 0 )
{
}
WSymmetricSphericalHarmonic
::
WSymmetricSphericalHarmonic
(
const
WV
alue
<
double
>
&
SHCoefficients
)
:
WSymmetricSphericalHarmonic
::
WSymmetricSphericalHarmonic
(
const
WV
ector_2
&
SHCoefficients
)
:
m_SHCoefficients
(
SHCoefficients
)
{
// calculate order from SHCoefficients.size:
...
...
@@ -57,28 +58,28 @@ WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
double
WSymmetricSphericalHarmonic
::
getValue
(
double
theta
,
double
phi
)
const
{
double
result
=
0.0
;
int
j
=
0
;
const
double
rootOf2
=
std
::
sqrt
(
2.0
);
for
(
int
k
=
0
;
k
<=
static_cast
<
int
>
(
m_order
);
k
+=
2
)
{
// m = 1 .. k
for
(
int
m
=
1
;
m
<=
k
;
m
++
)
double
result
=
0.0
;
int
j
=
0
;
const
double
rootOf2
=
std
::
sqrt
(
2.0
);
for
(
int
k
=
0
;
k
<=
static_cast
<
int
>
(
m_order
);
k
+=
2
)
{
j
=
(
k
*
k
+
k
+
2
)
/
2
+
m
;
result
+=
m_SHCoefficients
[
j
-
1
]
*
rootOf2
*
static_cast
<
double
>
(
std
::
pow
(
-
1.0
,
m
+
1
)
)
*
boost
::
math
::
spherical_harmonic_i
(
k
,
m
,
theta
,
phi
);
}
// m = 0
result
+=
m_SHCoefficients
[
(
k
*
k
+
k
+
2
)
/
2
-
1
]
*
boost
::
math
::
spherical_harmonic_r
(
k
,
0
,
theta
,
phi
);
// m = -k .. -1
for
(
int
m
=
-
k
;
m
<
0
;
m
++
)
{
j
=
(
k
*
k
+
k
+
2
)
/
2
+
m
;
result
+=
m_SHCoefficients
[
j
-
1
]
*
rootOf2
*
boost
::
math
::
spherical_harmonic_r
(
k
,
-
m
,
theta
,
phi
);
// m = 1 .. k
for
(
int
m
=
1
;
m
<=
k
;
m
++
)
{
j
=
(
k
*
k
+
k
+
2
)
/
2
+
m
;
result
+=
m_SHCoefficients
[
j
-
1
]
*
rootOf2
*
static_cast
<
double
>
(
std
::
pow
(
-
1.0
,
m
+
1
)
)
*
boost
::
math
::
spherical_harmonic_i
(
k
,
m
,
theta
,
phi
);
}
// m = 0
result
+=
m_SHCoefficients
[
(
k
*
k
+
k
+
2
)
/
2
-
1
]
*
boost
::
math
::
spherical_harmonic_r
(
k
,
0
,
theta
,
phi
);
// m = -k .. -1
for
(
int
m
=
-
k
;
m
<
0
;
m
++
)
{
j
=
(
k
*
k
+
k
+
2
)
/
2
+
m
;
result
+=
m_SHCoefficients
[
j
-
1
]
*
rootOf2
*
boost
::
math
::
spherical_harmonic_r
(
k
,
-
m
,
theta
,
phi
);
}
}
}
return
result
;
return
result
;
}
double
WSymmetricSphericalHarmonic
::
getValue
(
const
WUnitSphereCoordinates
&
coordinates
)
const
...
...
@@ -86,14 +87,14 @@ double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coor
return
getValue
(
coordinates
.
getTheta
(),
coordinates
.
getPhi
()
);
}
const
WV
alue
<
double
>
&
WSymmetricSphericalHarmonic
::
getCoefficients
()
const
const
WV
ector_2
&
WSymmetricSphericalHarmonic
::
getCoefficients
()
const
{
return
m_SHCoefficients
;
}
WV
alue
<
double
>
WSymmetricSphericalHarmonic
::
getCoefficientsSchultz
()
const
WV
ector_2
WSymmetricSphericalHarmonic
::
getCoefficientsSchultz
()
const
{
WV
alue
<
double
>
res
(
m_SHCoefficients
.
size
()
);
WV
ector_2
res
(
m_SHCoefficients
.
size
()
);
size_t
k
=
0
;
for
(
int
l
=
0
;
l
<=
static_cast
<
int
>
(
m_order
);
l
+=
2
)
{
...
...
@@ -114,9 +115,9 @@ WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
return
res
;
}
WV
alue
<
std
::
complex
<
double
>
>
WSymmetricSphericalHarmonic
::
getCoefficientsComplex
()
const
WV
ectorComplex_2
WSymmetricSphericalHarmonic
::
getCoefficientsComplex
()
const
{
WV
alue
<
std
::
complex
<
double
>
>
res
(
m_SHCoefficients
.
size
()
);
WV
ectorComplex_2
res
(
m_SHCoefficients
.
size
()
);
size_t
k
=
0
;
double
r
=
1.0
/
sqrt
(
2.0
);
std
::
complex
<
double
>
i
(
0.0
,
-
1.0
);
...
...
@@ -187,30 +188,30 @@ double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates
return
gfa
;
}
double
WSymmetricSphericalHarmonic
::
calcGFA
(
WMatrix
<
double
>
const
&
B
)
const
double
WSymmetricSphericalHarmonic
::
calcGFA
(
const
WMatrix
_2
&
B
)
const
{
// sh coeffs
WMatrix
<
double
>
s
(
B
.
getNbC
ols
(),
1
);
WAssert
(
B
.
getNbC
ols
()
==
m_SHCoefficients
.
size
(),
""
);
for
(
std
::
size_
t
k
=
0
;
k
<
B
.
getNbC
ols
();
++
k
)
WMatrix
_2
s
(
B
.
c
ols
(),
1
);
WAssert
(
B
.
c
ols
()
==
m_SHCoefficients
.
size
(),
""
);
for
(
in
t
k
=
0
;
k
<
B
.
c
ols
();
++
k
)
{
s
(
k
,
0
)
=
m_SHCoefficients
[
k
];
}
s
=
B
*
s
;
WAssert
(
s
.
getNbR
ows
()
==
B
.
getNbR
ows
(),
""
);
WAssert
(
s
.
getNbC
ols
()
==
1u
,
""
);
WAssert
(
s
.
r
ows
()
==
B
.
r
ows
(),
""
);
WAssert
(
s
.
c
ols
()
==
1u
,
""
);
double
n
=
static_cast
<
double
>
(
s
.
getNbR
ows
()
);
double
n
=
static_cast
<
double
>
(
s
.
r
ows
()
);
double
d
=
0.0
;
double
gfa
=
0.0
;
double
mean
=
0.0
;
for
(
std
::
size_
t
i
=
0
;
i
<
s
.
getNbR
ows
();
++
i
)
for
(
in
t
i
=
0
;
i
<
s
.
r
ows
();
++
i
)
{
mean
+=
s
(
i
,
0
);
}
mean
/=
n
;
for
(
std
::
size_
t
i
=
0
;
i
<
s
.
getNbR
ows
();
++
i
)
for
(
in
t
i
=
0
;
i
<
s
.
r
ows
();
++
i
)
{
double
f
=
s
(
i
,
0
);
// we use gfa as a temporary here
...
...
@@ -238,12 +239,12 @@ double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
return
gfa
;
}
void
WSymmetricSphericalHarmonic
::
applyFunkRadonTransformation
(
WMatrix
<
double
>
const
&
frtMat
)
void
WSymmetricSphericalHarmonic
::
applyFunkRadonTransformation
(
const
WMatrix
_2
&
frtMat
)
{
WAssert
(
frtMat
.
getNbC
ols
()
==
m_SHCoefficients
.
size
(),
""
);
WAssert
(
frtMat
.
getNbR
ows
()
==
m_SHCoefficients
.
size
(),
""
);
WAssert
(
frtMat
.
c
ols
()
==
m_SHCoefficients
.
size
(),
""
);
WAssert
(
frtMat
.
r
ows
()
==
m_SHCoefficients
.
size
(),
""
);
// Funk-Radon-Transformation as in Descoteaux's thesis
for
(
size_
t
j
=
0
;
j
<
m_SHCoefficients
.
size
();
j
++
)
for
(
in
t
j
=
0
;
j
<
m_SHCoefficients
.
size
();
j
++
)
{
m_SHCoefficients
[
j
]
*=
frtMat
(
j
,
j
);
}
...
...
@@ -254,7 +255,7 @@ size_t WSymmetricSphericalHarmonic::getOrder() const
return
m_order
;
}
WMatrix
<
double
>
WSymmetricSphericalHarmonic
::
getSHFittingMatrix
(
const
std
::
vector
<
WVector3D
>&
orientations
,
WMatrix
_2
WSymmetricSphericalHarmonic
::
getSHFittingMatrix
(
const
std
::
vector
<
WVector3D
>&
orientations
,
int
order
,
double
lambda
,
bool
withFRT
)
...
...
@@ -268,33 +269,34 @@ WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vect
return
WSymmetricSphericalHarmonic
::
getSHFittingMatrix
(
sphereCoordinates
,
order
,
lambda
,
withFRT
);
}
WMatrix
<
double
>
WSymmetricSphericalHarmonic
::
getSHFittingMatrix
(
const
std
::
vector
<
WUnitSphereCoordinates
>&
orientations
,
int
order
,
double
lambda
,
bool
withFRT
)
WMatrix
_2
WSymmetricSphericalHarmonic
::
getSHFittingMatrix
(
const
std
::
vector
<
WUnitSphereCoordinates
>&
orientations
,
int
order
,
double
lambda
,
bool
withFRT
)
{
WMatrix
<
double
>
B
(
WSymmetricSphericalHarmonic
::
calcBaseMatrix
(
orientations
,
order
)
);
WMatrix
<
double
>
Bt
(
B
.
transpose
d
()
);
WMatrix
<
double
>
result
(
Bt
*
B
);
WMatrix
_2
B
(
WSymmetricSphericalHarmonic
::
calcBaseMatrix
(
orientations
,
order
)
);
WMatrix
_2
Bt
(
B
.
transpose
()
);
WMatrix
_2
result
(
Bt
*
B
);
if
(
lambda
!=
0.0
)
{
WMatrix
<
double
>
L
(
WSymmetricSphericalHarmonic
::
calcSmoothingMatrix
(
order
)
);
WMatrix
_2
L
(
WSymmetricSphericalHarmonic
::
calcSmoothingMatrix
(
order
)
);
result
+=
lambda
*
L
;
}
result
=
pseudoInverse
(
result
)
*
Bt
;
// result *= Bt;
if
(
withFRT
)
{
WMatrix
<
double
>
P
(
WSymmetricSphericalHarmonic
::
calcFRTMatrix
(
order
)
);
WMatrix
_2
P
(
WSymmetricSphericalHarmonic
::
calcFRTMatrix
(
order
)
);
return
(
P
*
result
);
}
return
result
;
}
WMatrix
<
double
>
WSymmetricSphericalHarmonic
::
calcBaseMatrix
(
const
std
::
vector
<
WUnitSphereCoordinates
>&
orientations
,
WMatrix
_2
WSymmetricSphericalHarmonic
::
calcBaseMatrix
(
const
std
::
vector
<
WUnitSphereCoordinates
>&
orientations
,
int
order
)
{
WMatrix
<
double
>
B
(
orientations
.
size
(),
(
(
order
+
1
)
*
(
order
+
2
)
)
/
2
);
WMatrix
_2
B
(
orientations
.
size
(),
(
(
order
+
1
)
*
(
order
+
2
)
)
/
2
);
// calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
int
j
=
0
;
...
...
@@ -325,10 +327,10 @@ WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector<
return
B
;
}
WMatrix
<
std
::
complex
<
double
>
>
WMatrix
Complex_2
WSymmetricSphericalHarmonic
::
calcComplexBaseMatrix
(
std
::
vector
<
WUnitSphereCoordinates
>
const
&
orientations
,
int
order
)
{
WMatrix
<
std
::
complex
<
double
>
>
B
(
orientations
.
size
(),
(
(
order
+
1
)
*
(
order
+
2
)
)
/
2
);
WMatrix
Complex_2
B
(
orientations
.
size
(),
(
(
order
+
1
)
*
(
order
+
2
)
)
/
2
);
for
(
std
::
size_t
row
=
0
;
row
<
orientations
.
size
();
row
++
)
{
...
...
@@ -355,11 +357,11 @@ WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoor
return
B
;
}
WMatrix
<
double
>
WSymmetricSphericalHarmonic
::
calcSmoothingMatrix
(
size_t
order
)