diff --git a/six/modules/c++/six.sicd/include/six/sicd/ImageData.h b/six/modules/c++/six.sicd/include/six/sicd/ImageData.h index ee30f220e..13af4affd 100644 --- a/six/modules/c++/six.sicd/include/six/sicd/ImageData.h +++ b/six/modules/c++/six.sicd/include/six/sicd/ImageData.h @@ -99,11 +99,11 @@ struct ImageData bool validate(const GeoData& geoData, logging::Logger& log) const; - static void testing_fromComplex_(std::span, std::span); // for unit-tests + static void testing_fromComplex_(std::span, std::span); // for unit-tests - static void toComplex(const six::Amp8iPhs8iLookup_t& lookup, std::span, std::span); - void toComplex(std::span, std::span) const; - void fromComplex(std::span, std::span) const; + static void toComplex(const six::Amp8iPhs8iLookup_t& lookup, std::span, std::span); + void toComplex(std::span, std::span) const; + void fromComplex(std::span, std::span) const; /*! * Create a lookup table for converting from AMP8I_PHS8I to complex. diff --git a/six/modules/c++/six.sicd/source/ComplexToAMP8IPHS8I.cpp b/six/modules/c++/six.sicd/source/ComplexToAMP8IPHS8I.cpp index a6bb542ce..c67e09476 100644 --- a/six/modules/c++/six.sicd/source/ComplexToAMP8IPHS8I.cpp +++ b/six/modules/c++/six.sicd/source/ComplexToAMP8IPHS8I.cpp @@ -25,10 +25,16 @@ #include #include -#include +#include +#include +#include + +#include +#if CODA_OSS_cpp17 +#include +#endif #include -#include #include #include @@ -61,17 +67,17 @@ The resulting green point is then what's used to find the nearest magnitude via * @param v complex value * @return phase between [0, 2PI] */ -inline long double GetPhase(const std::complex& v) +inline auto GetPhase(const std::complex& v) { - auto phase = std::arg(v); - if (phase < 0.0) phase += M_PI * 2.0; // Wrap from [0, 2PI] + double phase = std::arg(v); + if (phase < 0.0) phase += std::numbers::pi * 2.0; // Wrap from [0, 2PI] return phase; } template -static std::vector make_magnitudes_(TToComplexFunc toComplex) +static std::vector make_magnitudes_(TToComplexFunc toComplex) { - std::vector retval; + std::vector retval; retval.reserve(UINT8_MAX + 1); for (const auto amplitude : six::sicd::Utilities::iota_0_256()) { @@ -106,8 +112,8 @@ static inline auto make_magnitudes() return make_magnitudes_(toComplex); } -static const std::vector& get_magnitudes(const six::AmplitudeTable* pAmplitudeTable, - std::vector& uncached_magnitudes) +static const std::vector& get_magnitudes(const six::AmplitudeTable* pAmplitudeTable, + std::vector& uncached_magnitudes) { if (pAmplitudeTable == nullptr) { @@ -130,8 +136,8 @@ six::sicd::details::ComplexToAMP8IPHS8I::ComplexToAMP8IPHS8I(const six::Amplitud size_t i = 0; for(const auto value : six::sicd::Utilities::iota_0_256()) { - const units::Radians angle{ p0 + value * phase_delta }; - long double y, x; + const units::Radians angle{ static_cast(p0) + value * phase_delta }; + float y, x; SinCos(angle, y, x); phase_directions[i] = { x, y }; i++; @@ -143,7 +149,7 @@ six::sicd::details::ComplexToAMP8IPHS8I::ComplexToAMP8IPHS8I(const six::Amplitud * @param value query value * @return index of nearest value within the iterator range. */ -static inline uint8_t nearest(const std::vector& magnitudes, long double value) +static inline uint8_t nearest(const std::vector& magnitudes, float value) { const auto begin = magnitudes.begin(); const auto end = magnitudes.end(); @@ -179,6 +185,59 @@ six::AMP8I_PHS8I_t six::sicd::details::ComplexToAMP8IPHS8I::nearest_neighbor(con return retval; } + // Yes, this is duplicated code :-( 1) hopefully it will go away someday "soon," + // that is, we'll be at C++17; 2) the cutoff/dimension values may be different. + // + // First of all, C++11's std::async() is now (in 2023) thought of as maybe a + // bit "half baked," and perhaps shouldn't be emulated. Then, C++17 added + // parallel algorithms which might be a better way of satisfying our immediate + // needs (below) ... although we're still at C++14. +template +static inline OutputIt transform_async(const InputIt first1, const InputIt last1, OutputIt d_first, TFunc f, + typename std::iterator_traits::difference_type cutoff) +{ + // https://en.cppreference.com/w/cpp/thread/async + const auto len = std::distance(first1, last1); + if (len < cutoff) + { + return std::transform(first1, last1, d_first, f); + } + + constexpr auto policy = std::launch::async; + + const auto mid1 = first1 + len / 2; + const auto d_mid = d_first + len / 2; + auto handle = std::async(policy, transform_async, mid1, last1, d_mid, f, cutoff); + transform_async(first1, mid1, d_first, f, cutoff); + return handle.get(); +} +template +static inline void transform(std::span inputs, std::span results, TFunc f) +{ +#if CODA_OSS_cpp17 + std::ignore = std::transform(std::execution::par, inputs.begin(), inputs.end(), results.begin(), f); +#else + constexpr ptrdiff_t cutoff_ = 0; // too slow w/o multi-threading + // The value of "default_cutoff" was determined by testing; there is nothing special about it, feel free to change it. + constexpr auto dimension = 128 * 8; + constexpr auto default_cutoff = dimension * dimension; + const auto cutoff = cutoff_ == 0 ? default_cutoff : cutoff_; + std::ignore = transform_async(inputs.begin(), inputs.end(), results.begin(), f, cutoff); +#endif // CODA_OSS_cpp17 +} + +void six::sicd::details::ComplexToAMP8IPHS8I::nearest_neighbors(std::span inputs, std::span results, + const six::AmplitudeTable* pAmplitudeTable) +{ + // make a structure to quickly find the nearest neighbor + auto& converter = make(pAmplitudeTable); + const auto fromComplex_ = [&converter](const auto& v) + { + return converter.nearest_neighbor(v); + }; + transform(inputs, results, fromComplex_); +} + const six::sicd::details::ComplexToAMP8IPHS8I& six::sicd::details::ComplexToAMP8IPHS8I::make(const six::AmplitudeTable* pAmplitudeTable) { if (pAmplitudeTable == nullptr) @@ -189,6 +248,7 @@ const six::sicd::details::ComplexToAMP8IPHS8I& six::sicd::details::ComplexToAMP8 } else { + // Might have already "cached" this on the AmplitudeTable instance. auto pFromComplex = pAmplitudeTable->getFromComplex(); if (pFromComplex != nullptr) { diff --git a/six/modules/c++/six.sicd/source/ImageData.cpp b/six/modules/c++/six.sicd/source/ImageData.cpp index e61fe4430..357db6db4 100644 --- a/six/modules/c++/six.sicd/source/ImageData.cpp +++ b/six/modules/c++/six.sicd/source/ImageData.cpp @@ -31,6 +31,10 @@ #include #include +#include +#if CODA_OSS_cpp17 +#include +#endif #include "six/AmplitudeTable.h" #include "six/sicd/GeoData.h" @@ -39,7 +43,7 @@ using namespace six; using namespace six::sicd; - // There was in coda-oss, but I removed it. + // This was in coda-oss, but I removed it. // // First of all, C++11's std::async() is now (in 2023) thought of as maybe a // bit "half baked," and perhaps shouldn't be emulated. Then, C++17 added @@ -67,12 +71,15 @@ static inline OutputIt transform_async(const InputIt first1, const InputIt last1 template static inline void transform(std::span inputs, std::span results, TFunc f) { +#if CODA_OSS_cpp17 + std::ignore = std::transform(std::execution::par, inputs.begin(), inputs.end(), results.begin(), f); +#else constexpr ptrdiff_t cutoff_ = 0; // too slow w/o multi-threading - if (cutoff_ < 0) - { - std::ignore = std::transform(inputs.begin(), inputs.end(), results.begin(), f); - } - else + //if (cutoff_ < 0) + //{ + // std::ignore = std::transform(inputs.begin(), inputs.end(), results.begin(), f); + //} + //else { // The value of "default_cutoff" was determined by testing; there is nothing special about it, feel free to change it. constexpr auto dimension = 128 * 8; @@ -81,6 +88,7 @@ static inline void transform(std::span inputs, std::span inputs, std::span inputs, std::span results) const +void ImageData::fromComplex(std::span inputs, std::span results) const { - // make a structure to quickly find the nearest neighbor - auto& converter = six::sicd::details::ComplexToAMP8IPHS8I::make(amplitudeTable.get()); - const auto fromComplex_ = [&converter](const auto& v) - { - return converter.nearest_neighbor(v); - }; - transform(inputs, results, fromComplex_); + six::sicd::details::ComplexToAMP8IPHS8I::nearest_neighbors(inputs, results, amplitudeTable.get()); } -void ImageData::testing_fromComplex_(std::span inputs, std::span results) +void ImageData::testing_fromComplex_(std::span inputs, std::span results) { static const ImageData imageData; assert(imageData.amplitudeTable.get() == nullptr); diff --git a/six/modules/c++/six.sicd/source/Utilities.cpp b/six/modules/c++/six.sicd/source/Utilities.cpp index 306b41c55..e118c4836 100644 --- a/six/modules/c++/six.sicd/source/Utilities.cpp +++ b/six/modules/c++/six.sicd/source/Utilities.cpp @@ -27,12 +27,13 @@ #include #include #include -#include +#include #include #include #include #include #include +#include #include #include @@ -83,16 +84,16 @@ six::Region buildRegion(const types::RowCol& offset, } } -static std::complex toComplex_(long double A, uint8_t phase) +static auto toComplex_(double A, uint8_t phase) { // The phase values should be read in (values 0 to 255) and converted to float by doing: // P = (1 / 256) * input_value - const long double P = (1.0 / 256.0) * phase; + const double P = (1.0 / 256.0) * phase; // To convert the amplitude and phase values to complex float (i.e. real and imaginary): // S = A * cos(2 * pi * P) + j * A * sin(2 * pi * P) - const auto angle = units::Radians{ 2 * M_PI * P }; - long double sin_angle, cos_angle; + const auto angle = units::Radians{ 2 * std::numbers::pi * P }; + double sin_angle, cos_angle; SinCos(angle, sin_angle, cos_angle); std::complex S(A * cos_angle, A * sin_angle); return S; @@ -100,12 +101,12 @@ static std::complex toComplex_(long double A, uint8_t phase) std::complex six::sicd::Utilities::toComplex(uint8_t amplitude, uint8_t phase) { // A = input_amplitude(i.e. 0 to 255) - const long double A = amplitude; + const double A = amplitude; return toComplex_(A, phase); } std::complex six::sicd::Utilities::toComplex(uint8_t amplitude, uint8_t phase, const six::AmplitudeTable& amplitudeTable) { - const long double A = amplitudeTable.index(amplitude); + const double A = amplitudeTable.index(amplitude); return toComplex_(A, phase); } std::complex six::sicd::Utilities::toComplex(uint8_t amplitude, uint8_t phase, const six::AmplitudeTable* pAmplitudeTable) diff --git a/six/modules/c++/six/include/six/AmplitudeTable.h b/six/modules/c++/six/include/six/AmplitudeTable.h index 91685cc0e..6724cd114 100644 --- a/six/modules/c++/six/include/six/AmplitudeTable.h +++ b/six/modules/c++/six/include/six/AmplitudeTable.h @@ -203,16 +203,17 @@ class ComplexToAMP8IPHS8I final * @return nearest amplitude and phase value */ AMP8I_PHS8I_t nearest_neighbor(const six::zfloat& v) const; + static void nearest_neighbors(std::span inputs, std::span results, const six::AmplitudeTable*); private: //! The sorted set of possible magnitudes order from small to large. - std::vector uncached_magnitudes; // Order is important! This must be ... - const std::vector& magnitudes; // ... before this. + std::vector uncached_magnitudes; // Order is important! This must be ... + const std::vector& magnitudes; // ... before this. //! The difference in phase angle between two UINT phase values. - long double phase_delta; + float phase_delta; //! Unit vector rays that represent each direction that phase can point. - std::array, UINT8_MAX + 1> phase_directions; + std::array phase_directions; }; } }