Index: celestia.vcproj =================================================================== --- celestia.vcproj (revision 4824) +++ celestia.vcproj (working copy) @@ -37,7 +37,7 @@ /> + + + + #include #include +#include #include #include #include @@ -1733,6 +1734,7 @@ } +#if 0 class OrbitSampler : public OrbitSampleProc { public: @@ -1747,6 +1749,22 @@ samples->push_back(samp); }; }; +#endif +class OrbitSampler : public OrbitSampleProc +{ +public: + CurvePlot* m_orbitPath; + + OrbitSampler(CurvePlot* orbitPath) : m_orbitPath(orbitPath) {}; + void sample(double t, const Vector3d& position, const Vector3d& velocity) + { + CurvePlotSample samp; + samp.position = position; + samp.velocity = velocity; + samp.t = t; + m_orbitPath->addSample(samp); + }; +}; void renderOrbitColor(const Body *body, bool selected, float opacity) @@ -1808,6 +1826,7 @@ } +#if 0 // Subdivide the orbit into sections and compute a bounding volume for each section. The bounding // volumes used are capsules, the set of all points less than some constant distance from a line // segment. @@ -1885,12 +1904,230 @@ ((3.0 * (p1 - p0) - 2.0 * v0 - v1) * (t * t)) + (v0 * t)); } +#endif -static int splinesRendered = 0; static int orbitsRendered = 0; static int orbitsSkipped = 0; static int sectionsCulled = 0; + +void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath, + double t, + const Quatf& cameraOrientationf, + const Frustum& frustum, + float nearDist, + float farDist) +{ + Body* body = orbitPath.body; + Quatd cameraOrientation(cameraOrientationf.w, cameraOrientationf.x, cameraOrientationf.y, cameraOrientationf.z); + double nearZ = -nearDist; // negate, becase z is into the screen in camera space + double farZ = -farDist; + + const Orbit* orbit = NULL; + if (body != NULL) + orbit = body->getOrbit(t); + else + orbit = orbitPath.star->getOrbit(); + + CurvePlot* cachedOrbit = NULL; + OrbitCache::iterator cached = orbitCache.find(orbit); + if (cached != orbitCache.end()) + { + cachedOrbit = cached->second; + cachedOrbit->setLastUsed(frameCount); + } + + // If it's not in the cache already + if (cachedOrbit == NULL) + { + double startTime = t; + int nSamples = detailOptions.orbitPathSamplePoints; + + // Adjust the number of samples used for aperiodic orbits--these aren't + // true orbits, but are sampled trajectories, generally of spacecraft. + // Better control is really needed--some sort of adaptive sampling would + // be ideal. + if (!orbit->isPeriodic()) + { + double begin = 0.0, end = 0.0; + orbit->getValidRange(begin, end); + + if (begin != end) + { + startTime = begin; + nSamples = (int) (orbit->getPeriod() * 100.0); + nSamples = max(min(nSamples, 1000), 100); + } + else + { + // If the orbit is aperiodic and doesn't have a + // finite duration, we don't render it. A compromise + // would be to pick some time window centered at the + // current time, but we'd have to pick some arbitrary + // duration. + nSamples = 0; + } + } + else + { + startTime = t - orbit->getPeriod() / 2.0; + } + + cachedOrbit = new CurvePlot(); + cachedOrbit->setLastUsed(frameCount); + + OrbitSampler sampler(cachedOrbit); + orbit->sample(startTime, + orbit->getPeriod(), + nSamples, + sampler); + + // If the orbit cache is full, first try and eliminate some old orbits + if (orbitCache.size() > OrbitCacheCullThreshold) + { + // Check for old orbits at most once per frame + if (lastOrbitCacheFlush != frameCount) + { + for (OrbitCache::iterator iter = orbitCache.begin(); iter != orbitCache.end();) + { + // Tricky code to eliminate a node in the orbit cache without screwing + // up the iterator. Should work in all STL implementations. + if (frameCount - iter->second->lastUsed() > OrbitCacheRetireAge) + orbitCache.erase(iter++); + else + ++iter; + } + lastOrbitCacheFlush = frameCount; + } + } + + orbitCache[orbit] = cachedOrbit; + } + + if (cachedOrbit->empty()) + return; + + // 'Periodic' orbits are generally not strictly periodic because of perturbations + // from other bodies. Here we update the trajectory samples to make sure that the + // orbit covers a time range centered at the current time and covering a full revolution. + if (orbit->isPeriodic()) + { + double startTime = t - orbit->getPeriod() / 2.0; + double endTime = t + orbit->getPeriod() / 2.0; + double dt = orbit->getPeriod() / detailOptions.orbitPathSamplePoints; + + if (startTime < cachedOrbit->startTime()) + { + cachedOrbit->removeSamplesAfter(endTime + dt); + + double orbitStartTime = cachedOrbit->empty() ? endTime : cachedOrbit->startTime() - dt; + do { + CurvePlotSample sample; + sample.t = orbitStartTime; + sample.position = orbit->positionAtTime(orbitStartTime); + sample.velocity = orbit->velocityAtTime(orbitStartTime); + cachedOrbit->addSample(sample); + orbitStartTime -= dt; + } while (orbitStartTime > startTime); + } + else if (endTime > cachedOrbit->endTime()) + { + cachedOrbit->removeSamplesBefore(startTime - dt); + + double orbitEndTime = cachedOrbit->empty() ? startTime : cachedOrbit->endTime() + dt; + do { + CurvePlotSample sample; + sample.t = orbitEndTime; + sample.position = orbit->positionAtTime(orbitEndTime); + sample.velocity = orbit->velocityAtTime(orbitEndTime); + cachedOrbit->addSample(sample); + orbitEndTime += dt; + } while (orbitEndTime < endTime); + } + } + + // We perform vertex tranformations on the CPU because double precision is necessary to + // render orbits properly. Start by computing the modelview matrix, to transform orbit + // vertices into camera space. + Mat4d modelview; + { + Quatd orientation(1.0); + if (body != NULL) + { + orientation = fromEigen(body->getOrbitFrame(t)->getOrientation(t)); + } + + // Equivalent to: + // glRotate(cameraOrientation); + // glTranslate(orbitPath.origin); + // glRotate(~orientation); + modelview = + (orientation).toMatrix4() * + Mat4d::translation(ptFromEigen(orbitPath.origin)) * + (~cameraOrientation).toMatrix4(); + } + + glPushMatrix(); + glLoadIdentity(); + + bool highlight; + if (body != NULL) + highlight = highlightObject.body() == body; + else + highlight = highlightObject.star() == orbitPath.star; + renderOrbitColor(body, highlight, orbitPath.opacity); + +#ifdef STIPPLED_LINES + glLineStipple(3, 0x5555); + glEnable(GL_LINE_STIPPLE); +#endif + + double subdivisionThreshold = pixelSize * 40.0; + Eigen::Transform3d eigen_modelview; + eigen_modelview.matrix() = Eigen::Matrix4d(&modelview[0][0]); + + Eigen::Vector3d viewFrustumPlaneNormals[4]; + for (int i = 0; i < 4; i++) + { + Vec3f n = frustum.getPlane(i).normal; + viewFrustumPlaneNormals[i] = Eigen::Vector3d(n.x, n.y, n.z); + } + + if (orbit->isPeriodic()) + { + cachedOrbit->render(eigen_modelview, + nearZ, farZ, viewFrustumPlaneNormals, + subdivisionThreshold, + t - orbit->getPeriod() / 2.0, t + orbit->getPeriod() / 2.0); + } + else + { + if (renderFlags & ShowPartialTrajectories) + { + // Show the trajectory from the start time until the current simulation time + cachedOrbit->render(eigen_modelview, + nearZ, farZ, viewFrustumPlaneNormals, + subdivisionThreshold, + cachedOrbit->startTime(), t); + } + else + { + // Show the entire trajectory + cachedOrbit->render(eigen_modelview, + nearZ, farZ, viewFrustumPlaneNormals, + subdivisionThreshold); + } + } + +#ifdef STIPPLED_LINES + glDisable(GL_LINE_STIPPLE); +#endif + + glPopMatrix(); +} + +#if 0 +static int splinesRendered = 0; static Point3d renderOrbitSplineSegment(const Renderer::OrbitSample& p0, const Renderer::OrbitSample& p1, const Renderer::OrbitSample& p2, @@ -2455,6 +2692,7 @@ glPopMatrix(); } +#endif // Convert a position in the universal coordinate system to astrocentric @@ -4016,19 +4254,18 @@ for (i = 0; i < (int) orbitPathList.size(); i++) { const OrbitPathListEntry& o = orbitPathList[i]; - float minNearDistance = min(-o.radius * 0.0001f, o.centerZ + o.radius); + float minNearDistance = min(-MinNearPlaneDistance, o.centerZ + o.radius); if (minNearDistance > zNearest) zNearest = minNearDistance; } - // Adjust the nearest interval to include the closest marker (if it's // closer to the observer than anything else if (!depthSortedAnnotations.empty()) { // Factor of 0.999 makes sure ensures that the near plane does not fall // exactly at the marker's z coordinate (in which case the marker - // would be susceptible to being clipped.) + // would be susceptible to getting clipped.) if (-depthSortedAnnotations[0].position.z() > zNearest) zNearest = -depthSortedAnnotations[0].position.z() * 0.999f; } @@ -4198,13 +4435,8 @@ float farZ = -orbitIter->centerZ + orbitIter->radius; // Don't render orbits when they're completely outside this - // depth interval. Also, don't render an orbit in this - // interval if it is vastly larger than the interval - // range; otherwise, the GPU will have precision troubles - // when clipping, producing visual artifacts. The factor - // of 1e5 may need some tuning. - if (nearZ < farPlaneDistance && farZ > nearPlaneDistance && - orbitIter->radius < 1.0e8f * (farPlaneDistance - nearPlaneDistance)) + // depth interval. + if (nearZ < farPlaneDistance && farZ > nearPlaneDistance) { #ifdef DEBUG_COALESCE switch (interval % 6) @@ -4259,12 +4491,10 @@ #if 0 // TODO: Debugging output for new orbit code; remove when development is complete clog << "orbits: " << orbitsRendered - << ", splines: " << splinesRendered << ", skipped: " << orbitsSkipped << ", sections culled: " << sectionsCulled << ", nIntervals: " << nIntervals << "\n"; #endif - splinesRendered = 0; orbitsRendered = 0; orbitsSkipped = 0; sectionsCulled = 0; @@ -8808,7 +9038,7 @@ path.star = NULL; path.centerZ = origf * viewMatZ; path.radius = (float) boundingRadius; - path.origin = toEigen(origf); + path.origin = toEigen(relOrigin); path.opacity = sizeFade(orbitRadiusInPixels, minOrbitSize, 2.0f); orbitPathList.push_back(path); } @@ -9032,7 +9262,7 @@ path.body = NULL; path.centerZ = origf * viewMatZ; path.radius = (float) boundingRadius; - path.origin = toEigen(origf); + path.origin = toEigen(origf).cast(); path.opacity = sizeFade(orbitRadiusInPixels, minOrbitSize, 2.0f); orbitPathList.push_back(path); } Index: src/celengine/render.h =================================================================== --- src/celengine/render.h (revision 4824) +++ src/celengine/render.h (working copy) @@ -28,6 +28,7 @@ class RendererWatcher; class FrameTree; class ReferenceMark; +class CurvePlot; struct LightSource { @@ -307,7 +308,7 @@ float radius; Body* body; const Star* star; - Eigen::Vector3f origin; + Eigen::Vector3d origin; float opacity; bool operator<(const OrbitPathListEntry&) const; @@ -683,6 +684,7 @@ public: +#if 0 struct OrbitSample { double t; @@ -704,9 +706,10 @@ std::vector sections; uint32 lastUsed; }; +#endif private: - typedef std::map OrbitCache; + typedef std::map OrbitCache; OrbitCache orbitCache; uint32 lastOrbitCacheFlush; Index: src/celengine/samporbit.cpp =================================================================== --- src/celengine/samporbit.cpp (revision 4824) +++ src/celengine/samporbit.cpp (working copy) @@ -395,12 +395,44 @@ template void SampledOrbit::sample(double start, double t, int, OrbitSampleProc& proc) const { + for (unsigned int i = 0; i < samples.size(); i++) + { + Vector3d v; + Vector3d p(samples[i].x, samples[i].y, samples[i].z); + + if (samples.size() == 1) + { + v = Vector3d::Zero(); + } + else if (i == 0) + { + double dt = samples[i + 1].t - samples[i].t; + v = (Vector3d(samples[i + 1].x, samples[i + 1].y, samples[i + 1].z) - p) / dt; + } + else if (i == samples.size() - 1) + { + double dt = samples[i].t - samples[i - 1].t; + v = (p - Vector3d(samples[i - 1].x, samples[i - 1].y, samples[i - 1].z)) / dt; + } + else + { + double dt0 = samples[i + 1].t - samples[i].t; + Vector3d v0 = (Vector3d(samples[i + 1].x, samples[i + 1].y, samples[i + 1].z) - p) / dt0; + double dt1 = samples[i].t - samples[i - 1].t; + Vector3d v1 = (p - Vector3d(samples[i - 1].x, samples[i - 1].y, samples[i - 1].z)) / dt1; + v = (v0 + v1) * 0.5; + } + + proc.sample(samples[i].t, Vector3d(p.x(), p.z(), -p.y()), Vector3d(v.x(), v.z(), -v.y())); + } + +#if 0 double cosThresholdAngle = cos(degToRad(SampleThresholdAngle)); double dt = MinSampleInterval; double end = start + t; double current = start; - proc.sample(current, positionAtTime(current)); + //proc.sample(samples[i].t, Vector3d(p.x(), p.z(), -p.y()), Vector3d(v.x(), v.z(), -v.y())); while (current < end) { @@ -433,6 +465,7 @@ } current += gooddt; } +#endif } @@ -677,6 +710,15 @@ template void SampledOrbitXYZV::sample(double start, double t, int, OrbitSampleProc& proc) const { + for (typename vector >::const_iterator iter = samples.begin(); + iter != samples.end(); iter++) + { + proc.sample(iter->t, + Vector3d(iter->position.x(), iter->position.z(), -iter->position.y()), + Vector3d(iter->velocity.x(), iter->velocity.z(), -iter->velocity.y())); + } + +#if 0 double cosThresholdAngle = cos(degToRad(SampleThresholdAngle)); double dt = MinSampleInterval; double end = start + t; @@ -720,6 +762,7 @@ } current += gooddt; } +#endif } Index: src/celestia.pro =================================================================== --- src/celestia.pro (revision 4824) +++ src/celestia.pro (working copy) @@ -377,9 +377,15 @@ ../thirdparty/glew/include/GL/glxew.h \ ../thirdparty/glew/include/GL/wglew.h -THIRDPARTY_SOURCES = $$GLEW_SOURCES -THIRDPARTY_HEADERS = $$GLEW_HEADERS +CURVEPLOT_SOURCES = \ + ../thirdparty/curveplot/src/curveplot.cpp +CURVEPLOT_HEADERS = \ + ../thirdparty/curveplot/include/curveplot.h + +THIRDPARTY_SOURCES = $$GLEW_SOURCES $$CURVEPLOT_SOURCES +THIRDPARTY_HEADERS = $$GLEW_HEADERS $$CURVEPLOT_HEADERS + DEFINES += GLEW_STATIC @@ -428,6 +434,7 @@ # Third party libraries INCLUDEPATH += ../thirdparty/glew/include INCLUDEPATH += ../thirdparty/Eigen +INCLUDEPATH += ../thirdparty/curveplot/include win32 { INCLUDEPATH += \ Index: thirdparty/curveplot/include/curveplot.h =================================================================== --- thirdparty/curveplot/include/curveplot.h (revision 0) +++ thirdparty/curveplot/include/curveplot.h (revision 0) @@ -0,0 +1,99 @@ +// curveplot.h +// +// Copyright (C) 2009 Chris Laurel . +// +// curveplot is a module for rendering curves in OpenGL at high precision. A +// plot is a series of cubic curves. The curves are transformed +// to camera space in software because double precision is absolutely +// required. The cubics are adaptively subdivided based on distance from +// the camera position. +// +// curveplot is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// curveplot is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// orbitpath. If not, see . + +#include +#include + + +class HighPrec_Frustum; + +class CurvePlotSample +{ +public: + Eigen::Vector3d position; + double t; + Eigen::Vector3d velocity; + double boundingRadius; +}; + + +class CurvePlot +{ + public: + CurvePlot(); + + double duration() const { return m_duration; } + void setDuration(double duration); + + double startTime() const + { + if (m_samples.empty()) + return 0.0; + else + return m_samples.front().t; + } + + double endTime() const + { + if (m_samples.empty()) + return 0.0; + else + return m_samples.back().t; + } + + void render(const Eigen::Transform3d& modelview, + double nearZ, + double farZ, + const Eigen::Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold) const; + void render(const Eigen::Transform3d& modelview, + double nearZ, + double farZ, + const Eigen::Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold, + double startTime, + double endTime) const; + + unsigned int lastUsed() const { return m_lastUsed; } + void setLastUsed(unsigned int lastUsed) { m_lastUsed = lastUsed; } + + void addSample(const CurvePlotSample& sample); + void removeSamplesBefore(double t); + void removeSamplesAfter(double t); + + bool empty() const { return m_samples.empty(); } + + private: + std::deque m_samples; + + double m_duration; + + unsigned int m_lastUsed; +}; + Index: thirdparty/curveplot/LICENSE =================================================================== --- thirdparty/curveplot/LICENSE (revision 0) +++ thirdparty/curveplot/LICENSE (revision 0) @@ -0,0 +1,165 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + + This version of the GNU Lesser General Public License incorporates +the terms and conditions of version 3 of the GNU General Public +License, supplemented by the additional permissions listed below. + + 0. Additional Definitions. + + As used herein, "this License" refers to version 3 of the GNU Lesser +General Public License, and the "GNU GPL" refers to version 3 of the GNU +General Public License. + + "The Library" refers to a covered work governed by this License, +other than an Application or a Combined Work as defined below. + + An "Application" is any work that makes use of an interface provided +by the Library, but which is not otherwise based on the Library. +Defining a subclass of a class defined by the Library is deemed a mode +of using an interface provided by the Library. + + A "Combined Work" is a work produced by combining or linking an +Application with the Library. The particular version of the Library +with which the Combined Work was made is also called the "Linked +Version". + + The "Minimal Corresponding Source" for a Combined Work means the +Corresponding Source for the Combined Work, excluding any source code +for portions of the Combined Work that, considered in isolation, are +based on the Application, and not on the Linked Version. + + The "Corresponding Application Code" for a Combined Work means the +object code and/or source code for the Application, including any data +and utility programs needed for reproducing the Combined Work from the +Application, but excluding the System Libraries of the Combined Work. + + 1. Exception to Section 3 of the GNU GPL. + + You may convey a covered work under sections 3 and 4 of this License +without being bound by section 3 of the GNU GPL. + + 2. Conveying Modified Versions. + + If you modify a copy of the Library, and, in your modifications, a +facility refers to a function or data to be supplied by an Application +that uses the facility (other than as an argument passed when the +facility is invoked), then you may convey a copy of the modified +version: + + a) under this License, provided that you make a good faith effort to + ensure that, in the event an Application does not supply the + function or data, the facility still operates, and performs + whatever part of its purpose remains meaningful, or + + b) under the GNU GPL, with none of the additional permissions of + this License applicable to that copy. + + 3. Object Code Incorporating Material from Library Header Files. + + The object code form of an Application may incorporate material from +a header file that is part of the Library. You may convey such object +code under terms of your choice, provided that, if the incorporated +material is not limited to numerical parameters, data structure +layouts and accessors, or small macros, inline functions and templates +(ten or fewer lines in length), you do both of the following: + + a) Give prominent notice with each copy of the object code that the + Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the object code with a copy of the GNU GPL and this license + document. + + 4. Combined Works. + + You may convey a Combined Work under terms of your choice that, +taken together, effectively do not restrict modification of the +portions of the Library contained in the Combined Work and reverse +engineering for debugging such modifications, if you also do each of +the following: + + a) Give prominent notice with each copy of the Combined Work that + the Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the Combined Work with a copy of the GNU GPL and this license + document. + + c) For a Combined Work that displays copyright notices during + execution, include the copyright notice for the Library among + these notices, as well as a reference directing the user to the + copies of the GNU GPL and this license document. + + d) Do one of the following: + + 0) Convey the Minimal Corresponding Source under the terms of this + License, and the Corresponding Application Code in a form + suitable for, and under terms that permit, the user to + recombine or relink the Application with a modified version of + the Linked Version to produce a modified Combined Work, in the + manner specified by section 6 of the GNU GPL for conveying + Corresponding Source. + + 1) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (a) uses at run time + a copy of the Library already present on the user's computer + system, and (b) will operate properly with a modified version + of the Library that is interface-compatible with the Linked + Version. + + e) Provide Installation Information, but only if you would otherwise + be required to provide such information under section 6 of the + GNU GPL, and only to the extent that such information is + necessary to install and execute a modified version of the + Combined Work produced by recombining or relinking the + Application with a modified version of the Linked Version. (If + you use option 4d0, the Installation Information must accompany + the Minimal Corresponding Source and Corresponding Application + Code. If you use option 4d1, you must provide the Installation + Information in the manner specified by section 6 of the GNU GPL + for conveying Corresponding Source.) + + 5. Combined Libraries. + + You may place library facilities that are a work based on the +Library side by side in a single library together with other library +facilities that are not Applications and are not covered by this +License, and convey such a combined library under terms of your +choice, if you do both of the following: + + a) Accompany the combined library with a copy of the same work based + on the Library, uncombined with any other library facilities, + conveyed under the terms of this License. + + b) Give prominent notice with the combined library that part of it + is a work based on the Library, and explaining where to find the + accompanying uncombined form of the same work. + + 6. Revised Versions of the GNU Lesser General Public License. + + The Free Software Foundation may publish revised and/or new versions +of the GNU Lesser General Public License from time to time. Such new +versions will be similar in spirit to the present version, but may +differ in detail to address new problems or concerns. + + Each version is given a distinguishing version number. If the +Library as you received it specifies that a certain numbered version +of the GNU Lesser General Public License "or any later version" +applies to it, you have the option of following the terms and +conditions either of that published version or of any later version +published by the Free Software Foundation. If the Library as you +received it does not specify a version number of the GNU Lesser +General Public License, you may choose any version of the GNU Lesser +General Public License ever published by the Free Software Foundation. + + If the Library as you received it specifies that a proxy can decide +whether future versions of the GNU Lesser General Public License shall +apply, that proxy's public statement of acceptance of any version is +permanent authorization for you to choose that version for the +Library. Index: thirdparty/curveplot/README =================================================================== --- thirdparty/curveplot/README (revision 0) +++ thirdparty/curveplot/README (revision 0) @@ -0,0 +1,6 @@ +Curveplot is a library for high precision plotting of piecewise cubic +curves in OpenGL. It is distributed under the terms of the GNU Lesser +Public License. + +Curveplot is Copyright (C) 2009 Chris Laurel + Index: thirdparty/curveplot/src/curveplot.cpp =================================================================== --- thirdparty/curveplot/src/curveplot.cpp (revision 0) +++ thirdparty/curveplot/src/curveplot.cpp (revision 0) @@ -0,0 +1,807 @@ +// curveplot.cpp +// +// Copyright (C) 2009 Chris Laurel . +// +// curveplot is a module for rendering curves in OpenGL at high precision. A +// plot is a series of cubic curves. The curves are transformed +// to camera space in software because double precision is absolutely +// required. The cubics are adaptively subdivided based on distance from +// the camera position. +// +// curveplot is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// curveplot is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// CurvePlot. If not, see . + +#define DEBUG_ADAPTIVE_SPLINE 0 +#if DEBUG_ADAPTIVE_SPLINE +#define USE_VERTEX_BUFFER 0 +#else +#define USE_VERTEX_BUFFER 0 +#endif + +#include "curveplot.h" +#include "GL/glew.h" +#include +#include + +using namespace std; +using namespace Eigen; + +static const unsigned int SubdivisionFactor = 8; +static const double InvSubdivisionFactor = 1.0 / (double) SubdivisionFactor; + + + +#if DEBUG_ADAPTIVE_SPLINE +static float SplineColors[10][3] = { + { 0, 0, 1 }, + { 0, 1, 1 }, + { 0, 1, 0 }, + { 1, 1, 0 }, + { 1, 0, 0 }, + { 1, 0, 1 }, + { 0.5f, 0.5f, 1.0f }, + { 0.5f, 1.0f, 1.0f }, + { 0.5f, 1.0f, 0.5f }, + { 1.0f, 1.0f, 0.5f }, +}; + +static unsigned int SegmentCounts[32]; +#endif + +#ifndef EIGEN_VECTORIZE +// Vectorization should be enabled for improved performance. +#endif + +// Convert a 3-vector to a 4-vector by adding a zero +static inline Vector4d zeroExtend(const Vector3d& v) +{ + return Vector4d(v.x(), v.y(), v.z(), 0.0); +} + + +class HighPrec_Frustum +{ +public: + HighPrec_Frustum(double nearZ, double farZ, const Vector3d planeNormals[]) : + m_nearZ(nearZ), + m_farZ(farZ) + { + for (unsigned int i = 0; i < 4; i++) + { + m_planeNormals[i] = zeroExtend(planeNormals[i]); + } + } + + inline bool cullSphere(const Vector3d& center, + double radius) const + { + return (center.z() - radius > m_nearZ || + center.z() + radius < m_farZ || + center.dot(m_planeNormals[0].start<3>()) < -radius || + center.dot(m_planeNormals[1].start<3>()) < -radius || + center.dot(m_planeNormals[2].start<3>()) < -radius || + center.dot(m_planeNormals[3].start<3>()) < -radius); + } + + inline bool cullSphere(const Vector4d& center, + double radius) const + { + return (center.z() - radius > m_nearZ || + center.z() + radius < m_farZ || + center.dot(m_planeNormals[0]) < -radius || + center.dot(m_planeNormals[1]) < -radius || + center.dot(m_planeNormals[2]) < -radius || + center.dot(m_planeNormals[3]) < -radius); + } + + inline double nearZ() const { return m_nearZ; } + inline double farZ() const { return m_farZ; } + +private: + double m_nearZ; + double m_farZ; + Vector4d m_planeNormals[4]; +}; + + +static inline Matrix4d cubicHermiteCoefficients(const Vector4d& p0, + const Vector4d& p1, + const Vector4d& v0, + const Vector4d& v1) +{ + Matrix4d coeff; + coeff.col(0) = p0; + coeff.col(1) = v0; + coeff.col(2) = 3.0 * (p1 - p0) - (2.0 * v0 + v1); + coeff.col(3) = 2.0 * (p0 - p1) + (v1 + v0); + + return coeff; +} + + +// Test a point to see if it lies within the frustum defined by +// planes z=nearZ, z=farZ, and the four side planes with specified +// normals. +static inline bool frustumCull(const Vector4d& curvePoint, + double curveBoundingRadius, + double nearZ, double farZ, + const Vector4d viewFrustumPlaneNormals[]) +{ + return (curvePoint.z() - curveBoundingRadius > nearZ || + curvePoint.z() + curveBoundingRadius < farZ || + curvePoint.dot(viewFrustumPlaneNormals[0]) < -curveBoundingRadius || + curvePoint.dot(viewFrustumPlaneNormals[1]) < -curveBoundingRadius || + curvePoint.dot(viewFrustumPlaneNormals[2]) < -curveBoundingRadius || + curvePoint.dot(viewFrustumPlaneNormals[3]) < -curveBoundingRadius); +} + + +class HighPrec_VertexBuffer +{ +public: + HighPrec_VertexBuffer() : + currentPosition(0), + capacity(4096), + data(NULL), + vbobj(0), + currentStripLength(0) + { + data = new Vector4f[capacity]; + } + + void setup() + { +#if USE_VERTEX_BUFFER + if (vbobj) + { + glx::glBindBufferARB(GL_ARRAY_BUFFER_ARB, vbobj); + } + + glEnableClientState(GL_VERTEX_ARRAY); + + mapBuffer(); + + Vector4f* vertexBase = vbobj ? (Vector4f*) NULL : data; + glVertexPointer(3, GL_FLOAT, sizeof(Vector4f), vertexBase); + + stripLengths.clear(); + currentStripLength = 0; + currentPosition = 0; +#endif + } + + void finish() + { + unmapBuffer(); + +#if USE_VERTEX_BUFFER + if (vbobj) + { + glDisableClientState(GL_VERTEX_ARRAY); + glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0); + } +#endif + } + + inline void vertex(const Vector3d& v) + { +#if USE_VERTEX_BUFFER + data[currentPosition++].segment<3>(0) = v.cast(); + ++currentStripLength; + if (currentPosition == capacity) + { + flush(); + + data[0].segment<3>(0) = v.cast(); + currentPosition = 1; + currentStripLength = 1; + } +#else + glVertex3dv(v.data()); +#endif + } + + inline void vertex(const Vector4d& v) + { +#if USE_VERTEX_BUFFER + data[currentPosition++] = v.cast(); + ++currentStripLength; + if (currentPosition == capacity) + { + flush(); + + data[0] = v.cast(); + currentPosition = 1; + currentStripLength = 1; + } +#else + glVertex3dv(v.data()); +#endif + } + + inline void begin() + { +#if !USE_VERTEX_BUFFER + glBegin(GL_LINE_STRIP); +#endif + } + + inline void end() + { +#if USE_VERTEX_BUFFER + stripLengths.push_back(currentStripLength); + currentStripLength = 0; +#else + glEnd(); +#endif + } + + inline void flush() + { +#if USE_VERTEX_BUFFER + if (currentPosition > 0) + { + unmapBuffer(); + + // Finish the current line strip + if (currentStripLength > 1) + end(); + + unsigned int startIndex = 0; + for (vector::const_iterator iter = stripLengths.begin(); iter != stripLengths.end(); ++iter) + { + glDrawArrays(GL_LINE_STRIP, startIndex, *iter); + startIndex += *iter; + } + + mapBuffer(); + + currentPosition = 0; + stripLengths.clear(); + } + + currentStripLength = 0; +#endif + } + + void createVertexBuffer() + { +#if USE_VERTEX_BUFFER + if (!vbobj) + { + glx::glGenBuffersARB(1, &vbobj); + glx::glBindBufferARB(GL_ARRAY_BUFFER_ARB, vbobj); + glx::glBufferDataARB(GL_ARRAY_BUFFER_ARB, + capacity * sizeof(Vector4f), + NULL, + GL_STREAM_DRAW_ARB); + } +#endif + } + + void mapBuffer() + { + if (vbobj) + { + // Calling glBufferDataARB() with NULL before mapping the buffer + // is a hint to OpenGL that previous contents of vertex buffer will + // be discarded and overwritten. It enables renaming in the driver, + // hopefully resulting in performance gains. + glBufferDataARB(GL_ARRAY_BUFFER_ARB, + capacity * sizeof(Vector4f), + NULL, + GL_STREAM_DRAW_ARB); + + data = reinterpret_cast(glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY_ARB)); + } + } + + void unmapBuffer() + { +#if USE_VERTEX_BUFFER + if (vbobj) + { + glUnmapBufferARB(GL_ARRAY_BUFFER_ARB); + data = NULL; + } +#endif + } + +private: + unsigned int currentPosition; + unsigned int capacity; + Vector4f* data; + GLuint vbobj; + unsigned int currentStripLength; + vector stripLengths; +}; + + +class HighPrec_RenderContext +{ +public: + HighPrec_RenderContext(HighPrec_VertexBuffer& vbuf, + HighPrec_Frustum& viewFrustum, + double subdivisionThreshold) : + m_vbuf(vbuf), + m_viewFrustum(viewFrustum), + m_subdivisionThreshold(subdivisionThreshold) + { + } + + ~HighPrec_RenderContext() + { + /* + vbuf.flush(); + vbuf.finish(); + */ + } + + // Return the GL restart status: true if the last segment of the + // curve was culled and we need to start a new primitive sequence + // with glBegin(). + bool renderCubic(bool restartCurve, + const Matrix4d& coeff, + double t0, double t1, + double curveBoundingRadius, + int depth) const + { + const double dt = (t1 - t0) * InvSubdivisionFactor; + double segmentBoundingRadius = curveBoundingRadius * InvSubdivisionFactor; + +#if DEBUG_ADAPTIVE_SPLINE + { + int c = depth % 10; + glColor4f(SplineColors[c][0], SplineColors[c][1], SplineColors[c][2], 1.0f); + ++SegmentCounts[depth]; + } +#endif + + Vector4d lastP = coeff * Vector4d(1.0, t0, t0 * t0, t0 * t0 * t0); + + for (unsigned int i = 1; i <= SubdivisionFactor; i++) + { + double t = t0 + dt * i; + Vector4d p = coeff * Vector4d(1.0, t, t * t, t * t * t); + + double minDistance = max(-m_viewFrustum.nearZ(), abs(p.z()) - segmentBoundingRadius); + + if (segmentBoundingRadius >= m_subdivisionThreshold * minDistance) + { + if (m_viewFrustum.cullSphere(p, segmentBoundingRadius)) + { + if (!restartCurve) + { + m_vbuf.end(); + restartCurve = true; + } + } + else + { + restartCurve = renderCubic(restartCurve, + coeff, t - dt, t, + segmentBoundingRadius, + depth + 1); + } + } + else + { +#if DEBUG_ADAPTIVE_SPLINE + { + int c = depth % 10; + glColor4f(SplineColors[c][0], SplineColors[c][1], SplineColors[c][2], i % 2 ? 0.25f : 1.0f); + } +#endif + + if (restartCurve) + { + m_vbuf.begin(); + m_vbuf.vertex(lastP); + restartCurve = false; + } + m_vbuf.vertex(p); + } + lastP = p; + } + + return restartCurve; + } + +private: + HighPrec_VertexBuffer& m_vbuf; + HighPrec_Frustum& m_viewFrustum; + double m_subdivisionThreshold; +}; + + + +static HighPrec_VertexBuffer vbuf; + + +CurvePlot::CurvePlot() +{ +} + + +/** Add a new sample to the path. If the sample time is less than the first time, + * it is added at the end. If it is greater than the last time, it is appended + * to the path. The sample is ignored if it has a time in between the first and + * last times of the path. + */ +void +CurvePlot::addSample(const CurvePlotSample& sample) +{ + bool addToBack = false; + + if (m_samples.empty() || sample.t > m_samples.back().t) + { + addToBack = true; + } + else if (sample.t < m_samples.front().t) + { + addToBack = false; + } + else + { + // Sample falls within range of current samples; discard it + return; + } + + if (addToBack) + m_samples.push_back(sample); + else + m_samples.push_front(sample); + + if (m_samples.size() > 1) + { + // Calculate a bounding radius for this segment. No point on the curve will + // be further from the start point than the bounding radius. + if (addToBack) + { + const CurvePlotSample& lastSample = m_samples[m_samples.size() - 2]; + double dt = sample.t - lastSample.t; + Matrix4d coeff = cubicHermiteCoefficients(zeroExtend(lastSample.position), + zeroExtend(sample.position), + zeroExtend(lastSample.velocity * dt), + zeroExtend(sample.velocity * dt)); + Vector4d extents = coeff.cwise().abs() * Vector4d(0.0, 1.0, 1.0, 1.0); + m_samples[m_samples.size() - 1].boundingRadius = extents.norm(); + } + else + { + const CurvePlotSample& nextSample = m_samples[1]; + double dt = nextSample.t - sample.t; + Matrix4d coeff = cubicHermiteCoefficients(zeroExtend(sample.position), + zeroExtend(nextSample.position), + zeroExtend(sample.velocity * dt), + zeroExtend(nextSample.velocity * dt)); + Vector4d extents = coeff.cwise().abs() * Vector4d(0.0, 1.0, 1.0, 1.0); + m_samples[1].boundingRadius = extents.norm(); + } + } +} + + +void +CurvePlot::removeSamplesBefore(double t) +{ + while (!m_samples.empty() && m_samples.front().t < t) + { + m_samples.pop_front(); + } +} + + +void +CurvePlot::removeSamplesAfter(double t) +{ + while (!m_samples.empty() && m_samples.back().t > t) + { + m_samples.pop_back(); + } +} + + +void +CurvePlot::setDuration(double duration) +{ + m_duration = duration; +} + + +// Trajectory consists of segments, each of which is a cubic +// polynomial. + +void +CurvePlot::render(const Transform3d& modelview, + double nearZ, + double farZ, + const Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold) const +{ + // Flag to indicate whether we need to issue a glBegin() + bool restartCurve = true; + + const Vector3d& p0_ = m_samples[0].position; + const Vector3d& v0_ = m_samples[0].velocity; + Vector4d p0 = modelview * Vector4d(p0_.x(), p0_.y(), p0_.z(), 1.0); + Vector4d v0 = modelview * Vector4d(v0_.x(), v0_.y(), v0_.z(), 0.0); + + HighPrec_Frustum viewFrustum(nearZ, farZ, viewFrustumPlaneNormals); + HighPrec_RenderContext rc(vbuf, viewFrustum, subdivisionThreshold); + +#if DEBUG_ADAPTIVE_SPLINE + for (unsigned int i = 0; i < sizeof(SegmentCounts) / sizeof(SegmentCounts[0]); i++) + SegmentCounts[i] = 0; +#endif + + vbuf.createVertexBuffer(); + vbuf.setup(); + + clog << "size: " << m_samples.size() << endl; + for (unsigned int i = 1; i < m_samples.size(); i++) + { + // Transform the points into camera space. + const Vector3d& p1_ = m_samples[i].position; + const Vector3d& v1_ = m_samples[i].velocity; + Vector4d p1 = modelview * Vector4d(p1_.x(), p1_.y(), p1_.z(), 1.0); + Vector4d v1 = modelview * Vector4d(v1_.x(), v1_.y(), v1_.z(), 0.0); + + // O(t) is an approximating function for this segment of + // the orbit, with 0 <= t <= 1 + // C is the viewer position + // d(t) = |O(t) - C|, the distance from viewer to the + // orbit segment. + + double curveBoundingRadius = m_samples[i].boundingRadius; + + // Estimate the minimum possible distance from the + // curve to the z=0 plane. If the curve is far enough + // away to be approximated as a straight line, we'll just + // render it. Otherwise, it should be a performance win + // to do a sphere-frustum cull test before subdividing and + // rendering segment. + double minDistance = abs(p0.z()) - curveBoundingRadius; + + // Render close segments as splines with adaptive subdivision. The + // subdivisions eliminates kinks between line segments and also + // prevents clipping precision problems that occur when a + // very long line is rendered with a relatively small view + // volume. + if (curveBoundingRadius >= subdivisionThreshold * minDistance) + { +#if DEBUG_ADAPTIVE_SPLINE + ++SegmentCounts[0]; +#endif + // Skip rendering this section if it lies outside the view + // frustum. + if (viewFrustum.cullSphere(p0, curveBoundingRadius)) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + double dt = m_samples[i].t - m_samples[i - 1].t; + Matrix4d coeff = cubicHermiteCoefficients(p0, p1, v0 * dt, v1 * dt); + + restartCurve = rc.renderCubic(restartCurve, coeff, 0.0, 1.0, curveBoundingRadius, 1); + } + } + else + { +#if DEBUG_ADAPTIVE_SPLINE + glColor4f(SplineColors[0][0], SplineColors[0][1], SplineColors[0][2], 1.0f); +#endif + + // Apparent size of curve is small enough that we can approximate + // it as a line. + + // Simple cull test--just check the far plane + if (p0.z() + curveBoundingRadius < farZ) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + if (restartCurve) + { + vbuf.begin(); + vbuf.vertex(p0); + restartCurve = false; + } + vbuf.vertex(p1); + } + } + + p0 = p1; + v0 = v1; + } + + if (!restartCurve) + { + vbuf.end(); + } + + vbuf.flush(); + vbuf.finish(); + +#if DEBUG_ADAPTIVE_SPLINE3 + for (unsigned int i = 0; SegmentCounts[i] != 0 || i < 3; i++) + { + clog << i << ":" << SegmentCounts[i] << ", "; + } + clog << endl; +#endif +} + + +void +CurvePlot::render(const Transform3d& modelview, + double nearZ, + double farZ, + const Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold, + double startTime, + double endTime) const +{ + // Flag to indicate whether we need to issue a glBegin() + bool restartCurve = true; + + if (m_samples.empty() || endTime <= m_samples.front().t || startTime >= m_samples.back().t) + return; + + // Linear search for the first sample + unsigned int startSample = 0; + while (startSample < m_samples.size() - 1 && startTime < m_samples[startSample].t) + startSample++; + + // Start at the first sample with time <= startTime + if (startSample > 0) + startSample--; + + const Vector3d& p0_ = m_samples[startSample].position; + const Vector3d& v0_ = m_samples[startSample].velocity; + Vector4d p0 = modelview * Vector4d(p0_.x(), p0_.y(), p0_.z(), 1.0); + Vector4d v0 = modelview * Vector4d(v0_.x(), v0_.y(), v0_.z(), 0.0); + + HighPrec_Frustum viewFrustum(nearZ, farZ, viewFrustumPlaneNormals); + HighPrec_RenderContext rc(vbuf, viewFrustum, subdivisionThreshold); + + vbuf.createVertexBuffer(); + vbuf.setup(); + + bool firstSegment = true; + bool lastSegment = false; + + for (unsigned int i = startSample + 1; i < m_samples.size() && !lastSegment; i++) + { + // Transform the points into camera space. + const Vector3d& p1_ = m_samples[i].position; + const Vector3d& v1_ = m_samples[i].velocity; + Vector4d p1 = modelview * Vector4d(p1_.x(), p1_.y(), p1_.z(), 1.0); + Vector4d v1 = modelview * Vector4d(v1_.x(), v1_.y(), v1_.z(), 0.0); + + if (endTime <= m_samples[i].t) + { + lastSegment = true; + } + + // O(t) is an approximating function for this segment of + // the orbit, with 0 <= t <= 1 + // C is the viewer position + // d(t) = |O(t) - C|, the distance from viewer to the + // orbit segment. + + double curveBoundingRadius = m_samples[i].boundingRadius; + + // Estimate the minimum possible distance from the + // curve to the z=0 plane. If the curve is far enough + // away to be approximated as a straight line, we'll just + // render it. Otherwise, it should be a performance win + // to do a sphere-frustum cull test before subdividing and + // rendering segment. + double minDistance = abs(p0.z()) - curveBoundingRadius; + + // Render close segments as splines with adaptive subdivision. The + // subdivisions eliminates kinks between line segments and also + // prevents clipping precision problems that occur when a + // very long line is rendered with a relatively small view + // volume. + if (curveBoundingRadius >= subdivisionThreshold * minDistance || lastSegment || firstSegment) + { + // Skip rendering this section if it lies outside the view + // frustum. + if (viewFrustum.cullSphere(p0, curveBoundingRadius)) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + double dt = m_samples[i].t - m_samples[i - 1].t; + double t0 = 0.0; + double t1 = 1.0; + + if (firstSegment) + { + t0 = (startTime - m_samples[i - 1].t) / dt; + t0 = std::max(0.0, std::min(1.0, t0)); + firstSegment = false; + } + + if (lastSegment) + { + t1 = (endTime - m_samples[i - 1].t) / dt; + } + + Matrix4d coeff = cubicHermiteCoefficients(p0, p1, v0 * dt, v1 * dt); + restartCurve = rc.renderCubic(restartCurve, coeff, t0, t1, curveBoundingRadius, 1); + } + } + else + { + // Apparent size of curve is small enough that we can approximate + // it as a line. + + // Simple cull test--just check the far plane. This is required because + // apparent clipping precision limitations can cause a GPU to draw lines + // that lie completely beyond the far plane. + if (p0.z() + curveBoundingRadius < farZ) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + if (restartCurve) + { + vbuf.begin(); + vbuf.vertex(p0); + restartCurve = false; + } + vbuf.vertex(p1); + } + } + + p0 = p1; + v0 = v1; + } + + if (!restartCurve) + { + vbuf.end(); + } + + vbuf.flush(); + vbuf.finish(); +}