Faster 8AMPI_PHSI -> complex (#549)

* Move typedefs to header.

* Expose the method for creating the 256x256 lookup table.

* Expose a method for converting 8-bit ampphase to complex.

* Move method into class.

* Replace sequential 8-bit to complex method with ImageData's lookup+parallel impl.

six/modules/c++/six.sicd/include/six/sicd/ImageData.h

@@ -41,6 +41,10 @@ namespace sicd
 {
 using cx_float = std::complex<float>;
 using AMP8I_PHS8I_t = std::pair<uint8_t, uint8_t>;
+//! Fixed size 256 element array of complex values.
+using input_values_t = std::array<std::complex<float>, UINT8_MAX + 1>;
+//! Fixed size 256 x 256 matrix of complex values.
+using input_amplitudes_t = std::array<input_values_t, UINT8_MAX + 1>;
 
 class GeoData;
 /*!
@@ -102,9 +106,20 @@ struct ImageData
     cx_float from_AMP8I_PHS8I(const AMP8I_PHS8I_t&) const; // for unit-tests
     static void to_AMP8I_PHS8I(const AmplitudeTable*, std::span<const cx_float>, std::span<AMP8I_PHS8I_t>, ptrdiff_t cutoff = -1); // for unit-tests
 
+    static void from_AMP8I_PHS8I(const input_amplitudes_t& lookup, std::span<const AMP8I_PHS8I_t>, std::span<cx_float>, ptrdiff_t cutoff = -1);
     void from_AMP8I_PHS8I(std::span<const AMP8I_PHS8I_t>, std::span<cx_float>, ptrdiff_t cutoff = -1) const;
     void to_AMP8I_PHS8I(std::span<const cx_float>, std::span<AMP8I_PHS8I_t>, ptrdiff_t cutoff = -1) const;
+
+    /*!
+     * Create a lookup table for converting from AMP8I_PHS8I to complex.
+     * @param pAmplitudeTable Input amplitude table. May be nullptr if no amplitude table is defined.
+     * @param pValues_ Output table's scope to keep it around past the function call. May be empty if there was no input amplitude table.
+     * @return reference to the output lookup table.
+     */
+    static const input_amplitudes_t& get_RE32F_IM32F_values(const six::AmplitudeTable* pAmplitudeTable,
+                                                            std::unique_ptr<input_amplitudes_t>& pValues_);
 };
+
 }
 }
 

six/modules/c++/six.sicd/source/ImageData.cpp

@@ -155,9 +155,6 @@ static std::vector<KDNode_t> make_nodes(const six::AmplitudeTable* pAmplitudeTab
     return retval;
 }
 
-using input_values_t = std::array<std::complex<float>, UINT8_MAX + 1>;
-using input_amplitudes_t = std::array<input_values_t, UINT8_MAX + 1>;
-
 // input_amplitudes_t is too big for the stack
 static std::unique_ptr<input_amplitudes_t> AMP8I_PHS8I_to_RE32F_IM32F_(const six::AmplitudeTable* pAmplitudeTable)
 {
@@ -175,7 +172,7 @@ static std::unique_ptr<input_amplitudes_t> AMP8I_PHS8I_to_RE32F_IM32F_(const six
 }
 
 // This is a non-templatized function so that there is copy of the "static" data with a NULL AmplutdeTable.
-static const input_amplitudes_t* get_RE32F_IM32F_values(const six::AmplitudeTable* pAmplitudeTable)
+static const input_amplitudes_t* get_cached_RE32F_IM32F_values(const six::AmplitudeTable* pAmplitudeTable)
 {
     if (pAmplitudeTable == nullptr)
     {
@@ -193,7 +190,7 @@ std::complex<float> ImageData::from_AMP8I_PHS8I(const AMP8I_PHS8I_t& input) cons
     }
 
     auto const pAmplitudeTable = amplitudeTable.get();
-    auto const pValues = get_RE32F_IM32F_values(pAmplitudeTable);
+    auto const pValues = get_cached_RE32F_IM32F_values(pAmplitudeTable);
 
     // Do we have a cahced result to use (no amplitude table)?
     // Or must it be recomputed (have an amplutude table)?
@@ -206,10 +203,10 @@ std::complex<float> ImageData::from_AMP8I_PHS8I(const AMP8I_PHS8I_t& input) cons
     return std::complex<float>(gsl::narrow_cast<float>(S.real()), gsl::narrow_cast<float>(S.imag()));
 }
 
-static const input_amplitudes_t& get_RE32F_IM32F_values(const six::AmplitudeTable* pAmplitudeTable,
+const input_amplitudes_t& ImageData::get_RE32F_IM32F_values(const six::AmplitudeTable* pAmplitudeTable,
     std::unique_ptr<input_amplitudes_t>& pValues_)
 {
-    const input_amplitudes_t* pValues = get_RE32F_IM32F_values(pAmplitudeTable);
+    const input_amplitudes_t* pValues = get_cached_RE32F_IM32F_values(pAmplitudeTable);
     if (pValues == nullptr)
     {
         assert(pAmplitudeTable != nullptr);
@@ -230,6 +227,12 @@ void ImageData::from_AMP8I_PHS8I(std::span<const AMP8I_PHS8I_t> inputs, std::spa
 
     std::unique_ptr<input_amplitudes_t> pValues_;
     const auto& values = get_RE32F_IM32F_values(amplitudeTable.get(), pValues_);
+    from_AMP8I_PHS8I(values, inputs, results, cutoff_);
+}
+
+void ImageData::from_AMP8I_PHS8I(const input_amplitudes_t& values, std::span<const AMP8I_PHS8I_t> inputs, std::span<std::complex<float>> results,
+    ptrdiff_t cutoff_)
+{
     const auto get_RE32F_IM32F_value_f = [&values](const six::sicd::AMP8I_PHS8I_t& v)
     {
         return values[v.first][v.second];

six/modules/c++/six.sicd/source/Utilities.cpp

@@ -169,27 +169,29 @@ class SICD_readerAndConverter final
     {
         auto bufferPtr = buffer + ((row - offset.row) * (elementsPerRow / 2));
 
-        // Take each (uint8_t, uint8_t) out of the temp buffer and put it into the real buffer as a std::complex<float>
-        for (size_t index = 0; index < elementsPerRow * rowsToRead; index += 2)
-        {
-            // "For amplitude and phase components, the amplitude component is  stored first."
-            const auto& input_amplitude = tempVector[index];
-            const auto& input_value = tempVector[index + 1];
-
-            *bufferPtr = six::sicd::Utilities::from_AMP8I_PHS8I(input_amplitude, input_value, pAmplitudeTable);
-            bufferPtr++;
-        }
+        // There's type mangling going on here.
+        // We're taking std::vector<uint8_t> and saying it's packed with std::pair<uint8_t, uint8_t>.
+        static_assert(sizeof(uint8_t) * 2 == sizeof(six::sicd::AMP8I_PHS8I_t), "expected packed layout in pair");
+        auto packed = reinterpret_cast<const six::sicd::AMP8I_PHS8I_t*>(tempVector.data());
+
+        // Reuse image data's conversion to complex.
+        static const ptrdiff_t kDefaultCutoff = 0;
+        size_t count = (elementsPerRow * rowsToRead) / 2;
+        std::span<const six::sicd::AMP8I_PHS8I_t> input(packed, count);
+        std::span<std::complex<float>> output(bufferPtr, input.size());
+        six::sicd::ImageData::from_AMP8I_PHS8I(lookup, input, output, kDefaultCutoff);
     }
     const types::RowCol<size_t>& offset;
     std::complex<float>* buffer;
-    const six::AmplitudeTable* pAmplitudeTable = nullptr;
+    std::unique_ptr<six::sicd::input_amplitudes_t> lookupScope;
+    const six::sicd::input_amplitudes_t& lookup;
 
 public:
     SICD_readerAndConverter(six::NITFReadControl& reader, size_t imageNumber,
 			    const types::RowCol<size_t>& offset, const types::RowCol<size_t>& extent,
                 size_t elementsPerRow,
 			    std::complex<float>* buffer,  const six::AmplitudeTable* pAmplitudeTable = nullptr)
-      : offset(offset), buffer(buffer), pAmplitudeTable(pAmplitudeTable)
+      : offset(offset), buffer(buffer), lookupScope(nullptr), lookup(six::sicd::ImageData::get_RE32F_IM32F_values(pAmplitudeTable, lookupScope))
     {
         SICDreader<T>(reader, imageNumber, offset, extent, elementsPerRow,
             [&](size_t elementsPerRow, size_t row, size_t rowsToRead, const std::vector<T>& tempVector)