helper.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: helper functions
 * Author: Joseph Meyer
 * Date: September, 2025
 *
 ******************************************************************************/

 
#pragma once
 
#include <iostream>
#include <filesystem>
#include <mutex>
 
#include <boost/unordered/unordered_flat_map.hpp>
#include <xoshiro.h>
#include <gdal_priv.h>
#include <ogrsf_frmts.h>
#include <ogr_core.h>
 
#define MAXINT8     127
#define MAXINT16    32767
 
namespace sgs {
namespace helper {

struct Index {
    int x = -1;
    int y = -1;
};

struct Field {
    std::string fname; //field name
    int fval; //field value
};

struct RasterBandMetaData {
    GDALRasterBand *p_band = nullptr;
    void *p_buffer = nullptr;
    GDALDataType type = GDT_Unknown;
    size_t size = 0;
    std::string name = "";
    double nan = -1;
    int xBlockSize = -1;
    int yBlockSize = -1;
    std::mutex *p_mutex = nullptr;
};

struct VRTBandDatasetInfo {
    GDALDataset *p_dataset = nullptr;
    std::string filename = "";
 
};

inline void
setStratBandTypeAndSize(
    size_t maxStrata,
    GDALDataType *p_type,
    size_t *p_size
) {
    if (maxStrata <= MAXINT8) {
        *p_type = GDT_Int8;
        *p_size = sizeof(int8_t);
    }
    else if (maxStrata <= MAXINT16) {
        *p_type = GDT_Int16;
        *p_size = sizeof(int16_t);
    }
    else {
        *p_type = GDT_Int32;
        *p_size = sizeof(int32_t);
    }
}

template <typename T>
inline T
getPixelValueDependingOnType(
    GDALDataType type,
    void *p_data,
    size_t index
) {
    switch (type) {
        case GDT_Int8:
            return static_cast<T>(((int8_t *)p_data)[index]);
        case GDT_UInt16:
            return static_cast<T>(((uint16_t *)p_data)[index]);
        case GDT_Int16:
            return static_cast<T>(((int16_t *)p_data)[index]);
        case GDT_UInt32:
            return static_cast<T>(((uint32_t *)p_data)[index]);
        case GDT_Int32:
            return static_cast<T>(((int32_t *)p_data)[index]);
        case GDT_Float32:
            return static_cast<T>(((float *)p_data)[index]);
        case GDT_Float64:
            return static_cast<T>(((double *)p_data)[index]);
        default:
            throw std::runtime_error("raster pixel data type not supported.");
    }
}

inline void
setStrataPixelDependingOnType(
    GDALDataType type,
    void *p_data,
    size_t index,
    bool isNan,
    size_t strata
) {
    switch(type) {
        case GDT_Int8:
            reinterpret_cast<int8_t *>(p_data)[index] = isNan ?
                static_cast<int8_t>(-1) :
                static_cast<int8_t>(strata);
            break;
        case GDT_Int16: 
            reinterpret_cast<int16_t *>(p_data)[index] = isNan ?
                static_cast<int16_t>(-1) :
                static_cast<int16_t>(strata);
            break;
        case GDT_Int32: 
            reinterpret_cast<int32_t *>(p_data)[index] = isNan ?
                static_cast<int32_t>(-1) :
                static_cast<int32_t>(strata);
            break;
        default:
            throw std::runtime_error("strata pixel data type not supported.");
    }
}

inline void
printTypeWarningsForInt32Conversion(GDALDataType type) {
    switch(type) {
        case GDT_UInt32:
            std::cout << "**warning** the pixel type of one of the bands given is an unsigned 32 bit integer. This may result in undefined behavior if the value is not castable to a 32-bit signed integer type." << std::endl;
            break;
        case GDT_Float32:
            std::cout << "**warning** the pixel type of one of the bands given is a 32 bit floating point value. This may result in undefined behavior if the value is not castable to a 32-bit signed integer type." << std::endl;
            break;
        case GDT_Float64:
            std::cout << "**warning** the pixel type of one of the bands given is a 64 bit floating point value. This may result in undefined behavior if the value is not castable to a 32-bit signed integer type." << std::endl;
        default:
            //don't care if converting the type to int32 won't result in overflow problems
            break;
    }
}

inline GDALDataset *
createVirtualDataset(
    std::string driverName, 
    int width, 
    int height, 
    double *geotransform, 
    std::string projection) 
{
    GDALDriver *p_driver = GetGDALDriverManager()->GetDriverByName(driverName.c_str());
    if (!p_driver) {
        throw std::runtime_error("unable to find dataset driver.");
    }
 
    GDALDataset *p_dataset = p_driver->Create(
        "",
        width,
        height,
        0,
        GDT_Unknown,
        nullptr
    );
    if (!p_dataset) {
        throw std::runtime_error("unable to create dataset with driver.");
    }
 
    CPLErr err = p_dataset->SetGeoTransform(geotransform);
    if (err) {
        throw std::runtime_error("error setting geotransform.");
    }
 
    err = p_dataset->SetProjection(projection.c_str());
    if (err) {
        throw std::runtime_error("error setting projection.");
    }
 
    return p_dataset;
}

inline GDALDataset *
createDataset(
    std::string filename,
    std::string driverName, 
    int width, 
    int height, 
    double *geotransform, 
    std::string projection,
    RasterBandMetaData *bands,
    size_t bandCount,
    bool useTiles,
    std::map<std::string, std::string>& driverOptions) 
{
    GDALDriver *p_driver = GetGDALDriverManager()->GetDriverByName(driverName.c_str());
    if (!p_driver) {
        throw std::runtime_error("unable to find dataset driver.");
    }
 
    char **papszOptions = nullptr;
    if (useTiles) {
        if (driverOptions.find("TILED") == driverOptions.end()) {
            papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES");
        }
        if (driverOptions.find("BLOCKXSIZE") == driverOptions.end()) {
            const char *xBlockSizeOption = std::to_string(bands[0].xBlockSize).c_str();
            papszOptions = CSLSetNameValue(papszOptions, "BLOCKXSIZE", xBlockSizeOption);
        }
        if (driverOptions.find("BLOCKYSIZE") == driverOptions.end()) {
            const char *yBlockSizeOption = std::to_string(bands[0].yBlockSize).c_str();
            papszOptions = CSLSetNameValue(papszOptions, "BLOCKYSIZE", yBlockSizeOption);
        }
    }
 
    for (auto const& [key, val] : driverOptions) {
        papszOptions = CSLSetNameValue(papszOptions, key.c_str(), val.c_str());
    }
 
    GDALDataset *p_dataset = p_driver->Create(
        filename.c_str(),
        width,
        height,
        bandCount,
        bands[0].type,
        papszOptions
    );
    CSLDestroy(papszOptions);
 
    if (!p_dataset) {
        throw std::runtime_error("unable to create dataset with driver.");
    }
 
    CPLErr err = p_dataset->SetGeoTransform(geotransform);
    if (err) {
        throw std::runtime_error("error setting geotransform.");
    }
 
    err = p_dataset->SetProjection(projection.c_str());
    if (err) {
        throw std::runtime_error("error setting projection.");
    }
 
    for (size_t band = 0; band < bandCount; band++) {
        GDALRasterBand *p_band = p_dataset->GetRasterBand(band + 1);
        p_band->SetDescription(bands[band].name.c_str());
        p_band->SetNoDataValue(bands[band].nan);
        p_band->GetBlockSize(&bands[band].xBlockSize, &bands[band].yBlockSize);
        bands[band].p_band = p_band;
    }
 
    return p_dataset;
}
 

inline void
addBandToMEMDataset(
    GDALDataset *p_dataset,
    RasterBandMetaData& band)
{
    //allocate data buffer if it has not been allocated yet
    if (!band.p_buffer) {
        band.p_buffer = VSIMalloc3(
            p_dataset->GetRasterXSize(),
            p_dataset->GetRasterYSize(),
            band.size
        );
    }
 
    CPLErr err;
    char **papszOptions = nullptr;
    std::string datapointer = std::to_string((size_t)band.p_buffer);
    papszOptions = CSLSetNameValue(papszOptions, "DATAPOINTER", datapointer.c_str());
 
    err = p_dataset->AddBand(band.type, papszOptions);
    CSLDestroy(papszOptions);
    if (err) {
        throw std::runtime_error("unable to add band to dataset.");
    }
    
    GDALRasterBand *p_band = p_dataset->GetRasterBand(p_dataset->GetRasterCount());
    p_band->SetNoDataValue(band.nan);
    p_band->SetDescription(band.name.c_str());
    band.p_band = p_band;
}

inline void
createVRTBandDataset(
    GDALDataset *p_dataset,
    RasterBandMetaData& band,
    std::string tempFolder,
    std::string key,
    std::vector<VRTBandDatasetInfo>& VRTBandInfo,
    std::map<std::string, std::string>& driverOptions)
{
    std::filesystem::path tmpPath = tempFolder;
    std::filesystem::path tmpName = "strat_breaks_" + key + ".tif";
    tmpPath = tmpPath / tmpName;
 
    VRTBandDatasetInfo info;
    info.filename = tmpPath.string();
    
    bool useTiles = band.xBlockSize != p_dataset->GetRasterXSize() && 
            band.yBlockSize != p_dataset->GetRasterYSize();
 
    double geotransform[6];
    CPLErr err = p_dataset->GetGeoTransform(geotransform);
    if (err) {
        throw std::runtime_error("unable to get geotransform from dataset.");
    }
 
    info.p_dataset = createDataset(
        info.filename,
        "GTiff",
        p_dataset->GetRasterXSize(),
        p_dataset->GetRasterYSize(),
        geotransform,
        std::string(p_dataset->GetProjectionRef()),
        &band,
        1,
        useTiles,
        driverOptions
    );
 
    GDALRasterBand *p_band = info.p_dataset->GetRasterBand(1);
    p_band->SetDescription(band.name.c_str());
    p_band->SetNoDataValue(band.nan);
    p_band->GetBlockSize(&band.xBlockSize, &band.yBlockSize);
    band.p_band = p_band;
        
    VRTBandInfo.push_back(info);    
}
 

inline void
addBandToVRTDataset(
    GDALDataset *p_dataset,
    RasterBandMetaData& band,
    VRTBandDatasetInfo& info
) {
    char **papszOptions = nullptr;
    const char *filename = info.filename.c_str();
    papszOptions = CSLSetNameValue(papszOptions, "subclass", "VRTRawRasterBand");
    papszOptions = CSLSetNameValue(papszOptions, "SourceFilename", filename);
 
    CPLErr err = p_dataset->AddBand(band.type, papszOptions);
    CSLDestroy(papszOptions);
    if (err) {
        throw std::runtime_error("unable to add band to dataset.");
    }
 
    GDALRasterBand *p_VRTBand = p_dataset->GetRasterBand(p_dataset->GetRasterCount());
    p_VRTBand->SetDescription(band.name.c_str());
    p_VRTBand->SetNoDataValue(band.nan);
}

inline void
rasterBandIO(
    RasterBandMetaData& band,
    void *p_buffer,
    int xBlockSize,
    int yBlockSize,
    int xBlock,
    int yBlock,
    int xValid,
    int yValid,
    bool read,
    bool threaded = true)
{
    CPLErr err;
    bool useBlock = xBlockSize == band.xBlockSize &&
            yBlockSize == band.yBlockSize;
 
    if (threaded) {
        band.p_mutex->lock();
    }
    if (useBlock && read) {
        err = read ?
            band.p_band->ReadBlock(xBlock, yBlock, p_buffer) :
            band.p_band->WriteBlock(xBlock, yBlock, p_buffer);
    }
    else {
        err = band.p_band->RasterIO(
            read ? GF_Read : GF_Write,
            xBlock * xBlockSize,
            yBlock * yBlockSize,
            xValid,
            yValid,
            p_buffer,
            xBlockSize,
            yBlockSize,
            band.type,
            0,
            0
        );
    }
    if (threaded) {
        band.p_mutex->unlock();
    }
    if (err) {
        throw read ?
            std::runtime_error("unable to read block from raster.") :
            std::runtime_error("unable to write block to raster.");
    }
}

inline void
addPoint(OGRPoint *p_point, OGRLayer *p_layer) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}

inline void
addPoint(const OGRPoint *p_point, OGRLayer *p_layer) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}

inline void
addPoint(OGRPoint *p_point, OGRLayer *p_layer, Field *p_field) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    p_feature->SetField(p_field->fname.c_str(), p_field->fval);
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}

inline void
addPoint(const OGRPoint *p_point, OGRLayer *p_layer, Field *p_field) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    p_feature->SetField(p_field->fname.c_str(), p_field->fval);
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}

inline void
addPoint(OGRPoint *p_point, OGRLayer *p_layer, std::vector<Field *> *p_fields) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    for (size_t i = 0; i < p_fields->size(); i++) {
        p_feature->SetField((*p_fields)[i]->fname.c_str(), (*p_fields)[i]->fval);
    }
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}

inline void
addPoint(const OGRPoint *p_point, OGRLayer *p_layer, std::vector<Field *> *p_fields) {
    OGRFeature *p_feature = OGRFeature::CreateFeature(p_layer->GetLayerDefn());
    p_feature->SetGeometry(p_point);
    for (size_t i = 0; i < p_fields->size(); i++) {
        p_feature->SetField((*p_fields)[i]->fname.c_str(), (*p_fields)[i]->fval);
    }
    p_layer->CreateFeature(p_feature);
    OGRFeature::DestroyFeature(p_feature);  
}
 

template <typename T>
inline T
point2index(double xCoord, double yCoord, double *IGT, T width) {
    T x = static_cast<T>(IGT[0] + xCoord * IGT[1] + yCoord * IGT[2]);
    T y = static_cast<T>(IGT[3] + xCoord * IGT[4] + yCoord * IGT[5]);
 
    T index = y * width + x;
    return index;
}

class Variance {
    private:
    int64_t k = 0;  //count
    double M = 0;   //running mean
    double S = 0;   //sum of squares
    double oldM = 0;
    
    public:
    inline void
    update(double x) {
        k++;
        oldM = M;
 
        //update running mean
        M = M + (x - M) / static_cast<double>(k);
 
        //update sum of squares
        S = S + (x - M) * (x - oldM);
    }

    inline double 
    getMean() {
        return M;
    }

    inline double 
    getStdev() {
        double variance = S / static_cast<double>(k);
        return std::sqrt(variance);
    }

    inline int64_t
    getCount() {
        return this->k;
    }
};

inline uint64_t
getProbabilityMultiplier(double width, double height, double pixelWidth, double pixelHeight, int startMult, int numSamples, bool useMindist, double accessibleArea) {
    double numer = static_cast<double>(numSamples * startMult * (useMindist ? 3 : 1));
    double denom = height * width;
 
    if (accessibleArea != -1) {
        double totalArea = width * pixelWidth * height * pixelHeight;
        numer *= (totalArea / accessibleArea);
    }
 
    if (numer > denom) {
        return 0;
    }
 
    uint8_t bits = static_cast<uint8_t>(std::ceil(std::log2(denom) - std::log2(numer)));
    return (1 << bits) - 1;
}

class RandValController {
private:
    std::vector<bool> randVals;
    size_t randValIndex = 0;
    uint64_t multiplier = 0;
    xso::xoshiro_4x64_plus *p_rng = nullptr;
    bool alwaysTrue = false;
 
public:
    RandValController(int xBlockSize, int yBlockSize, uint64_t multiplier, xso::xoshiro_4x64_plus *p_rng) {
        if (multiplier == 0) {
            this->alwaysTrue = true;
        }
        else {
            this->randVals.resize(xBlockSize * yBlockSize);
            this->randValIndex = static_cast<size_t>(xBlockSize * yBlockSize);
            this->multiplier = multiplier;
            this->p_rng = p_rng;
        }
    }

    inline void 
    calculateRandValues(void) {
        if (alwaysTrue) {
            return;
        }
 
        for (size_t i = 0; i < randValIndex; i++) {
            randVals[i] = (((*p_rng)() >> 11) & multiplier) == multiplier;
        }
        randValIndex = 0;
    }

    inline bool 
    next(void) {
        if (alwaysTrue) {
            return true;
        }
 
        bool retval = randVals[randValIndex];
        randValIndex++;
        return retval;
    }
};

struct PointHash {
    inline std::size_t operator()(const std::pair<int, int> &v) const {
        std::size_t h1 = std::hash<int>{}(v.first);
        std::size_t h2 = std::hash<int>{}(v.second);
        return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2));
    }
};

typedef boost::unordered::unordered_flat_map<std::pair<int, int>, std::vector<std::pair<double, double>>, PointHash> NeighborMap;

inline bool is_valid_sample(double x, double y, NeighborMap& neighbor_map, float mindist, float mindist_sq) {
    int cx = static_cast<int>(std::floor(x / mindist));
    int cy = static_cast<int>(std::floor(y / mindist));
 
    for (int dx = -1; dx <= 1; dx++) {
        for (int dy = -1; dy <= 1; dy++) {
            std::pair<int, int> neighbor = {cx + dx, cy + dy};
 
            auto it = neighbor_map.find(neighbor);
            if (it == neighbor_map.end())
                continue;
 
            for (const auto &[nx, ny] : it->second) {
                float dxp = x - nx;
                float dyp = y - ny;
                float dist_sq = dxp * dxp + dyp * dyp;
 
                if (dist_sq < mindist_sq) {
                    return false;
                }
            }
        }
    }
 
    // add point to neighborhood
    neighbor_map[{cx, cy}].emplace_back(x, y);
    return true;
}

 
inline std::pair<double, double> sample_to_point(double *GT, Index &index) {
    double px = index.x + 0.5;
    double py = index.y + 0.5;
    double x = GT[0] + px * GT[1] + py * GT[2];
    double y = GT[3] + px * GT[4] + py * GT[5];
 
    return {x, y};
}

 
inline std::pair<double, double> sample_to_point(double *GT, int xs, int ys) {
    double px = xs + 0.5;
    double py = ys + 0.5;
    double x = GT[0] + px * GT[1] + py * GT[2];
    double y = GT[3] + px * GT[4] + py * GT[5];
 
    return {x, y};
}
 
} //namespace helper
} //namespace sgs