Tetra Star-ID algorithm

.gitignore

@@ -11,7 +11,7 @@
 *.out.*
 *.lisp
 *.org
-*.py
+# *.py
 *.json
 /.gdb_history
 

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
@@ -85,6 +86,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)

documentation/database.man

@@ -80,6 +80,19 @@ Sets the maximum distance \fImax\fP (degrees). Defaults to 15 if option is not s
 \fB--kvector-distance-bins\fP \fInum-bins\fP
 Sets the number of distance bins in the kvector building method to \fInum-bins\fP.  Defaults to 10000 if option is not selected, which is pretty reasonable for most cases.
 
+.SH TETRA DATABASE OPTIONS
+
+The Tetra database is needed for the Tetra star identification algorithm. Its primary payload is
+a pattern catalog , essentially a large hash table that stores 4-star patterns
+
+.TP
+\fB--tetra\fP
+Generate a Tetra database
+
+.TP
+\fB--tetra-max-angle \fImax\fP
+Sets the maximum allowable angle between stars in a Tetra pattern to \fImax\fP (degrees). Defaults to 12 if option is not selected.
+
 .SH OTHER OPTIONS
 
 .TP

documentation/pipeline.man

@@ -121,7 +121,7 @@ Runs the dummy centroiding algorithm (random centroid algorithm) with \fInum-sta
 
 .TP
 \fB--star-id-algo\fP \fIalgo\fP
-Runs the \fIalgo\fP star identification algorithm. Current options are "dummy", "gv", and "pyramid". Defaults to "dummy" if option is not selected.
+Runs the \fIalgo\fP star identification algorithm. Current options are "dummy", "gv", "tetra", and "pyramid". Defaults to "dummy" if option is not selected.
 
 .TP
 \fB--angular-tolerance\fP [\fItolerance\fP] Sets the estimated angular centroiding error tolerance,

reg.py

@@ -0,0 +1,65 @@
+# Script to check performance on images in given folder and prevent regression testing
+# Go over all images in some folder and check performance, log it to stdout (for now)
+# TODO: figure out some storage plan
+
+# CLI:
+# test directory name, output log name
+# (optional) attitude estimator
+# Example usage: python reg.py samples/ log.txt -att_estimator quest
+
+import subprocess
+import argparse
+
+# import os
+import glob
+import datetime
+
+parser = argparse.ArgumentParser()
+parser.add_argument("test_dir", type=str)
+parser.add_argument("log", type=str, default="log.txt")
+parser.add_argument("-att_estimator", type=str, nargs="?", default="dqm")
+args = parser.parse_args()
+
+print(f"Testing images in {args.test_dir}")
+print(f"attitude estimator: {args.att_estimator}")
+print(f"Logging outputs to {args.log}")
+
+
+def get_diff(expected, actual):
+    """Get element-wise angle difference between expected and actual (what we got)"""
+    return [expected[i] - actual[i] for i in range(len(actual))]
+
+
+output_log = open(args.log, "a+")  # append to end of log, don't overwrite
+for img_name in glob.glob(args.test_dir + "**/*.png", recursive=True):
+    print(f"===================={img_name}====================")
+    cmd = (
+        f"./lost pipeline \
+      --png {img_name} \
+      --fov 17 \
+      --centroid-algo cog \
+      --centroid-filter-brightest 8 \
+      --star-id-algo tetra \
+      --database tetra-algo-12.dat \
+      --false-stars 0 \
+      --attitude-algo {args.att_estimator} \
+      --print-attitude attitude.txt \
+      --plot-output annotated-res.png",
+    )
+    subprocess.call(cmd, shell=True)
+    # log attitude.txt
+    # Log:
+    # which image was tested
+    # Output from attitude.txt
+    dt = datetime.datetime.now().isoformat()
+    output_log.write("----------------------------\n")
+    output_log.write(f"{img_name}-{dt}-{args.att_estimator}\n")
+    output_log.write("----------------------------\n")
+    with open("attitude.txt") as att_log:
+        targets = ["attitude_ra", "attitude_de", "attitude_roll"]
+        for line in att_log:
+            line = line.split(" ")
+            if line[0] in targets:
+                output_log.write(line[1])
+
+output_log.close()

src/attitude-estimators.cpp

@@ -125,7 +125,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

@@ -169,6 +169,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));
@@ -200,12 +210,12 @@ 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 };
 }
@@ -225,6 +235,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 {

src/attitude-utils.hpp

@@ -1,7 +1,10 @@
 #ifndef ATTITUDE_UTILS_H
 #define ATTITUDE_UTILS_H
 
+#include <algorithm>
+#include <fstream>
 #include <memory>
+#include <numeric>  // iota
 #include <vector>
 
 #include "serialize-helpers.hpp"
@@ -26,6 +29,8 @@ 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
@@ -37,6 +42,9 @@ class Vec3 {
     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;
@@ -45,8 +53,11 @@ class Vec3 {
     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
@@ -183,8 +194,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.hpp

@@ -41,7 +41,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

@@ -11,4 +11,6 @@ LOST_CLI_OPTION("kvector-max-distance"   , float      , kvectorMaxDistance    ,
 LOST_CLI_OPTION("kvector-distance-bins"  , long       , kvectorNumDistanceBins, 10000 , atol(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("swap-integer-endianness", bool       , swapIntegerEndianness , false , atobool(optarg), true)
 LOST_CLI_OPTION("swap-float-endianness"  , bool       , swapFloatEndianness   , false , atobool(optarg), true)
+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)