SpatialUtil.h

Go to the documentation of this file.

// -*- lsst-c++ -*-

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 

#ifndef LSST_AP_SPATIAL_UTIL_H
#define LSST_AP_SPATIAL_UTIL_H

#include <cassert>
#include <cmath>

#include <algorithm>
#include <vector>

#include "Common.h"
#include "CircularRegion.h"
#include "RectangularRegion.h"


namespace lsst { namespace ap {

namespace {
    double const RA_DEC_SCALE = 1.19304647111111111111111111111e07; // 1073741824.0/90.0
}

class ZoneStripeChunkDecomposition {

public :

    ZoneStripeChunkDecomposition(
        int const zonesPerDegree,
        int const zonesPerStripe,
        int const maxEntriesPerZoneEstimate
    );

    ZoneStripeChunkDecomposition(ZoneStripeChunkDecomposition const & zsc);

    ZoneStripeChunkDecomposition & operator=(ZoneStripeChunkDecomposition const & zsc) {
        if (this != &zsc) {
            ZoneStripeChunkDecomposition copy(zsc);
            swap(copy);
        }
        return *this;
    }

    int getNumChunksPerStripe(int const stripeId, double const minWidth) const;

    double stripeAndCircleToRaRange(
        int    const stripeId,
        double const cenRa,
        double const cenDec,
        double const rad
    ) const;

    int getNumChunksPerStripe(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        return _chunksPerStripe[stripeId - _minStripe];
    }

    int getFirstChunkForStripe(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        return stripeId << 16;
    }

    int radecToChunk(double const ra, double const dec) const {
        int const s  = decToStripe(dec);
        int const nc = getNumChunksPerStripe(s);
        int c = static_cast<int>(std::floor((ra/360.0)*nc));
        if (c >= nc) {
            --c;
        }
        return (c + getFirstChunkForStripe(s));
    }

    int decToZone(double const dec) const {
        assert(dec >= -90.0 && dec <= 90.0 && "declination out of bounds");
        int const zoneId = static_cast<int>(std::floor(dec*_zonesPerDegree));
        return (zoneId > _maxZone) ? _maxZone : zoneId;
    }

    double getZoneDecMin(int const zone) const {
        assert(zone >= _minZone && zone <= _maxZone && "zone id out of bounds");
        double d = static_cast<double>(zone)/_zonesPerDegree;
        return d <= -90.0 ? -90.0 : d;
    }

    double getZoneDecMax(int const zone) const {
        assert(zone >= _minZone && zone <= _maxZone && "zone id out of bounds");
        double d = static_cast<double>(zone + 1)/_zonesPerDegree;
        return d >= 90.0 ? 90.0 : d;
    }

    int decToStripe(double const dec) const {
        assert(dec >= -90.0 && dec <= 90.0 && "declination out of bounds");
        double zoneId = std::floor(dec*_zonesPerDegree);
        return (zoneId > _maxZoneAsDouble) ? _maxStripe :
            static_cast<int>(std::floor(zoneId/_zonesPerStripe));
    }

    static int chunkToStripe(int const chunkId) {
        return chunkId >> 16;
    }

    static int chunkToSequence(int const chunkId) {
        return chunkId & 0x7fff;
    }

    double getStripeDecMin(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        double d = static_cast<double>(stripeId*_zonesPerStripe)/_zonesPerDegree;
        return d <= -90.0 ? -90.0 : d;
    }

    double getStripeDecMax(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        double d = static_cast<double>((stripeId + 1)*_zonesPerStripe)/_zonesPerDegree;
        return d >= 90.0 ? 90.0 : d;
    }

    int getStripeZoneMin(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        return stripeId*_zonesPerStripe;
    }

    int getStripeZoneMax(int const stripeId) const {
        assert(stripeId >= _minStripe && stripeId <= _maxStripe && "stripe id out of bounds");
        return stripeId*_zonesPerStripe + _zonesPerStripe - 1;
    }

    int getZonesPerStripe() const { return _zonesPerStripe; }

    int getMaxEntriesPerZoneEstimate() const {
        return _maxEntriesPerZoneEstimate;
    }

    void swap(ZoneStripeChunkDecomposition & zsc);

private :

    std::vector<int> _chunksPerStripe;

    double _zonesPerDegree;
    double _maxZoneAsDouble;
    int _zonesPerStripe;
    int _maxEntriesPerZoneEstimate;
    int _minZone;
    int _maxZone;
    int _minStripe;
    int _maxStripe;
};

inline void swap(ZoneStripeChunkDecomposition & a, ZoneStripeChunkDecomposition & b) {
    a.swap(b);
}


double alpha(double const theta, double const centerDec, double const dec);

double maxAlpha(double const theta, double const centerDec);

void computeChunkIds(
    std::vector<int>                   & chunkIds,
    CircularRegion               const & region,
    ZoneStripeChunkDecomposition const & zsc,
    int                          const   workerId    = 0,
    int                          const   numWorkers  = 1
);

void computeChunkIds(
    std::vector<int>                   & chunkIds,
    RectangularRegion            const & region,
    ZoneStripeChunkDecomposition const & zsc,
    int                          const   workerId    = 0,
    int                          const   numWorkers  = 1
);

inline boost::uint32_t raToScaledInteger(double const ra) {
    assert(ra >= 0.0 && ra < 360.0);
    double d = std::floor(ra*RA_DEC_SCALE);
    return d >= 4294967295.0 ? 4294967295u : static_cast<uint32_t>(d);
}

inline boost::uint32_t deltaRaToScaledInteger(double const delta) {
    assert(delta >= 0.0 && delta <= 180.0);
    return static_cast<int32_t>(std::ceil(delta*RA_DEC_SCALE));
}

inline boost::int32_t decToScaledInteger(double const dec) {
    assert(dec >= -90.0 && dec <= 90.0);
    return static_cast<int32_t>(std::floor(dec*RA_DEC_SCALE));
}

inline boost::int32_t deltaDecToScaledInteger(double const delta) {
    assert(delta >= 0.0 && delta <= 90.0);
    return static_cast<int32_t>(std::ceil(delta*RA_DEC_SCALE));
}

inline double clampDec(double const dec) {
    if (dec <= -90.0) {
        return -90.0;
    } else if (dec >= 90.0) {
        return 90.0;
    }
    return dec;
}


}} // end of namespace lsst::ap

#endif // LSST_AP_SPATIAL_UTIL_H