use a 3x3 eigendecomposition in gridvolume::lookupVector()
Wenzel Jakob
parent
0ca137b4b2
commit
1eabd9ce9a
|
|
@ -678,6 +678,17 @@ template <typename T, int M1, int N1, int M2, int N2> inline Matrix<M1, N2, T>
|
|||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* \brief Fast 3x3 eigenvalue decomposition
|
||||
* \param m
|
||||
* Matrix in question -- will be replaced with the eigenvectors
|
||||
* \param lambda
|
||||
* Parameter used to returns the eigenvalues
|
||||
* \return
|
||||
* \c true upon success.
|
||||
*/
|
||||
extern MTS_EXPORT_CORE bool eig3(Matrix3x3 &m, Float lambda[3]);
|
||||
|
||||
MTS_NAMESPACE_END
|
||||
|
||||
#include "matrix.inl"
|
||||
|
|
|
|||
|
|
@ -206,4 +206,150 @@ std::string Transform::toString() const {
|
|||
return m_transform.toString();
|
||||
}
|
||||
|
||||
/**
|
||||
* \brief Tridiagonalizes a matrix using a Householder transformation
|
||||
* and returns the result.
|
||||
*
|
||||
* Based on "Geometric Tools" by David Eberly.
|
||||
*/
|
||||
static inline void tred3(Matrix3x3 &m, Float *diag, Float *subd) {
|
||||
Float m00 = m(0, 0), m01 = m(0, 1), m02 = m(0, 2),
|
||||
m11 = m(1, 1), m12 = m(1, 2), m22 = m(2, 2);
|
||||
|
||||
diag[0] = m00;
|
||||
subd[2] = 0;
|
||||
if (std::abs(m02) > std::numeric_limits<Float>::epsilon()) {
|
||||
Float length = std::sqrt(m01*m01 + m02*m02),
|
||||
invLength = 1 / length;
|
||||
m01 *= invLength;
|
||||
m02 *= invLength;
|
||||
Float q = 2*m01*m12 + m02*(m22 - m11);
|
||||
diag[1] = m11 + m02*q;
|
||||
diag[2] = m22 - m02*q;
|
||||
subd[0] = length;
|
||||
subd[1] = m12 - m01*q;
|
||||
m = Matrix3x3(
|
||||
1, 0, 0,
|
||||
0, m01, m02,
|
||||
0, m02, -m01);
|
||||
} else {
|
||||
/* The matrix is already tridiagonal,
|
||||
return an identity transformation matrix */
|
||||
diag[1] = m11;
|
||||
diag[2] = m22;
|
||||
subd[0] = m01;
|
||||
subd[1] = m12;
|
||||
m = Matrix3x3(
|
||||
1, 0, 0,
|
||||
0, 1, 0,
|
||||
0, 0, 1);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the eigenvalues and eigenvectors of a 3x3 symmetric tridiagonal
|
||||
* matrix using the QL algorithm with implicit shifts
|
||||
*
|
||||
* Based on "Geometric Tools" by David Eberly and Numerical Recipes by
|
||||
* Teukolsky, et al. Originally, this code goes back to the Handbook
|
||||
* for Automatic Computation by Wilkionson and Reinsch, as well as
|
||||
* the corresponding EISPACK routine.
|
||||
*/
|
||||
|
||||
static inline bool ql3(Matrix3x3 &m, Float *diag, Float *subd) {
|
||||
const int maxIter = 32;
|
||||
|
||||
for (int i = 0; i < 3; ++i) {
|
||||
int k;
|
||||
for (k = 0; k < maxIter; ++k) {
|
||||
int j;
|
||||
for (j = i; j < 2; ++j) {
|
||||
Float tmp = std::abs(diag[j]) + std::abs(diag[j+1]);
|
||||
if (std::abs(subd[j]) + tmp == tmp)
|
||||
break;
|
||||
}
|
||||
if (j == i)
|
||||
break;
|
||||
|
||||
Float value0 = (diag[i + 1] - diag[i])/(2*subd[i]);
|
||||
Float value1 = std::sqrt(value0*value0 + 1);
|
||||
value0 = diag[j] - diag[i] + subd[i] /
|
||||
((value0 < 0) ? (value0 - value1) : (value0 + value1));
|
||||
|
||||
Float sn = 1, cs = 1, value2 = 0;
|
||||
for (int l = j - 1; l >= i; --l) {
|
||||
Float value3 = sn*subd[l], value4 = cs*subd[l];
|
||||
if (std::abs(value3) >= std::abs(value0)) {
|
||||
cs = value0 / value3;
|
||||
value1 = std::sqrt(cs*cs + 1);
|
||||
subd[l + 1] = value3*value1;
|
||||
sn = 1.0f / value1;
|
||||
cs *= sn;
|
||||
} else {
|
||||
sn = value3 / value0;
|
||||
value1 = std::sqrt(sn*sn + 1);
|
||||
subd[l + 1] = value0*value1;
|
||||
cs = 1.0f / value1;
|
||||
sn *= cs;
|
||||
}
|
||||
value0 = diag[l + 1] - value2;
|
||||
value1 = (diag[l] - value0)*sn + 2*value4*cs;
|
||||
value2 = sn*value1;
|
||||
diag[l + 1] = value0 + value2;
|
||||
value0 = cs*value1 - value4;
|
||||
|
||||
for (int t = 0; t < 3; ++t) {
|
||||
value3 = m(t, l + 1);
|
||||
m(t, l + 1) = sn*m(t, l) + cs*value3;
|
||||
m(t, l) = cs*m(t, l) - sn*value3;
|
||||
}
|
||||
}
|
||||
diag[i] -= value2;
|
||||
subd[i] = value0;
|
||||
subd[j] = 0;
|
||||
}
|
||||
if (k == maxIter)
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/// Fast 3x3 eigenvalue decomposition
|
||||
bool eig3(Matrix3x3 &m, Float lambda[3]) {
|
||||
Float subd[3];
|
||||
|
||||
/* Reduce to Hessenberg form */
|
||||
tred3(m, lambda, subd);
|
||||
|
||||
/* Iteratively find the eigenvalues and eigenvectors
|
||||
using implicitly shifted QL iterations */
|
||||
if (!ql3(m, lambda, subd))
|
||||
return false;
|
||||
|
||||
/* Sort the eigenvalues in decreasing order */
|
||||
for (int i=0; i<2; ++i) {
|
||||
/* Locate the maximum eigenvalue. */
|
||||
int largest = i;
|
||||
Float maxValue = lambda[largest];
|
||||
for (int j = i+1; j<3; ++j) {
|
||||
if (lambda[j] > maxValue) {
|
||||
largest = j;
|
||||
maxValue = lambda[largest];
|
||||
}
|
||||
}
|
||||
|
||||
if (largest != i) {
|
||||
// Swap the eigenvalues.
|
||||
std::swap(lambda[i], lambda[largest]);
|
||||
|
||||
// Swap the eigenvectors
|
||||
for (int j=0; j<3; ++j)
|
||||
std::swap(m(j, i), m(j, largest));
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
MTS_NAMESPACE_END
|
||||
|
|
|
|||
|
|
@ -289,6 +289,10 @@ public:
|
|||
inline Vector toVector() const {
|
||||
return Vector(value[0], value[1], value[2]);
|
||||
}
|
||||
|
||||
float operator[](int i) const {
|
||||
return value[i];
|
||||
}
|
||||
|
||||
inline Matrix3x3 tensor() const {
|
||||
return Matrix3x3(
|
||||
|
|
@ -466,30 +470,30 @@ public:
|
|||
}
|
||||
#else
|
||||
Float _fx = 1.0f - fx, _fy = 1.0f - fy, _fz = 1.0f - fz;
|
||||
float3 d000, d001, d010, d011, d100, d101, d110, d111;
|
||||
float3 d[8];
|
||||
|
||||
switch (m_volumeType) {
|
||||
case EFloat32: {
|
||||
const float3 *vectorData = (float3 *) m_data;
|
||||
d000 = vectorData[(z1*m_res.y + y1)*m_res.x + x1];
|
||||
d001 = vectorData[(z1*m_res.y + y1)*m_res.x + x2];
|
||||
d010 = vectorData[(z1*m_res.y + y2)*m_res.x + x1];
|
||||
d011 = vectorData[(z1*m_res.y + y2)*m_res.x + x2];
|
||||
d100 = vectorData[(z2*m_res.y + y1)*m_res.x + x1];
|
||||
d101 = vectorData[(z2*m_res.y + y1)*m_res.x + x2];
|
||||
d110 = vectorData[(z2*m_res.y + y2)*m_res.x + x1];
|
||||
d111 = vectorData[(z2*m_res.y + y2)*m_res.x + x2];
|
||||
d[0] = vectorData[(z1*m_res.y + y1)*m_res.x + x1];
|
||||
d[1] = vectorData[(z1*m_res.y + y1)*m_res.x + x2];
|
||||
d[2] = vectorData[(z1*m_res.y + y2)*m_res.x + x1];
|
||||
d[3] = vectorData[(z1*m_res.y + y2)*m_res.x + x2];
|
||||
d[4] = vectorData[(z2*m_res.y + y1)*m_res.x + x1];
|
||||
d[5] = vectorData[(z2*m_res.y + y1)*m_res.x + x2];
|
||||
d[6] = vectorData[(z2*m_res.y + y2)*m_res.x + x1];
|
||||
d[7] = vectorData[(z2*m_res.y + y2)*m_res.x + x2];
|
||||
}
|
||||
break;
|
||||
case EQuantizedDirections: {
|
||||
d000 = lookupQuantizedDirection((z1*m_res.y + y1)*m_res.x + x1);
|
||||
d001 = lookupQuantizedDirection((z1*m_res.y + y1)*m_res.x + x2);
|
||||
d010 = lookupQuantizedDirection((z1*m_res.y + y2)*m_res.x + x1);
|
||||
d011 = lookupQuantizedDirection((z1*m_res.y + y2)*m_res.x + x2);
|
||||
d100 = lookupQuantizedDirection((z2*m_res.y + y1)*m_res.x + x1);
|
||||
d101 = lookupQuantizedDirection((z2*m_res.y + y1)*m_res.x + x2);
|
||||
d110 = lookupQuantizedDirection((z2*m_res.y + y2)*m_res.x + x1);
|
||||
d111 = lookupQuantizedDirection((z2*m_res.y + y2)*m_res.x + x2);
|
||||
d[0] = lookupQuantizedDirection((z1*m_res.y + y1)*m_res.x + x1);
|
||||
d[1] = lookupQuantizedDirection((z1*m_res.y + y1)*m_res.x + x2);
|
||||
d[2] = lookupQuantizedDirection((z1*m_res.y + y2)*m_res.x + x1);
|
||||
d[3] = lookupQuantizedDirection((z1*m_res.y + y2)*m_res.x + x2);
|
||||
d[4] = lookupQuantizedDirection((z2*m_res.y + y1)*m_res.x + x1);
|
||||
d[5] = lookupQuantizedDirection((z2*m_res.y + y1)*m_res.x + x2);
|
||||
d[6] = lookupQuantizedDirection((z2*m_res.y + y2)*m_res.x + x1);
|
||||
d[7] = lookupQuantizedDirection((z2*m_res.y + y2)*m_res.x + x2);
|
||||
}
|
||||
break;
|
||||
default:
|
||||
|
|
@ -497,17 +501,29 @@ public:
|
|||
}
|
||||
|
||||
#if defined(VINTERP_LINEAR)
|
||||
value = (((d000*_fx + d001*fx)*_fy +
|
||||
(d010*_fx + d011*fx)*fy)*_fz +
|
||||
((d100*_fx + d101*fx)*_fy +
|
||||
(d110*_fx + d111*fx)*fy)*fz).toVector();
|
||||
value = (((d[0]*_fx + d[1]*fx)*_fy +
|
||||
(d[2]*_fx + d[3]*fx)* fy)*_fz +
|
||||
((d[4]*_fx + d[5]*fx)*_fy +
|
||||
(d[6]*_fx + d[7]*fx)* fy)* fz).toVector();
|
||||
#elif defined(VINTERP_STRUCTURE_TENSOR)
|
||||
Matrix3x3 tensor =
|
||||
((d000.tensor()*_fx + d001.tensor()*fx)*_fy +
|
||||
(d010.tensor()*_fx + d011.tensor()*fx)*fy)*_fz +
|
||||
((d100.tensor()*_fx + d101.tensor()*fx)*_fy +
|
||||
(d110.tensor()*_fx + d111.tensor()*fx)*fy)*fz;
|
||||
|
||||
Matrix3x3 tensor(0.0f);
|
||||
for (int k=0; k<8; ++k) {
|
||||
Float factor = ((k & 1) ? fx : _fx) * ((k & 2) ? fy : _fy)
|
||||
* ((k & 4) ? fz : _fz);
|
||||
for (int i=0; i<3; ++i)
|
||||
for (int j=0; j<3; ++j)
|
||||
tensor(i, j) += factor * d[k][i] * d[k][j];
|
||||
}
|
||||
|
||||
Float lambda[3];
|
||||
Matrix3x3 Q(tensor);
|
||||
if (!eig3(Q, lambda)) {
|
||||
Log(EWarn, "lookupVector(): Eigendecomposition failed!");
|
||||
return Vector(0.0f);
|
||||
}
|
||||
// Float specularity = 1-lambda[1]/lambda[0];
|
||||
value = Q.col(0);
|
||||
#if 0
|
||||
if (tensor.isZero())
|
||||
return Vector(0.0f);
|
||||
|
||||
|
|
@ -522,6 +538,7 @@ public:
|
|||
for (int i=0; i<POWER_ITERATION_STEPS-1; ++i)
|
||||
value = normalize(tensor * value);
|
||||
value = tensor * value;
|
||||
#endif
|
||||
#else
|
||||
#error Need to choose a vector interpolation method!
|
||||
#endif
|
||||
|
|
@ -532,7 +549,7 @@ public:
|
|||
else
|
||||
return Vector(0.0f);
|
||||
}
|
||||
|
||||
|
||||
bool supportsFloatLookups() const { return m_channels == 1; }
|
||||
bool supportsSpectrumLookups() const { return m_channels == 3; }
|
||||
bool supportsVectorLookups() const { return m_channels == 3; }
|
||||
|
|
|
|||
Loading…
Reference in New Issue