1 // This is bbas/vsph/vsph_spherical_coord.cxx
2 #include <complex>
3 #include <iostream>
4 #include <cmath>
5 #include "vsph_spherical_coord.h"
6 #ifdef _MSC_VER
7 # include "vcl_msvc_warnings.h"
8 #endif
9
10 #include <vgl/io/vgl_io_point_3d.h>
11
12 #define RADIUS_THRESH 0.0001
13
14
spherical_coord(vgl_point_3d<double> cp,vsph_sph_point_3d & sp)15 void vsph_spherical_coord::spherical_coord(vgl_point_3d<double> cp, vsph_sph_point_3d& sp)
16 {
17 // translate the point to the spherical coordinate system
18 double x = cp.x() - origin_.x();
19 double y = cp.y() - origin_.y();
20 double z = cp.z() - origin_.z();
21
22 double radius = std::sqrt(x*x+y*y+z*z);
23
24 // if radius is zero, the rsult does not make sense
25 if (std::abs(radius) < RADIUS_THRESH) {
26 sp.set(0.0,0.0,0.0);
27 return;
28 }
29
30 double phi = std::atan2(y,x);
31 double theta = std::acos(z/radius);
32 sp.set(radius, theta, phi);
33 }
34
cart_coord(vsph_sph_point_3d const & p) const35 vgl_point_3d<double> vsph_spherical_coord::cart_coord(vsph_sph_point_3d const& p) const
36 {
37 // Important note: does not use the "radius_" part of *this! Only origin_
38 double x = p.radius_*std::sin(p.theta_)*std::cos(p.phi_);
39 double y = p.radius_*std::sin(p.theta_)*std::sin(p.phi_);
40 double z = p.radius_*std::cos(p.theta_);
41 // translate the point based on the origin
42 return {x+origin_.x(),y+origin_.y(),z+origin_.z()};
43 }
44
move_point(vsph_sph_point_3d & p)45 bool vsph_spherical_coord::move_point(vsph_sph_point_3d& p)
46 {
47 // if it is already on the sphere
48 if (std::abs(p.radius_ - radius_) < 0.01)
49 return false;
50
51 // create a new point with the right radius
52 p.set(radius_, p.theta_, p.phi_);
53 return true;
54 }
radial_vector(vsph_sph_point_3d const & p)55 vgl_vector_3d<double> vsph_spherical_coord::radial_vector(vsph_sph_point_3d const& p){
56 vgl_point_3d<double> p_cart = cart_coord(p);
57 vgl_vector_3d<double> radial_dir = p_cart-origin_;
58 radial_dir /= radial_dir.length();
59 return radial_dir;
60 }
61
print(std::ostream & os) const62 void vsph_spherical_coord::print(std::ostream& os) const
63 {
64 os << " vsph_spherical_coord:[radius=" <<radius_ << ", origin=" << origin_ << "] ";
65 }
66
b_read(vsl_b_istream & is)67 void vsph_spherical_coord::b_read(vsl_b_istream& is)
68 {
69 short version;
70 vsl_b_read(is, version);
71 switch (version) {
72 case 1:
73 vsl_b_read(is, radius_);
74 vsl_b_read(is, origin_);
75 break;
76 default:
77 std::cerr << "I/O ERROR: vsl_b_read(vsl_b_istream&, vsph_spherical_coord&)\n"
78 << " Unknown version number "<< version << '\n';
79 is.is().clear(std::ios::badbit); // Set an unrecoverable IO error on stream
80 break;
81 }
82 }
83
b_write(vsl_b_ostream & os)84 void vsph_spherical_coord::b_write(vsl_b_ostream& os)
85 {
86 vsl_b_write(os, version());
87 vsl_b_write(os, radius_);
88 vsl_b_write(os, origin_);
89 }
90
operator <<(std::ostream & os,vsph_spherical_coord const & p)91 std::ostream& operator<<(std::ostream& os, vsph_spherical_coord const& p)
92 {
93 p.print(os);
94 return os;
95 }
96
operator <<(std::ostream & os,vsph_spherical_coord_sptr const & p)97 std::ostream& operator<<(std::ostream& os, vsph_spherical_coord_sptr const& p)
98 {
99 p->print(os);
100 return os;
101 }
102