From b259987e951af992ee3dcd6bf99cefcccfdb53d6 Mon Sep 17 00:00:00 2001 From: Hleb Valoshka <375gnu@gmail.com> Date: Fri, 8 Nov 2019 20:29:27 +0300 Subject: [PATCH] Move ellipsoidTangent to mathlib.h --- src/celengine/render.cpp | 51 --------------------------------- src/celengine/visibleregion.cpp | 50 -------------------------------- src/celmath/mathlib.h | 51 +++++++++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 101 deletions(-) diff --git a/src/celengine/render.cpp b/src/celengine/render.cpp index 6657da56..587a72da 100644 --- a/src/celengine/render.cpp +++ b/src/celengine/render.cpp @@ -3746,57 +3746,6 @@ struct LightIrradiancePredicate }; -template static Matrix -ellipsoidTangent(const Matrix& recipSemiAxes, - const Matrix& w, - const Matrix& e, - const Matrix& 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 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 v = -e * (1 - t) + w * t; - Matrix 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, diff --git a/src/celengine/visibleregion.cpp b/src/celengine/visibleregion.cpp index 39764125..9576c83d 100644 --- a/src/celengine/visibleregion.cpp +++ b/src/celengine/visibleregion.cpp @@ -81,56 +81,6 @@ VisibleRegion::setOpacity(float opacity) } -template static Matrix -ellipsoidTangent(const Matrix& recipSemiAxes, - const Matrix& w, - const Matrix& e, - const Matrix& 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 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 v = -e * (1 - t) + w * t; - Matrix 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 diff --git a/src/celmath/mathlib.h b/src/celmath/mathlib.h index d12bc007..a1f73d5c 100644 --- a/src/celmath/mathlib.h +++ b/src/celmath/mathlib.h @@ -12,6 +12,7 @@ #include #include +#include #define PI 3.14159265358979323846 @@ -99,5 +100,55 @@ template inline constexpr T sphereArea(T r) { return 4 * static_cast(PI) * r * r; } + +template static Eigen::Matrix +ellipsoidTangent(const Eigen::Matrix& recipSemiAxes, + const Eigen::Matrix& w, + const Eigen::Matrix& e, + const Eigen::Matrix& 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 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 v = -e * (1 - t) + w * t; + Eigen::Matrix 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_