Commit 89c4481b authored by dcookastrog's avatar dcookastrog
Browse files

Removed SurfacePoint::SetRadii and target body radii from SurfacePoint

parent 2d1238d4
Loading
Loading
Loading
Loading
+71 −160
Original line number Diff line number Diff line
@@ -17,9 +17,9 @@ namespace Isis {
   *
   */
  SurfacePoint::SurfacePoint() {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
  }

  /**
@@ -27,25 +27,11 @@ namespace Isis {
   *
   */
  SurfacePoint::SurfacePoint(const SurfacePoint &other) {
    if(other.p_majorAxis) {
      p_majorAxis = new Distance(*other.p_majorAxis);
    if(other.p_localRadius) {
      p_localRadius = new Distance(*other.p_localRadius);
    }
    else {
      p_majorAxis = NULL;
    }

    if(other.p_minorAxis) {
      p_minorAxis = new Distance(*other.p_minorAxis);
    }
    else {
      p_minorAxis = NULL;
    }

    if(other.p_polarAxis) {
      p_polarAxis = new Distance(*other.p_polarAxis);
    }
    else {
      p_polarAxis = NULL;
      p_localRadius = NULL;
    }

    if(other.p_x) {
@@ -94,9 +80,9 @@ namespace Isis {
   */
  SurfacePoint::SurfacePoint(const Latitude &lat, const Longitude &lon,
      const Distance &radius) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetSphericalPoint(lat, lon, radius);
  }

@@ -119,9 +105,9 @@ namespace Isis {
  SurfacePoint::SurfacePoint(const Latitude &lat, const Longitude &lon,
      const Distance &radius, const Angle &latSigma, const Angle &lonSigma,
      const Distance &radiusSigma) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetSpherical(lat, lon, radius, latSigma, lonSigma, radiusSigma);
  }

@@ -133,9 +119,9 @@ namespace Isis {
   */
  SurfacePoint::SurfacePoint(const Latitude &lat, const Longitude &lon,
      const Distance &radius, const symmetric_matrix<double, upper> &covar) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetSpherical(lat, lon, radius, covar);
  }

@@ -149,9 +135,9 @@ namespace Isis {
   */
  SurfacePoint::SurfacePoint(const Displacement &x, const Displacement &y,
      const Displacement &z) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetRectangular(x, y, z);
  }

@@ -173,9 +159,9 @@ namespace Isis {
  SurfacePoint::SurfacePoint(const Displacement &x, const Displacement &y,
      const Displacement &z, const Distance &xSigma, const Distance &ySigma,
      const Distance &zSigma) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetRectangular(x, y, z, xSigma, ySigma, zSigma);
  }

@@ -191,9 +177,9 @@ namespace Isis {
   */
  SurfacePoint::SurfacePoint(const Displacement &x, const Displacement &y,
      const Displacement &z, const symmetric_matrix<double, upper> &covar) {
    p_localRadius = NULL;
    InitCovariance();
    InitPoint();
    InitRadii();
    SetRectangular(x, y, z, covar);
  }

@@ -227,16 +213,6 @@ namespace Isis {
    p_z = NULL;
  }

  /**
   * Initialize the target surface radii
   *
   */
  void SurfacePoint::InitRadii() {
    p_majorAxis = NULL;
    p_minorAxis = NULL;
    p_polarAxis = NULL;
  }

  
  /**
   * This is a private method to set a surface point in rectangular, body-fixed
@@ -450,6 +426,12 @@ namespace Isis {
    SetRectangularPoint(Displacement(rect[0], Displacement::Kilometers),
                        Displacement(rect[1], Displacement::Kilometers),
                        Displacement(rect[2], Displacement::Kilometers));
    if(p_localRadius) {
      *p_localRadius = radius;
    }
    else {
      p_localRadius = new Distance(radius);
    }
  }


@@ -551,35 +533,46 @@ namespace Isis {
   *                  in meters
   * @param radiusSigma Radius sigma of body-fixed coordinate of surface point
   *                  in meters
   * @internal
   *   @history  2018-06-28 Debbie A. Cook  Revised to use the local radius of
   *                 the SurfacePoint to convert distance to angle instead of the
   *                 major equatorial axis.  Also corrected longitude conversion.
   *                 See note in SurfacePoint.h.
   */
  void SurfacePoint::SetSphericalSigmasDistance(const Distance &latSigma,
    const Distance &lonSigma, const Distance &radiusSigma) {

    if (!p_majorAxis || !p_minorAxis || !p_polarAxis || !p_majorAxis->isValid() ||
        !p_minorAxis->isValid() || !p_polarAxis->isValid()) {
      IString msg = "In order to use sigmas in meter units, the equitorial "
        "radius must be set with a call to SetRadii or an appropriate "
        "constructor";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (!Valid()) {
      IString msg = "Cannot set spherical sigmas on an invalid surface point";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    double scaledLatSig = latSigma / *p_majorAxis;
    double scaledLonSig = lonSigma * cos((double)GetLatitude().radians())
                                   / *p_majorAxis;
    SetSphericalSigmas( Angle(scaledLatSig, Angle::Radians),
                        Angle(scaledLonSig, Angle::Radians), radiusSigma);
    // Convert Latitude sigma to radians
    Distance scalingRadius = GetLocalRadius();
    double latSigRadians = latSigma / scalingRadius;

    // Convert Longitude sigma to radians
    double convFactor = cos((double)GetLatitude().radians());
    double lonSigRadians;
    
    if (convFactor > 0.0000000000000001) {             
      scalingRadius *= convFactor;
      lonSigRadians = lonSigma / scalingRadius;
    }
    else {
      //  Brent Archinal suggested setting sigma to pi in the case of a point near the pole
      lonSigRadians = PI;
    }

    SetSphericalSigmas( Angle(latSigRadians, Angle::Radians),
                        Angle(lonSigRadians, Angle::Radians), radiusSigma);
  }


  /**
   * Set spherical covariance matrix
   *
   * @param covar Spherical variance/covariance matrix (radians**2)
   * @param covar Spherical variance/covariance matrix (radians**2 for lat and lon)
   *
   * @return void
   */
@@ -700,48 +693,6 @@ namespace Isis {
  }


  /**
   * Reset the radii of the surface body of the surface point
   *
   * @param majorAxis  The semi-major axis of the surface model
   * @param minorAxis  The semi-minor axis of the surface model
   * @param polarAxis  The polar axis of the surface model
   */
  void SurfacePoint::SetRadii(const Distance &majorRadius,
                              const Distance &minorRadius,
                              const Distance &polarRadius) {

    if (!majorRadius.isValid()  ||
        !minorRadius.isValid()  ||
        !polarRadius.isValid()) {
      IString msg = "Radii must be set to valid distances.  One or more radii "
        "have been set to an invalid distance.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if(p_majorAxis) {
      *p_majorAxis = majorRadius;
    }
    else {
      p_majorAxis = new Distance(majorRadius);
    }

    if(p_minorAxis) {
      *p_minorAxis = minorRadius;
    }
    else {
      p_minorAxis = new Distance(minorRadius);
    }

    if(p_polarAxis) {
      *p_polarAxis = polarRadius;
    }
    else {
      p_polarAxis = new Distance(polarRadius);
    }
  }


  /**
   * This method resets the local radius of a SurfacePoint
   *
@@ -905,12 +856,17 @@ namespace Isis {
      if (!Valid())
        return Distance();

      if (p_localRadius) {
        return *p_localRadius;
      }
      else {
        double x = p_x->meters();
        double y = p_y->meters();
        double z = p_z->meters();

        return Distance(sqrt(x*x + y*y + z*z), Distance::Meters);
      }
    }


  /**
@@ -924,14 +880,10 @@ namespace Isis {
        Angle latSigma = GetLatSigma();

        if (latSigma.isValid()) {
          if (!p_majorAxis || !p_majorAxis->isValid()) {
            IString msg = "In order to calculate sigmas in meter units, the "
              "equitorial radius must be set with a call to SetRadii.";
            throw IException(IException::Programmer, msg, _FILEINFO_);
          }
          Distance scalingRadius = GetLocalRadius();

          // Convert from radians to meters
          latSigmaDistance = latSigma.radians() * *p_majorAxis;
          latSigmaDistance = latSigma.radians() * scalingRadius;
        }
      }

@@ -940,7 +892,7 @@ namespace Isis {


  /**
   * Return the longiitude sigma in meters
   * Return the longitude sigma in meters
   *
   */
  Distance SurfacePoint::GetLonSigmaDistance() const {
@@ -950,18 +902,12 @@ namespace Isis {
      Angle lonSigma = GetLonSigma();

      if (lonSigma.isValid()) {
        if (!p_majorAxis || !p_majorAxis->isValid()) {
          IString msg = "In order to calculate sigmas in meter units, the "
            "equitorial radius must be set with a call to SetRadii.";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }

        Latitude lat = GetLatitude();
        double scaler = cos(lat.radians());
        Distance scalingRadius = cos(GetLatitude().radians()) * GetLocalRadius();

        // Convert from radians to meters and return
        if (scaler != 0.)
          lonSigmaDistance = lonSigma.radians() * *p_majorAxis / scaler;
        // TODO What do we do when the scaling radius is 0 (at the pole)?
        if (scalingRadius.meters() != 0.)
          lonSigmaDistance = lonSigma.radians() * scalingRadius;
      }
    }

@@ -1038,18 +984,18 @@ namespace Isis {
      }

  /**
   * Computes and returns the distance between two surface points. This does
   *   not currently support ellipsoids and so any attempt with points with
   *   planetary radii will fail. The average of the local radii will be
   *   used.
   * Computes and returns the distance between two surface points.  The average of 
   * the local radii will be used.
   */
  Distance SurfacePoint::GetDistanceToPoint(const SurfacePoint &other) const {
    if(p_majorAxis || p_minorAxis || p_polarAxis ||
       other.p_majorAxis || other.p_minorAxis || other.p_polarAxis) {
      IString msg = "SurfacePoint::GetDistanceToPoint not yet implemented for "
          "ellipsoids";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    // TODO What do we do for the check? We no longer have the ability to check.
    
    // if(p_majorAxis || p_minorAxis || p_polarAxis ||
    //    other.p_majorAxis || other.p_minorAxis || other.p_polarAxis) {
    //   IString msg = "SurfacePoint::GetDistanceToPoint not yet implemented for "
    //       "ellipsoids";
    //   throw IException(IException::Programmer, msg, _FILEINFO_);
    // }

    if(!Valid() || !other.Valid())
      return Distance();
@@ -1112,17 +1058,6 @@ namespace Isis {
      equal = equal && p_z == NULL && other.p_z == NULL;
    }

    if(equal && p_majorAxis && p_minorAxis && p_polarAxis) {
      equal = equal && *p_majorAxis == *other.p_majorAxis;
      equal = equal && *p_minorAxis == *other.p_minorAxis;
      equal = equal && *p_polarAxis == *other.p_polarAxis;
    }
    else if(equal) {
      equal = equal && p_majorAxis == NULL && other.p_majorAxis == NULL;
      equal = equal && p_minorAxis == NULL && other.p_minorAxis == NULL;
      equal = equal && p_polarAxis == NULL && other.p_polarAxis == NULL;
    }

    if(equal && p_rectCovar) {
      equal = equal && (*p_rectCovar)(0, 0) == (*other.p_rectCovar)(0, 0);
      equal = equal && (*p_rectCovar)(0, 1) == (*other.p_rectCovar)(0, 1);
@@ -1156,9 +1091,6 @@ namespace Isis {
    if(p_x && other.p_x &&
       p_y && other.p_y &&
       p_z && other.p_z &&
       !p_majorAxis && !other.p_majorAxis &&
       !p_minorAxis && !other.p_minorAxis &&
       !p_polarAxis && !other.p_polarAxis &&
       !p_rectCovar && !other.p_rectCovar &&
       !p_sphereCovar && !other.p_sphereCovar) {
      *p_x = *other.p_x;
@@ -1167,17 +1099,6 @@ namespace Isis {
    }
    else {
      FreeAllocatedMemory();
      if(other.p_majorAxis) {
        p_majorAxis = new Distance(*other.p_majorAxis);
      }

      if(other.p_minorAxis) {
        p_minorAxis = new Distance(*other.p_minorAxis);
      }

      if(other.p_polarAxis) {
        p_polarAxis = new Distance(*other.p_polarAxis);
      }

      if(other.p_x) {
        p_x = new Displacement(*other.p_x);
@@ -1219,19 +1140,9 @@ namespace Isis {
      p_z = NULL;
    }

    if(p_majorAxis) {
      delete p_majorAxis;
      p_majorAxis = NULL;
    }

    if(p_minorAxis) {
      delete p_minorAxis;
      p_minorAxis = NULL;
    }

    if(p_polarAxis) {
      delete p_polarAxis;
      p_polarAxis = NULL;
    if(p_localRadius) {
      delete p_localRadius;
      p_localRadius = NULL;
    }

    if(p_rectCovar) {
+14 −8
Original line number Diff line number Diff line
@@ -81,12 +81,24 @@ namespace Isis {
   *   @history 2011-10-06 Steven Lambright - Get*SigmaDistance will no longer
   *                           throw an exception if there is no stored sigma
   *                           and there is no stored target radii.
   *   @history 2018-06-28 Debbie A. Cook - Removed target body radii and all
   *                           related methods.  Any reference to the major axis (equatorial
   *                           radius) was replaced with the local radius. Technically I think
   *                           we should be using the target body minor axis (polar) for the
   *                           conversion of latitude degrees to/from distance and the local 
   *                           radius of the surface point for conversion of longitude degrees 
   *                           to/from distance.  For large bodies the percent error will be small, 
   *                           so we will use the local radius for convenience.  For small bodies, 
   *                           any bundle adjustment will likely use rectangular coordinates 
   *                           where degree conversions will not be necessary.  Added new
   *                           member p_localRadius to avoid recalculating when coordinates
   *                           have not changed.  Also corrected the longitude conversion equation
   *                           in SetSphericalSigmasDistance and GetLonSigmaDistance.
   */

  class SurfacePoint {
    public:
      // Constructors
//      SurfacePoint(const std::vector <double> radii);
      SurfacePoint();
      SurfacePoint(const SurfacePoint &other);
      SurfacePoint(const Latitude &lat, const Longitude &lon,
@@ -152,9 +164,6 @@ namespace Isis {
                                      const Distance &lonSigma,
                                      const Distance &radiusSigma);

      void SetRadii(const Distance &majorRadius, const Distance &minorRadius,
                    const Distance &polarRadius);

      void ResetLocalRadius(const Distance &radius);
      bool Valid() const;

@@ -197,14 +206,11 @@ namespace Isis {
    private:
      void InitCovariance();
      void InitPoint();
      void InitRadii();
      void SetRectangularPoint(const Displacement &x, const Displacement &y, const Displacement &z);
      void SetSphericalPoint(const Latitude &lat, const Longitude &lon, const Distance &radius);
      void FreeAllocatedMemory();

      Distance *p_majorAxis;
      Distance *p_minorAxis;
      Distance *p_polarAxis;
      Distance *p_localRadius;
      Displacement *p_x;
      Displacement *p_y;
      Displacement *p_z;
+3 −5
Original line number Diff line number Diff line
@@ -20,11 +20,13 @@ using namespace boost::numeric::ublas;

/**
 *
 * @author 2010-07-30 Tracie Sucharski, Ken L. Edmunson, and Debbie A. Cook
 * @author 2010-07-30 Tracie Sucharski, Ken L. Edmundson, and Debbie A. Cook
 *
 * @internal
 *   @history 2015-02-20 Jeannie Backer - Added print statement for
 *            latitude, longitude and radius weight accessor methods.
 *   @history 2018-06-28 Debbie A.Cook - Removed test of obsolete method
 *            SetRadii
 */
int main(int argc, char *argv[]) {
  Isis::Preference::Preferences(true);
@@ -37,10 +39,6 @@ int main(int argc, char *argv[]) {
          "sigmaX=10. m, sigmaY=50. m, sigmaZ=20. m" << endl << endl;
    Isis::SurfacePoint spRec;

    spRec.SetRadii(Distance(1000., Distance::Meters),
                   Distance(1000., Distance::Meters),
                   Distance(1000., Distance::Meters));

    symmetric_matrix<double,upper> covar;
    covar.resize(3);
    covar.clear();
+0 −6
Original line number Diff line number Diff line
@@ -790,9 +790,6 @@ int main(int argc, char *argv[]) {
    qDebug() << "Create ConstrainedPoint from constrained point with adjusted surface point "
                "(32, 120, 1000)...";
    SurfacePoint aprioriSurfPt;
    aprioriSurfPt.SetRadii(Distance(1000.0, Distance::Meters),
                           Distance(1000.0, Distance::Meters),
                           Distance(1000.0, Distance::Meters));
    boost::numeric::ublas::symmetric_matrix<double, boost::numeric::ublas::upper> covar;
    covar.resize(3);
    covar.clear();
@@ -811,9 +808,6 @@ int main(int argc, char *argv[]) {
                                Angle(1.64192315,Angle::Degrees),
                                Angle(1.78752107, Angle::Degrees),
                                Distance(38.454887335682053718134171237789, Distance::Meters));
    adjustedSurfPt.SetRadii(Distance(1000.0, Distance::Meters),
                            Distance(1000.0, Distance::Meters),
                            Distance(1000.0, Distance::Meters));
    bcp5.setAdjustedSurfacePoint(adjustedSurfPt);
    qDebug().noquote() << bcp5.formatBundleOutputSummaryString(errorProp);
    qDebug().noquote() << bcp5.formatBundleOutputDetailString(errorProp, radiansToMeters);
+0 −14
Original line number Diff line number Diff line
@@ -366,10 +366,6 @@ namespace Isis {

          if ( !m_header.targetRadii.empty() ) {

            aprioriSurfacePoint.SetRadii( m_header.targetRadii[0],
                                          m_header.targetRadii[1],
                                          m_header.targetRadii[2] );

            if ( aprioriSurfacePoint.GetLatSigmaDistance().meters() != Isis::Null
                 && aprioriSurfacePoint.GetLonSigmaDistance().meters() != Isis::Null
                 && aprioriSurfacePoint.GetLocalRadiusSigma().meters() != Isis::Null ) {
@@ -449,10 +445,6 @@ namespace Isis {

          if ( !m_header.targetRadii.empty() ) {

            adjustedSurfacePoint.SetRadii( m_header.targetRadii[0],
                                           m_header.targetRadii[1],
                                           m_header.targetRadii[2] );

            if ( adjustedSurfacePoint.GetLatSigmaDistance().meters() != Isis::Null
                 && adjustedSurfacePoint.GetLonSigmaDistance().meters() != Isis::Null
                 && adjustedSurfacePoint.GetLocalRadiusSigma().meters() != Isis::Null ) {
@@ -1573,12 +1565,6 @@ namespace Isis {

      SurfacePoint aprioriSurfacePoint = controlPoint->GetAprioriSurfacePoint();
      SurfacePoint adjustedSurfacePoint = controlPoint->GetAdjustedSurfacePoint();
      aprioriSurfacePoint.SetRadii(m_header.targetRadii[0],
                                   m_header.targetRadii[1],
                                   m_header.targetRadii[2]);
      adjustedSurfacePoint.SetRadii(m_header.targetRadii[0],
                                    m_header.targetRadii[1],
                                    m_header.targetRadii[2]);
      controlPoint->SetAdjustedSurfacePoint(adjustedSurfacePoint);
      controlPoint->SetAprioriSurfacePoint(aprioriSurfacePoint);
    }
Loading