diff --git a/.github/workflows/build-test-lint.yml b/.github/workflows/build-test-lint.yml index 57f71511..93821aa2 100644 --- a/.github/workflows/build-test-lint.yml +++ b/.github/workflows/build-test-lint.yml @@ -12,10 +12,14 @@ jobs: steps: - uses: actions/checkout@v3 - - name: Build + - name: Build Double run: CXXFLAGS=-Werror make -j$(($(nproc)+1)) - - name: Test + - name: Test Double run: CXXFLAGS=-Werror make test -j$(($(nproc)+1)) + - name: Build Float + run: CXXFLAGS=-Werror make clean all LOST_FLOAT_MODE=1 -j$(($(nproc)+1)) + - name: Test Float + run: CXXFLAGS=-Werror make test LOST_FLOAT_MODE=1 -j$(($(nproc)+1)) lint: runs-on: ubuntu-latest diff --git a/Makefile b/Makefile index f20a4557..341dc9eb 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -# Copyright (c) 2020 Mark Polyakov, Karen Haining (If you edit the file, add your name here!) +# Copyright (c) 2020 Mark Polyakov, Karen Haining, Muki Kiboigo (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 @@ -51,6 +51,12 @@ ifndef LOST_DISABLE_ASAN LDFLAGS := $(LDFLAGS) -fsanitize=address endif +# Use Double Mode by default. +# If compiled with LOST_FLOAT_MODE=1, we will use floats. +ifdef LOST_FLOAT_MODE + CXXFLAGS := $(CXXFLAGS) -Wdouble-promotion -Werror=double-promotion -D LOST_FLOAT_MODE +endif + all: $(BIN) $(BSC) release: CXXFLAGS := $(RELEASE_CXXFLAGS) diff --git a/src/attitude-estimators.cpp b/src/attitude-estimators.cpp index 1e80f692..bf55d8f1 100644 --- a/src/attitude-estimators.cpp +++ b/src/attitude-estimators.cpp @@ -1,10 +1,13 @@ #include "attitude-estimators.hpp" + #include #include +#include "decimal.hpp" + namespace lost { -#define EPSILON 0.0001 // threshold for 0 for Newton-Raphson method +#define EPSILON DECIMAL(0.0001) // threshold for 0 for Newton-Raphson method Attitude DavenportQAlgorithm::Go(const Camera &camera, const Stars &stars, @@ -16,35 +19,63 @@ Attitude DavenportQAlgorithm::Go(const Camera &camera, assert(stars.size() >= 2); // attitude profile matrix - Eigen::Matrix3f B; + #ifdef LOST_FLOAT_MODE + Eigen::Matrix3f B; + #else + Eigen::Matrix3d B; + #endif + + B.setZero(); for (const StarIdentifier &s : starIdentifiers) { Star bStar = stars[s.starIndex]; Vec3 bStarSpatial = camera.CameraToSpatial(bStar.position); - Eigen::Vector3f bi; + + #ifdef LOST_FLOAT_MODE + Eigen::Vector3f bi; + Eigen::Vector3f ri; + #else + Eigen::Vector3d bi; + Eigen::Vector3d ri; + #endif + bi << bStarSpatial.x, bStarSpatial.y, bStarSpatial.z; CatalogStar rStar = catalog[s.catalogIndex]; - Eigen::Vector3f ri; ri << rStar.spatial.x, rStar.spatial.y, rStar.spatial.z; - //Weight = 1 (can be changed later, in which case we want to make a vector to hold all weights {ai}) // Calculate matrix B = sum({ai}{bi}{ri}T) B += ri * bi.transpose() * s.weight; } // S = B + Transpose(B) - Eigen::Matrix3f S = B + B.transpose(); + #ifdef LOST_FLOAT_MODE + Eigen::Matrix3f S = B + B.transpose(); + #else + Eigen::Matrix3d S = B + B.transpose(); + #endif + //sigma = B[0][0] + B[1][1] + B[2][2] - float sigma = B.trace(); + decimal sigma = B.trace(); + //Z = [[B[1][2] - B[2][1]], [B[2][0] - B[0][2]], [B[0][1] - B[1][0]]] - Eigen::Vector3f Z; + #ifdef LOST_FLOAT_MODE + Eigen::Vector3f Z; + #else + Eigen::Vector3d Z; + #endif + Z << B(1,2) - B(2,1), B(2,0) - B(0,2), B(0,1) - B(1,0); // K = [[[sigma], [Z[0]], [Z[1]], [Z[2]]], [[Z[0]], [S[0][0] - sigma], [S[0][1]], [S[0][2]]], [[Z[1]], [S[1][0]], [S[1][1] - sigma], [S[1][2]]], [[Z[2]], [S[2][0]], [S[2][1]], [S[2][2] - sigma]]] - Eigen::Matrix4f K; + #ifdef LOST_FLOAT_MODE + Eigen::Matrix4f K; + #else + Eigen::Matrix4d K; + #endif + K << sigma, Z(0), Z(1), Z(2), Z(0), S(0,0) - sigma, S(0,1), S(0,2), Z(1), S(1,0), S(1,1) - sigma, S(1,2), @@ -52,11 +83,18 @@ Attitude DavenportQAlgorithm::Go(const Camera &camera, //Find eigenvalues of K, store the largest one as lambda //find the maximum index - Eigen::EigenSolver solver(K); - Eigen::Vector4cf values = solver.eigenvalues(); - Eigen::Matrix4cf vectors = solver.eigenvectors(); + #ifdef LOST_FLOAT_MODE + Eigen::EigenSolver solver(K); + Eigen::Vector4cf values = solver.eigenvalues(); + Eigen::Matrix4cf vectors = solver.eigenvectors(); + #else + Eigen::EigenSolver solver(K); + Eigen::Vector4cd values = solver.eigenvalues(); + Eigen::Matrix4cd vectors = solver.eigenvectors(); + #endif + int maxIndex = 0; - float maxEigenvalue = values(0).real(); + decimal maxEigenvalue = values(0).real(); for (int i = 1; i < values.size(); i++) { if (values(i).real() > maxEigenvalue) { maxIndex = i; @@ -65,7 +103,12 @@ Attitude DavenportQAlgorithm::Go(const Camera &camera, } //The return quaternion components = eigenvector assocaited with lambda - Eigen::Vector4cf maxEigenvector = vectors.col(maxIndex); + #ifdef LOST_FLOAT_MODE + Eigen::Vector4cf maxEigenvector = vectors.col(maxIndex); + #else + Eigen::Vector4cd maxEigenvector = vectors.col(maxIndex); + #endif + // IMPORTANT: The matrix K is symmetric -- clearly first row and first column are equal. // Furthermore, S is symmetric because s_i,j = b_i,j + b_j,i and s_j,i=b_j,i + b_i,j, so the // bottom right 3x3 block of K is symmetric too. Thus all its eigenvalues and eigenvectors @@ -117,19 +160,19 @@ Attitude TriadAlgorithm::Go(const Camera &camera, * Characteristic polynomial of the quest K-matrix * @see equation 19b of https://arc.aiaa.org/doi/pdf/10.2514/1.62549 */ -float QuestCharPoly(float x, float a, float b, float c, float d, float s) {return (pow(x,2)-a) * (pow(x,2)-b) - (c*x) + (c*s) - d;} +decimal QuestCharPoly(decimal x, decimal a, decimal b, decimal c, decimal d, decimal s) {return (DECIMAL_POW(x,2)-a) * (DECIMAL_POW(x,2)-b) - (c*x) + (c*s) - d;} /** * Derivitive of the characteristic polynomial of the quest K-matrix */ -float QuestCharPolyPrime(float x, float a, float b, float c) {return 4*pow(x,3) - 2*(a+b)*x - c;} +decimal QuestCharPolyPrime(decimal x, decimal a, decimal b, decimal c) {return 4*DECIMAL_POW(x,3) - 2*(a+b)*x - c;} /** * 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) { - float height; +decimal QuestEigenvalueEstimator(decimal guess, decimal a, decimal b, decimal c, decimal d, decimal s) { + decimal height; do { height = QuestCharPoly(guess, a, b, c, d, s) / QuestCharPolyPrime(guess, a, b, c); guess -= height; @@ -149,7 +192,7 @@ Attitude QuestAlgorithm::Go(const Camera &camera, assert(stars.size() >= 2); // initial guess for eigenvalue (sum of the weights) - float guess = 0; + decimal guess = 0; // attitude profile matrix Mat3 B = {0, 0, 0, 0, 0, 0, 0, 0, 0}; @@ -171,7 +214,7 @@ Attitude QuestAlgorithm::Go(const Camera &camera, // S = B + Transpose(B) Mat3 S = B + B.Transpose(); //sigma = B[0][0] + B[1][1] + B[2][2] - float sigma = B.Trace(); + decimal sigma = B.Trace(); //Z = [[B[1][2] - B[2][1]], [B[2][0] - B[0][2]], [B[0][1] - B[1][0]]] Vec3 Z = { B.At(1,2) - B.At(2,1), @@ -180,23 +223,23 @@ Attitude QuestAlgorithm::Go(const Camera &camera, }; // calculate coefficients for characteristic polynomial - float delta = S.Det(); - float kappa = (S.Inverse() * delta).Trace(); - float a = pow(sigma,2) - kappa; - float b = pow(sigma,2) + (Z * Z); - float c = delta + (Z * S * Z); - float d = Z * (S * S) * Z; + decimal delta = S.Det(); + decimal kappa = (S.Inverse() * delta).Trace(); + decimal a = DECIMAL_POW(sigma,2) - kappa; + decimal b = DECIMAL_POW(sigma,2) + (Z * Z); + decimal c = delta + (Z * S * Z); + decimal d = Z * (S * S) * Z; // Newton-Raphson method for estimating the largest eigenvalue - float eig = QuestEigenvalueEstimator(guess, a, b, c, d, sigma); + decimal eig = QuestEigenvalueEstimator(guess, a, b, c, d, sigma); // solve for the optimal quaternion: from https://ahrs.readthedocs.io/en/latest/filters/quest.html - float alpha = pow(eig,2) - pow(sigma, 2) + kappa; - float beta = eig - sigma; - float gamma = (eig + sigma) * alpha - delta; + decimal alpha = DECIMAL_POW(eig,2) - DECIMAL_POW(sigma, 2) + kappa; + decimal beta = eig - sigma; + decimal gamma = (eig + sigma) * alpha - delta; Vec3 X = ((kIdentityMat3 * alpha) + (S * beta) + (S * S)) * Z; - float scalar = 1 / sqrt(pow(gamma,2) + X.MagnitudeSq()); + decimal scalar = 1 / DECIMAL_SQRT(DECIMAL_POW(gamma,2) + X.MagnitudeSq()); X = X * scalar; gamma *= scalar; diff --git a/src/attitude-utils.cpp b/src/attitude-utils.cpp index fdc8f8c1..9921a0e4 100644 --- a/src/attitude-utils.cpp +++ b/src/attitude-utils.cpp @@ -2,8 +2,10 @@ #include #include +#include #include +#include "decimal.hpp" #include "serialize-helpers.hpp" namespace lost { @@ -41,12 +43,12 @@ Quaternion::Quaternion(const Vec3 &input) { } /// Create a quaternion which represents a rotation of theta around the axis input -Quaternion::Quaternion(const Vec3 &input, float theta) { - real = cos(theta/2); +Quaternion::Quaternion(const Vec3 &input, decimal theta) { + real = DECIMAL_COS(theta/2); // the compiler will optimize it. Right? - i = input.x * sin(theta/2); - j = input.y * sin(theta/2); - k = input.z * sin(theta/2); + i = input.x * DECIMAL_SIN(theta/2); + j = input.y * DECIMAL_SIN(theta/2); + k = input.z * DECIMAL_SIN(theta/2); } /// Rotate a 3d vector according to the rotation represented by the quaternion. @@ -56,24 +58,24 @@ Vec3 Quaternion::Rotate(const Vec3 &input) const { } /// How many radians the rotation represented by this quaternion has. -float Quaternion::Angle() const { +decimal Quaternion::Angle() const { if (real <= -1) { return 0; // 180*2=360=0 } // TODO: we shouldn't need this nonsense, right? how come acos sometimes gives nan? (same as in AngleUnit) - return (real >= 1 ? 0 : acos(real))*2; + return (real >= 1 ? 0 : DECIMAL_ACOS(real))*2; } -float Quaternion::SmallestAngle() const { - float rawAngle = Angle(); - return rawAngle > M_PI - ? 2*M_PI - rawAngle +decimal Quaternion::SmallestAngle() const { + decimal rawAngle = Angle(); + return rawAngle > DECIMAL_M_PI + ? 2*DECIMAL_M_PI - rawAngle : rawAngle; } -void Quaternion::SetAngle(float newAngle) { - real = cos(newAngle/2); - SetVector(Vector().Normalize() * sin(newAngle/2)); +void Quaternion::SetAngle(decimal newAngle) { + real = DECIMAL_COS(newAngle/2); + SetVector(Vector().Normalize() * DECIMAL_SIN(newAngle/2)); } EulerAngles Quaternion::ToSpherical() const { @@ -84,21 +86,21 @@ EulerAngles Quaternion::ToSpherical() const { // and 2, we store the conjugate of the quaternion (double check why?), which means we need to // invert the final de and roll terms, as well as negate all the terms involving a mix between // the real and imaginary parts. - float ra = atan2(2*(-real*k+i*j), 1-2*(j*j+k*k)); + decimal ra = DECIMAL_ATAN2(2*(-real*k+i*j), 1-2*(j*j+k*k)); if (ra < 0) - ra += 2*M_PI; - float de = -asin(2*(-real*j-i*k)); // allow de to be positive or negaive, as is convention - float roll = -atan2(2*(-real*i+j*k), 1-2*(i*i+j*j)); + ra += 2*DECIMAL_M_PI; + decimal de = -DECIMAL_ASIN(2*(-real*j-i*k)); // allow de to be positive or negaive, as is convention + decimal roll = -DECIMAL_ATAN2(2*(-real*i+j*k), 1-2*(i*i+j*j)); if (roll < 0) - roll += 2*M_PI; + roll += 2*DECIMAL_M_PI; return EulerAngles(ra, de, roll); } -Quaternion SphericalToQuaternion(float ra, float dec, float roll) { - assert(roll >= 0.0 && roll <= 2*M_PI); - assert(ra >= 0 && ra <= 2*M_PI); - assert(dec >= -M_PI && dec <= M_PI); +Quaternion SphericalToQuaternion(decimal ra, decimal dec, decimal roll) { + assert(roll >= DECIMAL(0.0) && roll <= 2*DECIMAL_M_PI); + assert(ra >= DECIMAL(0.0) && ra <= 2*DECIMAL_M_PI); + assert(dec >= -DECIMAL_M_PI && dec <= DECIMAL_M_PI); // when we are modifying the coordinate axes, the quaternion multiplication works so that the // rotations are applied from left to right. This is the opposite as for modifying vectors. @@ -115,7 +117,7 @@ Quaternion SphericalToQuaternion(float ra, float dec, float roll) { } /// Whether the quaternion is a unit quaternion. All quaternions representing rotations should be units. -bool Quaternion::IsUnit(float tolerance) const { +bool Quaternion::IsUnit(decimal tolerance) const { return abs(i*i+j*j+k*k+real*real - 1) < tolerance; } @@ -131,82 +133,82 @@ Quaternion Quaternion::Canonicalize() const { } /// Convert from right ascension & declination to a 3d point on the unit sphere. -Vec3 SphericalToSpatial(float ra, float de) { +Vec3 SphericalToSpatial(decimal ra, decimal de) { return { - cos(ra)*cos(de), - sin(ra)*cos(de), - sin(de), + DECIMAL_COS(ra)*DECIMAL_COS(de), + DECIMAL_SIN(ra)*DECIMAL_COS(de), + DECIMAL_SIN(de), }; } /// Convert from a 3d point on the unit sphere to right ascension & declination. -void SpatialToSpherical(const Vec3 &vec, float *ra, float *de) { - *ra = atan2(vec.y, vec.x); +void SpatialToSpherical(const Vec3 &vec, decimal *ra, decimal *de) { + *ra = DECIMAL_ATAN2(vec.y, vec.x); if (*ra < 0) - *ra += M_PI*2; - *de = asin(vec.z); + *ra += DECIMAL_M_PI*2; + *de = DECIMAL_ASIN(vec.z); } -float RadToDeg(float rad) { - return rad*180.0/M_PI; +decimal RadToDeg(decimal rad) { + return rad*DECIMAL(180.0)/DECIMAL_M_PI; } -float DegToRad(float deg) { - return deg/180.0*M_PI; +decimal DegToRad(decimal deg) { + return deg/DECIMAL(180.0)*DECIMAL_M_PI; } -float RadToArcSec(float rad) { - return RadToDeg(rad) * 3600.0; +decimal RadToArcSec(decimal rad) { + return RadToDeg(rad) * DECIMAL(3600.0); } -float ArcSecToRad(float arcSec) { - return DegToRad(arcSec / 3600.0); +decimal ArcSecToRad(decimal arcSec) { + return DegToRad(arcSec / DECIMAL(3600.0)); } -float FloatModulo(float x, float mod) { +decimal DecimalModulo(decimal x, decimal mod) { // first but not last chatgpt generated code in lost: - float result = x - mod * floor(x / mod); + decimal result = x - mod * DECIMAL_FLOOR(x / mod); return result >= 0 ? result : result + mod; } /// The square of the magnitude -float Vec3::MagnitudeSq() const { - return fma(x,x,fma(y,y, z*z)); +decimal Vec3::MagnitudeSq() const { + return DECIMAL_FMA(x,x,DECIMAL_FMA(y,y, z*z)); } /// The square of the magnitude -float Vec2::MagnitudeSq() const { - return fma(x,x, y*y); +decimal Vec2::MagnitudeSq() const { + return DECIMAL_FMA(x,x, y*y); } -float Vec3::Magnitude() const { - return hypot(hypot(x, y), z); // not sure if this is faster than a simple sqrt, but it does have less error? +decimal Vec3::Magnitude() const { + return DECIMAL_HYPOT(DECIMAL_HYPOT(x, y), z); // not sure if this is faster than a simple sqrt, but it does have less error? } -float Vec2::Magnitude() const { - return hypot(x, y); +decimal Vec2::Magnitude() const { + return DECIMAL_HYPOT(x, y); } /// Create a vector pointing in the same direction with magnitude 1 Vec3 Vec3::Normalize() const { - float mag = Magnitude(); + decimal mag = Magnitude(); return { x/mag, y/mag, z/mag, }; } /// Dot product -float Vec3::operator*(const Vec3 &other) const { - return fma(x,other.x, fma(y,other.y, z*other.z)); +decimal Vec3::operator*(const Vec3 &other) const { + return DECIMAL_FMA(x,other.x, DECIMAL_FMA(y,other.y, z*other.z)); } /// Dot product -Vec2 Vec2::operator*(const float &other) const { +Vec2 Vec2::operator*(const decimal &other) const { return { x*other, y*other }; } /// Vector-Scalar multiplication -Vec3 Vec3::operator*(const float &other) const { +Vec3 Vec3::operator*(const decimal &other) const { return { x*other, y*other, z*other }; } @@ -253,7 +255,7 @@ Vec3 Vec3::operator*(const Mat3 &other) const { } /// Access the i,j-th element of the matrix -float Mat3::At(int i, int j) const { +decimal Mat3::At(int i, int j) const { return x[3*i+j]; } @@ -297,7 +299,7 @@ Vec3 Mat3::operator*(const Vec3 &vec) const { } /// Matrix-Scalar multiplication -Mat3 Mat3::operator*(const float &s) const { +Mat3 Mat3::operator*(const decimal &s) const { return { s*At(0,0), s*At(0,1), s*At(0,2), s*At(1,0), s*At(1,1), s*At(1,2), @@ -315,19 +317,19 @@ Mat3 Mat3::Transpose() const { } /// Trace of a matrix -float Mat3::Trace() const { +decimal Mat3::Trace() const { return At(0,0) + At(1,1) + At(2,2); } /// Determinant of a matrix -float Mat3::Det() const { +decimal Mat3::Det() const { return (At(0,0) * (At(1,1)*At(2,2) - At(2,1)*At(1,2))) - (At(0,1) * (At(1,0)*At(2,2) - At(2,0)*At(1,2))) + (At(0,2) * (At(1,0)*At(2,1) - At(2,0)*At(1,1))); } /// Inverse of a matrix Mat3 Mat3::Inverse() const { // https://ardoris.wordpress.com/2008/07/18/general-formula-for-the-inverse-of-a-3x3-matrix/ - float scalar = 1 / Det(); + decimal scalar = 1 / Det(); Mat3 res = { At(1,1)*At(2,2) - At(1,2)*At(2,1), At(0,2)*At(2,1) - At(0,1)*At(2,2), At(0,1)*At(1,2) - At(0,2)*At(1,1), @@ -367,9 +369,9 @@ Quaternion DCMToQuaternion(const Mat3 &dcm) { // the DCM itself does Vec3 oldXAxis = Vec3({1, 0, 0}); Vec3 newXAxis = dcm.Column(0); // this is where oldXAxis is mapped to - assert(abs(newXAxis.Magnitude()-1) < 0.001); + assert(DECIMAL_ABS(newXAxis.Magnitude()-1) < DECIMAL(0.001)); Vec3 xAlignAxis = oldXAxis.CrossProduct(newXAxis).Normalize(); - float xAlignAngle = AngleUnit(oldXAxis, newXAxis); + decimal xAlignAngle = AngleUnit(oldXAxis, newXAxis); Quaternion xAlign(xAlignAxis, xAlignAngle); // Make a quaternion that will rotate the Y-axis into place @@ -450,22 +452,22 @@ bool Attitude::IsKnown() const { /// Serialize a Vec3 to buffer. Takes up space according to SerializeLengthVec3 void SerializeVec3(SerializeContext *ser, const Vec3 &vec) { - SerializePrimitive(ser, vec.x); - SerializePrimitive(ser, vec.y); - SerializePrimitive(ser, vec.z); + SerializePrimitive(ser, vec.x); + SerializePrimitive(ser, vec.y); + SerializePrimitive(ser, vec.z); } Vec3 DeserializeVec3(DeserializeContext *des) { Vec3 result = { - DeserializePrimitive(des), - DeserializePrimitive(des), - DeserializePrimitive(des), + DeserializePrimitive(des), + DeserializePrimitive(des), + DeserializePrimitive(des), }; return result; } /// Calculate the inner angle, in radians, between two vectors. -float Angle(const Vec3 &vec1, const Vec3 &vec2) { +decimal Angle(const Vec3 &vec1, const Vec3 &vec2) { return AngleUnit(vec1.Normalize(), vec2.Normalize()); } @@ -474,10 +476,10 @@ float Angle(const Vec3 &vec1, const Vec3 &vec2) { * Slightly faster than Angle() * @warn If the vectors are not already unit vectors, will return the wrong result! */ -float AngleUnit(const Vec3 &vec1, const Vec3 &vec2) { - float dot = vec1*vec2; +decimal AngleUnit(const Vec3 &vec1, const Vec3 &vec2) { + decimal dot = vec1*vec2; // TODO: we shouldn't need this nonsense, right? how come acos sometimes gives nan? - return dot >= 1 ? 0 : dot <= -1 ? M_PI-0.0000001 : acos(dot); + return dot >= 1 ? 0 : dot <= -1 ? DECIMAL_M_PI-DECIMAL(0.0000001) : DECIMAL_ACOS(dot); } } diff --git a/src/attitude-utils.hpp b/src/attitude-utils.hpp index 842c1220..0935afc8 100644 --- a/src/attitude-utils.hpp +++ b/src/attitude-utils.hpp @@ -5,6 +5,7 @@ #include #include "serialize-helpers.hpp" +#include "decimal.hpp" namespace lost { @@ -12,58 +13,58 @@ namespace lost { // to Quaterinon, and another storing as Quaternion and converting to Euler. But abstract classes // make everything more annoying, because you need vectors of pointers...ugh! -/// A two dimensional vector with floating point components +/// A two dimensional vector with decimaling point components struct Vec2 { - float x; - float y; + decimal x; + decimal y; - float Magnitude() const; - float MagnitudeSq() const; + decimal Magnitude() const; + decimal MagnitudeSq() const; Vec2 Normalize() const; - float operator*(const Vec2 &) const; - Vec2 operator*(const float &) const; + decimal operator*(const Vec2 &) const; + Vec2 operator*(const decimal &) const; Vec2 operator-(const Vec2 &) const; Vec2 operator+(const Vec2 &) const; }; class Mat3; // define above so we can use in Vec3 class -/// Three dimensional vector with floating point components +/// Three dimensional vector with decimaling point components class Vec3 { public: - float x; - float y; - float z; + decimal x; + decimal y; + decimal z; - float Magnitude() const; - float MagnitudeSq() const; + decimal Magnitude() const; + decimal MagnitudeSq() const; Vec3 Normalize() const; - float operator*(const Vec3 &) const; - Vec3 operator*(const float &) const; + decimal operator*(const Vec3 &) const; + Vec3 operator*(const decimal &) const; Vec3 operator*(const Mat3 &) const; Vec3 operator-(const Vec3 &) const; Vec3 CrossProduct(const Vec3 &) const; Mat3 OuterProduct(const Vec3 &) const; }; -/// 3x3 vector with floating point components +/// 3x3 vector with decimaling point components class Mat3 { public: - float x[9]; + decimal x[9]; - float At(int i, int j) const; + decimal At(int i, int j) const; Mat3 operator+(const Mat3 &) const; Mat3 operator*(const Mat3 &) const; Vec3 operator*(const Vec3 &) const; - Mat3 operator*(const float &) const; + Mat3 operator*(const decimal &) const; Mat3 Transpose() const; Vec3 Column(int) const; Vec3 Row(int) const; - float Trace() const; - float Det() const; + decimal Trace() const; + decimal Det() const; Mat3 Inverse() const; }; @@ -72,8 +73,8 @@ extern const Mat3 kIdentityMat3; void SerializeVec3(SerializeContext *, const Vec3 &); Vec3 DeserializeVec3(DeserializeContext *des); -float Distance(const Vec2 &, const Vec2 &); -float Distance(const Vec3 &, const Vec3 &); +decimal Distance(const Vec2 &, const Vec2 &); +decimal Distance(const Vec3 &, const Vec3 &); /** * A "human-readable" way to represent a 3d rotation or orientation. @@ -82,15 +83,15 @@ float Distance(const Vec3 &, const Vec3 &); */ class EulerAngles { public: - EulerAngles(float ra, float de, float roll) + EulerAngles(decimal ra, decimal de, decimal roll) : ra(ra), de(de), roll(roll) { }; /// Right ascension. How far we yaw left. Yaw is performed first. - float ra; + decimal ra; /// Declination. How far we pitch up (or down if negative). Pitch is performed second, after yaw. - float de; + decimal de; /// How far we roll counterclockwise. Roll is performed last (after yaw and pitch). - float roll; + decimal roll; }; /// A quaternion is a common way to represent a 3d rotation. @@ -98,9 +99,9 @@ class Quaternion { public: Quaternion() = default; explicit Quaternion(const Vec3 &); - Quaternion(const Vec3 &, float); + Quaternion(const Vec3 &, decimal); - Quaternion(float real, float i, float j, float k) + Quaternion(decimal real, decimal i, decimal j, decimal k) : real(real), i(i), j(j), k(k) { }; Quaternion operator*(const Quaternion &other) const; @@ -108,19 +109,19 @@ class Quaternion { Vec3 Vector() const; void SetVector(const Vec3 &); Vec3 Rotate(const Vec3 &) const; - float Angle() const; + decimal 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); + decimal SmallestAngle() const; + void SetAngle(decimal); EulerAngles ToSpherical() const; - bool IsUnit(float tolerance) const; + bool IsUnit(decimal tolerance) const; Quaternion Canonicalize() const; - float real; - float i; - float j; - float k; + decimal real; + decimal i; + decimal j; + decimal k; }; // @@ -165,23 +166,23 @@ 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. -Quaternion SphericalToQuaternion(float ra, float dec, float roll); +Quaternion SphericalToQuaternion(decimal ra, decimal dec, decimal roll); /// returns unit vector -Vec3 SphericalToSpatial(float ra, float de); -void SpatialToSpherical(const Vec3 &, float *ra, float *de); +Vec3 SphericalToSpatial(decimal ra, decimal de); +void SpatialToSpherical(const Vec3 &, decimal *ra, decimal *de); /// angle between two vectors, using dot product and magnitude division -float Angle(const Vec3 &, const Vec3 &); +decimal Angle(const Vec3 &, const Vec3 &); /// angle between two vectors, /assuming/ that they are already unit length -float AngleUnit(const Vec3 &, const Vec3 &); +decimal AngleUnit(const Vec3 &, const Vec3 &); -float RadToDeg(float); -float DegToRad(float); -float RadToArcSec(float); -float ArcSecToRad(float); -/// Given a float, find it "modulo" another float, in the true mathematical sense (not remainder). +decimal RadToDeg(decimal); +decimal DegToRad(decimal); +decimal RadToArcSec(decimal); +decimal ArcSecToRad(decimal); +/// Given a decimal, find it "modulo" another decimal, in the true mathematical sense (not remainder). /// Always returns something in [0,mod) Eg -0.8 mod 0.6 = 0.4 -float FloatModulo(float x, float mod); +decimal DecimalModulo(decimal x, decimal mod); // TODO: quaternion and euler angle conversion, conversion between ascension/declination to rec9tu diff --git a/src/camera.cpp b/src/camera.cpp index af45d953..956c2553 100644 --- a/src/camera.cpp +++ b/src/camera.cpp @@ -4,6 +4,7 @@ #include #include "attitude-utils.hpp" +#include "decimal.hpp" namespace lost { @@ -16,10 +17,10 @@ Vec2 Camera::SpatialToCamera(const Vec3 &vector) const { assert(vector.x > 0); // TODO: is there any sort of accuracy problem when vector.y and vector.z are small? - float focalFactor = focalLength/vector.x; + decimal focalFactor = focalLength/vector.x; - float yPixel = vector.y*focalFactor; - float zPixel = vector.z*focalFactor; + decimal yPixel = vector.y*focalFactor; + decimal zPixel = vector.z*focalFactor; return { -yPixel + xCenter, -zPixel + yCenter }; } @@ -36,8 +37,8 @@ Vec3 Camera::CameraToSpatial(const Vec2 &vector) const { // isn't it interesting: To convert from center-based to left-corner-based coordinates is the // same formula; f(x)=f^{-1}(x) ! - float xPixel = -vector.x + xCenter; - float yPixel = -vector.y + yCenter; + decimal xPixel = -vector.x + xCenter; + decimal yPixel = -vector.y + yCenter; return { 1, @@ -54,15 +55,15 @@ bool Camera::InSensor(const Vec2 &vector) const { && vector.y >= 0 && vector.y <= yResolution; } -float FovToFocalLength(float xFov, float xResolution) { - return xResolution / 2.0f / tan(xFov/2); +decimal FovToFocalLength(decimal xFov, decimal xResolution) { + return xResolution / DECIMAL(2.0) / DECIMAL_TAN(xFov/2); } -float FocalLengthToFov(float focalLength, float xResolution, float pixelSize) { - return atan(xResolution/2 * pixelSize / focalLength) * 2; +decimal FocalLengthToFov(decimal focalLength, decimal xResolution, decimal pixelSize) { + return DECIMAL_ATAN(xResolution/2 * pixelSize / focalLength) * 2; } -float Camera::Fov() const { +decimal Camera::Fov() const { return FocalLengthToFov(focalLength, xResolution, 1.0); } diff --git a/src/camera.hpp b/src/camera.hpp index 2c076cfd..703c0f87 100644 --- a/src/camera.hpp +++ b/src/camera.hpp @@ -13,16 +13,16 @@ class Camera { /** * @param xCenter,yCenter The "principal point" of the camera. In an ideal camera, just half the resolution, but physical cameras often have a bit of offset. */ - Camera(float focalLength, - float xCenter, float yCenter, + Camera(decimal focalLength, + decimal xCenter, decimal yCenter, int xResolution, int yResolution) : focalLength(focalLength), xCenter(xCenter), yCenter(yCenter), xResolution(xResolution), yResolution(yResolution) {}; - Camera(float focalLength, int xResolution, int yResolution) + Camera(decimal focalLength, int xResolution, int yResolution) : Camera(focalLength, - xResolution / (float) 2.0, yResolution / (float) 2.0, + xResolution / DECIMAL(2.0), yResolution / DECIMAL(2.0), xResolution, yResolution) {}; Vec2 SpatialToCamera(const Vec3 &) const; @@ -30,7 +30,7 @@ class Camera { // 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; + // void CoordinateAngles(const Vec2 &vector, decimal *ra, decimal *de) const; bool InSensor(const Vec2 &vector) const; @@ -39,20 +39,20 @@ class Camera { /// Height of the sensor in pixels int YResolution() const { return yResolution; }; /// Focal length in pixels - float FocalLength() const { return focalLength; }; + decimal FocalLength() const { return focalLength; }; /// Horizontal field of view in radians - float Fov() const; + decimal Fov() const; - void SetFocalLength(float focalLength) { this->focalLength = focalLength; } + void SetFocalLength(decimal focalLength) { this->focalLength = focalLength; } private: // TODO: distortion - float focalLength; - float xCenter; float yCenter; + decimal focalLength; + decimal xCenter; decimal yCenter; int xResolution; int yResolution; }; -float FovToFocalLength(float xFov, float xResolution); +decimal FovToFocalLength(decimal xFov, decimal xResolution); } diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 0ddb6103..e9c05302 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -2,14 +2,16 @@ #include #include -#include #include +#include #include #include #include #include +#include "decimal.hpp" + namespace lost { // DUMMY @@ -19,7 +21,7 @@ std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, in unsigned int randomSeed = 123456; for (int i = 0; i < numStars; i++) { - result.push_back(Star(rand_r(&randomSeed) % imageWidth, rand_r(&randomSeed) % imageHeight, 10.0)); + result.push_back(Star(rand_r(&randomSeed) % imageWidth, rand_r(&randomSeed) % imageHeight, DECIMAL(10.0))); } return result; @@ -39,11 +41,11 @@ int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) { int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { // code here, duh long total = imageWidth * imageHeight; - //float top = 255; - float sumB = 0; - float sum1 = 0; - float wB = 0; - float maximum = 0; + //decimal top = 255; + decimal sumB = 0; + decimal sum1 = 0; + decimal wB = 0; + decimal maximum = 0; int level = 0; // make the histogram (array length 256) int histogram[256]; @@ -57,12 +59,12 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { sum1 += i * histogram[i]; } for (int i = 0; i < 256; i ++) { - float wF = total - wB; + decimal wF = total - wB; //std::cout << "wF\n" << wB << "\n"; //std::cout << "wB\n" << wF << "\n"; if (wB > 0 && wF > 0) { - float mF = (sum1 - sumB) / wF; - float val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF); + decimal mF = (sum1 - sumB) / wF; + decimal val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF); //std::cout << val << "\n"; if (val >= maximum) { level = i; @@ -78,38 +80,38 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { // a simple, but well tested thresholding algorithm that works well with star images int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight) { unsigned long totalMag = 0; - float std = 0; + decimal std = 0; long totalPixels = imageHeight * imageWidth; for (long i = 0; i < totalPixels; i++) { totalMag += image[i]; } - float mean = totalMag / totalPixels; + decimal mean = totalMag / totalPixels; for (long i = 0; i < totalPixels; i++) { - std += std::pow(image[i] - mean, 2); + std += DECIMAL_POW(image[i] - mean, 2); } - std = std::sqrt(std / totalPixels); + std = DECIMAL_SQRT(std / totalPixels); return mean + (std * 5); } // basic thresholding, but do it faster (trade off of some accuracy?) int BasicThresholdOnePass(unsigned char *image, int imageWidth, int imageHeight) { unsigned long totalMag = 0; - float std = 0; - float sq_totalMag = 0; + decimal std = 0; + decimal sq_totalMag = 0; long totalPixels = imageHeight * imageWidth; for (long i = 0; i < totalPixels; i++) { totalMag += image[i]; sq_totalMag += image[i] * image[i]; } - float mean = totalMag / totalPixels; - float variance = (sq_totalMag / totalPixels) - (mean * mean); - std = std::sqrt(variance); + decimal mean = totalMag / totalPixels; + decimal variance = (sq_totalMag / totalPixels) - (mean * mean); + std = DECIMAL_SQRT(variance); return mean + (std * 5); } struct CentroidParams { - float yCoordMagSum; - float xCoordMagSum; + decimal yCoordMagSum; + decimal xCoordMagSum; long magSum; int xMin; int xMax; @@ -182,11 +184,11 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi yDiameter = (p.yMax - p.yMin) + 1; //use the sums to finish CoG equation and add stars to the result - float xCoord = (p.xCoordMagSum / (p.magSum * 1.0)); - float yCoord = (p.yCoordMagSum / (p.magSum * 1.0)); + decimal xCoord = (p.xCoordMagSum / (p.magSum * DECIMAL(1.0))); + decimal yCoord = (p.yCoordMagSum / (p.magSum * DECIMAL(1.0))); if (p.isValid) { - result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f, ((float)(xDiameter))/2.0f, ((float)(yDiameter))/2.0f, p.checkedIndices.size() - sizeBefore)); + result.push_back(Star(xCoord + DECIMAL(0.5), yCoord + DECIMAL(0.5), (xDiameter)/DECIMAL(2.0), (yDiameter)/DECIMAL(2.0), p.checkedIndices.size() - sizeBefore)); } } } @@ -195,7 +197,7 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi //Determines how accurate and how much iteration is done by the IWCoG algorithm, //smaller means more accurate and more iterations. -float iWCoGMinChange = 0.0002; +decimal iWCoGMinChange = DECIMAL(0.0002); struct IWCoGParams { int xMin; @@ -254,12 +256,12 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im p.maxIntensity = 0; int xDiameter = 0; int yDiameter = 0; - float yWeightedCoordMagSum = 0; - float xWeightedCoordMagSum = 0; - float weightedMagSum = 0; - float fwhm; //fwhm variable - float standardDeviation; - float w; //weight value + decimal yWeightedCoordMagSum = 0; + decimal xWeightedCoordMagSum = 0; + decimal weightedMagSum = 0; + decimal fwhm; //fwhm variable + decimal standardDeviation; + decimal w; //weight value p.xMax = i % imageWidth; p.xMin = i % imageWidth; @@ -274,20 +276,20 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im yDiameter = (p.yMax - p.yMin) + 1; //calculate fwhm - float count = 0; + decimal count = 0; for (int j = 0; j < (int) starIndices.size(); j++) { if (image[starIndices.at(j)] > p.maxIntensity / 2) { count++; } } - fwhm = sqrt(count); - standardDeviation = fwhm / (2.0 * sqrt(2.0 * log(2.0))); - float modifiedStdDev = 2.0 * pow(standardDeviation, 2); - // TODO: Why are these floats? --Mark - float guessXCoord = (float) (p.guess % imageWidth); - float guessYCoord = (float) (p.guess / imageWidth); + fwhm = DECIMAL_SQRT(count); + standardDeviation = fwhm / (DECIMAL(2.0) * DECIMAL_SQRT(DECIMAL(2.0) * DECIMAL_LOG(2.0))); + decimal modifiedStdDev = DECIMAL(2.0) * DECIMAL_POW(standardDeviation, 2); + // TODO: Why are these decimals? --Mark + decimal guessXCoord = (p.guess % imageWidth); + decimal guessYCoord = (p.guess / imageWidth); //how much our new centroid estimate changes w each iteration - float change = INFINITY; + decimal change = INFINITY; int stop = 0; //while we see some large enough change in estimated, maybe make it a global variable while (change > iWCoGMinChange && stop < 100000) { @@ -298,16 +300,16 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im stop++; for (long j = 0; j < (long)starIndices.size(); j++) { //calculate w - float currXCoord = (float) (starIndices.at(j) % imageWidth); - float currYCoord = (float) (starIndices.at(j) / imageWidth); - w = p.maxIntensity * exp(-1.0 * ((pow(currXCoord - guessXCoord, 2) / modifiedStdDev) + (pow(currYCoord - guessYCoord, 2) / modifiedStdDev))); + decimal currXCoord = starIndices.at(j) % imageWidth; + decimal currYCoord = starIndices.at(j) / imageWidth; + w = p.maxIntensity * DECIMAL_EXP(DECIMAL(-1.0) * ((DECIMAL_POW(currXCoord - guessXCoord, 2) / modifiedStdDev) + (DECIMAL_POW(currYCoord - guessYCoord, 2) / modifiedStdDev))); - xWeightedCoordMagSum += w * currXCoord * ((float) image[starIndices.at(j)]); - yWeightedCoordMagSum += w * currYCoord * ((float) image[starIndices.at(j)]); - weightedMagSum += w * ((float) image[starIndices.at(j)]); + xWeightedCoordMagSum += w * currXCoord * DECIMAL(image[starIndices.at(j)]); + yWeightedCoordMagSum += w * currYCoord * DECIMAL(image[starIndices.at(j)]); + weightedMagSum += w * DECIMAL(image[starIndices.at(j)]); } - float xTemp = xWeightedCoordMagSum / weightedMagSum; - float yTemp = yWeightedCoordMagSum / weightedMagSum; + decimal xTemp = xWeightedCoordMagSum / weightedMagSum; + decimal yTemp = yWeightedCoordMagSum / weightedMagSum; change = abs(guessXCoord - xTemp) + abs(guessYCoord - yTemp); @@ -315,7 +317,7 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im guessYCoord = yTemp; } if (p.isValid) { - result.push_back(Star(guessXCoord + 0.5f, guessYCoord + 0.5f, ((float)(xDiameter))/2.0f, ((float)(yDiameter))/2.0f, starIndices.size())); + result.push_back(Star(guessXCoord + DECIMAL(0.5), guessYCoord + DECIMAL(0.5), xDiameter/DECIMAL(2.0), yDiameter/DECIMAL(2.0), starIndices.size())); } } } diff --git a/src/database-options.hpp b/src/database-options.hpp index bbad9e8f..0f609004 100644 --- a/src/database-options.hpp +++ b/src/database-options.hpp @@ -1,14 +1,15 @@ // see pipeline-options.hpp for more information #include +#include "decimal.hpp" -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("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("swap-integer-endianness", bool , swapIntegerEndianness , false , atobool(optarg), true) -LOST_CLI_OPTION("swap-float-endianness" , bool , swapFloatEndianness , false , atobool(optarg), true) -LOST_CLI_OPTION("output" , std::string, outputPath , "-" , optarg , kNoDefaultArgument) +LOST_CLI_OPTION("min-mag" , decimal , minMag , 100 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("max-stars" , int , maxStars , 10000 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("min-separation" , decimal , minSeparation , 0.08 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("kvector" , bool , kvector , false , atobool(optarg), true) +LOST_CLI_OPTION("kvector-min-distance" , decimal , kvectorMinDistance , 0.5 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("kvector-max-distance" , decimal , kvectorMaxDistance , 15 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +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-decimal-endianness", bool , swapDecimalEndianness , false , atobool(optarg), true) +LOST_CLI_OPTION("output" , std::string, outputPath , "-" , optarg , kNoDefaultArgument) diff --git a/src/databases.cpp b/src/databases.cpp index 93302d61..150c7885 100644 --- a/src/databases.cpp +++ b/src/databases.cpp @@ -3,8 +3,8 @@ #include #include #include - #include +#include #include #include @@ -16,10 +16,14 @@ namespace lost { const int32_t PairDistanceKVectorDatabase::kMagicValue = 0x2536f009; +inline bool isFlagSet(uint32_t dbFlags, uint32_t flag) { + return (dbFlags & flag) != 0; +} + struct KVectorPair { int16_t index1; int16_t index2; - float distance; + decimal distance; }; bool CompareKVectorPairs(const KVectorPair &p1, const KVectorPair &p2) { @@ -33,8 +37,8 @@ bool CompareKVectorPairs(const KVectorPair &p1, const KVectorPair &p2) { | size | name | description | |---------------+------------+-------------------------------------------------------------| | 4 | numEntries | | - | sizeof float | min | minimum value contained in the database | - | sizeof float | max | max value contained in index | + | sizeof decimal | min | minimum value contained in the database | + | sizeof decimal | max | max value contained in index | | 4 | numBins | | | 4*(numBins+1) | bins | The `i'th bin (starting from zero) stores how many pairs of | | | | stars have a distance lesst han or equal to: | @@ -57,9 +61,9 @@ bool CompareKVectorPairs(const KVectorPair &p1, const KVectorPair &p2) { * @param numBins the number of "bins" the KVector should use. A higher number makes query results "tighter" but takes up more disk space. Usually should be set somewhat smaller than (max-min) divided by the "width" of the typical query. * @param buffer[out] index is written here. */ -void SerializeKVectorIndex(SerializeContext *ser, const std::vector &values, float min, float max, long numBins) { +void SerializeKVectorIndex(SerializeContext *ser, const std::vector &values, decimal min, decimal max, long numBins) { std::vector kVector(numBins+1); // We store sums before and after each bin - float binWidth = (max - min) / numBins; + decimal binWidth = (max - min) / numBins; // generate the k-vector part // Idea: When we find the first star that's across any bin boundary, we want to update all the newly sealed bins @@ -92,8 +96,8 @@ void SerializeKVectorIndex(SerializeContext *ser, const std::vector &valu // metadata fields SerializePrimitive(ser, values.size()); - SerializePrimitive(ser, min); - SerializePrimitive(ser, max); + SerializePrimitive(ser, min); + SerializePrimitive(ser, max); SerializePrimitive(ser, numBins); // kvector index field @@ -106,11 +110,11 @@ void SerializeKVectorIndex(SerializeContext *ser, const std::vector &valu KVectorIndex::KVectorIndex(DeserializeContext *des) { numValues = DeserializePrimitive(des); - min = DeserializePrimitive(des); - max = DeserializePrimitive(des); + min = DeserializePrimitive(des); + max = DeserializePrimitive(des); numBins = DeserializePrimitive(des); - assert(min >= 0.0f); + assert(min >= DECIMAL(0.0)); assert(max > min); binWidth = (max - min) / numBins; @@ -122,13 +126,13 @@ KVectorIndex::KVectorIndex(DeserializeContext *des) { * @param upperIndex[out] Is set to the index of the last returned value +1. * @return the index (starting from zero) of the first value matching the query */ -long KVectorIndex::QueryLiberal(float minQueryDistance, float maxQueryDistance, long *upperIndex) const { +long KVectorIndex::QueryLiberal(decimal minQueryDistance, decimal maxQueryDistance, long *upperIndex) const { assert(maxQueryDistance > minQueryDistance); if (maxQueryDistance >= max) { - maxQueryDistance = max - 0.00001; // TODO: better way to avoid hitting the bottom bin + maxQueryDistance = max - DECIMAL(0.00001); // TODO: better way to avoid hitting the bottom bin } if (minQueryDistance <= min) { - minQueryDistance = min + 0.00001; + minQueryDistance = min + DECIMAL(0.00001); } if (minQueryDistance > max || maxQueryDistance < min) { *upperIndex = 0; @@ -152,7 +156,7 @@ long KVectorIndex::QueryLiberal(float minQueryDistance, float maxQueryDistance, } /// return the lowest-indexed bin that contains the number of pairs with distance <= dist -long KVectorIndex::BinFor(float query) const { +long KVectorIndex::BinFor(decimal query) const { long result = (long)ceil((query - min) / binWidth); assert(result >= 0); assert(result <= numBins); @@ -168,7 +172,7 @@ long KVectorIndex::BinFor(float query) const { | sizeof kvectorIndex | kVectorIndex | Serialized KVector index | | 2*sizeof(int16)*numPairs | pairs | Bulk pair data | */ -std::vector CatalogToPairDistances(const Catalog &catalog, float minDistance, float maxDistance) { +std::vector CatalogToPairDistances(const Catalog &catalog, decimal minDistance, decimal maxDistance) { std::vector result; for (int16_t i = 0; i < (int16_t)catalog.size(); i++) { for (int16_t k = i+1; k < (int16_t)catalog.size(); k++) { @@ -176,7 +180,7 @@ std::vector CatalogToPairDistances(const Catalog &catalog, float mi KVectorPair pair = { i, k, AngleUnit(catalog[i].spatial, catalog[k].spatial) }; assert(isfinite(pair.distance)); assert(pair.distance >= 0); - assert(pair.distance <= M_PI); + assert(pair.distance <= DECIMAL_M_PI); if (pair.distance >= minDistance && pair.distance <= maxDistance) { // we'll sort later @@ -191,13 +195,13 @@ std::vector CatalogToPairDistances(const Catalog &catalog, float mi * Serialize a pair-distance KVector into buffer. * Use SerializeLengthPairDistanceKVector to determine how large the buffer needs to be. See command line documentation for other options. */ -void SerializePairDistanceKVector(SerializeContext *ser, const Catalog &catalog, float minDistance, float maxDistance, long numBins) { +void SerializePairDistanceKVector(SerializeContext *ser, const Catalog &catalog, decimal minDistance, decimal maxDistance, long numBins) { std::vector kVector(numBins+1); // numBins = length, all elements zero std::vector pairs = CatalogToPairDistances(catalog, minDistance, maxDistance); // sort pairs in increasing order. std::sort(pairs.begin(), pairs.end(), CompareKVectorPairs); - std::vector distances; + std::vector distances; for (const KVectorPair &pair : pairs) { distances.push_back(pair.distance); @@ -221,7 +225,7 @@ PairDistanceKVectorDatabase::PairDistanceKVectorDatabase(DeserializeContext *des } /// Return the value in the range [low,high] which is closest to num -float Clamp(float num, float low, float high) { +decimal Clamp(decimal num, decimal low, decimal high) { return num < low ? low : num > high ? high : num; } @@ -231,9 +235,9 @@ float Clamp(float num, float low, float high) { * @return A pointer to the start of the matched pairs. Each pair is stored as simply two 16-bit integers, each of which is a catalog index. (you must increment the pointer twice to get to the next pair). */ const int16_t *PairDistanceKVectorDatabase::FindPairsLiberal( - float minQueryDistance, float maxQueryDistance, const int16_t **end) const { + decimal minQueryDistance, decimal maxQueryDistance, const int16_t **end) const { - assert(maxQueryDistance <= M_PI); + assert(maxQueryDistance <= DECIMAL_M_PI); long upperIndex = -1; long lowerIndex = index.QueryLiberal(minQueryDistance, maxQueryDistance, &upperIndex); @@ -242,16 +246,16 @@ const int16_t *PairDistanceKVectorDatabase::FindPairsLiberal( } const int16_t *PairDistanceKVectorDatabase::FindPairsExact(const Catalog &catalog, - float minQueryDistance, float maxQueryDistance, const int16_t **end) const { + decimal minQueryDistance, decimal maxQueryDistance, const int16_t **end) const { // Instead of computing the angle for every pair in the database, we pre-compute the /cosines/ // of the min and max query distances so that we can compare against dot products directly! As - // angle increases, cosine decreases, up to M_PI (and queries larger than that don't really make + // angle increases, cosine decreases, up to DECIMAL_M_PI (and queries larger than that don't really make // sense anyway) - assert(maxQueryDistance <= M_PI); + assert(maxQueryDistance <= DECIMAL_M_PI); - float maxQueryCos = cos(minQueryDistance); - float minQueryCos = cos(maxQueryDistance); + decimal maxQueryCos = DECIMAL_COS(minQueryDistance); + decimal minQueryCos = DECIMAL_COS(maxQueryDistance); long liberalUpperIndex; long liberalLowerIndex = index.QueryLiberal(minQueryDistance, maxQueryDistance, &liberalUpperIndex); @@ -280,8 +284,8 @@ long PairDistanceKVectorDatabase::NumPairs() const { } /// Return the distances from the given star to each star it's paired with in the database (for debugging). -std::vector PairDistanceKVectorDatabase::StarDistances(int16_t star, const Catalog &catalog) const { - std::vector result; +std::vector PairDistanceKVectorDatabase::StarDistances(int16_t star, const Catalog &catalog) const { + std::vector result; for (int i = 0; i < NumPairs(); i++) { if (pairs[i*2] == star || pairs[i*2+1] == star) { result.push_back(AngleUnit(catalog[pairs[i*2]].spatial, catalog[pairs[i*2+1]].spatial)); @@ -296,6 +300,7 @@ std::vector PairDistanceKVectorDatabase::StarDistances(int16_t star, cons | size | name | description | |------+----------------+---------------------------------------------| | 4 | magicValue | unique database identifier | + | 4 | flags | [X, X, X, isDouble?] | | 4 | databaseLength | length in bytes (32-bit unsigned) | | n | database | the entire database. 8-byte aligned | | ... | ... | More databases (each has value, length, db) | @@ -318,6 +323,21 @@ const unsigned char *MultiDatabase::SubDatabasePointer(int32_t magicValue) const if (curMagicValue == 0) { return nullptr; } + uint32_t dbFlags = DeserializePrimitive(des); + + // Ensure that our database is using the same type as the runtime. + #ifdef LOST_FLOAT_MODE + if (!isFlagSet(dbFlags, MULTI_DB_FLOAT_FLAG)) { + std::cerr << "LOST was compiled in float mode. This database was serialized in double mode and is incompatible." << std::endl; + exit(1); + } + #else + if (isFlagSet(dbFlags, MULTI_DB_FLOAT_FLAG)) { + std::cerr << "LOST was compiled in double mode. This database was serialized in float mode and is incompatible." << std::endl; + exit(1); + } + #endif + uint32_t dbLength = DeserializePrimitive(des); assert(dbLength > 0); DeserializePadding(des); // align to an 8-byte boundary @@ -331,9 +351,11 @@ const unsigned char *MultiDatabase::SubDatabasePointer(int32_t magicValue) const } void SerializeMultiDatabase(SerializeContext *ser, - const MultiDatabaseDescriptor &dbs) { + const MultiDatabaseDescriptor &dbs, + uint32_t flags) { for (const MultiDatabaseEntry &multiDbEntry : dbs) { SerializePrimitive(ser, multiDbEntry.magicValue); + SerializePrimitive(ser, flags); SerializePrimitive(ser, multiDbEntry.bytes.size()); SerializePadding(ser); std::copy(multiDbEntry.bytes.cbegin(), multiDbEntry.bytes.cend(), std::back_inserter(ser->buffer)); diff --git a/src/databases.hpp b/src/databases.hpp index c8337c11..c510ab78 100644 --- a/src/databases.hpp +++ b/src/databases.hpp @@ -12,6 +12,8 @@ namespace lost { const int32_t kCatalogMagicValue = 0xF9A283BC; +inline bool isFlagSet(uint32_t dbFlags, uint32_t flag); + /** * A data structure enabling constant-time range queries into fixed numerical data. * @@ -22,27 +24,27 @@ class KVectorIndex { public: explicit KVectorIndex(DeserializeContext *des); - long QueryLiberal(float minQueryDistance, float maxQueryDistance, long *upperIndex) const; + long QueryLiberal(decimal minQueryDistance, decimal maxQueryDistance, long *upperIndex) const; /// The number of data points in the data referred to by the kvector long NumValues() const { return numValues; }; long NumBins() const { return numBins; }; /// Upper bound on elements - float Max() const { return max; }; + decimal Max() const { return max; }; // Lower bound on elements - float Min() const { return min; }; + decimal Min() const { return min; }; private: - long BinFor(float dist) const; + long BinFor(decimal dist) const; long numValues; - float min; - float max; - float binWidth; + decimal min; + decimal max; + decimal binWidth; long numBins; const int32_t *bins; }; -void SerializePairDistanceKVector(SerializeContext *, const Catalog &, float minDistance, float maxDistance, long numBins); +void SerializePairDistanceKVector(SerializeContext *, const Catalog &, decimal minDistance, decimal maxDistance, long numBins); /** * A database storing distances between pairs of stars. @@ -53,14 +55,14 @@ class PairDistanceKVectorDatabase { public: explicit PairDistanceKVectorDatabase(DeserializeContext *des); - const int16_t *FindPairsLiberal(float min, float max, const int16_t **end) const; - const int16_t *FindPairsExact(const Catalog &, float min, float max, const int16_t **end) const; - std::vector StarDistances(int16_t star, const Catalog &) const; + const int16_t *FindPairsLiberal(decimal min, decimal max, const int16_t **end) const; + const int16_t *FindPairsExact(const Catalog &, decimal min, decimal max, const int16_t **end) const; + std::vector StarDistances(int16_t star, const Catalog &) const; /// Upper bound on stored star pair distances - float MaxDistance() const { return index.Max(); }; + decimal MaxDistance() const { return index.Max(); }; /// Lower bound on stored star pair distances - float MinDistance() const { return index.Min(); }; + decimal MinDistance() const { return index.Min(); }; /// Exact number of stored pairs long NumPairs() const; @@ -85,7 +87,7 @@ class PairDistanceKVectorDatabase { // public: // explicit TripleInnerKVectorDatabase(const unsigned char *databaseBytes); -// void FindTriplesLiberal(float min, float max, long **begin, long **end) const; +// void FindTriplesLiberal(decimal min, decimal max, long **begin, long **end) const; // private: // KVectorIndex index; // int16_t *triples; @@ -96,6 +98,9 @@ class PairDistanceKVectorDatabase { * This is almost always the database that is actually passed to star-id algorithms in the real world, since you'll want to store at least the catalog plus one specific database. * Multi-databases are essentially a map from "magic values" to database buffers. */ + +#define MULTI_DB_FLOAT_FLAG 0x0001 // By default, our DB is in double mode. + class MultiDatabase { public: /// Create a multidatabase from a serialized multidatabase. @@ -111,12 +116,13 @@ class MultiDatabaseEntry { : magicValue(magicValue), bytes(bytes) { } int32_t magicValue; + uint32_t flags; std::vector bytes; }; typedef std::vector MultiDatabaseDescriptor; -void SerializeMultiDatabase(SerializeContext *, const MultiDatabaseDescriptor &dbs); +void SerializeMultiDatabase(SerializeContext *, const MultiDatabaseDescriptor &dbs, uint32_t flags); } diff --git a/src/decimal.hpp b/src/decimal.hpp new file mode 100644 index 00000000..643d72ef --- /dev/null +++ b/src/decimal.hpp @@ -0,0 +1,64 @@ +#ifndef DECIMAL_H +#define DECIMAL_H + +#include +#include + +#ifdef LOST_FLOAT_MODE + typedef float decimal; + #define STR_TO_DECIMAL(x) std::stof(x) +#else + typedef double decimal; + #define STR_TO_DECIMAL(x) std::stod(x) +#endif + +// This should only be used sparingly. +// It's better to verbosely typecast sometimes. Only use these to prevent promotions. +// The reason why this isn't used everywhere instead of the wrapped macros is +// because the code becomes hard to read when there are multiple layers of typecasting. +// With this method, we might have more preprocessing to do BUT the code remains readable +// as the methods remain relatively the same. +#define DECIMAL(x) ((decimal) x) + +// Math Constants wrapped with Decimal typecast +#define DECIMAL_M_E ((decimal) M_E) /* e */ +#define DECIMAL_M_LOG2E ((decimal) M_LOG2E) /* log_2 e */ +#define DECIMAL_M_LOG10E ((decimal) M_LOG10E) /* log_10 e */ +#define DECIMAL_M_LN2 ((decimal) M_LN2) /* log_e 2 */ +#define DECIMAL_M_LN10 ((decimal) M_LN10) /* log_e 10 */ +#define DECIMAL_M_PI ((decimal) M_PI) /* pi */ +#define DECIMAL_M_PI_2 ((decimal) M_PI_2) /* pi/2 */ +#define DECIMAL_M_PI_4 ((decimal) M_PI_4) /* pi/4 */ +#define DECIMAL_M_1_PI ((decimal) M_1_PI) /* 1/pi */ +#define DECIMAL_M_2_PI ((decimal) M_2_PI) /* 2/pi */ +#define DECIMAL_M_2_SQRTPI ((decimal) M_2_SQRTPI) /* 2/sqrt(pi) */ +#define DECIMAL_M_SQRT2 ((decimal) M_SQRT2) /* sqrt(2) */ +#define DECIMAL_M_SQRT1_2 ((decimal) M_SQRT1_2) /* 1/sqrt(2) */ + +// Math Functions wrapped with Decimal typecast +#define DECIMAL_POW(base,power) ((decimal) std::pow(base, power)) +#define DECIMAL_SQRT(x) ((decimal) std::sqrt(x)) +#define DECIMAL_LOG(x) ((decimal) std::log(x)) +#define DECIMAL_EXP(x) ((decimal) std::exp(x)) +#define DECIMAL_ERF(x) ((decimal) std::erf(x)) + +// Rouding methods wrapped with Decimal typecast) +#define DECIMAL_ROUND(x) ((decimal) std::round(x)) +#define DECIMAL_CEIL(x) ((decimal) std::ceil(x)) +#define DECIMAL_FLOOR(x) ((decimal) std::floor(x)) +#define DECIMAL_ABS(x) ((decimal) std::abs(x)) + +// Trig Methods wrapped with Decimal typecast) +#define DECIMAL_SIN(x) ((decimal) std::sin(x)) +#define DECIMAL_COS(x) ((decimal) std::cos(x)) +#define DECIMAL_TAN(x) ((decimal) std::tan(x)) +#define DECIMAL_ASIN(x) ((decimal) std::asin(x)) +#define DECIMAL_ACOS(x) ((decimal) std::acos(x)) +#define DECIMAL_ATAN(x) ((decimal) std::atan(x)) +#define DECIMAL_ATAN2(x,y) ((decimal) std::atan2(x,y)) + +// Float methods wrapped with Decimal typecast) +#define DECIMAL_FMA(x,y,z) ((decimal) std::fma(x,y,z)) +#define DECIMAL_HYPOT(x,y) ((decimal) std::hypot(x,y)) + +#endif // decimal.hpp diff --git a/src/io.cpp b/src/io.cpp index 32ef63a9..389f14cc 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -24,6 +24,7 @@ #include "attitude-estimators.hpp" #include "attitude-utils.hpp" #include "databases.hpp" +#include "decimal.hpp" #include "star-id.hpp" #include "star-utils.hpp" @@ -57,18 +58,25 @@ UserSpecifiedOutputStream::~UserSpecifiedOutputStream() { std::vector BscParse(std::string tsvPath) { std::vector result; FILE *file; - double raj2000, dej2000; + decimal raj2000, dej2000; int magnitudeHigh, magnitudeLow, name; char weird; file = fopen(tsvPath.c_str(), "r"); + if (file == NULL) { printf("Error opening file: %s\n", strerror(errno)); exit(1); // TODO: do we want any other error handling? return result; } - while (EOF != fscanf(file, "%lf|%lf|%d|%c|%d.%d", + #ifdef LOST_FLOAT_MODE + std::string format = "%f|%f|%d|%c|%d.%d"; + #else + std::string format = "%lf|%lf|%d|%c|%d.%d"; + #endif + + while (EOF != fscanf(file, format.c_str(), &raj2000, &dej2000, &name, &weird, &magnitudeHigh, &magnitudeLow)) { @@ -102,7 +110,7 @@ const Catalog &CatalogRead() { return a.spatial.x < b.spatial.x; }); for (int i = catalog.size(); i > 0; i--) { - if ((catalog[i].spatial - catalog[i-1].spatial).Magnitude() < 5e-5) { // 70 stars removed at this threshold. + if ((catalog[i].spatial - catalog[i-1].spatial).Magnitude() < DECIMAL(5e-5)) { // 70 stars removed at this threshold. if (catalog[i].magnitude > catalog[i-1].magnitude) { catalog.erase(catalog.begin() + i); } else { @@ -191,13 +199,13 @@ void SurfacePlot(std::string description, cairo_set_font_options(cairoCtx, cairoFontOptions); cairo_text_extents_t cairoTextExtents; cairo_text_extents(cairoCtx, "1234567890", &cairoTextExtents); - double textHeight = cairoTextExtents.height; + decimal textHeight = cairoTextExtents.height; for (const Star ¢roid : stars) { // plot the box around the star - if (centroid.radiusX > 0.0f) { - float radiusX = centroid.radiusX; - float radiusY = centroid.radiusY > 0.0f ? + if (centroid.radiusX > DECIMAL(0.0)) { + decimal radiusX = centroid.radiusX; + decimal radiusY = centroid.radiusY > DECIMAL(0.0) ? centroid.radiusY : radiusX; // Rectangles should be entirely /outside/ the radius of the star, so the star is @@ -210,8 +218,8 @@ void SurfacePlot(std::string description, cairo_stroke(cairoCtx); } else { cairo_rectangle(cairoCtx, - floor(centroid.position.x), - floor(centroid.position.y), + DECIMAL_FLOOR(centroid.position.x), + DECIMAL_FLOOR(centroid.position.y), 1, 1); cairo_fill(cairoCtx); } @@ -226,11 +234,11 @@ void SurfacePlot(std::string description, const Star ¢roid = stars[starId.starIndex]; cairo_move_to(cairoCtx, - centroid.radiusX > 0.0f + centroid.radiusX > DECIMAL(0.0) ? centroid.position.x + centroid.radiusX + 3 : centroid.position.x + 8, - centroid.radiusY > 0.0f + centroid.radiusY > DECIMAL(0.0) ? centroid.position.y - centroid.radiusY + textHeight : centroid.position.y + 10); @@ -266,7 +274,7 @@ typedef StarIdAlgorithm *(*StarIdAlgorithmFactory)(); typedef AttitudeEstimationAlgorithm *(*AttitudeEstimationAlgorithmFactory)(); SerializeContext serFromDbValues(const DatabaseOptions &values) { - return SerializeContext(values.swapIntegerEndianness, values.swapFloatEndianness); + return SerializeContext(values.swapIntegerEndianness, values.swapDecimalEndianness); } MultiDatabaseDescriptor GenerateDatabases(const Catalog &catalog, const DatabaseOptions &values) { @@ -278,8 +286,8 @@ MultiDatabaseDescriptor GenerateDatabases(const Catalog &catalog, const Database dbEntries.emplace_back(kCatalogMagicValue, catalogSer.buffer); if (values.kvector) { - float minDistance = DegToRad(values.kvectorMinDistance); - float maxDistance = DegToRad(values.kvectorMaxDistance); + decimal minDistance = DegToRad(values.kvectorMinDistance); + decimal maxDistance = DegToRad(values.kvectorMaxDistance); long numBins = values.kvectorNumDistanceBins; SerializeContext ser = serFromDbValues(values); SerializePairDistanceKVector(&ser, catalog, minDistance, maxDistance, numBins); @@ -308,7 +316,7 @@ std::ostream &operator<<(std::ostream &os, const Camera &camera) { * Calculate the focal length, in pixels, based on the given command line options. * This function exists because there are two ways to specify how "zoomed-in" the camera is. One way is using just FOV, which is useful when generating false images. Another is a combination of pixel size and focal length, which is useful for physical cameras. */ -float FocalLengthFromOptions(const PipelineOptions &values, int xResolution) { +decimal FocalLengthFromOptions(const PipelineOptions &values, int xResolution) { if ((values.pixelSize != -1) ^ (values.focalLength != 0)) { std::cerr << "ERROR: Exactly one of --pixel-size or --focal-length were set." << std::endl; exit(1); @@ -383,7 +391,7 @@ PipelineInputList GetPngPipelineInput(const PipelineOptions &values) { int xResolution = cairo_image_surface_get_width(cairoSurface); int yResolution = cairo_image_surface_get_height(cairoSurface); - float focalLengthPixels = FocalLengthFromOptions(values, xResolution); + decimal focalLengthPixels = FocalLengthFromOptions(values, xResolution); Camera cam = Camera(focalLengthPixels, xResolution, yResolution); result.push_back(std::unique_ptr(new PngPipelineInput(cairoSurface, cam, CatalogRead()))); @@ -406,11 +414,11 @@ PipelineInputList GetPngPipelineInput(const PipelineOptions &values) { /// A star used in simulated image generation. Contains extra data about how to simulate the star. class GeneratedStar : public Star { public: - GeneratedStar(Star star, float peakBrightness, Vec2 motionBlurDelta) + GeneratedStar(Star star, decimal peakBrightness, Vec2 motionBlurDelta) : Star(star), peakBrightness(peakBrightness), delta(motionBlurDelta) { }; /// the brightness density per time unit at the center of the star. 0.0 is black, 1.0 is white. - float peakBrightness; + decimal peakBrightness; /// (only meaningful with motion blur) Where the star will appear one time unit in the future. Vec2 delta; @@ -432,23 +440,23 @@ class GeneratedStar : public Star { * @param stddev The standard deviation of spread of the star. Higher values make stars more spread out. See command line documentation. * @return Indefinite integral of brightness density. */ -static float MotionBlurredPixelBrightness(const Vec2 &pixel, const GeneratedStar &generatedStar, - float t, float stddev) { +static decimal MotionBlurredPixelBrightness(const Vec2 &pixel, const GeneratedStar &generatedStar, + decimal t, decimal stddev) { const Vec2 &p0 = generatedStar.position; const Vec2 &delta = generatedStar.delta; const Vec2 d0 = p0 - pixel; return generatedStar.peakBrightness - * stddev*sqrt(M_PI) / (sqrt(2)*delta.Magnitude()) - * exp(pow(d0.x*delta.x + d0.y*delta.y, 2) / (2*stddev*stddev*delta.MagnitudeSq()) + * stddev*DECIMAL_SQRT(DECIMAL_M_PI) / (DECIMAL_SQRT(2)*delta.Magnitude()) + * DECIMAL_EXP(DECIMAL_POW(d0.x*delta.x + d0.y*delta.y, 2) / (2*stddev*stddev*delta.MagnitudeSq()) - d0.MagnitudeSq() / (2*stddev*stddev)) - * erf((t*delta.MagnitudeSq() + d0.x*delta.x + d0.y*delta.y) / (stddev*sqrt(2)*delta.Magnitude())); + * DECIMAL_ERF((t*delta.MagnitudeSq() + d0.x*delta.x + d0.y*delta.y) / (stddev*DECIMAL_SQRT(2)*delta.Magnitude())); } /// Like motionBlurredPixelBrightness, but for when motion blur is disabled. -static float StaticPixelBrightness(const Vec2 &pixel, const GeneratedStar &generatedStar, - float t, float stddev) { +static decimal StaticPixelBrightness(const Vec2 &pixel, const GeneratedStar &generatedStar, + decimal t, decimal stddev) { const Vec2 d0 = generatedStar.position - pixel; - return generatedStar.peakBrightness * t * exp(-d0.MagnitudeSq() / (2 * stddev * stddev)); + return generatedStar.peakBrightness * t * DECIMAL_EXP(-d0.MagnitudeSq() / (2 * stddev * stddev)); } /** @@ -460,12 +468,12 @@ static float StaticPixelBrightness(const Vec2 &pixel, const GeneratedStar &gener * deviation 1/5th of the cutoff brightness. We compute the probability that, taking read noise into * account, the observed energy would be less than the cutoff energy. */ -static float CentroidImagingProbability(float mag, float cutoffMag) { - float brightness = MagToBrightness(mag); - float cutoffBrightness = MagToBrightness(cutoffMag); - float stddev = cutoffBrightness/5.0; +static decimal CentroidImagingProbability(decimal mag, decimal cutoffMag) { + decimal brightness = MagToBrightness(mag); + decimal cutoffBrightness = MagToBrightness(cutoffMag); + decimal stddev = cutoffBrightness/DECIMAL(5.0); // CDF of Normal distribution with given mean and stddev - return 1 - (0.5 * (1 + erf((cutoffBrightness-brightness)/(stddev*sqrt(2.0))))); + return 1 - (DECIMAL(0.5) * (1 + DECIMAL_ERF((cutoffBrightness-brightness)/(stddev*DECIMAL_SQRT(2.0))))); } const int kMaxBrightness = 255; @@ -480,38 +488,38 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, std::default_random_engine *rng, bool centroidsOnly, - float zeroMagTotalPhotons, - float starSpreadStdDev, - float saturationPhotons, - float darkCurrent, - float readNoiseStdDev, + decimal zeroMagTotalPhotons, + decimal starSpreadStdDev, + decimal saturationPhotons, + decimal darkCurrent, + decimal readNoiseStdDev, Attitude motionBlurDirection, // applied on top of the attitude - float exposureTime, - float readoutTime, // zero for no rolling shutter + decimal exposureTime, + decimal readoutTime, // zero for no rolling shutter bool shotNoise, int oversampling, int numFalseStars, int falseStarMinMagnitude, int falseStarMaxMagnitude, int cutoffMag, - float perturbationStddev) + decimal perturbationStddev) : camera(camera), attitude(attitude), catalog(catalog) { assert(falseStarMaxMagnitude <= falseStarMinMagnitude); - assert(perturbationStddev >= 0.0); + assert(perturbationStddev >= DECIMAL(0.0)); image.width = camera.XResolution(); image.height = camera.YResolution(); // number of true photons each pixel receives. assert(oversampling >= 1); - int oversamplingPerAxis = ceil(sqrt(oversampling)); + int oversamplingPerAxis = DECIMAL_CEIL(DECIMAL_SQRT(oversampling)); if (oversamplingPerAxis*oversamplingPerAxis != oversampling) { std::cerr << "WARNING: oversampling was not a perfect square. Rounding up to " << oversamplingPerAxis*oversamplingPerAxis << "." << std::endl; } assert(exposureTime > 0); - bool motionBlurEnabled = abs(motionBlurDirection.GetQuaternion().Angle()) > 0.001; + bool motionBlurEnabled = abs(motionBlurDirection.GetQuaternion().Angle()) > DECIMAL(0.001); Quaternion motionBlurDirectionQ = motionBlurDirection.GetQuaternion(); // attitude at the middle of exposure time Quaternion currentAttitude = attitude.GetQuaternion(); @@ -521,20 +529,20 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, // a star with 1 photon has peak density 1/(2pi sigma^2), because 2d gaussian formula. Then just // multiply up proportionally! - float zeroMagPeakPhotonDensity = zeroMagTotalPhotons / (2*M_PI * starSpreadStdDev*starSpreadStdDev); + decimal zeroMagPeakPhotonDensity = zeroMagTotalPhotons / (2*DECIMAL_M_PI * starSpreadStdDev*starSpreadStdDev); // TODO: Is it 100% correct to just copy the standard deviation in both dimensions? - std::normal_distribution perturbation1DDistribution(0.0, perturbationStddev); + std::normal_distribution perturbation1DDistribution(DECIMAL(0.0), perturbationStddev); Catalog catalogWithFalse = catalog; - std::uniform_real_distribution uniformDistribution(0.0, 1.0); + std::uniform_real_distribution uniformDistribution(DECIMAL(0.0), DECIMAL(1.0)); std::uniform_int_distribution magnitudeDistribution(falseStarMaxMagnitude, falseStarMinMagnitude); for (int i = 0; i < numFalseStars; i++) { - float ra = uniformDistribution(*rng) * 2*M_PI; + decimal ra = uniformDistribution(*rng) * 2*DECIMAL_M_PI; // to be uniform around sphere. Borel-Kolmogorov paradox is calling - float de = asin(uniformDistribution(*rng)*2 - 1); - float magnitude = magnitudeDistribution(*rng); + decimal de = DECIMAL_ASIN(uniformDistribution(*rng)*2 - 1); + decimal magnitude = magnitudeDistribution(*rng); catalogWithFalse.push_back(CatalogStar(ra, de, magnitude, -1)); } @@ -553,15 +561,15 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, Vec3 futureSpatial = futureAttitude.Rotate(catalogWithFalse[i].spatial); Vec2 delta = camera.SpatialToCamera(futureSpatial) - camCoords; if (!motionBlurEnabled) { - delta = {0, 0}; // avoid floating point funny business + delta = {0, 0}; // avoid decimaling point funny business } // radiant intensity, in photons per time unit per pixel, at the center of the star. - float peakBrightnessPerTime = zeroMagPeakPhotonDensity * MagToBrightness(catalogStar.magnitude); - float interestingThreshold = 0.05; // we don't need to check pixels that are expected to + decimal peakBrightnessPerTime = zeroMagPeakPhotonDensity * MagToBrightness(catalogStar.magnitude); + decimal interestingThreshold = DECIMAL(0.05); // we don't need to check pixels that are expected to // receive this many photons or fewer. // inverse of the function defining the Gaussian distribution: Find out how far from the // mean we'll have to go until the number of photons is less than interestingThreshold - float radius = ceil(sqrt(-log(interestingThreshold/peakBrightnessPerTime/exposureTime)*2*M_PI*starSpreadStdDev*starSpreadStdDev)); + decimal radius = DECIMAL_CEIL(DECIMAL_SQRT(-DECIMAL_LOG(interestingThreshold/peakBrightnessPerTime/exposureTime)*2*DECIMAL_M_PI*starSpreadStdDev*starSpreadStdDev)); Star star = Star(camCoords.x, camCoords.y, radius, radius, // important to invert magnitude here, so that centroid magnitude becomes larger for brighter stars. @@ -584,7 +592,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, // for input, though, add perturbation and stuff. Star inputStar = star; - if (perturbationStddev > 0.0) { + if (perturbationStddev > DECIMAL(0.0)) { // clamp to within 2 standard deviations for some reason: inputStar.position.x += std::max(std::min(perturbation1DDistribution(*rng), 2*perturbationStddev), -2*perturbationStddev); inputStar.position.y += std::max(std::min(perturbation1DDistribution(*rng), 2*perturbationStddev), -2*perturbationStddev); @@ -608,12 +616,12 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, return; } - std::vector photonsBuffer(image.width*image.height, 0); + std::vector photonsBuffer(image.width*image.height, 0); for (const GeneratedStar &star : generatedStars) { // delta will be exactly (0,0) when motion blur disabled - Vec2 earliestPosition = star.position - star.delta*(exposureTime/2.0 + readoutTime/2.0); - Vec2 latestPosition = star.position + star.delta*(exposureTime/2.0 + readoutTime/2.0); + Vec2 earliestPosition = star.position - star.delta*(exposureTime/DECIMAL(2.0) + readoutTime/DECIMAL(2.0)); + Vec2 latestPosition = star.position + star.delta*(exposureTime/DECIMAL(2.0) + readoutTime/DECIMAL(2.0)); int xMin = std::max(0, (int)std::min(earliestPosition.x - star.radiusX, latestPosition.x - star.radiusX)); int xMax = std::min(image.width-1, (int)std::max(earliestPosition.x + star.radiusX, latestPosition.x + star.radiusX)); int yMin = std::max(0, (int)std::min(earliestPosition.y - star.radiusX, latestPosition.y - star.radiusX)); @@ -621,7 +629,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, // peak brightness is measured in photons per time unit per pixel, so if oversampling, we // need to convert units to photons per time unit per sample - float oversamplingBrightnessFactor = oversamplingPerAxis*oversamplingPerAxis; + decimal oversamplingBrightnessFactor = oversamplingPerAxis*oversamplingPerAxis; // the star.x and star.y refer to the pixel whose top left corner the star should appear at // (and fractional amounts are relative to the corner). When we color a pixel, we ideally @@ -631,17 +639,17 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, for (int yPixel = yMin; yPixel <= yMax; yPixel++) { // offset of beginning & end of readout compared to beginning & end of readout for // center row - float readoutOffset = readoutTime * (yPixel - image.height/2.0) / image.height; - float tStart = -exposureTime/2.0 + readoutOffset; - float tEnd = exposureTime/2.0 + readoutOffset; + decimal readoutOffset = readoutTime * (yPixel - image.height/DECIMAL(2.0)) / image.height; + decimal tStart = -exposureTime/DECIMAL(2.0) + readoutOffset; + decimal tEnd = exposureTime/DECIMAL(2.0) + readoutOffset; // loop through all samples in the current pixel for (int xSample = 0; xSample < oversamplingPerAxis; xSample++) { for (int ySample = 0; ySample < oversamplingPerAxis; ySample++) { - float x = xPixel + (xSample+0.5)/oversamplingPerAxis; - float y = yPixel + (ySample+0.5)/oversamplingPerAxis; + decimal x = xPixel + (xSample+DECIMAL(0.5))/oversamplingPerAxis; + decimal y = yPixel + (ySample+DECIMAL(0.5))/oversamplingPerAxis; - float curPhotons; + decimal curPhotons; if (motionBlurEnabled) { curPhotons = (MotionBlurredPixelBrightness({x, y}, star, tEnd, starSpreadStdDev) @@ -652,7 +660,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, / oversamplingBrightnessFactor; } - assert(0.0 <= curPhotons); + assert(DECIMAL(0.0) <= curPhotons); photonsBuffer[xPixel + yPixel*image.width] += curPhotons; } @@ -661,13 +669,13 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, } } - std::normal_distribution readNoiseDist(0.0, readNoiseStdDev); + std::normal_distribution readNoiseDist(DECIMAL(0.0), readNoiseStdDev); // convert from photon counts to observed pixel brightnesses, applying noise and such. imageData = std::vector(image.width*image.height); image.image = imageData.data(); for (int i = 0; i < image.width * image.height; i++) { - float curBrightness = 0; + decimal curBrightness = 0; // dark current (Constant) curBrightness += darkCurrent; @@ -682,8 +690,8 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, // range. This is problematic if the mean is far above the max long value, because then it // might have to sample many many times (and furthermore, the results won't be useful // anyway) - float photons = photonsBuffer[i]; - if (photons > (float)LONG_MAX - 3.0*sqrt(LONG_MAX)) { + decimal photons = photonsBuffer[i]; + if (photons > DECIMAL(LONG_MAX) - DECIMAL(3.0) * DECIMAL_SQRT(LONG_MAX)) { std::cout << "ERROR: One of the pixels had too many photons. Generated image would not be physically accurate, exiting." << std::endl; exit(1); } @@ -695,7 +703,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, curBrightness += quantizedPhotons / saturationPhotons; // std::clamp not introduced until C++17, so we avoid it. - float clampedBrightness = std::max(std::min(curBrightness, (float)1.0), (float)0.0); + decimal clampedBrightness = std::max(std::min(curBrightness, DECIMAL(1.0)), DECIMAL(0.0)); imageData[i] = floor(clampedBrightness * kMaxBrightness); // TODO: off-by-one, 256? } } @@ -705,7 +713,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, * Takes a random engine as a parameter. */ static Attitude RandomAttitude(std::default_random_engine* pReng) { - std::uniform_real_distribution randomAngleDistribution(0, 1); + std::uniform_real_distribution randomAngleDistribution(0, 1); // normally the ranges of the Ra and Dec are: // Dec: [-90 deg, 90 deg] --> [-pi/2 rad, pi/2 rad], where negative means south @@ -713,9 +721,9 @@ static Attitude RandomAttitude(std::default_random_engine* pReng) { // Ra: [0 deg, 360 deg] --> [0 rad, 2pi rad ] // Roll: [0 rad, 2 pi rad] - float randomRa = 2 * M_PI * randomAngleDistribution(*pReng); - float randomDec = (M_PI / 2) - acos(1 - 2 * randomAngleDistribution(*pReng)); //acos returns a float in range [0, pi] - float randomRoll = 2 * M_PI * randomAngleDistribution(*pReng); + decimal randomRa = 2 * DECIMAL_M_PI * randomAngleDistribution(*pReng); + decimal randomDec = (DECIMAL_M_PI / 2) - acos(1 - 2 * randomAngleDistribution(*pReng)); //acos returns a decimal in range [0, pi] + decimal randomRoll = 2 * DECIMAL_M_PI * randomAngleDistribution(*pReng); Attitude randAttitude = Attitude(SphericalToQuaternion(randomRa, randomDec, randomRoll)); @@ -748,7 +756,7 @@ PipelineInputList GetGeneratedPipelineInput(const PipelineOptions &values) { DegToRad(values.generateBlurRoll))); PipelineInputList result; - float focalLength = FocalLengthFromOptions(values, values.generateXRes); + decimal focalLength = FocalLengthFromOptions(values, values.generateXRes); for (int i = 0; i < values.generate; i++) { @@ -837,6 +845,7 @@ Pipeline SetPipeline(const PipelineOptions &values) { // TODO: more flexible or sth // TODO: don't allow setting star-id until database is set, and perhaps limit the star-id // choices to those compatible with the database? + // // centroid algorithm stage if (values.centroidAlgo == "dummy") { @@ -1025,18 +1034,18 @@ class CentroidComparison { * Average distance from actual to expected centroids (in pixels) * Only correct centroids are considered in this average. */ - float meanError; + decimal meanError; /** * Number of actual stars within the centroiding threshold of an expected star. */ - float numCorrectCentroids; + decimal numCorrectCentroids; /** * Stars in actual but not expected. Ideally 0 - * This is a float because we may average multiple centroid comparisons together. + * This is a decimal because we may average multiple centroid comparisons together. */ - float numExtraCentroids; + decimal numExtraCentroids; // We no longer have a num missing because often generated stars have too low of a signal-to-noise ratio and the centroid algo won't pick them up. }; @@ -1045,21 +1054,21 @@ class CentroidComparison { /// `two` whose distance to the current `one` star is <=threshold. In a (ordered) multimap, the /// insertion order is preserved for elements with the same key, and indeed we'll sort the elements /// corresponding to each key by distance from the corresponding `one` star. -static std::multimap FindClosestCentroids(float threshold, +static std::multimap FindClosestCentroids(decimal threshold, const Stars &one, const Stars &two) { std::multimap result; for (int i = 0; i < (int)one.size(); i++) { - std::vector> closest; + std::vector> closest; for (int k = 0; k < (int)two.size(); k++) { - float currDistance = (one[i].position - two[k].position).Magnitude(); + decimal currDistance = (one[i].position - two[k].position).Magnitude(); if (currDistance <= threshold) { closest.emplace_back(currDistance, k); } } std::sort(closest.begin(), closest.end()); - for (const std::pair &pair : closest) { + for (const std::pair &pair : closest) { result.emplace(i, pair.second); } } @@ -1072,7 +1081,7 @@ static std::multimap FindClosestCentroids(float threshold, * Useful for debugging and benchmarking. * @param threshold The maximum number of pixels apart two centroids can be to be considered the same. */ -CentroidComparison CentroidsCompare(float threshold, +CentroidComparison CentroidsCompare(decimal threshold, const Stars &expected, const Stars &actual) { @@ -1120,7 +1129,7 @@ CentroidComparison CentroidComparisonsCombine(std::vector co StarIdComparison StarIdsCompare(const StarIdentifiers &expected, const StarIdentifiers &actual, // use these to map indices to names for the respective lists of StarIdentifiers const Catalog &expectedCatalog, const Catalog &actualCatalog, - float centroidThreshold, + decimal centroidThreshold, const Stars &expectedStars, const Stars &inputStars) { StarIdComparison result = { @@ -1266,7 +1275,7 @@ static void PipelineComparatorCentroids(std::ostream &os, const PipelineOptions &values) { int size = (int)expected.size(); - float threshold = values.centroidCompareThreshold; + decimal threshold = values.centroidCompareThreshold; std::vector comparisons; for (int i = 0; i < size; i++) { @@ -1288,7 +1297,7 @@ static void PrintCentroids(const std::string &prefix, // May be NULL. Should be the only the first starId, because we don't have any reasonable aggregative action to perform. const StarIdentifiers *starIds) { assert(starses.size() > 0); - float avgNumStars = 0; + decimal avgNumStars = 0; for (const Stars &stars : starses) { avgNumStars += stars.size(); } @@ -1527,9 +1536,9 @@ static void PipelineComparatorAttitude(std::ostream &os, // TODO: use Wahba loss function (maybe average per star) instead of just angle. Also break // apart roll error from boresight error. This is just quick and dirty for testing - float angleThreshold = DegToRad(values.attitudeCompareThreshold); + decimal angleThreshold = DegToRad(values.attitudeCompareThreshold); - float attitudeErrorSum = 0.0f; + decimal attitudeErrorSum = 0.0f; int numCorrect = 0; int numIncorrect = 0; @@ -1537,7 +1546,7 @@ static void PipelineComparatorAttitude(std::ostream &os, if (actual[i].attitude->IsKnown()) { Quaternion expectedQuaternion = expected[i]->ExpectedAttitude()->GetQuaternion(); Quaternion actualQuaternion = actual[i].attitude->GetQuaternion(); - float attitudeError = (expectedQuaternion * actualQuaternion.Conjugate()).SmallestAngle(); + decimal attitudeError = (expectedQuaternion * actualQuaternion.Conjugate()).SmallestAngle(); assert(attitudeError >= 0); if (attitudeError <= angleThreshold) { @@ -1549,9 +1558,9 @@ static void PipelineComparatorAttitude(std::ostream &os, } } - float attitudeErrorMean = attitudeErrorSum / numCorrect; - float fractionCorrect = (float)numCorrect / expected.size(); - float fractionIncorrect = (float)numIncorrect / expected.size(); + decimal attitudeErrorMean = DECIMAL(attitudeErrorSum) / numCorrect; + decimal fractionCorrect = DECIMAL(numCorrect) / expected.size(); + decimal fractionIncorrect = DECIMAL(numIncorrect) / expected.size(); os << "attitude_error_mean " << attitudeErrorMean << std::endl; os << "attitude_availability " << fractionCorrect << std::endl; @@ -1830,17 +1839,17 @@ void PipelineComparison(const PipelineInputList &expected, // void InspectFindStar(const Catalog &catalog) { // std::string raStr = PromptLine("Right Ascension"); -// float raRadians; +// decimal raRadians; // int raHours, raMinutes; -// float raSeconds; +// decimal raSeconds; // int raFormatTime = sscanf(raStr.c_str(), "%dh %dm %fs", &raHours, &raMinutes, &raSeconds); -// float raDeg; +// decimal raDeg; // int raFormatDeg = sscanf(raStr.c_str(), "%f", &raDeg); // if (raFormatTime == 3) { -// raRadians = (raHours * 2*M_PI/24) + (raMinutes * 2*M_PI/24/60) + (raSeconds * 2*M_PI/24/60/60); +// raRadians = (raHours * 2*DECIMAL_M_PI/24) + (raMinutes * 2*DECIMAL_M_PI/24/60) + (raSeconds * 2*DECIMAL_M_PI/24/60/60); // } else if (raFormatDeg == 1) { // raRadians = DegToRad(raFormatDeg); // } else { @@ -1850,18 +1859,18 @@ void PipelineComparison(const PipelineInputList &expected, // std::string deStr = PromptLine("Declination"); -// float deRadians; +// decimal deRadians; // int deDegPart, deMinPart; -// float deSecPart; +// decimal deSecPart; // char dummy[8]; // int deFormatParts = sscanf(deStr.c_str(), "%d%s %d%s %f%s", &deDegPart, dummy, &deMinPart, dummy, &deSecPart, dummy); -// float deDeg; +// decimal deDeg; // int deFormatDeg = sscanf(deStr.c_str(), "%f", &deDeg); // if (deFormatParts == 6) { -// deRadians = DegToRad(deDegPart + (float)deMinPart/60 + (float)deSecPart/60/60); +// deRadians = DegToRad(deDegPart + (decimal)deMinPart/60 + (decimal)deSecPart/60/60); // } else if (deFormatDeg == 1) { // deRadians = DegToRad(deFormatDeg); // } else { @@ -1871,7 +1880,7 @@ void PipelineComparison(const PipelineInputList &expected, // // find the star -// float tolerance = 0.001; +// decimal tolerance = 0.001; // Vec3 userSpatial = SphericalToSpatial(raRadians, deRadians); // int i = 0; // for (const CatalogStar &curStar : catalog) { @@ -1888,7 +1897,7 @@ void PipelineComparison(const PipelineInputList &expected, // void InspectPrintStar(const Catalog &catalog) { // auto stars = PromptCatalogStars(catalog, 1); -// float ra, de; +// decimal ra, de; // SpatialToSpherical(stars[0]->spatial, &ra, &de); // std::cout << "star_ra " << RadToDeg(ra) << std::endl; diff --git a/src/io.hpp b/src/io.hpp index 0ce9dff6..c950541e 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -76,6 +76,7 @@ class Image { // PIPELINE INPUT // //////////////////// + /// The command line options passed when running a pipeline class PipelineOptions { public: @@ -126,13 +127,13 @@ class GeneratedPipelineInput : public PipelineInput { GeneratedPipelineInput(const Catalog &, Attitude, Camera, std::default_random_engine *, bool centroidsOnly, - float observedReferenceBrightness, float starSpreadStdDev, - float sensitivity, float darkCurrent, float readNoiseStdDev, - Attitude motionBlurDirection, float exposureTime, float readoutTime, + decimal observedReferenceBrightness, decimal starSpreadStdDev, + decimal sensitivity, decimal darkCurrent, decimal readNoiseStdDev, + Attitude motionBlurDirection, decimal exposureTime, decimal readoutTime, bool shotNoise, int oversampling, int numFalseStars, int falseMinMagnitude, int falseMaxMagnitude, int cutoffMag, - float perturbationStddev); + decimal perturbationStddev); const Image *InputImage() const override { return ℑ }; @@ -274,7 +275,7 @@ void PipelineComparison(const PipelineInputList &expected, StarIdComparison StarIdsCompare(const StarIdentifiers &expected, const StarIdentifiers &actual, // use these to map indices to names for the respective lists of StarIdentifiers const Catalog &expectedCatalog, const Catalog &actualCatalog, - float centroidThreshold, + decimal centroidThreshold, const Stars &expectedStars, const Stars &inputStars); //////////////// diff --git a/src/main.cpp b/src/main.cpp index 5d26fec4..86c20b77 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,9 +5,11 @@ */ #include +#include #include #include +#include #include #include #include @@ -17,6 +19,7 @@ #include "databases.hpp" #include "centroiders.hpp" +#include "decimal.hpp" #include "io.hpp" #include "man-database.h" #include "man-pipeline.h" @@ -30,9 +33,16 @@ static void DatabaseBuild(const DatabaseOptions &values) { MultiDatabaseDescriptor dbEntries = GenerateDatabases(narrowedCatalog, values); SerializeContext ser = serFromDbValues(values); - SerializeMultiDatabase(&ser, dbEntries); + + // Create & Set Flags. + uint32_t dbFlags = 0; + dbFlags |= typeid(decimal) == typeid(float) ? MULTI_DB_FLOAT_FLAG : 0; + + // Serialize Flags + SerializeMultiDatabase(&ser, dbEntries, dbFlags); std::cerr << "Generated database with " << ser.buffer.size() << " bytes" << std::endl; + std::cerr << "Database flagged with " << std::bitset<8*sizeof(dbFlags)>(dbFlags) << std::endl; UserSpecifiedOutputStream pos = UserSpecifiedOutputStream(values.outputPath, true); pos.Stream().write((char *) ser.buffer.data(), ser.buffer.size()); diff --git a/src/pipeline-options.hpp b/src/pipeline-options.hpp index df39ce54..fc364e4f 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -14,67 +14,69 @@ #include +#include "decimal.hpp" + // CAMERA -LOST_CLI_OPTION("png" , std::string, png , "" , optarg , kNoDefaultArgument) -LOST_CLI_OPTION("focal-length" , float , focalLength , 0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("pixel-size" , float , pixelSize , -1 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("fov" , float , fov , 20 , atof(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("png" , std::string , png , "" , optarg , kNoDefaultArgument) +LOST_CLI_OPTION("focal-length" , decimal , focalLength , 0 , atof(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("pixel-size" , decimal , pixelSize , -1 , atof(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("fov" , decimal , fov , 20 , atof(optarg) , kNoDefaultArgument) // PIPELINE STAGES -LOST_CLI_OPTION("centroid-algo" , std::string, centroidAlgo , "" , optarg , "cog") -LOST_CLI_OPTION("centroid-dummy-stars" , int , centroidDummyNumStars , 5 , atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("centroid-mag-filter" , float , centroidMagFilter , -1 , atof(optarg) , 5) -LOST_CLI_OPTION("centroid-filter-brightest", int , centroidFilterBrightest, -1 , atoi(optarg) , 10) -LOST_CLI_OPTION("database" , std::string, databasePath , "" , optarg , kNoDefaultArgument) -LOST_CLI_OPTION("star-id-algo" , std::string, idAlgo , "" , optarg , "pyramid") -LOST_CLI_OPTION("angular-tolerance" , float , angularTolerance , .04 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("false-stars-estimate" , int , estimatedNumFalseStars , 500 , atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("max-mismatch-probability" , float , maxMismatchProb , .001, atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("attitude-algo" , std::string, attitudeAlgo , "" , optarg , "dqm") +LOST_CLI_OPTION("centroid-algo" , std::string, centroidAlgo , "" , optarg , "cog") +LOST_CLI_OPTION("centroid-dummy-stars" , int , centroidDummyNumStars , 5 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("centroid-mag-filter" , decimal , centroidMagFilter , -1 , STR_TO_DECIMAL(optarg) , 5) +LOST_CLI_OPTION("centroid-filter-brightest", int , centroidFilterBrightest , -1 , atoi(optarg) , 10) +LOST_CLI_OPTION("database" , std::string, databasePath , "" , optarg , kNoDefaultArgument) +LOST_CLI_OPTION("star-id-algo" , std::string, idAlgo , "" , optarg , "pyramid") +LOST_CLI_OPTION("angular-tolerance" , decimal , angularTolerance , .04 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("false-stars-estimate" , int , estimatedNumFalseStars , 500 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("max-mismatch-probability" , decimal , maxMismatchProb , .001, STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("attitude-algo" , std::string, attitudeAlgo , "" , optarg , "dqm") // OUTPUT COMPARISON -LOST_CLI_OPTION("centroid-compare-threshold", float , centroidCompareThreshold, 2 , atof(optarg), kNoDefaultArgument) -LOST_CLI_OPTION("attitude-compare-threshold", float , attitudeCompareThreshold, 1 , atof(optarg), kNoDefaultArgument) -LOST_CLI_OPTION("plot-raw-input" , std::string, plotRawInput , "", optarg , "-") -LOST_CLI_OPTION("plot-input" , std::string, plotInput , "", optarg , "-") -LOST_CLI_OPTION("plot-expected" , std::string, plotExpected , "", optarg , "-") -LOST_CLI_OPTION("plot-centroid-indices" , std::string, plotCentroidIndices , "", optarg , "-") -LOST_CLI_OPTION("plot-output" , std::string, plotOutput , "", optarg , "-") -LOST_CLI_OPTION("print-expected-centroids" , std::string, printExpectedCentroids , "", optarg , "-") -LOST_CLI_OPTION("print-input-centroids" , std::string, printInputCentroids , "", optarg , "-") -LOST_CLI_OPTION("print-actual-centroids" , std::string, printActualCentroids , "", optarg , "-") -LOST_CLI_OPTION("print-attitude" , std::string, printAttitude , "", optarg , "-") -LOST_CLI_OPTION("print-expected-attitude" , std::string, printExpectedAttitude , "", optarg , "-") -LOST_CLI_OPTION("print-speed" , std::string, printSpeed , "", optarg , "-") -LOST_CLI_OPTION("compare-centroids" , std::string, compareCentroids , "", optarg , "-") -LOST_CLI_OPTION("compare-star-ids" , std::string, compareStarIds , "", optarg , "-") -LOST_CLI_OPTION("compare-attitudes" , std::string, compareAttitudes , "", optarg , "-") +LOST_CLI_OPTION("centroid-compare-threshold", decimal , centroidCompareThreshold, 2 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("attitude-compare-threshold", decimal , attitudeCompareThreshold, 1 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("plot-raw-input" , std::string, plotRawInput , "", optarg , "-") +LOST_CLI_OPTION("plot-input" , std::string, plotInput , "", optarg , "-") +LOST_CLI_OPTION("plot-expected" , std::string, plotExpected , "", optarg , "-") +LOST_CLI_OPTION("plot-centroid-indices" , std::string, plotCentroidIndices , "", optarg , "-") +LOST_CLI_OPTION("plot-output" , std::string, plotOutput , "", optarg , "-") +LOST_CLI_OPTION("print-expected-centroids" , std::string, printExpectedCentroids , "", optarg , "-") +LOST_CLI_OPTION("print-input-centroids" , std::string, printInputCentroids , "", optarg , "-") +LOST_CLI_OPTION("print-actual-centroids" , std::string, printActualCentroids , "", optarg , "-") +LOST_CLI_OPTION("print-attitude" , std::string, printAttitude , "", optarg , "-") +LOST_CLI_OPTION("print-expected-attitude" , std::string, printExpectedAttitude , "", optarg , "-") +LOST_CLI_OPTION("print-speed" , std::string, printSpeed , "", optarg , "-") +LOST_CLI_OPTION("compare-centroids" , std::string, compareCentroids , "", optarg , "-") +LOST_CLI_OPTION("compare-star-ids" , std::string, compareStarIds , "", optarg , "-") +LOST_CLI_OPTION("compare-attitudes" , std::string, compareAttitudes , "", optarg , "-") // IMAGE GENERATION -LOST_CLI_OPTION("generate" , int , generate , 0 , atoi(optarg) , 1) -LOST_CLI_OPTION("generate-x-resolution" , int , generateXRes , 1024 , atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-y-resolution" , int , generateYRes , 1024 , atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-centroids-only" , bool , generateCentroidsOnly , false , atobool(optarg) , true) -LOST_CLI_OPTION("generate-zero-mag-photons" , float , generateZeroMagPhotons , 20000 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-saturation-photons" , float , generateSaturationPhotons , 150 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-spread-stddev" , float , generateSpreadStdDev , 1 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-shot-noise" , bool , generateShotNoise , true , atobool(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-dark-current" , float , generateDarkCurrent , 0.1 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-read-noise-stddev" , float , generateReadNoiseStdDev , .05 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-ra" , float , generateRa , 88 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-de" , float , generateDe , 7 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-roll" , float , generateRoll , 0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-random-attitudes" , bool , generateRandomAttitudes , false , atobool(optarg) , true) -LOST_CLI_OPTION("generate-blur-ra" , float , generateBlurRa , 0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-blur-de" , float , generateBlurDe , 0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-blur-roll" , float , generateBlurRoll , 0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-exposure" , float , generateExposure , 0.2 , atof(optarg) , 0.1) -LOST_CLI_OPTION("generate-readout-time" , float , generateReadoutTime , 0 , atof(optarg) , 0.01) -LOST_CLI_OPTION("generate-oversampling" , int , generateOversampling , 4 , atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-false-stars" , int , generateNumFalseStars , 0 , atoi(optarg) , 50) -LOST_CLI_OPTION("generate-false-min-mag" , float , generateFalseMinMag , 8 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-false-max-mag" , float , generateFalseMaxMag , 1 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-perturb-centroids" , float , generatePerturbationStddev, 0 , atof(optarg) , 0.2) -LOST_CLI_OPTION("generate-cutoff-mag" , float , generateCutoffMag , 6.0 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-seed" , int , generateSeed , 394859, atoi(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-time-based-seed" , bool , timeSeed , false , atobool(optarg) , true) +LOST_CLI_OPTION("generate" , int , generate , 0 , atoi(optarg) , 1) +LOST_CLI_OPTION("generate-x-resolution" , int , generateXRes , 1024 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-y-resolution" , int , generateYRes , 1024 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-centroids-only" , bool , generateCentroidsOnly , false , atobool(optarg) , true) +LOST_CLI_OPTION("generate-zero-mag-photons" , decimal , generateZeroMagPhotons , 20000 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-saturation-photons" , decimal , generateSaturationPhotons , 150 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-spread-stddev" , decimal , generateSpreadStdDev , 1 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-shot-noise" , bool , generateShotNoise , true , atobool(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-dark-current" , decimal , generateDarkCurrent , 0.1 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-read-noise-stddev" , decimal , generateReadNoiseStdDev , .05 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-ra" , decimal , generateRa , 88 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-de" , decimal , generateDe , 7 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-roll" , decimal , generateRoll , 0 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-random-attitudes" , bool , generateRandomAttitudes , false , atobool(optarg) , true) +LOST_CLI_OPTION("generate-blur-ra" , decimal , generateBlurRa , 0 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-blur-de" , decimal , generateBlurDe , 0 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-blur-roll" , decimal , generateBlurRoll , 0 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-exposure" , decimal , generateExposure , 0.2 , STR_TO_DECIMAL(optarg) , 0.1) +LOST_CLI_OPTION("generate-readout-time" , decimal , generateReadoutTime , 0 , STR_TO_DECIMAL(optarg) , 0.01) +LOST_CLI_OPTION("generate-oversampling" , int , generateOversampling , 4 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-false-stars" , int , generateNumFalseStars , 0 , atoi(optarg) , 50) +LOST_CLI_OPTION("generate-false-min-mag" , decimal , generateFalseMinMag , 8 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-false-max-mag" , decimal , generateFalseMaxMag , 1 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-perturb-centroids" , decimal , generatePerturbationStddev, 0 , STR_TO_DECIMAL(optarg) , 0.2) +LOST_CLI_OPTION("generate-cutoff-mag" , decimal , generateCutoffMag , 6.0 , STR_TO_DECIMAL(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-seed" , int , generateSeed , 394859, atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("generate-time-based-seed" , bool , timeSeed , false , atobool(optarg) , true) diff --git a/src/serialize-helpers.hpp b/src/serialize-helpers.hpp index 13a0a3fb..2e57fb87 100644 --- a/src/serialize-helpers.hpp +++ b/src/serialize-helpers.hpp @@ -1,7 +1,7 @@ /** * Helpers to serialize and deserialize arbitrary data types to disk * - * The serialization and deserialization helpers here assume that (a) integers and floating point + * The serialization and deserialization helpers here assume that (a) integers and decimal point * numbers are stored in the same format on the source and target systems except for endianness, and * (b) the target system requires no greater than n-byte alignment for n-byte values (eg, an int32_t * only needs 4-byte alignment, not 8-byte), and (c) the database itself will be aligned at a @@ -23,16 +23,18 @@ #include #include +#include "decimal.hpp" + namespace lost { class SerializeContext { public: - SerializeContext(bool swapIntegerEndianness, bool swapFloatEndianness) - : swapIntegerEndianness(swapIntegerEndianness), swapFloatEndianness(swapFloatEndianness) { } + SerializeContext(bool swapIntegerEndianness, bool swapDecimalEndianness) + : swapIntegerEndianness(swapIntegerEndianness), swapDecimalEndianness(swapDecimalEndianness) { } SerializeContext() : SerializeContext(false, false) { } bool swapIntegerEndianness; - bool swapFloatEndianness; + bool swapDecimalEndianness; std::vector buffer; }; @@ -66,8 +68,8 @@ void SwapEndianness(unsigned char *buffer) { } /// Swap the endianness of a value if necessary. Uses -/// LOST_DATABASE_{SOURCE,TARGET}_INTEGER_ENDIANNESS to determine to switch all values but float and -/// double, which use LOST_DATABASE_{SOURCE,TARGET}_FLOAT_ENDIANNESS. +/// LOST_DATABASE_{SOURCE,TARGET}_INTEGER_ENDIANNESS to determine to switch all values but decimal +/// which uses LOST_DATABASE_{SOURCE,TARGET}_DECIMAL_ENDIANNESS. template void SwapEndiannessIfNecessary(unsigned char *buffer, SerializeContext *ser) { if (ser->swapIntegerEndianness) { @@ -78,16 +80,9 @@ void SwapEndiannessIfNecessary(unsigned char *buffer, SerializeContext *ser) { // template specializations template <> -inline void SwapEndiannessIfNecessary(unsigned char *buffer, SerializeContext *ser) { - if (ser->swapFloatEndianness) { - SwapEndianness(buffer); - } -} - -template <> -inline void SwapEndiannessIfNecessary(unsigned char *buffer, SerializeContext *ser) { - if (ser->swapFloatEndianness) { - SwapEndianness(buffer); +inline void SwapEndiannessIfNecessary(unsigned char *buffer, SerializeContext *ser) { + if (ser->swapDecimalEndianness) { + SwapEndianness(buffer); } } diff --git a/src/star-id-private.hpp b/src/star-id-private.hpp index c910f898..943997c9 100644 --- a/src/star-id-private.hpp +++ b/src/star-id-private.hpp @@ -17,7 +17,7 @@ namespace lost { class IRUnidentifiedCentroid { public: IRUnidentifiedCentroid(const Star &star, int16_t index) - : bestAngleFrom90(std::numeric_limits::max()), // should be infinity + : bestAngleFrom90(std::numeric_limits::max()), // should be infinity bestStar1(0,0), bestStar2(0,0), index(index), star(&star) { @@ -29,7 +29,7 @@ class IRUnidentifiedCentroid { : bestStar1(0,0), bestStar2(0,0), index(-1) { } - float bestAngleFrom90; /// For the pair of other centroids forming the triangular angle closest to 90 degrees, how far from 90 degrees it is (in radians) + decimal bestAngleFrom90; /// For the pair of other centroids forming the triangular angle closest to 90 degrees, how far from 90 degrees it is (in radians) StarIdentifier bestStar1; /// One star corresponding to bestAngleFrom90 StarIdentifier bestStar2; /// The other star corresponding to bestAngleFrom90 int16_t index; /// Index into list of all centroids @@ -37,10 +37,10 @@ class IRUnidentifiedCentroid { private: // possible improvement: Use a tree map here to allow binary search - std::vector> identifiedStarsInRange; + std::vector> identifiedStarsInRange; private: - float VerticalAnglesToAngleFrom90(float v1, float v2); + decimal VerticalAnglesToAngleFrom90(decimal v1, decimal v2); public: void AddIdentifiedStar(const StarIdentifier &starId, const Stars &stars); @@ -49,15 +49,15 @@ class IRUnidentifiedCentroid { std::vector IdentifyThirdStar(const PairDistanceKVectorDatabase &db, const Catalog &catalog, int16_t catalogIndex1, int16_t catalogIndex2, - float distance1, float distance2, - float tolerance); + decimal distance1, decimal distance2, + decimal tolerance); int IdentifyRemainingStarsPairDistance(StarIdentifiers *, const Stars &, const PairDistanceKVectorDatabase &, const Catalog &, const Camera &, - float tolerance); + decimal tolerance); } diff --git a/src/star-id.cpp b/src/star-id.cpp index b0787313..7de2c303 100644 --- a/src/star-id.cpp +++ b/src/star-id.cpp @@ -46,17 +46,17 @@ StarIdentifiers GeometricVotingStarIdAlgorithm::Go( // TODO: find a faster way to do this: std::vector votedInPair(catalog.size(), false); Vec3 jSpatial = camera.CameraToSpatial(stars[j].position).Normalize(); - float greatCircleDistance = AngleUnit(iSpatial, jSpatial); + decimal greatCircleDistance = AngleUnit(iSpatial, jSpatial); //give a greater range for min-max Query for bigger radius (GreatCircleDistance) - float lowerBoundRange = greatCircleDistance - tolerance; - float upperBoundRange = greatCircleDistance + tolerance; + decimal lowerBoundRange = greatCircleDistance - tolerance; + decimal upperBoundRange = greatCircleDistance + tolerance; const int16_t *upperBoundSearch; const int16_t *lowerBoundSearch = vectorDatabase.FindPairsLiberal( lowerBoundRange, upperBoundRange, &upperBoundSearch); //loop from lowerBoundSearch till numReturnedPairs, add one vote to each star in the pairs in the datastructure for (const int16_t *k = lowerBoundSearch; k != upperBoundSearch; k++) { if ((k - lowerBoundSearch) % 2 == 0) { - float actualAngle = AngleUnit(catalog[*k].spatial, catalog[*(k+1)].spatial); + decimal actualAngle = AngleUnit(catalog[*k].spatial, catalog[*(k+1)].spatial); assert(actualAngle <= greatCircleDistance + tolerance * 2); assert(actualAngle >= greatCircleDistance - tolerance * 2); } @@ -82,7 +82,7 @@ StarIdentifiers GeometricVotingStarIdAlgorithm::Go( } } // if (i == 542) { - // for (float dist : vectorDatabase.StarDistances(9085, catalog)) { + // for (decimal dist : vectorDatabase.StarDistances(9085, catalog)) { // printf("Actual 9085 distance: %f\n", dist); // } // puts("Debug star."); @@ -112,17 +112,17 @@ StarIdentifiers GeometricVotingStarIdAlgorithm::Go( // Calculate distance between catalog stars CatalogStar first = catalog[identified[i].catalogIndex]; CatalogStar second = catalog[identified[j].catalogIndex]; - float cDist = AngleUnit(first.spatial, second.spatial); + decimal cDist = AngleUnit(first.spatial, second.spatial); Star firstIdentified = stars[identified[i].starIndex]; Star secondIdentified = stars[identified[j].starIndex]; Vec3 firstSpatial = camera.CameraToSpatial(firstIdentified.position); Vec3 secondSpatial = camera.CameraToSpatial(secondIdentified.position); - float sDist = Angle(firstSpatial, secondSpatial); + decimal sDist = Angle(firstSpatial, secondSpatial); //if sDist is in the range of (distance between stars in the image +- R) //add a vote for the match - if (abs(sDist - cDist) < tolerance) { + if (DECIMAL_ABS(sDist - cDist) < tolerance) { verificationVotes[i]++; verificationVotes[j]++; } @@ -265,8 +265,8 @@ std::unordered_multimap PairDistanceQueryToMap(const int16_t * return result; } -float IRUnidentifiedCentroid::VerticalAnglesToAngleFrom90(float v1, float v2) { - return abs(FloatModulo(v1-v2, M_PI) - M_PI_2); +decimal IRUnidentifiedCentroid::VerticalAnglesToAngleFrom90(decimal v1, decimal v2) { + return DECIMAL_ABS(DecimalModulo(v1-v2, DECIMAL_M_PI) - DECIMAL_M_PI_2); } /** @@ -276,10 +276,10 @@ float IRUnidentifiedCentroid::VerticalAnglesToAngleFrom90(float v1, float v2) { void IRUnidentifiedCentroid::AddIdentifiedStar(const StarIdentifier &starId, const Stars &stars) { const Star &otherStar = stars[starId.starIndex]; Vec2 positionDifference = otherStar.position - star->position; - float angleFromVertical = atan2(positionDifference.y, positionDifference.x); + decimal angleFromVertical = DECIMAL_ATAN2(positionDifference.y, positionDifference.x); for (const auto &otherPair : identifiedStarsInRange) { - float curAngleFrom90 = VerticalAnglesToAngleFrom90(otherPair.first, angleFromVertical); + decimal curAngleFrom90 = VerticalAnglesToAngleFrom90(otherPair.first, angleFromVertical); if (curAngleFrom90 < bestAngleFrom90) { bestAngleFrom90 = curAngleFrom90; bestStar1 = starId; @@ -298,17 +298,17 @@ void IRUnidentifiedCentroid::AddIdentifiedStar(const StarIdentifier &starId, con */ std::vector::iterator> FindUnidentifiedCentroidsInRange( std::vector *centroids, const Star &star, const Camera &camera, - float minDistance, float maxDistance) { + decimal minDistance, decimal maxDistance) { Vec3 ourSpatial = camera.CameraToSpatial(star.position).Normalize(); - float minCos = cos(maxDistance); - float maxCos = cos(minDistance); + decimal minCos = DECIMAL_COS(maxDistance); + decimal maxCos = DECIMAL_COS(minDistance); std::vector::iterator> result; for (auto it = centroids->begin(); it != centroids->end(); ++it) { Vec3 theirSpatial = camera.CameraToSpatial((*it)->star->position).Normalize(); - float angleCos = ourSpatial * theirSpatial; + decimal angleCos = ourSpatial * theirSpatial; if (angleCos >= minCos && angleCos <= maxCos) { result.push_back(it); } @@ -320,19 +320,19 @@ std::vector::iterator> FindUnidentifiedCen // // Find the first centroid that is within range of the given centroid. // auto firstInRange = std::lower_bound(centroids->begin(), centroids->end(), star.position.x - maxDistance, - // [](const IRUnidentifiedCentroid ¢roid, float x) { + // [](const IRUnidentifiedCentroid ¢roid, decimal x) { // return centroid.star.position.x < x; // }); // // Find the first centroid that is not within range of the given centroid. // auto firstNotInRange = std::lower_bound(firstInRange, centroids->end(), star.position.x + maxDistance, - // [](const IRUnidentifiedCentroid ¢roid, float x) { + // [](const IRUnidentifiedCentroid ¢roid, decimal x) { // return centroid.star.position.x <= x; // }); // // Copy the pointers to the stars into the result vector. // for (auto it = firstInRange; it != firstNotInRange; ++it) { - // float distance = Distance(star.position, it->star.position); + // decimal distance = Distance(star.position, it->star.position); // if (distance >= minDistance && distance <= maxDistance) { // result.push_back(&*it); // } @@ -350,8 +350,8 @@ std::vector::iterator> FindUnidentifiedCen void AddToAllUnidentifiedCentroids(const StarIdentifier &starId, const Stars &stars, std::vector *aboveThresholdCentroids, std::vector *belowThresholdCentroids, - float minDistance, float maxDistance, - float angleFrom90Threshold, + decimal minDistance, decimal maxDistance, + decimal angleFrom90Threshold, const Camera &camera) { std::vector nowBelowThreshold; // centroid indices newly moved above the threshold @@ -384,8 +384,8 @@ void AddToAllUnidentifiedCentroids(const StarIdentifier &starId, const Stars &st std::vector IdentifyThirdStar(const PairDistanceKVectorDatabase &db, const Catalog &catalog, int16_t catalogIndex1, int16_t catalogIndex2, - float distance1, float distance2, - float tolerance) { + decimal distance1, decimal distance2, + decimal tolerance) { const int16_t *query1End; const int16_t *query1 = db.FindPairsExact(catalog, distance1-tolerance, distance1+tolerance, &query1End); @@ -404,7 +404,7 @@ std::vector IdentifyThirdStar(const PairDistanceKVectorDatabase &db, Vec3 candidateSpatial = catalog[*candidateIt].spatial; - float angle2 = AngleUnit(candidateSpatial, spatial2); + decimal angle2 = AngleUnit(candidateSpatial, spatial2); // check distance to second star if (!(angle2 >= distance2-tolerance && angle2 <= distance2+tolerance)) { @@ -412,7 +412,7 @@ std::vector IdentifyThirdStar(const PairDistanceKVectorDatabase &db, } // check spectrality - float spectralTorch = cross * candidateSpatial; + decimal spectralTorch = cross * candidateSpatial; // if they are nearly coplanar, don't need to check spectrality // TODO: Implement ^^. Not high priority, since always checking spectrality is conservative. if (spectralTorch <= 0) { @@ -440,7 +440,7 @@ IRUnidentifiedCentroid *SelectNextUnidentifiedCentroid(std::vectorbestAngleFrom90 < b->bestAngleFrom90; }); - // 10 is arbitrary; but really it should be less than M_PI_2 when set + // 10 is arbitrary; but really it should be less than DECIMAL_M_PI_2 when set if (bestAboveThreshold != aboveThresholdCentroids->end() && (*bestAboveThreshold)->bestAngleFrom90 < 10) { auto result = *bestAboveThreshold; aboveThresholdCentroids->erase(bestAboveThreshold); @@ -450,7 +450,7 @@ IRUnidentifiedCentroid *SelectNextUnidentifiedCentroid(std::vectorstar->position); Vec3 spatial1 = camera.CameraToSpatial(stars[nextUnidentifiedCentroid->bestStar1.starIndex].position); Vec3 spatial2 = camera.CameraToSpatial(stars[nextUnidentifiedCentroid->bestStar2.starIndex].position); - float d1 = Angle(spatial1, unidentifiedSpatial); - float d2 = Angle(spatial2, unidentifiedSpatial); - float spectralTorch = spatial1.CrossProduct(spatial2) * unidentifiedSpatial; + decimal d1 = Angle(spatial1, unidentifiedSpatial); + decimal d2 = Angle(spatial2, unidentifiedSpatial); + decimal spectralTorch = spatial1.CrossProduct(spatial2) * unidentifiedSpatial; // find all the catalog stars that are in both annuli // flip arguments for appropriate spectrality. @@ -581,9 +581,9 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( DeserializeContext des(databaseBuffer); PairDistanceKVectorDatabase vectorDatabase(&des); - // smallest normal single-precision float is around 10^-38 so we should be all good. See + // smallest normal single-precision decimal is around 10^-38 so we should be all good. See // Analytic_Star_Pattern_Probability on the HSL wiki for details. - float expectedMismatchesConstant = pow(numFalseStars, 4) * pow(tolerance, 5) / 2 / pow(M_PI, 2); + decimal expectedMismatchesConstant = DECIMAL_POW(numFalseStars, 4) * DECIMAL_POW(tolerance, 5) / 2 / DECIMAL_POW(DECIMAL_M_PI, 2); // this iteration technique is described in the Pyramid paper. Briefly: i will always be the // lowest index, then dj and dk are how many indexes ahead the j-th star is from the i-th, and @@ -628,17 +628,17 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( Vec3 jSpatial = camera.CameraToSpatial(stars[j].position).Normalize(); Vec3 kSpatial = camera.CameraToSpatial(stars[k].position).Normalize(); - float ijDist = AngleUnit(iSpatial, jSpatial); + decimal ijDist = AngleUnit(iSpatial, jSpatial); - float iSinInner = sin(Angle(jSpatial - iSpatial, kSpatial - iSpatial)); - float jSinInner = sin(Angle(iSpatial - jSpatial, kSpatial - jSpatial)); - float kSinInner = sin(Angle(iSpatial - kSpatial, jSpatial - kSpatial)); + decimal iSinInner = DECIMAL_SIN(Angle(jSpatial - iSpatial, kSpatial - iSpatial)); + decimal jSinInner = DECIMAL_SIN(Angle(iSpatial - jSpatial, kSpatial - jSpatial)); + decimal kSinInner = DECIMAL_SIN(Angle(iSpatial - kSpatial, jSpatial - kSpatial)); // if we made it this far, all 6 angles are confirmed! Now check // that this match would not often occur due to chance. // See Analytic_Star_Pattern_Probability on the HSL wiki for details - float expectedMismatches = expectedMismatchesConstant - * sin(ijDist) + decimal expectedMismatches = expectedMismatchesConstant + * DECIMAL_SIN(ijDist) / kSinInner / std::max(std::max(iSinInner, jSinInner), kSinInner); @@ -652,11 +652,11 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( // sign of determinant, to detect flipped patterns bool spectralTorch = iSpatial.CrossProduct(jSpatial)*kSpatial > 0; - float ikDist = AngleUnit(iSpatial, kSpatial); - float irDist = AngleUnit(iSpatial, rSpatial); - float jkDist = AngleUnit(jSpatial, kSpatial); - float jrDist = AngleUnit(jSpatial, rSpatial); - float krDist = AngleUnit(kSpatial, rSpatial); // TODO: we don't really need to + decimal ikDist = AngleUnit(iSpatial, kSpatial); + decimal irDist = AngleUnit(iSpatial, rSpatial); + decimal jkDist = AngleUnit(jSpatial, kSpatial); + decimal jrDist = AngleUnit(jSpatial, rSpatial); + decimal krDist = AngleUnit(kSpatial, rSpatial); // TODO: we don't really need to // check krDist, if k has been // verified by i and j it's fine. @@ -704,7 +704,7 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( } // small optimization: We can calculate jk before iterating through r, so we will! - float jkCandidateDist = AngleUnit(jCandidateSpatial, kCandidateSpatial); + decimal jkCandidateDist = AngleUnit(jCandidateSpatial, kCandidateSpatial); if (jkCandidateDist < jkDist - tolerance || jkCandidateDist > jkDist + tolerance) { continue; } @@ -715,8 +715,8 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( for (auto rCandidateIt = irMap.equal_range(iCandidate); rCandidateIt.first != rCandidateIt.second; rCandidateIt.first++) { int rCandidate = rCandidateIt.first->second; const Vec3 &rCandidateSpatial = catalog[rCandidate].spatial; - float jrCandidateDist = AngleUnit(jCandidateSpatial, rCandidateSpatial); - float krCandidateDist; + decimal jrCandidateDist = AngleUnit(jCandidateSpatial, rCandidateSpatial); + decimal krCandidateDist; if (jrCandidateDist < jrDist - tolerance || jrCandidateDist > jrDist + tolerance) { continue; } @@ -745,7 +745,8 @@ StarIdentifiers PyramidStarIdAlgorithm::Go( } if (iMatch != -1) { - printf("Matched unique pyramid!\nExpected mismatches: %e\n", expectedMismatches); + std::cout.precision(6); + std::cout << "Matched unique pyramid!" << std::endl << "Expected mismatches: " << std::scientific << expectedMismatches << std::endl << std::fixed; identified.push_back(StarIdentifier(i, iMatch)); identified.push_back(StarIdentifier(j, jMatch)); identified.push_back(StarIdentifier(k, kMatch)); diff --git a/src/star-id.hpp b/src/star-id.hpp index d82727c9..32c0dff5 100644 --- a/src/star-id.hpp +++ b/src/star-id.hpp @@ -39,9 +39,9 @@ class GeometricVotingStarIdAlgorithm : public StarIdAlgorithm { /** * @param tolerance Angular tolerance (Two inter-star distances are considered the same if within this many radians) */ - explicit GeometricVotingStarIdAlgorithm(float tolerance): tolerance(tolerance) { }; + explicit GeometricVotingStarIdAlgorithm(decimal tolerance): tolerance(tolerance) { }; private: - float tolerance; + decimal tolerance; }; @@ -60,13 +60,13 @@ class PyramidStarIdAlgorithm final : public StarIdAlgorithm { * @param maxMismatchProbability The maximum allowable probability for any star to be mis-id'd. * @param cutoff Maximum number of pyramids to iterate through before giving up. */ - PyramidStarIdAlgorithm(float tolerance, int numFalseStars, float maxMismatchProbability, long cutoff) + PyramidStarIdAlgorithm(decimal tolerance, int numFalseStars, decimal maxMismatchProbability, long cutoff) : tolerance(tolerance), numFalseStars(numFalseStars), maxMismatchProbability(maxMismatchProbability), cutoff(cutoff) { }; private: - float tolerance; + decimal tolerance; int numFalseStars; - float maxMismatchProbability; + decimal maxMismatchProbability; long cutoff; }; diff --git a/src/star-utils.cpp b/src/star-utils.cpp index af12dd9b..29e9e6ae 100644 --- a/src/star-utils.cpp +++ b/src/star-utils.cpp @@ -14,7 +14,7 @@ bool CatalogStarMagnitudeCompare(const CatalogStar &a, const CatalogStar &b) { return a.magnitude < b.magnitude; } -Catalog NarrowCatalog(const Catalog &catalog, int maxMagnitude, int maxStars, float minSeparation) { +Catalog NarrowCatalog(const Catalog &catalog, int maxMagnitude, int maxStars, decimal minSeparation) { Catalog result; for (int i = 0; i < (int)catalog.size(); i++) { if (catalog[i].magnitude <= maxMagnitude) { @@ -70,7 +70,7 @@ Catalog::const_iterator FindNamedStar(const Catalog &catalog, int name) { void SerializeCatalogStar(SerializeContext *ser, const CatalogStar &catalogStar, bool inclMagnitude, bool inclName) { SerializeVec3(ser, catalogStar.spatial); if (inclMagnitude) { - SerializePrimitive(ser, catalogStar.magnitude); + SerializePrimitive(ser, catalogStar.magnitude); } if (inclName) { // TODO: double check that bools aren't some special bitwise thing in C++ @@ -87,7 +87,7 @@ CatalogStar DeserializeCatalogStar(DeserializeContext *des, bool inclMagnitude, CatalogStar result; result.spatial = DeserializeVec3(des); if (inclMagnitude) { - result.magnitude = DeserializePrimitive(des); + result.magnitude = DeserializePrimitive(des); } else { result.magnitude = -424242; // TODO, what to do about special values, since there's no good ones for ints. } @@ -144,8 +144,8 @@ Catalog DeserializeCatalog(DeserializeContext *des, bool *inclMagnitudeReturn, b return result; } -float MagToBrightness(int mag) { - return pow(10.0, -mag/250.0); +decimal MagToBrightness(int mag) { + return DECIMAL_POW(DECIMAL(10.0), -mag/DECIMAL(250.0)); } } diff --git a/src/star-utils.hpp b/src/star-utils.hpp index 0af225bd..7d07ebbf 100644 --- a/src/star-utils.hpp +++ b/src/star-utils.hpp @@ -19,7 +19,7 @@ class CatalogStar { * @param magnitude See ::magnitude * @param name See ::name */ - CatalogStar(float raj2000, float dej2000, int magnitude, int name) : + CatalogStar(decimal raj2000, decimal dej2000, int magnitude, int name) : spatial(SphericalToSpatial(raj2000, dej2000)), magnitude(magnitude), name(name) {} CatalogStar(Vec3 spatial, int magnitude, int name) : @@ -48,21 +48,21 @@ class CatalogStar { */ class Star { public: - Star(float x, float y, float radiusX, float radiusY, int magnitude) : + Star(decimal x, decimal y, decimal radiusX, decimal radiusY, int magnitude) : position({x, y}), radiusX(radiusX), radiusY(radiusY), magnitude(magnitude) {}; /// Convenience constructor that sets Star.radiusY = radiusX and Star.magnitude = 0 - Star(float x, float y, float radiusX) : Star(x, y, radiusX, radiusX, 0) {}; + Star(decimal x, decimal y, decimal radiusX) : Star(x, y, radiusX, radiusX, 0) {}; /// Create a zeroed-out star. Fields should be set immediately after construction. - Star() : Star(0.0, 0.0, 0.0) {}; + Star() : Star(DECIMAL(0.0), DECIMAL(0.0), DECIMAL(0.0)) {}; /// The (x,y) pixel coordinates in the image (top left is 0,0) Vec2 position; /// Approximate horizontal radius of the bright area in pixels. - float radiusX; + decimal radiusX; /// Approximate vertical radius of the bright area in pixels. - float radiusY; + decimal radiusY; /** * A relative measure of magnitude of the star. Larger is brighter. * It's impossible to tell the true magnitude of the star from the image, without really good camera calibration. Anyway, this field is not meant to correspond to the usual measurement of magnitude. Instead, it's just some measure of brightness which may be specific to the centroiding algorithm. For example, it might be the total number of bright pixels in the star. @@ -81,7 +81,7 @@ class StarIdentifier { : starIndex(starIndex), catalogIndex(catalogIndex), weight(weight) { }; /// Sets StarIdentifier.weight = 1 StarIdentifier(int starIndex, int catalogIndex) - : StarIdentifier(starIndex, catalogIndex, 1.0f) { }; + : StarIdentifier(starIndex, catalogIndex, DECIMAL(1.0)) { }; // does not check weight bool operator==(const StarIdentifier& other) const { @@ -94,7 +94,7 @@ class StarIdentifier { /// An index into an array of CatalogStar objects. int catalogIndex; /// A weight indicating the confidence of this idenification. Often just 1. - float weight; + decimal weight; }; typedef std::vector Catalog; @@ -108,7 +108,7 @@ Catalog::const_iterator FindNamedStar(const Catalog &catalog, int name); /// returns some relative brightness measure, which is proportional to the total number of photons received from a star. /// As always, the magnitude is actually 100* the usual magnitude -float MagToBrightness(int magnitude); +decimal MagToBrightness(int magnitude); /** * Remove unwanted stars from an unfiltered catalog. @@ -129,7 +129,7 @@ float MagToBrightness(int magnitude); * mark, so that star-id algorithms can reason about them instead of thinking they are false stars, * but no such provision is made yet. */ -Catalog NarrowCatalog(const Catalog &, int maxMagnitude, int maxStars, float minSeparation); +Catalog NarrowCatalog(const Catalog &, int maxMagnitude, int maxStars, decimal minSeparation); } diff --git a/test/geometry.cpp b/test/geometry.cpp index 252434fe..7e4db6d7 100644 --- a/test/geometry.cpp +++ b/test/geometry.cpp @@ -11,8 +11,8 @@ using namespace lost; // NOLINT TEST_CASE("Convert coordinates: pixel -> spatial -> pixel", "[geometry]") { Camera camera(100, 512, 1024); - float expectedX = GENERATE(142, 90, 512, 255); - float expectedY = GENERATE(18, 512, 0, 800); + decimal expectedX = GENERATE(142, 90, 512, 255); + decimal expectedY = GENERATE(18, 512, 0, 800); Vec3 spatial = camera.CameraToSpatial({expectedX, expectedY}); Vec2 actualPixels = camera.SpatialToCamera(spatial); @@ -68,9 +68,9 @@ TEST_CASE("Angle from camera", "[geometry]") { TEST_CASE("spherical -> quaternion -> spherical", "[geometry]") { // 0.1 instead of 0, because at 0 it might sometimes return 2PI, which is fine for most // circumstances. Also, at 0, the epsilon=0, so there's no tolerance in Approx by default! - float ra = DegToRad(GENERATE(28.9, 83.2, 14.0, 0.1, 329.8)); - float de = DegToRad(GENERATE(7.82, 9.88, 88.8, 0.1, -72.0, -9.9)); - float roll = DegToRad(GENERATE(9.38, 300.9, 37.8, 199.9)); + decimal ra = DegToRad(GENERATE(28.9, 83.2, 14.0, 0.1, 329.8)); + decimal de = DegToRad(GENERATE(7.82, 9.88, 88.8, 0.1, -72.0, -9.9)); + decimal roll = DegToRad(GENERATE(9.38, 300.9, 37.8, 199.9)); Quaternion quat = SphericalToQuaternion(ra, de, roll); @@ -82,10 +82,10 @@ TEST_CASE("spherical -> quaternion -> spherical", "[geometry]") { } TEST_CASE("spherical -> spatial -> spherical", "[geometry]") { - float ra = DegToRad(GENERATE(28.9, 83.2, 14.0, 0.1, 329.8)); - float de = DegToRad(GENERATE(7.82, 9.88, 88.8, 0.1, -72.0, -9.9)); + decimal ra = DegToRad(GENERATE(28.9, 83.2, 14.0, 0.1, 329.8)); + decimal de = DegToRad(GENERATE(7.82, 9.88, 88.8, 0.1, -72.0, -9.9)); - float raOut, deOut; + decimal raOut, deOut; SpatialToSpherical(SphericalToSpatial(ra, de), &raOut, &deOut); CHECK(ra == Approx(raOut)); @@ -93,9 +93,9 @@ TEST_CASE("spherical -> spatial -> spherical", "[geometry]") { } TEST_CASE("quat -> dcm -> quat", "[geometry]") { - float ra = GENERATE(take(5, random(0.1, 3.14*2))); - float de = GENERATE(take(5, random(-3.14, 3.14))); - float roll = GENERATE(take(5, random(0.1, 3.14*2))); + decimal ra = GENERATE(take(5, random(0.1, 3.14*2))); + decimal de = GENERATE(take(5, random(-3.14, 3.14))); + decimal roll = GENERATE(take(5, random(0.1, 3.14*2))); Quaternion quat1 = SphericalToQuaternion(ra, de, roll).Canonicalize(); Mat3 dcm = QuaternionToDCM(quat1); diff --git a/test/identify-remaining-stars.cpp b/test/identify-remaining-stars.cpp index 20c2c468..21969ece 100644 --- a/test/identify-remaining-stars.cpp +++ b/test/identify-remaining-stars.cpp @@ -13,19 +13,19 @@ using namespace lost; // NOLINT TEST_CASE("IRUnidentifiedCentroid simple orthogonal", "[identify-remaining] [fast]") { IRUnidentifiedCentroid centroid(elevenStars[3], 0); - REQUIRE(centroid.bestAngleFrom90 > 9e9); + REQUIRE(centroid.bestAngleFrom90 > DECIMAL(9e9)); centroid.AddIdentifiedStar(elevenStarIds[1], elevenStars); // one star is not enough to get angle - REQUIRE(centroid.bestAngleFrom90 > 9e9); + REQUIRE(centroid.bestAngleFrom90 > DECIMAL(9e9)); centroid.AddIdentifiedStar(elevenStarIds[2], elevenStars); // we've set them up to be almost orthogonal - REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(1e-6)); + REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(DECIMAL(1e-6))); REQUIRE(((centroid.bestStar1 == elevenStarIds[1] && centroid.bestStar2 == elevenStarIds[2]) || (centroid.bestStar1 == elevenStarIds[2] && centroid.bestStar2 == elevenStarIds[1]))); // adding another, non-orthogonal one shouldn't break it centroid.AddIdentifiedStar(elevenStarIds[0], elevenStars); - REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(1e-6)); + REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(DECIMAL(1e-6))); REQUIRE(((centroid.bestStar1 == elevenStarIds[1] && centroid.bestStar2 == elevenStarIds[2]) || (centroid.bestStar1 == elevenStarIds[2] && centroid.bestStar2 == elevenStarIds[1]))); } @@ -34,26 +34,26 @@ TEST_CASE("IRUnidentifiedCentroid not orthogonal until they are", "[identify-rem IRUnidentifiedCentroid centroid(elevenStars[3], 0); centroid.AddIdentifiedStar(elevenStarIds[1], elevenStars); centroid.AddIdentifiedStar(elevenStarIds[0], elevenStars); - REQUIRE(centroid.bestAngleFrom90 == Approx(M_PI_4)); + REQUIRE(centroid.bestAngleFrom90 == Approx(DECIMAL_M_PI_4)); centroid.AddIdentifiedStar(elevenStarIds[6], elevenStars); - REQUIRE(centroid.bestAngleFrom90 == Approx(M_PI_4)); + REQUIRE(centroid.bestAngleFrom90 == Approx(DECIMAL_M_PI_4)); centroid.AddIdentifiedStar(elevenStarIds[8], elevenStars); - REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(1e-6)); + REQUIRE(centroid.bestAngleFrom90 == Approx(0).margin(DECIMAL(1e-6))); } TEST_CASE("IRUnidentifiedCentroid obtuse angle", "[identify-remaining] [fast]") { IRUnidentifiedCentroid centroid(elevenStars[3], 0); centroid.AddIdentifiedStar(elevenStarIds[1], elevenStars); centroid.AddIdentifiedStar(elevenStarIds[6], elevenStars); - REQUIRE(centroid.bestAngleFrom90 == Approx(M_PI_4)); + REQUIRE(centroid.bestAngleFrom90 == Approx(DECIMAL_M_PI_4)); } // TODO: Tests for FindAllInRange if we ever make the logic more complicated std::vector IdentifyThirdStarTest(const Catalog &catalog, int16_t catalogName1, int16_t catalogName2, - float dist1, float dist2, float tolerance) { + decimal dist1, decimal dist2, decimal tolerance) { SerializeContext ser; - SerializePairDistanceKVector(&ser, integralCatalog, 0, M_PI, 1000); + SerializePairDistanceKVector(&ser, integralCatalog, 0, DECIMAL_M_PI, 1000); DeserializeContext des(ser.buffer.data()); auto cs1 = FindNamedStar(catalog, catalogName1); auto cs2 = FindNamedStar(catalog, catalogName2); @@ -70,8 +70,8 @@ TEST_CASE("IdentifyThirdStar", "[identify-remaining] [fast]") { // TODO: does it std::vector stars = IdentifyThirdStarTest(integralCatalog, 42, 44, // (1,0,0), (0,1,0) - M_PI_2, M_PI_2, - 1e-6); + DECIMAL_M_PI_2, DECIMAL_M_PI_2, + DECIMAL(1e-6)); REQUIRE(stars.size() == 1); REQUIRE(integralCatalog[stars[0]].name == 50); @@ -82,8 +82,8 @@ TEST_CASE("IdentifyThirdStar with tolerance", "[identify-remaining] [fast]") { // try it again where we actually need the tolerance std::vector stars = IdentifyThirdStarTest(integralCatalog, 42, 44, // (1,0,0), (0,1,0) - M_PI_2 - DegToRad(1.0), M_PI_2 + DegToRad(1.0), - 0.1); + DECIMAL_M_PI_2 - DegToRad(1.0), DECIMAL_M_PI_2 + DegToRad(1.0), + DECIMAL(0.1)); REQUIRE(stars.size() == 1); REQUIRE(integralCatalog[stars[0]].name == 50); } @@ -91,8 +91,8 @@ TEST_CASE("IdentifyThirdStar with tolerance", "[identify-remaining] [fast]") { TEST_CASE("IdentifyThirdStar reversed spectrality", "[identify-remaining] [fast]") { std::vector stars = IdentifyThirdStarTest(integralCatalog, 44, 42, // (0,1,0), (1,0,0) - M_PI_2, M_PI_2, - 1e-6); + DECIMAL_M_PI_2, DECIMAL_M_PI_2, + DECIMAL(1e-6)); REQUIRE(stars.size() == 1); REQUIRE(integralCatalog[stars[0]].name == 58); } @@ -100,16 +100,16 @@ TEST_CASE("IdentifyThirdStar reversed spectrality", "[identify-remaining] [fast] TEST_CASE("IdentifyThirdStar no third star", "[identify-remaining] [fast]") { std::vector stars = IdentifyThirdStarTest(integralCatalog, 42, 44, // (1,0,0), (0,1,0) - 1, M_PI_2, - 1e-6); + 1, DECIMAL_M_PI_2, + DECIMAL(1e-6)); REQUIRE(stars.size() == 0); } TEST_CASE("IdentifyThirdStar just out of tolerance", "[identify-remaining] [fast]") { std::vector stars2 = IdentifyThirdStarTest(integralCatalog, 42, 44, // (1,0,0), (0,1,0) - M_PI_2 - 2e-6, M_PI_2, - 1e-6); + DECIMAL_M_PI_2 - DECIMAL(2e-6), DECIMAL_M_PI_2, + DECIMAL(1e-6)); REQUIRE(stars2.size() == 0); } @@ -140,7 +140,7 @@ TEST_CASE("IdentifyRemainingStars fuzz", "[identify-remaining] [fuzz]") { int numFakeStars = 100; std::default_random_engine rng(GENERATE(take(kIdentifyRemainingNumImages, random(0, 1000000)))); - std::uniform_real_distribution yDist(0.0, 256.0); + std::uniform_real_distribution yDist(DECIMAL(0.0), DECIMAL(256.0)); std::uniform_int_distribution starIndexDist(0, numFakeStars - 1); std::uniform_int_distribution moreStartingStars(0, 1); @@ -149,8 +149,8 @@ TEST_CASE("IdentifyRemainingStars fuzz", "[identify-remaining] [fuzz]") { Catalog fakeCatalog; for (int i = 0; i < numFakeStars; i++) { // smolCamera has width 256 - float x = i * 256.0 / numFakeStars; - float y = yDist(rng); + decimal x = i * DECIMAL(256.0) / numFakeStars; + decimal y = yDist(rng); fakeCentroids.emplace_back(x, y, 1); fakeCatalog.emplace_back(smolCamera.CameraToSpatial({x, y}).Normalize(), 1, i); fakeStarIds.emplace_back(i, i); @@ -166,11 +166,11 @@ TEST_CASE("IdentifyRemainingStars fuzz", "[identify-remaining] [fuzz]") { } SerializeContext ser; - SerializePairDistanceKVector(&ser, fakeCatalog, 0, M_PI, 1000); + SerializePairDistanceKVector(&ser, fakeCatalog, 0, DECIMAL_M_PI, 1000); DeserializeContext des(ser.buffer.data()); PairDistanceKVectorDatabase db(&des); - int numIdentified = IdentifyRemainingStarsPairDistance(&someFakeStarIds, fakeCentroids, db, fakeCatalog, smolCamera, 1e-5); + int numIdentified = IdentifyRemainingStarsPairDistance(&someFakeStarIds, fakeCentroids, db, fakeCatalog, smolCamera, DECIMAL(1e-5)); REQUIRE(numIdentified == numFakeStars - fakePatternSize); REQUIRE(AreStarIdentifiersEquivalent(fakeStarIds, someFakeStarIds)); diff --git a/test/kvector.cpp b/test/kvector.cpp index 40b05793..ff6e83d8 100644 --- a/test/kvector.cpp +++ b/test/kvector.cpp @@ -13,37 +13,37 @@ TEST_CASE("Kvector full database stuff", "[kvector]") { const Catalog &catalog = CatalogRead(); std::vector dbBytes; SerializeContext ser; - SerializePairDistanceKVector(&ser, catalog, DegToRad(1.0), DegToRad(2.0), 100); + SerializePairDistanceKVector(&ser, catalog, DegToRad(DECIMAL(1.0)), DegToRad(DECIMAL(2.0)), 100); DeserializeContext des(ser.buffer.data()); PairDistanceKVectorDatabase db(&des); SECTION("basic consistency checks") { long lastNumReturnedPairs = 999999; - for (float i = 1.1; i < 1.99; i += 0.1) { + for (decimal i = DECIMAL(1.1); i < DECIMAL(1.99); i += DECIMAL(0.1)) { const int16_t *end; - const int16_t *pairs = db.FindPairsLiberal(i * M_PI/180.0, 2.0 * M_PI/180.0, &end); - float shortestDistance = INFINITY; + const int16_t *pairs = db.FindPairsExact(catalog, i * DECIMAL_M_PI/DECIMAL(180.0), DECIMAL(2.0) * DECIMAL_M_PI/DECIMAL(180.0), &end); + decimal shortestDistance = INFINITY; for (const int16_t *pair = pairs; pair != end; pair += 2) { - float distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); + decimal distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); if (distance < shortestDistance) { shortestDistance = distance; } - CHECK(i * M_PI/180.0 <=distance); - CHECK(distance<= 2.01 * M_PI/180.0); + CHECK(i * DECIMAL_M_PI/DECIMAL(180.0) <= distance); + CHECK(distance <= DECIMAL(2.01) * DECIMAL_M_PI/DECIMAL(180.0)); } long numReturnedPairs = (end - pairs)/2; REQUIRE(0 < numReturnedPairs); REQUIRE(numReturnedPairs < lastNumReturnedPairs); - REQUIRE(shortestDistance < (i + 0.01) * M_PI/180.0); + REQUIRE(shortestDistance < (i + DECIMAL(0.01)) * DECIMAL_M_PI/DECIMAL(180.0)); lastNumReturnedPairs = numReturnedPairs; } } SECTION("form a partition") { long totalReturnedPairs = 0; - for (float i = 1.1; i < 2.01; i+= 0.1) { + for (decimal i = DECIMAL(1.1); i < DECIMAL(2.01); i+= DECIMAL(0.1)) { const int16_t *end; - const int16_t *pairs = db.FindPairsLiberal(DegToRad(i-0.1)+0.00001, DegToRad(i)-0.00001, &end); + const int16_t *pairs = db.FindPairsLiberal(DegToRad(i-DECIMAL(0.1))+DECIMAL(0.00001), DegToRad(i)-DECIMAL(0.00001), &end); long numReturnedPairs = (end-pairs)/2; totalReturnedPairs += numReturnedPairs; } @@ -54,24 +54,24 @@ TEST_CASE("Kvector full database stuff", "[kvector]") { TEST_CASE("Tighter tolerance test", "[kvector]") { const Catalog &catalog = CatalogRead(); SerializeContext ser; - SerializePairDistanceKVector(&ser, catalog, DegToRad(0.5), DegToRad(5.0), 1000); + SerializePairDistanceKVector(&ser, catalog, DegToRad(DECIMAL(0.5)), DegToRad(DECIMAL(5.0)), 1000); DeserializeContext des(ser.buffer.data()); PairDistanceKVectorDatabase db(&des); // radius we'll request - float delta = 0.0001; + decimal delta = DECIMAL(0.0001); // radius we expect back: radius we request + width of a bin - float epsilon = delta + DegToRad(5.0 - 0.5) / 1000; + decimal epsilon = delta + DegToRad(DECIMAL(5.0) - DECIMAL(0.5)) / 1000; // in the first test_case, the ends of each request pretty much line up with the ends of the // buckets (intentionally), so that we can do the "form a partition" test. Here, however, a // request may intersect a bucket, in which case things slightly outside the requested range should // be returned. SECTION("liberal") { bool outsideRangeReturned = false; - for (float i = DegToRad(0.6); i < DegToRad(4.9); i += DegToRad(0.1228)) { + for (decimal i = DegToRad(DECIMAL(0.6)); i < DegToRad(DECIMAL(4.9)); i += DegToRad(DECIMAL(0.1228))) { const int16_t *end; const int16_t *pairs = db.FindPairsLiberal(i - delta, i + delta, &end); for (const int16_t *pair = pairs; pair != end; pair += 2) { - float distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); + decimal distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); // only need to check one side, since we're only looking for one exception. if (i - delta > distance) { outsideRangeReturned = true; @@ -84,11 +84,11 @@ TEST_CASE("Tighter tolerance test", "[kvector]") { } SECTION("exact") { bool outsideRangeReturned = false; - for (float i = DegToRad(0.6); i < DegToRad(4.9); i += DegToRad(0.1228)) { + for (decimal i = DegToRad(DECIMAL(0.6)); i < DegToRad(DECIMAL(4.9)); i += DegToRad(DECIMAL(0.1228))) { const int16_t *end; const int16_t *pairs = db.FindPairsExact(catalog, i - delta, i + delta, &end); for (const int16_t *pair = pairs; pair != end; pair += 2) { - float distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); + decimal distance = AngleUnit(catalog[pair[0]].spatial, catalog[pair[1]].spatial); // only need to check one side, since we're only looking for one exception. if (i - delta > distance) { outsideRangeReturned = true; @@ -103,21 +103,21 @@ TEST_CASE("Tighter tolerance test", "[kvector]") { TEST_CASE("3-star database, check exact results", "[kvector] [fast]") { Catalog tripleCatalog = { - CatalogStar(DegToRad(2), DegToRad(-3), 3.0, 42), - CatalogStar(DegToRad(4), DegToRad(7), 2.0, 43), - CatalogStar(DegToRad(2), DegToRad(6), 4.0, 44), + CatalogStar(DegToRad(2), DegToRad(-3), DECIMAL(3.0), 42), + CatalogStar(DegToRad(4), DegToRad(7), DECIMAL(2.0), 43), + CatalogStar(DegToRad(2), DegToRad(6), DECIMAL(4.0), 44), }; SerializeContext ser; - SerializePairDistanceKVector(&ser, tripleCatalog, DegToRad(0.5), DegToRad(20.0), 1000); + SerializePairDistanceKVector(&ser, tripleCatalog, DegToRad(DECIMAL(0.5)), DegToRad(DECIMAL(20.0)), 1000); DeserializeContext des(ser.buffer.data()); PairDistanceKVectorDatabase db(&des); REQUIRE(db.NumPairs() == 3); - float distances[] = {0.038825754, 0.15707963, 0.177976474}; + decimal distances[] = {0.038825754, 0.15707963, 0.177976474}; SECTION("liberal") { - for (float distance : distances) { + for (decimal distance : distances) { const int16_t *end; - const int16_t *pairs = db.FindPairsLiberal(distance - 1e-6, distance + 1e-6, &end); + const int16_t *pairs = db.FindPairsLiberal(distance - DECIMAL(1e-6), distance + DECIMAL(1e-6), &end); REQUIRE(end - pairs == 2); CHECK(AngleUnit(tripleCatalog[pairs[0]].spatial, tripleCatalog[pairs[1]].spatial) == Approx(distance).epsilon(1e-4)); } @@ -125,9 +125,9 @@ TEST_CASE("3-star database, check exact results", "[kvector] [fast]") { // also serves as a regression test for an off-by-one error that used to be present in exact, where it assumed the end index was inclusive instead of "off-the-end" SECTION("exact") { - for (float distance : distances) { + for (decimal distance : distances) { const int16_t *end; - const int16_t *pairs = db.FindPairsExact(tripleCatalog, distance - 1e-4, distance + 1e-4, &end); + const int16_t *pairs = db.FindPairsExact(tripleCatalog, distance - DECIMAL(1e-4), distance + DECIMAL(1e-4), &end); REQUIRE(end - pairs == 2); CHECK(AngleUnit(tripleCatalog[pairs[0]].spatial, tripleCatalog[pairs[1]].spatial) == Approx(distance).epsilon(1e-4)); } diff --git a/test/serialize.cpp b/test/serialize.cpp index c8003875..b504a054 100644 --- a/test/serialize.cpp +++ b/test/serialize.cpp @@ -10,57 +10,57 @@ using namespace lost; // NOLINT TEST_CASE("Simple serialization, deserialization of primitives", "[fast] [serialize]") { int64_t val64 = 27837492938; - float valFloat = 23.71728; + decimal valDecimal = 23.71728; SerializeContext ser; SerializePrimitive(&ser, val64); - SerializePrimitive(&ser, valFloat); + SerializePrimitive(&ser, valDecimal); DeserializeContext des(ser.buffer.data()); int64_t deserializedVal64 = DeserializePrimitive(&des); - float deserializedFloat = DeserializePrimitive(&des); + decimal deserializedDecimal = DeserializePrimitive(&des); CHECK(val64 == deserializedVal64); - CHECK(valFloat == deserializedFloat); + CHECK(valDecimal == deserializedDecimal); } TEST_CASE("Endian-swapped serialization, deserialization of primitives", "[fast] [serialize]") { int64_t val64 = 27837492938; - float valFloat = 23.71728; + decimal valDecimal = 23.71728; SerializeContext ser1(true, true); SerializePrimitive(&ser1, val64); - SerializePrimitive(&ser1, valFloat); + SerializePrimitive(&ser1, valDecimal); DeserializeContext des(ser1.buffer.data()); int64_t deserializedVal64 = DeserializePrimitive(&des); - float deserializedValFloat = DeserializePrimitive(&des); + decimal deserializedValDecimal = DeserializePrimitive(&des); CHECK(val64 != deserializedVal64); - CHECK(valFloat != deserializedValFloat); + CHECK(valDecimal != deserializedValDecimal); // but if we serialize it again, it should be back to normal! SerializeContext ser2(true, true); SerializePrimitive(&ser2, deserializedVal64); - SerializePrimitive(&ser2, deserializedValFloat); + SerializePrimitive(&ser2, deserializedValDecimal); DeserializeContext des2(ser2.buffer.data()); int64_t redeserializedVal64 = DeserializePrimitive(&des2); - float redeserializedValFloat = DeserializePrimitive(&des2); + decimal redeserializedValDecimal = DeserializePrimitive(&des2); CHECK(val64 == redeserializedVal64); - CHECK(valFloat == redeserializedValFloat); + CHECK(valDecimal == redeserializedValDecimal); } -TEST_CASE("Endian-swapped floats only", "[fast] [serialize]") { +TEST_CASE("Endian-swapped decimals only", "[fast] [serialize]") { int64_t val64 = 27837492938; - float valFloat = 23.71728; + decimal valDecimal = 23.71728; SerializeContext ser1(false, true); SerializePrimitive(&ser1, val64); - SerializePrimitive(&ser1, valFloat); + SerializePrimitive(&ser1, valDecimal); DeserializeContext des(ser1.buffer.data()); int64_t deserializedVal64 = DeserializePrimitive(&des); - float deserializedValFloat = DeserializePrimitive(&des); + decimal deserializedValDecimal = DeserializePrimitive(&des); CHECK(val64 == deserializedVal64); - CHECK(valFloat != deserializedValFloat); + CHECK(valDecimal != deserializedValDecimal); SerializeContext ser2(false, true); - SerializePrimitive(&ser2, deserializedValFloat); + SerializePrimitive(&ser2, deserializedValDecimal); DeserializeContext des2(ser2.buffer.data()); - float redeserializedValFloat = DeserializePrimitive(&des2); - CHECK(valFloat == redeserializedValFloat); + decimal redeserializedValDecimal = DeserializePrimitive(&des2); + CHECK(valDecimal == redeserializedValDecimal); } TEST_CASE("Padding", "[fast] [serialize]") {