gridInterp.h

Go to the documentation of this file.

#pragma once

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <limits>
#include <cmath>

namespace GRID_INTERP{

    enum class METHOD{LINEAR, 
        NEAREST, 
        IDW, 
        INVALID_METHOD 
    };

    enum class TYPE{GRID, 
        LAYER, 
        INVALID_TYPE 
    };

    class axis{
    public:
        axis(){};
        void setAxis(double origin, double dx, int N);
        void setAxis(std::vector<double> &x_in);
        void findIIT(double x_in, int &i1, int &i2, double &t)const;
        
        void dataFromFile(std::string filename);

        void reset();

        double operator()(int i)const;

        int last()const;

    private:
        std::vector<double> x;
        int N = 0;
        double origin = 0.0;
        double dx = 0.0;
        bool bConst = false;
        void findCell(int &i, int &ii, const double x)const;
    };

    void axis::reset(){
        x.clear();
        origin = 0.0;
        dx = 0.0;
        bConst = false;
        N = 0;
    }
    double axis::operator()(int i)const{
        return x[i];
    }
    int axis::last()const{
        return x.size()-1;
    }

    void axis::setAxis(double origin_in, double dx_in, int n){
        x.clear();
        for(int i = 0; i < n; i++){
            x.push_back(origin_in + static_cast<double>(i)*dx_in);
        }
        bConst = true;
        origin = origin_in;
        dx = dx_in;
        N = x.size() - 1;
    }

    void axis::setAxis(std::vector<double> &x_in){
        x = x_in;
        N = x.size() - 1;
    }

    void axis::findCell(int &i, int &ii, double x_in)const{
        switch (bConst)
        {
        case true:
            {
                i = static_cast<int>((x_in - origin)/dx);
                ii = i + 1;
            }
            break;
        default:
            {
                int mid = (i + ii)/2;
                if (x_in <= x[mid])
                    ii = mid;
                else
                    i = mid;
                
                if (ii - i <= 1)
                    return;
                else
                    findCell(i, ii, x_in);
            }
            break;
        }
    }

    void axis::findIIT(double x_in, int &i1, int &i2, double &t)const{
        if (x_in <= x[0]){
            i1 = 0;
            i2 = 0;
            t = 0.0;
            return;
        }
        if (x_in >= x[N]){
            i1 = N;
            i2 = N;
            t = 1.0;
            return;
        }
        i1  = 0;
        i2 = N;
        findCell(i1, i2, x_in);
        t = (x_in - x[i1])/(x[i2] - x[i1]);
    }


    void axis::dataFromFile(std::string filename){
        std::ifstream datafile(filename.c_str());
        if (!datafile.good()) {
            std::cout << "Can't open the file: " << filename << std::endl;
        }
        else{
            int Nticks;
            std::string line, tmp;
            getline(datafile, line);
            {
                std::istringstream inp(line.c_str());
                inp >> tmp;
                
                if (tmp.compare("CONST") == 0)
                    bConst = true;
                else
                    bConst = false;

                inp >> Nticks;
            }

            if (bConst){
                getline(datafile, line);
                std::istringstream inp(line.c_str());
                inp >> origin;
                inp >> dx;
                setAxis(origin, dx, Nticks);
            }
            else{
                x.clear();
                double xx;
                for (int i = 0; i < Nticks; i++){
                    getline(datafile, line);
                    std::istringstream inp(line.c_str());
                    inp >> xx;
                    x.push_back(xx);
                }
                N = x.size() - 1;
            }
            datafile.close();
        }
    }


    class GridValues{
    public:
        GridValues(){};
        double operator()(int l, int r, int c)const;
        double operator()(int r, int c)const;
        double operator()(int c)const;
        void resize(int nl, int nr, int nc);
        // sets the value of the l layer, r row and c column with the value v_in
        void set(int l, int r, int c, double v_in);
        void reset();
        int nx()const;
        int ny()const;
        int nz()const;
    private:
        std::vector<std::vector<std::vector<double>>> v;
    };

    double GridValues::operator()(int l, int r, int c)const{
        return v[l][r][c];
    }
    double GridValues::operator()(int r, int c)const{
        return v[0][r][c];
    }
    double GridValues::operator()(int c)const{
        return v[0][0][c];
    }
    void GridValues::resize(int nl, int nr, int nc){
        v.clear();
        //std::vector<double>(nc,0);
        //std::vector<std::vector<double>>(nr, std::vector<double>(nc,0));
        v.resize(nl,std::vector<std::vector<double>>(nr, std::vector<double>(nc,0)));
    }
    void GridValues::set(int l, int r, int c, double v_in){
        v[l][r][c] = v_in;
    }

    void GridValues::reset(){
        v.clear();
    }

    int GridValues::nx()const{
        return v[0][0].size();
    }
    int GridValues::ny()const{
        return v[0].size();
    }
    int GridValues::nz()const{
        return v.size();
    }

    template<int dim>
    class interp{
    public:
        // The constructor resets everything and allocates space for the axis.
        interp();
        
        double interpolate(double x, double y = 0, double z = 0)const;
        void setAxis(int idim, std::vector<double> &x_in);
        void setAxis(int idim, double origin_in, double dx_in, int n);

        void setMethod(METHOD method_in);

        void getDataFromFile(std::string filename);

        void setValues(int l, int r, int c, double v_in);
        void setElevation(int l, int r, int c, double v_in);

        void reset();

    private:
        std::vector<axis> a;
        GridValues v;
        GridValues elev;
        METHOD method = METHOD::INVALID_METHOD;
        double interp1D(double x)const;
        double interp2D(double x, double y)const;
        double interp3D(double x, double y, double z)const;
        //void cellLinearcorrect(int idim, double x, int &i1, int &i2, double &x1, double &x2)const;
        double translateX = 0.0;
        double translateY = 0.0;
        double rotate = 0.0;
        double cosTH;
        double sinTH;
        bool bPreTransform = false;
    };


    template<int dim>
    interp<dim>::interp(){
        a.clear();
        a.resize(dim);
    }

    template<int dim>
    void interp<dim>::reset(){
        for(int idim = 0; idim < dim; ++ idim)
            a[idim].reset();
        v.reset();
        elev.reset();
        //mode = MODE::INVALID_MODE;
        method = METHOD::INVALID_METHOD;
    }

    template <int dim>
    void interp<dim>::setAxis(int idim, std::vector<double> &x_in){
        if (idim < dim)
            a[idim].setAxis(x_in);
    }

    template <int dim>
    void interp<dim>::setAxis(int idim, double origin_in, double dx_in, int n){
        if (idim < dim)
            a[idim].setAxis(origin_in, dx_in, n);
    }

    template<int dim>
    void interp<dim>::setValues(int l, int r, int c, double v_in){
        v.set(l,r,c,v_in);
    }

    template <int dim>
    void interp<dim>::setElevation(int l, int r, int c, double v_in){
        elev.set(l, r, c, v_in);
    }

    template<int dim>
    void interp<dim>::setMethod(METHOD method_in){
        //mode = mode_in;
        method = method_in;
    }

    template <int dim>
    void interp<dim>::getDataFromFile(std::string filename){
        std::ifstream datafile(filename.c_str());
        if (!datafile.good()) {
            std::cout << "Can't open the file: " << filename << std::endl;
        }
        else{
            //MODE md;
            METHOD mthd;
            std::string line, tmp;
            std::vector<int> dims;
            
            getline(datafile, line);
            {// METHOD NX NY NZ
                std::istringstream inp(line.c_str());
                inp >> tmp;
                if (tmp.compare("GRIDDED") == 0){
                    // If the first line is GRIDDED ignore it and read the next line
                    getline(datafile, line);
                    inp.clear();
                    inp.str(line.c_str());
                    inp >> tmp;
                }

                if (tmp.compare("LINEAR") == 0)
                    mthd = METHOD::LINEAR;
                else if (tmp.compare("NEAREST") == 0)
                    mthd = METHOD::NEAREST;
                else
                    mthd = METHOD::INVALID_METHOD; 
                int d;
                for(int idim = 0; idim < dim; ++idim){
                    inp >> d;
                    dims.push_back(d);
                }

                {// Attempt to read transformation parameters
                    inp >> translateX;
                    inp >> translateY;
                    inp >> rotate;
                    if (std::abs(translateX) > 0.000001 || std::abs(translateY) > 0.000001 || std::abs(rotate) > 0.000001){
                        bPreTransform = true;
                        // Calculate the rotation values
                        double pi = 4*std::atan(1);
                        cosTH = std::cos(rotate * pi/180.0);
                        sinTH = std::sin(rotate * pi/180.0);
                    }
                }

                if (dim == 1)
                    v.resize(1,1,dims[0]);
                else if (dim == 2)
                    v.resize(1,dims[1],dims[0]);
                else if (dim == 3)
                    v.resize(dims[2], dims[1], dims[0]);
            }

            getline(datafile, line);
            {// FILEX, FILEY, FILEZ
                std::istringstream inp(line.c_str());
                for(int i = 0; i < dim; ++i){
                    inp >> tmp;
                    a[i].dataFromFile(tmp);
                }
            }
            if (dim == 1){
                double vv;
                for (int i = 0; i < dims[0]; ++i){
                    getline(datafile, line);
                    std::istringstream inp(line.c_str());
                    inp >> vv;
                    v.set(0,0,i,vv);
                }
            }
            else{
                int nz = 1;
                int ny, nx;
                if (dim == 2){
                    nx = dims[0];
                    ny = dims[1];
                }
                else if (dim == 3){
                    nx = dims[0];
                    ny = dims[1];
                    nz = dims[2];

                }
                double vv;
                for (int k = 0; k < nz; ++k)
                    for (int i = 0; i < ny; ++i){
                        getline(datafile, line);
                        std::istringstream inp(line.c_str());
                        for (int j = 0; j < nx; ++j){
                            inp >> vv;
                            v.set(k, i, j, vv);
                        }
                    }
            }
            setMethod(mthd);
            datafile.close();
        }
    }

    template<int dim>
    double interp<dim>::interpolate(double x, double y, double z)const{
        if (bPreTransform){
            x = x + translateX;
            y = y + translateY;
            x = x * cosTH - y * sinTH;
            x = x * sinTH + y * cosTH;
        }
        double outcome;
        switch (dim)
        {
        case 1:
            outcome = interp1D(x);
            break;
        case 2:
            outcome = interp2D(x,y);
            break;
        case 3:
            outcome = interp3D(x,y,z);
            break;
        default:
            outcome = std::numeric_limits<double>::quiet_NaN();
            break;
        }
        return outcome;
    }

    template<int dim>
    double interp<dim>::interp1D(double x)const{
        double outcome;
        int i1, i2;
        double t;
        a[0].findIIT(x,i1, i2, t);
        switch (method)
        {
        case METHOD::LINEAR:
            outcome = v(i1)*(1-t) + v(i2)*t;
            break;
        case METHOD::NEAREST:
            if (std::abs(a[0](i1) - x) <= std::abs(a[0](i2) - x))
                outcome = v(i1);
            else
                outcome = v(i2);
            break;
        default:
            outcome = std::numeric_limits<double>::quiet_NaN();
            break;
        }
        return outcome;
    }
  

    template<int dim>
    double interp<dim>::interp2D(double x, double y)const{
        double outcome;
        int i1, i2, j1, j2;
        double tx, ty;
        a[0].findIIT(x, j1, j2, tx);
        a[1].findIIT(y, i1, i2, ty);
        
        switch (method)
        {
        case METHOD::LINEAR:
            {
                double y1 = v(i1, j1)*(1-tx) + v(i1, j2)*tx;
                double y2 = v(i2, j1)*(1-tx) + v(i2, j2)*tx;
                outcome = y1*(1-ty) + y2*ty;
                
            }
            break;
        case METHOD::NEAREST:
            {
                int i, j;
                if (std::abs(a[0](j1) - x) <= std::abs(a[0](j2) - x))
                    j = j1;
                else
                    j = j2;

                if (std::abs(a[1](i1) - y) <= std::abs(a[1](i2) - y))
                    i = i1;
                else
                    i = i2;
                outcome = v(i,j);
            }
            break;
        default:
            break;
        }
        return outcome;
    }

    template<int dim>
    double interp<dim>::interp3D(double x, double y, double z)const{
        int i1, i2, j1, j2, k1, k2;
        double tx, ty, tz;
        a[0].findIIT(x, j1, j2, tx);
        a[1].findIIT(y, i1, i2, ty);
        a[2].findIIT(z, k1, k2, tz);
        
        switch (method)
        {
        case METHOD::LINEAR:
            {
                double z1y1 = v(k1, i1, j1)*(1-tx) + v(k1, i1, j2)*tx;
                double z1y2 = v(k1, i2, j1)*(1-tx) + v(k1, i2, j2)*tx;
                double z1 = z1y1*(1-ty) + z1y2*ty;

                double z2y1 = v(k2, i1, j1)*(1-tx) + v(k2, i1, j2)*tx;
                double z2y2 = v(k2, i2, j1)*(1-tx) + v(k2, i2, j2)*tx;
                double z2 = z2y1*(1-ty) + z2y2*ty;

                return z1*(1-tz) + z2*tz;
            }
            break;
        case METHOD::NEAREST:
            {
                int i, j, k;
                if (std::abs(a[0](j1) - x) <= std::abs(a[0](j2) - x))
                    j = j1;
                else
                    j = j2;

                if (std::abs(a[1](i1) - y) <= std::abs(a[1](i2) - y))
                    i = i1;
                else
                    i = i2;
                if (std::abs(a[2](k1) - y) <= std::abs(a[2](k2) - y))
                    k = k1;
                else
                    k = k2;
                return v(k, i, j);

            }
            break;
        default:
            {
                return -99999999.9;
                break;
            }
        }   
    }

    void writeCoordsValues(std::string filename, std::vector<std::vector<double>> &coords, std::vector<double> &v){
        std::ofstream out_file;
        out_file.open(filename.c_str());
        for (unsigned int i = 0; i < coords.size(); i++){
            for (unsigned j = 0; j < coords[i].size(); j++){
                out_file << coords[i][j] << " ";
            }
            out_file << v[i] << std::endl;
        }
        out_file.close();
    }
}