Past Summer Project: swissALTI3D

TIF3D_Switzerland.cpp

/*
*  Summer project I worked on to get radio map given dropped antenna tower, height and so on.
*  Calculates fresnel map and draws strength of signal given origin point taking into account
*  the elevation, LOS, and general environment of Switzerland. 
*  TIF Map: https://www.swisstopo.admin.ch/en/height-model-swissalti3d
*/

/* main.cpp */

#include <cstdio>
#include <cstdlib>
#include <errno.h>
#include <vector>
#include <string>
#include <filesystem>
#include <iostream>
#include <bits/ostream.tcc>
#include <gdal/gdal_priv.h>
#include <fstream>
#include "ops_transform.h"


using namespace std;


int             antenna_count = 0;
constexpr float CELL_TOWER_HEIGHT = 10.0f;
constexpr float CELL_FREQUENCY    = 5800000000;
constexpr int   LIGHT_SPEED       = 299792458;


struct XYCoordinate
{
    int x           = 0;
    int y           = 0;
    float elevation = 0;
    float antenna   = 0;
    float height_m  = 0;

    void print() const
    {
        printf("(%d, %d, %0.10f, %.2f, %.10f)\n", x, y, elevation, antenna, height_m);
    }
};



XYCoordinate create_antenna(const float* elevationData, const int x, const int y, float antenna, int nXsize);
int get_rf_strength(const float* elevationData, XYCoordinate alpha, XYCoordinate beta, const int nXsize);
bool is_signal_obstructed(const float actual_clear, const float fresnel_clear);
std::vector<uint8_t> get_rf_map(const float* elevationData, XYCoordinate tx, int width, int height, const float wavelength);
void draw_map(const std::vector<uint8_t>& data, int width, int height, const std::string& filename);
float get_distance(const XYCoordinate alpha, const XYCoordinate beta, const float pixelResolution);
int get_fresnel_strength(const float actual, const float fresnel);


XYCoordinate create_antenna(const float* elevationData, const int x, const int y, float antenna, int nXsize)
{
    XYCoordinate cell;
    cell.x = x, cell.y = y;
    cell.antenna = antenna;
    cell.elevation = elevationData[cell.y * nXsize + cell.x];
    cell.height_m  = cell.elevation + cell.antenna;

    antenna_count++;

    return cell;
}


int get_rf_strength(const float* elevationData, XYCoordinate tx, XYCoordinate rx, const int nXsize, const float wavelength)
{
    auto min       = 1000.0f;
    auto min_fres  = 0.0f;
    auto n_samples = 100;

    for ( int i = 0; i < n_samples; i++ )
    {
        XYCoordinate interpolated;
        auto t = (float)i / (n_samples - 1);


        interpolated.x = tx.x + t * (rx.x - tx.x);
        interpolated.y = tx.y + t * (rx.y - tx.y);


        float ground_elevation = elevationData[interpolated.y * nXsize + interpolated.x];
        float signal_elevation = tx.height_m + t * (rx.height_m - tx.height_m);
        float actual_clear   = signal_elevation - ground_elevation;
        interpolated.height_m = signal_elevation;


        float d_tx_sample    = get_distance(tx, interpolated, 0.5);
        float d_sample_rx    = get_distance(interpolated, rx, 0.5);
        float fresnel_radius = std::sqrt(wavelength * d_tx_sample * d_sample_rx / (d_tx_sample + d_sample_rx));
        float fresnel_clear  = fresnel_radius * 0.6;


        if ( is_signal_obstructed(actual_clear, fresnel_clear) )
        {
            min = actual_clear;
            min_fres = fresnel_clear;
        }

        if ( actual_clear < min )
        {
            min = actual_clear;
            min_fres = fresnel_clear;
        }
    }
    return get_fresnel_strength(min, min_fres);
}

/**
 * @brief checks if signal is obstructed
 * @param actual_clear vertical distance between ground and signal
 * @param fresnel_clear f1 fat beam
 * @return true if obstructed
 */
bool is_signal_obstructed(const float actual_clear, const float fresnel_clear)
{
    if (actual_clear < fresnel_clear)
        return true;

    return false;
}

std::vector<uint8_t> get_rf_map(const float* elevationData, XYCoordinate tx,
                                int width, int height, const float wavelength)
{
    std::vector<uint8_t> rf_map(width * height, 0);;

    for ( int y = 0; y < height; y++ )
    {
        for ( int x = 0; x < width; x++ )
        {
            // edge case: transmitter location
            if ( x == tx.x && y == tx.y ) { rf_map[y * width + x] = 255; continue; }

            XYCoordinate rx;
            rx.x = x;
            rx.y = y;
            rx.elevation = elevationData[y * width + x];
            rx.antenna = 1.5f;
            rx.height_m = rx.elevation + rx.antenna;

            int signal = get_rf_strength(elevationData, tx, rx, width, wavelength);
            rf_map[y * width + x] = static_cast<uint8_t>((signal * 255) / 100);
        }

    }
    return rf_map;
}

/**
 * @brief writes uint8 (0-255) values into a file for viewing
 * @param data      0-255 values
 * @param width     width of each TIF  (2,000)
 * @param height    height of each TIF (2,000)
 * @param filename  filename to write to
 */
void draw_map(const std::vector<uint8_t>& terrain_v,
              const std::vector<uint8_t>& signal_v,
              int width, int height, const std::string& filename)
{
    std::ofstream file(filename, std::ios::binary);
    file << "P6\n" << width << " " << height << "\n255\n";

    for ( int i = 0; i < width * height; i++ )
    {
        uint8_t terrain = terrain_v[i];
        uint8_t signal = signal_v[i];

        uint8_t r,g,b;

        if ( signal == 0 ) { r = g = b = terrain; }
        else
        {
            float alpha = 0.3f;

            if      ( signal < 64  ) { r = 255; g = 0;   b = 0;   } // poor signal
            else if ( signal < 128 ) { r = 255; g = 255; b = 0;   } // okay signal
            else if ( signal < 192 ) { r = 0;   g = 255; b = 0;   } // good signal
            else                     { r = 0;   g = 0;   b = 255; } // exel signal

            r = static_cast<uint8_t>(alpha * r + (1.0f - alpha) * terrain);
            g = static_cast<uint8_t>(alpha * g + (1.0f - alpha) * terrain);
            b = static_cast<uint8_t>(alpha * b + (1.0f - alpha) * terrain);
        }

        file.write(reinterpret_cast<const char*>(&r), 1);
        file.write(reinterpret_cast<const char*>(&g), 1);
        file.write(reinterpret_cast<const char*>(&b), 1);
    }
}

/**
 * @brief distance using 3D Pythagorean Theorem
 * @param alpha first antenna
 * @param beta  second antenna
 * @param pixelResolution 0.5m
 * @return distance between two antennas
 */
float get_distance(const XYCoordinate alpha, const XYCoordinate beta, const float pixelResolution)
{
    float x = (alpha.x - beta.x) * pixelResolution;
    float y = (alpha.y - beta.y) * pixelResolution;
    float z = (alpha.height_m - beta.height_m);
    float d = std::sqrt( ( x*x ) + ( y*y ) + (z*z) );

    return d;
}

/**
 * @brief signal strength
 * @param actual   actual clearance
 * @param fresnel  fresnel clearance
 * @return integer (0-100) representing strength of RF signal
 */
int get_fresnel_strength(const float actual, const float fresnel)
{
    if ( fresnel <= 0.0f ) return 100;

    float clear_pct = actual / fresnel;

    if (clear_pct < 0  ) return 0;
    if (clear_pct < 0.6) return (int)(clear_pct * 25);
    if (clear_pct < 1.0) return 15 + (int)((clear_pct-0.6)*50);
    if (clear_pct < 2.0) return 35 + (int)((clear_pct-1.0)*40);
    if (clear_pct < 3.0) return 75 + (int)((clear_pct-2.0)*10);
    return 85 + std::min((int)((clear_pct-3.0) * 5), 15);
}

int main()
{
    const auto FileName = "/home/ash/CLionProjects/lidarGeo/tif/18.tif";

    GDALAllRegister();
    const GDALAccess eAccess = GA_ReadOnly;
    GDALDatasetUniquePtr poDataset = GDALDatasetUniquePtr(GDALDataset::FromHandle(GDALOpen(FileName, eAccess)));

    if( !poDataset )
    {
        cerr << "Failed to open dataset " << FileName << endl;
    }

    // #============================================================================================================# //

    double adfGeoTransform[6];
    printf("Size is X:%d Y:%d Bands:%d\n", poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), poDataset->GetRasterCount());
    printf("Driver:%s/%s\n", poDataset->GetDriver()->GetDescription(), poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ));

    if ( poDataset->GetProjectionRef() != NULL )
    {
        printf("Projection is: %s\n", poDataset->GetProjectionRef() );
    }

    if ( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        // Origin: x, y coordinates on earth
        // Pixel: real-world distance of west-east (horizontal), north-south (vertical)
        printf("Origin = (%.0f, %.0f)\n", adfGeoTransform[0], adfGeoTransform[3]);
        printf("Pixel size = (%.1fm, %.1fm)\n", adfGeoTransform[1], adfGeoTransform[5]);
    }

    GDALRasterBand *poBand;
    int             nBlockXSize, nBlockYSize;
    int             bGotMin, bGotMax;
    double          adfMinMax[2];

    poBand = poDataset->GetRasterBand( 1 );
    poBand -> GetBlockSize( &nBlockXSize, &nBlockYSize );

    // #============================================================================================================# //

    printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
        nBlockXSize, nBlockYSize,
        GDALGetDataTypeName(poBand->GetRasterDataType()),
        GDALGetColorInterpretationName(
        poBand->GetColorInterpretation()) );

    adfMinMax[0] = poBand->GetMinimum( &bGotMin );
    adfMinMax[1] = poBand->GetMaximum( &bGotMax );

    if ( ! (bGotMin && bGotMax) )
        GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);

    if ( poBand->GetOverviewCount() > 0 )
        printf("Min=%.3f has %d overviews.\n", poBand->GetOverviewCount());

    if ( poBand->GetColorTable() != NULL )
        printf("Band has a color table with %d entries.\n",
                poBand->GetColorTable()->GetColorEntryCount());



    // #============================================================================================================# //

    // Loads entire elevation data (band 1) into memory
    int nXsize = poDataset->GetRasterXSize( );
    int nYsize = poDataset->GetRasterYSize( );
    float *elevationData = (float *)CPLMalloc(sizeof(float)*nXsize*nYsize);

    // GF_Read        -> Read not write
    // 0, 0           -> nXOff, NYOff
    // nXsize, nYsize -> size of window  raster
    // elevationData  -> buffer to read into
    // nXsize, nYsize -> buffer size, used for downscaling
    // GDT_Float32    -> data type, auto data conversion
    // 0,0            -> controls access to memory
    poBand->RasterIO(GF_Read, 0, 0, nXsize, nYsize, elevationData, nXsize, nYsize, GDT_Float32, 0, 0);

    size_t total_pixels = static_cast<size_t>(nXsize) * static_cast<size_t>(nYsize);
    printf("%zu pixels (%zu MB)\n", total_pixels, total_pixels * sizeof(float) / 1024 / 1024);

    XYCoordinate alpha = create_antenna(elevationData, 0, 0, CELL_TOWER_HEIGHT, nXsize);

    alpha.print();

    float wavelength = LIGHT_SPEED / CELL_FREQUENCY;

    printf("speed:%d (km/s) freq:%.0f (Hz) wave:%f (m)\n\n",
            LIGHT_SPEED, CELL_FREQUENCY, wavelength);


    // #============================================================================================================# //

    auto signal_map = get_rf_map(elevationData, alpha, nXsize, nYsize, wavelength);
    auto gradient = compute_grad(elevationData, nXsize, nYsize, adfGeoTransform[1]);
    transform(gradient, 0.6, 3.0);
    auto elevationDataNorm = normalize(gradient);

    draw_map(elevationDataNorm, signal_map, nXsize, nYsize, "norm.pgm");


    CPLFree(elevationData);

    return 0;
}

/* ops_transform.h */

#ifndef OPS_TRANSFORM_H
#define OPS_TRANSFORM_H

#include <algorithm>
#include <vector>
#include <cmath>
#include <array>

class elevation_view {
private:
    const float* data;
    const int width, height;

public:
    elevation_view(const float* data, int w, int h) : data(data), width(w), height(h) {}

    double operator()(int x, int y) const {
        return static_cast<double>(data[y * width + x]);
    }

    bool valid(int x, int y) const {
        return x >= 0 && x < width && y >= 0 && y < height;
    }

    int w() const { return width; }
    int h() const { return height; }
};

struct stencil {
    std::array<double, 3> coeffs;
    std::array<int, 3> offsets;
};

const stencil LEFT_EDGE  = { {{-1.5,  2.0, -0.5}}, {{0, 1, 2}} };
const stencil RIGHT_EDGE = { {{ 0.5, -2.0,  1.5}}, {{-2, -1, 0}} };
const stencil CENTER    = { {{ -0.5,  0.0,  0.5}}, {{-1, 0, 1}} };

inline double apply_stencil_x(const elevation_view& view, int x, int y, const stencil& stencil) {
    double result = 0.0;
    for (int i = 0; i < 3; ++i) {
        int px = x + stencil.offsets[i];
        if (view.valid(px, y)) {
            result += stencil.coeffs[i] * view(px, y);
        }
    }
    return result;
}

inline double apply_stencil_y(const elevation_view& view, int x, int y, const stencil& stencil) {
    double result = 0.0;
    for (int i = 0; i < 3; ++i) {
        int py = y + stencil.offsets[i];
        if (view.valid(x, py)) {
            result += stencil.coeffs[i] * view(x, py);
        }
    }
    return result;
}

inline const stencil& select_stencil(int pos, int size) {
    if (pos == 0) return LEFT_EDGE;
    if (pos == size - 1) return RIGHT_EDGE;
    return CENTER;
}

inline std::vector<float> compute_grad(
    const float* elevationData,
    int width,
    int height,
    double spacing = 0.5) {

    if (width < 3 || height < 3) {
        throw std::invalid_argument("Raster too small (need ≥3×3)");
    }

    elevation_view view(elevationData, width, height);
    std::vector<float> result(width * height);

    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            const auto& stencil_x = select_stencil(x, width);
            const auto& stencil_y = select_stencil(y, height);

            double gx = apply_stencil_x(view, x, y, stencil_x) / spacing;
            double gy = apply_stencil_y(view, x, y, stencil_y) / spacing;

            result[y * width + x] = static_cast<float>(std::sqrt(gx*gx + gy*gy));
        }
    }

    return result;
}

inline std::pair<std::vector<float>, std::vector<float>> compute_grad_components(
    const float* elevationData,
    int width,
    int height,
    double spacing = 0.5) {

    elevation_view view(elevationData, width, height);
    std::vector<float> grad_x(width * height), grad_y(width * height);

    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            int idx = y * width + x;

            grad_x[idx] = apply_stencil_x(view, x, y, select_stencil(x, width)) / spacing;
            grad_y[idx] = apply_stencil_y(view, x, y, select_stencil(y, height)) / spacing;
        }
    }

    return {grad_x, grad_y};
}

/**
 * @brief performs clamping, and non-linearity using tanh and log
 * @param gradient float vector of elevation data with derivative
 * @param clamp max value
 * @param scale scaling for tanh
 */
inline void transform(std::vector<float> &gradient, float clamp = 0.6, float scale = 3.0)
{
    for (float& value : gradient)
    {
        if ( value > clamp ) value = clamp;
        value = std::tanh(scale * value);
        value = std::log(value + 1.0f);
    }
}

/**
 * @brief normalizes into 0-255 using min and max values in gradient
 * @param gradient float vector of elevation data with derivative
 * @return normalized elevation data
 */
inline std::vector<uint8_t> normalize(const std::vector<float>& gradient)
{
    auto [min_it, max_it] = std::minmax_element(gradient.begin(), gradient.end());
    float min_val = *min_it;
    float max_val = *max_it;
    float range = max_val - min_val;

    std::vector<uint8_t> normalized(gradient.size());
    for (size_t i = 0; i < gradient.size(); i++)
        normalized[i] = static_cast<uint8_t> (255 * (gradient[i] - min_val) / range);

    return normalized;
}

#endif //OPS_TRANSFORM_H