TriangularPlate.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "TriangularPlate.h"
 
#include <string>
#include <vector>
#include <numeric>
 
#include <QtGlobal>
 
#include "AbstractPlate.h"
#include "Angle.h"
#include "Distance.h"
#include "Intercept.h"
#include "IException.h"
#include "IString.h"
#include "Latitude.h"
#include "Longitude.h"
#include "NaifDskApi.h"
#include "SurfacePoint.h"
 
using namespace std;
 
namespace Isis {
 
  // no need to implement this method since it is private and never called
  // within this class
  // TriangularPlate::TriangularPlate() : AbstractPlate(), m_plate(3, 3, 0.0) { }
  
  TriangularPlate::TriangularPlate(const NaifTriangle &plate, const int &plateId) : 
                                   AbstractPlate(), m_plate(plate.copy()), 
                                   m_plateId(plateId) { }
  
  TriangularPlate::~TriangularPlate() { }
  
  int TriangularPlate::id() const {
    return m_plateId;
  }
 
  QString TriangularPlate::name() const {
    return "TriangularPlate";
  }
 
  Distance TriangularPlate::maxRadius() const {
    double radius = qMax(qMax(vnorm_c(m_plate[0]), vnorm_c(m_plate[1])), 
                         vnorm_c(m_plate[2]));
    return ( Distance(radius, Distance::Kilometers) );
  }
  
  Distance TriangularPlate::minRadius() const {
    double radius = qMin(qMin(vnorm_c(m_plate[0]), vnorm_c(m_plate[1])), 
                         vnorm_c(m_plate[2]));
    return ( Distance(radius, Distance::Kilometers) );
  }
  
  double TriangularPlate::area() const { 
  
    //  Get the lengths of each side
    NaifVector edge(3);
    vsub_c(m_plate[1], m_plate[0], &edge[0]);
    double s1 = vnorm_c(&edge[0]);
  
    vsub_c(m_plate[2], m_plate[0], &edge[0]);
    double s2 = vnorm_c(&edge[0]);
  
    vsub_c(m_plate[2], m_plate[1], &edge[0]);
    double s3 = vnorm_c(&edge[0]);
  
    // Heron's formula for area
    double S = (s1 + s2 + s3) / 2.0;
    double p_area = std::sqrt(S * (S - s1) * (S - s2) * (S - s3));
  
    return (p_area);
  }
  
  NaifVector TriangularPlate::normal() const {
  
    //  Get the lengths of each side
    NaifVector edge1(3);
    vsub_c(m_plate[1], m_plate[0], &edge1[0]);
  
    NaifVector edge2(3);
    vsub_c(m_plate[2], m_plate[0], &edge2[0]);
  
    NaifVector norm(3);
    ucrss_c(&edge1[0], &edge2[0], &norm[0]);
    
    return (norm);
  }
  
  NaifVector TriangularPlate::center() const {
    double third(0.33333333333333331);
    NaifVector midPt(3);
    vlcom3_c(third, m_plate[0], third, m_plate[1], third, m_plate[2], &midPt[0]);
    return (midPt);
  }
  
  Angle TriangularPlate::separationAngle(const NaifVector &raydir) const {
    //  Get two sides
    NaifVector norm(normal());
    double sepang = vsep_c(&norm[0], &raydir[0]);
    return ( Angle(sepang, Angle::Radians) );
  }
  
  bool TriangularPlate::hasIntercept(const NaifVertex &vertex, 
                                     const NaifVector &raydir) const {
    NaifVertex point;
    return (findPlateIntercept(vertex, raydir, point));
  }
  
  bool TriangularPlate::hasPoint(const Latitude &lat, 
                                 const Longitude &lon) const {
  
    //  Extend the maximum height of the plate to a resonable distance
    double maxrad = maxRadius().kilometers() * 1.5;
  
    // Create a surface point above the highest plate vertex
    SurfacePoint point(lat, lon, Distance(maxrad, Distance::Kilometers));
    NaifVertex obs(3);
    point.ToNaifArray(&obs[0]);
  
    // Set the ray direction back toward the center of the body
    NaifVector raydir(3);
    vminus_c(&obs[0], &raydir[0]);
  
    // Determine where the point/ray intercepts the plate
    NaifVertex xpt;
    return (findPlateIntercept(obs, raydir, xpt));
  }
  
  SurfacePoint *TriangularPlate::point(const Latitude &lat, 
                                       const Longitude &lon) const {
    //  Extend the maximum height of the plate
    double maxrad = maxRadius().kilometers() * 1.5;
  
    // Create a surface point 1.5 times above the highest plate vertex
    SurfacePoint point(lat, lon, Distance(maxrad, Distance::Kilometers));
    NaifVertex obs(3);
    point.ToNaifArray(&obs[0]);
  
    // Set the ray direction back toward the center of the body
    NaifVector raydir(3);
    vminus_c(&obs[0], &raydir[0]);
  
    // Determine if the point/ray intercepts the plate
    NaifVertex xpt;
    if ( !findPlateIntercept(obs, raydir, xpt) ) return (0);
  
    // Construct the intercept point and return it
    SurfacePoint *ipoint(new SurfacePoint());
    ipoint->FromNaifArray(&xpt[0]);
    return (ipoint);
  }
  
  
  Intercept *TriangularPlate::intercept(const NaifVertex &vertex, 
                                        const NaifVector &raydir) const {
    NaifVertex point;
    if ( !findPlateIntercept(vertex, raydir, point) ) return (0);
  
    // Got a valid intercept.  Construct it and return
    SurfacePoint *xpt = new SurfacePoint();
    xpt->FromNaifArray(&point[0]);
    return (new Intercept(vertex, raydir, xpt, clone()));
  }
  
 
  NaifVertex TriangularPlate::vertex(int v) const {
    NaifVertex vec(3);
    if ( (v < 0) || (v > 2) ) {
      QString msg = "Unable to get TriangularPlate vertex for index ["
                    + toString(v) + "]. Valid index range is 0-2.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    for ( int i = 0 ; i < 3 ; i++ ) {
      vec[i] = m_plate[v][i];
    }
    return (vec);
  }
  
  
  AbstractPlate *TriangularPlate::clone() const {
    return (new TriangularPlate(m_plate, m_plateId));
  }
  
  bool TriangularPlate::findPlateIntercept(const NaifVertex &obs,
                                           const NaifVector &raydir, 
                                           NaifVertex &point) const {
  
    // Construct three edges of the solid tehtrahderal between plate and observer
    NaifVector e1(3), e2(3), e3(3);
    vsub_c(m_plate[0], &obs[0], &e1[0]);
    vsub_c(m_plate[1], &obs[0], &e2[0]);
    vsub_c(m_plate[2], &obs[0], &e3[0]);
  
    // Test to see if the ray direction and plate normal are perpendicular
    NaifVector tnorm12(3);
    vcrss_c(&e1[0], &e2[0], &tnorm12[0]);
    double tdot12 = vdot_c(&raydir[0], &tnorm12[0]);
    double en   = vdot_c(&e3[0], &tnorm12[0]);
  
    //  Check for e3 perpendicular to plate normal.  If true, e3 is a linear
    // combination of e1 and e2.
    if ( qFuzzyCompare(en+1.0, 1.0) ) return (false);
  
    // Check if raydir and e3 are in same half space
    if ( (en > 0.0) && (tdot12 < 0.0) ) return (false);
    if ( (en < 0.0) && (tdot12 > 0.0) ) return (false);
  
    // Check that raydir and e1 are on the same side of the plane spanned by e2 
    // and e3.
    NaifVector tnorm23(3);
    vcrss_c(&e2[0], &e3[0], &tnorm23[0]);
    double tdot23 = vdot_c(&raydir[0], &tnorm23[0]);
  
    // Check if raydir and e3 are in same halfspace
    if ( (en > 0.0) && (tdot23 < 0.0) ) return (false);
    if ( (en < 0.0) && (tdot23 > 0.0) ) return (false);
  
    // Finally check that raydir and e2 are in the same half space bounded by e3
    // and e2.
    NaifVector tnorm31(3);
    vcrss_c(&e3[0], &e1[0], &tnorm31[0]);
    double tdot31 = vdot_c(&raydir[0], &tnorm31[0]);
  
    // Check if raydir and e2 are in same halfspace
    if ( (en > 0.0) && (tdot31 < 0.0) ) return (false);
    if ( (en < 0.0) && (tdot31 > 0.0) ) return (false);
  
    // Ok, we know raydir intersects the plate.  Compute the intercept point if
    // the denominator is not 0.
    double denom = tdot12 + tdot23 + tdot31;
 
    // NOTE: If we have gotten this far in the method,
    // we have checked that en != 0
    //     if en < 0, then tdot12 <= 0, tdot23 <= 0 and tdot31 <= 0
    //     if en > 0, then tdot12 >= 0, tdot23 >= 0 and tdot31 >= 0
    // So, in order for the denominator to be 0, then it must be that
    //     tdot12 == tdot23 == tdot31 == 0
    if ( qFuzzyCompare(denom+1.0, 1.0) ) return (false);
  
    double scale = en / denom;
    NaifVertex xpt(3);
    vlcom_c(1.0, &obs[0], scale, &raydir[0], &xpt[0]);
    point = xpt;
    return (true);
  }
 
} // namespace Isis