PragueSkyModel.h

 
// Copyright 2022 Charles University
// SPDX-License-Identifier: Apache-2.0
 
#include <exception>
#include <string>
#include <vector>
#include <cmath>
 
double lerp(const double from, const double to, const double factor);
 
class PragueSkyModel {
 
 
    // Public types
public:
    class DatasetNotFoundException : public std::exception {
    private:
        const std::string message;
 
    public:
        DatasetNotFoundException(const std::string &filename) : message(std::string("Dataset file ") + filename + std::string(" not found")) {
        }
 
        virtual const char *what() const throw() {
            return message.c_str();
        }
    };
 
    class DatasetReadException : public std::exception {
    private:
        const std::string message;
 
    public:
        DatasetReadException(const std::string &parameterName) : message(std::string("Dataset reading failed at ") + parameterName) {
        }
 
        virtual const char *what() const throw() {
            return message.c_str();
        }
    };
 
    class NoPolarisationException : public std::exception {
    private:
        const std::string message;
 
    public:
        NoPolarisationException() : message(std::string("The supplied dataset does not contain polarisation data")) {
        }
 
        virtual const char *what() const throw() {
            return message.c_str();
        }
    };
 
    class NotInitializedException : public std::exception {
    private:
        const std::string message;
 
    public:
        NotInitializedException() : message(std::string("The model is not initialized")) {
        }
 
        virtual const char *what() const throw() {
            return message.c_str();
        }
    };
 
    class Vector3 {
    public:
        double x, y, z;
 
        Vector3() {
            this->x = 0.0;
            this->y = 0.0;
            this->z = 0.0;
        }
 
        Vector3(double x, double y, double z) {
            this->x = x;
            this->y = y;
            this->z = z;
        }
 
        Vector3 operator+(const Vector3 &other) const {
            return Vector3(x + other.x, y + other.y, z + other.z);
        }
        Vector3 operator-(const Vector3 &other) const {
            return Vector3(x - other.x, y - other.y, z - other.z);
        }
 
        Vector3 operator*(const double factor) const {
            return Vector3(x * factor, y * factor, z * factor);
        }
 
        Vector3 operator/(const double factor) const {
            return Vector3(x / factor, y / factor, z / factor);
        }
 
        bool isZero() const {
            return x == 0.0 && y == 0.0 && z == 0.0;
        }
    };
 
    struct Parameters {
        double theta;
 
        double gamma;
 
        double shadow;
 
        double zero;
 
        double elevation;
 
        double altitude;
 
        double visibility;
 
        double albedo;
    };
 
    struct AvailableData {
        double albedoMin;
        double albedoMax;
        double altitudeMin;
        double altitudeMax;
        double elevationMin;
        double elevationMax;
        double visibilityMin;
        double visibilityMax;
        bool polarisation;
        int channels;
        double channelStart;
        double channelWidth;
    };
 
 
    // Private types
private:
    struct InterpolationParameter {
        double factor;
        int index;
    };
 
    struct AngleParameters {
        InterpolationParameter gamma, alpha, zero;
    };
 
    struct ControlParameters {
        std::array<std::vector<float>::const_iterator, 16> coefficients;
        std::array<double, 4> interpolationFactor;
    };
 
    struct Metadata {
        int rank;
 
        int sunOffset;
        int sunStride;
        std::vector<double> sunBreaks;
 
        int zenithOffset;
        int zenithStride;
        std::vector<double> zenithBreaks;
 
        int emphOffset; // not used for polarisation
        std::vector<double> emphBreaks; // not used for polarisation
 
        int totalCoefsSingleConfig;
        int totalCoefsAllConfigs;
    };
 
    struct TransmittanceParameters {
        InterpolationParameter altitude;
        InterpolationParameter distance;
    };
 
 
    // Private data
    int channels;
    double channelStart;
    double channelWidth;
 
    bool initialized;
 
    // Total number of configurations
    int totalConfigs;
 
    // Number of configurations skipped from the beginning of the radiance/polarisation coefficients part
    // of the dataset file (if loading of just one visibility was requested)
    int skippedConfigsBegin;
 
    // Number of configurations skipped till the end of the radiance/polarisation coefficients part
    // of the dataset file (if loading of just one visibility was requested)
    int skippedConfigsEnd;
 
    // Metadata common for radiance and polarisation
 
    std::vector<double> visibilitiesRad;
    std::vector<double> albedosRad;
    std::vector<double> altitudesRad;
    std::vector<double> elevationsRad;
 
    // Radiance metadata
 
    Metadata metadataRad;
 
    // Radiance data
    //
    // Structure:
    // [[[[[[ sunCoefsRad       (sunBreaksCountRad * float),
    //        zenithCoefsRad (zenithBreaksCountRad * float) ] * rankRad,
    //        emphCoefsRad     (emphBreaksCountRad * float) ]
    //  * channels ] * elevationCount ] * altitudeCount ] * albedoCount ] * visibilityCount
 
    std::vector<float> dataRad;
 
    // Polarisation metadata
 
    Metadata metadataPol;
 
    // Polarisation data
    //
    // Struture:
    // [[[[[[ sunCoefsPol       (sunBreaksCountPol * float),
    //        zenithCoefsPol (zenithBreaksCountPol * float) ] * rankPol]
    // * channels ] * elevationCount ] * altitudeCount ] * albedoCount ] * visibilityCount
 
    std::vector<float> dataPol;
 
    // Transmittance metadata
 
    int aDim;
    int dDim;
    int rankTrans;
    std::vector<double> altitudesTrans;
    std::vector<double> visibilitiesTrans;
 
    // Transmittance data
 
    std::vector<float> dataTransU;
    std::vector<float> dataTransV;
 
 
    // Public methods
public:
    PragueSkyModel() : initialized(false) {};
 
    void initialize(const std::string &filename, const double singleVisibility = 0.0);
 
    bool isInitialized() const {
        return initialized;
    }
 
    AvailableData getAvailableData() const;
 
    Parameters computeParameters(const Vector3 &viewPoint, const Vector3 &viewDirection, const double groundLevelSolarElevationAtOrigin, const double groundLevelSolarAzimuthAtOrigin, const double visibility, const double albedo) const;
 
    double skyRadiance(const Parameters &params, const double wavelength) const;
 
    double sunRadiance(const Parameters &params, const double wavelength) const;
 
    double polarisation(const Parameters &params, const double wavelength) const;
 
    double transmittance(const Parameters &params, const double wavelength, const double distance) const;
 
 
    // Private methods
private:
    void readRadiance(FILE *handle, const double singleVisibility);
 
    void readTransmittance(FILE *handle);
 
    void readPolarisation(FILE *handle);
 
    // Sky radiance and polarisation
 
    std::vector<float>::const_iterator getCoefficients(const std::vector<float> &dataset, const int totalCoefsSingleConfig, const int elevation, const int altitude, const int visibility, const int albedo, const int wavelength) const;
 
    template<int TOffset, int TLevel>
    double interpolate(const AngleParameters &angleParameters, const ControlParameters &controlParameters, const Metadata &metadata) const {
        if constexpr (TLevel == 4) {
            return reconstruct(angleParameters, controlParameters.coefficients[TOffset], metadata);
        } else {
            // Compute the first value
            const double resultLow = interpolate<TOffset, TLevel + 1>(angleParameters, controlParameters, metadata);
 
            // Skip the second value if not useful or not available.
            if (controlParameters.interpolationFactor[TLevel] < 1e-6) {
                return resultLow;
            }
 
            // Compute the second value
            const double resultHigh = interpolate<TOffset + (1 << (3 - TLevel)), TLevel + 1>(angleParameters, controlParameters, metadata);
 
            return lerp(resultLow, resultHigh, controlParameters.interpolationFactor[TLevel]);
        }
    }
 
    double reconstruct(const AngleParameters &angleParameters, const std::vector<float>::const_iterator controlParameters, const Metadata &metadata) const;
 
    InterpolationParameter getInterpolationParameter(const double queryVal, const std::vector<double> &breaks) const;
 
    double evaluateModel(const Parameters &params, const double wavelength, const std::vector<float> &data, const Metadata &metadata) const;
 
    // Transmittance
 
    std::vector<float>::const_iterator getCoefficientsTransBase(const int altitude, const int a, const int d) const;
 
    std::vector<float>::const_iterator getCoefficientsTrans(const int visibility, const int altitude, const int wavelength) const;
 
    double interpolateTrans(const int visibilityIndex, const InterpolationParameter altitudeParam, const TransmittanceParameters transParams, const int channelIndex) const;
 
    double reconstructTrans(const int visibilityIndex, const int altitudeIndex, const TransmittanceParameters transParams, const int channelIndex) const;
 
    InterpolationParameter getInterpolationParameterTrans(const double value, const int paramCount, const int power) const;
 
    TransmittanceParameters toTransmittanceParams(const double theta, const double distance, const double altitude) const;
};