pca.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: C++ implementation of PCA
 * Author: Joseph Meyer
 * Date: October, 2025
 *
 ******************************************************************************/

 
#include <iostream>
#include <numeric>
 
#include "utils/helper.h"
#include "utils/raster.h"
 
#include "oneapi/dal.hpp"
 
namespace sgs {
namespace pca {
 
typedef oneapi::dal::homogen_table              DALHomogenTable;

template <typename T>
struct PCAResult {
    std::vector<std::vector<T>> eigenvectors;
    std::vector<T> eigenvalues;
    std::vector<double> means;
    std::vector<double> stdevs;
};

template <typename T>
PCAResult<T>
calculatePCA(
    std::vector<helper::RasterBandMetaData>& bands,
    GDALDataType type,
    size_t size,
    int width,
    int height,
    int nComp)
{
    int bandCount = static_cast<int>(bands.size());
    T *p_data = reinterpret_cast<T *>(VSIMalloc3(width * height, bandCount, size));
 
    std::vector<helper::Variance> bandVariances(bandCount);
    std::vector<T> noDataVals(bandCount);
    for (size_t i = 0; i < bands.size(); i++) {
        noDataVals[i] = static_cast<T>(bands[i].nan);
    }
 
    //read full bands into p_data
    for (size_t i = 0; i < bands.size(); i++) {
        bands[i].p_band->RasterIO(
            GF_Read,
            0,
            0,
            width,
            height,
            (void *)((size_t)p_data + i * size),    
            width,
            height,
            type,
            size * bands.size(),
            size * bands.size() * width
        );  
    }
 
    //remove nodata values from p_data, ensuring data pixels are consecutive
    int nFeatures = 0;
    for (int i = 0; i < height * width; i++) {
        bool isNan = false;
        for (int b = 0; b < bandCount; b++) { 
            T val = p_data[i * bandCount + b];
            isNan = val == noDataVals[b] || std::isnan(val);
            if (isNan) {
                break;
            }
            p_data[nFeatures * bandCount + b] = val;
        }
        nFeatures += !isNan;
    }
 
    //update variance calculations
    for (int i = 0; i < nFeatures; i++) {
        for (int b = 0; b < bandCount; b++) {
            T val = p_data[i * bandCount + b];
            bandVariances[b].update(static_cast<double>(val));  
        }
    }
 
    //calculate pca
    DALHomogenTable table = DALHomogenTable::wrap<T>(p_data, nFeatures, bandCount, oneapi::dal::data_layout::row_major);
    const auto desc = oneapi::dal::pca::descriptor<T, oneapi::dal::pca::method::cov>().set_component_count(nComp).set_deterministic(true);
        const auto result = oneapi::dal::train(desc, table);
 
    VSIFree(p_data);
 
    PCAResult<T> retval;
    auto eigenvectors = result.get_eigenvectors();
    auto eigenvalues = result.get_eigenvalues();
    int64_t eigRows = eigenvectors.get_row_count();
    int64_t eigCols = eigenvectors.get_column_count();
 
    oneapi::dal::row_accessor<const float> eigVecAcc {eigenvectors};
    auto eigVecBlock = eigVecAcc.pull({0, eigRows});
 
    oneapi::dal::row_accessor<const float> eigValAcc {eigenvalues};
    auto eigValBlock = eigValAcc.pull({0, 1});
 
    retval.eigenvectors.resize(eigRows);
    retval.eigenvalues.resize(eigRows);
    for (int64_t i = 0; i < eigRows; i++) {
        retval.eigenvalues[i] = static_cast<double>(eigValBlock[i]);
 
        retval.eigenvectors[i].resize(eigCols);
        for (int64_t j = 0; j < eigCols; j++) {
            retval.eigenvectors[i][j] = static_cast<double>(eigVecBlock[i * eigCols + j]);
        }
    }
 
    for (int b = 0; b < bandCount; b++) {
        retval.means.push_back(bandVariances[b].getMean());
        retval.stdevs.push_back(bandVariances[b].getStdev());
    }
 
    return retval;  
}

template <typename T>
PCAResult<T>
calculatePCA(
    std::vector<helper::RasterBandMetaData>& bands,
    GDALDataType type,
    size_t size,
    int xBlockSize,
    int yBlockSize,
    int xBlocks,
    int yBlocks,
    int nComp)
{
    int bandCount = static_cast<int>(bands.size());
    T *p_data = reinterpret_cast<T *>(VSIMalloc3(xBlockSize * yBlockSize, bandCount, size));
 
    std::vector<helper::Variance> bandVariances(bandCount);
    std::vector<T> noDataVals(bandCount);
    for (size_t i = 0; i < bands.size(); i++) {
        noDataVals[i] = static_cast<T>(bands[i].nan);
    }
 
    const auto desc = oneapi::dal::pca::descriptor<T, oneapi::dal::pca::method::cov>().set_component_count(nComp).set_deterministic(true);
    oneapi::dal::pca::partial_train_result<> partial_result;
 
    for (int yBlock = 0; yBlock < yBlocks; yBlock++) {
        for (int xBlock = 0; xBlock < xBlocks; xBlock++) {
            int xValid, yValid;
            bands[0].p_band->GetActualBlockSize(xBlock, yBlock, &xValid, &yValid);
        
            //read bands into memory    
            for (int i = 0; i < bandCount; i++) {
                bands[i].p_band->RasterIO(
                    GF_Read,
                    xBlock * xBlockSize,
                    yBlock * yBlockSize,
                    xValid,
                    yValid,
                    (void *)((size_t)p_data + i * size),
                    xValid,
                    yValid,
                    type,
                    size * bands.size(),
                    size * bands.size() * static_cast<size_t>(xBlockSize)
                );              
            }
 
            //remove nodata values
            int nFeatures = 0;
            for (int x = 0; x < xValid; x++) {
                for (int y = 0; y < yValid; y++) {
                    bool isNan = false;
                    for (int b = 0; b < bandCount; b++) {
                        T val = p_data[((y * xBlockSize) + x) * bandCount + b];
                        isNan = std::isnan(val) || val == noDataVals[b];
                        if (isNan) {
                            break;
                        }
                        p_data[nFeatures * bandCount + b] = val;
                    }
                    nFeatures += !isNan;
                }
            }
 
            //update variance calculations
            for (int i = 0; i < nFeatures; i++) {
                for (int b = 0; b < bandCount; b++) {
                    T val = p_data[i * bandCount + b];
                    bandVariances[b].update(static_cast<double>(val));  
                }
            }
 
            //calculate partial result
            DALHomogenTable table = DALHomogenTable(p_data, nFeatures, bandCount, [](const T*){}, oneapi::dal::data_layout::row_major);
            partial_result = oneapi::dal::partial_train(desc, partial_result, table);
        }
    }
 
    auto result = oneapi::dal::finalize_train(desc, partial_result);
    
    VSIFree(p_data);
 
    PCAResult<T> retval;
    auto eigenvectors = result.get_eigenvectors();
    auto eigenvalues = result.get_eigenvalues();
    int64_t eigRows = eigenvectors.get_row_count();
    int64_t eigCols = eigenvectors.get_column_count();
 
    oneapi::dal::row_accessor<const T> eigVecAcc {eigenvectors};
    auto eigVecBlock = eigVecAcc.pull({0, eigRows});
 
    oneapi::dal::row_accessor<const T> eigValAcc {eigenvalues};
    auto eigValBlock = eigValAcc.pull({0, 1});
 
    retval.eigenvectors.resize(eigRows);
    retval.eigenvalues.resize(eigRows);
    for (int64_t i = 0; i < eigRows; i++) {
        retval.eigenvalues[i] = static_cast<T>(eigValBlock[i]);
 
        retval.eigenvectors[i].resize(eigCols);
        for (int64_t j = 0; j < eigCols; j++) {
            retval.eigenvectors[i][j] = static_cast<T>(eigValBlock[i * eigCols + j]);
        }
    }
    
    for (int b = 0; b < bandCount; b++) {
        retval.means.push_back(bandVariances[b].getMean());
        retval.stdevs.push_back(bandVariances[b].getStdev());
    }
 
    return retval;  
}

template <typename T>
void 
writePCA(
    std::vector<helper::RasterBandMetaData>& bands,
    std::vector<helper::RasterBandMetaData>& PCABands,
    PCAResult<T>& result,
    GDALDataType type,
    size_t size,
    int height,
    int width)
{
    int bandCount = static_cast<int>(bands.size());
    int nComp = static_cast<int>(PCABands.size());
    std::vector<T> noDataVals(bandCount); 
    T resultNan = std::nan("");
    
    //set no data values and read input bands
    T *p_data = reinterpret_cast<T *>(VSIMalloc3(bandCount, height * width, size));
    for (int b = 0; b < bandCount; b++) {
        noDataVals[b] = static_cast<T>(bands[b].nan);
 
        bands[b].p_band->RasterIO(
            GF_Read,
            0,
            0,
            width,
            height,
            (void *)((size_t)p_data + b * size),    
            width,
            height,
            type,
            size * bandCount,
            size * bandCount * width
        );  
    }
 
    //scale and shift data pixels, set no data pixels to nan
    for (int i = 0; i < height * width; i++) {
        int bi = i * bandCount;
 
        for (int b = 0; b < bandCount; b++) {
            T val = p_data[bi + b];
            p_data[bi + b] = (val == noDataVals[b]) ?
                resultNan :
                (val - result.means[b]) / result.stdevs[b];
        }
    }
 
    //read result eigenvectors into matrix format
    T *p_comp = reinterpret_cast<T *>(VSIMalloc3(nComp, bandCount, size));
    for (int c = 0; c < nComp; c++) {
        int ci = c * bandCount;
        for (int b = 0; b < bandCount; b++) {
            p_comp[ci + b] = result.eigenvectors[c][b];
        }
    }

 
    //create DAL homogen table wrappers for input data  
    const auto dataTable = DALHomogenTable(p_data, height * width, bandCount, [](const T*){}, oneapi::dal::data_layout::row_major);
    const auto compTable = DALHomogenTable(p_comp, nComp, bandCount, [](const T*){}, oneapi::dal::data_layout::row_major); 
    
    //create DAL descriptor and compute the result
    const auto kernel_desc = oneapi::dal::linear_kernel::descriptor{}.set_scale(1.0).set_shift(0.0);
    const auto compute_result = oneapi::dal::compute(kernel_desc, dataTable, compTable);    
 
    //get the raw data from the result
    const auto values = compute_result.get_values();
    oneapi::dal::row_accessor<const T> valAcc {values};
    const T *p_result = valAcc.pull({0, height * width}).get_data();
 
    //write the raw result to the output
    for (int c = 0; c < nComp; c++) {
        PCABands[c].p_band->RasterIO(
            GF_Write,
            0,
            0,
            width,
            height,
            (void *)((size_t)p_result + c * size),
            width,
            height,
            type,
            size * nComp,
            size * nComp * width
        );
    }
 
    VSIFree(p_data);
    VSIFree(p_comp);
}

template <typename T>
void 
writePCA(
    std::vector<helper::RasterBandMetaData>& bands,
    std::vector<helper::RasterBandMetaData>& PCABands,
    PCAResult<T>& result,
    GDALDataType type,
    size_t size,
    int xBlockSize,
    int yBlockSize,
    int xBlocks,
    int yBlocks)
{
    int bandCount = static_cast<int>(bands.size());
    int nComp = static_cast<int>(PCABands.size());
    std::vector<T> noDataVals(bandCount);
    T resultNan = std::nan("");
    for (int i = 0; i < bandCount; i++) {
        noDataVals[i] = static_cast<T>(bands[i].nan);
    }
 
    T *p_data = reinterpret_cast<T *>(VSIMalloc3(xBlockSize * yBlockSize, size, bandCount));
    T *p_comp = reinterpret_cast<T *>(VSIMalloc3(nComp, bandCount, size));
 
    //create DAL homogen table wrappers for input data
    const auto dataTable = DALHomogenTable(p_data, xBlockSize * yBlockSize, bandCount, [](const T*){}, oneapi::dal::data_layout::row_major);
    const auto compTable = DALHomogenTable(p_comp, nComp, bandCount, [](const T*){}, oneapi::dal::data_layout::row_major);
 
    for (int yBlock = 0; yBlock < yBlocks; yBlock++) {
        for (int xBlock = 0; xBlock < xBlocks; xBlock++) {
            int xValid, yValid;
            bands[0].p_band->GetActualBlockSize(xBlock, yBlock, &xValid, &yValid);
        
            //read bands into memory, specifically into the memory location of the dal data homogen table
            for (int b = 0; b < bandCount; b++) {
                bands[b].p_band->RasterIO(
                    GF_Read,
                    xBlock * xBlockSize,
                    yBlock * yBlockSize,
                    xValid,
                    yValid,
                    (void *)((size_t)p_data + b * size),
                    xValid,
                    yValid,
                    type,
                    size * bands.size(),
                    size * bands.size() * static_cast<size_t>(xBlockSize)
                );              
            }
 
            //scale and shift data pixels, set no data pixels to nan
            for (int i = 0; i < xBlockSize * yBlockSize; i++) {
                int bi = i * bandCount;
 
                for (int b = 0; b < bandCount; b++) {
                    T val = p_data[bi + b];
                    p_data[bi + b] = (val == noDataVals[b]) ?
                        resultNan :
                        (val - result.means[b]) / result.stdevs[b];
                }
            }

 
            //create DAL descriptor and compute the result
            const auto kernel_desc = oneapi::dal::linear_kernel::descriptor{}.set_scale(1.0).set_shift(0.0);
            const auto compute_result = oneapi::dal::compute(kernel_desc, dataTable, compTable);
 
            //get the raw data from the result
            const auto values = compute_result.get_values();
            oneapi::dal::row_accessor<const T> valAcc {values};
            const T *p_result = valAcc.pull({0, xBlockSize * yBlockSize}).get_data();
 
            //write the raw result to the output
            for (int c = 0; c < nComp; c++) {
                PCABands[c].p_band->RasterIO(
                    GF_Write,
                    xBlock * xBlockSize,
                    yBlock * yBlockSize,
                    xValid,
                    yValid,
                    (void *)((size_t)p_comp + c * size),
                    xValid,
                    yValid,
                    type,
                    size * nComp,
                    size * nComp * xBlockSize
                );
            }
        }
    }
 
    VSIFree(p_data);
    VSIFree(p_comp);
}

std::tuple<
    raster::GDALRasterWrapper *, 
    std::vector<std::vector<double>>, 
    std::vector<double>, 
    std::vector<double>, 
    std::vector<double>
>
pca(
    raster::GDALRasterWrapper *p_raster,
    int nComp,
    bool largeRaster,
    std::string tempFolder,
    std::string filename,
    std::map<std::string, std::string> driverOptions)
{
    GDALAllRegister();
 
    int bandCount = p_raster->getBandCount();
    int height = p_raster->getHeight();
    int width = p_raster->getWidth();
    double *geotransform = p_raster->getGeotransform();
    std::string projection = std::string(p_raster->getDataset()->GetProjectionRef());
    
    bool isMEMDataset = !largeRaster && filename == "";
    bool isVRTDataset = largeRaster && filename == "";
    GDALDataset *p_dataset = nullptr;
    
    std::vector<helper::RasterBandMetaData> bands(bandCount);
    std::vector<helper::RasterBandMetaData> pcaBands(nComp);
    std::vector<helper::VRTBandDatasetInfo> VRTBandInfo;
 
    int xBlockSize, yBlockSize;
    p_raster->getRasterBand(0)->GetBlockSize(&xBlockSize, &yBlockSize);
 
    GDALDataType type = GDT_Float32;
    size_t size = sizeof(float);
    for (int i = 0; i < bandCount; i++) {
        bands[i].p_band = p_raster->getRasterBand(i);
        bands[i].nan = bands[i].p_band->GetNoDataValue();
 
        if (p_raster->getRasterBandType(i) == GDT_Float64) {
            type = GDT_Float64;
            size = sizeof(double);
        }
    }
 
    if (isMEMDataset) {
        p_dataset = helper::createVirtualDataset("MEM", width, height, geotransform, projection);
    
        for (int i = 0; i < nComp; i++) {
            pcaBands[i].type = type == GDT_Float64 ? GDT_Float64 : GDT_Float32;
            pcaBands[i].size = type == GDT_Float64 ? sizeof(double) : sizeof(float);
            pcaBands[i].name = "comp_" + std::to_string(i + 1);
            pcaBands[i].nan = std::nan("");
            helper::addBandToMEMDataset(p_dataset, pcaBands[i]);
        }
    }
    else if (isVRTDataset){
        p_dataset = helper::createVirtualDataset("VRT", width, height, geotransform, projection);
    
        for (int i = 0; i < nComp; i++) {
            pcaBands[i].type = type == GDT_Float64 ? GDT_Float64 : GDT_Float32;
            pcaBands[i].size = type == GDT_Float64 ? sizeof(double) : sizeof(float);
            pcaBands[i].name = "comp_" + std::to_string(i + 1);
            pcaBands[i].nan = std::nan("");
            helper::createVRTBandDataset(p_dataset, pcaBands[i], tempFolder, pcaBands[i].name, VRTBandInfo, driverOptions);
        }
    }
    else {
        std::filesystem::path filepath = filename;
        std::string extension = filepath.extension().string();
        std::string driver;
 
        if (extension == ".tif") {
            driver = "GTiff";
        }
        else {
            throw std::runtime_error("sgs only supports .tif files right now.");
        }
        
        bool useTiles = xBlockSize != width &&
                yBlockSize != height;
 
        for (int i = 0; i < nComp; i++) {
            pcaBands[i].type = type == GDT_Float64 ? GDT_Float64 : GDT_Float32;
            pcaBands[i].size = type == GDT_Float64 ? sizeof(double) : sizeof(float);
            pcaBands[i].name = "comp_" + std::to_string(i + 1);
            pcaBands[i].nan = std::nan("");
        
            if (useTiles) {
                pcaBands[i].xBlockSize = xBlockSize;
                pcaBands[i].yBlockSize = yBlockSize;
            }
        }
 
        p_dataset = helper::createDataset(
            filename,
            driver,
            width,
            height, 
            geotransform,
            projection,
            pcaBands.data(),
            pcaBands.size(),
            useTiles,
            driverOptions
        );
    }
 
    //get block size
    int xBlocks = (width + xBlockSize - 1) / xBlockSize;
    int yBlocks = (height + yBlockSize - 1) / yBlockSize;
 
    //these vectors are not used for calculation, but set and returned to the Python side of the application for reference
    std::vector<std::vector<double>> eigenvectors; 
    std::vector<double> eigenvalues;
    std::vector<double> means;
    std::vector<double> stdevs;
 
    //calculate PCA eigenvectors (and eigenvalues), and write values to PCA bands
    switch(type) {
        case GDT_Float32: {
            PCAResult<float> result;
            if (largeRaster) {
                result = calculatePCA<float>(bands, type, size, xBlockSize, yBlockSize, xBlocks, yBlocks, nComp);
                writePCA<float>(bands, pcaBands, result, type, size, xBlockSize, yBlockSize, xBlocks, yBlocks);
            }
            else {
                result = calculatePCA<float>(bands, type, size, width, height, nComp);
                writePCA<float>(bands, pcaBands, result, type, size, height, width);
            }
 
            eigenvectors.resize(result.eigenvectors.size());
            for (size_t i = 0; i < result.eigenvectors.size(); i++) {
                eigenvectors[i].resize(result.eigenvectors[i].size());
                for (size_t j = 0; j < result.eigenvectors[i].size(); j++) {
                    eigenvectors[i][j] = static_cast<double>(result.eigenvectors[i][j]);
                }
            }
 
            eigenvalues.resize(result.eigenvalues.size());
            for (size_t i = 0; i < result.eigenvalues.size(); i++) {
                eigenvalues[i] = static_cast<double>(result.eigenvalues[i]);
            }
        
            means = result.means;
            stdevs = result.stdevs;
 
            break;
        }
        case GDT_Float64: {
            PCAResult<double> result;
            if (largeRaster) {
                result = calculatePCA<double>(bands, type, size, xBlockSize, yBlockSize, xBlocks, yBlocks, nComp);
                writePCA<double>(bands, pcaBands, result, type, size, xBlockSize, yBlockSize, xBlocks, yBlocks);
            }
            else {
                result = calculatePCA<double>(bands, type, size, width, height, nComp);
                writePCA<double>(bands, pcaBands, result, type, size, height, width);
            }
 
            eigenvectors = result.eigenvectors;
            eigenvalues = result.eigenvalues;
            means = result.means;
            stdevs = result.stdevs;
 
            break;
        }
        default:
            throw std::runtime_error("should not be here! GDALDataType should be one of Float32/Float64!");
    }
 
    if (isVRTDataset) {
        for (int c = 0; c < nComp; c++) {
            GDALClose(VRTBandInfo[c].p_dataset);
            helper::addBandToVRTDataset(p_dataset, pcaBands[c], VRTBandInfo[c]);
        }
    }
 
    std::vector<void *> buffers(nComp);
    if (isMEMDataset) {
        for (int c = 0; c < nComp; c++) {
            buffers[c] = pcaBands[c].p_buffer;
        }
    }
 
    raster::GDALRasterWrapper *p_outrast = !isMEMDataset ?
        new raster::GDALRasterWrapper(p_dataset) :
        new raster::GDALRasterWrapper(p_dataset, buffers);
 
    return {p_outrast, eigenvectors, eigenvalues, means, stdevs};
}
 
} //namespace pca
} //namespace sgs