From debcca3a1c8effcc364f4d3a4a12d66e807a47c8 Mon Sep 17 00:00:00 2001 From: Andrew Tribick Date: Sat, 13 Nov 2021 22:42:26 +0100 Subject: [PATCH] Further celmath cleanup - Tidy up includes and namespace handling - Replace Ray with Eigen::ParametrizedLine --- src/celengine/body.cpp | 2 +- src/celengine/deepskyobj.cpp | 2 +- src/celengine/deepskyobj.h | 3 +- src/celengine/galaxy.cpp | 10 +-- src/celengine/galaxy.h | 22 ++++-- src/celengine/geometry.h | 10 ++- src/celengine/globular.cpp | 10 +-- src/celengine/globular.h | 22 +++--- src/celengine/hdrfuncrender.cpp | 2 +- src/celengine/lodspheremesh.cpp | 2 +- src/celengine/modelgeometry.cpp | 4 +- src/celengine/modelgeometry.h | 14 ++-- src/celengine/nebula.cpp | 2 +- src/celengine/nebula.h | 2 +- src/celengine/observer.cpp | 4 +- src/celengine/opencluster.cpp | 2 +- src/celengine/opencluster.h | 9 ++- src/celengine/overlayimage.cpp | 2 +- src/celengine/particlesystem.cpp | 2 +- src/celengine/particlesystem.h | 2 +- src/celengine/planetgrid.cpp | 3 +- src/celengine/render.cpp | 32 +++++---- src/celengine/universe.cpp | 33 ++++----- src/celestia/eclipsefinder.cpp | 5 +- src/celestia/qt/qteventfinder.cpp | 2 +- src/celmath/distance.h | 19 ++--- src/celmath/ellipsoid.h | 22 +++--- src/celmath/frustum.cpp | 84 ++++------------------ src/celmath/frustum.h | 9 +-- src/celmath/geomutil.h | 73 +++++++++++-------- src/celmath/intersect.h | 70 +++++++++--------- src/celmath/mathlib.h | 49 +++++++------ src/celmath/perlin.h | 5 +- src/celmath/ray.h | 66 +++++------------ src/celmath/solve.h | 36 ++++++---- src/celmath/sphere.h | 16 ++--- src/tools/atmosphere/scattersim.cpp | 107 +++++++++++++++------------- 37 files changed, 354 insertions(+), 405 deletions(-) diff --git a/src/celengine/body.cpp b/src/celengine/body.cpp index 255771bc5..5268ef596 100644 --- a/src/celengine/body.cpp +++ b/src/celengine/body.cpp @@ -1015,7 +1015,7 @@ void Body::computeLocations() v.normalize(); v *= (float) boundingRadius; - Ray3d ray(v.cast(), -v.cast()); + Eigen::ParametrizedLine ray(v.cast(), -v.cast()); double t = 0.0; if (g->pick(ray, t)) { diff --git a/src/celengine/deepskyobj.cpp b/src/celengine/deepskyobj.cpp index 245146937..fc3f90b36 100644 --- a/src/celengine/deepskyobj.cpp +++ b/src/celengine/deepskyobj.cpp @@ -78,7 +78,7 @@ void DeepSkyObject::setInfoURL(const string& s) } -bool DeepSkyObject::pick(const Ray3d& ray, +bool DeepSkyObject::pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const { diff --git a/src/celengine/deepskyobj.h b/src/celengine/deepskyobj.h index e42f6272d..27659d313 100644 --- a/src/celengine/deepskyobj.h +++ b/src/celengine/deepskyobj.h @@ -13,7 +13,6 @@ #include #include #include -#include #include #ifdef USE_GLCONTEXT #include @@ -83,7 +82,7 @@ class DeepSkyObject : public AstroObject virtual const char* getObjTypeName() const = 0; - virtual bool pick(const celmath::Ray3d& ray, + virtual bool pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const = 0; virtual bool load(AssociativeArray*, const fs::path& resPath); diff --git a/src/celengine/galaxy.cpp b/src/celengine/galaxy.cpp index 06e6d3aa1..9b10bbc35 100644 --- a/src/celengine/galaxy.cpp +++ b/src/celengine/galaxy.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -227,7 +228,7 @@ const char* Galaxy::getObjTypeName() const static const float RADIUS_CORRECTION = 0.025f; static const float MAX_SPIRAL_THICKNESS = 0.06f; -bool Galaxy::pick(const Ray3d& ray, +bool Galaxy::pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const { @@ -246,7 +247,8 @@ bool Galaxy::pick(const Ray3d& ray, getRadius()*(form->scale.z() + RADIUS_CORRECTION)); Matrix3d rotation = getOrientation().cast().toRotationMatrix(); - return testIntersection(Ray3d(ray.origin - getPosition(), ray.direction).transform(rotation), + return testIntersection(transformRay(Eigen::ParametrizedLine(ray.origin() - getPosition(), ray.direction()), + rotation), Ellipsoidd(ellipsoidAxes), distanceToPicker, cosAngleToBoundCenter); @@ -380,7 +382,7 @@ void Galaxy::renderGalaxyPointSprites(const Vector3f& offset, int pow2 = 1; BlobVector* points = form->blobs; - unsigned int nPoints = (unsigned int) (points->size() * clamp(getDetail())); + unsigned int nPoints = (unsigned int) (points->size() * celmath::clamp(getDetail())); // corrections to avoid excessive brightening if viewed e.g. edge-on float brightness_corr = 1.0f; @@ -576,7 +578,7 @@ float Galaxy::getLightGain() void Galaxy::setLightGain(float lg) { - lightGain = clamp(lg); + lightGain = celmath::clamp(lg); } diff --git a/src/celengine/galaxy.h b/src/celengine/galaxy.h index eb07fcad6..f76883755 100644 --- a/src/celengine/galaxy.h +++ b/src/celengine/galaxy.h @@ -8,10 +8,21 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _GALAXY_H_ -#define _GALAXY_H_ +#pragma once -#include +#include +#include + +#include +#include + +#include +#include "deepskyobj.h" + + +class AssociativeArray; +struct Matrices; +class Renderer; struct Blob @@ -39,7 +50,7 @@ class Galaxy : public DeepSkyObject float getDetail() const; void setDetail(float); - bool pick(const celmath::Ray3d& ray, + bool pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const override; bool load(AssociativeArray*, const fs::path&) override; @@ -57,7 +68,7 @@ class Galaxy : public DeepSkyObject static float getLightGain(); static void setLightGain(float); - uint64_t getRenderMask() const override; + std::uint64_t getRenderMask() const override; unsigned int getLabelMask() const override; const char* getObjTypeName() const override; @@ -103,4 +114,3 @@ class Galaxy : public DeepSkyObject static float lightGain; }; -#endif // _GALAXY_H_ diff --git a/src/celengine/geometry.h b/src/celengine/geometry.h index ea8a5b348..575043394 100644 --- a/src/celengine/geometry.h +++ b/src/celengine/geometry.h @@ -8,11 +8,11 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELENGINE_GEOMETRY_H_ -#define _CELENGINE_GEOMETRY_H_ +#pragma once + +#include #include -#include class RenderContext; @@ -31,7 +31,7 @@ public: * and set distance; otherwise return false and leave * distance unmodified. */ - virtual bool pick(const celmath::Ray3d& r, double& distance) const = 0; + virtual bool pick(const Eigen::ParametrizedLine& r, double& distance) const = 0; virtual bool isOpaque() const = 0; @@ -54,5 +54,3 @@ public: { } }; - -#endif // _CELENGINE_GEOMETRY_H_ diff --git a/src/celengine/globular.cpp b/src/celengine/globular.cpp index 4ed01f6d1..ebc666eca 100644 --- a/src/celengine/globular.cpp +++ b/src/celengine/globular.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include "astro.h" @@ -250,7 +251,7 @@ const char* Globular::getObjTypeName() const } constexpr const float RADIUS_CORRECTION = 0.025f; -bool Globular::pick(const Ray3d& ray, +bool Globular::pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const { @@ -266,7 +267,8 @@ bool Globular::pick(const Ray3d& ray, getRadius() * (form->scale.z() + RADIUS_CORRECTION)); Vector3d p = getPosition(); - return testIntersection(Ray3d(ray.origin - p, ray.direction).transform(getOrientation().cast().toRotationMatrix()), + return testIntersection(transformRay(Eigen::ParametrizedLine(ray.origin() - p, ray.direction()), + getOrientation().cast().toRotationMatrix()), Ellipsoidd(ellipsoidAxes), distanceToPicker, cosAngleToBoundCenter); @@ -487,9 +489,9 @@ void Globular::renderGlobularPointSprites( * or when distance from globular center decreases. */ - GLsizei count = (GLsizei) (form->gblobs->size() * clamp(getDetail())); + GLsizei count = (GLsizei) (form->gblobs->size() * celmath::clamp(getDetail())); float t = pow(2, 1 + log2(minimumFeatureSize / brightness) / log2(1/1.25f)); - count = min(count, (GLsizei) clamp(t, 128.0f, (float) max(count, 128))); + count = min(count, (GLsizei) celmath::clamp(t, 128.0f, (float) max(count, 128))); globProg->use(); diff --git a/src/celengine/globular.h b/src/celengine/globular.h index a82fdf885..ef418ed59 100644 --- a/src/celengine/globular.h +++ b/src/celengine/globular.h @@ -12,12 +12,20 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _GLOBULAR_H_ -#define _GLOBULAR_H_ +#pragma once +#include +#include +#include + +#include #include -#include -#include + +#include "deepskyobj.h" +#include "vertexobject.h" + +struct Matrices; +class Renderer; struct GBlob @@ -54,7 +62,7 @@ class Globular : public DeepSkyObject float getBoundingSphereRadius() const override { return tidalRadius; } - bool pick(const celmath::Ray3d& ray, + bool pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const override; bool load(AssociativeArray*, const fs::path&) override; @@ -67,7 +75,7 @@ class Globular : public DeepSkyObject GlobularForm* getForm() const; - uint64_t getRenderMask() const override; + std::uint64_t getRenderMask() const override; unsigned int getLabelMask() const override; const char* getObjTypeName() const override; @@ -93,5 +101,3 @@ class Globular : public DeepSkyObject celgl::VertexObject vo{ GL_ARRAY_BUFFER, 0, GL_STATIC_DRAW }; }; - -#endif // _GLOBULAR_H_ diff --git a/src/celengine/hdrfuncrender.cpp b/src/celengine/hdrfuncrender.cpp index f8d6b29e2..d9742c754 100644 --- a/src/celengine/hdrfuncrender.cpp +++ b/src/celengine/hdrfuncrender.cpp @@ -477,7 +477,7 @@ void Renderer::adjustEclipsedStarExposure(double now) if (lightToBodyDir * bodyToEyeDir > 0.0) { - double dist = distance(posEye, Ray3d(posBody, lightToBodyDir)); + double dist = distance(posEye, Eigen::ParametrizedLine(posBody, lightToBodyDir)); if (dist < body->getRadius()) eyeNotEclipsed = false; } diff --git a/src/celengine/lodspheremesh.cpp b/src/celengine/lodspheremesh.cpp index 307c22056..f95732549 100644 --- a/src/celengine/lodspheremesh.cpp +++ b/src/celengine/lodspheremesh.cpp @@ -177,7 +177,7 @@ void LODSphereMesh::render(unsigned int attributes, lod /= (1 << (-lodBias)); else if (lodBias > 0) lod *= (1 << lodBias); - lod = clamp(lod, 2, maxDivisions); + lod = celmath::clamp(lod, 2, maxDivisions); int step = maxDivisions / lod; int thetaExtent = maxDivisions; diff --git a/src/celengine/modelgeometry.cpp b/src/celengine/modelgeometry.cpp index 8267ccf30..25deac6f6 100644 --- a/src/celengine/modelgeometry.cpp +++ b/src/celengine/modelgeometry.cpp @@ -59,9 +59,9 @@ ModelGeometry::ModelGeometry(unique_ptr&& model) : bool -ModelGeometry::pick(const Ray3d& r, double& distance) const +ModelGeometry::pick(const Eigen::ParametrizedLine& r, double& distance) const { - return m_model->pick(r.origin, r.direction, distance); + return m_model->pick(r.origin(), r.direction(), distance); } diff --git a/src/celengine/modelgeometry.h b/src/celengine/modelgeometry.h index dd7157823..5e98df21a 100644 --- a/src/celengine/modelgeometry.h +++ b/src/celengine/modelgeometry.h @@ -8,13 +8,15 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELENGINE_MODEL_GEOMETRY_H_ -#define _CELENGINE_MODEL_GEOMETRY_H_ +#pragma once + +#include + +#include -#include "geometry.h" #include #include -#include +#include "geometry.h" class CelestiaTextureResource : public cmod::Material::TextureResource @@ -49,7 +51,7 @@ class ModelGeometry : public Geometry * and set distance; otherwise return false and leave * distance unmodified. */ - virtual bool pick(const celmath::Ray3d& r, double& distance) const; + virtual bool pick(const Eigen::ParametrizedLine& r, double& distance) const; //! Render the model in the current OpenGL context virtual void render(RenderContext&, double t = 0.0); @@ -65,5 +67,3 @@ class ModelGeometry : public Geometry bool m_vbInitialized{ false }; std::unique_ptr m_glData; }; - -#endif // !_CELENGINE_MODEL_H_ diff --git a/src/celengine/nebula.cpp b/src/celengine/nebula.cpp index 9742e7a3d..09f56f2e1 100644 --- a/src/celengine/nebula.cpp +++ b/src/celengine/nebula.cpp @@ -58,7 +58,7 @@ const char* Nebula::getObjTypeName() const } -bool Nebula::pick(const Ray3d& ray, +bool Nebula::pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const { diff --git a/src/celengine/nebula.h b/src/celengine/nebula.h index 0e58733ef..61dc517b0 100644 --- a/src/celengine/nebula.h +++ b/src/celengine/nebula.h @@ -24,7 +24,7 @@ class Nebula : public DeepSkyObject void setType(const std::string&) override; std::string getDescription() const override; - bool pick(const celmath::Ray3d& ray, + bool pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const override; bool load(AssociativeArray*, const fs::path&) override; diff --git a/src/celengine/observer.cpp b/src/celengine/observer.cpp index 40e16d548..01942a058 100644 --- a/src/celengine/observer.cpp +++ b/src/celengine/observer.cpp @@ -266,7 +266,7 @@ void Observer::update(double dt, double timeScale) // durations correctly by skipping directly to the destination. float t = 1.0; if (journey.duration > 0) - t = (float) clamp((realTime - journey.startTime) / journey.duration); + t = (float) celmath::clamp((realTime - journey.startTime) / journey.duration); Vector3d jv = journey.to.offsetFromKm(journey.from); UniversalCoord p; @@ -419,7 +419,7 @@ void Observer::update(double dt, double timeScale) if (getVelocity() != targetVelocity) { - double t = clamp((realTime - beginAccelTime) / VELOCITY_CHANGE_TIME); + double t = celmath::clamp((realTime - beginAccelTime) / VELOCITY_CHANGE_TIME); Vector3d v = getVelocity() * (1.0 - t) + targetVelocity * t; // At some threshold, we just set the velocity to zero; otherwise, diff --git a/src/celengine/opencluster.cpp b/src/celengine/opencluster.cpp index d6e7d4bea..d4acccd54 100644 --- a/src/celengine/opencluster.cpp +++ b/src/celengine/opencluster.cpp @@ -46,7 +46,7 @@ const char* OpenCluster::getObjTypeName() const } -bool OpenCluster::pick(const Ray3d& ray, +bool OpenCluster::pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const { diff --git a/src/celengine/opencluster.h b/src/celengine/opencluster.h index cfd9ba2d8..c7c59f984 100644 --- a/src/celengine/opencluster.h +++ b/src/celengine/opencluster.h @@ -7,8 +7,9 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef CELENGINE_OPENCLUSTER_H_ -#define CELENGINE_OPENCLUSTER_H_ +#pragma once + +#include #include #include @@ -25,7 +26,7 @@ class OpenCluster : public DeepSkyObject void setType(const std::string&) override; std::string getDescription() const override; - bool pick(const celmath::Ray3d& ray, + bool pick(const Eigen::ParametrizedLine& ray, double& distanceToPicker, double& cosAngleToBoundCenter) const override; bool load(AssociativeArray*, const fs::path&) override; @@ -53,5 +54,3 @@ class OpenCluster : public DeepSkyObject // TODO: It could be very useful to have a list of stars that are members // of the cluster. }; - -#endif // CELENGINE_OPENCLUSTER_H_ diff --git a/src/celengine/overlayimage.cpp b/src/celengine/overlayimage.cpp index 1780b89e6..5e0b1eee7 100644 --- a/src/celengine/overlayimage.cpp +++ b/src/celengine/overlayimage.cpp @@ -53,7 +53,7 @@ void OverlayImage::render(float curr_time, int width, int height) float alpha = 1.0f; if (curr_time > start + fadeafter) { - alpha = clamp(start + duration - curr_time); + alpha = celmath::clamp(start + duration - curr_time); } celestia::Rect r(left, bottom, xSize, ySize); diff --git a/src/celengine/particlesystem.cpp b/src/celengine/particlesystem.cpp index 4b776d1c5..e1bff2df1 100644 --- a/src/celengine/particlesystem.cpp +++ b/src/celengine/particlesystem.cpp @@ -488,7 +488,7 @@ ParticleSystem::render(RenderContext& rc, double tsec) bool -ParticleSystem::pick(const Ray3d& /* r */, double& /* distance */) const +ParticleSystem::pick(const Eigen::ParametrizedLine& /* r */, double& /* distance */) const { // Pick selection for particle systems not supported (because it's // not typically desirable.) diff --git a/src/celengine/particlesystem.h b/src/celengine/particlesystem.h index 499c8c8af..f876f9acb 100644 --- a/src/celengine/particlesystem.h +++ b/src/celengine/particlesystem.h @@ -98,7 +98,7 @@ class ParticleSystem : public Geometry virtual ~ParticleSystem(); virtual void render(RenderContext& rc, double tsec = 0.0); - virtual bool pick(const celmath::Ray3d& r, double& distance) const; + virtual bool pick(const Eigen::ParametrizedLine& r, double& distance) const; virtual bool isOpaque() const; void addEmitter(ParticleEmitter* emitter); diff --git a/src/celengine/planetgrid.cpp b/src/celengine/planetgrid.cpp index c2864a47a..61c533f0a 100644 --- a/src/celengine/planetgrid.cpp +++ b/src/celengine/planetgrid.cpp @@ -66,7 +66,8 @@ static void longLatLabel(const string& labelText, // Draw the label only if it isn't obscured by the body ellipsoid double t = 0.0; - if (testIntersection(Ray3d(viewRayOrigin, pos - viewRayOrigin), Ellipsoidd(semiAxes.cast()), t) && t >= 1.0) + if (testIntersection(Eigen::ParametrizedLine(viewRayOrigin, pos - viewRayOrigin), + Ellipsoidd(semiAxes.cast()), t) && t >= 1.0) { // Compute the position of the label Vector3d labelPos = bodyCenter + diff --git a/src/celengine/render.cpp b/src/celengine/render.cpp index b02e36ef3..f0ee3a95e 100644 --- a/src/celengine/render.cpp +++ b/src/celengine/render.cpp @@ -33,6 +33,8 @@ std::ofstream hdrlog; #define EXPOSURE_HALFLIFE 0.4f #endif +#include + #include #include "render.h" @@ -2122,7 +2124,7 @@ void Renderer::renderEllipsoidAtmosphere(const Atmosphere& atmosphere, // Gradually fade in the atmosphere if it's thickness on screen is just // over one pixel. - float fade = clamp(pixSize - 2); + float fade = celmath::clamp(pixSize - 2); Matrix3f rot = orientation.toRotationMatrix(); Matrix3f irot = orientation.conjugate().toRotationMatrix(); @@ -2508,7 +2510,7 @@ Renderer::locationsToAnnotations(const Body& body, // the label is obscured by the planet. An exact calculation // for irregular objects would be too expensive, and the // ellipsoid approximation works reasonably well for them. - Ray3d testRay(viewRayOrigin, pcLabelPos - viewRayOrigin); + Eigen::ParametrizedLine testRay(viewRayOrigin, pcLabelPos - viewRayOrigin); bool hit = testIntersection(testRay, bodyEllipsoid, t); if (!hit || t >= 1.0) @@ -3035,7 +3037,7 @@ void Renderer::renderObject(const Vector3f& pos, { thicknessInPixels = atmosphere->height / ((distance - radius) * pixelSize); - fade = clamp(thicknessInPixels - 2); + fade = celmath::clamp(thicknessInPixels - 2); } else { @@ -3220,7 +3222,7 @@ bool Renderer::testEclipse(const Body& receiver, Vector3d receiverToCasterDir = posReceiver - posCaster; double dist = distance(posReceiver, - Ray3d(posCaster, lightToCasterDir)); + Eigen::ParametrizedLine(posCaster, lightToCasterDir)); if (dist < R && lightToCasterDir.dot(receiverToCasterDir) > 0.0) { Vector3d sunDir = lightToCasterDir.normalized(); @@ -3933,8 +3935,8 @@ void Renderer::renderAsterisms(const Universe& universe, float dist, const Matri float opacity = 1.0f; if (dist > MaxAsterismLinesConstDist) { - opacity = clamp((MaxAsterismLinesConstDist - dist) / - (MaxAsterismLinesDist - MaxAsterismLinesConstDist) + 1); + opacity = celmath::clamp((MaxAsterismLinesConstDist - dist) / + (MaxAsterismLinesDist - MaxAsterismLinesConstDist) + 1); } enableSmoothLines(); @@ -3964,8 +3966,8 @@ void Renderer::renderBoundaries(const Universe& universe, float dist, const Matr float opacity = 1.0f; if (dist > MaxAsterismLabelsConstDist) { - opacity = clamp((MaxAsterismLabelsConstDist - dist) / - (MaxAsterismLabelsDist - MaxAsterismLabelsConstDist) + 1); + opacity = celmath::clamp((MaxAsterismLabelsConstDist - dist) / + (MaxAsterismLabelsDist - MaxAsterismLabelsConstDist) + 1); } enableSmoothLines(); @@ -4501,7 +4503,7 @@ void Renderer::buildLabelLists(const Frustum& viewFrustum, lastPrimary = primary; } - Ray3d testRay(Vector3d::Zero(), pos.cast()); + Eigen::ParametrizedLine testRay(Vector3d::Zero(), pos.cast()); // Test the viewer-to-labeled object ray against // the primary sphere (TODO: handle ellipsoids) @@ -4905,8 +4907,8 @@ void Renderer::labelConstellations(const AsterismList& asterisms, float dist = observerPos.norm(); if (dist > MaxAsterismLabelsConstDist) { - opacity = clamp((MaxAsterismLabelsConstDist - dist) / - (MaxAsterismLabelsDist - MaxAsterismLabelsConstDist) + 1); + opacity = celmath::clamp((MaxAsterismLabelsConstDist - dist) / + (MaxAsterismLabelsDist - MaxAsterismLabelsConstDist) + 1); } // Use the default label color unless the constellation has an @@ -5158,7 +5160,7 @@ Renderer::renderAnnotations(vector::iterator startIter, { // Compute normalized device z float z = fisheye ? (1.0f - (iter->position.z() - nearDist) / d0 * 2.0f) : (d1 + d2 / -iter->position.z()); - float ndc_z = clamp(z, -1.0f, 1.0f); + float ndc_z = celmath::clamp(z, -1.0f, 1.0f); // Offsets to left align label int labelHOffset = 0; @@ -5361,7 +5363,7 @@ void Renderer::updateBodyVisibilityMask() void Renderer::setSolarSystemMaxDistance(float t) { - SolarSystemMaxDistance = clamp(t, 1.0f, 10.0f); + SolarSystemMaxDistance = celmath::clamp(t, 1.0f, 10.0f); } void Renderer::getViewport(int* x, int* y, int* w, int* h) const @@ -5802,7 +5804,7 @@ Renderer::setShadowMapSize(unsigned size) return; GLint t = 0; glGetIntegerv(GL_MAX_TEXTURE_SIZE, &t); - m_shadowMapSize = clamp(size, 0u, static_cast(t)); + m_shadowMapSize = celmath::clamp(size, 0u, static_cast(t)); if (m_shadowFBO != nullptr && m_shadowMapSize == m_shadowFBO->width()) return; if (m_shadowMapSize == 0) @@ -6047,7 +6049,7 @@ Renderer::adjustMagnitudeInsideAtmosphere(float &faintestMag, maxBodyMag = maxBodyMagPrev; saturationMag = maxBodyMag; #endif - float illumination = clamp(sunDir.dot(normal) + 0.2f); + float illumination = celmath::clamp(sunDir.dot(normal) + 0.2f); float lightness = illumination * density; faintestMag = faintestMag - 15.0f * lightness; diff --git a/src/celengine/universe.cpp b/src/celengine/universe.cpp index a9ba63c12..d12ae848f 100644 --- a/src/celengine/universe.cpp +++ b/src/celengine/universe.cpp @@ -20,6 +20,7 @@ #include "frametree.h" #include #include +#include #include #include @@ -290,7 +291,7 @@ struct PlanetPickInfo double closestDistance; double closestApproxDistance; Body* closestBody; - Ray3d pickRay; + Eigen::ParametrizedLine pickRay; double jd; float atanTolerance; }; @@ -305,7 +306,7 @@ static bool ApproxPlanetPickTraversal(Body* body, void* info) return true; Vector3d bpos = body->getAstrocentricPosition(pickInfo->jd); - Vector3d bodyDir = bpos - pickInfo->pickRay.origin; + Vector3d bodyDir = bpos - pickInfo->pickRay.origin(); double distance = bodyDir.norm(); // Check the apparent radius of the orbit against our tolerance factor. @@ -319,7 +320,7 @@ static bool ApproxPlanetPickTraversal(Body* body, void* info) } bodyDir.normalize(); - Vector3d bodyMiss = bodyDir - pickInfo->pickRay.direction; + Vector3d bodyMiss = bodyDir - pickInfo->pickRay.direction(); double sinAngle2 = bodyMiss.norm() / 2.0; if (sinAngle2 <= pickInfo->sinAngle2Closest) @@ -358,8 +359,8 @@ static bool ExactPlanetPickTraversal(Body* body, void* info) // Transform rotate the pick ray into object coordinates Matrix3d m = body->getEclipticToEquatorial(pickInfo->jd).toRotationMatrix(); - Ray3d r(pickInfo->pickRay.origin - bpos, pickInfo->pickRay.direction); - r = r.transform(m); + Eigen::ParametrizedLine r(pickInfo->pickRay.origin() - bpos, pickInfo->pickRay.direction()); + r = transformRay(r, m); if (!testIntersection(r, Ellipsoidd(ellipsoidAxes), distance)) distance = -1.0; } @@ -369,8 +370,8 @@ static bool ExactPlanetPickTraversal(Body* body, void* info) // Transform rotate the pick ray into object coordinates Quaterniond qd = body->getGeometryOrientation().cast(); Matrix3d m = (qd * body->getEclipticToBodyFixed(pickInfo->jd)).toRotationMatrix(); - Ray3d r(pickInfo->pickRay.origin - bpos, pickInfo->pickRay.direction); - r = r.transform(m); + Eigen::ParametrizedLine r(pickInfo->pickRay.origin() - bpos, pickInfo->pickRay.direction()); + r = transformRay(r, m); Geometry* geometry = GetGeometryManager()->find(body->getGeometry()); float scaleFactor = body->getGeometryScale(); @@ -381,8 +382,8 @@ static bool ExactPlanetPickTraversal(Body* body, void* info) // factor. Thus, the ray needs to be multiplied by the inverse of // the mesh scale factor. double is = 1.0 / scaleFactor; - r.origin *= is; - r.direction *= is; + r.origin() *= is; + r.direction() *= is; if (geometry != nullptr) { @@ -393,9 +394,9 @@ static bool ExactPlanetPickTraversal(Body* body, void* info) // Make also sure that the pickRay does not intersect the body in the // opposite hemisphere! Hence, need again the "bodyMiss" angle - Vector3d bodyDir = bpos - pickInfo->pickRay.origin; + Vector3d bodyDir = bpos - pickInfo->pickRay.origin(); bodyDir.normalize(); - Vector3d bodyMiss = bodyDir - pickInfo->pickRay.direction; + Vector3d bodyMiss = bodyDir - pickInfo->pickRay.direction(); double sinAngle2 = bodyMiss.norm() / 2.0; @@ -459,7 +460,7 @@ Selection Universe::pickPlanet(SolarSystem& solarSystem, // Transform the pick ray origin into astrocentric coordinates Vector3d astrocentricOrigin = origin.offsetFromKm(star->getPosition(when)); - pickInfo.pickRay = Ray3d(astrocentricOrigin, direction.cast()); + pickInfo.pickRay = Eigen::ParametrizedLine(astrocentricOrigin, direction.cast()); pickInfo.sinAngle2Closest = 1.0; pickInfo.closestDistance = 1.0e50; pickInfo.closestApproxDistance = 1.0e50; @@ -555,7 +556,7 @@ void StarPicker::process(const Star& star, float /*unused*/, float /*unused*/) // no intersection, then just use normal calculation. We actually test // intersection with a larger sphere to make sure we don't miss a star // right on the edge of the sphere. - if (testIntersection(Ray3f(Vector3f::Zero(), pickRay), + if (testIntersection(Eigen::ParametrizedLine(Vector3f::Zero(), pickRay), Spheref(relativeStarPos, orbitalRadius * 2.0f), distance)) { @@ -627,7 +628,7 @@ void CloseStarPicker::process(const Star& star, float distance = 0.0f; - if (testIntersection(Ray3f(Vector3f::Zero(), pickDir), + if (testIntersection(Eigen::ParametrizedLine(Vector3f::Zero(), pickDir), Spheref(starDir, star.getRadius()), distance)) { if (distance > 0.0f) @@ -740,7 +741,7 @@ void DSOPicker::process(DeepSkyObject* const & dso, double /*unused*/, float /*u Vector3d dsoDir = relativeDSOPos; double distance2 = 0.0; - if (testIntersection(Ray3d(Vector3d::Zero(), pickDir), + if (testIntersection(Eigen::ParametrizedLine(Vector3d::Zero(), pickDir), Sphered(relativeDSOPos, (double) dso->getRadius()), distance2)) { Vector3d dsoPos = dso->getPosition(); @@ -806,7 +807,7 @@ void CloseDSOPicker::process(DeepSkyObject* const & dso, double distanceToPicker = 0.0; double cosAngleToBoundCenter = 0.0; - if (dso->pick(Ray3d(pickOrigin, pickDir), distanceToPicker, cosAngleToBoundCenter)) + if (dso->pick(Eigen::ParametrizedLine(pickOrigin, pickDir), distanceToPicker, cosAngleToBoundCenter)) { // Don't select the object the observer is currently in: if ((pickOrigin - dso->getPosition()).norm() > dso->getRadius() && diff --git a/src/celestia/eclipsefinder.cpp b/src/celestia/eclipsefinder.cpp index da0237181..e372cd810 100644 --- a/src/celestia/eclipsefinder.cpp +++ b/src/celestia/eclipsefinder.cpp @@ -12,6 +12,9 @@ #include #include + +#include + #include "eclipsefinder.h" #include "celmath/ray.h" #include "celmath/distance.h" @@ -85,7 +88,7 @@ bool testEclipse(const Body& receiver, const Body& caster, double now) // If the distance is less than the sum of the caster's and receiver's // radii, then we have an eclipse. float R = receiver.getRadius() + shadowRadius; - double dist = distance(posReceiver, Ray3d(posCaster, posCaster)); + double dist = distance(posReceiver, Eigen::ParametrizedLine(posCaster, posCaster)); if (dist < R) { // Ignore "eclipses" where the caster and receiver have diff --git a/src/celestia/qt/qteventfinder.cpp b/src/celestia/qt/qteventfinder.cpp index f0ed11212..671021b41 100644 --- a/src/celestia/qt/qteventfinder.cpp +++ b/src/celestia/qt/qteventfinder.cpp @@ -430,7 +430,7 @@ static Vector3d findMaxEclipsePoint(const Vector3d& toCasterDir, double eclipsedBodyRadius) { double distance = 0.0; - bool intersect = testIntersection(Ray3d(Vector3d::Zero(), toCasterDir), + bool intersect = testIntersection(Eigen::ParametrizedLine(Vector3d::Zero(), toCasterDir), Sphered(toReceiver, eclipsedBodyRadius), distance); diff --git a/src/celmath/distance.h b/src/celmath/distance.h index 439a00633..d112fb630 100644 --- a/src/celmath/distance.h +++ b/src/celmath/distance.h @@ -9,24 +9,19 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_DISTANCE_H_ -#define _CELMATH_DISTANCE_H_ +#pragma once -#include "mathlib.h" -#include "ray.h" -#include "sphere.h" -#include "ellipsoid.h" #include +#include namespace celmath { -template T distance(const Eigen::Matrix& p, const Ray3& r) +template T distance(const Eigen::Matrix& p, const Eigen::ParametrizedLine& r) { - T t = ((p - r.origin).dot(r.direction)) / r.direction.squaredNorm(); - if (t <= 0) - return (p - r.origin).norm(); + T t = ((p - r.origin()).dot(r.direction())) / r.direction().squaredNorm(); + if (t <= static_cast(0)) + return (p - r.origin()).norm(); else - return (p - r.point(t)).norm(); + return (p - r.pointAt(t)).norm(); } } -#endif // _CELMATH_DISTANCE_H_ diff --git a/src/celmath/ellipsoid.h b/src/celmath/ellipsoid.h index 50dd7466e..047d131aa 100644 --- a/src/celmath/ellipsoid.h +++ b/src/celmath/ellipsoid.h @@ -7,15 +7,14 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_ELLIPSOID_H_ -#define _CELMATH_ELLIPSOID_H_ +#pragma once #include namespace celmath { - -template class Ellipsoid +template +class Ellipsoid { public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW @@ -24,8 +23,8 @@ template class Ellipsoid * at the origin. */ Ellipsoid() : - center(0, 0, 0), - axes(1, 1, 1) + center(Eigen::Matrix::Zero()), + axes(Eigen::Matrix::Ones()) { } @@ -33,7 +32,7 @@ template class Ellipsoid * at the origin. */ Ellipsoid(const Eigen::Matrix& _axes) : - center(0, 0, 0), + center(Eigen::Matrix::Zero()), axes(_axes) { } @@ -53,7 +52,7 @@ template class Ellipsoid { Eigen::Matrix v = p - center; v = v.cwise() / axes; - return v * v <= (T) 1.0; + return v * v <= static_cast(1); } @@ -62,9 +61,6 @@ template class Ellipsoid Eigen::Matrix axes; }; -typedef Ellipsoid Ellipsoidf; -typedef Ellipsoid Ellipsoidd; - +using Ellipsoidf = Ellipsoid; +using Ellipsoidd = Ellipsoid; } // namespace celmath - -#endif // _CELMATH_ELLIPSOID_H_ diff --git a/src/celmath/frustum.cpp b/src/celmath/frustum.cpp index 50b6ee55f..546cf3347 100644 --- a/src/celmath/frustum.cpp +++ b/src/celmath/frustum.cpp @@ -7,11 +7,12 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#include "frustum.h" -#include #include -using namespace Eigen; +#include + +#include "frustum.h" + namespace celmath { @@ -35,18 +36,18 @@ void Frustum::init(float fov, float aspectRatio, float n, float f) float h = std::tan(fov / 2.0f); float w = h * aspectRatio; - Vector3f normals[4]; - normals[Bottom] = Vector3f( 0.0f, 1.0f, -h); - normals[Top] = Vector3f( 0.0f, -1.0f, -h); - normals[Left] = Vector3f( 1.0f, 0.0f, -w); - normals[Right] = Vector3f(-1.0f, 0.0f, -w); + Eigen::Vector3f normals[4]; + normals[Bottom] = Eigen::Vector3f( 0.0f, 1.0f, -h); + normals[Top] = Eigen::Vector3f( 0.0f, -1.0f, -h); + normals[Left] = Eigen::Vector3f( 1.0f, 0.0f, -w); + normals[Right] = Eigen::Vector3f(-1.0f, 0.0f, -w); for (unsigned int i = 0; i < 4; i++) { - planes[i] = Hyperplane(normals[i].normalized(), 0.0f); + planes[i] = Eigen::Hyperplane(normals[i].normalized(), 0.0f); } - planes[Near] = Hyperplane(Vector3f(0.0f, 0.0f, -1.0f), -n); - planes[Far] = Hyperplane(Vector3f(0.0f, 0.0f, 1.0f), f); + planes[Near] = Eigen::Hyperplane(Eigen::Vector3f(0.0f, 0.0f, -1.0f), -n); + planes[Far] = Eigen::Hyperplane(Eigen::Vector3f(0.0f, 0.0f, 1.0f), f); } @@ -108,54 +109,8 @@ Frustum::testSphere(const Eigen::Vector3d& center, double radius) const } -/* -Frustum::Aspect Frustum::testCapsule(const Capsulef& capsule) const -{ - int nPlanes = infinite ? 5 : 6; - int intersections = 0; - float r2 = capsule.radius * capsule.radius; - - // TODO: Unnecessary after Eigen conversion of Capsule class - Vector3f capsuleOrigin(capsule.origin.x, capsule.origin.y, capsule.origin.z); - Vector3f capsuleAxis(capsule.axis.x, capsule.axis.y, capsule.axis.z); - - for (int i = 0; i < nPlanes; i++) - { - float signedDist0 = planes[i].signedDistance(capsuleOrigin); - float signedDist1 = planes[i].signedDistance(capsuleOrigin + capsuleAxis); - //float signedDist1 = signedDist0 + planes[i].normal * capsule.axis; - if (signedDist0 * signedDist1 > r2) - { - // Endpoints of capsule are on same side of plane; test closest endpoint to see if it - // lies closer to the plane than radius - if (abs(signedDist0) <= abs(signedDist1)) - { - if (signedDist0 < -capsule.radius) - return Outside; - if (signedDist0 < capsule.radius) - intersections |= (1 << i); - } - else - { - if (signedDist1 < -capsule.radius) - return Outside; - if (signedDist1 < capsule.radius) - intersections |= (1 << i); - } - } - else - { - // Capsule endpoints are on different sides of the plane, so we have an intersection - intersections |= (1 << i); - } - } - - return (intersections == 0) ? Inside : Intersect; -} -*/ - void -Frustum::transform(const Matrix3f& m) +Frustum::transform(const Eigen::Matrix3f& m) { unsigned int nPlanes = infinite ? 5 : 6; @@ -167,24 +122,15 @@ Frustum::transform(const Matrix3f& m) void -Frustum::transform(const Matrix4f& m) +Frustum::transform(const Eigen::Matrix4f& m) { unsigned int nPlanes = infinite ? 5 : 6; - Matrix4f invTranspose = m.inverse().transpose(); + Eigen::Matrix4f invTranspose = m.inverse().transpose(); for (unsigned int i = 0; i < nPlanes; i++) { planes[i].coeffs() = invTranspose * planes[i].coeffs(); planes[i].normalize(); - - //float s = 1.0f / planes[i].normal().norm(); - //planes[i].normal() *= s; - //planes[i].offset() *= s; - - //planes[i] = planes[i] * invTranspose; - //float s = 1.0f / planes[i].normal().norm(); - //planes[i].normal = planes[i].normal * s; - //planes[i].d *= s; } } diff --git a/src/celmath/frustum.h b/src/celmath/frustum.h index d09edec7c..92fc130e2 100644 --- a/src/celmath/frustum.h +++ b/src/celmath/frustum.h @@ -7,13 +7,11 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_FRUSTUM_H_ -#define _CELMATH_FRUSTUM_H_ +#pragma once #include #include -//#include namespace celmath { @@ -23,7 +21,7 @@ class Frustum public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW - typedef Eigen::Hyperplane PlaneType; + using PlaneType = Eigen::Hyperplane; Frustum(float fov, float aspectRatio, float nearDist); Frustum(float fov, float aspectRatio, float nearDist, float farDist); @@ -54,7 +52,6 @@ class Frustum Aspect test(const Eigen::Vector3f& point) const; Aspect testSphere(const Eigen::Vector3f& center, float radius) const; Aspect testSphere(const Eigen::Vector3d& center, double radius) const; -// Aspect testCapsule(const Capsulef&) const; private: void init(float, float, float, float); @@ -64,5 +61,3 @@ class Frustum }; } // namespace celmath - -#endif // _CELMATH_FRUSTUM_H_ diff --git a/src/celmath/geomutil.h b/src/celmath/geomutil.h index b00443327..9b9f0ca03 100644 --- a/src/celmath/geomutil.h +++ b/src/celmath/geomutil.h @@ -8,12 +8,12 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_GEOMUTIL_H_ -#define _CELMATH_GEOMUTIL_H_ +#pragma once + +#include #include #include -#include namespace celmath { @@ -21,24 +21,30 @@ namespace celmath template Eigen::Quaternion XRotation(T radians) { - T halfAngle = radians * T(0.5); - return Eigen::Quaternion(std::cos(halfAngle), std::sin(halfAngle), 0, 0); + using std::cos; + using std::sin; + T halfAngle = radians * static_cast(0.5); + return Eigen::Quaternion(cos(halfAngle), sin(halfAngle), 0, 0); } template Eigen::Quaternion YRotation(T radians) { - T halfAngle = radians * T(0.5); - return Eigen::Quaternion(std::cos(halfAngle), 0, std::sin(halfAngle), 0); + using std::cos; + using std::sin; + T halfAngle = radians * static_cast(0.5); + return Eigen::Quaternion(cos(halfAngle), 0, sin(halfAngle), 0); } template Eigen::Quaternion ZRotation(T radians) { - T halfAngle = radians * T(0.5); - return Eigen::Quaternion(std::cos(halfAngle), 0, 0, std::sin(halfAngle)); + using std::cos; + using std::sin; + T halfAngle = radians * static_cast(0.5); + return Eigen::Quaternion(cos(halfAngle), 0, 0, sin(halfAngle)); } @@ -70,14 +76,14 @@ ProjectPerspective(const Eigen::Matrix& from, const int viewport[4], Eigen::Matrix& to) { - Eigen::Matrix in(from.x(), from.y(), from.z(), T(1.0)); + Eigen::Matrix in(from.x(), from.y(), from.z(), static_cast(1)); Eigen::Matrix out = modelViewProjectionMatrix * in; - if (out.w() == T(0.0)) + if (out.w() == static_cast(0)) return false; out = out.array() / out.w(); // Map x, y and z to range 0-1 - out = T(0.5) + out.array() * T(0.5); + out = static_cast(0.5) + out.array() * static_cast(0.5); // Map x,y to viewport out.x() = viewport[0] + out.x() * viewport[2]; out.y() = viewport[1] + out.y() * viewport[3]; @@ -104,23 +110,28 @@ ProjectFisheye(const Eigen::Matrix& from, const int viewport[4], Eigen::Matrix& to) { - Eigen::Matrix inPos = modelViewMatrix * Eigen::Matrix(from.x(), from.y(), from.z(), T(1.0)); + using std::atan2; + using std::hypot; + Eigen::Matrix inPos = modelViewMatrix * Eigen::Matrix(from.x(), + from.y(), + from.z(), + static_cast(1)); T l = hypot(inPos.x(), inPos.y()); - if (l != T(0.0)) + if (l != static_cast(0)) { T phi = atan2(l, -inPos.z()); - T ratio = phi / T(M_PI_2) / l; + T ratio = phi / static_cast(M_PI_2) / l; inPos.x() *= ratio; inPos.y() *= ratio; } Eigen::Matrix out = projMatrix * inPos; - if (out.w() == T(0.0)) + if (out.w() == static_cast(0)) return false; out = out.array() / out.w(); // Map x, y and z to range 0-1 - out = T(0.5) + out.array() * T(0.5); + out = static_cast(0.5) + out.array() * static_cast(0.5); // Map x,y to viewport out.x() = viewport[0] + out.x() * viewport[2]; out.y() = viewport[1] + out.y() * viewport[3]; @@ -135,16 +146,19 @@ ProjectFisheye(const Eigen::Matrix& from, template Eigen::Matrix Perspective(T fovy, T aspect, T nearZ, T farZ) { - if (aspect == T(0.0)) + using std::cos; + using std::sin; + + if (aspect == static_cast(0)) return Eigen::Matrix::Identity(); T deltaZ = farZ - nearZ; - if (deltaZ == T(0.0)) + if (deltaZ == static_cast(0)) return Eigen::Matrix::Identity(); - T angle = degToRad(fovy / 2); + T angle = degToRad(fovy / static_cast(2)); T sine = sin(angle); - if (sine == T(0.0)) + if (sine == static_cast(0)) return Eigen::Matrix::Identity(); T ctg = cos(angle) / sine; @@ -152,9 +166,9 @@ Perspective(T fovy, T aspect, T nearZ, T farZ) m(0, 0) = ctg / aspect; m(1, 1) = ctg; m(2, 2) = -(farZ + nearZ) / deltaZ; - m(2, 3) = T(-2.0) * nearZ * farZ / deltaZ; - m(3, 2) = T(-1.0); - m(3, 3) = T(0.0); + m(2, 3) = static_cast(-2) * nearZ * farZ / deltaZ; + m(3, 2) = static_cast(-1); + m(3, 3) = static_cast(0); return m; } @@ -167,18 +181,17 @@ Ortho(T left, T right, T bottom, T top, T nearZ, T farZ) T tb = top - bottom; T fn = farZ - nearZ; Eigen::Matrix m; - m << 2/rl, 0, 0, - (right + left) / rl, - 0, 2/tb, 0, - (top + bottom) / tb, - 0, 0, -2/fn, - (farZ + nearZ) / fn, - 0, 0, 0, 1; + m << static_cast(2)/rl, static_cast(0), static_cast(0), -(right + left)/rl, + static_cast(0), static_cast(2)/tb, static_cast(0), -(top + bottom)/tb, + static_cast(0), static_cast(0), static_cast(-2)/fn, -(farZ + nearZ)/fn, + static_cast(0), static_cast(0), static_cast(0), static_cast(1); return m; } template Eigen::Matrix Ortho2D(T left, T right, T bottom, T top) { - return Ortho(left, right, bottom, top, T(-1), T(1)); + return Ortho(left, right, bottom, top, static_cast(-1), static_cast(1)); } } // namespace celmath -#endif // _CELMATH_GEOMUTIL_H_ diff --git a/src/celmath/intersect.h b/src/celmath/intersect.h index 93135278a..b3ecc3b6c 100644 --- a/src/celmath/intersect.h +++ b/src/celmath/intersect.h @@ -9,44 +9,48 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_INTERSECT_H_ -#define _CELMATH_INTERSECT_H_ +#pragma once -#include "ray.h" +#include + +#include +#include + +#include "mathlib.h" #include "sphere.h" #include "ellipsoid.h" -#include "mathlib.h" -#include namespace celmath { -template bool testIntersection(const Ray3& ray, +template bool testIntersection(const Eigen::ParametrizedLine& ray, const Sphere& sphere, T& distance) { - Eigen::Matrix diff = ray.origin - sphere.center; - T s = (T) 1.0 / square(sphere.radius); - T a = ray.direction.squaredNorm() * s; - T b = ray.direction.dot(diff) * s; + using std::sqrt; + + Eigen::Matrix diff = ray.origin() - sphere.center; + T s = static_cast(1) / square(sphere.radius); + T a = ray.direction().squaredNorm() * s; + T b = ray.direction().dot(diff) * s; T c = diff.squaredNorm() * s - (T) 1.0; T disc = b * b - a * c; - if (disc < 0.0) + if (disc < static_cast(0)) return false; - disc = (T) sqrt(disc); + disc = sqrt(disc); T sol0 = (-b + disc) / a; T sol1 = (-b - disc) / a; - if (sol0 > 0) + if (sol0 > static_cast(0)) { - if (sol0 < sol1 || sol1 < 0) + if (sol0 < sol1 || sol1 < static_cast(0)) distance = sol0; else distance = sol1; return true; } - else if (sol1 > 0) + else if (sol1 > static_cast(0)) { distance = sol1; return true; @@ -58,49 +62,52 @@ template bool testIntersection(const Ray3& ray, } -template bool testIntersection(const Ray3& ray, +template bool testIntersection(const Eigen::ParametrizedLine& ray, const Sphere& sphere, T& distanceToTester, T& cosAngleToCenter) { if (testIntersection(ray, sphere, distanceToTester)) { - cosAngleToCenter = (sphere.center - ray.origin).dot(ray.direction) / (sphere.center - ray.origin).norm(); + cosAngleToCenter = (sphere.center - ray.origin()).dot(ray.direction()) + / (sphere.center - ray.origin()).norm(); return true; } return false; } -template bool testIntersection(const Ray3& ray, +template bool testIntersection(const Eigen::ParametrizedLine& ray, const Ellipsoid& e, T& distance) { - Eigen::Matrix diff = ray.origin - e.center; + using std::sqrt; + + Eigen::Matrix diff = ray.origin() - e.center; Eigen::Matrix s = e.axes.cwiseInverse().array().square(); - Eigen::Matrix sdir = ray.direction.cwiseProduct(s); + Eigen::Matrix sdir = ray.direction().cwiseProduct(s); Eigen::Matrix sdiff = diff.cwiseProduct(s); - T a = ray.direction.dot(sdir); - T b = ray.direction.dot(sdiff); - T c = diff.dot(sdiff) - (T) 1.0; + T a = ray.direction().dot(sdir); + T b = ray.direction().dot(sdiff); + T c = diff.dot(sdiff) - static_cast(1); T disc = b * b - a * c; - if (disc < 0.0) + if (disc < static_cast(0)) return false; - disc = (T) sqrt(disc); + disc = sqrt(disc); T sol0 = (-b + disc) / a; T sol1 = (-b - disc) / a; - if (sol0 > 0) + if (sol0 > static_cast(0)) { - if (sol0 < sol1 || sol1 < 0) + if (sol0 < sol1 || sol1 < static_cast(0)) distance = sol0; else distance = sol1; return true; } - else if (sol1 > 0) + else if (sol1 > static_cast(0)) { distance = sol1; return true; @@ -112,18 +119,17 @@ template bool testIntersection(const Ray3& ray, } -template bool testIntersection(const Ray3& ray, +template bool testIntersection(const Eigen::ParametrizedLine& ray, const Ellipsoid& ellipsoid, T& distanceToTester, T& cosAngleToCenter) { if (testIntersection(ray, ellipsoid, distanceToTester)) { - cosAngleToCenter = (ellipsoid.center - ray.origin).dot(ray.direction) / (ellipsoid.center - ray.origin).norm(); + cosAngleToCenter = (ellipsoid.center - ray.origin()).dot(ray.direction()) + / (ellipsoid.center - ray.origin()).norm(); return true; } return false; } } // namespace celmath - -#endif // _CELMATH_INTERSECT_H_ diff --git a/src/celmath/mathlib.h b/src/celmath/mathlib.h index 5dec8eabd..29a0e99eb 100644 --- a/src/celmath/mathlib.h +++ b/src/celmath/mathlib.h @@ -7,8 +7,7 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_MATHLIB_H_ -#define _CELMATH_MATHLIB_H_ +#pragma once #include #include @@ -20,6 +19,8 @@ namespace celmath { template inline void sincos(T angle, T& s, T& c) { + using std::cos; + using std::sin; s = sin(angle); c = cos(angle); } @@ -46,7 +47,9 @@ template constexpr T lerp(T t, T a, T b) // return t clamped to [0, 1] template constexpr T clamp(T t) { - return (t < 0) ? 0 : ((t > 1) ? 1 : t); + return t < static_cast(0) ? static_cast(0) + : t > static_cast(1) ? 1 + : t; } #if __cplusplus < 201703L @@ -55,16 +58,18 @@ template constexpr T clamp(T t, T low, T high) { return (t < low) ? low : ((t > high) ? high : t); } +#else +using std::clamp; #endif template inline constexpr T degToRad(T d) { - return d / 180 * static_cast(PI); + return d / static_cast(180) * static_cast(PI); } template inline constexpr T radToDeg(T r) { - return r * 180 / static_cast(PI); + return r * static_cast(180) / static_cast(PI); } template inline constexpr T square(T x) @@ -79,13 +84,18 @@ template inline constexpr T cube(T x) template inline constexpr T sign(T x) { - return (x < 0) ? -1 : ((x > 0) ? 1 : 0); + return x < static_cast(0) ? static_cast(-1) + : x > static_cast(0) ? static_cast(1) + : static_cast(0); } // This function is like fmod except that it always returns // a positive value in the range [ 0, y ) template T pfmod(T x, T y) { + using std::abs; + using std::floor; + T quotient = floor(abs(x / y)); if (x < 0.0) return x + (quotient + 1) * y; @@ -100,7 +110,7 @@ template inline constexpr T circleArea(T r) template inline constexpr T sphereArea(T r) { - return 4 * static_cast(PI) * r * r; + return static_cast(4) * static_cast(PI) * r * r; } template static Eigen::Matrix @@ -110,6 +120,8 @@ ellipsoidTangent(const Eigen::Matrix& recipSemiAxes, const Eigen::Matrix& e_, T ee) { + using std::sqrt; + // 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 @@ -128,9 +140,9 @@ ellipsoidTangent(const Eigen::Matrix& recipSemiAxes, // 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 a = static_cast(4) * (square(ew) - ee * ww + ee + static_cast(2) * ew + ww); + T b = static_cast(-8) * (ee + ew); + T c = static_cast(4) * ee; T t = 0; T discriminant = b * b - 4 * a * c; @@ -138,32 +150,23 @@ ellipsoidTangent(const Eigen::Matrix& recipSemiAxes, if (discriminant < 0) discriminant = -discriminant; // Bad! - t = (-b + (T) sqrt(discriminant)) / (2 * a); + t = (-b + sqrt(discriminant)) / (static_cast(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 = -e * (static_cast(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); + T b1 = static_cast(2) * v_.dot(e_); + T t1 = -b1 / (static_cast(2) * a1); return e + v * t1; } } // namespace celmath -#if __cplusplus < 201703L -namespace std -{ - using celmath::clamp; -} -#endif - constexpr long double operator"" _deg (long double deg) { return celmath::degToRad(deg); } - -#endif // _MATHLIB_H_ diff --git a/src/celmath/perlin.h b/src/celmath/perlin.h index fe172ab5c..aea464fda 100644 --- a/src/celmath/perlin.h +++ b/src/celmath/perlin.h @@ -1,7 +1,6 @@ // perlin.h -#ifndef _PERLIN_H_ -#define _PERLIN_H_ +#pragma once #include @@ -18,5 +17,3 @@ extern float turbulence(const Eigen::Vector3f& p, float freq); extern float fractalsum(const float v[], float freq); extern float fractalsum(const Eigen::Vector2f& p, float freq); extern float fractalsum(const Eigen::Vector3f& p, float freq); - -#endif // _PERLIN_H_ diff --git a/src/celmath/ray.h b/src/celmath/ray.h index fa688dd67..9e858c381 100644 --- a/src/celmath/ray.h +++ b/src/celmath/ray.h @@ -7,63 +7,31 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_RAY_H_ -#define _CELMATH_RAY_H_ +#pragma once #include +#include namespace celmath { +// Eigen 3.4 will support a transform method on ParametrizedLine, which may enable +// removing at least one of the below overloads. -template class Ray3 -{ - public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW - - Ray3(); - Ray3(const Eigen::Matrix& origin, const Eigen::Matrix& direction); - - Eigen::Matrix point(T) const; - - Ray3 transform(const Eigen::Matrix& m) const - { - return Ray3(m * origin, m * direction); - } - - Ray3 transform(const Eigen::Matrix& m) const - { - Eigen::Matrix o(Eigen::Matrix::Ones()); - o.head(3) = origin; - Eigen::Matrix d(Eigen::Matrix::Zero()); - d.head(3) = direction; - return Ray3((m * o).head(3), (m * d).head(3)); - } - - public: - Eigen::Matrix origin; - Eigen::Matrix direction; -}; - -typedef Ray3 Ray3f; -typedef Ray3 Ray3d; - - -template Ray3::Ray3() : - origin(0, 0, 0), direction(0, 0, -1) +template +Eigen::ParametrizedLine transformRay(const Eigen::ParametrizedLine& line, + const Eigen::Matrix& m) { + return Eigen::ParametrizedLine(m * line.origin(), m * line.direction()); } -template Ray3::Ray3(const Eigen::Matrix& _origin, - const Eigen::Matrix& _direction) : - origin(_origin), direction(_direction) +template +Eigen::ParametrizedLine transformRay(const Eigen::ParametrizedLine& line, + const Eigen::Matrix& m) { + Eigen::Matrix o(Eigen::Matrix::Ones()); + o.head(3) = line.origin(); + Eigen::Matrix d(Eigen::Matrix::Zero()); + d.head(3) = line.direction(); + return Eigen::ParametrizedLine((m * o).head(3), (m * d).head(3)); } - -template Eigen::Matrix Ray3::point(T t) const -{ - return origin + direction * t; -} - -} // namespace celmath - -#endif // _CELMATH_RAY_H_ +} // end namespace celmath diff --git a/src/celmath/solve.h b/src/celmath/solve.h index 9aa98f3ba..b43f05eec 100644 --- a/src/celmath/solve.h +++ b/src/celmath/solve.h @@ -7,6 +7,7 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. +#include #include #pragma once @@ -16,12 +17,13 @@ namespace celmath // Solve a function using the bisection method. Returns a pair // with the solution as the first element and the error as the second. -template std::pair solve_bisection(F f, - T lower, T upper, - T err, - int maxIter = 100) +template +std::pair solve_bisection(F f, + T lower, T upper, + T err, + int maxIter = 100) { - T x = 0.0; + T x = static_cast(0); for (int i = 0; i < maxIter; i++) { @@ -30,7 +32,7 @@ template std::pair solve_bisection(F f, break; T y = f(x); - if (y < 0) + if (y < static_cast(0)) lower = x; else upper = x; @@ -42,12 +44,15 @@ template std::pair solve_bisection(F f, // Solve using iteration; terminate when error is below err or the maximum // number of iterations is reached. -template std::pair solve_iteration(F f, - T x0, - T err, - int maxIter = 100) +template +std::pair solve_iteration(F f, + T x0, + T err, + int maxIter = 100) { - T x = 0; + using std::abs; + + T x = static_cast(0); T x2 = x0; for (int i = 0; i < maxIter; i++) @@ -63,11 +68,12 @@ template std::pair solve_iteration(F f, // Solve using iteration method and a fixed number of steps. -template std::pair solve_iteration_fixed(F f, - T x0, - int maxIter) +template +std::pair solve_iteration_fixed(F f, + T x0, + int maxIter) { - T x = 0; + T x = static_cast(0); T x2 = x0; for (int i = 0; i < maxIter; i++) diff --git a/src/celmath/sphere.h b/src/celmath/sphere.h index dbf76aa90..f3f07686b 100644 --- a/src/celmath/sphere.h +++ b/src/celmath/sphere.h @@ -7,8 +7,7 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. -#ifndef _CELMATH_SPHERE_H_ -#define _CELMATH_SPHERE_H_ +#pragma once #include @@ -21,13 +20,13 @@ template class Sphere EIGEN_MAKE_ALIGNED_OPERATOR_NEW Sphere() : - center(0, 0, 0), - radius(1) + center(Eigen::Matrix::Zero()), + radius(static_cast(1)) { } Sphere(T _radius) : - center(0, 0, 0), + center(Eigen::Matrix::Zero()), radius(_radius) { } @@ -43,9 +42,6 @@ template class Sphere T radius; }; -typedef Sphere Spheref; -typedef Sphere Sphered; - +using Spheref = Sphere; +using Sphered = Sphere; } // namespace celmath - -#endif // _CELMATH_SPHERE_H_ diff --git a/src/tools/atmosphere/scattersim.cpp b/src/tools/atmosphere/scattersim.cpp index 91ed5cab5..081ae473e 100644 --- a/src/tools/atmosphere/scattersim.cpp +++ b/src/tools/atmosphere/scattersim.cpp @@ -15,6 +15,10 @@ #include #include #include + +#include +#include + #include #include #include @@ -122,7 +126,8 @@ public: Camera() = default; - Ray3d getViewRay(double viewportX, double viewportY) const + Eigen::ParametrizedLine + getViewRay(double viewportX, double viewportY) const { Vector3d viewDir; @@ -144,8 +149,8 @@ public: viewDir.normalize(); } - Ray3d viewRay(Vector3d::Zero(), viewDir); - return viewRay.transform(transform); + Eigen::ParametrizedLine viewRay(Vector3d::Zero(), viewDir); + return transformRay(viewRay, transform); } double fov{PI / 2.0}; @@ -285,8 +290,8 @@ public: void setParameters(ParameterSet& params); - Color raytrace(const Ray3d& ray) const; - Color raytrace_LUT(const Ray3d& ray) const; + Color raytrace(const Eigen::ParametrizedLine& ray) const; + Color raytrace_LUT(const Eigen::ParametrizedLine& ray) const; Color background; Light light; @@ -654,15 +659,15 @@ LUT3::lookup(double x, double y, double z) const } -template bool raySphereIntersect(const Ray3& ray, +template bool raySphereIntersect(const Eigen::ParametrizedLine& ray, const Sphere& sphere, T& dist0, T& dist1) { - Matrix diff = ray.origin - sphere.center; + Matrix diff = ray.origin() - sphere.center; T s = (T) 1.0 / square(sphere.radius); - T a = ray.direction.dot(ray.direction) * s; - T b = ray.direction.dot(diff) * s; + T a = ray.direction().squaredNorm() * s; + T b = ray.direction().dot(diff) * s; T c = diff.dot(diff) * s - (T) 1.0; T disc = b * b - a * c; if (disc < 0.0) @@ -695,15 +700,15 @@ template bool raySphereIntersect(const Ray3& ray, } -template bool raySphereIntersect2(const Ray3& ray, +template bool raySphereIntersect2(const Eigen::ParametrizedLine& ray, const Sphere& sphere, T& dist0, T& dist1) { - Matrix diff = ray.origin - sphere.center; + Matrix diff = ray.origin() - sphere.center; T s = (T) 1.0 / square(sphere.radius); - T a = ray.direction.dot(ray.direction) * s; - T b = ray.direction.dot(diff) * s; + T a = ray.direction().squaredNorm() * s; + T b = ray.direction().dot(diff) * s; T c = diff.dot(diff) * s - (T) 1.0; T disc = b * b - a * c; if (disc < 0.0) @@ -913,12 +918,12 @@ Vector3d integrateInscattering(const Scene& scene, for (unsigned int i = 0; i < nSteps; i++) { - Ray3d sunRay(samplePoint, lightDir); + Eigen::ParametrizedLine sunRay(samplePoint, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); // Compute the optical depth along path from sample point to the sun - OpticalDepths sunDepth = integrateOpticalDepth(scene, samplePoint, sunRay.point(sunDist)); + OpticalDepths sunDepth = integrateOpticalDepth(scene, samplePoint, sunRay.pointAt(sunDist)); // Compute the optical depth along the path from the sample point to the eye OpticalDepths eyeDepth = integrateOpticalDepth(scene, samplePoint, atmStart); @@ -971,12 +976,12 @@ Vector4d integrateInscatteringFactors(const Scene& scene, for (unsigned int i = 0; i < nSteps; i++) { - Ray3d sunRay(samplePoint, lightDir); + Eigen::ParametrizedLine sunRay(samplePoint, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); // Compute the optical depth along path from sample point to the sun - OpticalDepths sunDepth = integrateOpticalDepth(scene, samplePoint, sunRay.point(sunDist)); + OpticalDepths sunDepth = integrateOpticalDepth(scene, samplePoint, sunRay.pointAt(sunDist)); // Compute the optical depth along the path from the sample point to the eye OpticalDepths eyeDepth = integrateOpticalDepth(scene, samplePoint, atmStart); @@ -1043,14 +1048,14 @@ buildExtinctionLUT(const Scene& scene) double sinAngle = sqrt(1.0 - min(1.0, cosAngle * cosAngle)); Vector3d viewDir(cosAngle, sinAngle, 0.0); - Ray3d ray(atmStart, viewDir); + Eigen::ParametrizedLine ray(atmStart, viewDir); double dist = 0.0; if (!testIntersection(ray, shell, dist)) dist = 0.0; OpticalDepths depth = integrateOpticalDepth(scene, atmStart, - ray.point(dist)); + ray.pointAt(dist)); depth.rayleigh *= 4.0 * PI; depth.mie *= 4.0 * PI; Vector3d ext = scene.atmosphere.computeExtinction(depth); @@ -1103,14 +1108,14 @@ buildOpticalDepthLUT(const Scene& scene) double sinAngle = sqrt(1.0 - min(1.0, cosAngle * cosAngle)); Vector3d dir(cosAngle, sinAngle, 0.0); - Ray3d ray(atmStart, dir); + Eigen::ParametrizedLine ray(atmStart, dir); double dist = 0.0; if (!testIntersection(ray, shell, dist)) dist = 0.0; OpticalDepths depth = integrateOpticalDepth(scene, atmStart, - ray.point(dist)); + ray.pointAt(dist)); depth.rayleigh *= 4.0 * PI; depth.mie *= 4.0 * PI; @@ -1172,11 +1177,11 @@ Vector3d integrateInscattering_LUT(const Scene& scene, for (unsigned int i = 0; i < nSteps; i++) { - Ray3d sunRay(samplePoint, lightDir); + Eigen::ParametrizedLine sunRay(samplePoint, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); - Vector3d sunExt = lookupExtinction(scene, samplePoint, sunRay.point(sunDist)); + Vector3d sunExt = lookupExtinction(scene, samplePoint, sunRay.pointAt(sunDist)); Vector3d eyeExt; if (!eyeInsideAtmosphere) { @@ -1254,11 +1259,11 @@ Vector4d integrateInscatteringFactors_LUT(const Scene& scene, for (unsigned int i = 0; i < nSteps; i++) { - Ray3d sunRay(samplePoint, lightDir); + Eigen::ParametrizedLine sunRay(samplePoint, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); - Vector3d sunExt = lookupExtinction(scene, samplePoint, sunRay.point(sunDist)); + Vector3d sunExt = lookupExtinction(scene, samplePoint, sunRay.pointAt(sunDist)); Vector3d eyeExt; Vector3d subExt; if (planetHit) @@ -1320,12 +1325,12 @@ buildScatteringLUT(const Scene& scene) double sinAngle = sqrt(1.0 - min(1.0, cosAngle * cosAngle)); Vector3d viewDir(cosAngle, sinAngle, 0.0); - Ray3d viewRay(atmStart, viewDir); + Eigen::ParametrizedLine viewRay(atmStart, viewDir); double dist = 0.0; if (!testIntersection(viewRay, shell, dist)) dist = 0.0; - Vector3d atmEnd = viewRay.point(dist); + Vector3d atmEnd = viewRay.pointAt(dist); for (unsigned int k = 0; k < ScatteringLUTLightAngleSteps; k++) { @@ -1392,7 +1397,7 @@ Color getPlanetColor(const Scene& scene, const Vector3d& p) } -Color Scene::raytrace(const Ray3d& ray) const +Color Scene::raytrace(const Eigen::ParametrizedLine& ray) const { double dist = 0.0; double atmEnter = 0.0; @@ -1401,7 +1406,7 @@ Color Scene::raytrace(const Ray3d& ray) const double shellRadius = planet.radius + atmosphereShellHeight; Color color = background; - if (ray.direction.dot(-light.direction) > cos(sunAngularDiameter / 2.0)) + if (ray.direction().dot(-light.direction) > cos(sunAngularDiameter / 2.0)) color = light.color; if (raySphereIntersect(ray, @@ -1410,12 +1415,12 @@ Color Scene::raytrace(const Ray3d& ray) const atmExit)) { Color baseColor = color; - Vector3d atmStart = ray.origin + atmEnter * ray.direction; - Vector3d atmEnd = ray.origin + atmExit * ray.direction; + Vector3d atmStart = ray.origin() + atmEnter * ray.direction(); + Vector3d atmEnd = ray.origin() + atmExit * ray.direction(); if (testIntersection(ray, planet, dist)) { - Vector3d intersectPoint = ray.point(dist); + Vector3d intersectPoint = ray.pointAt(dist); Vector3d normal = intersectPoint - planet.center; normal.normalize(); Vector3d lightDir = -light.direction; @@ -1428,13 +1433,13 @@ Color Scene::raytrace(const Ray3d& ray) const // Compute ray from surface point to edge of the atmosphere in the direction // of the sun. - Ray3d sunRay(surfacePt, lightDir); + Eigen::ParametrizedLine sunRay(surfacePt, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); // Compute color of sunlight filtered by the atmosphere; consider extinction // along both the sun-to-surface and surface-to-eye paths. - OpticalDepths sunDepth = integrateOpticalDepth(*this, surfacePt, sunRay.point(sunDist)); + OpticalDepths sunDepth = integrateOpticalDepth(*this, surfacePt, sunRay.pointAt(sunDist)); OpticalDepths eyeDepth = integrateOpticalDepth(*this, atmStart, surfacePt); OpticalDepths totalDepth = sumOpticalDepths(sunDepth, eyeDepth); totalDepth.rayleigh *= 4.0 * PI; @@ -1445,7 +1450,7 @@ Color Scene::raytrace(const Ray3d& ray) const // surface color * sun color * atmospheric extinction baseColor = (planetColor * extinction) * light.color * diffuse; - atmEnd = ray.origin + dist * ray.direction; + atmEnd = ray.origin() + dist * ray.direction(); } Vector3d inscatter = integrateInscattering(*this, atmStart, atmEnd) * 4.0 * PI; @@ -1459,21 +1464,21 @@ Color Scene::raytrace(const Ray3d& ray) const Color -Scene::raytrace_LUT(const Ray3d& ray) const +Scene::raytrace_LUT(const Eigen::ParametrizedLine& ray) const { double atmEnter = 0.0; double atmExit = 0.0; double shellRadius = planet.radius + atmosphereShellHeight; Sphered shell(shellRadius); - Vector3d eyePt = Vector3d::Zero() + (ray.origin - planet.center); + Vector3d eyePt = Vector3d::Zero() + (ray.origin() - planet.center); Color color = background; - if (ray.direction.dot(-light.direction) > cos(sunAngularDiameter / 2.0)) + if (ray.direction().dot(-light.direction) > cos(sunAngularDiameter / 2.0)) color = light.color; // Transform ray to model space - Ray3d mray(eyePt, ray.direction); + Eigen::ParametrizedLine mray(eyePt, ray.direction()); bool hit = raySphereIntersect2(mray, shell, atmEnter, atmExit); if (hit && atmExit > 0.0) @@ -1481,18 +1486,18 @@ Scene::raytrace_LUT(const Ray3d& ray) const Color baseColor = color; bool eyeInsideAtmosphere = atmEnter < 0.0; - Vector3d atmStart = mray.origin + atmEnter * mray.direction; - Vector3d atmEnd = mray.origin + atmExit * mray.direction; + Vector3d atmStart = mray.origin() + atmEnter * mray.direction(); + Vector3d atmEnd = mray.origin() + atmExit * mray.direction(); double planetEnter = 0.0; double planetExit = 0.0; hit = raySphereIntersect2(mray, - Sphered(planet.radius), - planetEnter, - planetExit); + Sphered(planet.radius), + planetEnter, + planetExit); if (hit && planetEnter > 0.0) { - Vector3d surfacePt = mray.point(planetEnter); + Vector3d surfacePt = mray.pointAt(planetEnter); // Lambert lighting Vector3d normal = surfacePt - Vector3d::Zero(); @@ -1504,13 +1509,13 @@ Scene::raytrace_LUT(const Ray3d& ray) const // Compute ray from surface point to edge of the atmosphere in the direction // of the sun. - Ray3d sunRay(surfacePt, lightDir); + Eigen::ParametrizedLine sunRay(surfacePt, lightDir); double sunDist = 0.0; testIntersection(sunRay, shell, sunDist); // Compute color of sunlight filtered by the atmosphere; consider extinction // along both the sun-to-surface and surface-to-eye paths. - Vector3d sunExt = lookupExtinction(*this, surfacePt, sunRay.point(sunDist)); + Vector3d sunExt = lookupExtinction(*this, surfacePt, sunRay.pointAt(sunDist)); Vector3d eyeExt = lookupExtinction(*this, surfacePt, atmStart); if (eyeInsideAtmosphere) { @@ -1524,7 +1529,7 @@ Scene::raytrace_LUT(const Ray3d& ray) const // surface color * sun color * atmospheric extinction baseColor = (planetColor * extinction) * light.color * diffuse; - atmEnd = mray.point(planetEnter); + atmEnd = mray.pointAt(planetEnter); } Vector3d inscatter; @@ -1550,7 +1555,7 @@ Scene::raytrace_LUT(const Ray3d& ray) const else { atmEnd = atmStart; - atmStart = mray.point(planetEnter); + atmStart = mray.pointAt(planetEnter); //cout << atmEnter << ", " << planetEnter << ", " << atmExit << "\n"; rayleighScatter = lookupScattering(*this, atmStart, atmEnd, -light.direction) - @@ -1572,7 +1577,7 @@ Scene::raytrace_LUT(const Ray3d& ray) const } const Vector3d& rayleigh = atmosphere.rayleighCoeff; - double cosSunAngle = mray.direction.dot(-light.direction); + double cosSunAngle = mray.direction().dot(-light.direction); inscatter = phaseRayleigh(cosSunAngle) * rayleighScatter.cwiseProduct(rayleigh); inscatter = inscatter * 4.0 * PI; } @@ -1613,7 +1618,7 @@ void render(const Scene& scene, double viewportX = ((double) (j - viewport.x) / (double) (viewport.width - 1) - 0.5) * aspectRatio; double viewportY ((double) (i - viewport.y) / (double) (viewport.height - 1) - 0.5); - Ray3d viewRay = camera.getViewRay(viewportX, viewportY); + Eigen::ParametrizedLine viewRay = camera.getViewRay(viewportX, viewportY); Color color; if (LUTUsage != NoLUT)