lsst::sphgeom::Box Class Reference

Box represents a rectangle in spherical coordinate space that contains its boundary. More...

#include <Box.h>

Inheritance diagram for lsst::sphgeom::Box:

Public Member Functions
	Box ()
	This constructor creates an empty box. More...
	Box (LonLat const &p)
	This constructor creates a box containing a single point. More...
	Box (LonLat const &p1, LonLat const &p2)
	This constructor creates a box spanning the longitude interval [p1.getLon(), p2.getLon()] and latitude interval [p1.getLat(), p2.getLat()]. More...
	Box (LonLat const &p, Angle w, Angle h)
	This constructor creates a box with center p, half-width (in longitude angle) w and half-height (in latitude angle) h. More...
	Box (NormalizedAngleInterval const &lon, AngleInterval const &lat)
	This constructor creates a box spanning the given longitude and latitude intervals. More...
bool	operator== (Box const &b) const
	Two boxes are equal if they contain the same points. More...
bool	operator!= (Box const &b) const
bool	operator== (LonLat const &p) const
	A box is equal to a point p if it contains only p. More...
bool	operator!= (LonLat const &p) const
NormalizedAngleInterval const &	getLon () const
	getLon returns the longitude interval of this box. More...
AngleInterval const &	getLat () const
	getLat returns the latitude interval of this box. More...
bool	isEmpty () const
	isEmpty returns true if this box does not contain any points. More...
bool	isFull () const
	isFull returns true if this box contains all points on the unit sphere. More...
LonLat	getCenter () const
	getCenter returns the center of this box. More...
NormalizedAngle	getWidth () const
	getWidth returns the width in longitude angle of this box. More...
Angle	getHeight () const
	getHeight returns the height in latitude angle of this box. More...
Box &	clipTo (LonLat const &x)
	clipTo shrinks this box until it contains only x. More...
Box &	clipTo (Box const &x)
	x.clipTo(y) sets x to the smallest box containing the intersection of x and y. More...
Box	clippedTo (LonLat const &x) const
	clippedTo returns the intersection of this box and x. More...
Box	clippedTo (Box const &x) const
	clippedTo returns the smallest box containing the intersection of this box and x. More...
Box &	dilateBy (Angle r)
	dilateBy minimally expands this Box to include all points within angular separation r of its boundary. More...
Box	dilatedBy (Angle r) const
Box &	dilateBy (Angle w, Angle h)
	dilateBy morphologically dilates or erodes the longitude interval of this box by w, and the latitude interval of this box by h. More...
Box	dilatedBy (Angle w, Angle h) const
Box &	erodeBy (Angle r)
Box &	erodeBy (Angle w, Angle h)
Box	erodedBy (Angle r) const
Box	erodedBy (Angle w, Angle h) const
Relationship	relate (LonLat const &p) const
double	getArea () const
	getArea returns the area of this box in steradians. More...
std::unique_ptr< Region >	clone () const override
	clone returns a deep copy of this region. More...
Box	getBoundingBox () const override
	getBoundingBox returns a bounding-box for this region. More...
Box3d	getBoundingBox3d () const override
	getBoundingBox3d returns a 3-dimensional bounding-box for this region. More...
Circle	getBoundingCircle () const override
	getBoundingCircle returns a bounding-circle for this region. More...
bool	contains (UnitVector3d const &v) const override
	contains tests whether the given unit vector is inside this region. More...
Relationship	relate (Region const &r) const override
Relationship	relate (Box const &b) const override
Relationship	relate (Circle const &) const override
Relationship	relate (ConvexPolygon const &) const override
Relationship	relate (Ellipse const &) const override
std::vector< uint8_t >	encode () const override
	encode serializes this region into an opaque byte string. More...
bool	contains (LonLat const &x) const
bool	contains (Box const &x) const
bool	isDisjointFrom (LonLat const &x) const
bool	isDisjointFrom (Box const &x) const
bool	intersects (LonLat const &x) const
bool	intersects (Box const &x) const
bool	isWithin (LonLat const &x) const
bool	isWithin (Box const &x) const
Box &	expandTo (LonLat const &x)
Box &	expandTo (Box const &x)
Box	expandedTo (LonLat const &x) const
Box	expandedTo (Box const &x) const

Static Public Member Functions
static Box	fromDegrees (double lon1, double lat1, double lon2, double lat2)
static Box	fromRadians (double lon1, double lat1, double lon2, double lat2)
static Box	empty ()
static Box	full ()
static NormalizedAngle	halfWidthForCircle (Angle r, Angle lat)
	halfWidthForCircle computes the half-width of bounding boxes for circles with radius r and centers at the given latitude. More...
static NormalizedAngleInterval	allLongitudes ()
	allLongitudes returns a normalized angle interval containing all valid longitude angles. More...
static AngleInterval	allLatitudes ()
	allLatitudes returns an angle interval containing all valid latitude angles. More...
static std::unique_ptr< Box >	decode (std::vector< uint8_t > const &s)
static std::unique_ptr< Box >	decode (uint8_t const *buffer, size_t n)

Static Public Attributes
static constexpr uint8_t	TYPE_CODE = 'b'

Detailed Description

Box represents a rectangle in spherical coordinate space that contains its boundary.

A box can be empty or full (equal to the entire unit sphere), and may contain just a single point. Besides the usual rectangular regions, a box can also represent polar caps or annuli (i.e. when the box spans all longitudes).

For any instance b of this class, the following properties hold:

• b.isEmpty() == b.getLat().isEmpty()
• b.getLat().isEmpty() == b.getLon().isEmpty()
• Box::allLatitudes().contains(b.getLat())
• Box::allLongitudes().contains(b.getLon())

Definition at line 54 of file Box.h.

Constructor & Destructor Documentation

Box() [1/5]

lsst::sphgeom::Box::Box ( )

inline

This constructor creates an empty box.

Definition at line 92 of file Box.h.

{}

Box() [2/5]

lsst::sphgeom::Box::Box ( LonLat const &  p )

inlineexplicit

This constructor creates a box containing a single point.

Definition at line 95 of file Box.h.

                                   :
        _lon(p.getLon()),
        _lat(p.getLat())
    {
        _enforceInvariants();
    }

Box() [3/5]

lsst::sphgeom::Box::Box ( LonLat const &  p1, LonLat const &  p2  )

inline

This constructor creates a box spanning the longitude interval [p1.getLon(), p2.getLon()] and latitude interval [p1.getLat(), p2.getLat()].

Definition at line 105 of file Box.h.

                                              :
        _lon(p1.getLon(), p2.getLon()),
        _lat(p1.getLat(), p2.getLat())
    {
        _enforceInvariants();
    }

Box() [4/5]

lsst::sphgeom::Box::Box ( LonLat const &  p, Angle  w, Angle  h  )

inline

This constructor creates a box with center p, half-width (in longitude angle) w and half-height (in latitude angle) h.

Definition at line 114 of file Box.h.

                                            :
        _lon(NormalizedAngleInterval(p.getLon()).dilatedBy(w)),
        _lat(AngleInterval(p.getLat()).dilatedBy(h))
    {
        _enforceInvariants();
    }

Box() [5/5]

lsst::sphgeom::Box::Box ( NormalizedAngleInterval const &  lon, AngleInterval const &  lat  )

inline

This constructor creates a box spanning the given longitude and latitude intervals.

Definition at line 123 of file Box.h.

                                                                        :
        _lon(lon),
        _lat(lat)
    {
        _enforceInvariants();
    }

Member Function Documentation

allLatitudes()

static AngleInterval lsst::sphgeom::Box::allLatitudes ( )

inlinestatic

allLatitudes returns an angle interval containing all valid latitude angles.

Definition at line 87 of file Box.h.

                                        {
        return AngleInterval(Angle(-0.5 * PI), Angle(0.5 * PI));
    }

allLongitudes()

static NormalizedAngleInterval lsst::sphgeom::Box::allLongitudes ( )

inlinestatic

allLongitudes returns a normalized angle interval containing all valid longitude angles.

Definition at line 81 of file Box.h.

                                                   {
        return NormalizedAngleInterval::full();
    }

clippedTo() [1/2]

Box lsst::sphgeom::Box::clippedTo ( LonLat const &  x ) const

inline

clippedTo returns the intersection of this box and x.

Definition at line 235 of file Box.h.

{ return Box(*this).clipTo(x); }

clippedTo() [2/2]

Box lsst::sphgeom::Box::clippedTo ( Box const &  x ) const

inline

clippedTo returns the smallest box containing the intersection of this box and x.

The result is not always unique, and x.clippedTo(y) is not guaranteed to equal y.clippedTo(x).

Definition at line 240 of file Box.h.

{ return Box(*this).clipTo(x); }

clipTo() [1/2]

Box& lsst::sphgeom::Box::clipTo ( LonLat const &  x )

inline

clipTo shrinks this box until it contains only x.

If this box does not contain x, it is emptied.

Definition at line 217 of file Box.h.

                                   {
        _lon.clipTo(x.getLon());
        _lat.clipTo(x.getLat());
        _enforceInvariants();
        return *this;
    }

clipTo() [2/2]

Box& lsst::sphgeom::Box::clipTo ( Box const &  x )

inline

x.clipTo(y) sets x to the smallest box containing the intersection of x and y.

The result is not always unique, and x.clipTo(y) is not guaranteed to equal y.clipTo(x).

Definition at line 227 of file Box.h.

                                {
        _lon.clipTo(x.getLon());
        _lat.clipTo(x.getLat());
        _enforceInvariants();
        return *this;
    }

clone()

std::unique_ptr<Region> lsst::sphgeom::Box::clone ( ) const

inlineoverridevirtual

clone returns a deep copy of this region.

Implements lsst::sphgeom::Region.

Definition at line 303 of file Box.h.

                                                 {
        return std::unique_ptr<Box>(new Box(*this));
    }

contains() [1/3]

bool lsst::sphgeom::Box::contains ( LonLat const &  x ) const

inline

contains returns true if the intersection of this box and x is equal to x.

Definition at line 174 of file Box.h.

                                          {
        return _lat.contains(x.getLat()) && _lon.contains(x.getLon());
    }

contains() [2/3]

bool lsst::sphgeom::Box::contains ( Box const &  x ) const

inline

contains returns true if the intersection of this box and x is equal to x.

Definition at line 178 of file Box.h.

                                       {
        return _lat.contains(x._lat) && _lon.contains(x._lon);
    }

contains() [3/3]

bool lsst::sphgeom::Box::contains ( UnitVector3d const &  ) const

inlineoverridevirtual

contains tests whether the given unit vector is inside this region.

Implements lsst::sphgeom::Region.

Definition at line 311 of file Box.h.

                                                         {
        return contains(LonLat(v));
    }

decode() [1/2]

static std::unique_ptr<Box> lsst::sphgeom::Box::decode ( std::vector< uint8_t > const &  s )

inlinestatic

decode deserializes a Box from a byte string produced by encode.

Definition at line 338 of file Box.h.

                                                                   {
        return decode(s.data(), s.size());
    }

decode() [2/2]

std::unique_ptr< Box > lsst::sphgeom::Box::decode ( uint8_t const *  buffer, size_t  n  )

static

decode deserializes a Box from a byte string produced by encode.

Definition at line 459 of file Box.cc.

                                                               {
    if (buffer == nullptr || n != ENCODED_SIZE || *buffer != TYPE_CODE) {
        throw std::runtime_error("Byte-string is not an encoded Box");
    }
    std::unique_ptr<Box> box(new Box);
    ++buffer;
    double a = decodeDouble(buffer); buffer += 8;
    double b = decodeDouble(buffer); buffer += 8;
    box->_lon = NormalizedAngleInterval::fromRadians(a, b);
    a = decodeDouble(buffer); buffer += 8;
    b = decodeDouble(buffer); buffer += 8;
    box->_lat = AngleInterval::fromRadians(a, b);
    box->_enforceInvariants();
    return box;
}

dilateBy() [1/2]

Box & lsst::sphgeom::Box::dilateBy

(

Angle

r

)

dilateBy minimally expands this Box to include all points within angular separation r of its boundary.

If this box is empty or full, or if r is non-positive, there is no effect.

Definition at line 78 of file Box.cc.

                           {
    // The basic idea is to compute the union of the bounding boxes for all
    // circles of opening angle r with centers inside this box.
    //
    // The bounding box for a circle of opening angle r with center latitude
    // |δ| ≤ π/2 - r has height 2r.
    //
    // Given fixed r, the width of the bounding box for the circle centered at
    // latitude δ grows monotonically with |δ| - for justification, see the
    // derivation in halfWidthForCircle(). The maximum width is therefore
    // attained when the circle is centered at one of the latitude angle
    // boundaries of this box. If max(|δ|) ≥ π/2 - r, it is 2π.
    //
    // Dilating the longitude interval of this box by the maximum width and
    // the latitude interval by r gives the desired result.
    if (isEmpty() || isFull() || r <= Angle(0.0)) {
        return *this;
    }
    Angle maxAbsLatitude = std::max(abs(_lat.getA()), abs(_lat.getB()));
    NormalizedAngle w = halfWidthForCircle(r, maxAbsLatitude);
    return dilateBy(w, r);
}

dilateBy() [2/2]

Box & lsst::sphgeom::Box::dilateBy	(	Angle	w,
		Angle	h
	)

dilateBy morphologically dilates or erodes the longitude interval of this box by w, and the latitude interval of this box by h.

If w is positive, the longitude interval is dilated by [-w,w]. If w is zero, the corresponding interval is not modified, and if it is negative, the longitude interval is eroded by [w,-w]. The action of h on the latitude interval is analogous.

If this box is empty or full, there is no effect. Furthermore, a box containing the north or south pole is not considered to have a latitude boundary there, so that eroding it will not necessarily remove the pole.

If the desired outcome is the bounding box of points within some angular separation r of this box, then dilateBy(r) must be called, not dilateBy(r, r).

Definition at line 101 of file Box.cc.

                                    {
    if (isEmpty() || isFull()) {
        return *this;
    }
    _lon.dilateBy(w);
    if (!h.isNan()) {
        Angle a = (_lat.getA() > Angle(-0.5 * PI)) ? _lat.getA() - h : _lat.getA();
        Angle b = (_lat.getB() < Angle(0.5 * PI)) ? _lat.getB() + h : _lat.getB();
        _lat = AngleInterval(a, b);
    }
    _enforceInvariants();
    return *this;
}

dilatedBy() [1/2]

Box lsst::sphgeom::Box::dilatedBy ( Angle  r ) const

inline

Definition at line 273 of file Box.h.

{ return Box(*this).dilateBy(r); }

dilatedBy() [2/2]

Box lsst::sphgeom::Box::dilatedBy ( Angle  w, Angle  h  ) const

inline

Definition at line 291 of file Box.h.

{ return Box(*this).dilateBy(w, h); }

empty()

static Box lsst::sphgeom::Box::empty ( )

inlinestatic

Definition at line 69 of file Box.h.

{ return Box(); }

encode()

std::vector< uint8_t > lsst::sphgeom::Box::encode ( ) const

overridevirtual

encode serializes this region into an opaque byte string.

Byte strings emitted by encode can be deserialized with decode.

Implements lsst::sphgeom::Region.

Definition at line 447 of file Box.cc.

                                     {
    std::vector<uint8_t> buffer;
    uint8_t tc = TYPE_CODE;
    buffer.reserve(ENCODED_SIZE);
    buffer.push_back(tc);
    encodeDouble(_lon.getA().asRadians(), buffer);
    encodeDouble(_lon.getB().asRadians(), buffer);
    encodeDouble(_lat.getA().asRadians(), buffer);
    encodeDouble(_lat.getB().asRadians(), buffer);
    return buffer;
}

erodeBy() [1/2]

Box& lsst::sphgeom::Box::erodeBy ( Angle  r )

inline

Definition at line 292 of file Box.h.

{ return dilateBy(-r); }

erodeBy() [2/2]

Box& lsst::sphgeom::Box::erodeBy ( Angle  w, Angle  h  )

inline

Definition at line 293 of file Box.h.

{ return dilateBy(-w, -h); }

erodedBy() [1/2]

Box lsst::sphgeom::Box::erodedBy ( Angle  r ) const

inline

Definition at line 294 of file Box.h.

{ return dilatedBy(-r);  }

erodedBy() [2/2]

Box lsst::sphgeom::Box::erodedBy ( Angle  w, Angle  h  ) const

inline

Definition at line 295 of file Box.h.

{ return dilatedBy(-w, -h);  }

expandedTo() [1/2]

Box lsst::sphgeom::Box::expandedTo ( LonLat const &  x ) const

inline

expandedTo returns the smallest box containing the union of this box and x. The result is not always unique, and x.expandedTo(y) is not guaranteed to equal y.expandedTo(x).

Definition at line 263 of file Box.h.

{ return Box(*this).expandTo(x); }

expandedTo() [2/2]

Box lsst::sphgeom::Box::expandedTo ( Box const &  x ) const

inline

expandedTo returns the smallest box containing the union of this box and x. The result is not always unique, and x.expandedTo(y) is not guaranteed to equal y.expandedTo(x).

Definition at line 264 of file Box.h.

{ return Box(*this).expandTo(x); }

expandTo() [1/2]

Box& lsst::sphgeom::Box::expandTo ( LonLat const &  x )

inline

expandTo minimally expands this box to contain x. The result is not always unique, and x.expandTo(y) is not guaranteed to equal y.expandTo(x).

Definition at line 246 of file Box.h.

                                     {
        _lon.expandTo(x.getLon());
        _lat.expandTo(x.getLat());
        return *this;
    }

expandTo() [2/2]

Box& lsst::sphgeom::Box::expandTo ( Box const &  x )

inline

expandTo minimally expands this box to contain x. The result is not always unique, and x.expandTo(y) is not guaranteed to equal y.expandTo(x).

Definition at line 252 of file Box.h.

                                  {
        _lon.expandTo(x.getLon());
        _lat.expandTo(x.getLat());
        return *this;
    }

fromDegrees()

static Box lsst::sphgeom::Box::fromDegrees ( double  lon1, double  lat1, double  lon2, double  lat2  )

inlinestatic

Definition at line 59 of file Box.h.

                                                                               {
        return Box(NormalizedAngleInterval::fromDegrees(lon1, lon2),
                   AngleInterval::fromDegrees(lat1, lat2));
    }

fromRadians()

static Box lsst::sphgeom::Box::fromRadians ( double  lon1, double  lat1, double  lon2, double  lat2  )

inlinestatic

Definition at line 64 of file Box.h.

                                                                               {
        return Box(NormalizedAngleInterval::fromRadians(lon1, lon2),
                   AngleInterval::fromRadians(lat1, lat2));
    }

full()

static Box lsst::sphgeom::Box::full ( )

inlinestatic

Definition at line 71 of file Box.h.

{ return Box(allLongitudes(), allLatitudes()); }

getArea()

double lsst::sphgeom::Box::getArea

(

)

const

getArea returns the area of this box in steradians.

Definition at line 115 of file Box.cc.

                          {
    if (isEmpty()) {
        return 0.0;
    }
    // Given a, b ∈ [-π/2, π/2] and a std::sin implementation that is not
    // correctly rounded, b > a does not imply that std::sin(b) > std::sin(a).
    // To avoid potentially returning a negative area, defensively take an
    // absolute value.
    double dz = sin(_lat.getB()) - sin(_lat.getA());
    return std::fabs(_lon.getSize().asRadians() * dz);
}

getBoundingBox()

Box lsst::sphgeom::Box::getBoundingBox ( ) const

inlineoverridevirtual

getBoundingBox returns a bounding-box for this region.

Implements lsst::sphgeom::Region.

Definition at line 307 of file Box.h.

{ return *this; }

getBoundingBox3d()

Box3d lsst::sphgeom::Box::getBoundingBox3d ( ) const

overridevirtual

getBoundingBox3d returns a 3-dimensional bounding-box for this region.

Implements lsst::sphgeom::Region.

Definition at line 127 of file Box.cc.

                                  {
    if (isEmpty()) {
        return Box3d();
    }
    if (isFull()) {
        return Box3d::aroundUnitSphere();
    }
    double slata = sin(_lat.getA()), clata = cos(_lat.getA());
    double slatb = sin(_lat.getB()), clatb = cos(_lat.getB());
    double slona = sin(_lon.getA()), clona = cos(_lon.getA());
    double slonb = sin(_lon.getB()), clonb = cos(_lon.getB());
    // Compute the minimum/maximum x/y values of the box vertices.
    double xmin = std::min(std::min(clona * clata, clonb * clata),
                           std::min(clona * clatb, clonb * clatb)) - 2.5 * EPSILON;
    double xmax = std::max(std::max(clona * clata, clonb * clata),
                           std::max(clona * clatb, clonb * clatb)) + 2.5 * EPSILON;
    double ymin = std::min(std::min(slona * clata, slonb * clata),
                           std::min(slona * clatb, slonb * clatb)) - 2.5 * EPSILON;
    double ymax = std::max(std::max(slona * clata, slonb * clata),
                           std::max(slona * clatb, slonb * clatb)) + 2.5 * EPSILON;
    // Compute the maximum latitude cosine of points in the box.
    double mlc;
    if (_lat.contains(Angle(0.0))) {
        mlc = 1.0;
        // The box intersects the equator - the x or y extrema of the box may be
        // at the intersection of the box edge meridians with the equator.
        xmin = std::min(xmin, std::min(clona, clonb) - EPSILON);
        xmax = std::max(xmax, std::max(clona, clonb) + EPSILON);
        ymin = std::min(ymin, std::min(slona, slonb) - EPSILON);
        ymax = std::max(ymax, std::max(slona, slonb) + EPSILON);
    } else {
        // Note that clata and clatb are positive.
        mlc = std::max(clata, clatb) + EPSILON;
    }
    // Check for extrema on the box edges parallel to the equator.
    if (_lon.contains(NormalizedAngle(0.0))) {
        xmax = std::max(xmax, mlc);
    }
    if (_lon.contains(NormalizedAngle(0.5 * PI))) {
        ymax = std::max(ymax, mlc);
    }
    if (_lon.contains(NormalizedAngle(PI))) {
        xmin = std::min(xmin, -mlc);
    }
    if (_lon.contains(NormalizedAngle(1.5 * PI))) {
        ymin = std::min(ymin, -mlc);
    }
    // Clamp x/y extrema to [-1, 1]
    xmin = std::max(-1.0, xmin);
    xmax = std::min(1.0, xmax);
    ymin = std::max(-1.0, ymin);
    ymax = std::min(1.0, ymax);
    // Compute z extrema.
    double zmin = std::max(-1.0, slata - EPSILON);
    double zmax = std::min(1.0, slatb + EPSILON);
    return Box3d(Interval1d(xmin, xmax),
                 Interval1d(ymin, ymax),
                 Interval1d(zmin, zmax));
}

getBoundingCircle()

Circle lsst::sphgeom::Box::getBoundingCircle ( ) const

overridevirtual

getBoundingCircle returns a bounding-circle for this region.

Implements lsst::sphgeom::Region.

Definition at line 187 of file Box.cc.

                                    {
    if (isEmpty()) {
        return Circle::empty();
    }
    if (isFull()) {
        return Circle::full();
    }
    NormalizedAngle w = getWidth();
    // The minimal bounding circle center p lies on the meridian bisecting
    // this box. Let δ₁ and δ₂ be the minimum and maximum box latitudes.
    if (w.asRadians() <= PI) {
        UnitVector3d p;
        UnitVector3d boxVerts[4] = {
            UnitVector3d(_lon.getA(), _lat.getA()),
            UnitVector3d(_lon.getA(), _lat.getB()),
            UnitVector3d(_lon.getB(), _lat.getA()),
            UnitVector3d(_lon.getB(), _lat.getB())
        };
        // We take advantage of rotational symmetry to fix the bisecting
        // meridian at a longitude of zero. The box vertices then have
        // coordinates (±w/2, δ₁), (±w/2, δ₂), and p = (0, ϕ). Converting
        // to Cartesian coordinates gives p = (cos ϕ, 0, sin ϕ), and box
        // vertices at (cos w/2 cos δ₁, ±sin w/2 cos δ₁, sin δ₁) and
        // (cos w/2 cos δ₂, ±sin w/2 cos δ₂, sin δ₂).
        //
        // The point p₁ on the meridian that has minimum angular separation
        // to the vertices with latitude δ₁ lies on the plane they define.
        // The sum of the two vertex vectors is on that plane and on the plane
        // containing the meridian. Normalizing to obtain p₁, we have
        //
        //     (cos ϕ₁, 0, sin ϕ₁) =
        //         λ ((cos w/2 cos δ₁,  sin w/2 cos δ₁, sin δ₁) +
        //            (cos w/2 cos δ₁, -sin w/2 cos δ₁, sin δ₁))
        //
        // for some scaling factor λ. Simplifying, we get:
        //
        //     cos ϕ₁ = λ cos w/2 cos δ₁
        //     sin ϕ₁ = λ sin δ₁
        //
        // so that
        //
        //     tan ϕ₁ = sec w/2 tan δ₁
        //
        // Similarly, the point p₂ on the meridian that has minimum angular
        // separation to the vertices with latitude δ₂ satisfies:
        //
        //     tan ϕ₂ = sec w/2 tan δ₂
        //
        // where ϕ₁ ≤ ϕ₂ (since δ₁ ≤ δ₂). Finally, consider the point p₃
        // separated from each box vertex by the same angle. The dot
        // products of p₃ with the box vertices are all identical, so
        //
        //     cos ϕ₃ cos w/2 cos δ₁ + sin ϕ₃ sin δ₁ =
        //     cos ϕ₃ cos w/2 cos δ₂ + sin ϕ₃ sin δ₂
        //
        // Rearranging gives:
        //
        //     tan ϕ₃ = - cos w/2 (cos δ₁ - cos δ₂)/(sin δ₁ - sin δ₂)
        //
        // which can be simplified further using a tangent half-angle identity,
        // yielding:
        //
        //     tan ϕ₃ = cos w/2 tan (δ₁ + δ₂)/2
        //
        // Consider now the function f₁(ϕ) that gives the angular separation
        // between p with latitude ϕ and the vertices at latitude δ₁. It has
        // a line of symmetry at ϕ = ϕ₁, and increases monotonically with
        // |ϕ - ϕ₁|. Similarly, f₂(ϕ) has a minimum at ϕ₂ and increases
        // monotonically with |ϕ - ϕ₂|. The two functions cross at ϕ₃. The
        // opening angle of the bounding circle centered at latitude ϕ is
        // given by g = max(f₁, f₂), which we seek to minimize.
        //
        // If ϕ₁ ≤ ϕ₃ ≤ ϕ₂, then g is minimized at ϕ = ϕ₃. Otherwise, it
        // is minimized at either ϕ₁ or ϕ₂.
        double phi1, phi2, phi3;
        double c = cos(0.5 * w);
        if (c == 0.0) {
            // This code should never execute. If it does, the implementation
            // of std::cos is broken.
            phi1 = ::copysign(0.5 * PI, _lat.getA().asRadians());
            phi2 = ::copysign(0.5 * PI, _lat.getB().asRadians());
            phi3 = 0.0;
        } else {
            phi1 = std::atan(tan(_lat.getA()) / c);
            phi2 = std::atan(tan(_lat.getB()) / c);
            phi3 = std::atan(c * tan(_lat.getCenter()));
        }
        if (phi1 <= phi3 && phi3 <= phi2) {
            p = UnitVector3d(_lon.getCenter(), Angle(phi3));
        } else {
            UnitVector3d p1 = UnitVector3d(_lon.getCenter(), Angle(phi1));
            UnitVector3d p2 = UnitVector3d(_lon.getCenter(), Angle(phi2));
            if (p1.dot(boxVerts[0]) > p2.dot(boxVerts[1])) {
                p = p2;
            } else {
                p = p1;
            }
        }
        // Compute the maximum squared chord length between p and the box
        // vertices, so that each one is guaranteed to lie in the bounding
        // circle, regardless of numerical error in the above.
        double cl2 = (p - boxVerts[0]).getSquaredNorm();
        for (int i = 1; i < 4; ++i) {
            cl2 = std::max(cl2, (p - boxVerts[i]).getSquaredNorm());
        }
        // Add double the maximum squared-chord-length error, so that the
        // bounding circle we return also reliably CONTAINS this box.
        return Circle(p, cl2 + 2.0 * MAX_SQUARED_CHORD_LENGTH_ERROR);
    }
    // The box spans more than π radians in longitude. First, pick the smaller
    // of the bounding circles centered at the north and south pole.
    Angle r;
    UnitVector3d v;
    if (abs(_lat.getA()) <= abs(_lat.getB())) {
        v = UnitVector3d::Z();
        r = Angle(0.5 * PI) - _lat.getA();
    } else {
        v = -UnitVector3d::Z();
        r = _lat.getB() + Angle(0.5 * PI);
    }
    // If the box does not span all longitude angles, we also consider the
    // equatorial bounding circle with center longitude equal to the longitude
    // of the box center. The smaller of the polar and equatorial bounding
    // circles is returned.
    if (!_lon.isFull() && 0.5 * w < r) {
        r = 0.5 * w;
        v = UnitVector3d(_lon.getCenter(), Angle(0.0));
    }
    return Circle(v, r + 4.0 * Angle(MAX_ASIN_ERROR));
}

getCenter()

LonLat lsst::sphgeom::Box::getCenter ( ) const

inline

getCenter returns the center of this box.

It is NaN for empty boxes and arbitrary for full boxes.

Definition at line 159 of file Box.h.

                             {
        return LonLat(_lon.getCenter(), _lat.getCenter());
    }

getHeight()

Angle lsst::sphgeom::Box::getHeight ( ) const

inline

getHeight returns the height in latitude angle of this box.

It is negative or NaN for empty boxes.

Definition at line 169 of file Box.h.

{ return _lat.getSize(); }

getLat()

AngleInterval const& lsst::sphgeom::Box::getLat ( ) const

inline

getLat returns the latitude interval of this box.

Definition at line 148 of file Box.h.

{ return _lat; }

getLon()

NormalizedAngleInterval const& lsst::sphgeom::Box::getLon ( ) const

inline

getLon returns the longitude interval of this box.

Definition at line 145 of file Box.h.

{ return _lon; }

getWidth()

NormalizedAngle lsst::sphgeom::Box::getWidth ( ) const

inline

getWidth returns the width in longitude angle of this box.

It is NaN for empty boxes.

Definition at line 165 of file Box.h.

{ return _lon.getSize(); }

halfWidthForCircle()

NormalizedAngle lsst::sphgeom::Box::halfWidthForCircle ( Angle  r, Angle  lat  )

static

halfWidthForCircle computes the half-width of bounding boxes for circles with radius r and centers at the given latitude.

If r is non-positive, the result is zero, and if |lat| + r >= PI/2, the result is PI.

Definition at line 43 of file Box.cc.

                                                          {
    if (r <= Angle(0.0)) {
        return NormalizedAngle(0.0);
    }
    // If a circle centered at the given latitude contains a pole, then
    // its bounding box contains all possible longitudes.
    if (abs(lat) + r >= Angle(0.5 * PI)) {
        return NormalizedAngle(PI);
    }
    // Now, consider the circle with opening angle r > 0 centered at (0,δ)
    // with r < π/2 and |δ| ≠ π/2. The circle center vector in ℝ³ is
    // c = (cos δ, 0, sin δ). Its bounding box spans longitudes [-α,α], where
    // α is the desired half-width. The plane corresponding to longitude α has
    // normal vector (-sin α, cos α, 0) and is tangent to the circle at point
    // p. The great circle segment between the center of the circle and the
    // plane normal passes through p and has arc length π/2 + r, so that
    //
    //    (cos δ, 0, sin δ) · (-sin α, cos α, 0) = cos (π/2 + r)
    //
    // Solving for α gives
    //
    //    α = arcsin (sin r / cos δ)
    //
    // In the actual computation, there is an absolute value and an explicit
    // arcsin domain check to cope with rounding errors. An alternate way to
    // compute this is:
    //
    //    α = arctan (sin r / √(cos(δ - r) cos(δ + r)))
    double s = std::fabs(sin(r) / cos(lat));
    if (s >= 1.0) {
        return NormalizedAngle(0.5 * PI);
    }
    return NormalizedAngle(std::asin(s));
}

intersects() [1/2]

bool lsst::sphgeom::Box::intersects ( LonLat const &  x ) const

inline

intersects returns true if the intersection of this box and x is non-empty.

Definition at line 194 of file Box.h.

                                            {
        return _lat.intersects(x.getLat()) && _lon.intersects(x.getLon());
    }

intersects() [2/2]

bool lsst::sphgeom::Box::intersects ( Box const &  x ) const

inline

intersects returns true if the intersection of this box and x is non-empty.

Definition at line 198 of file Box.h.

                                         {
        return _lat.intersects(x._lat) && _lon.intersects(x._lon);
    }

isDisjointFrom() [1/2]

bool lsst::sphgeom::Box::isDisjointFrom ( LonLat const &  x ) const

inline

isDisjointFrom returns true if the intersection of this box and x is empty.

Definition at line 186 of file Box.h.

{ return !intersects(x); }

isDisjointFrom() [2/2]

bool lsst::sphgeom::Box::isDisjointFrom ( Box const &  x ) const

inline

isDisjointFrom returns true if the intersection of this box and x is empty.

Definition at line 188 of file Box.h.

{ return !intersects(x); }

isEmpty()

bool lsst::sphgeom::Box::isEmpty ( ) const

inline

isEmpty returns true if this box does not contain any points.

Definition at line 151 of file Box.h.

{ return _lat.isEmpty(); }

isFull()

bool lsst::sphgeom::Box::isFull ( ) const

inline

isFull returns true if this box contains all points on the unit sphere.

Definition at line 155 of file Box.h.

{ return _lon.isFull() && _lat == allLatitudes(); }

isWithin() [1/2]

bool lsst::sphgeom::Box::isWithin ( LonLat const &  x ) const

inline

isWithin returns true if the intersection of this box and x is this box.

Definition at line 206 of file Box.h.

                                          {
        return _lat.isWithin(x.getLat()) && _lon.isWithin(x.getLon());
    }

isWithin() [2/2]

bool lsst::sphgeom::Box::isWithin ( Box const &  x ) const

inline

isWithin returns true if the intersection of this box and x is this box.

Definition at line 210 of file Box.h.

                                       {
        return _lat.isWithin(x._lat) && _lon.isWithin(x._lon);
    }

operator!=() [1/2]

bool lsst::sphgeom::Box::operator!= ( Box const &  b ) const

inline

Definition at line 135 of file Box.h.

{ return !(*this == b); }

operator!=() [2/2]

bool lsst::sphgeom::Box::operator!= ( LonLat const &  p ) const

inline

Definition at line 142 of file Box.h.

{ return !(*this == p); }

operator==() [1/2]

bool lsst::sphgeom::Box::operator== ( Box const &  b ) const

inline

Two boxes are equal if they contain the same points.

Definition at line 131 of file Box.h.

                                         {
        return _lon == b._lon && _lat == b._lat;
    }

operator==() [2/2]

bool lsst::sphgeom::Box::operator== ( LonLat const &  p ) const

inline

A box is equal to a point p if it contains only p.

Definition at line 138 of file Box.h.

                                            {
        return _lat == p.getLat() && _lon == p.getLon();
    }

relate() [1/6]

Relationship lsst::sphgeom::Box::relate ( LonLat const &  p ) const

inline

Definition at line 297 of file Box.h.

{ return relate(Box(p)); }

relate() [2/6]

Relationship lsst::sphgeom::Box::relate ( Region const &  ) const

inlineoverridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

• Bit S & DISJOINT is set only if A and B do not have any points in common.
• Bit S & CONTAINS is set only if A contains all points in B.
• Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 315 of file Box.h.

                                                         {
        // Dispatch on the type of r.
        return invert(r.relate(*this));
    }

relate() [3/6]

Relationship lsst::sphgeom::Box::relate ( Box const &  ) const

inlineoverridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

• Bit S & DISJOINT is set only if A and B do not have any points in common.
• Bit S & CONTAINS is set only if A contains all points in B.
• Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 320 of file Box.h.

                                                      {
        Relationship r1 = _lon.relate(b._lon);
        Relationship r2 = _lat.relate(b._lat);
        // If the box longitude or latitude intervals are disjoint, then the
        // boxes are disjoint. The other spatial relationships must hold for
        // both the longitude and latitude intervals in order to hold for the
        // boxes.
        return ((r1 & r2) & (CONTAINS | WITHIN)) | ((r1 | r2) & DISJOINT);
    }

relate() [4/6]

Relationship lsst::sphgeom::Box::relate ( Circle const &  ) const

overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

• Bit S & DISJOINT is set only if A and B do not have any points in common.
• Bit S & CONTAINS is set only if A contains all points in B.
• Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 318 of file Box.cc.

                                               {
    if (isEmpty()) {
        if (c.isEmpty()) {
            return CONTAINS | DISJOINT | WITHIN;
        }
        return DISJOINT | WITHIN;
    } else if (c.isEmpty()) {
        return CONTAINS | DISJOINT;
    }
    if (isFull()) {
        if (c.isFull()) {
            return CONTAINS | WITHIN;
        }
        return CONTAINS;
    } else if (c.isFull()) {
        return WITHIN;
    }
    // Neither region is empty or full. We now determine whether or not the
    // circle and box boundaries intersect.
    //
    // If the box vertices are not all inside or all outside of c, then the
    // boundaries cross.
    LonLat vertLonLat[4] = {
        LonLat(_lon.getA(), _lat.getA()),
        LonLat(_lon.getA(), _lat.getB()),
        LonLat(_lon.getB(), _lat.getA()),
        LonLat(_lon.getB(), _lat.getB())
    };
    UnitVector3d verts[4];
    bool inside = false;
    for (int i = 0; i < 4; ++i) {
        verts[i] = UnitVector3d(vertLonLat[i]);
        double d = (verts[i] - c.getCenter()).getSquaredNorm();
        if (std::fabs(d - c.getSquaredChordLength()) <
            MAX_SQUARED_CHORD_LENGTH_ERROR) {
            // A box vertex is close to the circle boundary.
            return INTERSECTS;
        }
        bool b = d < c.getSquaredChordLength();
        if (i == 0) {
            inside = b;
        } else if (inside != b) {
            // There are box vertices both inside and outside of c.
            return INTERSECTS;
        }
    }
    UnitVector3d norms[2] = {
        UnitVector3d::orthogonalTo(_lon.getA()),
        UnitVector3d::orthogonalTo(_lon.getB())
    };
    if (inside) {
        // All box vertices are inside c. Look for points in the box edge
        // interiors that are outside c.
        for (int i = 0; i < 2; ++i) {
            double d = getMaxSquaredChordLength(
                c.getCenter(), verts[2 * i + 1], verts[2 * i], norms[i]);
            if (d > c.getSquaredChordLength() -
                    MAX_SQUARED_CHORD_LENGTH_ERROR) {
                return INTERSECTS;
            }
        }
        LonLat cc(-c.getCenter());
        if (_lon.contains(cc.getLon())) {
            // The points furthest from the center of c on the small circles
            // defined by the box edges with constant latitude are in the box
            // edge interiors. Find the largest squared chord length to either.
            Angle a = std::min(getMinAngleToCircle(cc.getLat(), _lat.getA()),
                               getMinAngleToCircle(cc.getLat(), _lat.getB()));
            double d = Circle::squaredChordLengthFor(Angle(PI) - a);
            if (d > c.getSquaredChordLength() -
                    MAX_SQUARED_CHORD_LENGTH_ERROR) {
                return INTERSECTS;
            }
        }
        // The box boundary is completely inside c. However, the box is not
        // necessarily within c: consider a circle with opening angle equal to
        // π - ε. If a box contains the complement of such a circle, then
        // intersecting it with that circle will punch a hole in the box. In
        // this case each region contains the boundary of the other, but
        // neither region contains the other.
        //
        // To handle this case, check that the box does not contain the
        // complement of c - since the boundaries do not intersect, this is the
        // case iff the box contains the center of the complement of c.
        if (contains(cc)) {
            return INTERSECTS;
        }
        return WITHIN;
    }
    // All box vertices are outside c. Look for points in the box edge
    // interiors that are inside c.
    for (int i = 0; i < 2; ++i) {
        double d = getMinSquaredChordLength(
            c.getCenter(), verts[2 * i + 1], verts[2 * i], norms[i]);
        if (d < c.getSquaredChordLength() + MAX_SQUARED_CHORD_LENGTH_ERROR) {
            return INTERSECTS;
        }
    }
    LonLat cc(c.getCenter());
    if (_lon.contains(cc.getLon())) {
        // The points closest to the center of c on the small circles
        // defined by the box edges with constant latitude are in the box
        // edge interiors. Find the smallest squared chord length to either.
        Angle a = std::min(getMinAngleToCircle(cc.getLat(), _lat.getA()),
                           getMinAngleToCircle(cc.getLat(), _lat.getB()));
        double d = Circle::squaredChordLengthFor(a);
        if (d < c.getSquaredChordLength() + MAX_SQUARED_CHORD_LENGTH_ERROR) {
            return INTERSECTS;
        }
    }
    // The box boundary is completely outside of c. If the box contains the
    // circle center, then the box contains c. Otherwise, the box and circle
    // are disjoint.
    if (contains(cc)) {
        return CONTAINS;
    }
    return DISJOINT;
}

relate() [5/6]

Relationship lsst::sphgeom::Box::relate ( ConvexPolygon const &  ) const

overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

• Bit S & DISJOINT is set only if A and B do not have any points in common.
• Bit S & CONTAINS is set only if A contains all points in B.
• Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 437 of file Box.cc.

                                                      {
    // ConvexPolygon-Box relations are implemented by ConvexPolygon.
    return invert(p.relate(*this));
}

relate() [6/6]

Relationship lsst::sphgeom::Box::relate ( Ellipse const &  ) const

overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

• Bit S & DISJOINT is set only if A and B do not have any points in common.
• Bit S & CONTAINS is set only if A contains all points in B.
• Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 442 of file Box.cc.

                                                {
    // Ellipse-Box relations are implemented by Ellipse.
    return invert(e.relate(*this));
}

Member Data Documentation

TYPE_CODE

constexpr uint8_t lsst::sphgeom::Box::TYPE_CODE = 'b'

static

Definition at line 56 of file Box.h.

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/Linux64/sphgeom/17.0+14/include/lsst/sphgeom/Box.h
• /j/snowflake/release/lsstsw/stack/Linux64/sphgeom/17.0+14/src/Box.cc