SeaCurrentsExtension.cpp

#include "SeaCurrentsExtension.h"
#include "Quadtree.h"
 
#include "SeaCurrentsSystem.h"
#include "Utilities/VectorPermutation.h"
#include "ResourceManifest.h"
 
#include "Context.h"
#include "EntityStore.h"
#include "ExtensionRegistry.h"
#include "Components/Core/TransformComponent.h"
#include "Resources/ResourceStore.h"
#include "Serialization/EntityReader.h"
#include "Utilities/Parsing.h"
 
#include "Foundation/HashFunctions.h"
#include "Foundation/Geometry/GeodeticUtils.h"
#include "Foundation/Platform/IO.h"
#include "Foundation/Logging/Logger.h"
 
#ifndef EMSCRIPTEN
#include "hdf5.h"
#include "H5LTpublic.h"
#endif
 
#include <map>
#include <set>
 
namespace {
  Cogs::Logging::Log  gLogger = Cogs::Logging::getLogger("SeaCurrentsExtension");
 
  struct RawCurrentData {
    float mSpeed;
    float mDirection;
  };
 
  struct RawPositionData {
    float mLatitude;
    float mLongitude;
  };
 
  constexpr const char* cSeaCurrentsEntityDefinition = R"(
  { "name": "SeaCurrents",
    "components" : [
      "SeaCurrentsComponent",
      "TransformComponent",
      "SceneComponent"
    ]
  })";
 
}
 
namespace Cogs::Core {
 
  void ComputePriority(SeaCurrentsComponent& seaCurrentsComp, const std::pair<glm::vec2, glm::vec2>& bounds, float threshold, int lowestPrio);
 
  struct SeaCurrentsExtension : public Extension
  {
    static  constexpr const char* cExtensionKey = "SeaCurrents";
 
    SeaCurrentsModuleAPI api = {
      reinterpret_cast<LoadFromMemoryFn>(&SeaCurrentsExtension::loadFromMemory),
      reinterpret_cast<LoadFromFileFn>(&SeaCurrentsExtension::loadFromFile),
      nullptr,
      nullptr,
    };
 
    SeaCurrentsExtension() { ExtensionRegistry::add(this, COGS_CORE_VERSION_STRING); }
 
    virtual bool initializeStatic() override {
      SeaCurrentsComponent::registerType();
      return true;
    }
 
    virtual const char* getExtensionKey() const override { return cExtensionKey; }
    virtual const void* getPublicAPI() const override { return &api; }
 
    virtual bool initialize(Context* context) override {
 
        const std::string resourceArchive = "Cogs.Core.Extensions.SeaCurrents.zip";
 
 
#ifdef EMSCRIPTEN
 
        bool resourceArchiveAdded = false;
        auto manifest = getResourceManifest(context);
        for (auto& item : manifest) {
            if (item.find(resourceArchive) != std::string::npos) {
                context->resourceStore->addResourceArchive(item);
                resourceArchiveAdded = true;
                break;
            }
        }
        if (!resourceArchiveAdded) {
            context->resourceStore->addResourceArchive(resourceArchive);
        }
 
#else
 
        context->resourceStore->addSearchPath("../Extensions/SeaCurrents/Data/");
        context->resourceStore->addResourceArchive(resourceArchive);
 
#endif
 
      readEntityDefinition(cSeaCurrentsEntityDefinition, context->store);
 
      ExtensionRegistry::registerExtensionSystem<SeaCurrentsSystem>(context, SystemPriority::Geometry, 16);
      return true;
    }
 
    static bool loadFromMemory(void* ctx, EntityId entityId, const char* extension, void* data, size_t dataSize, uint64_t timestamp, int utmZone = -1) {
      Context* context = static_cast<Context*>(ctx);
      EntityPtr             entity = context->store->getEntity(entityId);
      SeaCurrentsComponent* component = entity ? entity->getComponent<SeaCurrentsComponent>() : nullptr;
      bool                  ret = false;
 
      if (component) {
#ifndef EMSCRIPTEN 
        // To stop build warning unused variable 'gLogger'
        (void)gLogger;
        hid_t       fileId = H5LTopen_file_image(data, dataSize, H5LT_FILE_IMAGE_DONT_COPY | H5LT_FILE_IMAGE_DONT_RELEASE);
        hid_t       groupId;
        uint32_t    numInstances;
        uint32_t    dataCodingFormat;
 
        groupId = H5Gopen(fileId, "/SurfaceCurrent", H5P_DEFAULT);
        dataCodingFormat = readAttribute<uint32_t>(groupId, "dataCodingFormat");
        numInstances = readAttribute<uint32_t>(groupId, "numInstances");
        H5Gclose(groupId);
 
        // TODO:  Figure out what multiple instances represent...  If there is such a thing.
        if (numInstances == 1) {
          groupId = H5Gopen(fileId, "/SurfaceCurrent/SurfaceCurrent.01", H5P_DEFAULT);
 
          uint64_t    firstTime = parseISO8601(readAttribute<std::string>(groupId, "dateTimeOfFirstRecord"));
          uint64_t    lastTime = parseISO8601(readAttribute<std::string>(groupId, "dateTimeOfLastRecord"));
          uint32_t    noOfRecords = readAttribute<uint32_t>(groupId, "numberOfTimes");
          uint32_t    closestRecord;
          float       eastBound = readAttribute<float>(groupId, "eastBoundLongitude");
          float       westBound = readAttribute<float>(groupId, "westBoundLongitude");
          float       southBound = readAttribute<float>(groupId, "southBoundLatitude");
          float       northBound = readAttribute<float>(groupId, "northBoundLatitude");
          hid_t       rawCurrentDataMemType = H5Tcreate(H5T_COMPOUND, sizeof(RawCurrentData));
          double      lowerEasting;
          double      lowerNorthing;
          double      upperEasting;
          double      upperNorthing;
 
          LLtoUTM(lowerEasting, lowerNorthing, southBound, westBound, utmZone, true);
          LLtoUTM(upperEasting, upperNorthing, northBound, eastBound, utmZone, true);
 
          glm::dvec2  centre((lowerEasting + upperEasting) * 0.5f, (lowerNorthing + upperNorthing) * 0.5f);
 
          H5Tinsert(rawCurrentDataMemType, "surfaceCurrentSpeed", HOFFSET(RawCurrentData, mSpeed), H5T_NATIVE_FLOAT);
          H5Tinsert(rawCurrentDataMemType, "surfaceCurrentDirection", HOFFSET(RawCurrentData, mDirection), H5T_NATIVE_FLOAT);
 
          if (timestamp <= firstTime) {
            closestRecord = 0;
          }
          else if (timestamp >= lastTime) {
            closestRecord = noOfRecords;
          }
          else {
            closestRecord = static_cast<uint32_t>((timestamp - firstTime) * noOfRecords / (lastTime - firstTime));
          }
 
          char  filename[500];
          hid_t datasetId;
 
          sprintf(filename, "/SurfaceCurrent/SurfaceCurrent.01/Group_%03d/values", std::min(std::max(1U, closestRecord), noOfRecords));
          datasetId = H5Dopen(fileId, filename, H5P_DEFAULT);
 
          switch (dataCodingFormat) {
          case 1: {   // Time series at fixed stations...
            assert(false);
            break;
          }
          case 2: {   // Regularly gridded arrays...
            uint64_t                    width = readAttribute<uint32_t>(groupId, "numPointsLongitudinal");
            uint64_t                    height = readAttribute<uint32_t>(groupId, "numPointsLatitudinal");
            std::vector<RawCurrentData> buffer(width * height);
 
            if (H5Dread(datasetId, rawCurrentDataMemType, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer.data()) >= 0) {
              float   y = southBound;
              float   dx = readAttribute<float>(groupId, "gridSpacingLongitudinal");
              float   dy = readAttribute<float>(groupId, "gridSpacingLatitudinal");
              auto    i = buffer.begin();
              double  easting;
              double  northing;
 
              component->positions.clear();
              component->speeds.clear();
              component->angles.clear();
              component->priorities.clear();
              component->setChanged();
 
              for (uint64_t cy = height; cy--; y += dy) {
                float x = westBound;
 
                for (uint64_t cx = width; cx--; x += dx, ++i) {
                  if (i->mSpeed >= 0.0f) {
                    LLtoUTM(easting, northing, y, x, utmZone, true);
 
                    easting -= centre.x;
                    northing -= centre.y;
 
                    component->positions.push_back(glm::vec2(static_cast<float>(easting), static_cast<float>(northing)));
                    component->speeds.push_back(i->mSpeed);
                    component->angles.push_back(i->mDirection);
                  }
                }
              }
 
              entity->getComponent<TransformComponent>()->coordinates = glm::vec3(centre, 0.0f);
              entity->getComponent<TransformComponent>()->setChanged();
              ret = true;
            }
            break;
          }
          case 3: {   // Ungeorectified gridded arrays...
            uint32_t                      noOfNodes = readAttribute<uint32_t>(groupId, "numberOfNodes");
            hid_t                         rawPositionDataMemType = H5Tcreate(H5T_COMPOUND, sizeof(RawPositionData));
            std::vector<RawPositionData>  positions(noOfNodes);
            hid_t                         positionsId = H5Dopen(fileId, "/SurfaceCurrent/SurfaceCurrent.01/Positioning/geometryValues", H5P_DEFAULT);
 
            H5Tinsert(rawPositionDataMemType, "latitude", HOFFSET(RawPositionData, mLatitude), H5T_NATIVE_FLOAT);
            H5Tinsert(rawPositionDataMemType, "longitude", HOFFSET(RawPositionData, mLongitude), H5T_NATIVE_FLOAT);
 
            if (H5Dread(positionsId, rawPositionDataMemType, H5S_ALL, H5S_ALL, H5P_DEFAULT, positions.data()) >= 0) {
              std::vector<RawCurrentData>   currents(noOfNodes);
 
              if (H5Dread(datasetId, rawCurrentDataMemType, H5S_ALL, H5S_ALL, H5P_DEFAULT, currents.data()) >= 0) {
                double  easting;
                double  northing;
                auto    pi = positions.begin();
                auto    ci = currents.begin();
 
                component->positions.clear();
                component->speeds.clear();
                component->angles.clear();
                component->priorities.clear();
                component->setChanged();
 
                for (; noOfNodes--; ++ci, ++pi) {
                  LLtoUTM(easting, northing, pi->mLatitude, pi->mLongitude, utmZone, true);
 
                  easting -= centre.x;
                  northing -= centre.y;
 
                  component->positions.push_back(glm::vec2(static_cast<float>(easting), static_cast<float>(northing)));
                  component->speeds.push_back(ci->mSpeed);
                  component->angles.push_back(ci->mDirection);
                }
                entity->getComponent<TransformComponent>()->coordinates = glm::vec3(centre, 0.0f);
                entity->getComponent<TransformComponent>()->setChanged();
                ret = true;
              }
            }
            H5Dclose(positionsId);
            H5Tclose(rawPositionDataMemType);
            break;
          }
          case 4: {   // moving platform...
            assert(false);
            break;
          }
          }
 
          // Some special cases
          // Terminating conditions. Single point input, malformed inputs, etc...
          if (component->positions.size() != 0) {
 
            //Guarantee size consistency of component data vectors
            component->priorities.resize(component->positions.size());
 
            if (component->positions.size() == 1) {
              component->priorities[0] = 1;
              component->priorityLevel = 1;
            }
            else {
              glm::vec2         minPos(std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
              glm::vec2         maxPos(std::numeric_limits<float>::lowest(), std::numeric_limits<float>::lowest());
 
              // compute bounds and ensure all variables have valid data
              std::vector<glm::vec2>::iterator minId, maxId;
              minId = std::min_element(component->positions.begin(), component->positions.end(), [](const glm::vec2& a, const glm::vec2& b) {
                return (a == glm::min(a, b)) ? true : false;
                });
              maxId = std::max_element(component->positions.begin(), component->positions.end(), [](const glm::vec2& a, const glm::vec2& b) {
                return (b == glm::max(a, b)) ? true : false;
                });
 
              if (minId != component->positions.end()) minPos = *minId;
              if (maxId != component->positions.end()) maxPos = *maxId;
 
              //non mutating data members
              const std::pair<glm::vec2, glm::vec2> bounds = { minPos, maxPos };
 
              // Only perform if the bounds are within a reasonable area.
              constexpr float       deltaBounds = 0.0001f;
              if (deltaBounds <= glm::distance(bounds.first, bounds.second)) {
 
                // Compute the priority of each data point
                //ComputePriority(*component, bounds, deltaBounds, (int)component->positions.size());
                // create the quadtree here priority will be value of depth
                SeaCurrentQuadtree tree(minPos, maxPos);
                for (const auto& pos : component->positions)
                {
                  tree.insert(tree.getRoot(), 0, pos);
                }
                // assign priorites by navigating the tree
                for (size_t i = 0; i < component->positions.size(); i++)
                {
                  int priority = tree.search(tree.getRoot(), component->positions[i], 0);
                  component->priorities[i] = static_cast<int>(std::round(priority));
                }
           
                // sort data by priority
                auto p = Permutation::get_permutation(component->priorities,
                  [](const int& a, const int& b) { return a < b; });
 
                // sort the component array order based on priorities->
                Permutation::apply_permutation(component->priorities, p);
                Permutation::apply_permutation(component->positions, p);
                Permutation::apply_permutation(component->speeds, p);
                Permutation::apply_permutation(component->angles, p);
                component->priorityLevel = component->priorities.back();
              }
            }
          }
 
          H5Dclose(datasetId);
          H5Tclose(rawCurrentDataMemType);
          H5Gclose(groupId);
        }
        H5Fclose(fileId);
 
        (void)extension;
#else
        (void)data;
        (void)dataSize;
        (void)timestamp;
        (void)utmZone;
        LOG_ERROR(gLogger, "%s: Unsupported extension .h5 in EMSCRIPTEN build", IO::extension(extension).c_str());
#endif  
      }
      return ret;
    }
 
    static bool loadFromFile(void* context, EntityId entityId, const char* filename, uint64_t timestamp, int utmZone = -1) {
      ResourceBuffer  buffer = reinterpret_cast<Context*>(context)->resourceStore->getResourceContents(filename);
 
      if (!buffer.empty()) {
        return loadFromMemory(context, entityId, filename, buffer.data(), buffer.size(), timestamp, utmZone);
      }
      return false;
    }
#ifndef EMSCRIPTEN
    template<typename T> static T readAttribute(hid_t ownerId, const StringView& name) {
      hid_t attributeId = H5Aopen(ownerId, name.data(), H5P_DEFAULT);
      T     value = {};
 
      if (attributeId != H5I_INVALID_HID) {
        hid_t typeId = getTypeId<T>();
 
        assert(typeId != H5I_INVALID_HID);
 
        H5Aread(attributeId, typeId, &value);
        H5Aclose(attributeId);
      }
      return value;
    }
 
    template<> std::string readAttribute<std::string>(hid_t ownerId, const StringView& name) {
      hid_t       attributeId = H5Aopen(ownerId, name.data(), H5P_DEFAULT);
      std::string value;
 
      if (attributeId != H5I_INVALID_HID) {
        H5A_info_t  info;
        char* ptr;
 
        H5Aget_info(attributeId, &info);
        if (H5Aread(attributeId, H5Aget_type(attributeId), &ptr) >= 0) {
          value = std::string(ptr, info.data_size);
          H5free_memory(ptr);
        }
        H5Aclose(attributeId);
      }
      return value;
    }
 
    template<typename T> static hid_t getTypeId()           { return H5I_INVALID_HID; }
    template<> hid_t                  getTypeId<int32_t>()  { return H5T_NATIVE_INT32; }
    template<> hid_t                  getTypeId<uint32_t>() { return H5T_NATIVE_UINT32; }
    template<> hid_t                  getTypeId<float>()    { return H5T_NATIVE_FLOAT; }
    template<> hid_t                  getTypeId<double>()   { return H5T_NATIVE_DOUBLE; }
#endif
  }
  seaCurrentsExtensionInstance;
 
 
  void ComputePriority(SeaCurrentsComponent& seaCurrentsComp, const std::pair<glm::vec2, glm::vec2>& bounds, float threshold, int lowestPrio) {
 
    // if all points have a priority we are finished.
    if (std::find(seaCurrentsComp.priorities.begin(), seaCurrentsComp.priorities.end(), 0) == seaCurrentsComp.priorities.end()) {
      return;
    }
 
    std::map<int, std::vector<int>> neighbours;
    int priorityModifier = 0;
 
    //for each point count the number of neighbours
    for (int idx = 0, pointCount = (int)seaCurrentsComp.positions.size(); idx < pointCount; ++idx) {
 
      if (seaCurrentsComp.priorities[idx] != 0) continue;
 
      for (int idy = 0; idy < pointCount; ++idy) {
        // if points are close enough create the relationship
 
        if (idx == idy) continue;
 
        if (glm::distance(seaCurrentsComp.positions[idx], seaCurrentsComp.positions[idy]) < threshold) {
          neighbours[idx].push_back(idy);
          if (seaCurrentsComp.priorities[idy] == 0) neighbours[idy].push_back(idx);
        }
      }
    }
 
    // if no neighbours are found we need to keep checking  with a bigger threshold, otherwise sort the data and assign priorities.
    if (!neighbours.empty())
    {
      std::vector<std::pair<int, int>> sizeKeyTable;
      std::set<int> toIgnore;
 
      // sort data by number of neighbours
      sizeKeyTable.reserve(neighbours.size());
      for (auto& k : neighbours) {
        if (seaCurrentsComp.priorities[k.first] == 0) {
          sizeKeyTable.emplace_back(static_cast<int>(k.second.size()), k.first);
        }
      }
 
      std::stable_sort(sizeKeyTable.begin(), sizeKeyTable.end(), [](const std::pair<int, int>& a, const std::pair<int, int>& b) { return a.first > b.first; });
 
      // assign priority to clustered neighbours and continue with the recursion
      for (const auto& point : sizeKeyTable) {
 
        // with the sorted dataset we have to assign a priority to the most crowded point, but we also have to skip it's immediate neighbours
        if (std::find(toIgnore.begin(), toIgnore.end(), point.second) != toIgnore.end())
          continue;
 
        int32_t val = int32_t(lowestPrio - (sizeKeyTable.begin()->first - (point.first - 1))); // decrement prio by (max - (current-1)) 
        seaCurrentsComp.priorities[point.second] = val;
 
        for (auto& k : neighbours[point.second]) { // ignore immediate neighbours
          toIgnore.insert(k);
        }
      }
 
      priorityModifier = sizeKeyTable.begin()->first;
      // process discarded points within this threshold
      ComputePriority(seaCurrentsComp, bounds, threshold, lowestPrio - priorityModifier);
    }
    else // increase threshold and continue
      ComputePriority(seaCurrentsComp, bounds, threshold * 2.0f, lowestPrio - priorityModifier);
  }
 
 
}
 
#if defined( _WIN32 )
#define SEA_CURRENTS_EXT_API __declspec(dllexport)
#else
#define SEA_CURRENTS_EXT_API
#endif
 
extern "C" SEA_CURRENTS_EXT_API Cogs::Core::Extension * getExtension() {
  return &Cogs::Core::seaCurrentsExtensionInstance;
}
 
extern "C" SEA_CURRENTS_EXT_API SeaCurrentsModuleAPI * getAPI() {
  return &Cogs::Core::seaCurrentsExtensionInstance.api;
}