spacecruft
/
CelestiaContent
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Move ellipsoidTangent to mathlib.h

pull/3/head
Hleb Valoshka 2019-11-08 20:29:27 +03:00
parent fecd61fcb1
commit b259987e95
3 changed files with 51 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
src/celengine/render.cpp
View File

@ -3746,57 +3746,6 @@ struct LightIrradiancePredicate
};
template <typename T> static Matrix<T, 3, 1>
ellipsoidTangent(const Matrix<T, 3, 1>& recipSemiAxes,
const Matrix<T, 3, 1>& w,
const Matrix<T, 3, 1>& e,
const Matrix<T, 3, 1>& e_,
T ee)
{
// We want to find t such that -E(1-t) + Wt is the direction of a ray
// tangent to the ellipsoid. A tangent ray will intersect the ellipsoid
// at exactly one point. Finding the intersection between a ray and an
// ellipsoid ultimately requires using the quadratic formula, which has
// one solution when the discriminant (b^2 - 4ac) is zero. The code below
// computes the value of t that results in a discriminant of zero.
Matrix<T, 3, 1> w_ = w.cwiseProduct(recipSemiAxes);//(w.x * recipSemiAxes.x, w.y * recipSemiAxes.y, w.z * recipSemiAxes.z);
T ww = w_.dot(w_);
T ew = w_.dot(e_);
// Before elimination of terms:
// double a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1.0f);
// double b = -8 * ee * (ee + ew) - 4 * (-2 * (ee + ew) * (ee - 1.0f));
// double c = 4 * ee * ee - 4 * (ee * (ee - 1.0f));
// Simplify the below expression and eliminate the ee^2 terms; this
// prevents precision errors, as ee tends to be a very large value.
//T a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1);
T a = 4 * (square(ew) - ee * ww + ee + 2 * ew + ww);
T b = -8 * (ee + ew);
T c = 4 * ee;
T t = 0;
T discriminant = b * b - 4 * a * c;
if (discriminant < 0)
t = (-b + (T) sqrt(-discriminant)) / (2 * a); // Bad!
else
t = (-b + (T) sqrt(discriminant)) / (2 * a);
// V is the direction vector. We now need the point of intersection,
// which we obtain by solving the quadratic equation for the ray-ellipse
// intersection. Since we already know that the discriminant is zero,
// the solution is just -b/2a
Matrix<T, 3, 1> v = -e * (1 - t) + w * t;
Matrix<T, 3, 1> v_ = v.cwiseProduct(recipSemiAxes);
T a1 = v_.dot(v_);
T b1 = (T) 2 * v_.dot(e_);
T t1 = -b1 / (2 * a1);
return e + v * t1;
}
void Renderer::renderEllipsoidAtmosphere(const Atmosphere& atmosphere,
const Vector3f& center,
const Quaternionf& orientation,

50
src/celengine/visibleregion.cpp
View File

@ -81,56 +81,6 @@ VisibleRegion::setOpacity(float opacity)
}
template <typename T> static Matrix<T, 3, 1>
ellipsoidTangent(const Matrix<T, 3, 1>& recipSemiAxes,
const Matrix<T, 3, 1>& w,
const Matrix<T, 3, 1>& e,
const Matrix<T, 3, 1>& e_,
T ee)
{
// We want to find t such that -E(1-t) + Wt is the direction of a ray
// tangent to the ellipsoid. A tangent ray will intersect the ellipsoid
// at exactly one point. Finding the intersection between a ray and an
// ellipsoid ultimately requires using the quadratic formula, which has
// one solution when the discriminant (b^2 - 4ac) is zero. The code below
// computes the value of t that results in a discriminant of zero.
Matrix<T, 3, 1> w_ = w.cwiseProduct(recipSemiAxes);//(w.x * recipSemiAxes.x, w.y * recipSemiAxes.y, w.z * recipSemiAxes.z);
T ww = w_.dot(w_);
T ew = w_.dot(e_);
// Before elimination of terms:
// double a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1.0f);
// double b = -8 * ee * (ee + ew) - 4 * (-2 * (ee + ew) * (ee - 1.0f));
// double c = 4 * ee * ee - 4 * (ee * (ee - 1.0f));
// Simplify the below expression and eliminate the ee^2 terms; this
// prevents precision errors, as ee tends to be a very large value.
//T a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1);
T a = 4 * (square(ew) - ee * ww + ee + 2 * ew + ww);
T b = -8 * (ee + ew);
T c = 4 * ee;
T t = 0;
T discriminant = b * b - 4 * a * c;
if (discriminant < 0)
t = (-b + (T) sqrt(-discriminant)) / (2 * a); // Bad!
else
t = (-b + (T) sqrt(discriminant)) / (2 * a);
// V is the direction vector. We now need the point of intersection,
// which we obtain by solving the quadratic equation for the ray-ellipse
// intersection. Since we already know that the discriminant is zero,
// the solution is just -b/2a
Matrix<T, 3, 1> v = -e * (1 - t) + w * t;
Matrix<T, 3, 1> v_ = v.cwiseProduct(recipSemiAxes);
T a1 = v_.dot(v_);
T b1 = (T) 2 * v_.dot(e_);
T t1 = -b1 / (2 * a1);
return e + v * t1;
}
constexpr const unsigned maxSections = 360;
static void

51
src/celmath/mathlib.h
View File

@ -12,6 +12,7 @@
#include <cmath>
#include <cstdlib>
#include <Eigen/Core>
#define PI 3.14159265358979323846
@ -99,5 +100,55 @@ template<typename T> inline constexpr T sphereArea(T r)
{
return 4 * static_cast<T>(PI) * r * r;
}
template <typename T> static Eigen::Matrix<T, 3, 1>
ellipsoidTangent(const Eigen::Matrix<T, 3, 1>& recipSemiAxes,
const Eigen::Matrix<T, 3, 1>& w,
const Eigen::Matrix<T, 3, 1>& e,
const Eigen::Matrix<T, 3, 1>& e_,
T ee)
{
// We want to find t such that -E(1-t) + Wt is the direction of a ray
// tangent to the ellipsoid. A tangent ray will intersect the ellipsoid
// at exactly one point. Finding the intersection between a ray and an
// ellipsoid ultimately requires using the quadratic formula, which has
// one solution when the discriminant (b^2 - 4ac) is zero. The code below
// computes the value of t that results in a discriminant of zero.
Eigen::Matrix<T, 3, 1> w_ = w.cwiseProduct(recipSemiAxes);
T ww = w_.dot(w_);
T ew = w_.dot(e_);
// Before elimination of terms:
// double a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1.0f);
// double b = -8 * ee * (ee + ew) - 4 * (-2 * (ee + ew) * (ee - 1.0f));
// double c = 4 * ee * ee - 4 * (ee * (ee - 1.0f));
// Simplify the below expression and eliminate the ee^2 terms; this
// prevents precision errors, as ee tends to be a very large value.
//T a = 4 * square(ee + ew) - 4 * (ee + 2 * ew + ww) * (ee - 1);
T a = 4 * (square(ew) - ee * ww + ee + 2 * ew + ww);
T b = -8 * (ee + ew);
T c = 4 * ee;
T t = 0;
T discriminant = b * b - 4 * a * c;
if (discriminant < 0)
discriminant = -discriminant; // Bad!
t = (-b + (T) sqrt(discriminant)) / (2 * a);
// V is the direction vector. We now need the point of intersection,
// which we obtain by solving the quadratic equation for the ray-ellipse
// intersection. Since we already know that the discriminant is zero,
// the solution is just -b/2a
Eigen::Matrix<T, 3, 1> v = -e * (1 - t) + w * t;
Eigen::Matrix<T, 3, 1> v_ = v.cwiseProduct(recipSemiAxes);
T a1 = v_.dot(v_);
T b1 = (T) 2 * v_.dot(e_);
T t1 = -b1 / (2 * a1);
return e + v * t1;
}
}; // namespace celmath
#endif // _MATHLIB_H_