From ae4c90daf46e58ec834ad682a97e10fa36118b1e Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Thu, 2 Mar 2023 21:54:28 -0800 Subject: [PATCH 01/34] Start implementing Least Squares Gaussian Fit 1D --- .gitignore | 3 + src/centroiders.cpp | 661 ++++++++++++++++++++++++++------------------ src/centroiders.hpp | 7 + 3 files changed, 400 insertions(+), 271 deletions(-) diff --git a/.gitignore b/.gitignore index d8ebea35..d08a4241 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,9 @@ *.png *.dat *.txt +*.bin + +BSC5 /documentation/man-*.h diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 0ddb6103..55971402 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -1,325 +1,444 @@ #include "centroiders.hpp" +#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 +#include namespace lost { -// DUMMY +int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight); -std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, int imageHeight) const { - std::vector result; +// 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; +}; + +// TODO: function prototypes for the helper functions, please! - 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)); +class LSGF1DFunctor { + public: + LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) + : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} + + int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &residuals) const { + float a = params(0); + int xb = params(1); + float sigma = params(2); + for (int i = 0; i < X.size(); i++) { + residuals(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); + } + return 0; + } + + private: + const Eigen::VectorXd X; + // Eigen::VectorXd Y; // Don't need, calculate Y_i = V_{m, x_i} + const int nb = 2; + const int x0, y0; // coordinates of centroid window + const int w, h; // w = number of pixels in width, h = number of pixels in height + const unsigned char *image; + + // Get value of pixel at universal coordinates (x, y) + // TODO: is this correct representation? + // TODO: will need to deal with case where on boundary of acceptable nb + unsigned char Get(int x, int y) const { + // int ind = i * w + j; + int ind = x + y * w; + return image[ind]; + } + + // Assuming 0, 0 is the center + int XMarginal(int x) const { + int res = 0; + for (int j = -nb; j <= nb; j++) { + res += Get(x0 + x, y0 + j); } + return res; + } + + float Model(int xi, float a, int xb, float sigma) const { + return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); + } +}; + +typedef std::vector Point; + +std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + std::vector result; + + const int cb = 2; + std::set checkedPoints; + + int cutoff = BasicThreshold(image, imageWidth, imageHeight); + for (int i = 0; i < imageHeight * imageWidth; i++) { + int x = i / imageWidth; + int y = i % imageHeight; + Point p{x, y}; + if(x-cb < 0 || x+cb >= imageWidth || y-cb < 0 || y+cb >= imageHeight) continue; + if (image[i] >= cutoff && checkedPoints.count(p) == 0) { + + checkedPoints.insert(p); + + std::vector xs; + std::vector ys; + Eigen::VectorXd X(2*cb + 1); + for(int i = -cb; i <= cb; i++){ + xs.push_back(x+cb); + X << (x+cb); + ys.push_back(y+cb); + } + + // LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, + // const unsigned char *image) + // : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - return result; + Eigen::VectorXd beta(3); // initial guess + // TODO: fill beta with described values + beta << 1, 1, 1; + + LSGF1DFunctor functor(X, x, y, imageWidth, imageHeight, image); + Eigen::NumericalDiff numDiff(functor); + Eigen::LevenbergMarquardt, float> lm(numDiff); + + lm.parameters.maxfev = 1000; + lm.parameters.xtol = 1.0e-10; + + lm.minimize(beta); + + float a = beta(0); + int xb = beta(1); + float sigma = beta(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)); + } + + return result; } // a poorly designed thresholding algorithm int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) { - //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; + // 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; } // 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 sumB = 0; - float sum1 = 0; - float wB = 0; - float maximum = 0; - int level = 0; - // make the histogram (array length 256) - int histogram[256]; - - memset(histogram, 0, sizeof(int)*256); - - for (long i = 0; i < total; i++) { - histogram[image[i]]++; - } - for (int i = 0; i < 256; i ++) { - sum1 += i * histogram[i]; + // code here, duh + long total = imageWidth * imageHeight; + // float top = 255; + float sumB = 0; + float sum1 = 0; + float wB = 0; + float maximum = 0; + int level = 0; + // make the histogram (array length 256) + int histogram[256]; + + memset(histogram, 0, sizeof(int) * 256); + + for (long i = 0; i < total; i++) { + histogram[image[i]]++; + } + for (int i = 0; i < 256; i++) { + sum1 += i * histogram[i]; + } + for (int i = 0; i < 256; i++) { + float 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); + // std::cout << val << "\n"; + if (val >= maximum) { + level = i; + maximum = val; + } } - for (int i = 0; i < 256; i ++) { - float 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); - //std::cout << val << "\n"; - if (val >= maximum) { - level = i; - maximum = val; - } - } - wB = wB + histogram[i]; - sumB = sumB + i * histogram[i]; - } - return level; + wB = wB + histogram[i]; + sumB = sumB + i * histogram[i]; + } + 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; - long totalPixels = imageHeight * imageWidth; - for (long i = 0; i < totalPixels; i++) { - totalMag += image[i]; - } - float mean = totalMag / totalPixels; - for (long i = 0; i < totalPixels; i++) { - std += std::pow(image[i] - mean, 2); - } - std = std::sqrt(std / totalPixels); - return mean + (std * 5); + unsigned long totalMag = 0; + float std = 0; + long totalPixels = imageHeight * imageWidth; + for (long i = 0; i < totalPixels; i++) { + totalMag += image[i]; + } + float mean = totalMag / totalPixels; + for (long i = 0; i < totalPixels; i++) { + std += std::pow(image[i] - mean, 2); + } + std = std::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; - 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); - return mean + (std * 5); + unsigned long totalMag = 0; + float std = 0; + float 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); + 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) { - p->isValid = false; - } - p->checkedIndices.insert(i); - if (i % imageWidth > p->xMax) { - p->xMax = i % imageWidth; - } else if (i % imageWidth < p->xMin) { - p->xMin = i % imageWidth; - } - if (i / imageWidth > p->yMax) { - p->yMax = i / imageWidth; - } else if (i / imageWidth < p->yMin) { - p->yMin = i / imageWidth; - } - p->magSum += image[i]; - p->xCoordMagSum += ((i % imageWidth)) * image[i]; - p->yCoordMagSum += ((i / imageWidth)) * image[i]; - if (i % imageWidth != imageWidth - 1) { - CogHelper(p, i + 1, image, imageWidth, imageHeight); - } - if (i % imageWidth != 0) { - CogHelper(p, i - 1, image, imageWidth, imageHeight); - } - CogHelper(p, i + imageWidth, image, imageWidth, imageHeight); - CogHelper(p, i - imageWidth, image, imageWidth, 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) { + p->isValid = false; + } + p->checkedIndices.insert(i); + if (i % imageWidth > p->xMax) { + p->xMax = i % imageWidth; + } else if (i % imageWidth < p->xMin) { + p->xMin = i % imageWidth; + } + if (i / imageWidth > p->yMax) { + p->yMax = i / imageWidth; + } else if (i / imageWidth < p->yMin) { + p->yMin = i / imageWidth; + } + p->magSum += image[i]; + p->xCoordMagSum += ((i % imageWidth)) * image[i]; + p->yCoordMagSum += ((i / imageWidth)) * image[i]; + if (i % imageWidth != imageWidth - 1) { + CogHelper(p, i + 1, image, imageWidth, imageHeight); + } + if (i % imageWidth != 0) { + CogHelper(p, i - 1, image, imageWidth, imageHeight); } + CogHelper(p, i + imageWidth, image, imageWidth, imageHeight); + CogHelper(p, i - imageWidth, image, imageWidth, imageHeight); + } } -std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { - CentroidParams p; - - std::vector result; - - 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 - 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.xMax = i % imageWidth; - p.xMin = i % imageWidth; - p.yMax = i / imageWidth; - p.yMin = i / imageWidth; - p.isValid = true; - - int sizeBefore = p.checkedIndices.size(); - - CogHelper(&p, i, image, imageWidth, imageHeight); - 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 - 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)); - } - } +std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + CentroidParams p; + + std::vector result; + + 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 + 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.xMax = i % imageWidth; + p.xMin = i % imageWidth; + p.yMax = i / imageWidth; + p.yMin = i / imageWidth; + p.isValid = true; + + int sizeBefore = p.checkedIndices.size(); + + CogHelper(&p, i, image, imageWidth, imageHeight); + 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 + 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)); + } } - return result; + } + 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 { - int xMin; - int xMax; - int yMin; - int yMax; - int cutoff; - int maxIntensity; - int guess; - bool isValid; - std::unordered_set checkedIndices; + int xMin; + int xMax; + int yMin; + int yMax; + int cutoff; + int maxIntensity; + int guess; + bool isValid; + 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) { - p->isValid = false; - } - p->checkedIndices.insert(i); - starIndices->push_back(i); - if (image[i] > p->maxIntensity) { - p->maxIntensity = image[i]; - p->guess = i; - } - if (i % imageWidth > p->xMax) { - p->xMax = i % imageWidth; - } else if (i % imageWidth < p->xMin) { - p->xMin = i % imageWidth; - } - if (i / imageWidth > p->yMax) { - p->yMax = i / imageWidth; - } else if (i / imageWidth < p->yMin) { - p->yMin = i / imageWidth; - } - if (i % imageWidth != imageWidth - 1) { - IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, starIndices); - } - if (i % imageWidth != 0) { - IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, starIndices); - } - IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, starIndices); - IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, starIndices); +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); + starIndices->push_back(i); + if (image[i] > p->maxIntensity) { + p->maxIntensity = image[i]; + p->guess = i; + } + if (i % imageWidth > p->xMax) { + p->xMax = i % imageWidth; + } else if (i % imageWidth < p->xMin) { + p->xMin = i % imageWidth; + } + if (i / imageWidth > p->yMax) { + p->yMax = i / imageWidth; + } else if (i / imageWidth < p->yMin) { + p->yMin = i / imageWidth; } + if (i % imageWidth != imageWidth - 1) { + IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, starIndices); + } + if (i % imageWidth != 0) { + IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, starIndices); + } + IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, starIndices); + } } -Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { - IWCoGParams p; - std::vector result; - 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 - if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { - // TODO: store longs --Mark - 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 standardDeviation; - float w; //weight value - - p.xMax = i % imageWidth; - p.xMin = i % imageWidth; - p.yMax = i / imageWidth; - 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 - float 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); - //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 (change > iWCoGMinChange && stop < 100000) { - //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)]); - } - float xTemp = xWeightedCoordMagSum / weightedMagSum; - float yTemp = yWeightedCoordMagSum / weightedMagSum; - - change = abs(guessXCoord - xTemp) + abs(guessYCoord - yTemp); - - guessXCoord = xTemp; - 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())); - } +Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + IWCoGParams p; + std::vector result; + 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 + if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { + // TODO: store longs --Mark + 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 standardDeviation; + float w; // weight value + + p.xMax = i % imageWidth; + p.xMin = i % imageWidth; + p.yMax = i / imageWidth; + 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 + float 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); + // 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 (change > iWCoGMinChange && stop < 100000) { + // 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)]); + } + float xTemp = xWeightedCoordMagSum / weightedMagSum; + float yTemp = yWeightedCoordMagSum / weightedMagSum; + + change = abs(guessXCoord - xTemp) + abs(guessYCoord - yTemp); + + guessXCoord = xTemp; + 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())); + } } - return result; + } + return result; } -} +} // namespace lost diff --git a/src/centroiders.hpp b/src/centroiders.hpp index d18683f7..6a6f5164 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -21,6 +21,13 @@ class CentroidAlgorithm { virtual ~CentroidAlgorithm() { }; }; +class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ +public: + explicit LeastSquaresGaussianFit1D() {}; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; +}; + + /// A centroid algorithm for debugging that returns random centroids. class DummyCentroidAlgorithm: public CentroidAlgorithm { public: From c3cf0fa105fcd5c770a5ee204098f6c052fe4d03 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Fri, 3 Mar 2023 10:46:06 -0800 Subject: [PATCH 02/34] Now using LM correctly...I think --- src/centroiders.cpp | 101 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 79 insertions(+), 22 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 55971402..1a9399d9 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -7,11 +7,11 @@ #include #include +#include #include #include // TODO: remove later, this is just lazy to implement hash for unordered_set #include #include -#include #include namespace lost { @@ -32,19 +32,35 @@ struct CentroidParams { std::unordered_set checkedIndices; }; -// TODO: function prototypes for the helper functions, please! +// 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, m_values; + + 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; } +}; -class LSGF1DFunctor { +struct LSGF1DFunctor : Functor { public: LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) - : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} + : Functor(3, 3), X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &residuals) const { - float a = params(0); + int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &fvec) const { + double a = params(0); int xb = params(1); - float sigma = params(2); + double sigma = params(2); for (int i = 0; i < X.size(); i++) { - residuals(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); + fvec(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); } return 0; } @@ -75,13 +91,50 @@ class LSGF1DFunctor { return res; } - float Model(int xi, float a, int xb, float sigma) const { + double Model(int xi, double a, int xb, double sigma) const { return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); } }; typedef std::vector Point; +// TODO: duplicate, factor out better +unsigned char Get(int x, int y, unsigned char *image, int w, int h) { + int ind = x + y * w; + return image[ind]; +} + +void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, int *xb, + double *sigma) { + // a is set to max intensity value in the window + // (xb, yb) = coordinates of pixel with max intensity value + const int nb = 2; + int max = -1; + for (int i = -nb; i <= nb; i++) { + for (int j = -nb; j <= nb; j++) { + int pixelValue = Get(x + i, y + j, image, w, h); + if (pixelValue > max) { + *xb = x + i; + max = pixelValue; + } + } + } + *a = max; + // Sigma is kinda 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(x + i, y + j, image, w, h) > *a / 2){ + halfCount++; + } + } + } + double fwhm = std::sqrt(halfCount); + *sigma = fwhm / (2 * std::sqrt(2 * std::log(2))); +} + std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageWidth, int imageHeight) const { std::vector result; @@ -94,41 +147,45 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int x = i / imageWidth; int y = i % imageHeight; Point p{x, y}; - if(x-cb < 0 || x+cb >= imageWidth || y-cb < 0 || y+cb >= imageHeight) continue; + if (x - cb < 0 || x + cb >= imageWidth || y - cb < 0 || y + cb >= imageHeight) continue; + // TODO: check that entire window is not in checkedPoints? - no, just check the point itself if (image[i] >= cutoff && checkedPoints.count(p) == 0) { - checkedPoints.insert(p); std::vector xs; std::vector ys; - Eigen::VectorXd X(2*cb + 1); - for(int i = -cb; i <= cb; i++){ - xs.push_back(x+cb); - X << (x+cb); - ys.push_back(y+cb); + Eigen::VectorXd X(2 * cb + 1); + for (int i = -cb; i <= cb; i++) { + xs.push_back(x + cb); + X << (x + cb); + ys.push_back(y + cb); } // LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, // const unsigned char *image) // : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - Eigen::VectorXd beta(3); // initial guess + Eigen::VectorXd beta(3); // initial guess + double a, sigma; + int xb; + InitialGuess(x, y, image, imageWidth, imageHeight, &a, &xb, &sigma); // TODO: fill beta with described values - beta << 1, 1, 1; + beta << a, xb, sigma; LSGF1DFunctor functor(X, x, y, imageWidth, imageHeight, image); Eigen::NumericalDiff numDiff(functor); - Eigen::LevenbergMarquardt, float> lm(numDiff); + Eigen::LevenbergMarquardt, double> lm(numDiff); lm.parameters.maxfev = 1000; lm.parameters.xtol = 1.0e-10; lm.minimize(beta); - float a = beta(0); - int xb = beta(1); - float sigma = beta(2); + a = beta(0); + xb = beta(1); + sigma = beta(2); + std::cout << a << ", " << xb << ", " << sigma << std::endl; } } From 47916231d95728997dff70c35c8df9e6d923fd3b Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Fri, 3 Mar 2023 12:18:52 -0800 Subject: [PATCH 03/34] Enable calculation for both x and y marginals --- src/centroiders.cpp | 101 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 86 insertions(+), 15 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 1a9399d9..5590b7f8 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -50,9 +50,9 @@ struct Functor { int values() const { return m_values; } }; -struct LSGF1DFunctor : Functor { +struct LSGF1DXFunctor : Functor { public: - LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) + LSGF1DXFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) : Functor(3, 3), X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &fvec) const { @@ -71,6 +71,7 @@ struct LSGF1DFunctor : Functor { const int nb = 2; const int x0, y0; // coordinates of centroid window const int w, h; // w = number of pixels in width, h = number of pixels in height + // TODO: h is unused here const unsigned char *image; // Get value of pixel at universal coordinates (x, y) @@ -96,6 +97,53 @@ struct LSGF1DFunctor : Functor { } }; +struct LSGF1DYFunctor : Functor { + public: + LSGF1DYFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) + : Functor(3, 3), X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} + + int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &fvec) const { + double a = params(0); + int yb = params(1); + double sigma = params(2); + for (int i = 0; i < X.size(); i++) { + fvec(i) = YMarginal(X(i)) - Model(X(i), a, yb, sigma); + } + return 0; + } + + private: + const Eigen::VectorXd X; + // Eigen::VectorXd Y; // Don't need, calculate Y_i = V_{m, x_i} + const int nb = 2; + const int x0, y0; // coordinates of centroid window + const int w, h; // w = number of pixels in width, h = number of pixels in height + const unsigned char *image; + + // Get value of pixel at universal coordinates (x, y) + // TODO: is this correct representation? + // TODO: will need to deal with case where on boundary of acceptable nb + unsigned char Get(int x, int y) const { + // int ind = i * w + j; + int ind = x + y * w; + return image[ind]; + } + + // Assuming 0, 0 is the center + int YMarginal(int y) const { + int res = 0; + for (int j = -nb; j <= nb; j++) { + res += Get(x0 + j, y0 + y); + } + return res; + } + + double Model(int yi, double a, int yb, double sigma) const { + // Technically sigma_y + return a * exp(-1 * (yi - yb) * (yi - yb) / (2 * sigma * sigma)); + } +}; + typedef std::vector Point; // TODO: duplicate, factor out better @@ -104,7 +152,7 @@ unsigned char Get(int x, int y, unsigned char *image, int w, int h) { return image[ind]; } -void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, int *xb, +void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, int *xb, int *yb, double *sigma) { // a is set to max intensity value in the window // (xb, yb) = coordinates of pixel with max intensity value @@ -115,6 +163,7 @@ void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, i int pixelValue = Get(x + i, y + j, image, w, h); if (pixelValue > max) { *xb = x + i; + *yb = y + j; max = pixelValue; } } @@ -155,37 +204,59 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW std::vector xs; std::vector ys; Eigen::VectorXd X(2 * cb + 1); + Eigen::VectorXd Y(2 * cb + 1); for (int i = -cb; i <= cb; i++) { xs.push_back(x + cb); X << (x + cb); ys.push_back(y + cb); + Y << (y + cb); } // LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, // const unsigned char *image) // : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - Eigen::VectorXd beta(3); // initial guess + Eigen::VectorXd betaX(3); // initial guess + Eigen::VectorXd betaY(3); // initial guess + double a, sigma; - int xb; - InitialGuess(x, y, image, imageWidth, imageHeight, &a, &xb, &sigma); + int xb, yb; + InitialGuess(x, y, image, imageWidth, imageHeight, &a, &xb, &yb, &sigma); // TODO: fill beta with described values - beta << a, xb, sigma; + betaX << a, xb, sigma; + betaY << a, yb, sigma; - LSGF1DFunctor functor(X, x, y, imageWidth, imageHeight, image); - Eigen::NumericalDiff numDiff(functor); - Eigen::LevenbergMarquardt, double> lm(numDiff); + LSGF1DXFunctor functorX(X, x, y, imageWidth, imageHeight, image); + Eigen::NumericalDiff numDiffX(functorX); + Eigen::LevenbergMarquardt, double> lm(numDiffX); lm.parameters.maxfev = 1000; lm.parameters.xtol = 1.0e-10; - lm.minimize(beta); + lm.minimize(betaX); + + xb = betaX(1); + + LSGF1DYFunctor functorY(X, x, y, imageWidth, imageHeight, image); + Eigen::NumericalDiff numDiffY(functorY); + Eigen::LevenbergMarquardt, double> lmY(numDiffY); + + lmY.parameters.maxfev = 1000; + lmY.parameters.xtol = 1.0e-10; + + lmY.minimize(betaY); + + yb = betaY(1); + + a = betaX(0); + sigma = betaX(2); - a = beta(0); - xb = beta(1); - sigma = beta(2); + std::cout << a << ", " << xb << ", " << yb << ", " << sigma << std::endl; - std::cout << a << ", " << xb << ", " << sigma << std::endl; + // Add new centroid to result + // 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(xb, yb, sigma, sigma, a)); } } From 92f94780b901252c4d48988c2d2dcacaf2693777 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Fri, 3 Mar 2023 22:37:59 -0800 Subject: [PATCH 04/34] merge --- Makefile | 6 +- test/kvector.off | 111 ++++++++++++++++++++++++++++++ test/scripts/pyramid-incorrect.sh | 4 +- 3 files changed, 116 insertions(+), 5 deletions(-) create mode 100644 test/kvector.off diff --git a/Makefile b/Makefile index c28483a6..9a0d7875 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 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 From 7dc98a7580dedf89ad92786e1c88fc961033e432 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 4 Mar 2023 11:35:22 -0800 Subject: [PATCH 05/34] Figured out how to use LM, finally --- Makefile | 5 +- src/centroiders.cpp | 334 +++++++++++++++++++++++--------------------- src/io.cpp | 12 +- 3 files changed, 185 insertions(+), 166 deletions(-) diff --git a/Makefile b/Makefile index 9a0d7875..8e364f87 100644 --- a/Makefile +++ b/Makefile @@ -85,8 +85,9 @@ lint: test: $(BIN) $(BSC) $(TEST_BIN) $(TEST_BIN) - # bash ./test/scripts/pyramid-incorrect.sh - bash ./test/scripts/readme-examples-test.sh + bash ./test/scripts/pyramid-incorrect.sh + # bash ./test/scripts/readme-examples-test.sh + # bash ./test/scripts/tetra.sh $(TEST_BIN): $(TEST_OBJS) $(CXX) $(LDFLAGS) -o $(TEST_BIN) $(TEST_OBJS) $(LIBS) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 5590b7f8..24f8a020 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -32,150 +33,74 @@ struct CentroidParams { std::unordered_set checkedIndices; }; -// 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, m_values; - - 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; } -}; - -struct LSGF1DXFunctor : Functor { - public: - LSGF1DXFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) - : Functor(3, 3), X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - - int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &fvec) const { - double a = params(0); - int xb = params(1); - double sigma = params(2); - for (int i = 0; i < X.size(); i++) { - fvec(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); - } - return 0; - } - - private: - const Eigen::VectorXd X; - // Eigen::VectorXd Y; // Don't need, calculate Y_i = V_{m, x_i} - const int nb = 2; - const int x0, y0; // coordinates of centroid window - const int w, h; // w = number of pixels in width, h = number of pixels in height - // TODO: h is unused here - const unsigned char *image; - - // Get value of pixel at universal coordinates (x, y) - // TODO: is this correct representation? - // TODO: will need to deal with case where on boundary of acceptable nb - unsigned char Get(int x, int y) const { - // int ind = i * w + j; - int ind = x + y * w; - return image[ind]; - } - - // Assuming 0, 0 is the center - int XMarginal(int x) const { - int res = 0; - for (int j = -nb; j <= nb; j++) { - res += Get(x0 + x, y0 + j); - } - return res; - } - - double Model(int xi, double a, int xb, double sigma) const { - return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); - } -}; - -struct LSGF1DYFunctor : Functor { - public: - LSGF1DYFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, const unsigned char *image) - : Functor(3, 3), X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - - int operator()(const Eigen::VectorXd ¶ms, Eigen::VectorXd &fvec) const { - double a = params(0); - int yb = params(1); - double sigma = params(2); - for (int i = 0; i < X.size(); i++) { - fvec(i) = YMarginal(X(i)) - Model(X(i), a, yb, sigma); - } - return 0; - } - - private: - const Eigen::VectorXd X; - // Eigen::VectorXd Y; // Don't need, calculate Y_i = V_{m, x_i} - const int nb = 2; - const int x0, y0; // coordinates of centroid window - const int w, h; // w = number of pixels in width, h = number of pixels in height - const unsigned char *image; - - // Get value of pixel at universal coordinates (x, y) - // TODO: is this correct representation? - // TODO: will need to deal with case where on boundary of acceptable nb - unsigned char Get(int x, int y) const { - // int ind = i * w + j; - int ind = x + y * w; - return image[ind]; - } +/* +Get pixel value at (x, y) where x and y are 0-based +Note: top left is (0, 0) +*/ +int Get(int x, int y, const unsigned char *image, int w) { + int ind = y * w + x; + return image[ind]; +} - // Assuming 0, 0 is the center - int YMarginal(int y) const { - int res = 0; - for (int j = -nb; j <= nb; j++) { - res += Get(x0 + j, y0 + y); - } - return res; +/* +Get XMarginal - 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) { + int sum = 0; + for (int j = -nb; j <= nb; j++) { + sum += Get(x, y0 + j, image, w); } + return sum; +} - double Model(int yi, double a, int yb, double sigma) const { - // Technically sigma_y - return a * exp(-1 * (yi - yb) * (yi - yb) / (2 * sigma * sigma)); +/* +Get YMarginal - 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) { + int sum = 0; + for (int i = -nb; i <= nb; i++) { + sum += Get(x0 + i, y, image, w); } -}; - -typedef std::vector Point; + return sum; +} -// TODO: duplicate, factor out better -unsigned char Get(int x, int y, unsigned char *image, int w, int h) { - int ind = x + y * w; - return image[ind]; +/* +Can be used as f(xi, beta) or f(yi, beta) +*/ +int FitModel(int xi, int a, int xb, float sigma) { + return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); } -void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, int *xb, int *yb, - double *sigma) { +/* +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, int *a, int *xb, + int *yb, double *sigma) { // a is set to max intensity value in the window // (xb, yb) = coordinates of pixel with max intensity value - const int nb = 2; int max = -1; for (int i = -nb; i <= nb; i++) { for (int j = -nb; j <= nb; j++) { - int pixelValue = Get(x + i, y + j, image, w, h); + int pixelValue = Get(x0 + i, y0 + j, image, w); if (pixelValue > max) { - *xb = x + i; - *yb = y + j; + *xb = x0 + i; + *yb = y0 + j; max = pixelValue; } } } *a = max; - // Sigma is kinda complicated + // 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(x + i, y + j, image, w, h) > *a / 2){ + if (Get(x0 + i, y0 + j, image, w) > *a / 2) { halfCount++; } } @@ -184,11 +109,108 @@ void InitialGuess(int x, int y, unsigned char *image, int w, int h, double *a, i *sigma = 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, m_values; + + 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; } +}; + +// struct my_functor : Functor { +// // First param = number of parameters in your model +// // Second param = number of data points you want to test over +// // for us, = 2 * nb + 1 +// my_functor(void) : Functor(2, 2) {} +// int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { +// // Implement y = 10*(x0+3)^2 + (x1-5)^2 +// fvec(0) = 10.0 * pow(x(0) + 3.0, 2) + pow(x(1) - 5.0, 2); +// fvec(1) = 0; + +// return 0; +// } +// }; + +// void foo(int (*f)(int, int, int, unsigned char *, int), ) { +// std::cout << f +// } + +struct LSGFFunctor : Functor { + // First param = number of parameters in your model + // Second param = number of data points you want to test over + // for us, = 2 * nb + 1 + // LSGFFunctor(void) : Functor(2, 2) {} + LSGFFunctor(const int nb, const int marg, Eigen::VectorXd 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) {} + + /* + x = parameters (a, xb, sigma) + fvec = residual (one for every data point, of size = 2*nb+1) + */ + int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { + // // Implement y = 10*(x0+3)^2 + (x1-5)^2 + // fvec(0) = 10.0 * pow(x(0) + 3.0, 2) + pow(x(1) - 5.0, 2); + // fvec(1) = 0; + int a = x(0); + int xb = x(1); + double 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) = marginal - FitModel(X(i), a, xb, sigma); + } + // double a = params(0); + // int xb = params(1); + // double sigma = params(2); + // for (int i = 0; i < X.size(); i++) { + // fvec(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); + // } + + return 0; + } + + const int nb; + const int marg; + Eigen::VectorXd X; + const unsigned char *image; + const int w; + const int x0, y0; + + // Flag for XMarginal or YMarginal functor + // Field for data X + // Field for image, imageWidth + // Field for center coordinates (x0, y0) of this window +}; + +typedef std::vector Point; + std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageWidth, int imageHeight) const { std::vector result; - const int cb = 2; + // const int cb = 2; + const int nb = 2; std::set checkedPoints; int cutoff = BasicThreshold(image, imageWidth, imageHeight); @@ -196,66 +218,56 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int x = i / imageWidth; int y = i % imageHeight; Point p{x, y}; - if (x - cb < 0 || x + cb >= imageWidth || y - cb < 0 || y + cb >= imageHeight) continue; - // TODO: check that entire window is not in checkedPoints? - no, just check the point itself + if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; + // // TODO: check that entire window is not in checkedPoints? - no, just check the point + // itself if (image[i] >= cutoff && checkedPoints.count(p) == 0) { checkedPoints.insert(p); - std::vector xs; - std::vector ys; - Eigen::VectorXd X(2 * cb + 1); - Eigen::VectorXd Y(2 * cb + 1); - for (int i = -cb; i <= cb; i++) { - xs.push_back(x + cb); - X << (x + cb); - ys.push_back(y + cb); - Y << (y + cb); + Eigen::VectorXd X(2 * nb + 1); + Eigen::VectorXd Y(2 * nb + 1); + int vInd = 0; + for (int i = -nb; i <= nb; i++) { + X(vInd) = x + i; + Y(vInd) = y + i; + vInd++; } - // LSGF1DFunctor(const Eigen::VectorXd X, int x0, int y0, int w, int h, - // const unsigned char *image) - // : X(X), x0(x0), y0(y0), w(w), h(h), image(image) {} - - Eigen::VectorXd betaX(3); // initial guess - Eigen::VectorXd betaY(3); // initial guess - - double a, sigma; + int a; int xb, yb; - InitialGuess(x, y, image, imageWidth, imageHeight, &a, &xb, &yb, &sigma); - // TODO: fill beta with described values - betaX << a, xb, sigma; - betaY << a, yb, sigma; + double sigma; + InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); - LSGF1DXFunctor functorX(X, x, y, imageWidth, imageHeight, image); - Eigen::NumericalDiff numDiffX(functorX); - Eigen::LevenbergMarquardt, double> lm(numDiffX); + Eigen::VectorXd betaX(3); + betaX << a, xb, sigma; - lm.parameters.maxfev = 1000; - lm.parameters.xtol = 1.0e-10; + Eigen::VectorXd betaY(3); + betaY << a, yb, sigma; - lm.minimize(betaX); + LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffX(functorX); + Eigen::LevenbergMarquardt, double> lmX(numDiffX); - xb = betaX(1); + lmX.parameters.maxfev = 2000; + lmX.parameters.xtol = 1.0e-10; - LSGF1DYFunctor functorY(X, x, y, imageWidth, imageHeight, image); - Eigen::NumericalDiff numDiffY(functorY); - Eigen::LevenbergMarquardt, double> lmY(numDiffY); + LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffY(functorY); + Eigen::LevenbergMarquardt, double> lmY(numDiffY); - lmY.parameters.maxfev = 1000; + lmY.parameters.maxfev = 2000; lmY.parameters.xtol = 1.0e-10; + lmX.minimize(betaX); lmY.minimize(betaY); + xb = betaX(1); yb = betaY(1); - a = betaX(0); sigma = betaX(2); - std::cout << a << ", " << xb << ", " << yb << ", " << sigma << std::endl; + // std::cout << a << ", " << xb << ", " << yb << ", " << sigma << std::endl; - // Add new centroid to result - // 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(xb, yb, sigma, sigma, a)); } } diff --git a/src/io.cpp b/src/io.cpp index c7634fa1..7c317fab 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -713,11 +713,17 @@ 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()); } else if (values.centroidAlgo != "") { std::cout << "Illegal centroid algorithm." << std::endl; exit(1); From c812f3fd4862f24d57bc01d835ced3e644a9f604 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 4 Mar 2023 11:58:07 -0800 Subject: [PATCH 06/34] Cleanup some bugs --- src/centroiders.cpp | 52 ++++++++++++--------------------------------- src/centroiders.hpp | 5 +++++ 2 files changed, 18 insertions(+), 39 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 24f8a020..e086749b 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -127,29 +127,11 @@ struct Functor { int values() const { return m_values; } }; -// struct my_functor : Functor { -// // First param = number of parameters in your model -// // Second param = number of data points you want to test over -// // for us, = 2 * nb + 1 -// my_functor(void) : Functor(2, 2) {} -// int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { -// // Implement y = 10*(x0+3)^2 + (x1-5)^2 -// fvec(0) = 10.0 * pow(x(0) + 3.0, 2) + pow(x(1) - 5.0, 2); -// fvec(1) = 0; - -// return 0; -// } -// }; - -// void foo(int (*f)(int, int, int, unsigned char *, int), ) { -// std::cout << f -// } struct LSGFFunctor : Functor { // First param = number of parameters in your model // Second param = number of data points you want to test over // for us, = 2 * nb + 1 - // LSGFFunctor(void) : Functor(2, 2) {} LSGFFunctor(const int nb, const int marg, Eigen::VectorXd X, const unsigned char *image, const int w, const int x0, const int y0) : Functor(3, 2 * nb + 1), @@ -166,9 +148,7 @@ struct LSGFFunctor : Functor { fvec = residual (one for every data point, of size = 2*nb+1) */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - // // Implement y = 10*(x0+3)^2 + (x1-5)^2 - // fvec(0) = 10.0 * pow(x(0) + 3.0, 2) + pow(x(1) - 5.0, 2); - // fvec(1) = 0; + int a = x(0); int xb = x(1); double sigma = x(2); @@ -180,27 +160,21 @@ struct LSGFFunctor : Functor { marginal = YMarginal(x0, X(i), nb, image, w); fvec(i) = marginal - FitModel(X(i), a, xb, sigma); } - // double a = params(0); - // int xb = params(1); - // double sigma = params(2); - // for (int i = 0; i < X.size(); i++) { - // fvec(i) = XMarginal(X(i)) - Model(X(i), a, xb, sigma); - // } 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::VectorXd X; const unsigned char *image; - const int w; + const int w; // image width in pixels + // Center coordinates (x0, y0) of this window const int x0, y0; - // Flag for XMarginal or YMarginal functor - // Field for data X - // Field for image, imageWidth - // Field for center coordinates (x0, y0) of this window }; typedef std::vector Point; @@ -214,13 +188,13 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW std::set checkedPoints; int cutoff = BasicThreshold(image, imageWidth, imageHeight); - for (int i = 0; i < imageHeight * imageWidth; i++) { - int x = i / imageWidth; - int y = i % imageHeight; + // TODO: increased step size to 10. Is this fine? + // If we just increment by 1, it's possible for a big star to contribute many centroids + for (int i = 0; i < imageHeight * imageWidth; i+=10) { + int x = i % imageWidth; + int y = i / imageWidth; Point p{x, y}; if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; - // // TODO: check that entire window is not in checkedPoints? - no, just check the point - // itself if (image[i] >= cutoff && checkedPoints.count(p) == 0) { checkedPoints.insert(p); @@ -266,12 +240,12 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW a = betaX(0); sigma = betaX(2); - // std::cout << a << ", " << xb << ", " << yb << ", " << sigma << std::endl; - result.push_back(Star(xb, yb, sigma, sigma, a)); } } + std::cout << "Number of centroids: " << result.size() << std::endl; + return result; } diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 6a6f5164..14cb2640 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -21,6 +21,11 @@ class CentroidAlgorithm { virtual ~CentroidAlgorithm() { }; }; +/** + * @brief Least Squares Gaussian Fit 1D centroiding algorithm + * Uses Levenberg-Marquardt to solve nonlinear least-squares + * + */ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ public: explicit LeastSquaresGaussianFit1D() {}; From 10c1b23fa98147da83302e7b1aa080ea1caf7770 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 4 Mar 2023 22:22:28 -0800 Subject: [PATCH 07/34] Test window sizes --- src/centroiders.cpp | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index e086749b..f1f74677 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -183,20 +183,21 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int imageHeight) const { std::vector result; - // const int cb = 2; - const int nb = 2; + const int nb = 4; + const int step = 5; + std::set checkedPoints; int cutoff = BasicThreshold(image, imageWidth, imageHeight); // TODO: increased step size to 10. Is this fine? // If we just increment by 1, it's possible for a big star to contribute many centroids - for (int i = 0; i < imageHeight * imageWidth; i+=10) { + for (int i = 0; i < imageHeight * imageWidth; i+=step) { int x = i % imageWidth; int y = i / imageWidth; Point p{x, y}; if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; if (image[i] >= cutoff && checkedPoints.count(p) == 0) { - checkedPoints.insert(p); + // checkedPoints.insert(p); Eigen::VectorXd X(2 * nb + 1); Eigen::VectorXd Y(2 * nb + 1); @@ -240,6 +241,15 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW a = betaX(0); sigma = betaX(2); + // TODO: place all points in the window to checkedPoints + for(int xr = xb - nb; xr <= xb + nb; xr++){ + for(int yr = yb - nb; yr <= yb + nb; yr++){ + if(xr < 0 || xr >= imageWidth || yr < 0 || yr >= imageHeight) continue; + Point cp{xr, yr}; + checkedPoints.insert(cp); + } + } + result.push_back(Star(xb, yb, sigma, sigma, a)); } } From ecc21477938d6203aeb0683c2ca687e65742b226 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Mon, 6 Mar 2023 00:40:44 -0800 Subject: [PATCH 08/34] Allow floating-point values for coordinates, my bad --- src/centroiders.cpp | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index f1f74677..f1127a0e 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -67,7 +67,7 @@ int YMarginal(int x0, int y, int nb, const unsigned char *image, int w) { /* Can be used as f(xi, beta) or f(yi, beta) */ -int FitModel(int xi, int a, int xb, float sigma) { +float FitModel(float xi, float a, float xb, float sigma) { return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); } @@ -78,8 +78,8 @@ 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, int *a, int *xb, - int *yb, double *sigma) { +void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, float *xb, + float *yb, double *sigma) { // a is set to max intensity value in the window // (xb, yb) = coordinates of pixel with max intensity value int max = -1; @@ -149,8 +149,8 @@ struct LSGFFunctor : Functor { */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - int a = x(0); - int xb = x(1); + float a = x(0); + float xb = x(1); double sigma = x(2); for (int i = 0; i < X.size(); i++) { int marginal; @@ -183,8 +183,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int imageHeight) const { std::vector result; - const int nb = 4; - const int step = 5; + const int nb = 2; + const int step = 1; std::set checkedPoints; @@ -208,8 +208,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW vInd++; } - int a; - int xb, yb; + float a; + float xb, yb; double sigma; InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); @@ -245,12 +245,18 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW for(int xr = xb - nb; xr <= xb + nb; xr++){ for(int yr = yb - nb; yr <= yb + nb; yr++){ if(xr < 0 || xr >= imageWidth || yr < 0 || yr >= imageHeight) continue; - Point cp{xr, yr}; - checkedPoints.insert(cp); + // std::cout << xr << ", " << yr << std::endl; + // Point cp{xr, yr}; + // checkedPoints.insert(cp); + checkedPoints.insert({xr, yr}); } } - result.push_back(Star(xb, yb, sigma, sigma, a)); + std::cout << "Original: " << x << ", " << y << std::endl; + std::cout << "final: " << xb << ", " << yb << std::endl; + // result.push_back(Star(xb, yb, sigma, sigma, a)); + result.push_back(Star(xb, yb, 0)); + // result.push_back(Star(x, y, 0)); } } From df77c31f9ca17421e99aa0e57f13edb6b979bec9 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Mon, 6 Mar 2023 01:06:24 -0800 Subject: [PATCH 09/34] Floodfill as naive preprocessing --- src/centroiders.cpp | 183 ++++++++++++++++++++++++++++---------------- 1 file changed, 115 insertions(+), 68 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index f1127a0e..372d7be4 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -17,8 +18,77 @@ namespace lost { +typedef std::vector Point; + +int Get(int x, int y, const unsigned char *image, int w); + +std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight); + int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight); +std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight) { + std::vector res; + std::set checkedPoints; + + int cutoff = BasicThreshold(image, imageWidth, imageHeight); + + 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; + + int maxMag = -1; + int x0 = x; + int y0 = y; + + queue.push_back(pCurr); + while (!queue.size() == 0) { + Point p = queue[0]; + queue.pop_front(); + + // TODO: make Point a struct probably + 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; + + // 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}); + } + } + + std::cout << "Floodfill preprocessing: found " << res.size() << " total points" << std::endl; + return res; +} + // TODO: documentation! struct CentroidParams { float yCoordMagSum; @@ -78,8 +148,8 @@ 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) { +void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, + float *xb, float *yb, double *sigma) { // a is set to max intensity value in the window // (xb, yb) = coordinates of pixel with max intensity value int max = -1; @@ -127,7 +197,6 @@ struct Functor { int values() const { return m_values; } }; - struct LSGFFunctor : Functor { // First param = number of parameters in your model // Second param = number of data points you want to test over @@ -148,7 +217,6 @@ struct LSGFFunctor : Functor { fvec = residual (one for every data point, of size = 2*nb+1) */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - float a = x(0); float xb = x(1); double sigma = x(2); @@ -160,6 +228,7 @@ struct LSGFFunctor : Functor { marginal = YMarginal(x0, X(i), nb, image, w); fvec(i) = marginal - FitModel(X(i), a, xb, sigma); } + std::cout << fvec(0) << std::endl; return 0; } @@ -171,93 +240,71 @@ struct LSGFFunctor : Functor { // Data points (set of x or y coordinates) Eigen::VectorXd X; const unsigned char *image; - const int w; // image width in pixels + const int w; // image width in pixels // Center coordinates (x0, y0) of this window const int x0, y0; - }; -typedef std::vector Point; - std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageWidth, int imageHeight) const { std::vector result; const int nb = 2; - const int step = 1; - std::set checkedPoints; + // std::set checkedPoints; - int cutoff = BasicThreshold(image, imageWidth, imageHeight); - // TODO: increased step size to 10. Is this fine? - // If we just increment by 1, it's possible for a big star to contribute many centroids - for (int i = 0; i < imageHeight * imageWidth; i+=step) { - int x = i % imageWidth; - int y = i / imageWidth; - Point p{x, y}; + 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; - if (image[i] >= cutoff && checkedPoints.count(p) == 0) { - // checkedPoints.insert(p); - - Eigen::VectorXd X(2 * nb + 1); - Eigen::VectorXd Y(2 * nb + 1); - int vInd = 0; - for (int i = -nb; i <= nb; i++) { - X(vInd) = x + i; - Y(vInd) = y + i; - vInd++; - } - float a; - float xb, yb; - double sigma; - InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); + Eigen::VectorXd X(2 * nb + 1); + Eigen::VectorXd Y(2 * nb + 1); + int vInd = 0; + for (int i = -nb; i <= nb; i++) { + X(vInd) = x + i; + Y(vInd) = y + i; + vInd++; + } - Eigen::VectorXd betaX(3); - betaX << a, xb, sigma; + float a; + float xb, yb; + double sigma; + InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); - Eigen::VectorXd betaY(3); - betaY << a, yb, sigma; + Eigen::VectorXd betaX(3); + betaX << a, xb, sigma; - LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffX(functorX); - Eigen::LevenbergMarquardt, double> lmX(numDiffX); + Eigen::VectorXd betaY(3); + betaY << a, yb, sigma; - lmX.parameters.maxfev = 2000; - lmX.parameters.xtol = 1.0e-10; + LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffX(functorX); + Eigen::LevenbergMarquardt, double> lmX(numDiffX); - LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffY(functorY); - Eigen::LevenbergMarquardt, double> lmY(numDiffY); + lmX.parameters.maxfev = 2000; + lmX.parameters.xtol = 1.0e-10; - lmY.parameters.maxfev = 2000; - lmY.parameters.xtol = 1.0e-10; + LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffY(functorY); + Eigen::LevenbergMarquardt, double> lmY(numDiffY); - lmX.minimize(betaX); - lmY.minimize(betaY); + lmY.parameters.maxfev = 2000; + lmY.parameters.xtol = 1.0e-10; - xb = betaX(1); - yb = betaY(1); - a = betaX(0); - sigma = betaX(2); + lmX.minimize(betaX); + lmY.minimize(betaY); - // TODO: place all points in the window to checkedPoints - for(int xr = xb - nb; xr <= xb + nb; xr++){ - for(int yr = yb - nb; yr <= yb + nb; yr++){ - if(xr < 0 || xr >= imageWidth || yr < 0 || yr >= imageHeight) continue; - // std::cout << xr << ", " << yr << std::endl; - // Point cp{xr, yr}; - // checkedPoints.insert(cp); - checkedPoints.insert({xr, yr}); - } - } + xb = betaX(1); + yb = betaY(1); + a = betaX(0); + sigma = betaX(2); - std::cout << "Original: " << x << ", " << y << std::endl; - std::cout << "final: " << xb << ", " << yb << std::endl; - // result.push_back(Star(xb, yb, sigma, sigma, a)); - result.push_back(Star(xb, yb, 0)); - // result.push_back(Star(x, y, 0)); - } + // std::cout << "Original: " << x << ", " << y << std::endl; + // std::cout << "final: " << xb << ", " << yb << std::endl; + result.push_back(Star(xb, yb, 0)); } std::cout << "Number of centroids: " << result.size() << std::endl; From 08d1d15f0d89666917d548ca84b56f83d44a48ef Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 11 Mar 2023 17:51:50 -0800 Subject: [PATCH 10/34] Unit test Gaussian helper functions --- Makefile | 2 +- src/centroiders.cpp | 941 ++++++++++++++++++------------------- src/centroiders.hpp | 23 +- test/gaussian-centroid.cpp | 78 +++ 4 files changed, 568 insertions(+), 476 deletions(-) create mode 100644 test/gaussian-centroid.cpp diff --git a/Makefile b/Makefile index 8e364f87..147f5630 100644 --- a/Makefile +++ b/Makefile @@ -85,7 +85,7 @@ lint: test: $(BIN) $(BSC) $(TEST_BIN) $(TEST_BIN) - bash ./test/scripts/pyramid-incorrect.sh + # bash ./test/scripts/pyramid-incorrect.sh # bash ./test/scripts/readme-examples-test.sh # bash ./test/scripts/tetra.sh diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 372d7be4..e8869736 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -20,87 +20,87 @@ namespace lost { typedef std::vector Point; -int Get(int x, int y, const unsigned char *image, int w); +// int Get(int x, int y, const unsigned char *image, int w); std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight); int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight); std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight) { - std::vector res; - std::set checkedPoints; - - int cutoff = BasicThreshold(image, imageWidth, imageHeight); - - 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; - - int maxMag = -1; - int x0 = x; - int y0 = y; - - queue.push_back(pCurr); - while (!queue.size() == 0) { - Point p = queue[0]; - queue.pop_front(); - - // TODO: make Point a struct probably - 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; - - // 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; + std::vector res; + std::set checkedPoints; + + int cutoff = BasicThreshold(image, imageWidth, imageHeight); + + 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; + + int maxMag = -1; + int x0 = x; + int y0 = y; + + queue.push_back(pCurr); + while (!queue.size() == 0) { + Point p = queue[0]; + queue.pop_front(); + + // TODO: make Point a struct probably + 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; + + // 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}); } - - // 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}); } - } - std::cout << "Floodfill preprocessing: found " << res.size() << " total points" << std::endl; - return res; + 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; + float yCoordMagSum; + float xCoordMagSum; + long magSum; + int xMin; + int xMax; + int yMin; + int yMax; + int cutoff; + bool isValid; + std::unordered_set checkedIndices; }; /* @@ -108,379 +108,370 @@ Get pixel value at (x, y) where x and y are 0-based Note: top left is (0, 0) */ int Get(int x, int y, const unsigned char *image, int w) { - int ind = y * w + x; - return image[ind]; + int ind = y * w + x; + return image[ind]; } -/* -Get XMarginal - 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) { - int sum = 0; - for (int j = -nb; j <= nb; j++) { - sum += Get(x, y0 + j, image, w); - } - return sum; + int sum = 0; + for (int j = -nb; j <= nb; j++) { + sum += Get(x, y0 + j, image, w); + } + return sum; } -/* -Get YMarginal - 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) { - int sum = 0; - for (int i = -nb; i <= nb; i++) { - sum += Get(x0 + i, y, image, w); - } - return sum; + int sum = 0; + for (int i = -nb; i <= nb; i++) { + sum += Get(x0 + i, y, image, w); + } + return sum; } /* Can be used as f(xi, beta) or f(yi, beta) */ -float FitModel(float xi, float a, float xb, float sigma) { - return a * exp(-1 * (xi - xb) * (xi - xb) / (2 * sigma * sigma)); +float FitModel(float x, float a, float xb, float sigma) { + return a * exp(-1 * (x - xb) * (x - xb) / (2 * sigma * sigma)); } -/* -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) { - // 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 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++; - } + *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++; + } + } } - } - double fwhm = std::sqrt(halfCount); - *sigma = fwhm / (2 * std::sqrt(2 * std::log(2))); + double fwhm = std::sqrt(halfCount); + *sigma = 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; + typedef _Scalar Scalar; + enum { InputsAtCompileTime = NX, ValuesAtCompileTime = NY }; + typedef Eigen::Matrix InputType; + typedef Eigen::Matrix ValueType; + typedef Eigen::Matrix JacobianType; - int m_inputs, m_values; + int m_inputs, m_values; - Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} - Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} + 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; } + int inputs() const { return m_inputs; } + int values() const { return m_values; } }; struct LSGFFunctor : Functor { - // First param = number of parameters in your model - // Second param = number of data points you want to test over - // for us, = 2 * nb + 1 - LSGFFunctor(const int nb, const int marg, Eigen::VectorXd 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) {} - - /* - x = parameters (a, xb, sigma) - fvec = residual (one for every data point, of size = 2*nb+1) - */ - int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - float a = x(0); - float xb = x(1); - double 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) = marginal - FitModel(X(i), a, xb, sigma); + // First param = number of parameters in your model + // Second param = number of data points you want to test over + // for us, = 2 * nb + 1 + LSGFFunctor(const int nb, const int marg, Eigen::VectorXd 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) {} + + /* + x = parameters (a, xb, sigma) + fvec = residual (one for every data point, of size = 2*nb+1) + TODO: this isn't great?? We're just fitting 5 data points + */ + int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { + float a = x(0); + float xb = x(1); + double 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) = marginal - FitModel(X(i), a, xb, sigma); + } + std::cout << fvec(0) << std::endl; + + return 0; } - std::cout << fvec(0) << std::endl; - - 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::VectorXd X; - const unsigned char *image; - const int w; // image width in pixels - // Center coordinates (x0, y0) of this window - const int x0, y0; + + // 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::VectorXd X; + 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::vector result; - const int nb = 2; + const int nb = 2; - // std::set checkedPoints; + // std::set checkedPoints; - std::vector candidatePts = FloodfillPreproc(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; + 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; - Eigen::VectorXd X(2 * nb + 1); - Eigen::VectorXd Y(2 * nb + 1); - int vInd = 0; - for (int i = -nb; i <= nb; i++) { - X(vInd) = x + i; - Y(vInd) = y + i; - vInd++; - } + Eigen::VectorXd X(2 * nb + 1); + Eigen::VectorXd Y(2 * nb + 1); + int vInd = 0; + for (int i = -nb; i <= nb; i++) { + X(vInd) = x + i; + Y(vInd) = y + i; + vInd++; + } - float a; - float xb, yb; - double sigma; - InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); + float a; + float xb, yb; + double sigma; + InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); - Eigen::VectorXd betaX(3); - betaX << a, xb, sigma; + Eigen::VectorXd betaX(3); + betaX << a, xb, sigma; - Eigen::VectorXd betaY(3); - betaY << a, yb, sigma; + Eigen::VectorXd betaY(3); + betaY << a, yb, sigma; - LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffX(functorX); - Eigen::LevenbergMarquardt, double> lmX(numDiffX); + LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffX(functorX); + Eigen::LevenbergMarquardt, double> lmX(numDiffX); - lmX.parameters.maxfev = 2000; - lmX.parameters.xtol = 1.0e-10; + lmX.parameters.maxfev = 2000; + lmX.parameters.xtol = 1.0e-10; - LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffY(functorY); - Eigen::LevenbergMarquardt, double> lmY(numDiffY); + LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); + Eigen::NumericalDiff numDiffY(functorY); + Eigen::LevenbergMarquardt, double> lmY(numDiffY); - lmY.parameters.maxfev = 2000; - lmY.parameters.xtol = 1.0e-10; + lmY.parameters.maxfev = 2000; + lmY.parameters.xtol = 1.0e-10; - lmX.minimize(betaX); - lmY.minimize(betaY); + lmX.minimize(betaX); + lmY.minimize(betaY); - xb = betaX(1); - yb = betaY(1); - a = betaX(0); - sigma = betaX(2); + xb = betaX(1); + yb = betaY(1); + a = betaX(0); + sigma = betaX(2); - // std::cout << "Original: " << x << ", " << y << std::endl; - // std::cout << "final: " << xb << ", " << yb << std::endl; - result.push_back(Star(xb, yb, 0)); - } + // std::cout << "Original: " << x << ", " << y << std::endl; + // std::cout << "final: " << xb << ", " << yb << std::endl; + result.push_back(Star(xb, yb, 0)); + // result.push_back(Star(x, y, 0)); + } - std::cout << "Number of centroids: " << result.size() << std::endl; + std::cout << "Number of centroids: " << result.size() << std::endl; - return result; + return result; } // DUMMY std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, int imageHeight) const { - std::vector result; + 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)); - } + 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)); + } - return result; + return result; } // a poorly designed thresholding algorithm int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) { - // 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; + // 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; } // 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 sumB = 0; - float sum1 = 0; - float wB = 0; - float maximum = 0; - int level = 0; - // make the histogram (array length 256) - int histogram[256]; - - memset(histogram, 0, sizeof(int) * 256); - - for (long i = 0; i < total; i++) { - histogram[image[i]]++; - } - for (int i = 0; i < 256; i++) { - sum1 += i * histogram[i]; - } - for (int i = 0; i < 256; i++) { - float 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); - // std::cout << val << "\n"; - if (val >= maximum) { - level = i; - maximum = val; - } + // code here, duh + long total = imageWidth * imageHeight; + // float top = 255; + float sumB = 0; + float sum1 = 0; + float wB = 0; + float maximum = 0; + int level = 0; + // make the histogram (array length 256) + int histogram[256]; + + memset(histogram, 0, sizeof(int) * 256); + + for (long i = 0; i < total; i++) { + histogram[image[i]]++; } - wB = wB + histogram[i]; - sumB = sumB + i * histogram[i]; - } - return level; + for (int i = 0; i < 256; i++) { + sum1 += i * histogram[i]; + } + for (int i = 0; i < 256; i++) { + float 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); + // std::cout << val << "\n"; + if (val >= maximum) { + level = i; + maximum = val; + } + } + wB = wB + histogram[i]; + sumB = sumB + i * histogram[i]; + } + 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; - long totalPixels = imageHeight * imageWidth; - for (long i = 0; i < totalPixels; i++) { - totalMag += image[i]; - } - float mean = totalMag / totalPixels; - for (long i = 0; i < totalPixels; i++) { - std += std::pow(image[i] - mean, 2); - } - std = std::sqrt(std / totalPixels); - return mean + (std * 5); + unsigned long totalMag = 0; + float std = 0; + long totalPixels = imageHeight * imageWidth; + for (long i = 0; i < totalPixels; i++) { + totalMag += image[i]; + } + float mean = totalMag / totalPixels; + for (long i = 0; i < totalPixels; i++) { + std += std::pow(image[i] - mean, 2); + } + std = std::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; - 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); - return mean + (std * 5); + unsigned long totalMag = 0; + float std = 0; + float 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); + return mean + (std * 5); } // 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) { - p->isValid = false; - } - p->checkedIndices.insert(i); - if (i % imageWidth > p->xMax) { - p->xMax = i % imageWidth; - } else if (i % imageWidth < p->xMin) { - p->xMin = i % imageWidth; - } - if (i / imageWidth > p->yMax) { - p->yMax = i / imageWidth; - } else if (i / imageWidth < p->yMin) { - p->yMin = i / imageWidth; - } - p->magSum += image[i]; - p->xCoordMagSum += ((i % imageWidth)) * image[i]; - p->yCoordMagSum += ((i / imageWidth)) * image[i]; - if (i % imageWidth != imageWidth - 1) { - CogHelper(p, i + 1, image, imageWidth, imageHeight); - } - if (i % imageWidth != 0) { - CogHelper(p, i - 1, image, imageWidth, 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) { + p->isValid = false; + } + p->checkedIndices.insert(i); + if (i % imageWidth > p->xMax) { + p->xMax = i % imageWidth; + } else if (i % imageWidth < p->xMin) { + p->xMin = i % imageWidth; + } + if (i / imageWidth > p->yMax) { + p->yMax = i / imageWidth; + } else if (i / imageWidth < p->yMin) { + p->yMin = i / imageWidth; + } + p->magSum += image[i]; + p->xCoordMagSum += ((i % imageWidth)) * image[i]; + p->yCoordMagSum += ((i / imageWidth)) * image[i]; + if (i % imageWidth != imageWidth - 1) { + CogHelper(p, i + 1, image, imageWidth, imageHeight); + } + if (i % imageWidth != 0) { + CogHelper(p, i - 1, image, imageWidth, imageHeight); + } + CogHelper(p, i + imageWidth, image, imageWidth, imageHeight); + CogHelper(p, i - imageWidth, image, imageWidth, imageHeight); } - CogHelper(p, i + imageWidth, image, imageWidth, imageHeight); - CogHelper(p, i - imageWidth, image, imageWidth, imageHeight); - } } std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { - CentroidParams p; - - std::vector result; - - 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 - 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.xMax = i % imageWidth; - p.xMin = i % imageWidth; - p.yMax = i / imageWidth; - p.yMin = i / imageWidth; - p.isValid = true; - - int sizeBefore = p.checkedIndices.size(); - - CogHelper(&p, i, image, imageWidth, imageHeight); - 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 - 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)); - } + CentroidParams p; + + std::vector result; + + 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 + 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.xMax = i % imageWidth; + p.xMin = i % imageWidth; + p.yMax = i / imageWidth; + p.yMin = i / imageWidth; + p.isValid = true; + + int sizeBefore = p.checkedIndices.size(); + + CogHelper(&p, i, image, imageWidth, imageHeight); + 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 + 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)); + } + } } - } - return result; + return result; } // Determines how accurate and how much iteration is done by the IWCoG algorithm, @@ -488,134 +479,136 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi float iWCoGMinChange = 0.0002; struct IWCoGParams { - int xMin; - int xMax; - int yMin; - int yMax; - int cutoff; - int maxIntensity; - int guess; - bool isValid; - std::unordered_set checkedIndices; + int xMin; + int xMax; + int yMin; + int yMax; + int cutoff; + int maxIntensity; + int guess; + bool isValid; + 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) { - p->isValid = false; - } - p->checkedIndices.insert(i); - starIndices->push_back(i); - if (image[i] > p->maxIntensity) { - p->maxIntensity = image[i]; - p->guess = i; - } - if (i % imageWidth > p->xMax) { - p->xMax = i % imageWidth; - } else if (i % imageWidth < p->xMin) { - p->xMin = i % imageWidth; - } - if (i / imageWidth > p->yMax) { - p->yMax = i / imageWidth; - } else if (i / imageWidth < p->yMin) { - p->yMin = i / imageWidth; - } - if (i % imageWidth != imageWidth - 1) { - IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, starIndices); - } - if (i % imageWidth != 0) { - IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, 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); + starIndices->push_back(i); + if (image[i] > p->maxIntensity) { + p->maxIntensity = image[i]; + p->guess = i; + } + if (i % imageWidth > p->xMax) { + p->xMax = i % imageWidth; + } else if (i % imageWidth < p->xMin) { + p->xMin = i % imageWidth; + } + if (i / imageWidth > p->yMax) { + p->yMax = i / imageWidth; + } else if (i / imageWidth < p->yMin) { + p->yMin = i / imageWidth; + } + if (i % imageWidth != imageWidth - 1) { + IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, starIndices); + } + if (i % imageWidth != 0) { + IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, starIndices); + } + IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, starIndices); } - IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, starIndices); - IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, starIndices); - } } Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { - IWCoGParams p; - std::vector result; - 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 - if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { - // TODO: store longs --Mark - 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 standardDeviation; - float w; // weight value - - p.xMax = i % imageWidth; - p.xMin = i % imageWidth; - p.yMax = i / imageWidth; - 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 - float 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); - // 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 (change > iWCoGMinChange && stop < 100000) { - // 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)]); + IWCoGParams p; + std::vector result; + 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 + if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { + // TODO: store longs --Mark + 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 standardDeviation; + float w; // weight value + + p.xMax = i % imageWidth; + p.xMin = i % imageWidth; + p.yMax = i / imageWidth; + 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 + float 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); + // 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 (change > iWCoGMinChange && stop < 100000) { + // 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)]); + } + float xTemp = xWeightedCoordMagSum / weightedMagSum; + float yTemp = yWeightedCoordMagSum / weightedMagSum; + + change = abs(guessXCoord - xTemp) + abs(guessYCoord - yTemp); + + guessXCoord = xTemp; + 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())); + } } - float xTemp = xWeightedCoordMagSum / weightedMagSum; - float yTemp = yWeightedCoordMagSum / weightedMagSum; - - change = abs(guessXCoord - xTemp) + abs(guessYCoord - yTemp); - - guessXCoord = xTemp; - 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())); - } } - } - return result; + return result; } } // namespace lost diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 14cb2640..7721b7e6 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -62,6 +62,27 @@ class IterativeWeightedCenterOfGravityAlgorithm : public CentroidAlgorithm { Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; }; -} +// TODO: a bunch of functions that really should be private, but tests need to access + +/// 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 XMarginal - 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 YMarginal - 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); + +/* +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/test/gaussian-centroid.cpp b/test/gaussian-centroid.cpp new file mode 100644 index 00000000..a91cdb81 --- /dev/null +++ b/test/gaussian-centroid.cpp @@ -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 From 11b48bb1c0a604a4c12c22e3e55e767df5504256 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 11 Mar 2023 19:09:27 -0800 Subject: [PATCH 11/34] First try of 2D gaussian fit - still not good enough --- src/centroiders.cpp | 131 ++++++++++++++++++++++++++++++++++++++++++-- src/centroiders.hpp | 9 +++ src/io.cpp | 3 + 3 files changed, 139 insertions(+), 4 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index e8869736..eba18e93 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -135,6 +135,13 @@ float FitModel(float x, float a, float xb, float sigma) { return a * exp(-1 * (x - xb) * (x - xb) / (2 * sigma * sigma)); } +/// +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; +} + void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, float *xb, float *yb, double *sigma) { @@ -206,8 +213,8 @@ struct LSGFFunctor : Functor { TODO: this isn't great?? We're just fitting 5 data points */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - float a = x(0); - float xb = x(1); + double a = x(0); + double xb = x(1); double sigma = x(2); for (int i = 0; i < X.size(); i++) { int marginal; @@ -217,7 +224,6 @@ struct LSGFFunctor : Functor { marginal = YMarginal(x0, X(i), nb, image, w); fvec(i) = marginal - FitModel(X(i), a, xb, sigma); } - std::cout << fvec(0) << std::endl; return 0; } @@ -234,10 +240,70 @@ struct LSGFFunctor : Functor { const int x0, y0; }; +/// Functor for 2D Least-squares Gaussian Fit algo +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::VectorXd &x, Eigen::VectorXd &fvec) const { + double a = x(0); + double xb = x(1); + double yb = x(2); + double sigmaX = x(3); + double sigmaY = x(4); + + // 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) = marginal - FitModel(X(i), a, xb, sigma); + // } + + 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 = Get(xi, yi, image, w); + float modelPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); + fvec(ind) = yPred - modelPred; + + 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" << std::endl; + const int nb = 2; // std::set checkedPoints; @@ -286,9 +352,11 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW lmX.minimize(betaX); lmY.minimize(betaY); + a = betaX(0); + xb = betaX(1); yb = betaY(1); - a = betaX(0); + sigma = betaX(2); // std::cout << "Original: " << x << ", " << y << std::endl; @@ -302,6 +370,61 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW return result; } +std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageWidth, + int imageHeight) const { + std::vector result; + + std::cout << "2D GAUSSIAN FIT" << std::endl; + + const int nb = 2; + const int np = 2 * nb + 1; + + // std::set checkedPoints; + + 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 a; + float xb, yb; + double sigmaX, sigmaY; + InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigmaX); + + sigmaY = sigmaX; + + Eigen::VectorXd beta(5); + beta << a, xb, yb, sigmaX, sigmaY; + + // LSGF2DFunctor functor(nb, 0, X, image, imageWidth, x, y); + LSGF2DFunctor functor(nb, image, imageWidth, x, y); + Eigen::NumericalDiff numDiff(functor); + Eigen::LevenbergMarquardt, double> lm(numDiff); + + lm.parameters.maxfev = 2000; + lm.parameters.xtol = 1.0e-10; + + lm.minimize(beta); + + a = beta(0); + xb = beta(1); + yb = beta(2); + sigmaX = beta(3); + sigmaY = beta(4); + + // std::cout << "Original: " << x << ", " << y << std::endl; + // std::cout << "final: " << xb << ", " << yb << std::endl; + result.push_back(Star(xb, yb, 0)); + // result.push_back(Star(x, y, 0)); + } + + std::cout << "Number of centroids: " << result.size() << std::endl; + + return result; +} + // DUMMY std::vector DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, int imageHeight) const { diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 7721b7e6..4fe7e321 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -32,6 +32,15 @@ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; }; +/** + * @brief Least Squares Gaussian Fit 2D centroiding algorithm + * + */ +class LeastSquaresGaussianFit2D : public CentroidAlgorithm { + public: + explicit LeastSquaresGaussianFit2D(){}; + Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; +}; /// A centroid algorithm for debugging that returns random centroids. class DummyCentroidAlgorithm: public CentroidAlgorithm { diff --git a/src/io.cpp b/src/io.cpp index 7c317fab..527d7698 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -724,6 +724,9 @@ Pipeline SetPipeline(const PipelineOptions &values) { } else if (values.centroidAlgo == "lsgf1d") { result.centroidAlgorithm = std::unique_ptr(new LeastSquaresGaussianFit1D()); + } else if (values.centroidAlgo == "lsgf2d") { + result.centroidAlgorithm = + std::unique_ptr(new LeastSquaresGaussianFit2D()); } else if (values.centroidAlgo != "") { std::cout << "Illegal centroid algorithm." << std::endl; exit(1); From 7d308c6361879f559b828b72ab3a0540ed6426f7 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Wed, 15 Mar 2023 10:28:06 -0700 Subject: [PATCH 12/34] smol syntax change --- src/centroiders.cpp | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index eba18e93..ba098f03 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -31,6 +31,7 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im std::set checkedPoints; int cutoff = BasicThreshold(image, imageWidth, imageHeight); + std::cout << "cutoff: " << cutoff << std::endl; for (long i = 0; i < imageHeight * imageWidth; i++) { int x = i % imageWidth; @@ -47,7 +48,7 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im int y0 = y; queue.push_back(pCurr); - while (!queue.size() == 0) { + while (queue.size() != 0) { Point p = queue[0]; queue.pop_front(); @@ -264,20 +265,12 @@ struct LSGF2DFunctor : Functor { double sigmaX = x(3); double sigmaY = x(4); - // 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) = marginal - FitModel(X(i), a, xb, sigma); - // } - 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; + // TODO: other way around, yPred is our modelPred, yActual is our pixel intensity float yPred = Get(xi, yi, image, w); float modelPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); fvec(ind) = yPred - modelPred; @@ -383,6 +376,9 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); + // TODO: remove, testing + int same = 0; + for (const Point &pt : candidatePts) { int x = pt[0]; int y = pt[1]; @@ -415,12 +411,16 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW sigmaY = beta(4); // std::cout << "Original: " << x << ", " << y << std::endl; - // std::cout << "final: " << xb << ", " << yb << std::endl; - result.push_back(Star(xb, yb, 0)); - // result.push_back(Star(x, y, 0)); + std::cout << xb << ", " << yb << std::endl; + if(x == xb && y == yb) same++; + result.push_back(Star(xb, yb, nb)); + + + // result.push_back(Star(x, y, nb)); } std::cout << "Number of centroids: " << result.size() << std::endl; + std::cout << "same: " << same << std::endl; return result; } @@ -591,6 +591,8 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f, ((float)(xDiameter)) / 2.0f, ((float)(yDiameter)) / 2.0f, p.checkedIndices.size() - sizeBefore)); + + // std::cout << result.back().position.x << ", " << result.back().position.y << std::endl; } } } From 14e9dc27c20dbe6ad60e0761918b4437d169fce5 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 18 Mar 2023 11:54:49 -0700 Subject: [PATCH 13/34] more debugging LM --- src/centroiders.cpp | 46 ++++++++++++++++++++++++++------------------- src/centroiders.hpp | 4 ++++ 2 files changed, 31 insertions(+), 19 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index ba098f03..9771be77 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -20,8 +20,6 @@ namespace lost { typedef std::vector Point; -// int Get(int x, int y, const unsigned char *image, int w); - std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight); int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight); @@ -175,6 +173,19 @@ void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int *sigma = fwhm / (2 * std::sqrt(2 * std::log(2))); } +float InitialGuess2(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 { @@ -271,9 +282,9 @@ struct LSGF2DFunctor : Functor { int xi = x0 + i; int yi = y0 + j; // TODO: other way around, yPred is our modelPred, yActual is our pixel intensity - float yPred = Get(xi, yi, image, w); - float modelPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); - fvec(ind) = yPred - modelPred; + float yPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); + float yActual = Get(xi, yi, image, w); + fvec(ind) = yPred - yActual; ind++; } @@ -384,15 +395,12 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW int y = pt[1]; if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; - float a; - float xb, yb; - double sigmaX, sigmaY; - InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigmaX); - - sigmaY = sigmaX; + float a = Get(x, y, image, imageWidth); + double sigma = InitialGuess2(x, y, a, nb, image, imageWidth); + // InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigmaX); Eigen::VectorXd beta(5); - beta << a, xb, yb, sigmaX, sigmaY; + beta << a, x, y, sigma, sigma; // LSGF2DFunctor functor(nb, 0, X, image, imageWidth, x, y); LSGF2DFunctor functor(nb, image, imageWidth, x, y); @@ -405,15 +413,15 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW lm.minimize(beta); a = beta(0); - xb = beta(1); - yb = beta(2); - sigmaX = beta(3); - sigmaY = beta(4); + float xRes = beta(1); + float yRes = beta(2); + float sigmaX = beta(3); + float sigmaY = beta(4); // std::cout << "Original: " << x << ", " << y << std::endl; - std::cout << xb << ", " << yb << std::endl; - if(x == xb && y == yb) same++; - result.push_back(Star(xb, yb, nb)); + std::cout << xRes << ", " << yRes << ": " << a << " " << sigmaX << std::endl; + if(x == xRes && y == yRes) same++; + result.push_back(Star(xRes, yRes, nb)); // result.push_back(Star(x, y, nb)); diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 4fe7e321..5fc24bb4 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -92,6 +92,10 @@ 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); +// TODO: rename +// Basically initial guess for window, done after Floodfill preprocessing +float InitialGuess2(int x0, int y0, int maxMag, const int nb, const unsigned char *image, int w); + } // namespace lost #endif From fa29b4bb8be6c8457538a3c6dcf6e5bdf2babaee Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Tue, 21 Mar 2023 22:44:24 -0700 Subject: [PATCH 14/34] sus Jacobian math --- src/centroiders.cpp | 43 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 9771be77..eec22ca8 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -222,7 +222,6 @@ struct LSGFFunctor : Functor { /* x = parameters (a, xb, sigma) fvec = residual (one for every data point, of size = 2*nb+1) - TODO: this isn't great?? We're just fitting 5 data points */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { double a = x(0); @@ -293,6 +292,37 @@ struct LSGF2DFunctor : Functor { return 0; } + int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const{ + double a = x(0); + double xb = x(1); + double yb = x(2); + double sigmaX = x(3); + double sigmaY = x(4); + + for (int i = -nb; i <= nb; i++) { + for (int j = -nb; j <= nb; j++) { + int xi = x0 + i; + int yj = y0 + j; + + float yActual = Get(xi, yj, image, w); + float yPred = FitModel2D(xi, yj, a, xb, yb, sigmaX, sigmaY); + float ei = yActual - yPred; + + float expX = exp(-1 * (xi - xb) * (xi - xb) / (2 * sigmaX * sigmaX)); + float expY = exp(-1 * (yj - yb) * (yj - yb) / (2 * sigmaY * sigmaY)); + + int row = i + nb; + fjac(row, 0) = -expX * expY; + fjac(row, 1) = -a * expY * (xi - xb) / (sigmaX * sigmaX) * expX; + fjac(row, 2) = -a * expX * (yj - yb) / (sigmaY * sigmaY) * expY; + fjac(row, 3) = -a * expX * expY * std::pow(xi - xb, 2) / (std::pow(sigmaX, 3)); + fjac(row, 4) = -a * expX * expY * std::pow(yj - yb, 2) / (std::pow(sigmaY, 3)); + } + } + + return 0; + } + // Window size is (2*nb + 1) x (2*nb + 1) const int nb; @@ -404,8 +434,9 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW // LSGF2DFunctor functor(nb, 0, X, image, imageWidth, x, y); LSGF2DFunctor functor(nb, image, imageWidth, x, y); - Eigen::NumericalDiff numDiff(functor); - Eigen::LevenbergMarquardt, double> lm(numDiff); + // Eigen::NumericalDiff numDiff(functor); + // Eigen::LevenbergMarquardt, double> lm(numDiff); + Eigen::LevenbergMarquardt lm(functor); lm.parameters.maxfev = 2000; lm.parameters.xtol = 1.0e-10; @@ -419,9 +450,9 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW float sigmaY = beta(4); // std::cout << "Original: " << x << ", " << y << std::endl; - std::cout << xRes << ", " << yRes << ": " << a << " " << sigmaX << std::endl; - if(x == xRes && y == yRes) same++; - result.push_back(Star(xRes, yRes, nb)); + std::cout << xRes << ", " << yRes << std::endl; + // if(x == xRes && y == yRes) same++; + result.push_back(Star(xRes+0.5, yRes+0.5, 0)); // result.push_back(Star(x, y, nb)); From 78cba87a24aa42865d8281ce59a4508bf50f8a86 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Thu, 23 Mar 2023 14:44:04 -0700 Subject: [PATCH 15/34] LM working for 2D fit, implemented calculation of Jacobian --- src/centroiders.cpp | 43 ++++++++++++++++++------------------------- 1 file changed, 18 insertions(+), 25 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index eec22ca8..981026c2 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -299,24 +299,23 @@ struct LSGF2DFunctor : Functor { double sigmaX = x(3); double 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 yActual = Get(xi, yj, image, w); - float yPred = FitModel2D(xi, yj, a, xb, yb, sigmaX, sigmaY); - float ei = yActual - yPred; - float expX = exp(-1 * (xi - xb) * (xi - xb) / (2 * sigmaX * sigmaX)); float expY = exp(-1 * (yj - yb) * (yj - yb) / (2 * sigmaY * sigmaY)); - int row = i + nb; - fjac(row, 0) = -expX * expY; - fjac(row, 1) = -a * expY * (xi - xb) / (sigmaX * sigmaX) * expX; - fjac(row, 2) = -a * expX * (yj - yb) / (sigmaY * sigmaY) * expY; - fjac(row, 3) = -a * expX * expY * std::pow(xi - xb, 2) / (std::pow(sigmaX, 3)); - fjac(row, 4) = -a * expX * expY * std::pow(yj - yb, 2) / (std::pow(sigmaY, 3)); + 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++; } } @@ -417,9 +416,6 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); - // TODO: remove, testing - int same = 0; - for (const Point &pt : candidatePts) { int x = pt[0]; int y = pt[1]; @@ -443,23 +439,20 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW lm.minimize(beta); - a = beta(0); - float xRes = beta(1); - float yRes = beta(2); - float sigmaX = beta(3); - float sigmaY = beta(4); - - // std::cout << "Original: " << x << ", " << y << std::endl; - std::cout << xRes << ", " << yRes << std::endl; - // if(x == xRes && y == yRes) same++; - result.push_back(Star(xRes+0.5, yRes+0.5, 0)); + double aRes = beta(0); + double xRes = beta(1); + double yRes = beta(2); + double sigmaX = beta(3); + double sigmaY = beta(4); + // std::cout << xRes+0.5 << ", " << yRes+0.5 << std::endl; + // result.push_back(Star(xRes+0.5, yRes+0.5, 0)); + // TODO: not sure how much we care about making the radius/brightness accurate + result.push_back(Star(xRes+0.5, yRes+0.5, sigmaX, sigmaY, aRes)); - // result.push_back(Star(x, y, nb)); } std::cout << "Number of centroids: " << result.size() << std::endl; - std::cout << "same: " << same << std::endl; return result; } From 6ae68d2c2170b9a2684a24d2cd1355cef0dda365 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Thu, 23 Mar 2023 19:18:47 -0700 Subject: [PATCH 16/34] Calculate Jacobian for 1D fit --- src/centroiders.cpp | 41 +++++++++++++++++++++++++++++++++++------ 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 981026c2..a42fe4d9 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -233,7 +233,29 @@ struct LSGFFunctor : Functor { marginal = XMarginal(X(i), y0, nb, image, w); else marginal = YMarginal(x0, X(i), nb, image, w); - fvec(i) = marginal - FitModel(X(i), a, xb, sigma); + fvec(i) = FitModel(X(i), a, xb, sigma) - marginal; + } + + return 0; + } + + int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const { + double a = x(0); + double kb = x(1); + double 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; @@ -368,16 +390,23 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW Eigen::VectorXd betaY(3); betaY << a, yb, sigma; + // LSGF2DFunctor functor(nb, image, imageWidth, x, y); + // // Eigen::NumericalDiff numDiff(functor); + // // Eigen::LevenbergMarquardt, double> lm(numDiff); + // Eigen::LevenbergMarquardt lm(functor); + LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffX(functorX); - Eigen::LevenbergMarquardt, double> lmX(numDiffX); + // Eigen::NumericalDiff numDiffX(functorX); + Eigen::LevenbergMarquardt lmX(functorX); + lmX.parameters.maxfev = 2000; lmX.parameters.xtol = 1.0e-10; LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); - Eigen::NumericalDiff numDiffY(functorY); - Eigen::LevenbergMarquardt, double> lmY(numDiffY); + // Eigen::NumericalDiff numDiffY(functorY); + // Eigen::LevenbergMarquardt, double> lmY(numDiffY); + Eigen::LevenbergMarquardt lmY(functorY); lmY.parameters.maxfev = 2000; lmY.parameters.xtol = 1.0e-10; @@ -394,7 +423,7 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW // std::cout << "Original: " << x << ", " << y << std::endl; // std::cout << "final: " << xb << ", " << yb << std::endl; - result.push_back(Star(xb, yb, 0)); + result.push_back(Star(xb+0.5, yb+0.5, 0)); // result.push_back(Star(x, y, 0)); } From bbcc314fe2a4f7974c7a3794a62e4e74102618e8 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Thu, 23 Mar 2023 19:30:32 -0700 Subject: [PATCH 17/34] make 1D and 2D fit use same preprocessing --- src/centroiders.cpp | 19 +++++++++++-------- src/centroiders.hpp | 19 ++++++++++++------- 2 files changed, 23 insertions(+), 15 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index a42fe4d9..f25e95a7 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -379,16 +379,19 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW vInd++; } - float a; - float xb, yb; - double sigma; - InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); + // float a; + // float xb, yb; + // double sigma; + // InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); + + double a = Get(x, y, image, imageWidth); + double sigma = InitialGuess2(x, y, a, nb, image, imageWidth); Eigen::VectorXd betaX(3); - betaX << a, xb, sigma; + betaX << a, x, sigma; Eigen::VectorXd betaY(3); - betaY << a, yb, sigma; + betaY << a, y, sigma; // LSGF2DFunctor functor(nb, image, imageWidth, x, y); // // Eigen::NumericalDiff numDiff(functor); @@ -416,8 +419,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW a = betaX(0); - xb = betaX(1); - yb = betaY(1); + double xb = betaX(1); + double yb = betaY(1); sigma = betaX(2); diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 5fc24bb4..e4b7c21d 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -22,9 +22,11 @@ class CentroidAlgorithm { }; /** - * @brief Least Squares Gaussian Fit 1D centroiding algorithm - * Uses Levenberg-Marquardt to solve nonlinear least-squares + * @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: @@ -33,8 +35,11 @@ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ }; /** - * @brief Least Squares Gaussian Fit 2D centroiding algorithm + * @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: @@ -76,10 +81,10 @@ class IterativeWeightedCenterOfGravityAlgorithm : public CentroidAlgorithm { /// 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 XMarginal - keep x fixed, sum pixel values from [y0-nb, y0+nb] +/// 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 YMarginal - keep y fixed, sum pixel values from [x0-nb, x0+nb] +/// 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); /* @@ -89,8 +94,8 @@ 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); +// void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, +// float *xb, float *yb, double *sigma); // TODO: rename // Basically initial guess for window, done after Floodfill preprocessing From c814be85d17e9124b3f737ec83eb7c7489fb8218 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Fri, 24 Mar 2023 11:21:20 -0700 Subject: [PATCH 18/34] code cleanup --- src/centroiders.cpp | 162 +++++++++++++++++++------------------------- src/centroiders.hpp | 15 ++-- 2 files changed, 79 insertions(+), 98 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index f25e95a7..3817a07a 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -18,10 +18,25 @@ namespace lost { +/** +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); + 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); std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight) { @@ -29,7 +44,6 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im std::set checkedPoints; int cutoff = BasicThreshold(image, imageWidth, imageHeight); - std::cout << "cutoff: " << cutoff << std::endl; for (long i = 0; i < imageHeight * imageWidth; i++) { int x = i % imageWidth; @@ -50,7 +64,6 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im Point p = queue[0]; queue.pop_front(); - // TODO: make Point a struct probably int px = p[0]; int py = p[1]; if (px < 0 || px >= imageWidth || py < 0 || py >= imageHeight) continue; @@ -84,7 +97,7 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im } } - std::cout << "Floodfill preprocessing: found " << res.size() << " total points" << std::endl; + // std::cout << "Floodfill preprocessing: found " << res.size() << " total points" << std::endl; return res; } @@ -102,10 +115,6 @@ struct CentroidParams { std::unordered_set checkedIndices; }; -/* -Get pixel value at (x, y) where x and y are 0-based -Note: top left is (0, 0) -*/ int Get(int x, int y, const unsigned char *image, int w) { int ind = y * w + x; return image[ind]; @@ -127,53 +136,49 @@ int YMarginal(int x0, int y, int nb, const unsigned char *image, int w) { return sum; } -/* -Can be used as f(xi, beta) or f(yi, beta) -*/ -float FitModel(float x, float a, float xb, float sigma) { +static float FitModel(float x, float a, float xb, float sigma) { return a * exp(-1 * (x - xb) * (x - xb) / (2 * sigma * sigma)); } -/// -float FitModel2D(float x, float y, float a, float xb, float yb, float sigmaX, float sigmaY){ +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; } - -void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, - float *xb, float *yb, double *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++; - } - } - } - double fwhm = std::sqrt(halfCount); - *sigma = fwhm / (2 * std::sqrt(2 * std::log(2))); -} - -float InitialGuess2(int x0, int y0, int maxMag, const int nb, const unsigned char *image, int w){ +// DO NO DELETE +// void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, +// float *xb, float *yb, double *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++; +// } +// } +// } +// double 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++){ @@ -195,7 +200,8 @@ struct Functor { typedef Eigen::Matrix ValueType; typedef Eigen::Matrix JacobianType; - int m_inputs, m_values; + 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) {} @@ -204,11 +210,9 @@ struct Functor { int values() const { return m_values; } }; -struct LSGFFunctor : Functor { - // First param = number of parameters in your model - // Second param = number of data points you want to test over - // for us, = 2 * nb + 1 - LSGFFunctor(const int nb, const int marg, Eigen::VectorXd X, const unsigned char *image, +/// Functor for 1D Least Squares Gaussian Fit +struct LSGF1DFunctor : Functor { + LSGF1DFunctor(const int nb, const int marg, Eigen::VectorXd X, const unsigned char *image, const int w, const int x0, const int y0) : Functor(3, 2 * nb + 1), nb(nb), @@ -220,8 +224,8 @@ struct LSGFFunctor : Functor { y0(y0) {} /* + Calculate residuals (error = prediction - actual) x = parameters (a, xb, sigma) - fvec = residual (one for every data point, of size = 2*nb+1) */ int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { double a = x(0); @@ -239,6 +243,7 @@ struct LSGFFunctor : Functor { return 0; } + // Calculate Jacobian int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const { double a = x(0); double kb = x(1); @@ -273,7 +278,7 @@ struct LSGFFunctor : Functor { const int x0, y0; }; -/// Functor for 2D Least-squares Gaussian Fit algo +/// 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 @@ -302,7 +307,6 @@ struct LSGF2DFunctor : Functor { for(int j = -nb; j <= nb; j++){ int xi = x0 + i; int yi = y0 + j; - // TODO: other way around, yPred is our modelPred, yActual is our pixel intensity float yPred = FitModel2D(xi, yi, a, xb, yb, sigmaX, sigmaY); float yActual = Get(xi, yi, image, w); fvec(ind) = yPred - yActual; @@ -361,8 +365,6 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW const int nb = 2; - // std::set checkedPoints; - std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); for (const Point &pt : candidatePts) { @@ -379,13 +381,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW vInd++; } - // float a; - // float xb, yb; - // double sigma; - // InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigma); - double a = Get(x, y, image, imageWidth); - double sigma = InitialGuess2(x, y, a, nb, image, imageWidth); + double sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); Eigen::VectorXd betaX(3); betaX << a, x, sigma; @@ -393,23 +390,15 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW Eigen::VectorXd betaY(3); betaY << a, y, sigma; - // LSGF2DFunctor functor(nb, image, imageWidth, x, y); - // // Eigen::NumericalDiff numDiff(functor); - // // Eigen::LevenbergMarquardt, double> lm(numDiff); - // Eigen::LevenbergMarquardt lm(functor); - - LSGFFunctor functorX(nb, 0, X, image, imageWidth, x, y); - // Eigen::NumericalDiff numDiffX(functorX); - Eigen::LevenbergMarquardt lmX(functorX); + LSGF1DFunctor functorX(nb, 0, X, image, imageWidth, x, y); + Eigen::LevenbergMarquardt lmX(functorX); lmX.parameters.maxfev = 2000; lmX.parameters.xtol = 1.0e-10; - LSGFFunctor functorY(nb, 1, Y, image, imageWidth, x, y); - // Eigen::NumericalDiff numDiffY(functorY); - // Eigen::LevenbergMarquardt, double> lmY(numDiffY); - Eigen::LevenbergMarquardt lmY(functorY); + LSGF1DFunctor functorY(nb, 1, Y, image, imageWidth, x, y); + Eigen::LevenbergMarquardt lmY(functorY); lmY.parameters.maxfev = 2000; lmY.parameters.xtol = 1.0e-10; @@ -424,13 +413,11 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW sigma = betaX(2); - // std::cout << "Original: " << x << ", " << y << std::endl; - // std::cout << "final: " << xb << ", " << yb << std::endl; + // TODO: idk how much we care about making radius accurate result.push_back(Star(xb+0.5, yb+0.5, 0)); - // result.push_back(Star(x, y, 0)); } - std::cout << "Number of centroids: " << result.size() << std::endl; + // std::cout << "Number of centroids: " << result.size() << std::endl; return result; } @@ -439,12 +426,10 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW int imageHeight) const { std::vector result; - std::cout << "2D GAUSSIAN FIT" << std::endl; + std::cout << "2D Gaussian Fit" << std::endl; const int nb = 2; - const int np = 2 * nb + 1; - - // std::set checkedPoints; + // const int np = 2 * nb + 1; std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); @@ -454,16 +439,12 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; float a = Get(x, y, image, imageWidth); - double sigma = InitialGuess2(x, y, a, nb, image, imageWidth); - // InitialGuess(x, y, nb, image, imageWidth, &a, &xb, &yb, &sigmaX); + double sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); Eigen::VectorXd beta(5); beta << a, x, y, sigma, sigma; - // LSGF2DFunctor functor(nb, 0, X, image, imageWidth, x, y); LSGF2DFunctor functor(nb, image, imageWidth, x, y); - // Eigen::NumericalDiff numDiff(functor); - // Eigen::LevenbergMarquardt, double> lm(numDiff); Eigen::LevenbergMarquardt lm(functor); lm.parameters.maxfev = 2000; @@ -477,14 +458,12 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW double sigmaX = beta(3); double sigmaY = beta(4); - // std::cout << xRes+0.5 << ", " << yRes+0.5 << std::endl; - // result.push_back(Star(xRes+0.5, yRes+0.5, 0)); // TODO: not sure how much we care about making the radius/brightness accurate result.push_back(Star(xRes+0.5, yRes+0.5, sigmaX, sigmaY, aRes)); } - std::cout << "Number of centroids: " << result.size() << std::endl; + // std::cout << "Number of centroids: " << result.size() << std::endl; return result; } @@ -553,7 +532,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; diff --git a/src/centroiders.hpp b/src/centroiders.hpp index e4b7c21d..8e165d50 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -76,8 +76,6 @@ class IterativeWeightedCenterOfGravityAlgorithm : public CentroidAlgorithm { Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; }; -// TODO: a bunch of functions that really should be private, but tests need to access - /// Get value of pixel at (x, y) in image with width=w int Get(int x, int y, const unsigned char *image, int w); @@ -87,6 +85,15 @@ 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); +/** + * @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: @@ -97,10 +104,6 @@ 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); -// TODO: rename -// Basically initial guess for window, done after Floodfill preprocessing -float InitialGuess2(int x0, int y0, int maxMag, const int nb, const unsigned char *image, int w); - } // namespace lost #endif From 60ed8b2d6423d375385f0d8a1bfbce301c83c1cb Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Fri, 24 Mar 2023 11:26:56 -0700 Subject: [PATCH 19/34] replace double with float --- src/centroiders.cpp | 128 ++++++++++++++++++++++---------------------- 1 file changed, 63 insertions(+), 65 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 3817a07a..326b7ffd 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -140,7 +140,7 @@ 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){ +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; @@ -148,7 +148,7 @@ static float FitModel2D(float x, float y, float a, float xb, float yb, float sig // DO NO DELETE // void InitialGuess(int x0, int y0, const int nb, const unsigned char *image, int w, float *a, -// float *xb, float *yb, double *sigma) { +// 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; @@ -174,15 +174,16 @@ static float FitModel2D(float x, float y, float a, float xb, float yb, float sig // } // } // } -// double fwhm = std::sqrt(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 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){ + 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++; } } @@ -200,8 +201,8 @@ struct Functor { typedef Eigen::Matrix ValueType; typedef Eigen::Matrix JacobianType; - int m_inputs; // Number of parameters in your model - int m_values; // Number of data points + 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) {} @@ -211,10 +212,10 @@ struct Functor { }; /// Functor for 1D Least Squares Gaussian Fit -struct LSGF1DFunctor : Functor { - LSGF1DFunctor(const int nb, const int marg, Eigen::VectorXd X, const unsigned char *image, - const int w, const int x0, const int y0) - : Functor(3, 2 * nb + 1), +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), @@ -227,10 +228,10 @@ struct LSGF1DFunctor : Functor { Calculate residuals (error = prediction - actual) x = parameters (a, xb, sigma) */ - int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - double a = x(0); - double xb = x(1); - double sigma = x(2); + 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) @@ -244,21 +245,21 @@ struct LSGF1DFunctor : Functor { } // Calculate Jacobian - int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const { - double a = x(0); - double kb = x(1); - double sK = x(2); + 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); + 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); + fjac(ind, 1) = a * expK * (k - kb) / (sK * sK); + fjac(ind, 2) = a * expK * (k - kb) * (k - kb) / std::pow(sK, 3); ind++; } @@ -271,7 +272,7 @@ struct LSGF1DFunctor : Functor { // Flag for which marginal to use (0 = X, 1 = Y) const int marg; // Data points (set of x or y coordinates) - Eigen::VectorXd X; + Eigen::VectorXf X; const unsigned char *image; const int w; // image width in pixels // Center coordinates (x0, y0) of this window @@ -279,13 +280,12 @@ struct LSGF1DFunctor : Functor { }; /// Functor for 2D Least Squares Gaussian Fit -struct LSGF2DFunctor : Functor { +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)), + 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), @@ -295,16 +295,16 @@ struct LSGF2DFunctor : Functor { /* x = parameters (a, xb, yb, sigmaX, sigmaY) */ - int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const { - double a = x(0); - double xb = x(1); - double yb = x(2); - double sigmaX = x(3); - double sigmaY = x(4); + 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++){ + 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); @@ -318,12 +318,12 @@ struct LSGF2DFunctor : Functor { return 0; } - int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const{ - double a = x(0); - double xb = x(1); - double yb = x(2); - double sigmaX = x(3); - double sigmaY = x(4); + 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; @@ -372,8 +372,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int y = pt[1]; if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; - Eigen::VectorXd X(2 * nb + 1); - Eigen::VectorXd Y(2 * nb + 1); + Eigen::VectorXf X(2 * nb + 1); + Eigen::VectorXf Y(2 * nb + 1); int vInd = 0; for (int i = -nb; i <= nb; i++) { X(vInd) = x + i; @@ -381,24 +381,23 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW vInd++; } - double a = Get(x, y, image, imageWidth); - double sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); + float a = Get(x, y, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); - Eigen::VectorXd betaX(3); + Eigen::VectorXf betaX(3); betaX << a, x, sigma; - Eigen::VectorXd betaY(3); + Eigen::VectorXf betaY(3); betaY << a, y, sigma; LSGF1DFunctor functorX(nb, 0, X, image, imageWidth, x, y); - Eigen::LevenbergMarquardt lmX(functorX); - + Eigen::LevenbergMarquardt lmX(functorX); lmX.parameters.maxfev = 2000; lmX.parameters.xtol = 1.0e-10; LSGF1DFunctor functorY(nb, 1, Y, image, imageWidth, x, y); - Eigen::LevenbergMarquardt lmY(functorY); + Eigen::LevenbergMarquardt lmY(functorY); lmY.parameters.maxfev = 2000; lmY.parameters.xtol = 1.0e-10; @@ -408,13 +407,13 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW a = betaX(0); - double xb = betaX(1); - double yb = betaY(1); + float xb = betaX(1); + float yb = betaY(1); sigma = betaX(2); // TODO: idk how much we care about making radius accurate - result.push_back(Star(xb+0.5, yb+0.5, 0)); + result.push_back(Star(xb + 0.5, yb + 0.5, 0)); } // std::cout << "Number of centroids: " << result.size() << std::endl; @@ -439,28 +438,27 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW if (x - nb < 0 || x + nb >= imageWidth || y - nb < 0 || y + nb >= imageHeight) continue; float a = Get(x, y, image, imageWidth); - double sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, nb, image, imageWidth); - Eigen::VectorXd beta(5); + Eigen::VectorXf beta(5); beta << a, x, y, sigma, sigma; LSGF2DFunctor functor(nb, image, imageWidth, x, y); - Eigen::LevenbergMarquardt lm(functor); + Eigen::LevenbergMarquardt lm(functor); lm.parameters.maxfev = 2000; lm.parameters.xtol = 1.0e-10; lm.minimize(beta); - double aRes = beta(0); - double xRes = beta(1); - double yRes = beta(2); - double sigmaX = beta(3); - double sigmaY = beta(4); + float aRes = beta(0); + float xRes = beta(1); + float yRes = beta(2); + float sigmaX = beta(3); + float sigmaY = beta(4); // TODO: not sure how much we care about making the radius/brightness accurate - result.push_back(Star(xRes+0.5, yRes+0.5, sigmaX, sigmaY, aRes)); - + result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, aRes)); } // std::cout << "Number of centroids: " << result.size() << std::endl; From 82da27490e0d83153d76fc3b9fc8dba08c6d3904 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sun, 26 Mar 2023 09:47:36 -0700 Subject: [PATCH 20/34] label stars with guess of radius/magnitude --- src/centroiders.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 326b7ffd..b8a22069 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -412,8 +412,7 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW sigma = betaX(2); - // TODO: idk how much we care about making radius accurate - result.push_back(Star(xb + 0.5, yb + 0.5, 0)); + result.push_back(Star(xb + 0.5, yb + 0.5, sigma, sigma, a)); } // std::cout << "Number of centroids: " << result.size() << std::endl; @@ -457,7 +456,6 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW float sigmaX = beta(3); float sigmaY = beta(4); - // TODO: not sure how much we care about making the radius/brightness accurate result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, aRes)); } From 9cd89cfc9b190e715e73a2aad3f5d9e6a5c5c19f Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sun, 26 Mar 2023 10:19:51 -0700 Subject: [PATCH 21/34] Allow user to specify Gaussian fit window size through CLI --- src/centroiders.cpp | 11 ++--------- src/centroiders.hpp | 14 +++++++++++--- src/io.cpp | 4 ++-- src/pipeline-options.hpp | 1 + 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index b8a22069..197f453d 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -361,9 +361,7 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW int imageHeight) const { std::vector result; - std::cout << "1D Gaussian Fit" << std::endl; - - const int nb = 2; + std::cout << "1D Gaussian Fit (" << np << "x" << np << ")" << std::endl; std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); @@ -415,8 +413,6 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW result.push_back(Star(xb + 0.5, yb + 0.5, sigma, sigma, a)); } - // std::cout << "Number of centroids: " << result.size() << std::endl; - return result; } @@ -424,10 +420,7 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW int imageHeight) const { std::vector result; - std::cout << "2D Gaussian Fit" << std::endl; - - const int nb = 2; - // const int np = 2 * nb + 1; + std::cout << "2D Gaussian Fit (" << np << "x" << np << ")" << std::endl; std::vector candidatePts = FloodfillPreproc(image, imageWidth, imageHeight); diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 8e165d50..6802599c 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -30,8 +30,12 @@ class CentroidAlgorithm { */ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ public: - explicit LeastSquaresGaussianFit1D() {}; + explicit LeastSquaresGaussianFit1D(int nb) : nb(nb), np(nb*2+1) { }; Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + +private: + const int nb; + const int np; }; /** @@ -42,9 +46,13 @@ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ * Also more computationally expensive than other methods, scales exponentially with increasing window size */ class LeastSquaresGaussianFit2D : public CentroidAlgorithm { - public: - explicit LeastSquaresGaussianFit2D(){}; +public: + explicit LeastSquaresGaussianFit2D(int nb) : nb(nb), np(nb*2+1) { }; Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + +private: + const int nb; + const int np; }; /// A centroid algorithm for debugging that returns random centroids. diff --git a/src/io.cpp b/src/io.cpp index ac1a6a66..67b23a41 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -734,10 +734,10 @@ Pipeline SetPipeline(const PipelineOptions &values) { std::unique_ptr(new IterativeWeightedCenterOfGravityAlgorithm()); } else if (values.centroidAlgo == "lsgf1d") { result.centroidAlgorithm = - std::unique_ptr(new LeastSquaresGaussianFit1D()); + std::unique_ptr(new LeastSquaresGaussianFit1D(values.centroidFitRadius)); } else if (values.centroidAlgo == "lsgf2d") { result.centroidAlgorithm = - std::unique_ptr(new LeastSquaresGaussianFit2D()); + std::unique_ptr(new LeastSquaresGaussianFit2D(values.centroidFitRadius)); } else if (values.centroidAlgo != "") { std::cout << "Illegal centroid algorithm." << std::endl; exit(1); diff --git a/src/pipeline-options.hpp b/src/pipeline-options.hpp index a90d1c42..0fe77216 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -22,6 +22,7 @@ LOST_CLI_OPTION("fov" , float , fov , 20 , atof(optarg) , // PIPELINE STAGES 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-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) From 05f1fc3fd05cb8717dabba14504261d49cba9524 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sun, 26 Mar 2023 12:48:00 -0700 Subject: [PATCH 22/34] Implement 5x5 Gaussian Grid --- src/centroiders.cpp | 74 +++++++++++++++++++++++++++++++++++++++++++-- src/centroiders.hpp | 17 +++++++++-- src/io.cpp | 3 ++ 3 files changed, 88 insertions(+), 6 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 197f453d..dabc32a3 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -452,7 +452,77 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, aRes)); } - // std::cout << "Number of centroids: " << result.size() << std::endl; + return result; +} + +static void GetGGridCoeff(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; + + 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(Get(x + k, y + j, image, imageWidth)); + } + std::vector a(np); + std::vector b(np); + GetGGridCoeff(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); + GetGGridCoeff(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, Get(x, y, image, imageWidth))); + } return result; } @@ -622,8 +692,6 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f, ((float)(xDiameter)) / 2.0f, ((float)(yDiameter)) / 2.0f, p.checkedIndices.size() - sizeBefore)); - - // std::cout << result.back().position.x << ", " << result.back().position.y << std::endl; } } } diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 6802599c..c04d5d3c 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -79,9 +79,20 @@ 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 diff --git a/src/io.cpp b/src/io.cpp index 67b23a41..82d1e65d 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -738,6 +738,9 @@ Pipeline SetPipeline(const PipelineOptions &values) { } else if (values.centroidAlgo == "lsgf2d") { result.centroidAlgorithm = std::unique_ptr(new LeastSquaresGaussianFit2D(values.centroidFitRadius)); + } 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); From 5ee93ce13e8400e5b213d5f46668ebe71c45aa58 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Mon, 27 Mar 2023 15:33:57 -0700 Subject: [PATCH 23/34] smol --- src/centroiders.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index dabc32a3..ca0d3d93 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -455,7 +455,7 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW return result; } -static void GetGGridCoeff(const std::vector& w, std::vector* const a, std::vector* const b){ +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]; @@ -492,7 +492,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, } std::vector a(np); std::vector b(np); - GetGGridCoeff(w, &a, &b); + 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)); @@ -511,7 +511,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, } std::vector a(np); std::vector b(np); - GetGGridCoeff(w, &a, &b); + 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)); From 6694705093f217562f669ff90dfbb698e254a59a Mon Sep 17 00:00:00 2001 From: Mark Polyakov Date: Sun, 2 Apr 2023 17:21:06 -0700 Subject: [PATCH 24/34] Print total times, and fix 95% logic. --- src/io.cpp | 32 +++++++++++++++++++++++++------- test/scripts/random-crap.sh | 10 ++++++++++ 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 056ad483..136149af 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -1474,7 +1474,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 +1495,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! 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 ' From d3b1b3386b03cd5a47cc47c1f222037f445f5a37 Mon Sep 17 00:00:00 2001 From: Mark Polyakov Date: Tue, 4 Apr 2023 14:31:21 -0700 Subject: [PATCH 25/34] change parameters for image generation to be more actionable We are now talking more in terms of photoelectrons and less in terms of """brightness""" --- .gitignore | 7 +++++ Makefile | 3 --- documentation/pipeline.man | 14 +++++----- src/io.cpp | 31 +++++++++++----------- src/pipeline-options.hpp | 54 +++++++++++++++++++------------------- 5 files changed, 57 insertions(+), 52 deletions(-) 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 0df7c810..f20a4557 100644 --- a/Makefile +++ b/Makefile @@ -57,9 +57,6 @@ release: CXXFLAGS := $(RELEASE_CXXFLAGS) release: LDFLAGS := $(RELEASE_LDFLAGS) release: all -$(BSC): download-bsc.sh - ./download-bsc.sh - $(BIN): $(OBJS) $(CXX) $(LDFLAGS) -o $(BIN) $(OBJS) $(LIBS) diff --git a/documentation/pipeline.man b/documentation/pipeline.man index 969b6285..9a08a384 100644 --- a/documentation/pipeline.man +++ b/documentation/pipeline.man @@ -149,8 +149,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 +170,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/io.cpp b/src/io.cpp index 136149af..c82e46e9 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -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,12 +558,13 @@ 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 peakBrightness = peakBrightnessPerTime * exposureTime; + 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/peakBrightness)*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. @@ -695,7 +694,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 +771,9 @@ PipelineInputList GetGeneratedPipelineInput(const PipelineOptions &values) { &rng, values.generateCentroidsOnly, - values.generateRefBrightness, + values.generateZeroMagPhotons, values.generateSpreadStdDev, - values.generateSensitivity, + values.generateSaturationPhotons, values.generateDarkCurrent, values.generateReadNoiseStdDev, motionBlurDirection, diff --git a/src/pipeline-options.hpp b/src/pipeline-options.hpp index 23e4763c..713da148 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -48,30 +48,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) From 59255ac3c2049f8fcb91405237c60b3684982da8 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Tue, 4 Apr 2023 15:54:47 -0700 Subject: [PATCH 26/34] test weight for ggrid --- src/centroiders.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index ca0d3d93..f936d6fd 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -488,7 +488,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, for (int j = -nb; j <= nb; j++) { std::vector w; for (int k = -nb; k <= nb; k++) { - w.push_back(Get(x + k, y + j, image, imageWidth)); + w.push_back(std::pow(Get(x + k, y + j, image, imageWidth), 1)); } std::vector a(np); std::vector b(np); From c7908f0ef749d1fe73adbca482e48569588f7bc0 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Wed, 5 Apr 2023 08:41:33 -0700 Subject: [PATCH 27/34] Change magnitude of Gaussian methods to be floodfill size --- src/centroiders.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index f936d6fd..80d38895 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -9,7 +9,7 @@ #include #include #include -#include +// #include #include #include // TODO: remove later, this is just lazy to implement hash for unordered_set #include @@ -27,6 +27,7 @@ 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; /** @@ -59,6 +60,8 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im int x0 = x; int y0 = y; + int floodSize = 0; + queue.push_back(pCurr); while (queue.size() != 0) { Point p = queue[0]; @@ -72,6 +75,8 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im int mag = Get(px, py, image, imageWidth); if (mag < cutoff) continue; + floodSize++; + // 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); @@ -93,7 +98,7 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im } // Cool, our flood is done - res.push_back({x0, y0}); + res.push_back({x0, y0, floodSize}); } } @@ -410,7 +415,7 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW sigma = betaX(2); - result.push_back(Star(xb + 0.5, yb + 0.5, sigma, sigma, a)); + result.push_back(Star(xb + 0.5, yb + 0.5, sigma, sigma, pt[2])); } return result; @@ -449,7 +454,7 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW float sigmaX = beta(3); float sigmaY = beta(4); - result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, aRes)); + result.push_back(Star(xRes + 0.5, yRes + 0.5, sigmaX, sigmaY, pt[2])); } return result; @@ -521,7 +526,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, float yRes = y + nom/denom; - result.push_back(Star(xRes + 0.5, yRes + 0.5, nb, nb, Get(x, y, image, imageWidth))); + result.push_back(Star(xRes + 0.5, yRes + 0.5, nb, nb, pt[2])); } return result; From 420637f5a4928149b9750c3b9b96512cca2226a7 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Wed, 5 Apr 2023 09:07:07 -0700 Subject: [PATCH 28/34] Subtract noise from image BEFORE any other preprocessing --- src/centroiders.cpp | 68 ++++++++++++++++++++++++++++++++++----------- src/centroiders.hpp | 3 ++ 2 files changed, 55 insertions(+), 16 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 80d38895..77e67775 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -40,10 +40,25 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im /// 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)); + } +} + std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int imageHeight) { std::vector res; std::set checkedPoints; + SubtractNoise(image, imageWidth, imageHeight); + int cutoff = BasicThreshold(image, imageWidth, imageHeight); for (long i = 0; i < imageHeight * imageWidth; i++) { @@ -368,6 +383,8 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW 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) { @@ -427,6 +444,8 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW 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) { @@ -460,26 +479,40 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW 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]; +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 { +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) { @@ -505,7 +538,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, } } - float xRes = x + nom/denom; + float xRes = x + nom / denom; nom = 0; denom = 0; @@ -524,7 +557,7 @@ Stars GaussianGrid::Go(unsigned char *image, int imageWidth, } } - float yRes = y + nom/denom; + float yRes = y + nom / denom; result.push_back(Star(xRes + 0.5, yRes + 0.5, nb, nb, pt[2])); } @@ -667,6 +700,8 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi 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) { @@ -759,6 +794,7 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im 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 diff --git a/src/centroiders.hpp b/src/centroiders.hpp index c04d5d3c..e4bfe3fa 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -104,6 +104,9 @@ 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 From 648a681acad932be233d2ac5026bb76f7d2d224b Mon Sep 17 00:00:00 2001 From: Mark Polyakov Date: Wed, 5 Apr 2023 09:37:05 -0700 Subject: [PATCH 29/34] filter brightest centroids by number, not just magnitude --- documentation/pipeline.man | 7 +++++++ src/io.cpp | 25 ++++++++++++++++++------- src/io.hpp | 4 ++++ src/pipeline-options.hpp | 19 ++++++++++--------- 4 files changed, 39 insertions(+), 16 deletions(-) diff --git a/documentation/pipeline.man b/documentation/pipeline.man index 9a08a384..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. diff --git a/src/io.cpp b/src/io.cpp index c82e46e9..4a40f2d4 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -835,10 +835,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 @@ -858,7 +854,8 @@ Pipeline SetPipeline(const PipelineOptions &values) { } // 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 != "") { @@ -938,9 +935,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); } } @@ -1588,7 +1599,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 713da148..5c5c66ef 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -21,15 +21,16 @@ 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-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) From e0256f2a0125c67c8f7578f0173f66dfc8cf49a6 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Wed, 5 Apr 2023 09:39:31 -0700 Subject: [PATCH 30/34] edge case temporary fix for cutoff=0 --- src/centroiders.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 77e67775..563dd5c5 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -57,9 +57,13 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im std::vector res; std::set checkedPoints; - SubtractNoise(image, imageWidth, imageHeight); - 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; From 99b65e3c68162b69c3250956d26e83f8399b2473 Mon Sep 17 00:00:00 2001 From: Mark Polyakov Date: Wed, 5 Apr 2023 09:42:26 -0700 Subject: [PATCH 31/34] Keep track of correct centroids (regression?) --- src/io.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/io.cpp b/src/io.cpp index 4a40f2d4..5601daf9 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -1107,11 +1107,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; } From be293be915bec96891d930cd91e2a6d77cac275a Mon Sep 17 00:00:00 2001 From: Mark Polyakov Date: Thu, 6 Apr 2023 10:30:29 -0700 Subject: [PATCH 32/34] Fix motion blur after changing the generate params. --- src/io.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 5601daf9..bcaa3375 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)); } /** @@ -559,19 +559,18 @@ GeneratedPipelineInput::GeneratedPipelineInput(const Catalog &catalog, } // radiant intensity, in photons per time unit per pixel, at the center of the star. float peakBrightnessPerTime = zeroMagPeakPhotonDensity * MagToBrightness(catalogStar.magnitude); - float peakBrightness = peakBrightnessPerTime * exposureTime; 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)*2*M_PI*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: @@ -650,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; } From 3d297443873379e836c2862c65affe6aa8f063b8 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 8 Apr 2023 00:38:19 -0700 Subject: [PATCH 33/34] Add dynamic window resizing option for fitting algos --- src/centroiders.cpp | 53 +++++++++++++++++++++++++++++++--------- src/centroiders.hpp | 6 +++-- src/io.cpp | 4 +-- src/pipeline-options.hpp | 3 ++- 4 files changed, 50 insertions(+), 16 deletions(-) diff --git a/src/centroiders.cpp b/src/centroiders.cpp index 563dd5c5..7c113e41 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -53,6 +53,13 @@ void SubtractNoise(unsigned char *image, int imageWidth, int imageHeight){ } } +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; @@ -75,6 +82,8 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im std::vector pts; std::deque queue; + FloodParams fp{x, x, y, y}; + int maxMag = -1; int x0 = x; int y0 = y; @@ -96,6 +105,11 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im 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); @@ -117,7 +131,7 @@ std::vector FloodfillPreproc(unsigned char *image, int imageWidth, int im } // Cool, our flood is done - res.push_back({x0, y0, floodSize}); + res.push_back({x0, y0, floodSize, fp.xMax - fp.xMin, fp.yMax - fp.yMin}); } } @@ -394,19 +408,27 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW 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; + // TODO: just taking largest diameter right now + int dynamicWinSize = std::max(pt[3], pt[4]); + + int winRadius = (dynamic) ? dynamicWinSize : nb; // /2? + // int winRadius = nb; - Eigen::VectorXf X(2 * nb + 1); - Eigen::VectorXf Y(2 * nb + 1); + 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 = -nb; i <= nb; i++) { + 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, nb, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, winRadius, image, imageWidth); Eigen::VectorXf betaX(3); betaX << a, x, sigma; @@ -414,13 +436,13 @@ std::vector LeastSquaresGaussianFit1D::Go(unsigned char *image, int imageW Eigen::VectorXf betaY(3); betaY << a, y, sigma; - LSGF1DFunctor functorX(nb, 0, X, image, imageWidth, x, y); + 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(nb, 1, Y, image, imageWidth, x, y); + LSGF1DFunctor functorY(winRadius, 1, Y, image, imageWidth, x, y); Eigen::LevenbergMarquardt lmY(functorY); lmY.parameters.maxfev = 2000; @@ -455,15 +477,24 @@ std::vector LeastSquaresGaussianFit2D::Go(unsigned char *image, int imageW 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; + + // 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, nb, image, imageWidth); + float sigma = FitInitialGuessSigma(x, y, a, winRadius, image, imageWidth); Eigen::VectorXf beta(5); beta << a, x, y, sigma, sigma; - LSGF2DFunctor functor(nb, image, imageWidth, x, y); + LSGF2DFunctor functor(winRadius, image, imageWidth, x, y); Eigen::LevenbergMarquardt lm(functor); lm.parameters.maxfev = 2000; diff --git a/src/centroiders.hpp b/src/centroiders.hpp index e4bfe3fa..ec7226a9 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -30,11 +30,12 @@ class CentroidAlgorithm { */ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ public: - explicit LeastSquaresGaussianFit1D(int nb) : nb(nb), np(nb*2+1) { }; + 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; }; @@ -47,11 +48,12 @@ class LeastSquaresGaussianFit1D : public CentroidAlgorithm{ */ class LeastSquaresGaussianFit2D : public CentroidAlgorithm { public: - explicit LeastSquaresGaussianFit2D(int nb) : nb(nb), np(nb*2+1) { }; + 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; }; diff --git a/src/io.cpp b/src/io.cpp index 26385ff1..1d36adb8 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -852,10 +852,10 @@ Pipeline SetPipeline(const PipelineOptions &values) { std::unique_ptr(new IterativeWeightedCenterOfGravityAlgorithm()); } else if (values.centroidAlgo == "lsgf1d") { result.centroidAlgorithm = - std::unique_ptr(new LeastSquaresGaussianFit1D(values.centroidFitRadius)); + std::unique_ptr(new LeastSquaresGaussianFit1D(values.centroidFitRadius, values.centroidFitDynamic)); } else if (values.centroidAlgo == "lsgf2d") { result.centroidAlgorithm = - std::unique_ptr(new LeastSquaresGaussianFit2D(values.centroidFitRadius)); + std::unique_ptr(new LeastSquaresGaussianFit2D(values.centroidFitRadius, values.centroidFitDynamic)); } else if (values.centroidAlgo == "ggrid"){ result.centroidAlgorithm = std::unique_ptr(new GaussianGrid()); diff --git a/src/pipeline-options.hpp b/src/pipeline-options.hpp index c8fc07dd..6c76f682 100644 --- a/src/pipeline-options.hpp +++ b/src/pipeline-options.hpp @@ -22,7 +22,8 @@ LOST_CLI_OPTION("fov" , float , fov , 20 , atof(optarg) , // PIPELINE STAGES 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-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) From 890efb9f7aa539b2c3e2ecbda4e9c02f17bdc988 Mon Sep 17 00:00:00 2001 From: Edward Zhang Date: Sat, 8 Apr 2023 00:41:59 -0700 Subject: [PATCH 34/34] disable tests since I got rid of function --- test/{gaussian-centroid.cpp => gaussian-centroid.cpp.off} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{gaussian-centroid.cpp => gaussian-centroid.cpp.off} (100%) diff --git a/test/gaussian-centroid.cpp b/test/gaussian-centroid.cpp.off similarity index 100% rename from test/gaussian-centroid.cpp rename to test/gaussian-centroid.cpp.off