EmbreeShapeModel.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include "EmbreeShapeModel.h"
 
#include <numeric>
#include <float.h>
 
#include <QtGlobal>
#include <QList>
 
#include "IException.h"
#include "IString.h"
#include "Latitude.h"
#include "Longitude.h"
#include "NaifStatus.h"
#include "Pvl.h"
#include "ShapeModel.h"
#include "Target.h"
 
 
using namespace std;
 
namespace Isis {
 
  EmbreeShapeModel::EmbreeShapeModel() 
      : ShapeModel(),
        m_targetShape(0),
        m_targetManager(0),
        m_tolerance(DBL_MAX),
        m_shapeFile("") {
    // defaults for ShapeModel parent class include:
    //     name = empty string
    //     surfacePoint = null sp
    //     hasIntersection = false
    //     hasNormal = false
    //     normal = (0,0,0)
    //     hasEllipsoidIntersection = false
    setName("Embree");
  }
 
 
  EmbreeShapeModel::EmbreeShapeModel(Target *target, Pvl &pvl,
                                     EmbreeTargetManager *targetManager) 
      : ShapeModel(target),
        m_targetShape(0),
        m_targetManager(targetManager),
        m_tolerance(DBL_MAX),
        m_shapeFile("") {
 
    // defaults for ShapeModel parent class include:
    //     name = empty string
    //     surfacePoint = null sp
    //     hasIntersection = false
    //     hasNormal = false
    //     normal = (0,0,0)
    //     hasEllipsoidIntersection = false
 
    setName("Embree");  // Really is used as type in the system at present!
 
    PvlGroup &kernels = pvl.findGroup("Kernels", Pvl::Traverse);
 
    if (kernels.hasKeyword("ElevationModel")) {
      m_shapeFile = (QString) kernels["ElevationModel"];
    }
    else { // if (kernels.hasKeyword("ShapeModel")) {
      m_shapeFile = (QString) kernels["ShapeModel"];
    }
 
    try {
      // Request the EmbreeTargetShape from the manager
      // If the shapefile is being used by something else this will get a pointer
      // to the same target shape, otherwise it creates a new one.
      m_targetShape = m_targetManager->create(m_shapeFile);
    }
    catch(IException &e) {
      QString msg = "Cannot create a EmbreeShape from " + m_shapeFile;
      throw IException(e, IException::User, msg, _FILEINFO_);
    }
  }
 
 
  EmbreeShapeModel::EmbreeShapeModel(Target *target, const QString &shapefile,
                                     EmbreeTargetManager *targetManager) 
      : ShapeModel(target),
        m_targetShape(0),
        m_targetManager(targetManager),
        m_tolerance(DBL_MAX),
        m_shapeFile(shapefile) {
 
    // defaults for ShapeModel parent class include:
    //     name = empty string
    //     surfacePoint = null sp
    //     hasIntersection = false
    //     hasNormal = false
    //     normal = (0,0,0)
    //     hasEllipsoidIntersection = false
 
    setName("Embree");  // Really is used as type in the system at present!
 
    try {
      // Request the EmbreeTargetShape from the manager
      // If the shapefile is being used by something else this will get a pointer
      // to the same target shape, otherwise it creates a new one.
      m_targetShape = m_targetManager->create(m_shapeFile);
    }
    catch(IException &e) {
      QString msg = "Cannot create a EmbreeShape from " + m_shapeFile;
      throw IException(e, IException::User, msg, _FILEINFO_);
    }
  }
 
 
  EmbreeShapeModel::~EmbreeShapeModel() {
    if ( m_targetManager && !m_shapeFile.isEmpty() ) {
      m_targetManager->free(m_shapeFile);
    }
  }
 
 
  bool EmbreeShapeModel::intersectSurface(std::vector<double> observerPos,
                                          std::vector<double> lookDirection) {
    // Remove any previous intersection
    clearSurfacePoint();
 
    // Create a ray from the observer in the look direction
    RTCMultiHitRay ray(observerPos, lookDirection);
 
    m_targetShape->intersectRay(ray);
 
    // If nothing was hit
    if (ray.lastHit < 0) {
      setHasIntersection(false);
    }
    else {
      // Get the intersection point and the surface normal
      RayHitInformation hitInfo = m_targetShape->getHitInformation(ray, 0);
 
      // Update the surface point and surface normal
      updateIntersection(hitInfo);
    }
 
    return hasIntersection();
  }
 
 
  bool EmbreeShapeModel::intersectSurface(const Latitude &lat, const Longitude &lon,
                                          const std::vector<double> &observerPos,
                                          const bool &backCheck) {
    // Remove any previous intersection
    clearSurfacePoint();
 
    // Create a ray from the origin through the surface point
    RTCMultiHitRay ray = latlonToRay(lat,lon);
 
    m_targetShape->intersectRay(ray);
 
    // If no intersections (unlikely for this case), we are done!
    if ( ray.lastHit < 0 ) {
      return ( false );
    }
    LinearAlgebra::Vector observer = LinearAlgebra::vector(observerPos[0],
                                                           observerPos[1],
                                                           observerPos[2]);
 
    // Sorts hits based on distance to the observer
    QVector<RayHitInformation> hits = sortHits(ray, observer);
 
    // If desired, check occlusion
    if ( backCheck ) {
      for (int i = 0 ; i < hits.size() ; i++) {
        LinearAlgebra::Vector obsToIntersection = hits[i].intersection - observer;
        LinearAlgebra::Vector lookVector = LinearAlgebra::normalize(obsToIntersection);
 
        // Cast a ray from the observer to the intersection and check if it is occluded
        RTCOcclusionRay obsRay;
        obsRay.ray.org_x = observerPos[0];
        obsRay.ray.org_y = observerPos[1];
        obsRay.ray.org_z = observerPos[2];
        obsRay.ray.dir_x = lookVector[0];
        obsRay.ray.dir_y = lookVector[1];
        obsRay.ray.dir_z = lookVector[2];
        obsRay.ray.tnear = 0.0;
        obsRay.ray.tfar = LinearAlgebra::magnitude(obsToIntersection) - 0.0005;
        obsRay.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
        obsRay.hit.geomID = RTC_INVALID_GEOMETRY_ID;
        obsRay.hit.primID = RTC_INVALID_GEOMETRY_ID;
        obsRay.ray.mask = 0xFFFFFFFF;
        obsRay.ignorePrimID = hits[i].primID;
 
        // If the intersection point is not occluded,
        // then it is the closest intersection to the oberserver.
        if ( !m_targetShape->isOccluded(obsRay) ) {
          updateIntersection( hits[i] );
          break;
        }
      }
     }
    else {
      // If not testing for occlusion, take the hit closest to the observer
      updateIntersection( hits[0] );
    }
 
    return ( hasIntersection() );
  }
 
 
  bool EmbreeShapeModel::intersectSurface(const SurfacePoint &surfpt, 
                                          const std::vector<double> &observerPos,
                                          const bool &backCheck) {
    // Remove any previous intersection
    clearSurfacePoint();
 
    // Set up for finding all rays along origin vector through lat/lon surface point
    RTCMultiHitRay ray = pointToRay(surfpt);
 
    // Extend the ray to be 1.5 times the length of the SurfacePoint's radius
    ray.ray.tfar *= 1.5;
 
    m_targetShape->intersectRay(ray);
 
    // If no intersections (unlikely for this case), we are done!
    if ( ray.lastHit < 0 ) {
      return ( false );
    }
    // Convert the observer to a LinearAlgebra vector for occlusion testing
    LinearAlgebra::Vector observer = LinearAlgebra::vector(observerPos[0],
                                                           observerPos[1],
                                                           observerPos[2]);
 
    // Convert the surface point to a LinearAlgebra vector for sorting hits
    LinearAlgebra::Vector surfPoint(3);
    surfpt.ToNaifArray( &surfPoint[0] );
 
    // Sorts hits based on distance to the surface point
    QVector< RayHitInformation > hits = sortHits(ray, surfPoint);
 
    // If desired, check occlusion
    if ( backCheck ) {
      for (int i = 0 ; i < hits.size() ; i++) {
        LinearAlgebra::Vector obsToIntersection = hits[i].intersection - observer;
        LinearAlgebra::Vector lookVector = LinearAlgebra::normalize(obsToIntersection);
 
        // Cast a ray from the observer to the intersection and check if it is occluded
        RTCOcclusionRay obsRay;
        obsRay.ray.org_x = observerPos[0];
        obsRay.ray.org_y = observerPos[1];
        obsRay.ray.org_z = observerPos[2];
        obsRay.ray.dir_x = lookVector[0];
        obsRay.ray.dir_y = lookVector[1];
        obsRay.ray.dir_z = lookVector[2];
        obsRay.ray.tnear = 0.0;
        obsRay.ray.tfar = LinearAlgebra::magnitude(obsToIntersection);
        obsRay.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
        obsRay.hit.geomID = RTC_INVALID_GEOMETRY_ID;
        obsRay.hit.primID = RTC_INVALID_GEOMETRY_ID;
        obsRay.ray.mask = 0xFFFFFFFF;
        obsRay.ignorePrimID = hits[i].primID; 
 
        // If the intersection point is no occluded,
        // then it is the closest intersection to the oberserver.
        if ( !m_targetShape->isOccluded(obsRay) ) {
          updateIntersection( hits[i] );
          break;
        }
      }
     }
    else {
      // If not testing for occlusion, take the hit closest to the surface point
      updateIntersection( hits[0] );
    }
 
    return ( hasIntersection() );
  }
 
 
  void EmbreeShapeModel::updateIntersection(const RayHitInformation hitInfo) {
    // Flag that there is an intersection
    setHasIntersection(true);
 
    // Create the surfacepoint
    SurfacePoint intersectPoint;
    std::vector<double> intersectArray(3);
    std::copy( hitInfo.intersection.data().begin(),
               hitInfo.intersection.data().end(),
               intersectArray.begin() );
    intersectPoint.FromNaifArray( &intersectArray[0] );
    setSurfacePoint(intersectPoint);
 
    // Save the surface normal
    setNormal( hitInfo.surfaceNormal[0],
               hitInfo.surfaceNormal[1],
               hitInfo.surfaceNormal[2] );
  }
 
 
  void EmbreeShapeModel::clearSurfacePoint() {
    ShapeModel::clearSurfacePoint();
    setHasNormal(false);
    setHasLocalNormal(false);
    return;
  }
 
 
  Distance EmbreeShapeModel::localRadius(const Latitude &lat,
                                         const Longitude &lon) {
 
    // Create a ray from the origin to the surface point
    RTCMultiHitRay ray = latlonToRay(lat,lon);
 
    // Extend the ray to 2.5 times the maximum radius
    ray.ray.tfar *= 2.5;
 
    m_targetShape->intersectRay(ray);
 
    // If no intersections (unlikely for this case), we are done!
    if ( ray.lastHit < 0 ) {
      return ( Distance() );
    }
    // Otherwise, get the first intersection
    RayHitInformation hitInfo = m_targetShape->getHitInformation( ray, 0 );
 
    // Return the distance to the intersection
    return ( Distance( LinearAlgebra::magnitude(hitInfo.intersection),
                       Distance::Kilometers                            ) );
  }
 
 
  bool EmbreeShapeModel::isDEM() const {
    return false;
  }
 
 
  bool EmbreeShapeModel::isVisibleFrom(const std::vector<double> observerPos,
                                       const std::vector<double> lookDirection)  {
    //TODO check if there is a saved intersection
    // Create a ray from the observer in the look direction
    RTCMultiHitRay ray(observerPos, lookDirection);
 
    m_targetShape->intersectRay(ray);
 
    // If nothing was hit, something went really wrong. Just return false.
    if (ray.lastHit < 0) {
      return false;
    }
    else {
      // Get the new intersection point
      RayHitInformation hitInfo = m_targetShape->getHitInformation(ray, 0);
 
      // Check the distance between the new intersection and the saved intersection
      std::vector<double> intersectVect(3);
      surfaceIntersection()->ToNaifArray( &intersectVect[0] );
      LinearAlgebra::Vector oldIntersection = LinearAlgebra::vector( intersectVect[0],
                                                                     intersectVect[1],
                                                                     intersectVect[2] );
      return ( LinearAlgebra::magnitude(oldIntersection - hitInfo.intersection) < getTolerance() );
    }
  }
 
 
  void EmbreeShapeModel::calculateLocalNormal(QVector<double *> neighborPoints) {
    // Sanity check
    if ( !hasIntersection() ) { // hasIntersection()  <==>  hasNormall()
      QString mess = "Intercept point does not exist - cannot provide normal vector";
      throw IException(IException::Programmer, mess, _FILEINFO_);
    }
 
    return;
  }
 
 
  void EmbreeShapeModel::calculateDefaultNormal() {
    // ShapeModel (parent class) throws error if no intersection
     calculateSurfaceNormal();
  }
 
 
  void EmbreeShapeModel::calculateSurfaceNormal() {
    // ShapeModel (parent class) throws error if no intersection
    setNormal(ellipsoidNormal().toStdVector());// this takes care of setHasNormal(true);
    return;
  }
 
 
  QVector<double> EmbreeShapeModel::ellipsoidNormal()  {
 
    // Sanity check on state
    if ( !hasIntersection() ) {
       QString msg = "An intersection must be defined before computing the surface normal.";
       throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if ( !surfaceIntersection()->Valid() ) {
       QString msg = "The surface point intersection must be valid to compute the surface normal.";
       throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (!hasValidTarget()) {
       QString msg = "A valid target must be defined before computing the surface normal.";
       throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // Get the coordinates of the current surface point
    SpiceDouble pB[3];
    surfaceIntersection()->ToNaifArray(pB);
 
    // Get the body radii and compute the true normal of the ellipsoid
    QVector<double> norm(3);
    // need a case for target == NULL
    QVector<Distance> radii = QVector<Distance>::fromStdVector(targetRadii());
    NaifStatus::CheckErrors();
    surfnm_c(radii[0].kilometers(), radii[1].kilometers(), radii[2].kilometers(),
             pB, &norm[0]);
    NaifStatus::CheckErrors();
 
    return (norm);
  }
 
 
  double EmbreeShapeModel::incidenceAngle(const std::vector<double> &illuminatorBodyFixedPosition) {
 
    // If there is already a normal save it, because it's probably the local normal
    std::vector<double> localNormal;
    bool hadNormal = hasNormal();
    if ( hadNormal ) {
      localNormal = normal();
    }
 
    // Calculate the ellipsoid surface normal
    calculateDefaultNormal();
    
    // Use ShapeModel to calculate the ellipsoid incidence angle
    double ellipsoidEmission = ShapeModel::incidenceAngle(illuminatorBodyFixedPosition);
 
    // If there's a saved normal, reset it
    if ( hadNormal ) {
      setNormal(localNormal);
    }
 
    // Return the ellipsoid incidence angle
    return ellipsoidEmission;
  }
 
 
  RTCMultiHitRay EmbreeShapeModel::latlonToRay(const Latitude &lat, const Longitude &lon) const {
    // Initialize ray
    RTCMultiHitRay ray;
    ray.ray.org_x = 0.0;
    ray.ray.org_y = 0.0;
    ray.ray.org_z = 0.0;
 
    // Convert the lat, lon to a unit look direction
    double latAngle = lat.radians();
    double lonAngle = lon.radians();
    ray.ray.dir_x = cos(latAngle) * cos(lonAngle);
    ray.ray.dir_y = cos(latAngle) * sin(lonAngle);
    ray.ray.dir_z = sin(latAngle);
 
    // Set the ray's length to extend to the scene boundary
    ray.ray.tfar = m_targetShape->maximumSceneDistance();
 
    return (ray);
  }
 
 
  RTCMultiHitRay EmbreeShapeModel::pointToRay(const SurfacePoint &surfpt) const {
    // Setup everything but the direction component
    RTCMultiHitRay ray;
    ray.ray.org_x = 0.0;
    ray.ray.org_y = 0.0;
    ray.ray.org_z = 0.0;
    ray.ray.tnear = 0.0;
    ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID; // instance
    ray.hit.geomID = RTC_INVALID_GEOMETRY_ID;
    ray.hit.primID = RTC_INVALID_GEOMETRY_ID; // primitive id (triangle id)
    ray.ray.mask = 0xFFFFFFFF;
    ray.lastHit = -1;
 
    // Get the vector from the origin to the surface point
    std::vector<double> surfVect(3);
    surfpt.ToNaifArray( &surfVect[0] );
    LinearAlgebra::Vector direction = LinearAlgebra::vector( surfVect[0],
                                                             surfVect[1],
                                                             surfVect[2] );
 
    // Normalize the vector and store it in the ray
    direction = LinearAlgebra::normalize(direction);
    ray.ray.dir_x = direction[0];
    ray.ray.dir_y = direction[1];
    ray.ray.dir_z = direction[2];
 
    // Extend the ray to the surface point
    ray.ray.tfar = surfpt.GetLocalRadius().kilometers();
    return (ray);
  }
 
 
  QVector<RayHitInformation> EmbreeShapeModel::sortHits(RTCMultiHitRay &ray,
                                                        LinearAlgebra::Vector &observer) {
    int hitCount = ray.lastHit + 1;
    // Container that holds the hit (x, y, z) and distance from the origin
    // QMap sorts based upon the key, distance, so this handles sorting hits.
    QMap< double , RayHitInformation > hitsMap;
    for (int i = 0 ; i < hitCount ; i++) {
      RayHitInformation hitInfo = m_targetShape->getHitInformation( ray, i );
      hitsMap.insert( LinearAlgebra::magnitude(observer - hitInfo.intersection), hitInfo );
    }
 
    // Return sorted vector of hits
    return ( QVector< RayHitInformation >::fromList( hitsMap.values() ) );
  }
 
 
  double EmbreeShapeModel::getTolerance() const {
    return m_tolerance;
  }
 
 
  void EmbreeShapeModel::setTolerance(const double &tolerance) {
    m_tolerance = tolerance;
  }
}; // namespace Isis