Update OpenAL-soft to 1.23.1-bc7cb17.

libs/openal-soft/utils/makemhr/loadsofa.cpp

@@ -33,10 +33,14 @@
 #include <iterator>
 #include <memory>
 #include <numeric>
+#include <optional>
 #include <string>
+#include <string_view>
 #include <thread>
 #include <vector>
 
+#include "alspan.h"
+#include "alnumeric.h"
 #include "makemhr.h"
 #include "polyphase_resampler.h"
 #include "sofa-support.h"
@@ -44,6 +48,9 @@
 #include "mysofa.h"
 
 
+namespace {
+
+using namespace std::string_view_literals;
 using uint = unsigned int;
 
 /* Attempts to produce a compatible layout.  Most data sets tend to be
@@ -52,7 +59,7 @@ using uint = unsigned int;
  * possible.  Those sets that contain purely random measurements or use
  * different major axes will fail.
  */
-static bool PrepareLayout(const uint m, const float *xyzs, HrirDataT *hData)
+bool PrepareLayout(const uint m, const float *xyzs, HrirDataT *hData)
 {
     fprintf(stdout, "Detecting compatible layout...\n");
 
@@ -63,9 +70,10 @@ static bool PrepareLayout(const uint m, const float *xyzs, HrirDataT *hData)
         return false;
     }
 
-    double distances[MAX_FD_COUNT]{};
-    uint evCounts[MAX_FD_COUNT]{};
-    auto azCounts = std::vector<uint>(MAX_FD_COUNT*MAX_EV_COUNT, 0u);
+    std::array<double,MAX_FD_COUNT> distances{};
+    std::array<uint,MAX_FD_COUNT> evCounts{};
+    auto azCounts = std::vector<std::array<uint,MAX_EV_COUNT>>(MAX_FD_COUNT);
+    for(auto &azs : azCounts) azs.fill(0u);
 
     uint fi{0u}, ir_total{0u};
     for(const auto &field : fds)
@@ -74,21 +82,21 @@ static bool PrepareLayout(const uint m, const float *xyzs, HrirDataT *hData)
         evCounts[fi] = field.mEvCount;
 
         for(uint ei{0u};ei < field.mEvStart;ei++)
-            azCounts[fi*MAX_EV_COUNT + ei] = field.mAzCounts[field.mEvCount-ei-1];
+            azCounts[fi][ei] = field.mAzCounts[field.mEvCount-ei-1];
         for(uint ei{field.mEvStart};ei < field.mEvCount;ei++)
         {
-            azCounts[fi*MAX_EV_COUNT + ei] = field.mAzCounts[ei];
+            azCounts[fi][ei] = field.mAzCounts[ei];
             ir_total += field.mAzCounts[ei];
         }
 
         ++fi;
     }
     fprintf(stdout, "Using %u of %u IRs.\n", ir_total, m);
-    return PrepareHrirData(fi, distances, evCounts, azCounts.data(), hData) != 0;
+    const auto azs = al::span{azCounts}.first<MAX_FD_COUNT>();
+    return PrepareHrirData(al::span{distances}.first(fi), evCounts, azs, hData);
 }
 
-
-bool PrepareSampleRate(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
+float GetSampleRate(MYSOFA_HRTF *sofaHrtf)
 {
     const char *srate_dim{nullptr};
     const char *srate_units{nullptr};
@@ -96,21 +104,21 @@ bool PrepareSampleRate(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
     MYSOFA_ATTRIBUTE *srate_attrs{srate_array->attributes};
     while(srate_attrs)
     {
-        if(std::string{"DIMENSION_LIST"} == srate_attrs->name)
+        if("DIMENSION_LIST"sv == srate_attrs->name)
         {
             if(srate_dim)
             {
                 fprintf(stderr, "Duplicate SampleRate.DIMENSION_LIST\n");
-                return false;
+                return 0.0f;
             }
             srate_dim = srate_attrs->value;
         }
-        else if(std::string{"Units"} == srate_attrs->name)
+        else if("Units"sv == srate_attrs->name)
         {
             if(srate_units)
             {
                 fprintf(stderr, "Duplicate SampleRate.Units\n");
-                return false;
+                return 0.0f;
             }
             srate_units = srate_attrs->value;
         }
@@ -122,78 +130,72 @@ bool PrepareSampleRate(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
     if(!srate_dim)
     {
         fprintf(stderr, "Missing sample rate dimensions\n");
-        return false;
+        return 0.0f;
     }
-    if(srate_dim != std::string{"I"})
+    if(srate_dim != "I"sv)
     {
         fprintf(stderr, "Unsupported sample rate dimensions: %s\n", srate_dim);
-        return false;
+        return 0.0f;
     }
     if(!srate_units)
     {
         fprintf(stderr, "Missing sample rate unit type\n");
-        return false;
+        return 0.0f;
     }
-    if(srate_units != std::string{"hertz"})
+    if(srate_units != "hertz"sv)
     {
         fprintf(stderr, "Unsupported sample rate unit type: %s\n", srate_units);
-        return false;
+        return 0.0f;
     }
     /* I dimensions guarantees 1 element, so just extract it. */
-    hData->mIrRate = static_cast<uint>(srate_array->values[0] + 0.5f);
-    if(hData->mIrRate < MIN_RATE || hData->mIrRate > MAX_RATE)
+    if(srate_array->values[0] < float{MIN_RATE} || srate_array->values[0] > float{MAX_RATE})
     {
-        fprintf(stderr, "Sample rate out of range: %u (expected %u to %u)", hData->mIrRate,
+        fprintf(stderr, "Sample rate out of range: %f (expected %u to %u)", srate_array->values[0],
             MIN_RATE, MAX_RATE);
-        return false;
+        return 0.0f;
     }
-    return true;
+    return srate_array->values[0];
 }
 
-bool PrepareDelay(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
+enum class DelayType : uint8_t {
+    None,
+    I_R, /* [1][Channels] */
+    M_R, /* [HRIRs][Channels] */
+    Invalid,
+};
+DelayType PrepareDelay(MYSOFA_HRTF *sofaHrtf)
 {
     const char *delay_dim{nullptr};
     MYSOFA_ARRAY *delay_array{&sofaHrtf->DataDelay};
     MYSOFA_ATTRIBUTE *delay_attrs{delay_array->attributes};
     while(delay_attrs)
     {
-        if(std::string{"DIMENSION_LIST"} == delay_attrs->name)
+        if("DIMENSION_LIST"sv == delay_attrs->name)
         {
             if(delay_dim)
             {
                 fprintf(stderr, "Duplicate Delay.DIMENSION_LIST\n");
-                return false;
+                return DelayType::Invalid;
             }
             delay_dim = delay_attrs->value;
         }
         else
             fprintf(stderr, "Unexpected delay attribute: %s = %s\n", delay_attrs->name,
-                delay_attrs->value);
+                delay_attrs->value ? delay_attrs->value : "<null>");
         delay_attrs = delay_attrs->next;
     }
     if(!delay_dim)
     {
         fprintf(stderr, "Missing delay dimensions\n");
-        /*return false;*/
+        return DelayType::None;
     }
-    else if(delay_dim != std::string{"I,R"})
-    {
-        fprintf(stderr, "Unsupported delay dimensions: %s\n", delay_dim);
-        return false;
-    }
-    else if(hData->mChannelType == CT_STEREO)
-    {
-        /* I,R is 1xChannelCount. Makemhr currently removes any delay constant,
-         * so we can ignore this as long as it's equal.
-         */
-        if(delay_array->values[0] != delay_array->values[1])
-        {
-            fprintf(stderr, "Mismatched delays not supported: %f, %f\n", delay_array->values[0],
-                delay_array->values[1]);
-            return false;
-        }
-    }
-    return true;
+    if(delay_dim == "I,R"sv)
+        return DelayType::I_R;
+    if(delay_dim == "M,R"sv)
+        return DelayType::M_R;
+
+    fprintf(stderr, "Unsupported delay dimensions: %s\n", delay_dim);
+    return DelayType::Invalid;
 }
 
 bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
@@ -203,7 +205,7 @@ bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
     MYSOFA_ATTRIBUTE *ir_attrs{ir_array->attributes};
     while(ir_attrs)
     {
-        if(std::string{"DIMENSION_LIST"} == ir_attrs->name)
+        if("DIMENSION_LIST"sv == ir_attrs->name)
         {
             if(ir_dim)
             {
@@ -214,7 +216,7 @@ bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
         }
         else
             fprintf(stderr, "Unexpected IR attribute: %s = %s\n", ir_attrs->name,
-                ir_attrs->value);
+                ir_attrs->value ? ir_attrs->value : "<null>");
         ir_attrs = ir_attrs->next;
     }
     if(!ir_dim)
@@ -222,7 +224,7 @@ bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
         fprintf(stderr, "Missing IR dimensions\n");
         return false;
     }
-    if(ir_dim != std::string{"M,R,N"})
+    if(ir_dim != "M,R,N"sv)
     {
         fprintf(stderr, "Unsupported IR dimensions: %s\n", ir_dim);
         return false;
@@ -232,11 +234,11 @@ bool CheckIrData(MYSOFA_HRTF *sofaHrtf)
 
 
 /* Calculate the onset time of a HRIR. */
-static constexpr int OnsetRateMultiple{10};
-static double CalcHrirOnset(PPhaseResampler &rs, const uint rate, const uint n,
-    std::vector<double> &upsampled, const double *hrir)
+constexpr int OnsetRateMultiple{10};
+double CalcHrirOnset(PPhaseResampler &rs, const uint rate, const uint n,
+    al::span<double> upsampled, const double *hrir)
 {
-    rs.process(n, hrir, static_cast<uint>(upsampled.size()), upsampled.data());
+    rs.process({hrir, n}, upsampled);
 
     auto abs_lt = [](const double &lhs, const double &rhs) -> bool
     { return std::abs(lhs) < std::abs(rhs); };
@@ -246,8 +248,7 @@ static double CalcHrirOnset(PPhaseResampler &rs, const uint rate, const uint n,
 }
 
 /* Calculate the magnitude response of a HRIR. */
-static void CalcHrirMagnitude(const uint points, const uint n, std::vector<complex_d> &h,
-    double *hrir)
+void CalcHrirMagnitude(const uint points, const uint n, al::span<complex_d> h, double *hrir)
 {
     auto iter = std::copy_n(hrir, points, h.begin());
     std::fill(iter, h.end(), complex_d{0.0, 0.0});
@@ -256,26 +257,35 @@ static void CalcHrirMagnitude(const uint points, const uint n, std::vector<compl
     MagnitudeResponse(n, h.data(), hrir);
 }
 
-static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
+bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData, const DelayType delayType,
+    const uint outRate)
 {
     std::atomic<uint> loaded_count{0u};
 
-    auto load_proc = [sofaHrtf,hData,&loaded_count]() -> bool
+    auto load_proc = [sofaHrtf,hData,delayType,outRate,&loaded_count]() -> bool
     {
         const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
-        hData->mHrirsBase.resize(channels * hData->mIrCount * hData->mIrSize, 0.0);
+        hData->mHrirsBase.resize(channels * size_t{hData->mIrCount} * hData->mIrSize, 0.0);
         double *hrirs = hData->mHrirsBase.data();
 
+        std::vector<double> restmp;
+        std::optional<PPhaseResampler> resampler;
+        if(outRate && outRate != hData->mIrRate)
+        {
+            resampler.emplace().init(hData->mIrRate, outRate);
+            restmp.resize(sofaHrtf->N);
+        }
+
         for(uint si{0u};si < sofaHrtf->M;++si)
         {
             loaded_count.fetch_add(1u);
 
-            float aer[3]{
-                sofaHrtf->SourcePosition.values[3*si],
-                sofaHrtf->SourcePosition.values[3*si + 1],
-                sofaHrtf->SourcePosition.values[3*si + 2]
+            std::array aer{
+                sofaHrtf->SourcePosition.values[3_uz*si],
+                sofaHrtf->SourcePosition.values[3_uz*si + 1],
+                sofaHrtf->SourcePosition.values[3_uz*si + 2]
             };
-            mysofa_c2s(aer);
+            mysofa_c2s(aer.data());
 
             if(std::abs(aer[1]) >= 89.999f)
                 aer[0] = 0.0f;
@@ -284,22 +294,21 @@ static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
 
             auto field = std::find_if(hData->mFds.cbegin(), hData->mFds.cend(),
                 [&aer](const HrirFdT &fld) -> bool
-                {
-                    double delta = aer[2] - fld.mDistance;
-                    return (std::abs(delta) < 0.001);
-                });
+                { return (std::abs(aer[2] - fld.mDistance) < 0.001); });
             if(field == hData->mFds.cend())
                 continue;
 
-            double ef{(90.0+aer[1]) / 180.0 * (field->mEvCount-1)};
-            auto ei = static_cast<int>(std::round(ef));
-            ef = (ef-ei) * 180.0 / (field->mEvCount-1);
+            const double evscale{180.0 / static_cast<double>(field->mEvs.size()-1)};
+            double ef{(90.0 + aer[1]) / evscale};
+            auto ei = static_cast<uint>(std::round(ef));
+            ef = (ef - ei) * evscale;
             if(std::abs(ef) >= 0.1) continue;
 
-            double af{aer[0] / 360.0 * field->mEvs[ei].mAzCount};
-            auto ai = static_cast<int>(std::round(af));
-            af = (af-ai) * 360.0 / field->mEvs[ei].mAzCount;
-            ai %= field->mEvs[ei].mAzCount;
+            const double azscale{360.0 / static_cast<double>(field->mEvs[ei].mAzs.size())};
+            double af{aer[0] / azscale};
+            auto ai = static_cast<uint>(std::round(af));
+            af = (af-ai) * azscale;
+            ai %= static_cast<uint>(field->mEvs[ei].mAzs.size());
             if(std::abs(af) >= 0.1) continue;
 
             HrirAzT *azd = &field->mEvs[ei].mAzs[ai];
@@ -312,15 +321,40 @@ static bool LoadResponses(MYSOFA_HRTF *sofaHrtf, HrirDataT *hData)
 
             for(uint ti{0u};ti < channels;++ti)
             {
-                azd->mIrs[ti] = &hrirs[hData->mIrSize * (hData->mIrCount*ti + azd->mIndex)];
-                std::copy_n(&sofaHrtf->DataIR.values[(si*sofaHrtf->R + ti)*sofaHrtf->N],
-                    hData->mIrPoints, azd->mIrs[ti]);
+                azd->mIrs[ti] = &hrirs[(size_t{hData->mIrCount}*ti + azd->mIndex)*hData->mIrSize];
+                if(!resampler)
+                    std::copy_n(&sofaHrtf->DataIR.values[(size_t{si}*sofaHrtf->R + ti)*sofaHrtf->N],
+                        sofaHrtf->N, azd->mIrs[ti]);
+                else
+                {
+                    std::copy_n(&sofaHrtf->DataIR.values[(size_t{si}*sofaHrtf->R + ti)*sofaHrtf->N],
+                        sofaHrtf->N, restmp.begin());
+                    resampler->process(restmp, {azd->mIrs[ti], hData->mIrSize});
+                }
             }
 
-            /* TODO: Since some SOFA files contain minimum phase HRIRs,
-             * it would be beneficial to check for per-measurement delays
-             * (when available) to reconstruct the HRTDs.
-             */
+            /* Include any per-channel or per-HRIR delays. */
+            if(delayType == DelayType::I_R)
+            {
+                const float *delayValues{sofaHrtf->DataDelay.values};
+                for(uint ti{0u};ti < channels;++ti)
+                    azd->mDelays[ti] = delayValues[ti] / static_cast<float>(hData->mIrRate);
+            }
+            else if(delayType == DelayType::M_R)
+            {
+                const float *delayValues{sofaHrtf->DataDelay.values};
+                for(uint ti{0u};ti < channels;++ti)
+                    azd->mDelays[ti] = delayValues[si*sofaHrtf->R + ti] /
+                        static_cast<float>(hData->mIrRate);
+            }
+        }
+
+        if(outRate && outRate != hData->mIrRate)
+        {
+            const double scale{static_cast<double>(outRate) / hData->mIrRate};
+            hData->mIrRate = outRate;
+            hData->mIrPoints = std::min(static_cast<uint>(std::ceil(hData->mIrPoints*scale)),
+                hData->mIrSize);
         }
         return true;
     };
@@ -352,7 +386,7 @@ struct MagCalculator {
     {
         auto htemp = std::vector<complex_d>(mFftSize);
 
-        while(1)
+        while(true)
         {
             /* Load the current index to process. */
             size_t idx{mCurrent.load()};
@@ -375,8 +409,10 @@ struct MagCalculator {
     }
 };
 
+} // namespace
+
 bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSize,
-    const uint truncSize, const ChannelModeT chanMode, HrirDataT *hData)
+    const uint truncSize, const uint outRate, const ChannelModeT chanMode, HrirDataT *hData)
 {
     int err;
     MySofaHrtfPtr sofaHrtf{mysofa_load(filename, &err)};
@@ -429,39 +465,45 @@ bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSiz
     /* Assume a default head radius of 9cm. */
     hData->mRadius = 0.09;
 
-    if(!PrepareSampleRate(sofaHrtf.get(), hData) || !PrepareDelay(sofaHrtf.get(), hData)
-        || !CheckIrData(sofaHrtf.get()))
+    hData->mIrRate = static_cast<uint>(std::lround(GetSampleRate(sofaHrtf.get())));
+    if(!hData->mIrRate)
+        return false;
+
+    DelayType delayType = PrepareDelay(sofaHrtf.get());
+    if(delayType == DelayType::Invalid)
+        return false;
+
+    if(!CheckIrData(sofaHrtf.get()))
         return false;
     if(!PrepareLayout(sofaHrtf->M, sofaHrtf->SourcePosition.values, hData))
         return false;
-
-    if(!LoadResponses(sofaHrtf.get(), hData))
+    if(!LoadResponses(sofaHrtf.get(), hData, delayType, outRate))
         return false;
     sofaHrtf = nullptr;
 
-    for(uint fi{0u};fi < hData->mFdCount;fi++)
+    for(uint fi{0u};fi < hData->mFds.size();fi++)
     {
         uint ei{0u};
-        for(;ei < hData->mFds[fi].mEvCount;ei++)
+        for(;ei < hData->mFds[fi].mEvs.size();ei++)
         {
             uint ai{0u};
-            for(;ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
+            for(;ai < hData->mFds[fi].mEvs[ei].mAzs.size();ai++)
             {
                 HrirAzT &azd = hData->mFds[fi].mEvs[ei].mAzs[ai];
                 if(azd.mIrs[0] != nullptr) break;
             }
-            if(ai < hData->mFds[fi].mEvs[ei].mAzCount)
+            if(ai < hData->mFds[fi].mEvs[ei].mAzs.size())
                 break;
         }
-        if(ei >= hData->mFds[fi].mEvCount)
+        if(ei >= hData->mFds[fi].mEvs.size())
         {
             fprintf(stderr, "Missing source references [ %d, *, * ].\n", fi);
             return false;
         }
         hData->mFds[fi].mEvStart = ei;
-        for(;ei < hData->mFds[fi].mEvCount;ei++)
+        for(;ei < hData->mFds[fi].mEvs.size();ei++)
         {
-            for(uint ai{0u};ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
+            for(uint ai{0u};ai < hData->mFds[fi].mEvs[ei].mAzs.size();ai++)
             {
                 HrirAzT &azd = hData->mFds[fi].mEvs[ei].mAzs[ai];
                 if(azd.mIrs[0] == nullptr)
@@ -477,42 +519,41 @@ bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSiz
     size_t hrir_total{0};
     const uint channels{(hData->mChannelType == CT_STEREO) ? 2u : 1u};
     double *hrirs = hData->mHrirsBase.data();
-    for(uint fi{0u};fi < hData->mFdCount;fi++)
+    for(uint fi{0u};fi < hData->mFds.size();fi++)
     {
         for(uint ei{0u};ei < hData->mFds[fi].mEvStart;ei++)
         {
-            for(uint ai{0u};ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
+            for(uint ai{0u};ai < hData->mFds[fi].mEvs[ei].mAzs.size();ai++)
             {
                 HrirAzT &azd = hData->mFds[fi].mEvs[ei].mAzs[ai];
-                for(uint ti{0u};ti < channels;ti++)
+                for(size_t ti{0u};ti < channels;ti++)
                     azd.mIrs[ti] = &hrirs[hData->mIrSize * (hData->mIrCount*ti + azd.mIndex)];
             }
         }
 
-        for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;ei++)
-            hrir_total += hData->mFds[fi].mEvs[ei].mAzCount * channels;
+        for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvs.size();ei++)
+            hrir_total += hData->mFds[fi].mEvs[ei].mAzs.size() * channels;
     }
 
     std::atomic<size_t> hrir_done{0};
     auto onset_proc = [hData,channels,&hrir_done]() -> bool
     {
         /* Temporary buffer used to calculate the IR's onset. */
-        auto upsampled = std::vector<double>(OnsetRateMultiple * hData->mIrPoints);
+        auto upsampled = std::vector<double>(size_t{OnsetRateMultiple} * hData->mIrPoints);
         /* This resampler is used to help detect the response onset. */
         PPhaseResampler rs;
         rs.init(hData->mIrRate, OnsetRateMultiple*hData->mIrRate);
 
-        for(uint fi{0u};fi < hData->mFdCount;fi++)
+        for(auto &field : hData->mFds)
         {
-            for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;ei++)
+            for(auto &elev : field.mEvs.subspan(field.mEvStart))
             {
-                for(uint ai{0};ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
+                for(auto &azd : elev.mAzs)
                 {
-                    HrirAzT &azd = hData->mFds[fi].mEvs[ei].mAzs[ai];
                     for(uint ti{0};ti < channels;ti++)
                     {
                         hrir_done.fetch_add(1u, std::memory_order_acq_rel);
-                        azd.mDelays[ti] = CalcHrirOnset(rs, hData->mIrRate, hData->mIrPoints,
+                        azd.mDelays[ti] += CalcHrirOnset(rs, hData->mIrRate, hData->mIrPoints,
                             upsampled, azd.mIrs[ti]);
                     }
                 }
@@ -533,13 +574,12 @@ bool LoadSofaFile(const char *filename, const uint numThreads, const uint fftSiz
         return false;
 
     MagCalculator calculator{hData->mFftSize, hData->mIrPoints};
-    for(uint fi{0u};fi < hData->mFdCount;fi++)
+    for(auto &field : hData->mFds)
     {
-        for(uint ei{hData->mFds[fi].mEvStart};ei < hData->mFds[fi].mEvCount;ei++)
+        for(auto &elev : field.mEvs.subspan(field.mEvStart))
         {
-            for(uint ai{0};ai < hData->mFds[fi].mEvs[ei].mAzCount;ai++)
+            for(auto &azd : elev.mAzs)
             {
-                HrirAzT &azd = hData->mFds[fi].mEvs[ei].mAzs[ai];
                 for(uint ti{0};ti < channels;ti++)
                     calculator.mIrs.push_back(azd.mIrs[ti]);
             }