|
| 1 | +/* |
| 2 | + * Example demonstrating how to parse SBP data from a TCP port and compute the |
| 3 | + * ellipse representing the 2 sigma confidence level. |
| 4 | + * |
| 5 | + * Requires libsbp (https://github.com/swift-nav/libsbp) and Eigen 3 (https://eigen.tuxfamily.org/). |
| 6 | + * Eigen 3 is provided by the libeigen3-dev package on Debian-like Linux distributions. |
| 7 | + * |
| 8 | + * Compile using: |
| 9 | + * g++ tcp_2sigma_example.cc -I ${LIBSBP}/c/include/ -L ${LIBSBP}/c/build/src/ -lsbp |
| 10 | + * |
| 11 | + * If compiled with -DSHOW_PLOT then the ellipses will be graphically displayed. |
| 12 | + * This functionality requires OpenCV3 (https://opencv.org/) which is provided |
| 13 | + * by the libopencv-dev package on Debian-like Linux distributions. |
| 14 | + * |
| 15 | + * Compile using: |
| 16 | + * g++ tcp_2sigma_example.cc -DSHOW_PLOT -I ${LIBSBP}/c/include/ -L ${LIBSBP}/c/build/src/ -lsbp -lopencv_core -lopencv_imgproc -lopencv_highgui |
| 17 | + * |
| 18 | + * Originally adapted from https://github.com/swift-nav/libsbp/blob/master/c/example/tcp_example/tcp_example.c |
| 19 | + */ |
| 20 | + |
| 21 | +#ifndef SHOW_PLOT |
| 22 | +#define SHOW_PLOT 0 |
| 23 | +#endif /* SHOW_PLOT */ |
| 24 | + |
| 25 | +#include <stdio.h> |
| 26 | +#include <sys/socket.h> |
| 27 | +#include <arpa/inet.h> |
| 28 | +#include <stdlib.h> |
| 29 | +#include <string.h> |
| 30 | +#include <unistd.h> |
| 31 | +#include <cmath> |
| 32 | +#include <vector> |
| 33 | +#include <eigen3/Eigen/Core> |
| 34 | +#include <eigen3/Eigen/Eigenvalues> |
| 35 | +#if SHOW_PLOT |
| 36 | +#include <opencv2/core.hpp> |
| 37 | +#include <opencv2/core/eigen.hpp> |
| 38 | +#include <opencv2/imgproc.hpp> |
| 39 | +#include <opencv2/highgui.hpp> |
| 40 | +#endif /* SHOW_PLOT */ |
| 41 | + |
| 42 | +#include <libsbp/sbp.h> |
| 43 | + |
| 44 | +// width and height of the plot window in pixels |
| 45 | +#define SIZE_PIXELS 300 |
| 46 | +// width and height of the plot window in metres |
| 47 | +#define SIZE_METRES 10 |
| 48 | + |
| 49 | + |
| 50 | +void usage(const char* prog_name) { |
| 51 | + fprintf(stderr, "usage: %s [-a address -p port]\n", prog_name); |
| 52 | +} |
| 53 | + |
| 54 | + |
| 55 | +int open_socket(const char* host, int port) { |
| 56 | + struct sockaddr_in server; |
| 57 | + int fd = socket(AF_INET , SOCK_STREAM , 0); |
| 58 | + |
| 59 | + if(fd < 0) { |
| 60 | + fprintf(stderr, "Could not create socket\n"); |
| 61 | + return -1; |
| 62 | + } |
| 63 | + |
| 64 | + memset(&server, '0', sizeof(server)); |
| 65 | + server.sin_addr.s_addr = inet_addr(host); |
| 66 | + server.sin_family = AF_INET; |
| 67 | + server.sin_port = htons(port); |
| 68 | + |
| 69 | + if(connect(fd, (struct sockaddr *)&server , sizeof(server)) < 0) { |
| 70 | + fprintf(stderr, "Connection error\n"); |
| 71 | + close(fd); |
| 72 | + return -1; |
| 73 | + } |
| 74 | + |
| 75 | + return fd; |
| 76 | +} |
| 77 | + |
| 78 | + |
| 79 | +/* |
| 80 | + * Calculate the half major axis length, half minor axis length and rotation |
| 81 | + * angle for an error ellipse corresponding to the confidence interval |
| 82 | + * 'chi_squared' with the covariance matrix 'covmat'. |
| 83 | + * |
| 84 | + * Assumes normal distribution. |
| 85 | + * |
| 86 | + * See https://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/ |
| 87 | + * for an explanation of the algorithm. |
| 88 | + */ |
| 89 | +bool get_ellipse_parameters(const double& chi_squared, const Eigen::Matrix2d& covmat, double& half_major_axis, double& half_minor_axis, double& angle) { |
| 90 | + std::vector<std::tuple<double, Eigen::Vector2d>> eigen_vectors_and_values; |
| 91 | + |
| 92 | + // compute eigenvalues and eigenvectors |
| 93 | + Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eigensolver(covmat); |
| 94 | + if(eigensolver.info() != Eigen::Success) { |
| 95 | + return false; |
| 96 | + } |
| 97 | + |
| 98 | + // copy results into eigen_vectors_and_values |
| 99 | + Eigen::Vector2d eigen_values = eigensolver.eigenvalues(); |
| 100 | + Eigen::Matrix2d eigen_vectors = eigensolver.eigenvectors(); |
| 101 | + for(int i = 0; i < eigen_values.size(); i++) { |
| 102 | + std::tuple<double, Eigen::Vector2d> vec_and_val(eigen_values[i], eigen_vectors.col(i)); |
| 103 | + eigen_vectors_and_values.push_back(vec_and_val); |
| 104 | + } |
| 105 | + |
| 106 | + // sort eigen_vectors_and_values in descending order of eigenvalue |
| 107 | + std::sort(eigen_vectors_and_values.begin(), eigen_vectors_and_values.end(), |
| 108 | + [&](const std::tuple<double, Eigen::Vector2d>& a, const std::tuple<double, Eigen::Vector2d>& b) -> bool{ |
| 109 | + return std::get<double>(b) < std::get<double>(a); |
| 110 | + }); |
| 111 | + |
| 112 | + // calculate the angle between the largest eigenvector and the x-axis |
| 113 | + Eigen::Vector2d largest_eigenvector = std::get<Eigen::Vector2d>(eigen_vectors_and_values[0]); |
| 114 | + angle = atan2(largest_eigenvector[1], largest_eigenvector[0]); |
| 115 | + |
| 116 | + // shift the angle to the [0, 2pi] interval instead of [-pi, pi] |
| 117 | + if(angle < 0) angle += 2 * M_PI; |
| 118 | + |
| 119 | + // convert to degrees |
| 120 | + angle = 180. * angle / M_PI; |
| 121 | + |
| 122 | + // calculate the length of the minor and major axes |
| 123 | + half_major_axis = chi_squared * sqrt(std::get<double>(eigen_vectors_and_values[0])); |
| 124 | + half_minor_axis = chi_squared * sqrt(std::get<double>(eigen_vectors_and_values[1])); |
| 125 | + |
| 126 | + return true; |
| 127 | +} |
| 128 | + |
| 129 | + |
| 130 | +#if SHOW_PLOT |
| 131 | +void plot_ellipse(const double& half_major_axis, const double& half_minor_axis, const double& angle) { |
| 132 | + // the mean is nominally at the origin (i.e. reported lat/lon) but |
| 133 | + // for graphing purposes we shift it to the middle of our plot window |
| 134 | + cv::Point2f mean(SIZE_PIXELS/2, SIZE_PIXELS/2); |
| 135 | + |
| 136 | + // scale factor from metres to pixels |
| 137 | + double scale = SIZE_PIXELS / SIZE_METRES; |
| 138 | + |
| 139 | + // negative angle because OpenCV defines the angle clockwise instead of |
| 140 | + // anti-clockwise |
| 141 | + cv::RotatedRect ellipse(mean, cv::Size2f(half_major_axis * scale, half_minor_axis * scale), -angle); |
| 142 | + |
| 143 | + //Show the result |
| 144 | + cv::Mat visualizeimage(SIZE_PIXELS, SIZE_PIXELS, CV_8UC1, cv::Scalar::all(0)); |
| 145 | + cv::ellipse(visualizeimage, ellipse, cv::Scalar::all(255), 2); |
| 146 | + cv::imshow("2 Sigma Ellipse", visualizeimage); |
| 147 | +} |
| 148 | +#endif /* SHOW_PLOT */ |
| 149 | + |
| 150 | + |
| 151 | +void pos_llh_cov_callback(u16 /*sender_id*/, sbp_msg_type_t /*msg_type*/, const sbp_msg_t *msg, void* /*context*/) { |
| 152 | + const sbp_msg_pos_llh_cov_t* msg_pos_llh_cov = &msg->pos_llh_cov; |
| 153 | + if((msg_pos_llh_cov->flags & 0x7) != 0) { |
| 154 | + printf("TOW %u ms, Lat: %.15f deg, Lon: %.15f deg, Hgt: %.15f m\n", msg_pos_llh_cov->tow, msg_pos_llh_cov->lat, msg_pos_llh_cov->lon, msg_pos_llh_cov->height); |
| 155 | + |
| 156 | + // north/east submatrix |
| 157 | + Eigen::Matrix2d covmat; |
| 158 | + covmat << msg_pos_llh_cov->cov_n_n, msg_pos_llh_cov->cov_n_e, msg_pos_llh_cov->cov_n_e, msg_pos_llh_cov->cov_e_e; |
| 159 | + |
| 160 | + /* |
| 161 | + * Starling computes horizontal accuracy (as reported in |
| 162 | + * MSG_POS_LLH.h_accuracy) as the square root of the L2 norm of |
| 163 | + * the north/east submatrix. So the naive (but inaccurate) approach |
| 164 | + * is just to double it to get the worst-case 2 \Sigma value |
| 165 | + */ |
| 166 | + double hrms = sqrt(covmat.operatorNorm()) * 2.; |
| 167 | + printf("\t2DRMS (naive): %.3f m -> area %.3f m^2\n", hrms, M_PI * hrms * hrms); |
| 168 | + |
| 169 | + /* |
| 170 | + * A more accurate approach is to compute the error ellipse corresponding to |
| 171 | + * the 95% confidence interval |
| 172 | + */ |
| 173 | + double half_major_axis, half_minor_axis, angle; |
| 174 | + |
| 175 | + /* |
| 176 | + * From https://people.richland.edu/james/lecture/m170/tbl-chi.html |
| 177 | + * 95% confidence interval => area to the right = 0.05, 2 degrees of freedom => |
| 178 | + * sqrt(5.991) = 2.4477 |
| 179 | + * |
| 180 | + * For 99% confidence interval this would be sqrt(9.210) = 3.0348 |
| 181 | + */ |
| 182 | + const double chi_squared = 2.4477; |
| 183 | + |
| 184 | + if(get_ellipse_parameters(chi_squared, covmat, half_major_axis, half_minor_axis, angle)) { |
| 185 | + double area = M_PI * half_major_axis * half_minor_axis; |
| 186 | + printf("\tEllipse parameters: half major axis %g m, half minor axis %g m, angle %g deg -> area %.3f m^2\n", half_major_axis, half_minor_axis, angle, area); |
| 187 | + |
| 188 | +#if SHOW_PLOT |
| 189 | + plot_ellipse(half_major_axis, half_minor_axis, angle); |
| 190 | +#endif /* SHOW_PLOT */ |
| 191 | + } |
| 192 | + } else { |
| 193 | + printf("No fix\n"); |
| 194 | + } |
| 195 | +} |
| 196 | + |
| 197 | + |
| 198 | +s32 socket_read(u8* buff, u32 n, void* context) { |
| 199 | + int fd = *(int*)context; |
| 200 | + return read(fd, buff, n); |
| 201 | +} |
| 202 | + |
| 203 | + |
| 204 | +int main(int argc, char* const argv[]) { |
| 205 | + int opt; |
| 206 | + char* host = NULL; |
| 207 | + char* port = NULL; |
| 208 | + sbp_state_t sbp_state; |
| 209 | + sbp_msg_callbacks_node_t callback_node[1]; // note: one node is required per callback |
| 210 | + int fd; |
| 211 | + |
| 212 | + if(argc <= 2) { |
| 213 | + usage(argv[0]); |
| 214 | + exit(EXIT_FAILURE); |
| 215 | + } |
| 216 | + |
| 217 | + while((opt = getopt(argc, argv, "a:p:")) != -1) { |
| 218 | + switch (opt) { |
| 219 | + case 'a': |
| 220 | + host = optarg; |
| 221 | + break; |
| 222 | + case 'p': |
| 223 | + port = optarg; |
| 224 | + break; |
| 225 | + case 'h': |
| 226 | + usage(argv[0]); |
| 227 | + exit(EXIT_FAILURE); |
| 228 | + default: |
| 229 | + break; |
| 230 | + } |
| 231 | + } |
| 232 | + |
| 233 | + if(!host || !port) { |
| 234 | + fprintf(stderr, "Please supply the address and port of the SBP data stream!\n"); |
| 235 | + exit(EXIT_FAILURE); |
| 236 | + } |
| 237 | + |
| 238 | + fd = open_socket(host, atoi(port)); |
| 239 | + if(fd < 0) { |
| 240 | + return EXIT_FAILURE; |
| 241 | + } |
| 242 | + |
| 243 | +#if SHOW_PLOT |
| 244 | + cv::namedWindow("2 Sigma Ellipse"); |
| 245 | + cv::startWindowThread(); |
| 246 | +#endif /* SHOW_PLOT */ |
| 247 | + |
| 248 | + sbp_state_init(&sbp_state); |
| 249 | + sbp_state_set_io_context(&sbp_state, (void*)&fd); |
| 250 | + sbp_callback_register(&sbp_state, SbpMsgPosLlhCov, &pos_llh_cov_callback, NULL, &callback_node[0]); |
| 251 | + |
| 252 | + while(1) { |
| 253 | + switch(sbp_process(&sbp_state, &socket_read)) { |
| 254 | + case SBP_OK: |
| 255 | + case SBP_OK_CALLBACK_EXECUTED: |
| 256 | + case SBP_OK_CALLBACK_UNDEFINED: |
| 257 | + break; |
| 258 | + case SBP_CRC_ERROR: |
| 259 | + printf("CRC Error\n"); |
| 260 | + break; |
| 261 | + } |
| 262 | + } |
| 263 | + |
| 264 | + close(fd); |
| 265 | + return EXIT_SUCCESS; |
| 266 | +} |
0 commit comments