Tetra Star-ID algorithm

Makefile

@@ -1,15 +1,16 @@
-# Copyright (c) 2020 Mark Polyakov, Karen Haining (If you edit the file, add your name here!)
-# 
+# Copyright (c) 2020 Mark Polyakov, Karen Haining, Edward Zhang
+# (If you edit the file, add your name here!)
+#
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
 # in the Software without restriction, including without limitation the rights
 # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 # copies of the Software, and to permit persons to whom the Software is
 # furnished to do so, subject to the following conditions:
-# 
+#
 # The above copyright notice and this permission notice shall be included in all
 # copies or substantial portions of the Software.
-# 
+#
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -88,6 +89,7 @@ test: $(BIN) $(BSC) $(TEST_BIN)
 	# bash ./test/scripts/pyramid-incorrect.sh
 	bash ./test/scripts/readme-examples-test.sh
 	bash ./test/scripts/random-crap.sh
+	bash ./test/scripts/tetra.sh
 
 $(TEST_BIN): $(TEST_OBJS)
 	$(CXX) $(LDFLAGS) -o $(TEST_BIN) $(TEST_OBJS) $(LIBS)

README.md

@@ -32,14 +32,14 @@ You need Linux or macOS. On Windows, we recommend installing the Windows Subsyst
   libeigen3-dev`. Elsewhere, download the latest stable release from https://eigen.tuxfamily.org/.
   You can install it system-wide, or just extract it so that the main folder is in
   `vendor/eigen3/Eigen` under the LOST repository.
-  
+
 ### macOS Prerequisites
 
 - Download [Homebrew](https://brew.sh/), run `/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"` in terminal
 - Install [cairo](https://formulae.brew.sh/formula/cairo#default) via homebrew `brew install cairo`
 - Install [Eigen]("https://formulae.brew.sh/formula/eigen#default") via homebrew `brew install eigen`
     - Locate Eigen files and move them to `vendor/eigen3/Eigen` under the LOST repository
-	
+
 ### Building
 
 Clone this repository (`git clone https://github.com/uwcubesat/lost`), then `cd lost`, then

download-bsc.sh

@@ -2,4 +2,4 @@
 
 # copied from Firefox dev tools. Alternately, follow documentation/downloading-stars.md
 curl 'https://vizier.cds.unistra.fr/viz-bin/asu-tsv' --compressed -H 'Content-Type: application/x-www-form-urlencoded' -H 'Origin: https://vizier.u-strasbg.fr' --data-raw '-ref=VIZ5f9b96c42e6d&-to=4&-from=-3&-this=-3&%2F%2Fsource=V%2F50%2Fcatalog&%2F%2Ftables=V%2F50%2Fcatalog&-out.max=unlimited&%2F%2FCDSportal=http%3A%2F%2Fcdsportal.u-strasbg.fr%2FStoreVizierData.html&-out.form=%7C+-Separated-Values&-out.add=_RAJ%2C_DEJ&%2F%2Foutaddvalue=default&-oc.form=dec&-nav=cat%3AV%2F50%26tab%3A%7BV%2F50%2Fcatalog%7D%26key%3Asource%3DV%2F50%2Fcatalog%26HTTPPRM%3A%26&-c=&-c.eq=J2000&-c.r=++2&-c.u=arcmin&-c.geom=r&-source=V%2F50%2Fcatalog&-order=I&recno=&-out=HR&HR=&Name=&DM=&HD=&SAO=&FK5=&IRflag=&r_IRflag=&-out=Multiple&Multiple=&ADS=&ADScomp=&VarID=&RAJ2000=&DEJ2000=&GLON=&GLAT=&-out=Vmag&Vmag=&n_Vmag=&u_Vmag=&B-V=&u_B-V=&U-B=&u_U-B=&R-I=&n_R-I=&SpType=&n_SpType=&pmRA=&pmDE=&n_Parallax=&Parallax=&RadVel=&n_RadVel=&l_RotVel=&RotVel=&u_RotVel=&Dmag=&Sep=&MultID=&MultCnt=&NoteFlag=&%2F%2Fnoneucd1p=on&-file=.&-meta.ucd=2&-meta=1&-meta.foot=1&-usenav=1&-bmark=POST' | \
-    grep -E '^[0-9]{3}\.[0-9]{6}' > bright-star-catalog.tsv
+    grep -E '^[0-9]{3}\.[0-9]{6}' > bright-star-catalog.tsv

src/attitude-estimators.cpp

@@ -120,7 +120,7 @@ float QuestCharPoly(float x, float a, float b, float c, float d, float s) {retur
 float QuestCharPolyPrime(float x, float a, float b, float c) {return 4*pow(x,3) - 2*(a+b)*x - c;}
 
 /**
- * Approximates roots of a real function using the Newton-Raphson algorithm 
+ * Approximates roots of a real function using the Newton-Raphson algorithm
  * @see https://www.geeksforgeeks.org/program-for-newton-raphson-method/
  */
 float QuestEigenvalueEstimator(float guess, float a, float b, float c, float d, float s) {

src/attitude-utils.cpp

@@ -161,6 +161,16 @@ float FloatModulo(float x, float mod) {
     return result >= 0 ? result : result + mod;
 }
 
+std::ostream& operator<< (std::ostream& output, const Vec2 &vec){
+    output << "Vec2(x: " << vec.x << ", y: " << vec.y << ")";
+    return output;
+}
+
+std::ostream& operator<< (std::ostream& output, const Vec3 &vec){
+    output << "Vec3(x: " << vec.x << ", y: " << vec.y << ", z: " << vec.z << ")";
+    return output;
+}
+
 /// The square of the magnitude
 float Vec3::MagnitudeSq() const {
     return fma(x,x,fma(y,y, z*z));
@@ -187,17 +197,24 @@ Vec3 Vec3::Normalize() const {
     };
 }
 
+// Allows Vec3 to be used as key in associative containers
+bool Vec3::operator<(const Vec3 &other) const{
+    if (x != other.x) return (x < other.x);
+    if (y != other.y) return (y < other.y);
+    return z < other.z;
+}
+
 /// Dot product
 float Vec3::operator*(const Vec3 &other) const {
     return fma(x,other.x, fma(y,other.y, z*other.z));
 }
 
-/// Dot product
+/// Vector-scalar multiplication
 Vec2 Vec2::operator*(const float &other) const {
     return { x*other, y*other };
 }
 
-/// Vector-Scalar multiplication
+/// Vector-scalar multiplication
 Vec3 Vec3::operator*(const float &other) const {
     return { x*other, y*other, z*other };
 }
@@ -217,6 +234,10 @@ Vec3 Vec3::operator-(const Vec3 &other) const {
     return { x - other.x, y - other.y, z - other.z };
 }
 
+Vec3 Vec3::operator+(const Vec3 &other) const{
+    return {x + other.x, y + other.y, z + other.z};
+}
+
 /// Usual vector cross product
 Vec3 Vec3::CrossProduct(const Vec3 &other) const {
     return {
@@ -473,4 +494,5 @@ float Distance(const Vec2 &v1, const Vec2 &v2) {
     return sqrt(pow(v1.x-v2.x, 2) + pow(v1.y-v2.y, 2));
 }
 
+
 }

src/attitude-utils.hpp

@@ -1,7 +1,11 @@
 #ifndef ATTITUDE_UTILS_H
 #define ATTITUDE_UTILS_H
 
+#include <algorithm>
+#include <fstream>
 #include <memory>
+#include <numeric>  // iota
+#include <vector>
 
 namespace lost {
 
@@ -23,32 +27,42 @@ struct Vec2 {
     Vec2 operator*(const float &) const;
     Vec2 operator-(const Vec2 &) const;
     Vec2 operator+(const Vec2 &) const;
+
+    friend std::ostream &operator<<(std::ostream &output, const Vec2 &vec);
 };
 
-class Mat3; // define above so we can use in Vec3 class
+class Mat3;  // define above so we can use in Vec3 class
 
 /// Three dimensional vector with floating point components
 class Vec3 {
-public:
+   public:
     float x;
     float y;
     float z;
 
+    Vec3(){};
+    Vec3(float x, float y, float z) : x(x), y(y), z(z){};
+
     float Magnitude() const;
     float MagnitudeSq() const;
     Vec3 Normalize() const;
 
+    bool operator<(const Vec3 &other) const;
+
     float operator*(const Vec3 &) const;
     Vec3 operator*(const float &) const;
     Vec3 operator*(const Mat3 &) const;
     Vec3 operator-(const Vec3 &) const;
+    Vec3 operator+(const Vec3 &) const;
     Vec3 CrossProduct(const Vec3 &) const;
     Mat3 OuterProduct(const Vec3 &) const;
+
+    friend std::ostream &operator<<(std::ostream &output, const Vec3 &vec);
 };
 
 /// 3x3 vector with floating point components
 class Mat3 {
-public:
+   public:
     float x[9];
 
     float At(int i, int j) const;
@@ -75,31 +89,31 @@ float Distance(const Vec3 &, const Vec3 &);
 
 /**
  * A "human-readable" way to represent a 3d rotation or orientation.
- * Euler angles roughly correspond to yaw, pitch, and roll of an airplane, which are easy for humans to understand.
- * There's no one single way to store Euler angles. We use z-y'-x'' angles, according to the notation used on the wikipedia page for euler angles.
+ * Euler angles roughly correspond to yaw, pitch, and roll of an airplane, which are easy for humans
+ * to understand. There's no one single way to store Euler angles. We use z-y'-x'' angles, according
+ * to the notation used on the wikipedia page for euler angles.
  */
 class EulerAngles {
-public:
-    EulerAngles(float ra, float de, float roll)
-        : ra(ra), de(de), roll(roll) { };
+   public:
+    EulerAngles(float ra, float de, float roll) : ra(ra), de(de), roll(roll){};
 
     /// Right ascension. How far we yaw left. Yaw is performed first.
     float ra;
-    /// Declination. How far we pitch up (or down if negative). Pitch is performed second, after yaw.
+    /// Declination. How far we pitch up (or down if negative). Pitch is performed second, after
+    /// yaw.
     float de;
     /// How far we roll counterclockwise. Roll is performed last (after yaw and pitch).
     float roll;
 };
 
 /// A quaternion is a common way to represent a 3d rotation.
 class Quaternion {
-public:
+   public:
     Quaternion() = default;
     explicit Quaternion(const Vec3 &);
     Quaternion(const Vec3 &, float);
 
-    Quaternion(float real, float i, float j, float k)
-        : real(real), i(i), j(j), k(k) { };
+    Quaternion(float real, float i, float j, float k) : real(real), i(i), j(j), k(k){};
 
     Quaternion operator*(const Quaternion &other) const;
     Quaternion Conjugate() const;
@@ -121,31 +135,32 @@ class Quaternion {
 //
 /**
  * The attitude (orientation) of a spacecraft.
- * The Attitude object stores either a rotation matrix (direction cosine matrix) or a quaternion, and converts automatically to the other format when needed.
- * @note When porting to an embedded device, you'll probably want to get rid of this class and adapt to
- * either quaternions or DCMs exclusively, depending on the natural output format of whatever
+ * The Attitude object stores either a rotation matrix (direction cosine matrix) or a quaternion,
+ * and converts automatically to the other format when needed.
+ * @note When porting to an embedded device, you'll probably want to get rid of this class and adapt
+ * to either quaternions or DCMs exclusively, depending on the natural output format of whatever
  * attitude estimation algorithm you're using.
  */
 class Attitude {
-public:
+   public:
     Attitude() = default;
-    explicit Attitude(const Quaternion &); // NOLINT
+    explicit Attitude(const Quaternion &);  // NOLINT
     explicit Attitude(const Mat3 &dcm);
 
     Quaternion GetQuaternion() const;
     Mat3 GetDCM() const;
     EulerAngles ToSpherical() const;
     Vec3 Rotate(const Vec3 &) const;
 
-private:
+   private:
     enum AttitudeType {
         NullType,
         QuaternionType,
         DCMType,
     };
 
     Quaternion quaternion;
-    Mat3 dcm; // direction cosine matrix
+    Mat3 dcm;  // direction cosine matrix
     AttitudeType type;
 };
 
@@ -173,8 +188,24 @@ float ArcSecToRad(float);
 /// Always returns something in [0,mod) Eg -0.8 mod 0.6 = 0.4
 float FloatModulo(float x, float mod);
 
+// Argsort function - Tetra
+// Sort first vector based on values of second vector (asc)
+template <class T, class U>
+std::vector<T> ArgsortVector(std::vector<T> arr, std::vector<U> cmp) {
+    std::vector<T> res;
+    std::vector<int> indices(arr.size());
+    std::iota(indices.begin(), indices.end(), 0);
+    std::sort(indices.begin(), indices.end(),
+              [&](int a, int b) -> bool { return cmp[a] < cmp[b]; });
+
+    for (int ind : indices) {
+        res.push_back(arr[ind]);
+    }
+    return res;
+}
+
 // TODO: quaternion and euler angle conversion, conversion between ascension/declination to rec9tu
 
-}
+}  // namespace lost
 
 #endif

src/camera.cpp

@@ -46,6 +46,23 @@ Vec3 Camera::CameraToSpatial(const Vec2 &vector) const {
     };
 }
 
+Vec3 Camera::CameraToSpatialFov(const Vec2 &vector) const{
+    // float fov = DegToRad(this->Fov());
+    float fov = this->Fov();
+    float centerX = this->XResolution() / 2.0;
+    float centerY = this->YResolution() / 2.0;
+    float scaleFactor = std::tan(fov / 2) / centerX;
+
+    float jOverI = (centerX - vector.x) * scaleFactor;
+    float kOverI = (centerY - vector.y) * scaleFactor;
+    float i = 1.0 / std::sqrt(1 + jOverI * jOverI + kOverI * kOverI);
+    float j = jOverI * i;
+    float k = kOverI * i;
+    Vec3 res(i, j, k);
+
+    return res;
+}
+
 /// Returns whether a given pixel is actually in the camera's field of view
 bool Camera::InSensor(const Vec2 &vector) const {
     // if vector.x == xResolution, then it is at the leftmost point of the pixel that's "hanging

src/camera.hpp

@@ -28,6 +28,8 @@ class Camera {
     Vec2 SpatialToCamera(const Vec3 &) const;
     Vec3 CameraToSpatial(const Vec2 &) const;
 
+    Vec3 CameraToSpatialFov(const Vec2 &) const;
+
     // converts from a 2d point in the camera sensor to right ascension and declination relative to
     // the center of the camera.
     // void CoordinateAngles(const Vec2 &vector, float *ra, float *de) const;
@@ -41,7 +43,7 @@ class Camera {
     /// Focal length in pixels
     float FocalLength() const { return focalLength; };
     /// Horizontal field of view in radians
-    float Fov() const;
+    float Fov() const; // in radians
 
     void SetFocalLength(float focalLength) { this->focalLength = focalLength; }
 

src/database-options.hpp

@@ -4,9 +4,11 @@
 
 LOST_CLI_OPTION("min-mag"              , float      , minMag                , 100   , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("max-stars"            , int        , maxStars              , 10000 , atoi(optarg)   , kNoDefaultArgument)
-LOST_CLI_OPTION("min-separation"       , float      , minSeparation         , 0.08  , atof(optarg)   , kNoDefaultArgument)
+LOST_CLI_OPTION("min-separation"       , float      , minSeparation         , 0.05  , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("kvector"              , bool       , kvector               , false , atobool(optarg), true)
 LOST_CLI_OPTION("kvector-min-distance" , float      , kvectorMinDistance    , 0.5   , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("kvector-max-distance" , float      , kvectorMaxDistance    , 15    , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("kvector-distance-bins", long       , kvectorNumDistanceBins, 10000 , atol(optarg)   , kNoDefaultArgument)
+LOST_CLI_OPTION("tetra"                , bool       , tetra                 , false , atobool(optarg), true)
+LOST_CLI_OPTION("tetra-max-angle"      , float      , tetraMaxAngle         , 12    , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("output"               , std::string, outputPath            , "-"   , optarg         , kNoDefaultArgument)