Isis 3 Programmer Reference
SurfacePoint.cpp
1
6/* SPDX-License-Identifier: CC0-1.0 */
7#include "SurfacePoint.h"
8
9#include <SpiceUsr.h>
10
11#include "IException.h"
12#include "IString.h"
13#include "Latitude.h"
14#include "Longitude.h"
15#include "SpecialPixel.h"
16
17using namespace boost::numeric::ublas;
18using namespace std;
19
20namespace Isis {
21
30
36 p_localRadius = other.p_localRadius;
37
38 if(other.p_x) {
39 p_x = new Displacement(*other.p_x);
40 }
41 else {
42 p_x = NULL;
43 }
44
45 if(other.p_y) {
46 p_y = new Displacement(*other.p_y);
47 }
48 else {
49 p_y = NULL;
50 }
51
52 if(other.p_z) {
53 p_z = new Displacement(*other.p_z);
54 }
55 else {
56 p_z = NULL;
57 }
58
59 if(other.p_rectCovar) {
60 p_rectCovar = new symmetric_matrix<double, upper>(*other.p_rectCovar);
61 }
62 else {
63 p_rectCovar = NULL;
64 }
65
66 if(other.p_sphereCovar) {
67 p_sphereCovar = new symmetric_matrix<double, upper>(*other.p_sphereCovar);
68 }
69 else {
70 p_sphereCovar = NULL;
71 }
72 }
73
74
83 const Distance &radius) {
85 InitPoint();
86 SetSphericalPoint(lat, lon, radius);
87 }
88
89
90
106 const Distance &radius, const Angle &latSigma, const Angle &lonSigma,
107 const Distance &radiusSigma) {
109 InitPoint();
110 SetSpherical(lat, lon, radius, latSigma, lonSigma, radiusSigma);
111 }
112
113
119 SurfacePoint::SurfacePoint(const Latitude &lat, const Longitude &lon,
120 const Distance &radius, const symmetric_matrix<double, upper> &covar) {
122 InitPoint();
123 SetSpherical(lat, lon, radius, covar);
124 }
125
126
135 const Displacement &z) {
137 InitPoint();
138 SetRectangular(x, y, z);
139 }
140
141
157 const Displacement &z, const Distance &xSigma, const Distance &ySigma,
158 const Distance &zSigma) {
160 InitPoint();
161 SetRectangular(x, y, z, xSigma, ySigma, zSigma);
162 }
163
164
175 const Displacement &z, const symmetric_matrix<double, upper> &covar) {
177 InitPoint();
178 SetRectangular(x, y, z, covar);
179 }
180
181
187 FreeAllocatedMemory();
188 }
189
190
196 p_rectCovar = NULL;
197 p_sphereCovar = NULL;
198 }
199
200
206 p_x = NULL;
207 p_y = NULL;
208 p_z = NULL;
209 p_localRadius = Distance();
210 }
211
212
226 const Displacement &y, const Displacement &z) {
227
228 if (!x.isValid() || !y.isValid() || !z.isValid()) {
229 IString msg = "x, y, and z must be set to valid displacements. One or "
230 "more coordinates have been set to an invalid displacement.";
231 throw IException(IException::User, msg, _FILEINFO_);
232 }
233
234 if(p_x) {
235 *p_x = x;
236 }
237 else {
238 p_x = new Displacement(x);
239 }
240
241 if(p_y) {
242 *p_y = y;
243 }
244 else {
245 p_y = new Displacement(y);
246 }
247
248 if(p_z) {
249 *p_z = z;
250 }
251 else {
252 p_z = new Displacement(z);
253 }
254
255 // Added 07-30-2018 to avoid trying to set an invalid SurfacePoint. This breaks the ringspt and cnetnewradii tests. Also the pole...mean test.
256 // if (!(this->Valid())) {
257 // IString msg = "All coodinates of the point are zero. At least one coordinate"
258 // " must be nonzero to be a valid SurfacePoint.";
259 // throw IException(IException::User, msg, _FILEINFO_);
260 // }
261
262 if (!p_localRadius.isValid()) {
264 p_localRadius = GetLocalRadius();
265 }
266 }
267
268
284 const Displacement &y, const Displacement &z, const Distance &xSigma,
285 const Distance &ySigma, const Distance &zSigma) {
286
287 // Wipe out current local radius to ensure it will be recalculated in SetRectangularPoint
288 p_localRadius = Distance();
289
290 SetRectangularPoint(x, y, z);
291
292 if (xSigma.isValid() && ySigma.isValid() && zSigma.isValid())
293 SetRectangularSigmas(xSigma, ySigma, zSigma);
294 }
295
296
309 const symmetric_matrix<double,upper>& covar) {
310 // Wipe out current local radius to ensure it will be recalulated in SetRectangularPoint
311 p_localRadius = Distance();
312
313 SetRectangularPoint(x, y, z);
315 }
316
317
328 const Displacement &y,
329 const Displacement &z) {
330 // Wipe out current local radius to ensure it will be recalculated in SetRectangularPoint
331 p_localRadius = Distance();
332
333 SetRectangularPoint(x, y, z);
334 }
335
336
351 const Distance &ySigma,
352 const Distance &zSigma) {
353 if (!xSigma.isValid() || !ySigma.isValid() || !zSigma.isValid()) {
354 IString msg = "x sigma, y sigma , and z sigma must be set to valid "
355 "distances. One or more sigmas have been set to an invalid distance.";
356 throw IException(IException::User, msg, _FILEINFO_);
357 }
358
359 symmetric_matrix<double,upper> covar(3);
360 covar.clear();
361 covar(0,0) = xSigma.meters() * xSigma.meters();
362 covar(1,1) = ySigma.meters() * ySigma.meters();
363 covar(2,2) = zSigma.meters() * zSigma.meters();
364 SetRectangularMatrix(covar, Meters);
365 }
366
367
376 void SurfacePoint::SetRectangularMatrix(const symmetric_matrix<double, upper> &covar,
377 SurfacePoint::CoordUnits units) {
378
379 // Make sure the point is set first
380 if (!Valid()) {
381 IString msg = "A point must be set before a variance/covariance matrix "
382 "can be set.";
383 throw IException(IException::Programmer, msg, _FILEINFO_);
384 }
385
386 // Make sure the units are valid
387 if (units != Kilometers && units != Meters) {
388 IString msg = "Units must be Kilometers or Meters";
389 throw IException(IException::Programmer, msg, _FILEINFO_);
390 }
391
392 // covar units are km**2
393 if (p_rectCovar) {
394 *p_rectCovar = covar;
395 }
396 else {
397 p_rectCovar = new symmetric_matrix<double, upper>(covar);
398 }
399
400 if (units == Meters) {
401 // Convert input matrix to km to hold in memory
402 (*p_rectCovar)(0,0) = covar(0,0)/1.0e6;
403 (*p_rectCovar)(0,1) = covar(0,1)/1.0e6;
404 (*p_rectCovar)(0,2) = covar(0,2)/1.0e6;
405 (*p_rectCovar)(1,1) = covar(1,1)/1.0e6;
406 (*p_rectCovar)(1,2) = covar(1,2)/1.0e6;
407 (*p_rectCovar)(2,2) = covar(2,2)/1.0e6;
408 }
409
410 SpiceDouble rectMat[3][3];
411
412 // Compute the local radius of the surface point in meters. We will convert to km before saving the matrix.
413 double x2 = p_x->meters() * p_x->meters();
414 double y2 = p_y->meters() * p_y->meters();
415 double z = p_z->meters();
416 double radius = sqrt(x2 + y2 + z*z);
417
418 // *** TODO *** Replace this section with LinearAlgebra multiply calls and avoid having to create a Spice matrix
419 rectMat[0][0] = covar(0,0);
420 rectMat[0][1] = rectMat[1][0] = covar(0,1);
421 rectMat[0][2] = rectMat[2][0] = covar(0,2);
422 rectMat[1][1] = covar(1,1);
423 rectMat[1][2] = rectMat[2][1] = covar(1,2);
424 rectMat[2][2] = covar(2,2);
425
426 // Compute the Jacobian in meters. Don't deal with unit mismatch yet to preserve precision.
427 SpiceDouble J[3][3];
428 double zOverR = p_z->meters() / radius;
429 double r2 = radius*radius;
430 double denom = r2*radius*sqrt(1.0 - (zOverR*zOverR));
431 J[0][0] = -p_x->meters() * p_z->meters() / denom;
432 J[0][1] = -p_y->meters() * p_z->meters() / denom;
433 J[0][2] = (r2 - p_z->meters() * p_z->meters()) / denom;
434 J[1][0] = -p_y->meters() / (x2 + y2);
435 J[1][1] = p_x->meters() / (x2 + y2);
436 J[1][2] = 0.0;
437 J[2][0] = p_x->meters() / radius; // This row is unitless
438 J[2][1] = p_y->meters() / radius;
439 J[2][2] = p_z->meters() / radius;
440
441 if(!p_sphereCovar)
442 p_sphereCovar = new symmetric_matrix<double, upper>(3);
443
444 SpiceDouble mat[3][3];
445 mxm_c (J, rectMat, mat);
446 mxmt_c (mat, J, mat);
447 if (units == Kilometers) {
448 // Now take care of unit mismatch between rect matrix in km and Jacobian in m
449 (*p_sphereCovar)(0,0) = mat[0][0] * 1.0e6;
450 (*p_sphereCovar)(0,1) = mat[0][1] * 1.0e6;
451 (*p_sphereCovar)(0,2) = mat[0][2] * 1000.0;
452 (*p_sphereCovar)(1,1) = mat[1][1] * 1.0e6;
453
454 (*p_sphereCovar)(1,2) = mat[1][2] * 1000.0;
455 (*p_sphereCovar)(2,2) = mat[2][2];
456 }
457 else { // (units == Meters)
458 // Convert matrix lengths from m to km
459 (*p_sphereCovar)(0,0) = mat[0][0];
460 (*p_sphereCovar)(0,1) = mat[0][1];
461 (*p_sphereCovar)(0,2) = mat[0][2] / 1000.0;
462 (*p_sphereCovar)(1,1) = mat[1][1];
463 (*p_sphereCovar)(1,2) = mat[1][2] / 1000.0;
464 (*p_sphereCovar)(2,2) = mat[2][2] / 1.0e6;
465 }
466 }
467
468
481 const Longitude &lon,
482 const Distance &radius) {
483// Is error checking necessary or does Latitude, Longitude, and Distance handle it?????
484 if (!lat.isValid() || !lon.isValid() || !radius.isValid()) {
485 IString msg = "Latitude, longitude, or radius is an invalid value.";
486 throw IException(IException::User, msg, _FILEINFO_);
487 }
488
489 SpiceDouble dlat = (double) lat.radians();
490 SpiceDouble dlon = (double) lon.radians();
491 SpiceDouble dradius = radius.kilometers();
492
493 SpiceDouble rect[3];
494 latrec_c ( dradius, dlon, dlat, rect);
495
496 // Set local radius now since we have it and avoid calculating it later
497 p_localRadius = radius;
498
502 }
503
504
518 void SurfacePoint::SetSpherical(const Latitude &lat, const Longitude &lon,
519 const Distance &radius, const Angle &latSigma, const Angle &lonSigma,
520 const Distance &radiusSigma) {
521
522 SetSphericalPoint(lat, lon, radius);
523
524 if (latSigma.isValid() && lonSigma.isValid() && radiusSigma.isValid())
525 SetSphericalSigmas(latSigma, lonSigma, radiusSigma);
526 }
527
528
538 void SurfacePoint::SetSpherical(const Latitude &lat, const Longitude &lon,
539 const Distance &radius, const symmetric_matrix<double, upper> &covar) {
540 SetSphericalPoint(lat, lon, radius);
541 SetSphericalMatrix(covar);
542 }
543
544
554 const Longitude &lon, const Distance &radius) {
555 SetSphericalPoint(lat, lon, radius);
556 }
557
558
574 const Angle &lonSigma,
575 const Distance &radiusSigma) {
576 if (latSigma.isValid() && lonSigma.isValid() && radiusSigma.isValid()) {
577 symmetric_matrix<double,upper> covar(3);
578 covar.clear();
579
580 double sphericalCoordinate;
581 sphericalCoordinate = (double) latSigma.radians();
582 covar(0,0) = sphericalCoordinate*sphericalCoordinate;
583 sphericalCoordinate = (double) lonSigma.radians();
584 covar(1,1) = sphericalCoordinate*sphericalCoordinate;
585 sphericalCoordinate = (double) radiusSigma.meters();
586 covar(2,2) = sphericalCoordinate*sphericalCoordinate;
587
588 // Use default set units for radius of meters*meters
589 SetSphericalMatrix(covar);
590 }
591 else {
592 delete p_sphereCovar;
593 p_sphereCovar = NULL;
594
595 delete p_rectCovar;
596 p_rectCovar = NULL;
597 }
598 }
599
600
618 const Distance &lonSigma, const Distance &radiusSigma) {
619
620 if (!Valid()) {
621 IString msg = "Cannot set spherical sigmas on an invalid surface point";
622 throw IException(IException::Programmer, msg, _FILEINFO_);
623 }
624
626 Angle(MetersToLongitude(lonSigma.meters()), Angle::Radians), radiusSigma);
627 }
628
629
638 const symmetric_matrix<double, upper> & covar,
639 SurfacePoint::CoordUnits units) {
640
641 // Make sure the point is set first
642 if (!Valid()) {
643 IString msg = "A point must be set before a variance/covariance matrix "
644 "can be set.";
645 throw IException(IException::Programmer, msg, _FILEINFO_);
646 }
647
648 // Make sure the units are valid
649 if (units != Kilometers && units != Meters) {
650 IString msg = "Units must be Kilometers or Meters";
651 throw IException(IException::Programmer, msg, _FILEINFO_);
652 }
653
654 // Get the radius of the point in the same units as the input matrix when saving the input matrix
655 double radius = (double) GetLocalRadius().kilometers();
656
657 // Save the spherical matrix in km and km**2
658 if (p_sphereCovar) {
659 *p_sphereCovar = covar;
660 }
661 else {
662 p_sphereCovar = new symmetric_matrix<double, upper>(covar);
663 }
664
665 if (units == Meters) {
666 // Convert input matrix to km to store
667 (*p_sphereCovar)(0,0) = covar(0,0);
668 (*p_sphereCovar)(0,1) = covar(0,1);
669 (*p_sphereCovar)(0,2) = covar(0,2) / 1000.;
670 (*p_sphereCovar)(1,1) = covar(1,1);
671 (*p_sphereCovar)(1,2) = covar(1,2) / 1000.;
672 (*p_sphereCovar)(2,2) = covar(2,2) / 1.0e6;
673 radius = (double) GetLocalRadius().meters();
674 }
675
676 // ***TODO*** Consider using LinearAlgebra matrix multiply and avoid creating SpiceDouble matrix.
677 SpiceDouble sphereMat[3][3];
678
679 sphereMat[0][0] = covar(0,0);
680 sphereMat[0][1] = sphereMat[1][0] = covar(0,1);
681 sphereMat[0][2] = sphereMat[2][0] = covar(0,2);
682 sphereMat[1][1] = covar(1,1);
683 sphereMat[1][2] = sphereMat[2][1] = covar(1,2);
684 sphereMat[2][2] = covar(2,2);
685
686// std::cout<<"Ocovar = "<<sphereMat[0][0]<<" "<<sphereMat[0][1]<<" "<<sphereMat[0][2]<<std::endl
687// <<" "<<sphereMat[1][0]<<" "<<sphereMat[1][1]<<" "<<sphereMat[1][2]<<std::endl
688// <<" "<<sphereMat[2][0]<<" "<<sphereMat[2][1]<<" "<<sphereMat[2][2]<<std::endl;
689
690 // Get the lat/lon/radius of the point
691 double lat = (double) GetLatitude().radians();
692 double lon = (double) GetLongitude().radians();
693
694 // Compute the Jacobian in same units as input matrix.
695 SpiceDouble J[3][3];
696 double cosPhi = cos(lat);
697 double sinPhi = sin(lat);
698 double cosLamda = cos(lon);
699 double sinLamda = sin(lon);
700 double rcosPhi = radius*cosPhi;
701 double rsinPhi = radius*sinPhi;
702 J[0][0] = -rsinPhi * cosLamda;
703 J[0][1] = -rcosPhi * sinLamda;
704 J[0][2] = cosPhi * cosLamda;
705 J[1][0] = -rsinPhi * sinLamda;
706 J[1][1] = rcosPhi * cosLamda;
707 J[1][2] = cosPhi * sinLamda;
708 J[2][0] = rcosPhi;
709 J[2][1] = 0.0;
710 J[2][2] = sinPhi;
711
712 if(!p_rectCovar)
713 p_rectCovar = new symmetric_matrix<double, upper>(3);
714
715 SpiceDouble mat[3][3];
716 mxm_c (J, sphereMat, mat);
717 mxmt_c (mat, J, mat);
718
719 if (units == Kilometers) {
720 (*p_rectCovar)(0,0) = mat[0][0];
721 (*p_rectCovar)(0,1) = mat[0][1];
722 (*p_rectCovar)(0,2) = mat[0][2];
723 (*p_rectCovar)(1,1) = mat[1][1];
724 (*p_rectCovar)(1,2) = mat[1][2];
725 (*p_rectCovar)(2,2) = mat[2][2];
726 }
727 else { // (units == Meters)
728 // Convert to km
729 (*p_rectCovar)(0,0) = mat[0][0]/1.0e6;
730 (*p_rectCovar)(0,1) = mat[0][1]/1.0e6;
731 (*p_rectCovar)(0,2) = mat[0][2]/1.0e6;
732 (*p_rectCovar)(1,1) = mat[1][1]/1.0e6;
733 (*p_rectCovar)(1,2) = mat[1][2]/1.0e6;
734 (*p_rectCovar)(2,2) = mat[2][2]/1.0e6;
735 }
736// std::cout<<"Rcovar = "<<p_rectCovar(0,0)<<" "<<p_rectCovar(0,1)<<" "<<p_rectCovar(0,2)<<std::endl
737// <<" "<<p_rectCovar(1,0)<<" "<<p_rectCovar(1,1)<<" "<<p_rectCovar(1,2)<<std::endl
738// <<" "<<p_rectCovar(2,0)<<" "<<p_rectCovar(2,1)<<" "<<p_rectCovar(2,2)<<std::endl;
739 }
740
741
749 void SurfacePoint::SetMatrix(CoordinateType type, const symmetric_matrix<double, upper>& covar) {
750
751 switch (type) {
752 case Latitudinal:
753 SetSphericalMatrix(covar, SurfacePoint::Kilometers);
754 break;
755 case Rectangular:
756 SetRectangularMatrix(covar, SurfacePoint::Kilometers);
757 break;
758 default:
759 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
760 throw IException(IException::Programmer, msg, _FILEINFO_);
761 }
762 }
763
764
773 std::vector<double> SurfacePoint::Partial(CoordinateType type, CoordIndex index) {
774 std::vector<double> derivative(3);
775 switch (type) {
776 case Latitudinal:
777 derivative = LatitudinalDerivative(index);
778 break;
779 case Rectangular:
780 derivative = RectangularDerivative(index);
781 break;
782 default:
783 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
784 throw IException(IException::Programmer, msg, _FILEINFO_);
785 }
786 return derivative;
787 }
788
789
798 std::vector<double> SurfacePoint::LatitudinalDerivative(CoordIndex index) {
799 std::vector<double> derivative(3);
800 double rlat = GetLatitude().radians();
801 double rlon = GetLongitude().radians();
802 double sinLon = sin(rlon);
803 double cosLon = cos(rlon);
804 double sinLat = sin(rlat);
805 double cosLat = cos(rlat);
806 double radkm = GetLocalRadius().kilometers();
807
808 switch (index) {
809 case One:
810 derivative[0] = -radkm * sinLat * cosLon;
811 derivative[1] = -radkm * sinLon * sinLat;
812 derivative[2] = radkm * cosLat;
813 break;
814 case Two:
815 derivative[0] = -radkm * cosLat * sinLon;
816 derivative[1] = radkm * cosLat * cosLon;
817 derivative[2] = 0.0;
818 break;
819 case Three:
820 derivative[0] = cosLon * cosLat;
821 derivative[1] = sinLon * cosLat;
822 derivative[2] = sinLat;
823 break;
824 default:
825 IString msg = "Invalid coordinate index (must be less than 3).";
826 throw IException(IException::Programmer, msg, _FILEINFO_);
827 }
828 return derivative;
829 }
830
831
840 std::vector<double> SurfacePoint::RectangularDerivative(CoordIndex index) {
841 std::vector<double> derivative(3,0.0);
842
843 switch (index) {
844 case One:
845 derivative[0] = 1.;
846 break;
847 case Two:
848 derivative[1] = 1.;
849 break;
850 case Three:
851 derivative[2] = 1.;
852 break;
853 default:
854 IString msg = "Invalid coordinate index (must be less than 3).";
855 throw IException(IException::Programmer, msg, _FILEINFO_);
856 }
857 return derivative;
858 }
859
860
870 void SurfacePoint::ToNaifArray(double naifOutput[3]) const {
871 if(Valid()) {
872 naifOutput[0] = p_x->kilometers();
873 naifOutput[1] = p_y->kilometers();
874 naifOutput[2] = p_z->kilometers();
875 }
876 else {
877 IString msg = "Cannot convert an invalid surface point to a naif array";
878 throw IException(IException::Programmer, msg, _FILEINFO_);
879 }
880 }
881
891 void SurfacePoint::FromNaifArray(const double naifValues[3]) {
892 if(p_x && p_y && p_z) {
893 p_x->setKilometers(naifValues[0]);
894 p_y->setKilometers(naifValues[1]);
895 p_z->setKilometers(naifValues[2]);
896 }
897 else {
898 p_x = new Displacement(naifValues[0], Displacement::Kilometers);
899 p_y = new Displacement(naifValues[1], Displacement::Kilometers);
900 p_z = new Displacement(naifValues[2], Displacement::Kilometers);
901 }
902
904 p_localRadius = GetLocalRadius();
905 }
906
907
915
916 if (!radius.isValid()) {
917 IString msg = "Radius value must be a valid Distance.";
918 throw IException(IException::Programmer, msg, _FILEINFO_);
919 }
920
921 if (!p_x || !p_y || !p_z || !p_x->isValid() || !p_y->isValid() ||
922 !p_z->isValid()) {
923 IString msg = "In order to reset the local radius, a Surface Point must "
924 "already be set.";
925 throw IException(IException::Programmer, msg, _FILEINFO_);
926 }
927
928 // Set latitudinal coordinates
929 SpiceDouble lat = (double) GetLatitude().radians();
930 SpiceDouble lon = (double) GetLongitude().radians();
931
932 p_localRadius = radius;
933
934 // Set rectangular coordinates
935 SpiceDouble rect[3];
936 latrec_c ((SpiceDouble) radius.kilometers(), lon, lat, rect);
937 p_x->setKilometers(rect[0]);
938 p_y->setKilometers(rect[1]);
939 p_z->setKilometers(rect[2]);
940
941 // TODO What should be done to the variance/covariance matrix when the
942 // radius is reset??? With Bundle updates will this functionality be
943 // obsolete??? Ask Ken
944 }
945
946
947 bool SurfacePoint::Valid() const {
948 static const Displacement zero(0, Displacement::Meters);
949 return p_x && p_y && p_z && p_x->isValid() && p_y->isValid() && p_z->isValid() &&
950 (*p_x != zero || *p_y != zero || *p_z != zero);
951 }
952
953
962 double SurfacePoint::GetCoord(CoordinateType type, CoordIndex index, CoordUnits units) {
963 // TODO *** Is there a better way to satisfy the compiler that a value will be initialized?
964 // Don't the enums take care of preventing any other possibilities? no
965 double value = 0.;
966
967 switch (type) {
968
969 case Latitudinal:
970
971 switch (index) {
972
973 case One: // Latitude
974 value = LatToDouble(GetLatitude(), units);
975 break;
976
977 case Two: // Longitude
978 value = LonToDouble(GetLongitude(), units);
979 break;
980
981 case Three: // LocalRadius
982 value = DistanceToDouble(GetLocalRadius(), units);
983 break;
984
985 default:
986 IString msg = "Invalid coordinate index (must be less than 3).";
987 throw IException(IException::Programmer, msg, _FILEINFO_);
988 }
989 break;
990
991 case Rectangular:
992
993 switch (index) {
994
995 case One: // X
996 value = DisplacementToDouble(GetX(), units);
997 break;
998
999 case Two: // Y
1000 value = DisplacementToDouble(GetY(), units);
1001 break;
1002
1003 case Three: // Z
1004 value = DisplacementToDouble(GetZ(), units);
1005 break;
1006
1007 default:
1008 IString msg = "Invalid coordinate index enum (must be less than 3).";
1009 throw IException(IException::Programmer, msg, _FILEINFO_);
1010 }
1011 break;
1012
1013 default:
1014 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
1015 throw IException(IException::Programmer, msg, _FILEINFO_);
1016 }
1017 return value;
1018 }
1019
1020
1029 double SurfacePoint::GetSigma(CoordinateType type, CoordIndex index, CoordUnits units) {
1030 double value = 0; // See first TODO in GetCoord
1031
1032 switch (type) {
1033
1034 case Latitudinal:
1035
1036 switch (index) {
1037
1038 case One:
1039 value = DistanceToDouble(GetLatSigmaDistance(), units);
1040 break;
1041
1042 case Two:
1043 value = DistanceToDouble(GetLonSigmaDistance(), units);
1044 break;
1045
1046 case Three:
1047 value = DistanceToDouble(GetLocalRadiusSigma(), units);
1048 break;
1049
1050 default:
1051 IString msg = "Invalid coordinate index (must be less than 3).";
1052 throw IException(IException::Programmer, msg, _FILEINFO_);
1053 }
1054 break;
1055
1056 case Rectangular:
1057
1058 switch (index) {
1059
1060 case One:
1061 value = DistanceToDouble(GetXSigma(), units);
1062 break;
1063
1064 case Two:
1065 value = DistanceToDouble(GetYSigma(), units);
1066 break;
1067
1068 case Three:
1069 value = DistanceToDouble(GetZSigma(), units);
1070 break;
1071
1072 default:
1073 IString msg = "Invalid coordinate index (must be less than 3).";
1074 throw IException(IException::Programmer, msg, _FILEINFO_);
1075 }
1076 break;
1077
1078 default:
1079 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
1080 throw IException(IException::Programmer, msg, _FILEINFO_);
1081 }
1082
1083 return value;
1084 }
1085
1086
1095 CoordIndex index) {
1096 Distance dist = Distance(); // See first TODO in GetCoord
1097
1098 switch (type) {
1099
1100 case Latitudinal:
1101
1102 switch (index) {
1103
1104 case One:
1105 dist = GetLatSigmaDistance();
1106 break;
1107
1108 case Two:
1109 dist = GetLonSigmaDistance();
1110 break;
1111
1112 case Three:
1113 dist = GetLocalRadiusSigma();
1114 break;
1115
1116 default:
1117 IString msg = "Invalid coordinate index (must be less than 3).";
1118 throw IException(IException::Programmer, msg, _FILEINFO_);
1119 }
1120 break;
1121
1122 case Rectangular:
1123
1124 switch (index) {
1125
1126 case One:
1127 dist = GetXSigma();
1128 break;
1129
1130 case Two:
1131 dist = GetYSigma();
1132 break;
1133
1134 case Three:
1135 dist = GetZSigma();
1136 break;
1137
1138 default:
1139 IString msg = "Invalid coordinate index (must be less than 3).";
1140 throw IException(IException::Programmer, msg, _FILEINFO_);
1141 }
1142 break;
1143
1144 default:
1145 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
1146 throw IException(IException::Programmer, msg, _FILEINFO_);
1147 }
1148
1149 return dist;
1150 }
1151
1152
1160 double SurfacePoint::DisplacementToDouble(Displacement disp, CoordUnits units) {
1161 double value;
1162
1163 switch (units) {
1164
1165 case Meters:
1166 value = disp.meters();
1167 break;
1168
1169 case Kilometers:
1170 value = disp.kilometers();
1171 break;
1172
1173 default:
1174 IString msg = "Unrecognized unit for a Displacement (must be meters or kilometers).";
1175 throw IException(IException::Programmer, msg, _FILEINFO_);
1176 }
1177 return value;
1178 }
1179
1180
1188 double SurfacePoint::DistanceToDouble(Distance dist, CoordUnits units) {
1189 double value;
1190
1191 switch (units) {
1192
1193 case Meters:
1194 value = dist.meters();
1195 break;
1196
1197 case Kilometers:
1198 value = dist.kilometers();
1199 break;
1200
1201 default:
1202 IString msg = "Unrecognized unit for a Distance (not meters or kilometers).";
1203 throw IException(IException::Programmer, msg, _FILEINFO_);
1204 }
1205 return value;
1206 }
1207
1208
1216 double SurfacePoint::LatToDouble(Latitude lat, CoordUnits units) {
1217 double value;
1218
1219 switch (units) {
1220
1221 case Radians:
1222 value = GetLatitude().radians();
1223 break;
1224
1225 case Degrees:
1226 value = GetLatitude().degrees();
1227 break;
1228
1229 default:
1230 IString msg = "Invalid unit for latitude (must be Radians or Degrees).";
1231 throw IException(IException::Programmer, msg, _FILEINFO_);
1232 }
1233 return value;
1234 }
1235
1236
1248 double SurfacePoint::MetersToLatitude(double latLength) {
1249 if (Valid() && !IsSpecial(latLength)) {
1250 // Convert angle measure in meters to radians relative to latitude of SurfacePoint.
1251 return latLength / GetLocalRadius().meters();
1252 }
1253 else {
1254 // Return Null to be consistent with the bundle classes
1255 return Isis::Null;
1256 }
1257 }
1258
1259
1276 double SurfacePoint::MetersToLongitude(double deltaLonMeters) {
1277
1278 if (Valid() && !IsSpecial(deltaLonMeters)) {
1279 double convFactor = cos((double)GetLatitude().radians());
1280 double deltaLonRadians;
1281
1282 // Convert angle displacement to radians relative to longitude of SurfacePoint.
1283 if (convFactor > 1.0e-50) {
1284 deltaLonRadians = deltaLonMeters / (convFactor*GetLocalRadius().meters());
1285 }
1286 else {
1287 // Brent Archinal suggested setting sigma to pi in the case of a point near the pole
1288 deltaLonRadians = PI;
1289 }
1290 return deltaLonRadians;
1291 }
1292 else {
1293 // Return Null to be consistent with the bundle classes
1294 return Isis::Null;
1295 }
1296 }
1297
1298
1310 double SurfacePoint::LatitudeToMeters(double relativeLat) const {
1311 // Returns are consistent with the bundle classes
1312 if (relativeLat == 0.) {
1313 return 0.;
1314 }
1315 else if (Valid() && !IsSpecial(relativeLat) && GetLocalRadius().isValid() ) {
1316 return relativeLat * GetLocalRadius().meters();
1317 }
1318 else {
1319 return Isis::Null;
1320 }
1321 }
1322
1323
1336 double SurfacePoint::LongitudeToMeters(double deltaLonRadians) const{
1337 // Returns are consistent with the bundle classes
1338 double deltaLonMeters = Isis::Null;
1339
1340 if (deltaLonRadians == 0.) {
1341 deltaLonMeters = 0.;
1342 }
1343 else if (Valid() && !IsSpecial(deltaLonRadians) && GetLocalRadius().isValid() ) {
1344 // Convert from radians to meters and return
1345 double scalingRadius = cos(GetLatitude().radians()) * GetLocalRadius().meters();
1346 deltaLonMeters = deltaLonRadians * scalingRadius;
1347 }
1348
1349 return deltaLonMeters;
1350 }
1351
1352
1364 if (type.compare("LATITUDINAL", Qt::CaseInsensitive) == 0) {
1366 }
1367 else if (type.compare("RECTANGULAR", Qt::CaseInsensitive) == 0) {
1369 }
1370 else {
1372 "Unknown coordinate type for a SurfacePoint [" + type + "].",
1373 _FILEINFO_);
1374 }
1375 }
1376
1377
1391 switch (type) {
1392
1393 case Latitudinal:
1394 return "Latitudinal";
1395 break;
1396
1397 case Rectangular:
1398 return "Rectangular";
1399 break;
1400
1401 default:
1402 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
1403 throw IException(IException::Programmer, msg, _FILEINFO_);
1404 }
1405 }
1406
1407
1415 double SurfacePoint::LonToDouble(Longitude lon, CoordUnits units) {
1416 double value;
1417
1418 switch (units) {
1419
1420 case Radians:
1421 value = GetLongitude().radians();
1422 break;
1423
1424 case Degrees:
1425 value = GetLongitude().degrees();
1426 break;
1427
1428 default:
1429 IString msg = "Invalid unit for longitude (must be Radians of Degrees).";
1430 throw IException(IException::Programmer, msg, _FILEINFO_);
1431 }
1432 return value;
1433 }
1434
1435
1436 Displacement SurfacePoint::GetX() const {
1437 if(!p_x) return Displacement();
1438
1439 return *p_x;
1440 }
1441
1442
1443 Displacement SurfacePoint::GetY() const {
1444 if(!p_y) return Displacement();
1445
1446 return *p_y;
1447 }
1448
1449
1450 Displacement SurfacePoint::GetZ() const {
1451 if(!p_z) return Displacement();
1452
1453 return *p_z;
1454 }
1455
1456
1457 Distance SurfacePoint::GetXSigma() const {
1458 if(!p_rectCovar) return Distance();
1459
1460 return Distance(sqrt((*p_rectCovar)(0, 0)), Distance::Kilometers);
1461 }
1462
1463
1464 Distance SurfacePoint::GetYSigma() const {
1465 if(!p_rectCovar) return Distance();
1466
1467 return Distance(sqrt((*p_rectCovar)(1, 1)), Distance::Kilometers);
1468 }
1469
1470
1471 Distance SurfacePoint::GetZSigma() const {
1472 if(!p_rectCovar) return Distance();
1473
1474 return Distance(sqrt((*p_rectCovar)(2, 2)), Distance::Kilometers);
1475 }
1476
1477
1491 double SurfacePoint::GetWeight(CoordinateType type, CoordIndex index) {
1492 double value = 0; // See first TODO in GetCoord
1493
1494 switch (type) {
1495
1496 case Latitudinal:
1497
1498 switch (index) {
1499
1500 case One:
1501 value = GetLatWeight(); // 1/radians**2
1502 break;
1503
1504 case Two:
1505 value = GetLonWeight(); // 1/radians**2
1506 break;
1507
1508 case Three:
1509 value = GetLocalRadiusWeight(); // 1/km**2
1510 break;
1511
1512 default:
1513 IString msg = "Invalid coordinate index (must be less than 3).";
1514 throw IException(IException::Programmer, msg, _FILEINFO_);
1515 }
1516 break;
1517
1518 case Rectangular:
1519
1520 switch (index) {
1521
1522 case One:
1523 value = GetXWeight(); // 1/km**2
1524 break;
1525
1526 case Two:
1527 value = GetYWeight(); // 1/km**2
1528 break;
1529
1530 case Three:
1531 value = GetZWeight(); // 1/km**2
1532 break;
1533
1534 default:
1535 IString msg = "Invalid coordinate index (must be less than 3).";
1536 throw IException(IException::Programmer, msg, _FILEINFO_);
1537 }
1538 break;
1539
1540 default:
1541 IString msg ="Unknown surface point coordinate type enum [" + toString(type) + "]." ;
1542 throw IException(IException::Programmer, msg, _FILEINFO_);
1543 }
1544 return value;
1545 }
1546
1547
1554
1555 double dXSigma = GetXSigma().kilometers();
1556
1557 if (dXSigma <= 0.0 ) {
1558 IString msg = "SurfacePoint::GetXWeight(): Sigma <= 0.0";
1559 throw IException(IException::Programmer, msg, _FILEINFO_);
1560 }
1561
1562 return 1.0/(dXSigma*dXSigma);
1563 }
1564
1565
1572
1573 double dYSigma = GetYSigma().kilometers();
1574
1575 if (dYSigma <= 0.0 ) {
1576 IString msg = "SurfacePoint::GetYWeight(): Sigma <= 0.0";
1577 throw IException(IException::Programmer, msg, _FILEINFO_);
1578 }
1579
1580 return 1.0/(dYSigma*dYSigma);
1581 }
1582
1583
1590
1591 double dZSigma = GetZSigma().kilometers();
1592
1593 if (dZSigma <= 0.0 ) {
1594 IString msg = "SurfacePoint::GetZWeight(): Sigma <= 0.0";
1595 throw IException(IException::Programmer, msg, _FILEINFO_);
1596 }
1597
1598 return 1.0/(dZSigma*dZSigma);
1599 }
1600
1601
1602 symmetric_matrix<double, upper> SurfacePoint::GetRectangularMatrix
1603 (SurfacePoint::CoordUnits units) const {
1604 if(!p_rectCovar) {
1605 symmetric_matrix<double, upper> tmp(3);
1606 tmp.clear();
1607 return tmp;
1608 }
1609
1610 // Make sure the units are valid
1611 if (units != Kilometers && units != Meters) {
1612 IString msg = "Units must be Kilometers or Meters";
1613 throw IException(IException::Programmer, msg, _FILEINFO_);
1614 }
1615
1616 symmetric_matrix<double, upper> covar(3);
1617 covar.clear();
1618
1619 switch (units) {
1620
1621 case Meters:
1622 // Convert member matrix to Meters to return
1623 covar(0,0) = (*p_rectCovar)(0,0)*1.0e6;
1624 covar(0,1) = (*p_rectCovar)(0,1)*1.0e6;
1625 covar(0,2) = (*p_rectCovar)(0,2)*1.0e6;
1626 covar(1,1) = (*p_rectCovar)(1,1)*1.0e6;
1627 covar(1,2) = (*p_rectCovar)(1,2)*1.0e6;
1628 covar(2,2) = (*p_rectCovar)(2,2)*1.0e6;
1629 return covar;
1630 break;
1631
1632 case Kilometers:
1633 return *p_rectCovar;
1634 break;
1635
1636 default:
1637 IString msg = "Unrecognized unit for a Displacement (must be meters or kilometers).";
1638 throw IException(IException::Programmer, msg, _FILEINFO_);
1639 }
1640
1641 return *p_rectCovar;
1642 }
1643
1644
1645 Angle SurfacePoint::GetLatSigma() const {
1646 if(!p_sphereCovar)
1647 return Angle();
1648
1649 return Angle(sqrt((*p_sphereCovar)(0, 0)), Angle::Radians);
1650 }
1651
1652
1653 Angle SurfacePoint::GetLonSigma() const {
1654 if(!p_sphereCovar)
1655 return Angle();
1656
1657 return Angle(sqrt((*p_sphereCovar)(1, 1)), Angle::Radians);
1658 }
1659
1660
1666 if (!Valid())
1667 return Latitude();
1668
1669 // TODO Scale for accuracy with coordinate of largest magnitude
1670 double x = p_x->meters();
1671 double y = p_y->meters();
1672 double z = p_z->meters();
1673
1674 if (x != 0. || y != 0. || z != 0.)
1675 return Latitude(atan2(z, sqrt(x*x + y*y) ), Angle::Radians);
1676 else
1677 return Latitude();
1678 }
1679
1680
1686 if (!Valid())
1687 return Longitude();
1688
1689 double x = p_x->meters();
1690 double y = p_y->meters();
1691
1692 if(x == 0.0 && y == 0.0) {
1693 return Longitude(0, Angle::Radians);
1694 }
1695
1696 double lon = atan2(y, x);
1697 if (lon < 0) {
1698 lon += 2 * PI;
1699 }
1700
1701 return Longitude(lon, Angle::Radians);
1702 }
1703
1704
1710 static const Displacement zero(0, Displacement::Meters);
1711 if (Valid()) {
1712 double x = p_x->meters();
1713 double y = p_y->meters();
1714 double z = p_z->meters();
1715
1716 p_localRadius = Distance(sqrt(x*x + y*y + z*z), Distance::Meters);
1717 }
1718 else if (*p_x == zero && *p_y == zero && *p_z == zero) { // for backwards compatability
1719 p_localRadius = Distance(0., Distance::Meters);
1720 }
1721 else { // Invalid point
1722 IString msg = "SurfacePoint::Can't compute local radius on invalid point";
1723 throw IException(IException::Programmer, msg, _FILEINFO_);
1724 }
1725 }
1726
1727
1733 return p_localRadius;
1734 }
1735
1736
1748 double d = LatitudeToMeters(GetLatSigma().radians());
1749
1750 if (d > 1.0e-50) {
1751 return Distance(d, Distance::Meters);
1752 }
1753 else {
1754 return Distance();
1755 }
1756 }
1757
1758
1764 // return lonSigmaDistance;
1765 double d = LongitudeToMeters(GetLonSigma().radians());
1766// TODO What do we do when the scaling radius is 0 (at the pole)?
1767 // if (d > DBL_EPSILON) {
1768 return Distance(d, Distance::Meters);
1769 // }
1770 // else { // Too close to the pole
1771 // return Distance();
1772 // }
1773 }
1774
1775
1776 Distance SurfacePoint::GetLocalRadiusSigma() const {
1777 if(!p_sphereCovar)
1778 return Distance();
1779
1780 return Distance(sqrt((*p_sphereCovar)(2, 2)), Distance::Kilometers);
1781 }
1782
1783
1784 symmetric_matrix<double, upper> SurfacePoint::GetSphericalMatrix
1785 (SurfacePoint::CoordUnits units) const {
1786 if(!p_sphereCovar) {
1787 symmetric_matrix<double, upper> tmp(3);
1788 tmp.clear();
1789 return tmp;
1790 }
1791
1792 // Make sure the units are valid
1793 if (units != Kilometers && units != Meters) {
1794 IString msg = "Units must be Kilometers or Meters";
1795 throw IException(IException::Programmer, msg, _FILEINFO_);
1796 }
1797
1798 symmetric_matrix<double, upper> covar(3);
1799 covar.clear();
1800
1801 switch (units) {
1802
1803 case Meters:
1804 // Convert member matrix to Meters to return
1805 covar(0,0) = (*p_sphereCovar)(0,0);
1806 covar(0,1) = (*p_sphereCovar)(0,1);
1807 covar(0,2) = (*p_sphereCovar)(0,2)*1000.;
1808 covar(1,1) = (*p_sphereCovar)(1,1);
1809 covar(1,2) = (*p_sphereCovar)(1,2)*1000.;
1810 covar(2,2) = (*p_sphereCovar)(2,2)*1.0e6;
1811 return covar;
1812 break;
1813
1814 case Kilometers:
1815 return *p_sphereCovar;
1816 break;
1817
1818 default:
1819 IString msg = "Unrecognized unit for a Displacement (must be meters or kilometers).";
1820 throw IException(IException::Programmer, msg, _FILEINFO_);
1821 }
1822
1823 return *p_sphereCovar;
1824 }
1825
1826
1833 double dlatSigma = GetLatSigma().radians();
1834
1835 if( dlatSigma <= 0.0 ) {
1836 IString msg = "SurfacePoint::GetLatWeight(): Sigma <= 0.0";
1837 throw IException(IException::Programmer, msg, _FILEINFO_);
1838 }
1839
1840 return 1.0/(dlatSigma*dlatSigma);
1841 }
1842
1849 double dlonSigma = GetLonSigma().radians();
1850
1851 if( dlonSigma <= 0.0 ) {
1852 IString msg = "SurfacePoint::GetLonWeight(): Sigma <= 0.0";
1853 throw IException(IException::Programmer, msg, _FILEINFO_);
1854 }
1855
1856 return 1.0/(dlonSigma*dlonSigma);
1857 }
1858
1865
1866 double dlocalRadiusSigma = GetLocalRadiusSigma().kilometers();
1867
1868 if (dlocalRadiusSigma <= 0.0 ) {
1869 IString msg = "SurfacePoint::GetRadWeight(): Sigma <= 0.0";
1870 throw IException(IException::Programmer, msg, _FILEINFO_);
1871 }
1872
1873 return 1.0/(dlocalRadiusSigma*dlocalRadiusSigma);
1874 }
1875
1881
1882 if(!Valid() || !other.Valid())
1883 return Distance();
1884
1885 return GetDistanceToPoint(other,
1886 ((GetLocalRadius() + other.GetLocalRadius()) / 2.0));
1887 }
1888
1889
1898 const Distance &sphereRadius) const {
1899 if(!Valid() || !other.Valid())
1900 return Distance();
1901
1902 // Convert lat/lon values to radians
1903 const Angle &latitude = GetLatitude();
1904 const Angle &longitude = GetLongitude();
1905 const Angle &otherLatitude = other.GetLatitude();
1906 const Angle &otherLongitude = other.GetLongitude();
1907
1908 // The harvestine method:
1909 // http://en.wikipedia.org/wiki/Haversine_formula
1910 Angle deltaLat = latitude - otherLatitude;
1911 Angle deltaLon = longitude - otherLongitude;
1912
1913 double haversinLat = sin(deltaLat.radians() / 2.0);
1914 haversinLat *= haversinLat;
1915
1916 double haversinLon = sin(deltaLon.radians() / 2.0);
1917 haversinLon *= haversinLon;
1918
1919 double a = haversinLat + cos(latitude.radians()) *
1920 cos(otherLatitude.radians()) *
1921 haversinLon;
1922
1923 double c = 2 * atan(sqrt(a) / sqrt(1 - a));
1924
1925 return sphereRadius * c;
1926 }
1927
1928
1929 bool SurfacePoint::operator==(const SurfacePoint &other) const {
1930 bool equal = true;
1931
1932 if(equal && p_x && p_y && p_z && other.p_x && other.p_y && other.p_z) {
1933 equal = equal && *p_x == *other.p_x;
1934 equal = equal && *p_y == *other.p_y;
1935 equal = equal && *p_z == *other.p_z;
1936 }
1937 else if(equal) {
1938 equal = equal && p_x == NULL && other.p_x == NULL;
1939 equal = equal && p_y == NULL && other.p_y == NULL;
1940 equal = equal && p_z == NULL && other.p_z == NULL;
1941 }
1942
1943 if(equal && p_rectCovar) {
1944 equal = equal && (*p_rectCovar)(0, 0) == (*other.p_rectCovar)(0, 0);
1945 equal = equal && (*p_rectCovar)(0, 1) == (*other.p_rectCovar)(0, 1);
1946 equal = equal && (*p_rectCovar)(0, 2) == (*other.p_rectCovar)(0, 2);
1947 equal = equal && (*p_rectCovar)(1, 1) == (*other.p_rectCovar)(1, 1);
1948 equal = equal && (*p_rectCovar)(1, 2) == (*other.p_rectCovar)(1, 2);
1949 equal = equal && (*p_rectCovar)(2, 2) == (*other.p_rectCovar)(2, 2);
1950 }
1951 else if(equal) {
1952 equal = equal && p_rectCovar == NULL && other.p_rectCovar == NULL;
1953 }
1954
1955 if(equal && p_sphereCovar) {
1956 equal = equal && (*p_sphereCovar)(0, 0) == (*other.p_sphereCovar)(0, 0);
1957 equal = equal && (*p_sphereCovar)(0, 1) == (*other.p_sphereCovar)(0, 1);
1958 equal = equal && (*p_sphereCovar)(0, 2) == (*other.p_sphereCovar)(0, 2);
1959 equal = equal && (*p_sphereCovar)(1, 1) == (*other.p_sphereCovar)(1, 1);
1960 equal = equal && (*p_sphereCovar)(1, 2) == (*other.p_sphereCovar)(1, 2);
1961 equal = equal && (*p_sphereCovar)(2, 2) == (*other.p_sphereCovar)(2, 2);
1962 }
1963 else if(equal) {
1964 equal = equal && p_sphereCovar == NULL && other.p_sphereCovar == NULL;
1965 }
1966
1967 return equal;
1968 }
1969
1970 SurfacePoint &SurfacePoint::operator=(const SurfacePoint &other) {
1971 // The lazy way of doing this (free all memory and copy) is too expensive
1972 // in the default case!
1973 if(p_x && other.p_x &&
1974 p_y && other.p_y &&
1975 p_z && other.p_z &&
1976 !p_rectCovar && !other.p_rectCovar &&
1977 !p_sphereCovar && !other.p_sphereCovar) {
1978 *p_x = *other.p_x;
1979 *p_y = *other.p_y;
1980 *p_z = *other.p_z;
1981 }
1982 else {
1983 FreeAllocatedMemory();
1984
1985 if(other.p_x) {
1986 p_x = new Displacement(*other.p_x);
1987 }
1988
1989 if(other.p_y) {
1990 p_y = new Displacement(*other.p_y);
1991 }
1992
1993 if(other.p_z) {
1994 p_z = new Displacement(*other.p_z);
1995 }
1996
1997 if(other.p_rectCovar) {
1998 p_rectCovar = new symmetric_matrix<double, upper>(*other.p_rectCovar);
1999 }
2000
2001 if(other.p_sphereCovar) {
2002 p_sphereCovar = new symmetric_matrix<double, upper>(*other.p_sphereCovar);
2003 }
2004 }
2005
2006 // Finally initialize local radius to avoid using a previous value
2007 p_localRadius = other.GetLocalRadius();
2008
2009 return *this;
2010 }
2011
2012 void SurfacePoint::FreeAllocatedMemory() {
2013 if(p_x) {
2014 delete p_x;
2015 p_x = NULL;
2016 }
2017
2018 if(p_y) {
2019 delete p_y;
2020 p_y = NULL;
2021 }
2022
2023 if(p_z) {
2024 delete p_z;
2025 p_z = NULL;
2026 }
2027
2028 if(p_rectCovar) {
2029 delete p_rectCovar;
2030 p_rectCovar = NULL;
2031 }
2032
2033 if(p_sphereCovar) {
2034 delete p_sphereCovar;
2035 p_sphereCovar = NULL;
2036 }
2037 }
2038}
Defines an angle and provides unit conversions.
Definition Angle.h:45
double radians() const
Convert an angle to a double.
Definition Angle.h:226
double degrees() const
Get the angle in units of Degrees.
Definition Angle.h:232
@ Radians
Radians are generally used in mathematical equations, 0-2*PI is one circle, however these are more di...
Definition Angle.h:63
Displacement is a signed length, usually in meters.
double kilometers() const
Get the displacement in kilometers.
bool isValid() const
Test if this displacement has been initialized or not.
@ Kilometers
The distance is being specified in kilometers.
@ Meters
The distance is being specified in meters.
void setKilometers(double displacementInKilometers)
Set the displacement in kilometers.
double meters() const
Get the displacement in meters.
Distance measurement, usually in meters.
Definition Distance.h:34
bool isValid() const
Test if this distance has been initialized or not.
Definition Distance.cpp:192
double kilometers() const
Get the distance in kilometers.
Definition Distance.cpp:106
@ Kilometers
The distance is being specified in kilometers.
Definition Distance.h:45
@ Meters
The distance is being specified in meters.
Definition Distance.h:43
double meters() const
Get the distance in meters.
Definition Distance.cpp:85
Isis exception class.
Definition IException.h:91
@ User
A type of error that could only have occurred due to a mistake on the user's part (e....
Definition IException.h:126
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
Adds specific functionality to C++ strings.
Definition IString.h:165
This class is designed to encapsulate the concept of a Latitude.
Definition Latitude.h:51
This class is designed to encapsulate the concept of a Longitude.
Definition Longitude.h:40
This class defines a body-fixed surface point.
void SetRectangularCoordinates(const Displacement &x, const Displacement &y, const Displacement &z)
Set surface point in rectangular coordinates.
CoordinateType
Defines the coordinate typ, units, and coordinate index for some of the output methods.
@ Latitudinal
Planetocentric latitudinal (lat/lon/rad) coordinates.
@ Rectangular
Body-fixed rectangular x/y/z coordinates.
std::vector< double > LatitudinalDerivative(CoordIndex index)
Compute partial derivative of the conversion of the latitudinal coordinates to body-fixed rectangular...
void SetMatrix(CoordinateType type, const boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > &covar)
Set the covariance matrix.
void SetSphericalCoordinates(const Latitude &lat, const Longitude &lon, const Distance &radius)
Update spherical coordinates (lat/lon/radius)
static QString coordinateTypeToString(CoordinateType type)
Converts the given SurfacePoint::CoordinateType enumeration to a string.
double LatToDouble(Latitude lat, CoordUnits units)
This method returns a double version of a Latitude in the specified units.
Latitude GetLatitude() const
Return the body-fixed latitude for the surface point.
void SetRectangularSigmas(const Distance &xSigma, const Distance &ySigma, const Distance &zSigma)
Set surface point and sigmas in rectangular coordinates and convert to planetocentric.
double GetLonWeight() const
Return longitude weight for bundle adjustment Units are 1/(radians)^2.
double DistanceToDouble(Distance dist, CoordUnits units)
This method returns a double version of a Distance in the specified units.
void SetSphericalSigmasDistance(const Distance &latSigma, const Distance &lonSigma, const Distance &radiusSigma)
Set the spherical sigmas (in Distance units) into the spherical variance/covariance matrix.
Distance GetLonSigmaDistance() const
Return the longitude sigma in meters.
double GetWeight(CoordinateType type, CoordIndex index)
This method returns the weight of a SurfacePoint coordinate Note: At this time a units argument is no...
void ToNaifArray(double naifOutput[3]) const
A naif array is a c-style array of size 3.
void SetRectangularPoint(const Displacement &x, const Displacement &y, const Displacement &z)
This is a private method to set a surface point in rectangular, body-fixed coordinates.
boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > * p_rectCovar
3x3 upper triangular covariance matrix rectangular coordinates
void ComputeLocalRadius()
Compute the local radius of the surface point.
Distance GetDistanceToPoint(const SurfacePoint &other) const
Computes and returns the distance between two surface points.
void SetSphericalSigmas(const Angle &latSigma, const Angle &lonSigma, const Distance &radiusSigma)
Set the spherical sigmas into the spherical variance/covariance matrix in diagonal units of radians^2...
void SetSpherical(const Latitude &lat, const Longitude &lon, const Distance &radius, const Angle &latSigma=Angle(), const Angle &lonSigma=Angle(), const Distance &radiusSigma=Distance())
Set surface point and covariance matrix in planetocentric coordinates and convert to rectangular (Lat...
double MetersToLongitude(double lonLength)
This method returns an angular measure in radians of a distance in the direction of and relative to t...
boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > * p_sphereCovar
3x3 upper triangular covariance matrix ocentric coordinates
~SurfacePoint()
Destroys a SurfacePoint object/.
void InitCovariance()
Initialize the variance/covariance matrices.
void SetRectangular(const Displacement &x, const Displacement &y, const Displacement &z, const Distance &xSigma=Distance(), const Distance &ySigma=Distance(), const Distance &zSigma=Distance())
Set surface point in rectangular body-fixed coordinates wtih optional sigmas.
double LongitudeToMeters(double longitude) const
This method returns a length in meters version of a delta longitude angle in radians relative to the ...
double GetZWeight() const
Return Z weight for bundle adjustment Units are 1/(kilometers)^2.
Longitude GetLongitude() const
Return the body-fixed longitude for the surface point.
void SetSphericalMatrix(const boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > &covar, SurfacePoint::CoordUnits units=SurfacePoint::Meters)
Set spherical covariance matrix.
static CoordinateType stringToCoordinateType(QString type)
This method converts the given string value to a SurfacePoint::CoordinateType enumeration.
std::vector< double > RectangularDerivative(CoordIndex index)
Compute partial derivative of the body-fixed rectangular coordinates with respect to the indicated co...
double GetCoord(CoordinateType type, CoordIndex index, CoordUnits units)
This method returns a coordinate of a SurfacePoint.
double GetLocalRadiusWeight() const
Return radius weight for bundle adjustment Units are 1/(kilometers)^2.
double LatitudeToMeters(double latitude) const
This method returns a Displacement of an Angle relative to the current SurfacePoint latitude.
SurfacePoint()
Constructs an empty SurfacePoint object.
double MetersToLatitude(double latLength)
This method returns an angular measure of a distance in the direction of and relative to the latitude...
double DisplacementToDouble(Displacement disp, CoordUnits units)
This method returns a double version of a Displacement in the specified units.
double GetYWeight() const
Return Y weight for bundle adjustment Units are 1/(kilometers)^2.
double GetSigma(CoordinateType type, CoordIndex index, CoordUnits units)
This method returns a sigma of a SurfacePoint coordinate.
double GetLatWeight() const
Return latitude weight for bundle adjustment Units are 1/(radians)^2.
Distance GetSigmaDistance(CoordinateType type, CoordIndex index)
This method returns a sigma of a SurfacePoint coordinate as a Distance.
std::vector< double > Partial(CoordinateType type, CoordIndex index)
Compute partial derivative of conversion from body-fixed coordinates to the specified.
void SetRectangularMatrix(const boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > &covar, SurfacePoint::CoordUnits units=SurfacePoint::Meters)
Set rectangular covariance matrix and store in units of km**2.
void ResetLocalRadius(const Distance &radius)
This method resets the local radius of a SurfacePoint.
void FromNaifArray(const double naifValues[3])
A naif array is a c-style array of size 3.
double GetXWeight() const
Return X weight for bundle adjustment Units are 1/(kilometers)^2.
Distance GetLatSigmaDistance() const
Return the latitude sigma as a Distance.
void SetSphericalPoint(const Latitude &lat, const Longitude &lon, const Distance &radius)
This is a private method to set a surface point in spherical (lat/lon/radius), body-fixed coordinates...
Distance GetLocalRadius() const
Return the radius of the surface point.
double LonToDouble(Longitude lon, CoordUnits units)
This method returns a double version of a Longitude in the specified units.
void InitPoint()
Initialize a surface point.
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition IString.cpp:211
const double Null
Value for an Isis Null pixel.
bool IsSpecial(const double d)
Returns if the input pixel is special.
const double PI
The mathematical constant PI.
Definition Constants.h:40
Namespace for the standard library.