StelSphereGeometry.cpp   StelSphereGeometry.cpp 
skipping to change at line 20 skipping to change at line 20
* This program is distributed in the hope that it will be useful, * This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of * but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details. * GNU General Public License for more details.
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software * along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, U SA. * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, U SA.
*/ */
#include "StelSphereGeometry.hpp"
#include <QDebug> #include <QDebug>
#include <QBuffer>
#include <stdexcept>
#include "StelSphereGeometry.hpp"
#include "StelUtils.hpp"
#include "StelJsonParser.hpp"
using namespace StelGeom; // Definition of static constants.
const QVariant::Type SphericalRegionP::qVariantType = (QVariant::Type)(QVar
iant::UserType+1);
int SphericalRegionP::metaTypeId = SphericalRegionP::initialize();
//! Special constructor for 3 halfspaces convex int SphericalRegionP::initialize()
ConvexS::ConvexS(const Vec3d &e0,const Vec3d &e1,const Vec3d &e2)
{ {
reserve(3); int id = qRegisterMetaType<SphericalRegionP>();
push_back(e1^e0); qRegisterMetaTypeStreamOperators<SphericalRegionP>("SphericalRegionP
push_back(e2^e1); ");
push_back(e0^e2); return id;
}
// Warning: vectors not normalized while they should be QDataStream& operator<<(QDataStream& out, const SphericalRegionP& region)
// In this case it works because d==0 for each HalfSpace {
out << (quint8)region->getType();
region->serialize(out);
return out;
} }
//! Special constructor for 4 halfspaces convex QDataStream& operator>>(QDataStream& in, SphericalRegionP& region)
ConvexS::ConvexS(const Vec3d &e0,const Vec3d &e1,const Vec3d &e2, const Vec
3d &e3)
{ {
reserve(4); quint8 regType;
const double d = e3*((e2-e3)^(e1-e3)); in >> regType;
if (d > 0) switch (regType)
{ {
push_back(e1^e0); case SphericalRegion::Empty:
push_back(e2^e1); region = EmptySphericalRegion::staticInstance;
push_back(e3^e2); return in;
push_back(e0^e3); case SphericalRegion::AllSky:
region = AllSkySphericalRegion::staticInstance;
return in;
case SphericalRegion::Cap:
region = SphericalCap::deserialize(in);
return in;
case SphericalRegion::ConvexPolygon:
region = SphericalConvexPolygon::deserialize(in);
return in;
case SphericalRegion::Polygon:
region = SphericalPolygon::deserialize(in);
return in;
case SphericalRegion::Point:
region = SphericalPoint::deserialize(in);
return in;
default:
Q_ASSERT(0); // Unknown region type
}
return in;
}
///////////////////////////////////////////////////////////////////////////
////////////////////
// Default implementations of methods for SphericalRegion
///////////////////////////////////////////////////////////////////////////
////////////////////
QByteArray SphericalRegion::toJSON() const
{
QByteArray res;
QBuffer buf1(&res);
buf1.open(QIODevice::WriteOnly);
StelJsonParser::write(toQVariant(), &buf1);
buf1.close();
return res;
}
SphericalRegionP SphericalRegion::getEnlarged(double margin) const
{
Q_ASSERT(margin>=0);
if (margin>=M_PI)
return SphericalRegionP(new AllSkySphericalRegion());
const SphericalCap& cap = getBoundingCap();
double newRadius = std::acos(cap.d)+margin;
if (newRadius>=M_PI)
return SphericalRegionP(new AllSkySphericalRegion());
return SphericalRegionP(new SphericalCap(cap.n, std::cos(newRadius))
);
}
QVector<SphericalCap> SphericalRegion::getBoundingSphericalCaps() const
{
QVector<SphericalCap> res;
res << getBoundingCap();
return res;
}
// Warning: vectors not normalized while they should be // Default slow implementation o(n^2).
// In this case it works because d==0 for each HalfSpace SphericalCap SphericalRegion::getBoundingCap() const
{
SphericalCap res;
getOctahedronPolygon().getBoundingCap(res.n, res.d);
return res;
}
bool SphericalRegion::contains(const SphericalPolygon& r) const {return con
tainsDefault(&r);}
bool SphericalRegion::contains(const SphericalConvexPolygon& r) const {retu
rn containsDefault(&r);}
bool SphericalRegion::contains(const SphericalCap& r) const {return contain
sDefault(&r);}
bool SphericalRegion::contains(const SphericalPoint& r) const {return conta
ins(r.n);}
bool SphericalRegion::contains(const AllSkySphericalRegion& r) const {retur
n containsDefault(&r);}
bool SphericalRegion::contains(const SphericalRegion* r) const
{
switch (r->getType())
{
case SphericalRegion::Point:
return contains(static_cast<const SphericalPoint*>(r
)->n);
case SphericalRegion::Cap:
return contains(*static_cast<const SphericalCap*>(r)
);
case SphericalRegion::Polygon:
return contains(*static_cast<const SphericalPolygon*
>(r));
case SphericalRegion::ConvexPolygon:
return contains(*static_cast<const SphericalConvexPo
lygon*>(r));
case SphericalRegion::AllSky:
return contains(*static_cast<const AllSkySphericalRe
gion*>(r));
case SphericalRegion::Empty:
return false;
default:
return containsDefault(r);
} }
else Q_ASSERT(0);
return false;
}
bool SphericalRegion::intersects(const SphericalPolygon& r) const {return i
ntersectsDefault(&r);}
bool SphericalRegion::intersects(const SphericalConvexPolygon& r) const {re
turn intersectsDefault(&r);}
bool SphericalRegion::intersects(const SphericalCap& r) const {return inter
sectsDefault(&r);}
bool SphericalRegion::intersects(const SphericalPoint& r) const {return con
tains(r.n);}
bool SphericalRegion::intersects(const AllSkySphericalRegion& r) const {ret
urn getType()==SphericalRegion::Empty ? false : true;}
bool SphericalRegion::intersects(const SphericalRegion* r) const
{
switch (r->getType())
{ {
push_back((e2-e3)^(e1-e3)); case SphericalRegion::Point:
(*begin()).d = d; return intersects(static_cast<const SphericalPoint*>
(*begin()).n.normalize(); (r)->n);
case SphericalRegion::Cap:
return intersects(*static_cast<const SphericalCap*>(
r));
case SphericalRegion::Polygon:
return intersects(*static_cast<const SphericalPolygo
n*>(r));
case SphericalRegion::ConvexPolygon:
return intersects(*static_cast<const SphericalConvex
Polygon*>(r));
case SphericalRegion::AllSky:
return intersects(*static_cast<const AllSkySpherical
Region*>(r));
case SphericalRegion::Empty:
return false;
default:
return intersectsDefault(r);
} }
Q_ASSERT(0);
return false;
} }
//! Special constructor for 3 points polygon SphericalRegionP SphericalRegion::getIntersection(const SphericalRegion* r)
Polygon::Polygon(const Vec3d &e0,const Vec3d &e1,const Vec3d &e2) const
{ {
reserve(3); switch (r->getType())
push_back(e0); {
push_back(e1); case SphericalRegion::Point:
push_back(e2); return getIntersection(static_cast<const SphericalPo
int*>(r)->n);
case SphericalRegion::Cap:
return getIntersection(*static_cast<const SphericalC
ap*>(r));
case SphericalRegion::Polygon:
return getIntersection(*static_cast<const SphericalP
olygon*>(r));
case SphericalRegion::ConvexPolygon:
return getIntersection(*static_cast<const SphericalC
onvexPolygon*>(r));
case SphericalRegion::AllSky:
return getIntersection(*static_cast<const AllSkySphe
ricalRegion*>(r));
case SphericalRegion::Empty:
return false;
default:
return getIntersectionDefault(r);
}
Q_ASSERT(0);
return SphericalRegionP();
} }
SphericalRegionP SphericalRegion::getIntersection(const SphericalPolygon& r
) const {return getIntersectionDefault(&r);}
SphericalRegionP SphericalRegion::getIntersection(const SphericalConvexPoly
gon& r) const {return getIntersectionDefault(&r);}
SphericalRegionP SphericalRegion::getIntersection(const SphericalCap& r) co
nst {return getIntersectionDefault(&r);}
SphericalRegionP SphericalRegion::getIntersection(const SphericalPoint& r)
const {return getIntersectionDefault(&r);}
SphericalRegionP SphericalRegion::getIntersection(const AllSkySphericalRegi
on& r) const {return getIntersectionDefault(&r);}
SphericalRegionP SphericalRegion::getIntersection(const EmptySphericalRegio
n& r) const {return SphericalRegionP(new EmptySphericalRegion());}
//! Special constructor for 4 points polygon SphericalRegionP SphericalRegion::getUnion(const SphericalRegion* r) const
Polygon::Polygon(const Vec3d &e0,const Vec3d &e1,const Vec3d &e2, const Vec
3d &e3)
{ {
reserve(4); switch (r->getType())
push_back(e0); {
push_back(e1); case SphericalRegion::Point:
push_back(e2); return getUnion(static_cast<const SphericalPoint*>(r
push_back(e3); )->n);
case SphericalRegion::Cap:
return getUnion(*static_cast<const SphericalCap*>(r)
);
case SphericalRegion::Polygon:
return getUnion(*static_cast<const SphericalPolygon*
>(r));
case SphericalRegion::ConvexPolygon:
return getUnion(*static_cast<const SphericalConvexPo
lygon*>(r));
case SphericalRegion::AllSky:
return getUnion(*static_cast<const AllSkySphericalRe
gion*>(r));
case SphericalRegion::Empty:
return false;
default:
return getUnionDefault(r);
}
Q_ASSERT(0);
return SphericalRegionP();
} }
SphericalRegionP SphericalRegion::getUnion(const SphericalPolygon& r) const
{return getUnionDefault(&r);}
SphericalRegionP SphericalRegion::getUnion(const SphericalConvexPolygon& r)
const {return getUnionDefault(&r);}
SphericalRegionP SphericalRegion::getUnion(const SphericalCap& r) const {re
turn getUnionDefault(&r);}
SphericalRegionP SphericalRegion::getUnion(const SphericalPoint& r) const {
return getUnionDefault(&r);}
SphericalRegionP SphericalRegion::getUnion(const AllSkySphericalRegion& r)
const {return SphericalRegionP(new AllSkySphericalRegion());}
SphericalRegionP SphericalRegion::getUnion(const EmptySphericalRegion& r) c
onst {return getUnionDefault(&r);}
ConvexPolygon ConvexPolygon::fullSky() SphericalRegionP SphericalRegion::getSubtraction(const SphericalRegion* r) const
{ {
ConvexPolygon poly; switch (r->getType())
poly.asConvex().push_back(HalfSpace(Vec3d(1.,0.,0.), -1.)); {
return poly; case SphericalRegion::Point:
return getSubtraction(static_cast<const SphericalPoi
nt*>(r)->n);
case SphericalRegion::Cap:
return getSubtraction(*static_cast<const SphericalCa
p*>(r));
case SphericalRegion::Polygon:
return getSubtraction(*static_cast<const SphericalPo
lygon*>(r));
case SphericalRegion::ConvexPolygon:
return getSubtraction(*static_cast<const SphericalCo
nvexPolygon*>(r));
case SphericalRegion::AllSky:
return getSubtraction(*static_cast<const AllSkySpher
icalRegion*>(r));
case SphericalRegion::Empty:
return false;
default:
return getSubtractionDefault(r);
}
Q_ASSERT(0);
return SphericalRegionP();
} }
SphericalRegionP SphericalRegion::getSubtraction(const SphericalPolygon& r)
const {return getSubtractionDefault(&r);}
SphericalRegionP SphericalRegion::getSubtraction(const SphericalConvexPolyg
on& r) const {return getSubtractionDefault(&r);}
SphericalRegionP SphericalRegion::getSubtraction(const SphericalCap& r) con
st {return getSubtractionDefault(&r);}
SphericalRegionP SphericalRegion::getSubtraction(const SphericalPoint& r) c
onst {return getSubtractionDefault(&r);}
SphericalRegionP SphericalRegion::getSubtraction(const AllSkySphericalRegio
n& r) const {return SphericalRegionP(new EmptySphericalRegion());}
SphericalRegionP SphericalRegion::getSubtraction(const EmptySphericalRegion
& r) const {return getSubtractionDefault(&r);}
//! Check if the polygon is valid, i.e. it has no side >180 // Returns whether another SphericalPolygon intersects with the SphericalPo
bool ConvexPolygon::checkValid() const lygon.
bool SphericalRegion::containsDefault(const SphericalRegion* r) const
{ {
const ConvexS& cvx = asConvex(); if (!getBoundingCap().contains(r->getBoundingCap()))
const Polygon& poly = asPolygon();
if (cvx.size()<3)
return false; return false;
bool res=true; return getOctahedronPolygon().contains(r->getOctahedronPolygon());
for (size_t i=0;i<cvx.size();++i) }
// Returns whether another SphericalPolygon intersects with the SphericalPo
lygon.
bool SphericalRegion::intersectsDefault(const SphericalRegion* r) const
{
if (!getBoundingCap().intersects(r->getBoundingCap()))
return false;
return getOctahedronPolygon().intersects(r->getOctahedronPolygon());
}
// Return a new SphericalPolygon consisting of the intersection of this and
the given SphericalPolygon.
SphericalRegionP SphericalRegion::getIntersectionDefault(const SphericalReg
ion* r) const
{
if (!getBoundingCap().intersects(r->getBoundingCap()))
return SphericalRegionP(new EmptySphericalRegion());
OctahedronPolygon resOct(getOctahedronPolygon());
resOct.inPlaceIntersection(r->getOctahedronPolygon());
return SphericalRegionP(new SphericalPolygon(resOct));
}
// Return a new SphericalPolygon consisting of the union of this and the gi
ven SphericalPolygon.
SphericalRegionP SphericalRegion::getUnionDefault(const SphericalRegion* r)
const
{
OctahedronPolygon resOct(getOctahedronPolygon());
resOct.inPlaceUnion(r->getOctahedronPolygon());
return SphericalRegionP(new SphericalPolygon(resOct));
}
// Return a new SphericalPolygon consisting of the subtraction of the given
SphericalPolygon from this.
SphericalRegionP SphericalRegion::getSubtractionDefault(const SphericalRegi
on* r) const
{
OctahedronPolygon resOct(getOctahedronPolygon());
resOct.inPlaceSubtraction(r->getOctahedronPolygon());
return SphericalRegionP(new SphericalPolygon(resOct));
}
///////////////////////////////////////////////////////////////////////////
/
// Methods for SphericalCap
///////////////////////////////////////////////////////////////////////////
/
// Returns whether a SphericalPolygon is contained into the region.
bool SphericalCap::contains(const SphericalConvexPolygon& cvx) const
{
foreach (const Vec3d& v, cvx.getConvexContour())
{ {
// Check that all points not on the current convex plane are if (!contains(v))
included in it return false;
for (size_t p=0;p<cvx.size()-2;++p)
res &= cvx[i].contains(poly[(p+i+2)%poly.size()]);
} }
return res; return true;
}
bool SphericalCap::containsTriangle(const Vec3d* v) const
{
return contains(*(v++)) && contains(*(v++)) && contains(*v);
} }
//! Return the convex polygon area in steradians bool SphericalCap::intersectsTriangle(const Vec3d* v) const
// TODO Optimize using add oc formulas from http://en.wikipedia.org/wiki/So
lid_angle
double ConvexPolygon::getArea() const
{ {
// Use Girard's theorem if (contains(*v) || contains(*(v+1)) || contains(*(v+2)))
double angleSum=0.; return true;
const ConvexS& cvx = asConvex(); // No points of the triangle are inside the cap
const int size = cvx.size(); if (d<=0)
return false;
if (!sideHalfSpaceIntersects(*v, *(v+1), *this) || !sideHalfSpaceInt
ersects(*(v+1), *(v+2), *this) || !sideHalfSpaceIntersects(*(v+2), *v, *thi
s))
return false;
// Warning!!!! There is a last case which is not managed!
// When all the points of the polygon are outside the circle but the
halfspace of the corner the closest to the
// circle intersects the circle halfspace..
return true;
}
if (size==1) bool SphericalCap::intersectsConvexContour(const Vec3d* vertice, int nbVert
ice) const
{
for (int i=0;i<nbVertice;++i)
{ {
// Handle special case for > 180 degree polygons if (contains(vertice[i]))
return cvx[0].getArea(); return true;
} }
// No points of the convex polygon are inside the cap
if (d<=0)
return false;
// Sum the angles at each corner of the polygon for (int i=0;i<nbVertice-1;++i)
// (the non cartesian angle is found from the plan normals)
for (int i=0;i<size-1;++i)
{ {
angleSum += M_PI-cvx[i].n.angle(cvx[i+1].n); if (!sideHalfSpaceIntersects(vertice[i], vertice[i+1], *this
))
return false;
} }
// Last angle if (!sideHalfSpaceIntersects(vertice[nbVertice-1], vertice[0], *this
angleSum += M_PI-cvx[size-1].n.angle(cvx[0].n); ))
return angleSum - M_PI*(size-2); return false;
// Warning!!!! There is a last case which is not managed!
// When all the points of the polygon are outside the circle but the
halfspace of the corner the closest to the
// circle intersects the circle halfspace..
return true;
} }
//! Return the convex polygon barycenter // Returns whether a SphericalPolygon intersects the region.
// TODO this code is quite wrong but good for triangles bool SphericalCap::intersects(const SphericalConvexPolygon& cvx) const
Vec3d ConvexPolygon::getBarycenter() const
{ {
if (ConvexS::size()==1) // TODO This algo returns sometimes false positives!!
{ return intersectsConvexContour(cvx.getConvexContour().constData(), c
// Handle special case for > 180 degree polygons vx.getConvexContour().size());
return asConvex()[0].n; }
}
Vec3d barycenter(0.); bool SphericalCap::intersects(const SphericalPolygon& polyBase) const
for (unsigned int i=0;i<Polygon::size();++i) {
// Go through the full list of triangle
const QVector<Vec3d>& vArray = polyBase.getFillVertexArray().vertex;
for (int i=0;i<vArray.size()/3;++i)
{ {
barycenter += Polygon::operator[](i); if (intersectsConvexContour(vArray.constData()+i*3, 3))
return true;
} }
barycenter.normalize(); return false;
return barycenter;
} }
namespace StelGeom bool SphericalCap::clipGreatCircle(Vec3d& v1, Vec3d& v2) const
{ {
const bool b1 = contains(v1);
const bool b2 = contains(v2);
if (b1)
{
if (b2)
{
// Both points inside, nothing to do
return true;
}
else
{
// v1 inside, v2 outside. Return segment v1 -- inter
section
Vec3d v = v1^v2;
v.normalize();
Vec3d vv;
SphericalCap::intersectionPoints(*this, SphericalCap
(v, 0), v, vv);
const float cosDist = v1*v2;
v2 = (v1*v >= cosDist && v2*v >= cosDist) ? v : vv;
return true;
}
}
else
{
if (b2)
{
// v2 inside, v1 outside. Return segment v2 -- inter
section
Vec3d v = v1^v2;
v.normalize();
Vec3d vv;
SphericalCap::intersectionPoints(*this, SphericalCap
(v, 0), v, vv);
const float cosDist = v1*v2;
v1 = (v1*v >= cosDist && v2*v >= cosDist) ? v : vv;
return true;
}
else
{
// Both points outside
Vec3d v = v1^v2;
v.normalize();
Vec3d vv;
if (!SphericalCap::intersectionPoints(*this, Spheric
alCap(v, 0), v, vv))
return false;
const float cosDist = v1*v2;
if (v1*v >= cosDist && v2*v >= cosDist && v1*vv >= c
osDist && v2*vv >= cosDist)
{
v1 = v;
v2 = vv;
return true;
}
return false;
}
}
Q_ASSERT(0);
return false;
}
//! Compute the intersection of the planes defined by the 2 halfspaces on t //! Compute the intersection of the circles defined by the 2 caps on the sp
he sphere (usually on 2 points) and return it in p1 and p2. here (usually on 2 points) and return it in p1 and p2.
//! If the 2 HalfSpaces don't interesect or intersect only at 1 point, fals //! If the 2 SphericalCaps don't interesect or intersect only at 1 point, f
e is returned and p1 and p2 are undefined alse is returned and p1 and p2 are undefined.
bool planeIntersect2(const HalfSpace& h1, const HalfSpace& h2, Vec3d& p1, V bool SphericalCap::intersectionPoints(const SphericalCap& h1, const Spheric
ec3d& p2) alCap& h2, Vec3d& p1, Vec3d& p2)
{ {
if (!h1.intersects(h2))
return false;
const Vec3d& n1 = h1.n; const Vec3d& n1 = h1.n;
const Vec3d& n2 = h2.n; const Vec3d& n2 = h2.n;
const double& d1 = -h1.d; const double& d1 = -h1.d;
const double& d2 = -h2.d; const double& d2 = -h2.d;
const double& a1 = n1[0]; const double& a1 = n1[0];
const double& b1 = n1[1]; const double& b1 = n1[1];
const double& c1 = n1[2]; const double& c1 = n1[2];
const double& a2 = n2[0]; const double& a2 = n2[0];
const double& b2 = n2[1]; const double& b2 = n2[1];
const double& c2 = n2[2]; const double& c2 = n2[2];
skipping to change at line 234 skipping to change at line 525
const double t2 = (-b-sqrtD)/2.; const double t2 = (-b-sqrtD)/2.;
p1 = p0+u*t1; p1 = p0+u*t1;
p2 = p0+u*t2; p2 = p0+u*t2;
Q_ASSERT(fabs(p1.lengthSquared()-1.)<0.000001); Q_ASSERT(fabs(p1.lengthSquared()-1.)<0.000001);
Q_ASSERT(fabs(p2.lengthSquared()-1.)<0.000001); Q_ASSERT(fabs(p2.lengthSquared()-1.)<0.000001);
return true; return true;
} }
double SphericalCap::relativeAreaOverlap(const SphericalCap& c1, const Sphe
ricalCap& c2)
{
// TODO, this is wrong for many reasons
Q_ASSERT(0);
if (!c1.intersects(c2))
return 0.;
if (c1.contains(c2))
return c2.getArea()/c1.getArea();
if (c2.contains(c1))
return c1.getArea()/c2.getArea();
Vec3d p1, p2;
double area1=c1.getArea();
double area2=c2.getArea();
#ifndef NDEBUG
bool ok = SphericalCap::intersectionPoints(c1, c2, p1, p2);
Q_ASSERT(ok);
#else
SphericalCap::intersectionPoints(c1, c2, p1, p2);
#endif
Vec3d c(c1.n);
c*=c1.d;
const double a1 = std::acos((p1-c)*(p2-c)/(1.-fabs(c1.d)))/(2.*M_PI)
*area1 - OctahedronPolygon::sphericalTriangleArea(p1,p2,c1.n);
c=c2.n;
c*=c2.d;
const double a2 = std::acos((p1-c)*(p2-c)/(1.-fabs(c2.d)))/(2.*M_PI)
*area2 - OctahedronPolygon::sphericalTriangleArea(p2,p1,c2.n);
const double overlapArea = a1+a2;
return qMin(overlapArea/area1, overlapArea/area2);
}
double SphericalCap::relativeDiameterOverlap(const SphericalCap& c1, const
SphericalCap& c2)
{
if (!c1.intersects(c2))
return 0.;
if (c1.contains(c2))
return c2.getRadius()/c1.getRadius();
if (c2.contains(c1))
return c1.getRadius()/c2.getRadius();
const double r1 = c1.getRadius();
const double r2 = c2.getRadius();
const double a = c1.n.angle(c2.n);
const double overlapDist = (a-r1-r2)/(std::fabs(r1-r2)-r1-r2);
Q_ASSERT(overlapDist>=0);
return overlapDist*qMin(r1/r2, r2/r1);
}
QVector<Vec3d> SphericalCap::getClosedOutlineContour() const
{
static const int nbStep = 40;
QVector<Vec3d> contour;
Vec3d p(n);
Vec3d axis = n^Vec3d(1,0,0);
if (axis.lengthSquared()<0.1)
axis = n^Vec3d(0,1,0); // Improve precision
p.transfo4d(Mat4d::rotation(axis, std::acos(d)));
const Mat4d& rot = Mat4d::rotation(n, -2.*M_PI/nbStep);
for (int i=0;i<nbStep;++i)
{
contour.append(p);
p.transfo4d(rot);
}
return contour;
}
OctahedronPolygon SphericalCap::getOctahedronPolygon() const
{
return OctahedronPolygon(getClosedOutlineContour());
}
QVariantMap SphericalCap::toQVariant() const
{
QVariantMap res;
res.insert("type", "CAP");
double ra, dec;
StelUtils::rectToSphe(&ra, &dec, n);
QVariantList l;
l << ra*180./M_PI << dec*180./M_PI;
res.insert("center", l);
res.insert("radius", std::acos(d)*180./M_PI);
return res;
}
SphericalRegionP SphericalCap::deserialize(QDataStream& in)
{
Vec3d nn;
double dd;
in >> nn >> dd;
return SphericalRegionP(new SphericalCap(nn, dd));
}
///////////////////////////////////////////////////////////////////////////
/
// Methods for SphericalPoint
///////////////////////////////////////////////////////////////////////////
/
bool SphericalPoint::intersects(const SphericalPolygon& r) const
{
return r.contains(n);
}
bool SphericalPoint::intersects(const SphericalConvexPolygon& r) const
{
return r.contains(n);
}
OctahedronPolygon SphericalPoint::getOctahedronPolygon() const
{
QVector<Vec3d> contour;
contour << n << n << n;
return OctahedronPolygon(contour);
}
QVariantMap SphericalPoint::toQVariant() const
{
QVariantMap res;
res.insert("type", "POINT");
double ra, dec;
StelUtils::rectToSphe(&ra, &dec, n);
QVariantList l;
l << ra*180./M_PI << dec*180./M_PI;
res.insert("pos", l);
return res;
}
SphericalRegionP SphericalPoint::deserialize(QDataStream& in)
{
Vec3d nn;
in >> nn;
return SphericalRegionP(new SphericalPoint(nn));
}
///////////////////////////////////////////////////////////////////////////
/
// Methods for AllSkySphericalRegion
///////////////////////////////////////////////////////////////////////////
/
const SphericalRegionP AllSkySphericalRegion::staticInstance = SphericalReg
ionP(new AllSkySphericalRegion());
QVariantMap AllSkySphericalRegion::toQVariant() const
{
QVariantMap res;
res.insert("type", "ALLSKY");
return res;
}
///////////////////////////////////////////////////////////////////////////
/
// Methods for EmptySphericalRegion
///////////////////////////////////////////////////////////////////////////
/
const SphericalRegionP EmptySphericalRegion::staticInstance = SphericalRegi
onP(new EmptySphericalRegion());
QVariantMap EmptySphericalRegion::toQVariant() const
{
QVariantMap res;
res.insert("type", "EMPTY");
return res;
}
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalPolygon
///////////////////////////////////////////////////////////////////////////
////
SphericalCap SphericalPolygon::getBoundingCap() const
{
SphericalCap res;
octahedronPolygon.getBoundingCap(res.n, res.d);
return res;
}
QVariantMap SphericalPolygon::toQVariant() const
{
QVariantMap res;
// QVariantList worldCoordinates;
// double ra, dec;
Q_ASSERT(0);
// foreach (const QVector<Vec3d>& contour, getFillVertexArray().vertex)
// {
// QVariantList cv;
// foreach (const Vec3d& v, contour)
// {
// StelUtils::rectToSphe(&ra, &dec, v);
// QVariantList vv;
// vv << ra*180./M_PI << dec*180./M_PI;
// cv.append((QVariant)vv);
// }
// worldCoordinates.append((QVariant)cv);
// }
// res.insert("worldCoords", worldCoordinates);
return res;
}
void SphericalPolygon::serialize(QDataStream& out) const
{
out << octahedronPolygon;
}
SphericalRegionP SphericalPolygon::deserialize(QDataStream& in)
{
OctahedronPolygon p;
in >> p;
return SphericalRegionP(new SphericalPolygon(p));
}
bool SphericalPolygon::contains(const SphericalConvexPolygon& r) const {ret
urn octahedronPolygon.contains(r.getOctahedronPolygon());}
bool SphericalPolygon::intersects(const SphericalConvexPolygon& r) const {r
eturn r.intersects(*this);}
SphericalRegionP SphericalPolygon::multiUnion(const QList<SphericalRegionP>
& regions)
{
QList<OctahedronPolygon> l;
foreach (const SphericalRegionP& r, regions)
l.append(r->getOctahedronPolygon());
return SphericalRegionP(new SphericalPolygon(l));
}
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalTexturedPolygon
///////////////////////////////////////////////////////////////////////////
////
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalConvexPolygon
///////////////////////////////////////////////////////////////////////////
////
// Check if the polygon is valid, i.e. it has no side >180
bool SphericalConvexPolygon::checkValid() const
{
return SphericalConvexPolygon::checkValidContour(contour);
}
bool SphericalConvexPolygon::checkValidContour(const QVector<Vec3d>& contou
r)
{
if (contour.size()<3)
return false;
bool res=true;
for (int i=0;i<contour.size()-1;++i)
{
// Check that all points not on the current convex plane are
included in it
for (int p=0;p<contour.size()-2;++p)
res &= sideHalfSpaceContains(contour.at(i+1), contou
r.at(i), contour[(p+i+2)%contour.size()]);
}
for (int p=0;p<contour.size()-2;++p)
res &= sideHalfSpaceContains(contour.first(), contour.last()
, contour[(p+contour.size()+1)%contour.size()]);
return res;
}
//! Return the area of the region in steradians.
double SphericalConvexPolygon::getArea() const
{
double area = 0.;
Vec3d ar[3];
Vec3d v1, v2, v3;
ar[0]=contour.at(0);
for (int i=1;i<contour.size()-1;++i)
{
// Use Girard's theorem for each subtriangles
ar[1]=contour.at(i);
ar[2]=contour.at(i+1);
v1 = ar[0] ^ ar[1];
v2 = ar[1] ^ ar[2];
v3 = ar[2] ^ ar[0];
area += 2.*M_PI - v1.angle(v2) - v2.angle(v3) - v3.angle(v1)
;
}
return area;
}
//! Return a point located inside the region.
Vec3d SphericalConvexPolygon::getPointInside() const
{
Q_ASSERT(!isEmpty());
Vec3d p(contour.at(0));
p+=contour.at(1);
p+=contour.at(2);
p.normalize();
return p;
}
// Return the list of halfspace bounding the ConvexPolygon.
QVector<SphericalCap> SphericalConvexPolygon::getBoundingSphericalCaps() co
nst
{
QVector<SphericalCap> res;
for (int i=0;i<contour.size()-1;++i)
res << SphericalCap(contour.at(i+1)^contour.at(i), 0);
res << SphericalCap(contour.first()^contour.last(), 0);
return res;
}
// Returns whether a point is contained into the region.
bool SphericalConvexPolygon::contains(const Vec3d& p) const
{
if (!cachedBoundingCap.contains(p))
return false;
for (int i=0;i<contour.size()-1;++i)
{
if (!sideHalfSpaceContains(contour.at(i+1), contour.at(i), p
))
return false;
}
return sideHalfSpaceContains(contour.first(), contour.last(), p);
}
bool SphericalConvexPolygon::contains(const SphericalCap& c) const
{
if (!cachedBoundingCap.contains(c))
return false;
for (int i=0;i<contour.size()-1;++i)
{
if (!sideHalfSpaceContains(contour.at(i+1), contour.at(i), c
))
return false;
}
return sideHalfSpaceContains(contour.first(), contour.last(), c);
}
bool SphericalConvexPolygon::containsConvexContour(const Vec3d* vertice, in
t nbVertex) const
{
for (int i=0;i<nbVertex;++i)
{
if (!contains(vertice[i]))
return false;
}
return true;
}
bool SphericalConvexPolygon::contains(const SphericalConvexPolygon& cvx) co
nst
{
if (!cachedBoundingCap.contains(cvx.cachedBoundingCap))
return false;
return containsConvexContour(cvx.getConvexContour().constData(), cvx
.getConvexContour().size());
}
bool SphericalConvexPolygon::contains(const SphericalPolygon& poly) const
{
if (!cachedBoundingCap.contains(poly.getBoundingCap()))
return false;
// For standard polygons, go through the full list of triangles
const QVector<Vec3d>& vArray = poly.getFillVertexArray().vertex;
for (int i=0;i<vArray.size()/3;++i)
{
if (!containsConvexContour(vArray.constData()+i*3, 3))
return false;
}
return true;
}
bool SphericalConvexPolygon::areAllPointsOutsideOneSide(const Vec3d* thisCo
ntour, int nbThisContour, const Vec3d* points, int nbPoints)
{
for (int i=0;i<nbThisContour-1;++i)
{
bool allOutside = true;
for (int j=0;j<nbPoints&& allOutside==true;++j)
{
allOutside = allOutside && !sideHalfSpaceContains(th
isContour[i+1], thisContour[i], points[j]);
}
if (allOutside)
return true;
}
// Last iteration
bool allOutside = true;
for (int j=0;j<nbPoints&& allOutside==true;++j)
{
allOutside = allOutside && !sideHalfSpaceContains(thisContou
r[0], thisContour[nbThisContour-1], points[j]);
}
if (allOutside)
return true;
// Else
return false;
}
bool SphericalConvexPolygon::intersects(const SphericalConvexPolygon& cvx)
const
{
if (!cachedBoundingCap.intersects(cvx.cachedBoundingCap))
return false;
return !areAllPointsOutsideOneSide(cvx.contour) && !cvx.areAllPoints
OutsideOneSide(contour);
}
bool SphericalConvexPolygon::intersects(const SphericalPolygon& poly) const
{
if (!cachedBoundingCap.intersects(poly.getBoundingCap()))
return false;
// For standard polygons, go through the full list of triangles
const QVector<Vec3d>& vArray = poly.getFillVertexArray().vertex;
for (int i=0;i<vArray.size()/3;++i)
{
if (!areAllPointsOutsideOneSide(contour.constData(), contour
.size(), vArray.constData()+i*3, 3) && !areAllPointsOutsideOneSide(vArray.c
onstData()+i*3, 3, contour.constData(), contour.size()))
return true;
}
return false;
}
// This algo is wrong
void SphericalConvexPolygon::updateBoundingCap()
{
Q_ASSERT(contour.size()>2);
// Use this crapy algorithm instead
cachedBoundingCap.n.set(0,0,0);
foreach (const Vec3d& v, contour)
cachedBoundingCap.n+=v;
cachedBoundingCap.n.normalize();
cachedBoundingCap.d = 1.;
foreach (const Vec3d& v, contour)
{
if (cachedBoundingCap.n*v<cachedBoundingCap.d)
cachedBoundingCap.d = cachedBoundingCap.n*v;
}
cachedBoundingCap.d*=cachedBoundingCap.d>0 ? 0.9999999 : 1.0000001;
#ifndef NDEBUG
foreach (const Vec3d& v, contour)
Q_ASSERT(cachedBoundingCap.contains(v));
#endif
}
QVariantMap SphericalConvexPolygon::toQVariant() const
{
QVariantMap res;
res.insert("type", "CVXPOLYGON");
QVariantList cv;
double ra, dec;
foreach (const Vec3d& v, contour)
{
StelUtils::rectToSphe(&ra, &dec, v);
QVariantList vv;
vv << ra*180./M_PI << dec*180./M_PI;
cv.append((QVariant)vv);
}
res.insert("worldCoords", cv);
return res;
}
SphericalRegionP SphericalConvexPolygon::deserialize(QDataStream& in)
{
QVector<Vec3d> contour;
in >> contour;
return SphericalRegionP(new SphericalConvexPolygon(contour));
}
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalTexturedConvexPolygon
///////////////////////////////////////////////////////////////////////////
////
QVariantMap SphericalTexturedConvexPolygon::toQVariant() const
{
QVariantMap res = SphericalConvexPolygon::toQVariant();
QVariantList cv;
foreach (const Vec2f& v, textureCoords)
{
QVariantList vv;
vv << v[0] << v[1];
cv.append((QVariant)vv);
}
res.insert("textureCoords", cv);
return res;
}
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalTexturedPolygon
///////////////////////////////////////////////////////////////////////////
////
QVariantMap SphericalTexturedPolygon::toQVariant() const
{
Q_ASSERT(0);
// TODO store a tesselated polygon?, including edge flags?
return QVariantMap();
}
///////////////////////////////////////////////////////////////////////////
////
// Utility methods
///////////////////////////////////////////////////////////////////////////
////
Vec3d greatCircleIntersection(const Vec3d& p1, const Vec3d& p2, const Vec3d
& p3, const Vec3d& p4, bool& ok)
{
if (p3*p4>1.-1E-10)
{
// p3 and p4 are too close to avoid floating point problems
ok=false;
return p3;
}
Vec3d n2 = p3^p4;
n2.normalize();
return greatCircleIntersection(p1, p2, n2, ok);
}
Vec3d greatCircleIntersection(const Vec3d& p1, const Vec3d& p2, const Vec3d
& n2, bool& ok)
{
if (p1*p2>1.-1E-10)
{
// p1 and p2 are too close to avoid floating point problems
ok=false;
return p1;
}
Vec3d n1 = p1^p2;
Q_ASSERT(std::fabs(n2.lengthSquared()-1.)<0.00000001);
n1.normalize();
// Compute the parametric equation of the line at the intersection o
f the 2 planes
Vec3d u = n1^n2;
if (u.length()<1e-7)
{
// The planes are parallel
ok = false;
return u;
}
u.normalize();
// The 2 candidates are u and -u. Now need to find which point is th
e correct one.
ok = true;
n1 = p1; n1+=p2;
n1.normalize();
if (n1*u>0.)
return u;
else
return -u;
}
///////////////////////////////////////////////////////////////////////////
////
// Methods for SphericalRegionP
///////////////////////////////////////////////////////////////////////////
////
SphericalRegionP SphericalRegionP::loadFromJson(QIODevice* in)
{
StelJsonParser parser;
return loadFromQVariant(parser.parse(in).toMap());
}
SphericalRegionP SphericalRegionP::loadFromJson(const QByteArray& a)
{
QBuffer buf;
buf.setData(a);
buf.open(QIODevice::ReadOnly);
return loadFromJson(&buf);
}
inline void parseRaDec(const QVariant& vRaDec, Vec3d& v)
{
const QVariantList& vl = vRaDec.toList();
bool ok;
if (vl.size()!=2)
throw std::runtime_error(qPrintable(QString("invalid Ra,Dec
pair: \"%1\" (expect 2 double values in degree, got %2)").arg(vRaDec.toStri
ng()).arg(vl.size())));
StelUtils::spheToRect(vl.at(0).toDouble(&ok)*M_PI/180., vl.at(1).toD
ouble(&ok)*M_PI/180., v);
if (!ok)
throw std::runtime_error(qPrintable(QString("invalid Ra,Dec
pair: \"%1\" (expect 2 double values in degree)").arg(vRaDec.toString())));
}
QVector<QVector<Vec3d> > SphericalRegionP::loadContourFromQVariant(const QV
ariantList& contoursList)
{
QVector<QVector<Vec3d> > contours;
QVector<Vec3d> vertices;
bool ok;
for (int i=0;i<contoursList.size();++i)
{
const QVariantList& contourToList = contoursList.at(i).toLis
t();
if (contourToList.size()<1)
throw std::runtime_error(qPrintable(QString("invalid
contour definition: %1").arg(contoursList.at(i).toString())));
if (contourToList.at(0).toString()=="CAP")
{
// We now parse a cap, the format is "CAP",[ra, dec]
,aperture
if (contourToList.size()!=3)
throw std::runtime_error(qPrintable(QString(
"invalid CAP description: %1 (expect \"CAP\",[ra, dec],aperture)").arg(cont
oursList.at(i).toString())));
Vec3d v;
parseRaDec(contourToList.at(1), v);
double d = contourToList.at(2).toDouble(&ok)*M_PI/18
0.;
if (!ok)
throw std::runtime_error(qPrintable(QString(
"invalid aperture angle: \"%1\" (expect a double value in degree)").arg(con
tourToList.at(2).toString())));
SphericalCap cap(v,std::cos(d));
contours.append(cap.getClosedOutlineContour());
continue;
}
if (contourToList.at(0).toString()=="PATH")
{
// We now parse a path, the format is
// "PATH",[ra, dec],[["greatCircleTo", [ra, dec]], [
"smallCircle", [raAxis, decAxis], angle], [etc..]]
Q_ASSERT(vertices.isEmpty());
Vec3d v;
parseRaDec(contourToList.at(1), v);
vertices.append(v); // Starting point
foreach (const QVariant& elem, contourToList.at(2).t
oList())
{
const QVariantList& elemList = elem.toList()
;
if (elemList.size()<1)
throw std::runtime_error(qPrintable(
QString("invalid PATH description: \"%1\" (expect a list of greatCircleTo o
r smallCircle").arg(contourToList.at(2).toString())));
if (elemList.at(0)=="greatCircleTo")
{
parseRaDec(elemList.at(1), v);
vertices.append(v);
continue;
}
if (elemList.at(0)=="smallCircle")
{
Vec3d axis;
parseRaDec(elemList.at(1), axis);
double angle = elemList.at(2).toDoub
le(&ok)*M_PI/180.;
if (!ok || std::fabs(angle)>2.*M_PI)
throw std::runtime_error(qPr
intable(QString("invalid small circle rotation angle: \"%1\" (expect a doub
le value in degree betwwen -2pi and 2pi)").arg(elemList.at(2).toString())))
;
int nbStep = 1+(int)(std::fabs(angle
)/(2.*M_PI)*50);
Q_ASSERT(nbStep>0);
v = vertices.last();
const Mat4d& rotMat = Mat4d::rotatio
n(axis, angle/nbStep);
for (int step=0; step<nbStep;++step)
{
v.transfo4d(rotMat);
vertices.append(v);
}
continue;
}
throw std::runtime_error(qPrintable(QString(
"invalid PATH description: \"%1\" (expect a list of greatCircleTo or smallC
ircle").arg(contourToList.at(2).toString())));
}
Q_ASSERT(vertices.size()>2);
contours.append(vertices);
vertices.clear();
continue;
}
// If no type is provided, assume a polygon
if (contourToList.size()<3)
throw std::runtime_error("a polygon contour must hav
e at least 3 vertices");
Q_ASSERT(vertices.isEmpty());
Vec3d v;
foreach (const QVariant& vRaDec, contourToList)
{
parseRaDec(vRaDec, v);
vertices.append(v);
}
Q_ASSERT(vertices.size()>2);
contours.append(vertices);
vertices.clear();
}
return contours;
}
SphericalRegionP SphericalRegionP::loadFromQVariant(const QVariantList& con
toursList)
{
// It can only be a pure shape definition, without texture coords
const QVector<QVector<Vec3d> >& contours = loadContourFromQVariant(c
ontoursList);
return SphericalRegionP(new SphericalPolygon(contours));
}
SphericalRegionP SphericalRegionP::loadFromQVariant(const QVariantMap& map)
{
QVariantList contoursList = map.value("skyConvexPolygons").toList();
if (contoursList.empty())
contoursList = map.value("worldCoords").toList();
else
qWarning() << "skyConvexPolygons in preview JSON files is de
precated. Replace with worldCoords.";
if (contoursList.empty())
throw std::runtime_error("missing sky contours description r
equired for Spherical Geometry elements.");
// Load the matching textures positions (if any)
QVariantList texCoordList = map.value("textureCoords").toList();
if (!texCoordList.isEmpty() && contoursList.size()!=texCoordList.siz
e())
throw std::runtime_error(qPrintable(QString("the number of s
ky contours (%1) does not match the number of texture space contours (%2)")
.arg( contoursList.size()).arg(texCoordList.size())));
bool ok;
if (texCoordList.isEmpty())
{
// No texture coordinates
const QVector<QVector<Vec3d> >& contours = loadContourFromQV
ariant(contoursList);
return SphericalRegionP(new SphericalPolygon(contours));
}
else
{
// With texture coordinates
QVector<QVector<SphericalTexturedPolygon::TextureVertex> > c
ontours;
QVector<SphericalTexturedPolygon::TextureVertex> vertices;
for (int i=0;i<contoursList.size();++i)
{
// Load vertices
const QVariantList& polyRaDecToList = contoursList.a
t(i).toList();
if (polyRaDecToList.size()<3)
throw std::runtime_error("a polygon contour
must have at least 3 vertices");
SphericalTexturedPolygon::TextureVertex v;
foreach (const QVariant& vRaDec, polyRaDecToList)
{
parseRaDec(vRaDec, v.vertex);
vertices.append(v);
}
Q_ASSERT(vertices.size()>2);
// Add the texture coordinates
const QVariantList& polyXYToList = texCoordList.at(i
).toList();
if (polyXYToList.size()!=vertices.size())
throw std::runtime_error("texture coordinate
and vertices number mismatch for contour");
for (int n=0;n<polyXYToList.size();++n)
{
const QVariantList& vl = polyXYToList.at(n).
toList();
if (vl.size()!=2)
throw std::runtime_error("invalid te
xture coordinate pair (expect 2 double values in degree)");
vertices[n].texCoord.set(vl.at(0).toDouble(&
ok), vl.at(1).toDouble(&ok));
if (!ok)
throw std::runtime_error("invalid te
xture coordinate pair (expect 2 double values in degree)");
}
contours.append(vertices);
vertices.clear();
}
return SphericalRegionP(new SphericalTexturedPolygon(contour
s));
}
Q_ASSERT(0);
return SphericalRegionP(new SphericalCap());
} }
 End of changes. 45 change blocks. 
94 lines changed or deleted 1245 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/