Merge branch 'master' into pyramid-smart-pyramid-selection

.gitignore

@@ -8,6 +8,13 @@
 *.txt
 *.bin
 
+*.out.*
+*.lisp
+*.org
+*.py
+*.json
+/.gdb_history
+
 BSC5
 
 /documentation/man-*.h

Makefile

@@ -57,9 +57,6 @@ release: CXXFLAGS := $(RELEASE_CXXFLAGS)
 release: LDFLAGS := $(RELEASE_LDFLAGS)
 release: all
 
-$(BSC): download-bsc.sh
-	./download-bsc.sh
-
 $(BIN): $(OBJS)
 	$(CXX) $(LDFLAGS) -o $(BIN) $(OBJS) $(LIBS)
 

documentation/pipeline.man

@@ -92,6 +92,13 @@ The field of view of the camera that took the picture (in degrees). Defaults to
 \fB--centroid-mag-filter\fP \fImin-mag\fP
 Will not consider centroids with magnitude below \fImin-mag\fP.
 
+.TP
+\fB--centroid-filter-brightest\fP \fInum-stars\fP
+Remove all but the brightest \fInum-stars\fP many stars from the list of centroids before sending to
+star-id. Often a better choice than \fB--centroid-mag-filter\fP, because you can ensure that you
+keep enough stars to do star-id. If both this option and \fB--centroid-mag-filter\fP are provided,
+then all stars satisfying both criteria are kept (intersection).
+
 .TP
 \fB--database\fP \fIfilename\fP
 Chooses \fIfilename\fP as the database to use during star identification.
@@ -149,8 +156,14 @@ Sets the horizontal resolution of the generated image(s) to \fIpixels\fP. Defaul
 Sets the vertical resolution of the generated image(s) to \fIpixels\fP. Defaults to 1024.
 
 .TP
-\fB--generate-reference-brightness\fP \fIobserved-brightness\fP
-A star with magnitude 0 in the generated image will have have brightness \fIobserved-brightness\fP. An observed brightness of 1 fully saturates the pixel at the center of the star (completely white). Defaults to 100.
+\fB--generate-zero-mag-photons\fP \fInum-photoelectrons\fP
+A star with magnitude 0 will cause \fInum-photoelectrons\fP many photoelectrons to be received by
+the sensor. A default value of 20,000 is chosen. See Liebe's, "tutorial on star tracker accuracy" for theoretical information on how to calculate this.
+
+.TP
+\fB--generate-saturation-photons\fP \fIsaturation-photoelectrons\fP
+
+When a pixel receives at least this many photoelectrons, it will appear completely white (at least before other noise is applied). Note that, because of noise, a pixel may still appear completely white if it receives less than \fIsaturation-photoelectrons\fP many photoelectrons.
 
 .TP
 \fB--generate-spread-stddev\fP \fIstddev\fP
@@ -164,10 +177,6 @@ Enables or disables shot noise simulation in generated images. Defaults to true.
 \fB--generate-dark-current\fP \fInoise-level\fP
 Set observed brightness of dark current in the image, from 0 (no dark noise) to 1 (whole image pure white). Defaults to 0.1. cf \fB--generate-sensitivity\fP to control shot noise intensity.
 
-.TP
-\fB--generate-sensitivity\fP \fIsensitivity\fP
-Controls the simulated camera sensitivity for generated images. This only has an observable effect when shot noise is enabled (cf \fB--generate-shot-noise\fP). A higher sensitivity means that fewer photons can cause the same observed brightness on the sensor, so variation in number of photons causes more shot noise. Brightness is normalized according to \fB--generate-reference-brightness\fP so that adjusting sensitivity does not change observed brightness.
-
 .TP
 \fB--generate-read-noise-stddev\fP \fIstddev\fP
 Sets the standard deviation of Gaussian noise in the generated image(s) to \fIstddev\fP. Noise is measured in observed brightness, where 1 is the difference between pure white and pure black. Defaults to 0.05.
@@ -267,16 +276,28 @@ perturbations and magnitude cutoff is applied).
 
 .TP
 \fB--print-actual-centroids\fP [\fIpath\fP]
-Mainly for debugging. Argument is optional. Print list of centroids in the pipeline output to
-\fIpath\fP. If star identification was also performed, print which catalog star each centroid was
-identified as. Defaults to stdout.
+Argument is optional. Print list of centroids in the pipeline output to \fIpath\fP. If star
+identification was also performed, print which catalog star each centroid was identified as. If
+input is from a generated image, also print the expected id for each centroid (note that there can
+sometimes be multiple acceptable star-ids for a star, if it's within the centroid comparison
+threshold of multiple catalog stars. Usually setting --min-separation sufficiently high prevents
+this). Defaults to stdout.
+
+.TP
+\fB--print-input-centroids\fP [\fIpath\fP]
+Argument is optional. Print list of centroids in the pipeline input to \fIpath\fP. This only works
+for generated images. Prints whatever centroids would normally be sent to the star-id algorithm if
+provided in the pipeline. Also prints star-ids, but don't rely on them, because if the input
+centroids were perturbed the generated star-id may not actually be the closest catalog star anymore!
+Defaults to stdout.
 
 .TP
 \fB--print-expected-centroids\fP [\fIpath\fP]
-Mainly for debugging. Argument is optional. Print list of "expected" centroids from the pipeline
-input to \fIpath\fP. Also prints out expected star identifications for those stars, if present. This
-only works for generated images. Generally you should use \fB--print-actual-centroids\fP instead,
-and use \fB--compare-centroids\fP if you need to compare actual and expected centroids.
+Argument is optional. Print list of "expected" centroids from the pipeline input to \fIpath\fP. The
+"expected" centroids are all centroids used to generate the real image, including false stars and
+stars below the cutoff magnitude. Also prints out expected star identifications for those stars.
+This only works for generated images. Generally you should use \fB--print-actual-centroids\fP
+instead, and use \fB--compare-centroids\fP if you need to compare actual and expected centroids.
 
 .TP
 \fB--print-attitude\fP [\fIpath\fP]

src/attitude-estimators.cpp

@@ -10,6 +10,10 @@ Attitude DavenportQAlgorithm::Go(const Camera &camera,
                                  const Stars &stars,
                                  const Catalog &catalog,
                                  const StarIdentifiers &starIdentifiers) {
+    if (starIdentifiers.size() < 2) {
+        return Attitude();
+    }
+    assert(stars.size() >= 2);
 
     // attitude profile matrix
     Eigen::Matrix3f B;
@@ -88,9 +92,10 @@ Attitude TriadAlgorithm::Go(const Camera &camera,
                             const Stars &stars,
                             const Catalog &catalog,
                             const StarIdentifiers &starIds) {
-    if ((int)stars.size() < 2 || (int)starIds.size() < 2) {
-        return Attitude(Quaternion(1,0,0,0));
+    if (starIds.size() < 2) {
+        return Attitude();
     }
+    assert(stars.size() >= 2);
 
     // TODO: Better way of picking the two stars
     StarIdentifier
@@ -138,6 +143,11 @@ Attitude QuestAlgorithm::Go(const Camera &camera,
                             const Catalog &catalog,
                             const StarIdentifiers &starIdentifiers) {
 
+    if (starIdentifiers.size() < 2) {
+        return Attitude();
+    }
+    assert(stars.size() >= 2);
+
     // initial guess for eigenvalue (sum of the weights)
     float guess = 0;
 

src/attitude-utils.cpp

@@ -62,7 +62,13 @@ float Quaternion::Angle() const {
     return (real >= 1 ? 0 : acos(real))*2;
 }
 
-/// Change the amount of rotation in a quaternion, keeping that rotation around the same axis.
+float Quaternion::SmallestAngle() const {
+    float rawAngle = Angle();
+    return rawAngle > M_PI
+        ? 2*M_PI - rawAngle
+        : rawAngle;
+}
+
 void Quaternion::SetAngle(float newAngle) {
     real = cos(newAngle/2);
     SetVector(Vector().Normalize() * sin(newAngle/2));
@@ -336,10 +342,10 @@ Mat3 Mat3::Inverse() const {
                      0,0,1};
 
 Attitude::Attitude(const Quaternion &quat)
-    : quaternion(quat), type(QuaternionType) {}
+    : type(AttitudeType::QuaternionType), quaternion(quat) {}
 
 Attitude::Attitude(const Mat3 &matrix)
-    : dcm(matrix), type(DCMType) {}
+    : type(AttitudeType::DCMType), dcm(matrix) {}
 
 /// Convert a quaternion to a rotation matrix (Direction Cosine Matrix)
 Mat3 QuaternionToDCM(const Quaternion &quat) {
@@ -385,48 +391,58 @@ Quaternion DCMToQuaternion(const Mat3 &dcm) {
 /// Get the quaternion representing the attitude, converting from whatever format is stored.
 Quaternion Attitude::GetQuaternion() const {
     switch (type) {
-        case QuaternionType:
-            return quaternion;
-        case DCMType:
-            return DCMToQuaternion(dcm);
-        default:
-            assert(false);
+    case AttitudeType::QuaternionType:
+        return quaternion;
+    case AttitudeType::DCMType:
+        return DCMToQuaternion(dcm);
+    default:
+        assert(false);
     }
 }
 
 /// Get the rotation matrix (direction cosine matrix) representing the attitude, converting from whatever format is stored.
 Mat3 Attitude::GetDCM() const {
     switch (type) {
-        case DCMType:
-            return dcm;
-        case QuaternionType:
-            return QuaternionToDCM(quaternion);
-        default:
-            assert(false);
+    case AttitudeType::DCMType:
+        return dcm;
+    case AttitudeType::QuaternionType:
+        return QuaternionToDCM(quaternion);
+    default:
+        assert(false);
     }
 }
 
 /// Convert a vector from the reference frame to the body frame.
 Vec3 Attitude::Rotate(const Vec3 &vec) const {
     switch (type) {
-        case DCMType:
-            return dcm*vec;
-        case QuaternionType:
-            return quaternion.Rotate(vec);
-        default:
-            assert(false);
+    case AttitudeType::DCMType:
+        return dcm*vec;
+    case AttitudeType::QuaternionType:
+        return quaternion.Rotate(vec);
+    default:
+        assert(false);
     }
 }
 
 /// Get the euler angles from the attitude, converting from whatever format is stored.
 EulerAngles Attitude::ToSpherical() const {
     switch (type) {
-        case DCMType:
-            return GetQuaternion().ToSpherical();
-        case QuaternionType:
-            return quaternion.ToSpherical();
-        default:
-            assert(false);
+    case AttitudeType::DCMType:
+        return GetQuaternion().ToSpherical();
+    case AttitudeType::QuaternionType:
+        return quaternion.ToSpherical();
+    default:
+        assert(false);
+    }
+}
+
+bool Attitude::IsKnown() const {
+    switch (type) {
+    case AttitudeType::DCMType:
+    case AttitudeType::QuaternionType:
+        return true;
+    default:
+        return false;
     }
 }
 

src/attitude-utils.hpp

@@ -107,6 +107,9 @@ class Quaternion {
     void SetVector(const Vec3 &);
     Vec3 Rotate(const Vec3 &) const;
     float Angle() const;
+    /// Returns the smallest angle that can be used to represent the rotation represented by the
+    /// quaternion. I.e, min(Angle, 2pi-Angle).
+    float SmallestAngle() const;
     void SetAngle(float);
     EulerAngles ToSpherical() const;
     bool IsUnit(float tolerance) const;
@@ -121,13 +124,17 @@ 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
+ * The Attitude object stores either a rotation matrix (direction cosine matrix) or a quaternion,
+ * and converts automatically to the other format when needed. Importantly, an attitude may also be
+ * "unknown", representing for example if there weren't enough stars to determine an attitude. When
+ * set to unknown, it is illegal to try and convert it to anything, so check IsKnown first!
+ * @note When porting to an embedded device, you may 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:
+    /// constructs unknown attitude:
     Attitude() = default;
     explicit Attitude(const Quaternion &); // NOLINT
     explicit Attitude(const Mat3 &dcm);
@@ -136,33 +143,34 @@ class Attitude {
     Mat3 GetDCM() const;
     EulerAngles ToSpherical() const;
     Vec3 Rotate(const Vec3 &) const;
+    bool IsKnown() const;
 
 private:
-    enum AttitudeType {
-        NullType,
+    enum class AttitudeType {
+        UnknownType, /// Doesn't mean "unknown type", but rather "we have no fucking clue where we're pointing"
         QuaternionType,
         DCMType,
     };
 
+    AttitudeType type;
     Quaternion quaternion;
     Mat3 dcm; // direction cosine matrix
-    AttitudeType type;
 };
 
 Mat3 QuaternionToDCM(const Quaternion &);
 Quaternion DCMToQuaternion(const Mat3 &);
 
-// Return a quaternion that will reorient the coordinate axes so that the x-axis points at the given
-// right ascension and declination, then roll the coordinate axes counterclockwise (i.e., the stars
-// will appear to rotate clockwise). This is an "improper" z-y'-x' Euler rotation.
+/// Return a quaternion that will reorient the coordinate axes so that the x-axis points at the given
+/// right ascension and declination, then roll the coordinate axes counterclockwise (i.e., the stars
+/// will appear to rotate clockwise). This is an "improper" z-y'-x' Euler rotation.
 Quaternion SphericalToQuaternion(float ra, float dec, float roll);
 
-// returns unit vector
+/// returns unit vector
 Vec3 SphericalToSpatial(float ra, float de);
 void SpatialToSpherical(const Vec3 &, float *ra, float *de);
-// angle between two vectors, using dot product and magnitude division
+/// angle between two vectors, using dot product and magnitude division
 float Angle(const Vec3 &, const Vec3 &);
-// angle between two vectors, /assuming/ that they are already unit length
+/// angle between two vectors, /assuming/ that they are already unit length
 float AngleUnit(const Vec3 &, const Vec3 &);
 
 float RadToDeg(float);