LUT1d.icc

#if !defined(ISCE_CORE_LUT1D_ICC)
#error "LUT1d.icc is an implementation detail of class LUT1d"
#endif
 
#include <pyre/journal.h>

template <typename T>
T isce3::core::LUT1d<T>::
eval(double x) const {
 
    // Check if data are available; if not, return ref value
    T value = _refValue;
    if (!_haveData) {
        return value;
    }
 
    // Check bounds to see if we need to perform linear extrapolation
    const int n = _coords.size();
    if (x < _coords[0]) {
        if (_extrapolate) {
            const double dx = _coords[0] - _coords[1];
            const double dy = _values[0] - _values[1];
            const double d = x - _coords[1];
            T result = (dy / dx) * d + _values[1];
            return result;
        } else {
            pyre::journal::error_t errorChannel("isce.core.LUT1d");
            errorChannel
                << pyre::journal::at(__HERE__)
                << "Out of bounds evaluation for LUT1d."
                << pyre::journal::newline
                << pyre::journal::endl;
            return 0;
        }
    } else if (x > _coords[n-1]) {
        if (_extrapolate) {
            const double dx = _coords[n-1] - _coords[n-2];
            const double dy = _values[n-1] - _values[n-2];
            const double d = x - _coords[n-2];
            T result = (dy / dx) * d + _values[n-2];
            return result;
        } else {
            pyre::journal::error_t errorChannel("isce.core.LUT1d");
            errorChannel
                << pyre::journal::at(__HERE__)
                << "Out of bounds evaluation for LUT1d."
                << pyre::journal::newline
                << pyre::journal::endl;
            return 0;
        }
    }
 
    // Otherwise, proceed with interpolation
    // Binary search to find leftmost coordinate
    int low = 0;
    int high = _coords.size();
    while (low < high) {
        const int midpoint = (low + high) / 2;
        if (_coords[midpoint] < x) {
            low = midpoint + 1;
        } else {
            high = midpoint;
        }
    }
 
    // Check search convergence
    if (low != high) {
        pyre::journal::error_t errorChannel("isce.core.LUT1d");
        errorChannel
            << pyre::journal::at(__HERE__)
            << "Binary search did not converge."
            << pyre::journal::newline
            << pyre::journal::endl;
        return 0;
    }
 
    // Check if right on top of a coordinate
    if (std::abs(_coords[high] - x) < 1.0e-12) {
        return _values[high];
    }
 
    // The indices of the x bounds
    const int j0 = high - 1;
    const int j1 = high;
 
    // Get coordinates at bounds
    double x1 = _coords[j0];
    double x2 = _coords[j1];
 
    // Interpolate
    T result = (x2 - x) / (x2 - x1) * _values[j0] + (x - x1) / (x2 - x1) * _values[j1];
    return result;
}
 
template<typename T>
typename isce3::core::LUT1d<T>::ArrayXt
isce3::core::LUT1d<T>::eval(const Eigen::Ref<const Eigen::ArrayXd> & x) const {
    auto out = ArrayXt(x.size());
    #pragma omp parallel for
    for(Eigen::Index n=0; n < x.size(); ++n)
      out(n) = eval(x(n));
    return out;
}
 
template <typename T>
isce3::core::LUT1d<T>
isce3::core::avgLUT2dToLUT1d(const isce3::core::LUT2d<T> & lut2d,
                          const int axis) {
    if (axis != 0 && axis != 1) {
        std::string error_msg = "ERROR axis not equal to 0 or 1";
        throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
    }
 
    // Check if LUT2d has actual data; if not, just return LUT1d with reference value
    if (!lut2d.haveData()) {
        return isce3::core::LUT1d<T>(lut2d.refValue());
    }
 
    // Determine lut1d size and number of elements to sum
    const auto lut1d_size = (axis == 0) ? lut2d.width() : lut2d.length();
    const double n_to_sum_f = (axis == 0) ? static_cast<double>(lut2d.length()) :
                                            static_cast<double>(lut2d.width());
 
    // Get a reference to the LUT2d data
    const Matrix<T> & data = lut2d.data();
 
    // Compute sum and average
    isce3::core::EArray2D<double> ea_values;
    if (axis == 0)
        // Sum along rows (x-direction)
        ea_values = data.map().colwise().sum();
    else
        // Sum along columns (y-direction)
        ea_values = data.map().rowwise().sum();
    ea_values /= n_to_sum_f;
 
    // Initialize working valarrays for computing mean along y-direction
    std::valarray<double> values(0.0, lut1d_size);
    std::copy(ea_values.data(), ea_values.data() + ea_values.size(), begin(values));
 
    // Compute final coordinates and values
    std::valarray<double> coords(lut1d_size);
    for (size_t j = 0; j < lut1d_size; ++j) {
        if (axis == 0)
            coords[j] = lut2d.xStart() + j * lut2d.xSpacing();
        else
            coords[j] = lut2d.yStart() + j * lut2d.ySpacing();
    }
 
    return isce3::core::LUT1d<T>(coords, values, true);
}