diff --git a/.gitignore b/.gitignore index d08a4241..24f784a7 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,13 @@ *.txt *.bin +*.out.* +*.lisp +*.org +*.py +*.json +/.gdb_history + BSC5 /documentation/man-*.h diff --git a/Makefile b/Makefile index c4d1daf3..2591f297 100644 --- a/Makefile +++ b/Makefile @@ -1,15 +1,15 @@ # Copyright (c) 2020 Mark Polyakov, Karen Haining (If you edit the file, add your name here!) -# +# # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: -# +# # The above copyright notice and this permission notice shall be included in all # copies or substantial portions of the Software. -# +# # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE @@ -57,9 +57,6 @@ release: CXXFLAGS := $(RELEASE_CXXFLAGS) release: LDFLAGS := $(RELEASE_LDFLAGS) release: all -$(BSC): download-bsc.sh - ./download-bsc.sh - $(BIN): $(OBJS) $(CXX) $(LDFLAGS) -o $(BIN) $(OBJS) $(LIBS) diff --git a/documentation/pipeline.man b/documentation/pipeline.man index 969b6285..e47a661a 100644 --- a/documentation/pipeline.man +++ b/documentation/pipeline.man @@ -92,6 +92,13 @@ The field of view of the camera that took the picture (in degrees). Defaults to \fB--centroid-mag-filter\fP \fImin-mag\fP Will not consider centroids with magnitude below \fImin-mag\fP. +.TP +\fB--centroid-filter-brightest\fP \fInum-stars\fP +Remove all but the brightest \fInum-stars\fP many stars from the list of centroids before sending to +star-id. Often a better choice than \fB--centroid-mag-filter\fP, because you can ensure that you +keep enough stars to do star-id. If both this option and \fB--centroid-mag-filter\fP are provided, +then all stars satisfying both criteria are kept (intersection). + .TP \fB--database\fP \fIfilename\fP Chooses \fIfilename\fP as the database to use during star identification. @@ -149,8 +156,14 @@ Sets the horizontal resolution of the generated image(s) to \fIpixels\fP. Defaul Sets the vertical resolution of the generated image(s) to \fIpixels\fP. Defaults to 1024. .TP -\fB--generate-reference-brightness\fP \fIobserved-brightness\fP -A star with magnitude 0 in the generated image will have have brightness \fIobserved-brightness\fP. An observed brightness of 1 fully saturates the pixel at the center of the star (completely white). Defaults to 100. +\fB--generate-zero-mag-photons\fP \fInum-photoelectrons\fP +A star with magnitude 0 will cause \fInum-photoelectrons\fP many photoelectrons to be received by +the sensor. A default value of 20,000 is chosen. See Liebe's, "tutorial on star tracker accuracy" for theoretical information on how to calculate this. + +.TP +\fB--generate-saturation-photons\fP \fIsaturation-photoelectrons\fP + +When a pixel receives at least this many photoelectrons, it will appear completely white (at least before other noise is applied). Note that, because of noise, a pixel may still appear completely white if it receives less than \fIsaturation-photoelectrons\fP many photoelectrons. .TP \fB--generate-spread-stddev\fP \fIstddev\fP @@ -164,10 +177,6 @@ Enables or disables shot noise simulation in generated images. Defaults to true. \fB--generate-dark-current\fP \fInoise-level\fP Set observed brightness of dark current in the image, from 0 (no dark noise) to 1 (whole image pure white). Defaults to 0.1. cf \fB--generate-sensitivity\fP to control shot noise intensity. -.TP -\fB--generate-sensitivity\fP \fIsensitivity\fP -Controls the simulated camera sensitivity for generated images. This only has an observable effect when shot noise is enabled (cf \fB--generate-shot-noise\fP). A higher sensitivity means that fewer photons can cause the same observed brightness on the sensor, so variation in number of photons causes more shot noise. Brightness is normalized according to \fB--generate-reference-brightness\fP so that adjusting sensitivity does not change observed brightness. - .TP \fB--generate-read-noise-stddev\fP \fIstddev\fP Sets the standard deviation of Gaussian noise in the generated image(s) to \fIstddev\fP. Noise is measured in observed brightness, where 1 is the difference between pure white and pure black. Defaults to 0.05. diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 0ddb6103..7c113e41 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -1,25 +1,614 @@ #include "centroiders.hpp" +#include #include #include -#include #include -#include +#include +#include +#include +#include +// #include #include +#include // TODO: remove later, this is just lazy to implement hash for unordered_set #include #include +#include namespace lost { -// DUMMY +/** +Prediction of 1D Gaussian model for 1D LSGF method +Note: can be used as f(xi, beta) or f(yi, beta) +*/ +static float FitModel(float x, float a, float xb, float sigma); + +/// Prediction of 2D Gaussian model for 2D LSGF method +static float FitModel2D(float x, float y, float a, float xb, float yb, float sigmaX, float sigmaY); + +// Structure: (x, y, size of floodfill around (x, y)) +typedef std::vector Point; + +/** + * @brief Get a candidate point for each possible star cluster by running a floodfill through image + * + * Eliminates possibility of detecting multiple centroids for a single star + */ +std::vector FloodfillPreproc(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); + +void SubtractNoise(unsigned char *image, int imageWidth, int imageHeight){ + float sum = 0; + for (long i = 0; i < imageHeight * imageWidth; i++) { + int x = i % imageWidth; + int y = i / imageWidth; + sum += Get(x, y, image, imageWidth); + } + float noise = sum / (imageWidth * imageHeight); + for (long i = 0; i < imageHeight * imageWidth; i++) { + image[i] = std::max(0, int(image[i] - noise)); + } +} + +struct FloodParams{ + int xMin; + int xMax; + int yMin; + int yMax; +}; + +std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight) { + std::vector res; + std::set checkedPoints; + + int cutoff = BasicThreshold(image, imageWidth, imageHeight); + std::cout << "Cutoff: " << cutoff << std::endl; + if(cutoff == 0){ + // TODO: scuffed, floodfill will never stop if cutoff=0, so stop it here + std::cerr << "No stars detected in image, killing process" << std::endl; + exit(EXIT_FAILURE); + } + + for (long i = 0; i < imageHeight * imageWidth; i++) { + int x = i % imageWidth; + int y = i / imageWidth; + Point pCurr{x, y}; + if (image[i] >= cutoff && checkedPoints.count(pCurr) == 0) { + // Floodfill from this point (x, y) + // Obtain coordinates (x0, y0) of pixel with largest magnitude in the floodfill + std::vector pts; + std::deque queue; + + FloodParams fp{x, x, y, y}; + + int maxMag = -1; + int x0 = x; + int y0 = y; + + int floodSize = 0; + + queue.push_back(pCurr); + while (queue.size() != 0) { + Point p = queue[0]; + queue.pop_front(); + + int px = p[0]; + int py = p[1]; + if (px < 0 || px >= imageWidth || py < 0 || py >= imageHeight) continue; + if (checkedPoints.count(p) != 0) continue; + + int mag = Get(px, py, image, imageWidth); + if (mag < cutoff) continue; + + floodSize++; + + fp.xMin = std::min(fp.xMin, px); + fp.xMax = std::max(fp.xMax, px); + fp.yMin = std::min(fp.yMin, py); + fp.yMax = std::max(fp.yMax, py); + + // Add this point to pts and checkedPoints + // We can add to checkedPoints since cutoff is global - ensure no 2 fills collide + pts.push_back(p); + checkedPoints.insert(p); + + // Update max pixel value in fill + if (mag > maxMag) { + maxMag = mag; + x0 = px; + y0 = py; + } + + // Add all 8 adjacent points to the queue + for (int dx = -1; dx <= 1; dx++) { + for (int dy = -1; dy <= 1; dy++) { + queue.push_back({px + dx, py + dy}); + } + } + } + + // Cool, our flood is done + res.push_back({x0, y0, floodSize, fp.xMax - fp.xMin, fp.yMax - fp.yMin}); + } + } + + // std::cout << "Floodfill preprocessing: found " << res.size() << " total points" << std::endl; + return res; +} + +// TODO: documentation! +struct CentroidParams { + float yCoordMagSum; + float xCoordMagSum; + long magSum; + int xMin; + int xMax; + int yMin; + int yMax; + int cutoff; + bool isValid; + std::unordered_set checkedIndices; +}; + +int Get(int x, int y, const unsigned char *image, int w) { + int ind = y * w + x; + return image[ind]; +} + +int XMarginal(int x, int y0, int nb, const unsigned char *image, int w) { + int sum = 0; + for (int j = -nb; j <= nb; j++) { + sum += Get(x, y0 + j, image, w); + } + return sum; +} + +int YMarginal(int x0, int y, int nb, const unsigned char *image, int w) { + int sum = 0; + for (int i = -nb; i <= nb; i++) { + sum += Get(x0 + i, y, image, w); + } + return sum; +} + +static float FitModel(float x, float a, float xb, float sigma) { + return a * exp(-1 * (x - xb) * (x - xb) / (2 * sigma * sigma)); +} + +static float FitModel2D(float x, float y, float a, float xb, float yb, float sigmaX, float sigmaY) { + float term1 = exp(-1 * (x - xb) * (x - xb) / (2 * sigmaX * sigmaX)); + float term2 = exp(-1 * (y - yb) * (y - yb) / (2 * sigmaY * sigmaY)); + return a * term1 * term2; +} + +// DO NO DELETE +// void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, +// float *xb, float *yb, float *sigma) { +// // a is set to max intensity value in the window +// // (xb, yb) = coordinates of pixel with max intensity value +// int max = -1; +// for (int i = -nb; i <= nb; i++) { +// for (int j = -nb; j <= nb; j++) { +// int pixelValue = Get(x0 + i, y0 + j, image, w); +// if (pixelValue > max) { +// *xb = x0 + i; +// *yb = y0 + j; +// max = pixelValue; +// } +// } +// } +// *a = max; +// // Sigma is a little complicated +// // sigma = f_whm / (2 * \sqrt{2log(2)}) +// // f_whm \approx sqrt of number of pixels with intensity > 0.5 * a +// int halfCount = 0; +// for (int i = -nb; i <= nb; i++) { +// for (int j = -nb; j <= nb; j++) { +// if (Get(x0 + i, y0 + j, image, w) > *a / 2) { +// halfCount++; +// } +// } +// } +// float fwhm = std::sqrt(halfCount); +// *sigma = fwhm / (2 * std::sqrt(2 * std::log(2))); +// } + +float FitInitialGuessSigma(int x0, int y0, int maxMag, const int nb, const unsigned char *image, + int w) { + float halfCount = 0; + for (int i = -nb; i <= nb; i++) { + for (int j = -nb; j < +nb; j++) { + if (Get(x0 + i, y0 + j, image, w) > maxMag / 2) { + halfCount++; + } + } + } + float fwhm = std::sqrt(halfCount); + return fwhm / (2 * std::sqrt(2 * std::log(2))); +} + +// Generic functor +template +struct Functor { + typedef _Scalar Scalar; + enum { InputsAtCompileTime = NX, ValuesAtCompileTime = NY }; + typedef Eigen::Matrix InputType; + typedef Eigen::Matrix ValueType; + typedef Eigen::Matrix JacobianType; + + int m_inputs; // Number of parameters in your model + int m_values; // Number of data points + + Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} + Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} + + int inputs() const { return m_inputs; } + int values() const { return m_values; } +}; + +/// Functor for 1D Least Squares Gaussian Fit +struct LSGF1DFunctor : Functor { + LSGF1DFunctor(const int nb, const int marg, Eigen::VectorXf X, const unsigned char *image, + const int w, const int x0, const int y0) + : Functor(3, 2 * nb + 1), + nb(nb), + marg(marg), + X(X), + image(image), + w(w), + x0(x0), + y0(y0) {} + + /* + Calculate residuals (error = prediction - actual) + x = parameters (a, xb, sigma) + */ + int operator()(const Eigen::VectorXf &x, Eigen::VectorXf &fvec) const { + float a = x(0); + float xb = x(1); + float sigma = x(2); + for (int i = 0; i < X.size(); i++) { + int marginal; + if (marg == 0) + marginal = XMarginal(X(i), y0, nb, image, w); + else + marginal = YMarginal(x0, X(i), nb, image, w); + fvec(i) = FitModel(X(i), a, xb, sigma) - marginal; + } + + return 0; + } + + // Calculate Jacobian + int df(const Eigen::VectorXf &x, Eigen::MatrixXf &fjac) const { + float a = x(0); + float kb = x(1); + float sK = x(2); + + int ind = 0; + + for (int i = -nb; i <= nb; i++) { + int k = (marg == 0) ? (x0 + i) : (y0 + i); + + float expK = exp(-1 * (k - kb) * (k - kb) / (2 * sK * sK)); + + fjac(ind, 0) = expK; + fjac(ind, 1) = a * expK * (k - kb) / (sK * sK); + fjac(ind, 2) = a * expK * (k - kb) * (k - kb) / std::pow(sK, 3); + + ind++; + } + + return 0; + } + + // Window size is (2*nb + 1) x (2*nb + 1) + const int nb; + // Flag for which marginal to use (0 = X, 1 = Y) + const int marg; + // Data points (set of x or y coordinates) + Eigen::VectorXf X; + const unsigned char *image; + const int w; // image width in pixels + // Center coordinates (x0, y0) of this window + const int x0, y0; +}; + +/// Functor for 2D Least Squares Gaussian Fit +struct LSGF2DFunctor : Functor { + // We now have 5 params, beta = (a, xb, yb, sigmaX, sigmaY) + // Let entire window be (np x np) pixels + // We have np^2 data points + LSGF2DFunctor(const int nb, const unsigned char *image, const int w, const int x0, const int y0) + : Functor(5, (2 * nb + 1) * (2 * nb + 1)), + nb(nb), + image(image), + w(w), + x0(x0), + y0(y0) {} + + /* + x = parameters (a, xb, yb, sigmaX, sigmaY) + */ + int operator()(const Eigen::VectorXf &x, Eigen::VectorXf &fvec) const { + float a = x(0); + float xb = x(1); + float yb = x(2); + float sigmaX = x(3); + float sigmaY = x(4); + + int ind = 0; + for (int i = -nb; i <= nb; i++) { + for (int j = -nb; j <= nb; j++) { + int xi = x0 + i; + int yi = y0 + j; + float yPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); + float yActual = Get(xi, yi, image, w); + fvec(ind) = yPred - yActual; + + ind++; + } + } + + return 0; + } + + int df(const Eigen::VectorXf &x, Eigen::MatrixXf &fjac) const { + float a = x(0); + float xb = x(1); + float yb = x(2); + float sigmaX = x(3); + float sigmaY = x(4); + + int ind = 0; + + for (int i = -nb; i <= nb; i++) { + for (int j = -nb; j <= nb; j++) { + int xi = x0 + i; + int yj = y0 + j; + + float expX = exp(-1 * (xi - xb) * (xi - xb) / (2 * sigmaX * sigmaX)); + float expY = exp(-1 * (yj - yb) * (yj - yb) / (2 * sigmaY * sigmaY)); + + fjac(ind, 0) = expX * expY; + fjac(ind, 1) = a * expX * expY * (xi - xb) / (sigmaX * sigmaX); + fjac(ind, 2) = a * expX * expY * (yj - yb) / (sigmaY * sigmaY); + fjac(ind, 3) = a * expX * expY * std::pow(xi - xb, 2) / (std::pow(sigmaX, 3)); + fjac(ind, 4) = a * expX * expY * std::pow(yj - yb, 2) / (std::pow(sigmaY, 3)); + + ind++; + } + } + + return 0; + } + + // Window size is (2*nb + 1) x (2*nb + 1) + const int nb; + + const unsigned char *image; + const int w; // image width in pixels + // Center coordinates (x0, y0) of this window + const int x0, y0; +}; + +std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + std::vector result; + + std::cout << "1D Gaussian Fit (" << np << "x" << np << ")" << std::endl; + + SubtractNoise(image, imageWidth, imageHeight); + + std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); + + for (const Point &pt : candidatePts) { + int x = pt[0]; + int y = pt[1]; + // TODO: just taking largest diameter right now + int dynamicWinSize = std::max(pt[3], pt[4]); + + int winRadius = (dynamic) ? dynamicWinSize : nb; // /2? + // int winRadius = nb; + + if (x - winRadius < 0 || x + winRadius >= imageWidth || y - winRadius < 0 || + y + winRadius >= imageHeight) + continue; + + Eigen::VectorXf X(2 * winRadius + 1); + Eigen::VectorXf Y(2 * winRadius + 1); + int vInd = 0; + for (int i = -winRadius; i <= winRadius; i++) { + X(vInd) = x + i; + Y(vInd) = y + i; + vInd++; + } + + float a = Get(x, y, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, winRadius, image, imageWidth); + + Eigen::VectorXf betaX(3); + betaX << a, x, sigma; + + Eigen::VectorXf betaY(3); + betaY << a, y, sigma; + + LSGF1DFunctor functorX(winRadius, 0, X, image, imageWidth, x, y); + Eigen::LevenbergMarquardt lmX(functorX); + + lmX.parameters.maxfev = 2000; + lmX.parameters.xtol = 1.0e-10; + + LSGF1DFunctor functorY(winRadius, 1, Y, image, imageWidth, x, y); + Eigen::LevenbergMarquardt lmY(functorY); + + lmY.parameters.maxfev = 2000; + lmY.parameters.xtol = 1.0e-10; + + lmX.minimize(betaX); + lmY.minimize(betaY); + + a = betaX(0); + + float xb = betaX(1); + float yb = betaY(1); + + sigma = betaX(2); + + result.push_back(Star(xb + 0.5, yb + 0.5, sigma, sigma, pt[2])); + } + + return result; +} + +std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + std::vector result; + + std::cout << "2D Gaussian Fit (" << np << "x" << np << ")" << std::endl; + + SubtractNoise(image, imageWidth, imageHeight); + + std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); + + for (const Point &pt : candidatePts) { + int x = pt[0]; + int y = pt[1]; + + // TODO: just taking largest diameter right now + int dynamicWinSize = std::max(pt[3], pt[4]); + + int winRadius = (dynamic) ? dynamicWinSize : nb; // /2? + // int winRadius = nb; + + if (x - winRadius < 0 || x + winRadius >= imageWidth || y - winRadius < 0 || + y + winRadius >= imageHeight) + continue; + + float a = Get(x, y, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, winRadius, image, imageWidth); + + Eigen::VectorXf beta(5); + beta << a, x, y, sigma, sigma; + + LSGF2DFunctor functor(winRadius, image, imageWidth, x, y); + Eigen::LevenbergMarquardt lm(functor); + + lm.parameters.maxfev = 2000; + lm.parameters.xtol = 1.0e-10; + + lm.minimize(beta); + + float aRes = beta(0); + float xRes = beta(1); + float yRes = beta(2); + float sigmaX = beta(3); + float sigmaY = beta(4); -std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, int imageHeight) const { + result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, pt[2])); + } + + return result; +} + +static void GetGGridCoeffs(const std::vector &w, std::vector *const a, + std::vector *const b) { + (*a)[0] = 2 * w[2] * w[1] * w[0] + 6 * w[2] * w[0] * w[3] + + 12 * (w[3] * w[4] * w[0] + w[3] * w[1] * w[0]) + 32 * w[2] * w[4] * w[0] + + 36 * w[4] * w[1] * w[0]; + (*a)[1] = 2 * w[2] * w[3] * w[1] + 6 * w[3] * w[4] * w[1] - 4 * w[2] * w[1] * w[0] + + 12 * w[2] * w[4] * w[1] - 18 * w[3] * w[1] * w[0] - 48 * w[4] * w[1] * w[0]; + (*a)[2] = 2 * (w[2] * w[3] * w[4] - w[2] * w[1] * w[0]) - 4 * w[1] * w[2] * w[3] - + 18 * (w[0] * w[2] * w[3] + w[1] * w[2] * w[4]) - 64 * w[0] * w[2] * w[4]; + (*a)[3] = 2 * w[2] * w[1] * w[3] + 6 * w[3] * w[1] * w[0] - 4 * w[2] * w[3] * w[4] + + 12 * w[2] * w[3] * w[0] - 18 * w[3] * w[4] * w[1] - 48 * w[3] * w[4] * w[0]; + (*a)[4] = 2 * w[2] * w[3] * w[4] + 6 * w[2] * w[4] * w[1] + + 12 * (w[3] * w[4] * w[1] + w[4] * w[1] * w[0]) + 32 * w[2] * w[4] * w[0] + + 36 * w[3] * w[4] * w[0]; + + (*b)[0] = -w[2] * w[1] * w[0] + 3 * w[2] * w[3] * w[0] + + 18 * (w[3] * w[4] * w[0] + w[4] * w[1] * w[0]) + 32 * w[2] * w[4] * w[0]; + (*b)[1] = w[2] * w[3] * w[1] + 4 * w[2] * w[1] * w[0] + + 9 * (w[3] * w[4] * w[1] + w[3] * w[1] * w[0]) + 12 * w[2] * w[4] * w[1]; + (*b)[2] = 3 * (w[2] * w[3] * w[4] - w[2] * w[1] * w[0]) + + 9 * (w[2] * w[3] * w[0] - w[2] * w[1] * w[4]); + (*b)[3] = -w[2] * w[3] * w[1] - 4 * w[2] * w[3] * w[4] - + 9 * (w[3] * w[4] * w[1] + w[3] * w[0] * w[1]) - 12 * w[2] * w[3] * w[0]; + (*b)[4] = w[2] * w[3] * w[4] - 3 * w[2] * w[1] * w[4] - + 18 * (w[3] * w[0] * w[4] + w[1] * w[0] * w[4]) - 32 * w[2] * w[4] * w[0]; +} + +Stars GaussianGrid::Go(unsigned char *image, int imageWidth, int imageHeight) const { + std::vector result; + + std::cout << "Gaussian Grid (" << np << "x" << np << ")" << std::endl; + + SubtractNoise(image, imageWidth, imageHeight); + + std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); + + for (const Point &pt : candidatePts) { + int x = pt[0]; + int y = pt[1]; + if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; + + float nom = 0; + float denom = 0; + + for (int j = -nb; j <= nb; j++) { + std::vector w; + for (int k = -nb; k <= nb; k++) { + w.push_back(std::pow(Get(x + k, y + j, image, imageWidth), 1)); + } + std::vector a(np); + std::vector b(np); + GetGGridCoeffs(w, &a, &b); + for (int i = -nb; i <= nb; i++) { + float v = Get(x + i, y + j, image, imageWidth); + nom += b[i + nb] * ((v == 0) ? -1e6 : std::log(v)); + denom += a[i + nb] * ((v == 0) ? -1e6 : std::log(v)); + } + } + + float xRes = x + nom / denom; + + nom = 0; + denom = 0; + for (int i = -nb; i <= nb; i++) { + std::vector w; + for (int k = -nb; k <= nb; k++) { + w.push_back(Get(x + i, y + k, image, imageWidth)); + } + std::vector a(np); + std::vector b(np); + GetGGridCoeffs(w, &a, &b); + for (int j = -nb; j <= nb; j++) { + float v = Get(x + i, y + j, image, imageWidth); + nom += b[j + nb] * ((v == 0) ? -1e6 : std::log(v)); + denom += a[j + nb] * ((v == 0) ? -1e6 : std::log(v)); + } + } + + float yRes = y + nom / denom; + + result.push_back(Star(xRes + 0.5, yRes + 0.5, nb, nb, pt[2])); + } + + return result; +} + +// DUMMY +std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, + int imageHeight) const { std::vector result; 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, 10.0)); } return result; @@ -27,19 +616,19 @@ std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, in // a poorly designed thresholding algorithm int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) { - //loop through entire array, find sum of magnitudes + // loop through entire array, find sum of magnitudes long totalMag = 0; for (long i = 0; i < imageHeight * imageWidth; i++) { totalMag += image[i]; } - return (((totalMag/(imageHeight * imageWidth)) + 1) * 15) / 10; + return (((totalMag / (imageHeight * imageWidth)) + 1) * 15) / 10; } // a more sophisticated thresholding algorithm, not tailored to star images int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { // code here, duh long total = imageWidth * imageHeight; - //float top = 255; + // float top = 255; float sumB = 0; float sum1 = 0; float wB = 0; @@ -48,22 +637,22 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { // make the histogram (array length 256) int histogram[256]; - memset(histogram, 0, sizeof(int)*256); + memset(histogram, 0, sizeof(int) * 256); for (long i = 0; i < total; i++) { histogram[image[i]]++; } - for (int i = 0; i < 256; i ++) { + for (int i = 0; i < 256; i++) { sum1 += i * histogram[i]; } - for (int i = 0; i < 256; i ++) { + for (int i = 0; i < 256; i++) { float wF = total - wB; - //std::cout << "wF\n" << wB << "\n"; - //std::cout << "wB\n" << wF << "\n"; + // 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); - //std::cout << val << "\n"; + // std::cout << val << "\n"; if (val >= maximum) { level = i; maximum = val; @@ -75,7 +664,6 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { return level; } -// 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; @@ -107,25 +695,13 @@ int BasicThresholdOnePass(unsigned char *image, int imageWidth, int imageHeight) return mean + (std * 5); } -struct CentroidParams { - float yCoordMagSum; - float xCoordMagSum; - long magSum; - int xMin; - int xMax; - int yMin; - int yMax; - int cutoff; - bool isValid; - std::unordered_set checkedIndices; -}; - -//recursive helper here +// recursive helper here void CogHelper(CentroidParams *p, long i, unsigned char *image, int imageWidth, int imageHeight) { - - if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && p->checkedIndices.count(i) == 0) { - //check if pixel is on the edge of the image, if it is, we dont want to centroid this star - if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || i / imageWidth == imageHeight - 1) { + if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && + p->checkedIndices.count(i) == 0) { + // check if pixel is on the edge of the image, if it is, we dont want to centroid this star + if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || + i / imageWidth == imageHeight - 1) { p->isValid = false; } p->checkedIndices.insert(i); @@ -153,21 +729,23 @@ void CogHelper(CentroidParams *p, long i, unsigned char *image, int imageWidth, } } -std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { +std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, + int imageHeight) const { CentroidParams p; std::vector result; + SubtractNoise(image, imageWidth, imageHeight); + p.cutoff = BasicThreshold(image, imageWidth, imageHeight); for (long i = 0; i < imageHeight * imageWidth; i++) { if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { - - //iterate over pixels that are part of the star - int xDiameter = 0; //radius of current star + // iterate over pixels that are part of the star + int xDiameter = 0; // radius of current star int yDiameter = 0; - p.yCoordMagSum = 0; //y coordinate of current star - p.xCoordMagSum = 0; //x coordinate of current star - p.magSum = 0; //sum of magnitudes of current star + p.yCoordMagSum = 0; // y coordinate of current star + p.xCoordMagSum = 0; // x coordinate of current star + p.magSum = 0; // sum of magnitudes of current star p.xMax = i % imageWidth; p.xMin = i % imageWidth; @@ -181,20 +759,22 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi xDiameter = (p.xMax - p.xMin) + 1; yDiameter = (p.yMax - p.yMin) + 1; - //use the sums to finish CoG equation and add stars to the result + // 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)); 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 + 0.5f, yCoord + 0.5f, ((float)(xDiameter)) / 2.0f, + ((float)(yDiameter)) / 2.0f, + p.checkedIndices.size() - sizeBefore)); } } } return result; } -//Determines how accurate and how much iteration is done by the IWCoG algorithm, -//smaller means more accurate and more iterations. +// Determines how accurate and how much iteration is done by the IWCoG algorithm, +// smaller means more accurate and more iterations. float iWCoGMinChange = 0.0002; struct IWCoGParams { @@ -209,10 +789,13 @@ struct IWCoGParams { std::unordered_set checkedIndices; }; -void IWCoGHelper(IWCoGParams *p, long i, unsigned char *image, int imageWidth, int imageHeight, std::vector *starIndices) { - if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && p->checkedIndices.count(i) == 0) { - //check if pixel is on the edge of the image, if it is, we dont want to centroid this star - if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || i / imageWidth == imageHeight - 1) { +void IWCoGHelper(IWCoGParams *p, long i, unsigned char *image, int imageWidth, int imageHeight, + std::vector *starIndices) { + if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && + p->checkedIndices.count(i) == 0) { + // check if pixel is on the edge of the image, if it is, we dont want to centroid this star + if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || + i / imageWidth == imageHeight - 1) { p->isValid = false; } p->checkedIndices.insert(i); @@ -242,24 +825,26 @@ void IWCoGHelper(IWCoGParams *p, long i, unsigned char *image, int imageWidth, i } } -Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { +Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, + int imageHeight) const { IWCoGParams p; std::vector result; + SubtractNoise(image, imageWidth, imageHeight); p.cutoff = BasicThreshold(image, imageWidth, imageHeight); for (long i = 0; i < imageHeight * imageWidth; i++) { - //check if pixel is part of a "star" and has not been iterated over + // check if pixel is part of a "star" and has not been iterated over if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { // TODO: store longs --Mark - std::vector starIndices; //indices of the current star + std::vector starIndices; // indices of the current star p.maxIntensity = 0; int xDiameter = 0; int yDiameter = 0; float yWeightedCoordMagSum = 0; float xWeightedCoordMagSum = 0; float weightedMagSum = 0; - float fwhm; //fwhm variable + float fwhm; // fwhm variable float standardDeviation; - float w; //weight value + float w; // weight value p.xMax = i % imageWidth; p.xMin = i % imageWidth; @@ -267,15 +852,14 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im p.yMin = i / imageWidth; p.isValid = true; - IWCoGHelper(&p, i, image, imageWidth, imageHeight, &starIndices); xDiameter = (p.xMax - p.xMin) + 1; yDiameter = (p.yMax - p.yMin) + 1; - //calculate fwhm + // calculate fwhm float count = 0; - for (int j = 0; j < (int) starIndices.size(); j++) { + for (int j = 0; j < (int)starIndices.size(); j++) { if (image[starIndices.at(j)] > p.maxIntensity / 2) { count++; } @@ -284,27 +868,30 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im 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); - //how much our new centroid estimate changes w each iteration + float guessXCoord = (float)(p.guess % imageWidth); + float guessYCoord = (float)(p.guess / imageWidth); + // how much our new centroid estimate changes w each iteration float change = INFINITY; int stop = 0; - //while we see some large enough change in estimated, maybe make it a global variable + // while we see some large enough change in estimated, maybe make it a global variable while (change > iWCoGMinChange && stop < 100000) { - //traverse through star indices, calculate W at each coordinate, add to final coordinate sums + // traverse through star indices, calculate W at each coordinate, add to final + // coordinate sums yWeightedCoordMagSum = 0; xWeightedCoordMagSum = 0; weightedMagSum = 0; 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))); - - xWeightedCoordMagSum += w * currXCoord * ((float) image[starIndices.at(j)]); - yWeightedCoordMagSum += w * currYCoord * ((float) image[starIndices.at(j)]); - weightedMagSum += w * ((float) image[starIndices.at(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))); + + xWeightedCoordMagSum += w * currXCoord * ((float)image[starIndices.at(j)]); + yWeightedCoordMagSum += w * currYCoord * ((float)image[starIndices.at(j)]); + weightedMagSum += w * ((float)image[starIndices.at(j)]); } float xTemp = xWeightedCoordMagSum / weightedMagSum; float yTemp = yWeightedCoordMagSum / weightedMagSum; @@ -315,11 +902,13 @@ 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 + 0.5f, guessYCoord + 0.5f, + ((float)(xDiameter)) / 2.0f, ((float)(yDiameter)) / 2.0f, + starIndices.size())); } } } return result; } -} +} // namespace lost diff --git a/src/centroiders.hpp b/src/centroiders.hpp index d18683f7..ec7226a9 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -21,6 +21,42 @@ class CentroidAlgorithm { virtual ~CentroidAlgorithm() { }; }; +/** + * @brief Least Squares Gaussian Fit (1D) centroiding algorithm + * + * Detect centroids by fitting a 1D Gaussian to the marginals of each window + * Slightly less accurate than 2D fit + * Note: increasing window size may introduce noise, making this method less accurate on small stars + */ +class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ +public: + explicit LeastSquaresGaussianFit1D(int nb, bool dyn) : nb(nb), dynamic(dyn), np(nb*2+1) { }; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + +private: + const int nb; + bool dynamic; + const int np; +}; + +/** + * @brief Least Squares Gaussian Fit (2D) centroiding algorithm + * + * Detect centroids by fitting a 2D Gaussian to all pixels in each window + * This is the most accurate centroiding function to date + * Also more computationally expensive than other methods, scales exponentially with increasing window size + */ +class LeastSquaresGaussianFit2D : public CentroidAlgorithm { +public: + explicit LeastSquaresGaussianFit2D(int nb, bool dyn) : nb(nb), dynamic(dyn), np(nb*2+1) { }; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + +private: + const int nb; + bool dynamic; + const int np; +}; + /// A centroid algorithm for debugging that returns random centroids. class DummyCentroidAlgorithm: public CentroidAlgorithm { public: @@ -45,11 +81,53 @@ class CenterOfGravityAlgorithm : public CentroidAlgorithm { * Iteratively estimates the center of the centroid. Some papers report that it is slightly more precise than CenterOfGravityAlgorithm, but that has not been our experience. */ class IterativeWeightedCenterOfGravityAlgorithm : public CentroidAlgorithm { - public: - IterativeWeightedCenterOfGravityAlgorithm() { }; - Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; +public: + IterativeWeightedCenterOfGravityAlgorithm() { }; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; +}; + +class GaussianGrid : public CentroidAlgorithm { +public: + GaussianGrid(){ }; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + +private: + // Only support 5x5 Gaussian Grid for now + const int nb = 2; + const int np = nb*2 + 1; }; -} +/// Get value of pixel at (x, y) in image with width=w +int Get(int x, int y, const unsigned char *image, int w); + +/// Get value of x marginal - keep x fixed, sum pixel values from [y0-nb, y0+nb] +int XMarginal(int x, int y0, int nb, const unsigned char *image, int w); + +/// Get value of y marginal - keep y fixed, sum pixel values from [x0-nb, x0+nb] +int YMarginal(int x0, int y, int nb, const unsigned char *image, int w); + +/// Compute and subtract noise = mean of all pixels from image +void SubtractNoise(unsigned char *image, int imageWidth, int imageHeight); + +/** + * @brief Given window centered at (x0, y0) in image with width=w, calculate initial guess for sigma + * parameter of Gaussian function + * Only used with Least Squares Gaussian Fit algorithms + */ +float FitInitialGuessSigma(int x0, int y0, int maxMag, const int nb, const unsigned char *image, + int w); + +// DO NOT DELETE +/* +Get initial guess for window centered at (x0, y0) with given nb +Output: +a = max intensity value +(xb, yb) = coordinates of pixel with max intensity +sigma = standard deviation (sigmaX = sigmaY) +*/ +// void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, +// float *xb, float *yb, double *sigma); + +} // namespace lost #endif diff --git a/src/io.cpp b/src/io.cpp index 664b04e0..9ae9fc9b 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -411,7 +411,7 @@ class GeneratedStar : public Star { GeneratedStar(Star star, float peakBrightness, Vec2 motionBlurDelta) : Star(star), peakBrightness(peakBrightness), delta(motionBlurDelta) { }; - /// the brightness density at the center of the star. 0.0 is black, 1.0 is white. + /// the brightness density per time unit at the center of the star. 0.0 is black, 1.0 is white. float peakBrightness; /// (only meaningful with motion blur) Where the star will appear one time unit in the future. @@ -448,9 +448,9 @@ static float MotionBlurredPixelBrightness(const Vec2 &pixel, const GeneratedStar /// Like motionBlurredPixelBrightness, but for when motion blur is disabled. static float StaticPixelBrightness(const Vec2 &pixel, const GeneratedStar &generatedStar, - float stddev) { + float t, float stddev) { const Vec2 d0 = generatedStar.position - pixel; - return generatedStar.peakBrightness * exp(-d0.MagnitudeSq() / (2*stddev*stddev)); + return generatedStar.peakBrightness * t * exp(-d0.MagnitudeSq() / (2 * stddev * stddev)); } /** @@ -482,9 +482,9 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, std::default_random_engine *rng, bool centroidsOnly, - float observedReferenceBrightness, + float zeroMagTotalPhotons, float starSpreadStdDev, - float sensitivity, + float saturationPhotons, float darkCurrent, float readNoiseStdDev, Attitude motionBlurDirection, // applied on top of the attitude @@ -502,9 +502,6 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, assert(falseStarMaxMagnitude <= falseStarMinMagnitude); assert(perturbationStddev >= 0.0); - // in photons - float referenceBrightness = observedReferenceBrightness / sensitivity; - image.width = camera.XResolution(); image.height = camera.YResolution(); // number of true photons each pixel receives. @@ -515,11 +512,8 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, std::cerr << "WARNING: oversampling was not a perfect square. Rounding up to " << oversamplingPerAxis*oversamplingPerAxis << "." << std::endl; } - bool motionBlurEnabled = exposureTime > 0 && abs(motionBlurDirection.GetQuaternion().Angle()) > 0.001; - if (!motionBlurEnabled) { - exposureTime = 1.0; - motionBlurDirection = Attitude(Quaternion(0, 1, 0, 0)); - } + assert(exposureTime > 0); + bool motionBlurEnabled = abs(motionBlurDirection.GetQuaternion().Angle()) > 0.001; Quaternion motionBlurDirectionQ = motionBlurDirection.GetQuaternion(); // attitude at the middle of exposure time Quaternion currentAttitude = attitude.GetQuaternion(); @@ -527,6 +521,10 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, Quaternion futureAttitude = motionBlurDirectionQ*currentAttitude; std::vector generatedStars; + // 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); + // TODO: Is it 100% correct to just copy the standard deviation in both dimensions? std::normal_distribution perturbation1DDistribution(0.0, perturbationStddev); @@ -560,19 +558,19 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, delta = {0, 0}; // avoid floating point funny business } // radiant intensity, in photons per time unit per pixel, at the center of the star. - float peakBrightness = referenceBrightness * MagToBrightness(catalogStar.magnitude); - float interestingThreshold = 0.01; // we don't need to check pixels that are expected to + float peakBrightnessPerTime = zeroMagPeakPhotonDensity * MagToBrightness(catalogStar.magnitude); + float interestingThreshold = 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/peakBrightness/exposureTime)*2*starSpreadStdDev*starSpreadStdDev)); + float radius = ceil(sqrt(-log(interestingThreshold/peakBrightnessPerTime/exposureTime)*2*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. // It's possible to make it so that the magnitude is always positive too, but allowing weirder magnitudes helps keep star-id algos honest about their assumptions on magnitude. // we don't use its magnitude anywhere else in generation; peakBrightness was already calculated. -catalogStar.magnitude); - generatedStars.push_back(GeneratedStar(star, peakBrightness, delta)); + generatedStars.push_back(GeneratedStar(star, peakBrightnessPerTime, delta)); // Now add the star to the input and expected lists. // We do actually want to add false stars as well, because: @@ -651,7 +649,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, - MotionBlurredPixelBrightness({x, y}, star, tStart, starSpreadStdDev)) / oversamplingBrightnessFactor; } else { - curPhotons = StaticPixelBrightness({x, y}, star, starSpreadStdDev) + curPhotons = StaticPixelBrightness({x, y}, star, exposureTime, starSpreadStdDev) / oversamplingBrightnessFactor; } @@ -695,7 +693,7 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, } else { quantizedPhotons = round(photonsBuffer[i]); } - curBrightness += quantizedPhotons * sensitivity; + 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); @@ -772,9 +770,9 @@ PipelineInputList GetGeneratedPipelineInput(const PipelineOptions &values) { &rng, values.generateCentroidsOnly, - values.generateRefBrightness, + values.generateZeroMagPhotons, values.generateSpreadStdDev, - values.generateSensitivity, + values.generateSaturationPhotons, values.generateDarkCurrent, values.generateReadNoiseStdDev, motionBlurDirection, @@ -836,10 +834,6 @@ Pipeline::Pipeline(CentroidAlgorithm *centroidAlgorithm, /// Create a pipeline from command line options. Pipeline SetPipeline(const PipelineOptions &values) { - enum class PipelineStage { - Centroid, CentroidMagnitudeFilter, Database, StarId, AttitudeEstimation, Done - }; - Pipeline result; // TODO: more flexible or sth @@ -848,18 +842,31 @@ Pipeline SetPipeline(const PipelineOptions &values) { // centroid algorithm stage if (values.centroidAlgo == "dummy") { - result.centroidAlgorithm = std::unique_ptr(new DummyCentroidAlgorithm(values.centroidDummyNumStars)); + result.centroidAlgorithm = std::unique_ptr( + new DummyCentroidAlgorithm(values.centroidDummyNumStars)); } else if (values.centroidAlgo == "cog") { - result.centroidAlgorithm = std::unique_ptr(new CenterOfGravityAlgorithm()); + result.centroidAlgorithm = + std::unique_ptr(new CenterOfGravityAlgorithm()); } else if (values.centroidAlgo == "iwcog") { - result.centroidAlgorithm = std::unique_ptr(new IterativeWeightedCenterOfGravityAlgorithm()); + result.centroidAlgorithm = + std::unique_ptr(new IterativeWeightedCenterOfGravityAlgorithm()); + } else if (values.centroidAlgo == "lsgf1d") { + result.centroidAlgorithm = + std::unique_ptr(new LeastSquaresGaussianFit1D(values.centroidFitRadius, values.centroidFitDynamic)); + } else if (values.centroidAlgo == "lsgf2d") { + result.centroidAlgorithm = + std::unique_ptr(new LeastSquaresGaussianFit2D(values.centroidFitRadius, values.centroidFitDynamic)); + } else if (values.centroidAlgo == "ggrid"){ + result.centroidAlgorithm = + std::unique_ptr(new GaussianGrid()); } else if (values.centroidAlgo != "") { std::cout << "Illegal centroid algorithm." << std::endl; exit(1); } // centroid magnitude filter stage - if (values.centroidMagFilter != -1) result.centroidMinMagnitude = values.centroidMagFilter; + if (values.centroidMagFilter > 0) result.centroidMinMagnitude = values.centroidMagFilter; + if (values.centroidFilterBrightest > 0) result.centroidMinStars = values.centroidFilterBrightest; // database stage if (values.databasePath != "") { @@ -939,9 +946,23 @@ PipelineOutput Pipeline::Go(const PipelineInput &input) { std::chrono::time_point end = std::chrono::steady_clock::now(); result.centroidingTimeNs = std::chrono::duration_cast(end - start).count(); + // MAGNITUDE FILTERING + int minMagnitude = centroidMinMagnitude; + if (centroidMinStars > 0 + // don't need to filter if we don't even have that many stars + && centroidMinStars < (int)unfilteredStars.size()) { + + Stars magSortedStars = unfilteredStars; + // sort descending + std::sort(magSortedStars.begin(), magSortedStars.end(), [](const Star &a, const Star &b) { return a.magnitude > b.magnitude; }); + minMagnitude = std::max(minMagnitude, magSortedStars[centroidMinStars - 1].magnitude); + } + // determine the minimum magnitude according to sorted stars Stars *filteredStars = new std::vector(); for (const Star &star : unfilteredStars) { - if (star.magnitude >= centroidMinMagnitude) { + assert(star.magnitude >= 0); // catalog stars can have negative magnitude, but by our + // conventions, centroids shouldn't. + if (star.magnitude >= minMagnitude) { filteredStars->push_back(star); } } @@ -1097,11 +1118,13 @@ CentroidComparison CentroidComparisonsCombine(std::vector co for (const CentroidComparison &comparison : comparisons) { result.meanError += comparison.meanError; + result.numCorrectCentroids += comparison.numCorrectCentroids; result.numExtraCentroids += comparison.numExtraCentroids; } result.meanError /= comparisons.size(); result.numExtraCentroids /= comparisons.size(); + result.numCorrectCentroids /= comparisons.size(); return result; } @@ -1474,7 +1497,12 @@ static void PrintTimeStats(std::ostream &os, const std::string &prefix, const st long average = sum / times.size(); std::vector sortedTimes = times; std::sort(sortedTimes.begin(), sortedTimes.end()); - long ninetyFifthPercentile = sortedTimes[(int)(sortedTimes.size() * 0.95)]; + // what really is the 95th percentile? Being conservative, we want to pick a value that at least + // 95% of the times are less than. This means: (1) finding the number of times, (2) Finding + // Math.ceil(0.95 * numTimes), and (3) subtracting 1 to get the index. + int ninetyFiveIndex = (int)std::ceil(0.95 * times.size()) - 1; + assert(ninetyFiveIndex >= 0); + long ninetyFifthPercentile = sortedTimes[ninetyFiveIndex]; os << prefix << "_average_ns " << average << std::endl; os << prefix << "_min_ns " << min << std::endl; @@ -1490,20 +1518,33 @@ static void PipelineComparatorPrintSpeed(std::ostream &os, std::vector centroidingTimes; std::vector starIdTimes; std::vector attitudeTimes; + std::vector totalTimes; for (int i = 0; i < (int)actual.size(); i++) { - centroidingTimes.push_back(actual[i].centroidingTimeNs); - starIdTimes.push_back(actual[i].starIdTimeNs); - attitudeTimes.push_back(actual[i].attitudeEstimationTimeNs); + long totalTime = 0; + if (actual[i].centroidingTimeNs > 0) { + centroidingTimes.push_back(actual[i].centroidingTimeNs); + totalTime += actual[i].centroidingTimeNs; + } + if (actual[i].starIdTimeNs > 0) { + starIdTimes.push_back(actual[i].starIdTimeNs); + totalTime += actual[i].starIdTimeNs; + } + if (actual[i].attitudeEstimationTimeNs > 0) { + attitudeTimes.push_back(actual[i].attitudeEstimationTimeNs); + totalTime += actual[i].attitudeEstimationTimeNs; + } + totalTimes.push_back(totalTime); } - if (centroidingTimes[0] > 0) { + if (centroidingTimes.size() > 0) { PrintTimeStats(os, "centroiding", centroidingTimes); } - if (starIdTimes[0] > 0) { + if (starIdTimes.size() > 0) { PrintTimeStats(os, "starid", starIdTimes); } - if (attitudeTimes[0] > 0) { + if (attitudeTimes.size() > 0) { PrintTimeStats(os, "attitude", attitudeTimes); } + PrintTimeStats(os, "total", totalTimes); } // TODO: add these debug comparators back in! @@ -1571,7 +1612,7 @@ void PipelineComparison(const PipelineInputList &expected, const std::vector &actual, const PipelineOptions &values) { if (actual.size() == 0) { - std::cerr << "ERROR: No \"comparator\"/output action was specified. I.e., the star identification is all done, but you didn't specify how to return the results to you! Try --plot-output , perhaps." << std::endl; + std::cerr << "ERROR: No output! Did you specify any input images? Try --png or --generate." << std::endl; exit(1); } diff --git a/src/io.hpp b/src/io.hpp index 03072f99..9ce422ce 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -240,7 +240,11 @@ class Pipeline { private: std::unique_ptr centroidAlgorithm; + + // next two options are for magnitude filter: int centroidMinMagnitude = 0; + int centroidMinStars = 0; + std::unique_ptr starIdAlgorithm; std::unique_ptr attitudeEstimationAlgorithm; std::unique_ptr database; diff --git a/src/pipeline-options.hpp b/src/pipeline-options.hpp index 23e4763c..6c76f682 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -21,15 +21,18 @@ LOST_CLI_OPTION("pixel-size" , float , pixelSize , -1 , atof(optarg) , LOST_CLI_OPTION("fov" , float , 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("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-fit-radius" , int , centroidFitRadius , 2 , atoi(optarg) , kNoDefaultArgument) +LOST_CLI_OPTION("centroid-fit-dynamic" , bool , centroidFitDynamic , false , atobool(optarg) , kNoDefaultArgument) +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") // OUTPUT COMPARISON LOST_CLI_OPTION("centroid-compare-threshold", float , centroidCompareThreshold, 2 , atof(optarg), kNoDefaultArgument) @@ -48,30 +51,30 @@ LOST_CLI_OPTION("compare-star-ids" , std::string, compareStarIds 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-reference-brightness", float , generateRefBrightness , 100 , 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 , 4 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-blur-de" , float , generateBlurDe , 1 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-blur-roll" , float , generateBlurRoll , 20 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-sensitivity" , float , generateSensitivity , 0.01 , atof(optarg) , kNoDefaultArgument) -LOST_CLI_OPTION("generate-exposure" , float , generateExposure , 0 , 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" , 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) diff --git a/test/gaussian-centroid.cpp.off b/test/gaussian-centroid.cpp.off new file mode 100644 index 00000000..a91cdb81 --- /dev/null +++ b/test/gaussian-centroid.cpp.off @@ -0,0 +1,78 @@ +#include +#include + +#include "centroiders.hpp" + +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); + +// Vec3 spatial = camera.CameraToSpatial({expectedX, expectedY}); +// Vec2 actualPixels = camera.SpatialToCamera(spatial); + +// CHECK(actualPixels.x == Approx(expectedX).margin(0.000001)); +// CHECK(actualPixels.y == Approx(expectedY).margin(0.000001)); +// } + +TEST_CASE("Get pixel value from image", "[img]"){ + const int w = 5; + const int h = 5; + unsigned char image[w*h] = {0, 10, 30, 40, 10, + 10, 10, 50, 255, 255, + 0, 50, 255, 255, 255, + 40, 255, 255, 255, 255, + 35, 255, 255, 255, 255}; + + CHECK(Get(3, 0, image, w) == 40); + CHECK(Get(0, 4, image, w) == 35); + CHECK(Get(4, 4, image, w) == 255); +} + +TEST_CASE("Calculate X-Marginal", "[lsgf1]"){ + const int w = 5; + const int h = 5; + unsigned char image[w*h] = {0, 10, 30, 40, 10, + 10, 10, 50, 255, 255, + 0, 50, 255, 255, 255, + 40, 255, 255, 255, 255, + 35, 255, 255, 255, 255}; + + CHECK(XMarginal(2, 1, 1, image, w) == 30 + 50 + 255); + CHECK(XMarginal(2, 2, 2, image, w) == 30 + 50 + 255 * 3); + CHECK(XMarginal(3, 3, 1, image, w) == 255 * 3); +} + +TEST_CASE("Calculate Y-Marginal", "[lsgf1]"){ + const int w = 5; + const int h = 5; + unsigned char image[w*h] = {0, 10, 30, 40, 10, + 10, 10, 50, 255, 255, + 0, 50, 255, 255, 255, + 40, 255, 255, 255, 255, + 35, 255, 255, 255, 255}; + + CHECK(YMarginal(2, 1, 2, image, w) == 10 + 10 + 50 + 255 + 255); + CHECK(YMarginal(1, 4, 1, image, w) == 35 + 255 + 255); +} + +TEST_CASE("Calculate initial guess for Gaussian Fit params", "[lsgf]"){ + const int w = 5; + const int h = 5; + /////////////////////////// 0 1 2 3 4 + unsigned char image[w*h] = {0, 10, 30, 40, 10, // 0 + 10, 10, 50, 255, 255, // 1 + 0, 50, 255, 255, 255, // 2 + 40, 255, 255, 255, 255, // 3 + 35, 255, 255, 255, 255}; // 4 + + float a; + float xb, yb; + double sigma; + InitialGuess(2, 1, 1, image, w, &a, &xb, &yb, &sigma); + CHECK(a == 255); + std::cout << "xb, yb = " << xb << ", " << yb << std::endl; +} \ No newline at end of file diff --git a/test/kvector.off b/test/kvector.off new file mode 100644 index 00000000..7b02f54a --- /dev/null +++ b/test/kvector.off @@ -0,0 +1,111 @@ +#include + +#include "databases.hpp" +#include "io.hpp" +#include "attitude-utils.hpp" + +using namespace lost; // NOLINT + +static unsigned char *BuildPairDistanceKVectorDatabase( + const Catalog &catalog, long *length, float minDistance, float maxDistance, long numBins) { + + long dummyLength; + if (length == NULL) + length = &dummyLength; + + *length = SerializeLengthPairDistanceKVector(catalog, minDistance, maxDistance, numBins); + unsigned char *result = new unsigned char[*length]; + SerializePairDistanceKVector(catalog, minDistance, maxDistance, numBins, result); + return result; +} + +TEST_CASE("Kvector full database stuff", "[kvector]") { + long length; + Catalog &catalog = CatalogRead(); + unsigned char *dbBytes = BuildPairDistanceKVectorDatabase(catalog, &length, DegToRad(1.0), DegToRad(2.0), 100); + REQUIRE(length < 999999); + PairDistanceKVectorDatabase db(dbBytes); + + SECTION("basic consistency checks") { + long lastNumReturnedPairs = 999999; + for (float i = 1.1; i < 1.99; i += 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; + for (const int16_t *pair = pairs; pair != end; pair += 2) { + float 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); + } + long numReturnedPairs = (end - pairs)/2; + REQUIRE(0 < numReturnedPairs); + REQUIRE(numReturnedPairs < lastNumReturnedPairs); + REQUIRE(shortestDistance < (i + 0.01) * M_PI/180.0); + lastNumReturnedPairs = numReturnedPairs; + } + } + + SECTION("form a partition") { + long totalReturnedPairs = 0; + for (float i = 1.1; i < 2.01; i+= 0.1) { + const int16_t *end; + const int16_t *pairs = db.FindPairsLiberal(DegToRad(i-0.1)+0.00001, DegToRad(i)-0.00001, &end); + long numReturnedPairs = (end-pairs)/2; + totalReturnedPairs += numReturnedPairs; + } + REQUIRE(totalReturnedPairs == db.NumPairs()); + } +} + +TEST_CASE("Tighter tolerance test", "[kvector]") { + long length; + Catalog &catalog = CatalogRead(); + unsigned char *dbBytes = BuildPairDistanceKVectorDatabase(catalog, &length, DegToRad(0.5), DegToRad(5.0), 1000); + REQUIRE(length < 999999); + PairDistanceKVectorDatabase db(dbBytes); + // radius we'll request + float delta = 0.0001; + // radius we expect back: radius we request + width of a bin + float epsilon = delta + DegToRad(5.0 - 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. + bool outsideRangeReturned = false; + for (float i = DegToRad(0.6); i < DegToRad(4.9); i += DegToRad(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); + // only need to check one side, since we're only looking for one exception. + if (i - delta > distance) { + outsideRangeReturned = true; + } + CHECK(i - epsilon <=distance); + CHECK(distance<= i + epsilon); + } + } + CHECK(outsideRangeReturned); +} + +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), + }; + unsigned char *dbBytes = BuildPairDistanceKVectorDatabase(tripleCatalog, NULL, DegToRad(0.5), DegToRad(20.0), 1000); + PairDistanceKVectorDatabase db(dbBytes); + REQUIRE(db.NumPairs() == 3); + + float distances[] = {0.038823101, 0.157079488, 0.177976221}; + for (float distance : distances) { + const int16_t *end; + const int16_t *pairs = db.FindPairsLiberal(distance - 0.001, distance + 0.001, &end); + REQUIRE(end - pairs == 2); + CHECK(AngleUnit(tripleCatalog[pairs[0]].spatial, tripleCatalog[pairs[1]].spatial) == Approx(distance)); + } +} diff --git a/test/scripts/pyramid-incorrect.sh b/test/scripts/pyramid-incorrect.sh index 889c3ce9..928049c3 100755 --- a/test/scripts/pyramid-incorrect.sh +++ b/test/scripts/pyramid-incorrect.sh @@ -32,7 +32,7 @@ for _ in $(seq "${1:-10}"); do de=$(rand_int -89 89) roll=$(rand_int 0 359) tolerance=$(awk "BEGIN{print $RANDOM%10/100+.002}") # floating point numbers are hard, alright? - set -x + # set -x lost_output=$( ./lost pipeline \ --generate 1 \ @@ -52,7 +52,7 @@ for _ in $(seq "${1:-10}"); do --max-mismatch-prob 0.001 \ --compare-star-ids ) - set +x + # set +x num_incorrect_stars=$(grep -oP "(?<=starid_num_incorrect )\\d+" <<<"$lost_output") num_correct_stars=$(grep -oP "(?<=starid_num_correct )\\d+" <<<"$lost_output") if ((num_correct_stars > 0 && num_incorrect_stars == 0)); then diff --git a/test/scripts/random-crap.sh b/test/scripts/random-crap.sh index 57e0d014..71596c1b 100755 --- a/test/scripts/random-crap.sh +++ b/test/scripts/random-crap.sh @@ -31,6 +31,16 @@ echo 'Issue #36: Cog and Attitude without Star-ID' echo 'Run the generator without centroids a whole bunch and make sure no assertions go off' ./lost pipeline --generate 200 --generate-perturb-centroids 5 --generate-centroids-only +echo 'Speed 95-th percentile should be different than max for 20 but not 19 trials' +nineteen_out=$(./lost pipeline --generate 19 --generate-centroids-only --attitude-algo quest --print-speed -) +nineteen_max_ns=$(echo "$nineteen_out" | grep total_max_ns | cut -d' ' -f2) +nineteen_95_ns=$(echo "$nineteen_out" | grep 'total_95%_ns' | cut -d' ' -f2) +(( nineteen_max_ns == nineteen_95_ns )) || exit 1 # somehow single equals sign doesn't work? I am confusion. Probably because they actually have that do assignment, for for-loop purposes? +twenty_out=$(./lost pipeline --generate 20 --generate-centroids-only --attitude-algo quest --print-speed -) +twenty_max_ns=$(echo "$twenty_out" | grep total_max_ns | cut -d' ' -f2) +twenty_95_ns=$(echo "$twenty_out" | grep 'total_95%_ns' | cut -d' ' -f2) +(( twenty_max_ns > twenty_95_ns )) || exit 1 + set +x echo '