openscenegraph
CoordinateSystemNode
Go to the documentation of this file.
1/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2006 Robert Osfield
2 *
3 * This library is open source and may be redistributed and/or modified under
4 * the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
5 * (at your option) any later version. The full license is in LICENSE file
6 * included with this distribution, and on the openscenegraph.org website.
7 *
8 * This library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * OpenSceneGraph Public License for more details.
12*/
13
14#ifndef OSG_COORDINATESYSTEMNODE
15#define OSG_COORDINATESYSTEMNODE 1
16
17#include <osg/Group>
18#include <osg/Matrixd>
19
20namespace osg
21{
22
23const double WGS_84_RADIUS_EQUATOR = 6378137.0;
24const double WGS_84_RADIUS_POLAR = 6356752.3142;
25
26/** EllipsoidModel encapsulates the ellipsoid used to model astronomical bodies,
27 * such as sun, planets, moon etc.
28 * All distance quantities (i.e. heights + radius) are in meters,
29 * and latitude and longitude are in radians.*/
30class EllipsoidModel : public Object
31{
32 public:
33
34 /** WGS_84 is a common representation of the earth's spheroid */
35 EllipsoidModel(double radiusEquator = WGS_84_RADIUS_EQUATOR,
36 double radiusPolar = WGS_84_RADIUS_POLAR):
37 _radiusEquator(radiusEquator),
38 _radiusPolar(radiusPolar) { computeCoefficients(); }
39
40 EllipsoidModel(const EllipsoidModel& et,const CopyOp& copyop=CopyOp::SHALLOW_COPY):
41 Object(et,copyop),
42 _radiusEquator(et._radiusEquator),
43 _radiusPolar(et._radiusPolar) { computeCoefficients(); }
44
45 META_Object(osg,EllipsoidModel);
46
47 void setRadiusEquator(double radius) { _radiusEquator = radius; computeCoefficients(); }
48 double getRadiusEquator() const { return _radiusEquator; }
49
50 void setRadiusPolar(double radius) { _radiusPolar = radius; computeCoefficients(); }
51 double getRadiusPolar() const { return _radiusPolar; }
52
53 inline void convertLatLongHeightToXYZ(double latitude, double longitude, double height,
54 double& X, double& Y, double& Z) const;
55
56 inline void convertXYZToLatLongHeight(double X, double Y, double Z,
57 double& latitude, double& longitude, double& height) const;
58
59 inline void computeLocalToWorldTransformFromLatLongHeight(double latitude, double longitude, double height, osg::Matrixd& localToWorld) const;
60
61 inline void computeLocalToWorldTransformFromXYZ(double X, double Y, double Z, osg::Matrixd& localToWorld) const;
62
63 inline void computeCoordinateFrame(double latitude, double longitude, osg::Matrixd& localToWorld) const;
64
65 inline osg::Vec3d computeLocalUpVector(double X, double Y, double Z) const;
66
67 // Convenience method for determining if EllipsoidModel is a stock WGS84 ellipsoid
68 inline bool isWGS84() const {return(_radiusEquator == WGS_84_RADIUS_EQUATOR && _radiusPolar == WGS_84_RADIUS_POLAR);}
69
70 // Compares two EllipsoidModel by comparing critical internal values.
71 // Ignores _eccentricitySquared since it's just a cached value derived from
72 // the _radiusEquator and _radiusPolar members.
73 friend bool operator == ( const EllipsoidModel & e1, const EllipsoidModel& e2) {return(e1._radiusEquator == e2._radiusEquator && e1._radiusPolar == e2._radiusPolar);}
74
75
76 protected:
77
78 void computeCoefficients()
79 {
80 double flattening = (_radiusEquator-_radiusPolar)/_radiusEquator;
81 _eccentricitySquared = 2*flattening - flattening*flattening;
82 }
83
84 double _radiusEquator;
85 double _radiusPolar;
86 double _eccentricitySquared;
87
88};
89
90/** CoordinateFrame encapsulates the orientation of east, north and up.*/
91typedef Matrixd CoordinateFrame;
92
93/** CoordinateSystem encapsulate the coordinate system that is associated with objects in a scene.
94 For an overview of common earth bases coordinate systems see http://www.colorado.edu/geography/gcraft/notes/coordsys/coordsys_f.html */
95class OSG_EXPORT CoordinateSystemNode : public Group
96{
97 public:
98
99 CoordinateSystemNode();
100
101 CoordinateSystemNode(const std::string& format, const std::string& cs);
102
103 /** Copy constructor using CopyOp to manage deep vs shallow copy.*/
104 CoordinateSystemNode(const CoordinateSystemNode&,const osg::CopyOp& copyop=osg::CopyOp::SHALLOW_COPY);
105
106 META_Node(osg,CoordinateSystemNode);
107
108
109 /** Set the coordinate system node up by copying the format, coordinate system string, and ellipsoid model of another coordinate system node.*/
110 void set(const CoordinateSystemNode& csn);
111
112 /** Set the coordinate system format string. Typical values would be WKT, PROJ4, USGS etc.*/
113 void setFormat(const std::string& format) { _format = format; }
114
115 /** Get the coordinate system format string.*/
116 const std::string& getFormat() const { return _format; }
117
118 /** Set the CoordinateSystem reference string, should be stored in a form consistent with the Format.*/
119 void setCoordinateSystem(const std::string& cs) { _cs = cs; }
120
121 /** Get the CoordinateSystem reference string.*/
122 const std::string& getCoordinateSystem() const { return _cs; }
123
124
125 /** Set EllipsoidModel to describe the model used to map lat, long and height into geocentric XYZ and back. */
126 void setEllipsoidModel(EllipsoidModel* ellipsode) { _ellipsoidModel = ellipsode; }
127
128 /** Get the EllipsoidModel.*/
129 EllipsoidModel* getEllipsoidModel() { return _ellipsoidModel.get(); }
130
131 /** Get the const EllipsoidModel.*/
132 const EllipsoidModel* getEllipsoidModel() const { return _ellipsoidModel.get(); }
133
134 /** Compute the local coordinate frame for specified point.*/
135 CoordinateFrame computeLocalCoordinateFrame(const Vec3d& position) const;
136
137 /** Compute the local up-vector for specified point.*/
138 osg::Vec3d computeLocalUpVector(const Vec3d& position) const;
139
140 protected:
141
142 virtual ~CoordinateSystemNode() {}
143
144 std::string _format;
145 std::string _cs;
146 ref_ptr<EllipsoidModel> _ellipsoidModel;
147
148};
149
150
151
152////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
153// implement inline methods.
154//
155inline void EllipsoidModel::convertLatLongHeightToXYZ(double latitude, double longitude, double height,
156 double& X, double& Y, double& Z) const
157{
158 // for details on maths see http://www.colorado.edu/geography/gcraft/notes/datum/gif/llhxyz.gif
159 double sin_latitude = sin(latitude);
160 double cos_latitude = cos(latitude);
161 double N = _radiusEquator / sqrt( 1.0 - _eccentricitySquared*sin_latitude*sin_latitude);
162 X = (N+height)*cos_latitude*cos(longitude);
163 Y = (N+height)*cos_latitude*sin(longitude);
164 Z = (N*(1-_eccentricitySquared)+height)*sin_latitude;
165}
166
167
168inline void EllipsoidModel::convertXYZToLatLongHeight(double X, double Y, double Z,
169 double& latitude, double& longitude, double& height) const
170{
171 // handle polar and center-of-earth cases directly.
172 if (X != 0.0)
173 longitude = atan2(Y,X);
174 else
175 {
176 if (Y > 0.0)
177 longitude = PI_2;
178 else if (Y < 0.0)
179 longitude = -PI_2;
180 else
181 {
182 // at pole or at center of the earth
183 longitude = 0.0;
184 if (Z > 0.0)
185 { // north pole.
186 latitude = PI_2;
187 height = Z - _radiusPolar;
188 }
189 else if (Z < 0.0)
190 { // south pole.
191 latitude = -PI_2;
192 height = -Z - _radiusPolar;
193 }
194 else
195 { // center of earth.
196 latitude = PI_2;
197 height = -_radiusPolar;
198 }
199 return;
200 }
201 }
202
203 // http://www.colorado.edu/geography/gcraft/notes/datum/gif/xyzllh.gif
204 double p = sqrt(X*X + Y*Y);
205 double theta = atan2(Z*_radiusEquator , (p*_radiusPolar));
206 double eDashSquared = (_radiusEquator*_radiusEquator - _radiusPolar*_radiusPolar)/
207 (_radiusPolar*_radiusPolar);
208
209 double sin_theta = sin(theta);
210 double cos_theta = cos(theta);
211
212 latitude = atan( (Z + eDashSquared*_radiusPolar*sin_theta*sin_theta*sin_theta) /
213 (p - _eccentricitySquared*_radiusEquator*cos_theta*cos_theta*cos_theta) );
214
215 double sin_latitude = sin(latitude);
216 double N = _radiusEquator / sqrt( 1.0 - _eccentricitySquared*sin_latitude*sin_latitude);
217
218 height = p/cos(latitude) - N;
219}
220
221inline void EllipsoidModel::computeLocalToWorldTransformFromLatLongHeight(double latitude, double longitude, double height, osg::Matrixd& localToWorld) const
222{
223 double X, Y, Z;
224 convertLatLongHeightToXYZ(latitude,longitude,height,X,Y,Z);
225
226 localToWorld.makeTranslate(X,Y,Z);
227 computeCoordinateFrame(latitude, longitude, localToWorld);
228}
229
230inline void EllipsoidModel::computeLocalToWorldTransformFromXYZ(double X, double Y, double Z, osg::Matrixd& localToWorld) const
231{
232 double latitude, longitude, height;
233 convertXYZToLatLongHeight(X,Y,Z,latitude,longitude,height);
234
235 localToWorld.makeTranslate(X,Y,Z);
236 computeCoordinateFrame(latitude, longitude, localToWorld);
237}
238
239inline void EllipsoidModel::computeCoordinateFrame(double latitude, double longitude, osg::Matrixd& localToWorld) const
240{
241 // Compute up vector
242 osg::Vec3d up ( cos(longitude)*cos(latitude), sin(longitude)*cos(latitude), sin(latitude));
243
244 // Compute east vector
245 osg::Vec3d east (-sin(longitude), cos(longitude), 0);
246
247 // Compute north vector = outer product up x east
248 osg::Vec3d north = up ^ east;
249
250 // set matrix
251 localToWorld(0,0) = east[0];
252 localToWorld(0,1) = east[1];
253 localToWorld(0,2) = east[2];
254
255 localToWorld(1,0) = north[0];
256 localToWorld(1,1) = north[1];
257 localToWorld(1,2) = north[2];
258
259 localToWorld(2,0) = up[0];
260 localToWorld(2,1) = up[1];
261 localToWorld(2,2) = up[2];
262}
263
264inline osg::Vec3d EllipsoidModel::computeLocalUpVector(double X, double Y, double Z) const
265{
266 // Note latitude is angle between normal to ellipsoid surface and XY-plane
267 double latitude;
268 double longitude;
269 double altitude;
270 convertXYZToLatLongHeight(X,Y,Z,latitude,longitude,altitude);
271
272 // Compute up vector
273 return osg::Vec3d( cos(longitude) * cos(latitude),
274 sin(longitude) * cos(latitude),
275 sin(latitude));
276}
277
278}
279#endif