diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 00000000..182a0258 --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,17 @@ +{ + "configurations": [ + { + "name": "Linux", + "includePath": [ + "${workspaceFolder}/**" + ], + "defines": [], + "compilerPath": "/usr/bin/gcc", + "cStandard": "gnu17", + "cppStandard": "gnu++14", + "intelliSenseMode": "linux-gcc-x64", + "configurationProvider": "ms-vscode.makefile-tools" + } + ], + "version": 4 +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index 9b5fdd53..c5c0ad1c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -47,6 +47,21 @@ "streambuf": "cpp", "cinttypes": "cpp", "typeinfo": "cpp", - "files.eol": "\n" + "files.eol": "\n", + "regex": "cpp", + "bitset": "cpp", + "chrono": "cpp", + "complex": "cpp", + "condition_variable": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "set": "cpp", + "ratio": "cpp", + "future": "cpp", + "iomanip": "cpp", + "mutex": "cpp", + "shared_mutex": "cpp", + "thread": "cpp", + "variant": "cpp" } } \ No newline at end of file diff --git a/src/centroiders.cpp b/src/centroiders.cpp index aec56f50..b4d76f47 100644 --- a/src/centroiders.cpp +++ b/src/centroiders.cpp @@ -11,11 +11,11 @@ #include #include +#include namespace lost { -// DUMMY - +// DUMMYS std::vector DummyCentroidAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { std::vector result; @@ -26,7 +26,7 @@ std::vector DummyCentroidAlgorithm::Go(unsigned char *image, int imageWidt return result; } -// a poorly designed thresholding algorithm +// UNUSED: 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; @@ -36,7 +36,7 @@ int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) { return (((totalMag/(imageHeight * imageWidth)) + 1) * 15) / 10; } -// a more sophisticated thresholding algorithm, not tailored to star images +// UNUSED: 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; @@ -76,7 +76,7 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) { return level; } -// a simple, but well tested thresholding algorithm that works well with star images +// UNUSED: 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; @@ -92,7 +92,9 @@ int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight) { return mean + (std * 5); } -// basic thresholding, but do it faster (trade off of some accuracy?) +// UNUSED: Accepts an array of brightnesses and image dimensions, and returns the +// threshold brightness for the entire image, which is defined as 5 standard deviations +// above the mean. int BasicThresholdOnePass(unsigned char *image, int imageWidth, int imageHeight) { unsigned long totalMag = 0; float std = 0; @@ -105,9 +107,98 @@ int BasicThresholdOnePass(unsigned char *image, int imageWidth, int imageHeight) float mean = totalMag / totalPixels; float variance = (sq_totalMag / totalPixels) - (mean * mean); std = std::sqrt(variance); - return mean + (std * 5); + return std::round(mean + (std * 5)); +} + +// *****START OF COG ALGORITHM***** + +// This algorithm takes an image and finds all "centroids" or potential stars in it. It divides +// the image into a square grid, calculates the threshold brightness of each subsection, or "box", +// and uses that information to find the centroids. A threshold is defined as 5 standard deviations +// above the mean brightness of pixels. + +// Commonly used terms: + // box: A subsection in the subdivision scheme + // subdivision scheme: The square grid that contains the boxes the image is divided into. Each + // box can be assigned an index in this scheme, which follows row-major order. + // subdivisions: The number of slices the image is sliced into on each side. The square of + // this quantity indicates the number of boxes. + // leftover: When doing this division, we often cannot give each box the same number of + // rows/columns to each box. This quantity tell us how many rows/columns in the subdivision + // scheme will have 1 extra row/column, and these always start with the first rows/columns. + // div: The minimum rows/columns each box will have. + // For leftover and div, "horizontal" refers to rows and "vertical" refers to columns + +// By default, this algorithm ensures that regardless of the subdivisions entered, the number of +// pixels in each subdivision is imageWidth*imageHeight <= size <= 100 to ensure that enough +// pixles will be sampled in each subdivision + +// Accepts an index in the subdivision scheme, and a horizontal/vertical leftover and div +// This method uses these quantities to find the first row/column index in the actual image +// that "i" references. +int StartOfSubdivision(int i, int leftover, int div) { + if(i < leftover) { + return i * (div + 1); + } else { + return leftover * (div + 1) + (i - leftover) * div; + } +} + +// Refer to definition for details +long FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions); + + +// Accepts a subdivsion, or "box" number, the image's brightness array and associated dimensions, +// and the subdivisions. +// Returns the threshold for the corresponding box. +int LocalBasicThreshold(int box, unsigned char *image, int imageWidth, int imageHeight, + int subdivisions) { + + int horizontalDiv = imageHeight / subdivisions; + int verticalDiv = imageWidth / subdivisions; + int horizontalLeftover = imageHeight % subdivisions; + int verticalLeftover = imageWidth % subdivisions; + int row = box / subdivisions; // Finds the current row in the subdivision scheme + int col = box % subdivisions; // Finds the current column in the subdivision scheme + double average = 0; + double squareSum = 0; + long count = 0; + + // Runs through "box" in row-major order + for(long i = StartOfSubdivision(row, horizontalLeftover, horizontalDiv); + i < StartOfSubdivision(row + 1, horizontalLeftover, horizontalDiv); i++) { + for(long j = StartOfSubdivision(col, verticalLeftover, verticalDiv); + j < StartOfSubdivision(col + 1, verticalLeftover, verticalDiv); j++) { + average += image[i * imageWidth + j]; + squareSum += image[i * imageWidth + j] * image[i * imageWidth + j]; + count++; + + // Checks for Error + assert(FindSubdivision(i*imageWidth+j, imageWidth, imageHeight, subdivisions) == box); + } + } + + average /= count; + return std::round(average + (5 * std::sqrt((squareSum - count * average * average) / (count - 1)))); +} + +// Accepts the image's brightness array and dimensions, and the subdivisions, and returns +// a vector with the threshold for each "box" in row-major order +std::vector LocalThresholding(unsigned char *image, int imageWidth, int imageHeight, + int subdivisions) { + + std::vector thresholds; + for(int i = 0; i < subdivisions * subdivisions; i++) { + thresholds.push_back(LocalBasicThreshold(i, image, imageWidth, imageHeight, + subdivisions)); + } + + return thresholds; } +// Used in function CenterOfGravityAlgorithm:Go, and is useful for keeping track of +// star dimensions and characteristics, as well as the database of thresholds among +// other quantities relating to the image struct CentroidParams { float yCoordMagSum; float xCoordMagSum; @@ -116,17 +207,47 @@ struct CentroidParams { int xMax; int yMin; int yMax; - int cutoff; + std::vector localCutoff; bool isValid; std::unordered_set checkedIndices; }; -//recursive helper here -void CogHelper(CentroidParams &p, long i, unsigned char *image, int imageWidth, int imageHeight) { +// Accepts the row/column of an index on the image, the imageWidth/imageHeight, the subdivisions +// and the horizontal/vertical div and leftover +// Returns the index of the row/column in the subdivision scheme for the corresponding "i". +// NOTE: This function is the inverse of the StartOfSubdivision function in a way that +// RowOrColumn(StartOfSubdivision(x)) = x for any given row/column index "x" in the subdivision +// scheme, but StartOfSubdivision(RowOrColumn(x)) only equals x if x = StartOfSubdivision(y), +// for any given row/column index "y" in the subdivision scheme +long RowOrColumn(long i, int leftover, int div) { + if(i < (div + 1) * leftover) { + return i / (div + 1); + } else { + return leftover + (i - (div + 1) * leftover) / (div); + } +} + +// Accepts an index in the image, its dimensions and the subdivisions, and returns the +// corresponding "box" number that the index "i" is in +long FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions) { + long row = RowOrColumn(i / imageWidth, imageHeight % subdivisions, imageHeight / subdivisions); + long col = RowOrColumn(i % imageWidth, imageWidth % subdivisions, imageWidth / subdivisions); + return row * subdivisions + col; +} + +// Accepts a CentroidParams struct, the current index, image's array of brightnesses and dimensions, +// and the subdivisions. This is a recursive helper method that is used with +// CenterOfGravityAlgorithm:Go and finds the dimensions of an individual star in the context of the +// image +void CogHelper(CentroidParams &p, long i, unsigned char *image, int imageWidth, int imageHeight, + int subdivisions) { - if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { + if (i >= 0 && i < imageWidth * imageHeight && + image[i] >= p.localCutoff.at(FindSubdivision(i, imageWidth, imageHeight, subdivisions)) + && p.checkedIndices.count(i) == 0) { //check if pixel is on the edge of the image, if it is, we dont want to centroid this star - if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || i / imageWidth == imageHeight - 1) { + if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || + i / imageWidth == imageHeight - 1) { p.isValid = false; } p.checkedIndices.insert(i); @@ -144,25 +265,36 @@ void CogHelper(CentroidParams &p, long i, unsigned char *image, int imageWidth, p.xCoordMagSum += ((i % imageWidth)) * image[i]; p.yCoordMagSum += ((i / imageWidth)) * image[i]; if(i % imageWidth != imageWidth - 1) { - CogHelper(p, i + 1, image, imageWidth, imageHeight); + CogHelper(p, i + 1, image, imageWidth, imageHeight, subdivisions); } if (i % imageWidth != 0) { - CogHelper(p, i - 1, image, imageWidth, imageHeight); + CogHelper(p, i - 1, image, imageWidth, imageHeight, subdivisions); } - CogHelper(p, i + imageWidth, image, imageWidth, imageHeight); - CogHelper(p, i - imageWidth, image, imageWidth, imageHeight); + CogHelper(p, i + imageWidth, image, imageWidth, imageHeight, subdivisions); + CogHelper(p, i - imageWidth, image, imageWidth, imageHeight, subdivisions); } } +// Accepts an array of the image's brightnesses, and the image's dimensions, and finds all +// stars in the image and returns the stars as an array std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { CentroidParams p; - + // Program will use divisions to represent the subdivisions + int divisions = subdivisions; + int min = 0; + if(imageWidth > imageHeight) { + min = imageWidth; + } else { + min = imageHeight; + } + if(min / subdivisions < 10) { + divisions = min / 10; + } std::vector result; - - p.cutoff = BasicThreshold(image, imageWidth, imageHeight); + p.localCutoff = LocalThresholding(image, imageWidth, imageHeight, divisions); for (long i = 0; i < imageHeight * imageWidth; i++) { - if (image[i] >= p.cutoff && p.checkedIndices.count(i) == 0) { - + if (image[i] >= p.localCutoff.at(FindSubdivision(i, imageWidth, imageHeight, divisions)) + && p.checkedIndices.count(i) == 0) { //iterate over pixels that are part of the star int xDiameter = 0; //radius of current star int yDiameter = 0; @@ -178,7 +310,7 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi int sizeBefore = p.checkedIndices.size(); - CogHelper(p, i, image, imageWidth, imageHeight); + CogHelper(p, i, image, imageWidth, imageHeight, divisions); xDiameter = (p.xMax - p.xMin) + 1; yDiameter = (p.yMax - p.yMin) + 1; @@ -187,13 +319,16 @@ std::vector CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWi float yCoord = (p.yCoordMagSum / (p.magSum * 1.0)); if (p.isValid) { - result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f, ((float)(xDiameter))/2.0f, ((float)(yDiameter))/2.0f, p.checkedIndices.size() - sizeBefore)); + result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f, ((float)(xDiameter))/2.0f, + ((float)(yDiameter))/2.0f, p.checkedIndices.size() - sizeBefore)); } } } return result; } +// *****END OF COG ALGORITHM***** + //Determines how accurate and how much iteration is done by the IWCoG algorithm, //smaller means more accurate and more iterations. float iWCoGMinChange = 0.0002; @@ -203,17 +338,19 @@ struct IWCoGParams { int xMax; int yMin; int yMax; - int cutoff; + std::vector localCutoff; 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) { +void IWCoGHelper(IWCoGParams &p, long i, unsigned char *image, int imageWidth, int imageHeight, int subdivisions, std::vector &starIndices) { + if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p.localCutoff.at(FindSubdivision(i, imageWidth, imageHeight, subdivisions)) + && p.checkedIndices.count(i) == 0) { //check if pixel is on the edge of the image, if it is, we dont want to centroid this star - if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 || i / imageWidth == imageHeight - 1) { + if (i % imageWidth == 0 || i % imageWidth == imageWidth - 1 || i / imageWidth == 0 + || i / imageWidth == imageHeight - 1) { p.isValid = false; } p.checkedIndices.insert(i); @@ -233,23 +370,34 @@ void IWCoGHelper(IWCoGParams &p, long i, unsigned char *image, int imageWidth, i p.yMin = i / imageWidth; } if(i % imageWidth != imageWidth - 1) { - IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i + 1, image, imageWidth, imageHeight, subdivisions, starIndices); } if (i % imageWidth != 0) { - IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i - 1, image, imageWidth, imageHeight, subdivisions, starIndices); } - IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, starIndices); - IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i + imageWidth, image, imageWidth, imageHeight, subdivisions, starIndices); + IWCoGHelper(p, i - imageWidth, image, imageWidth, imageHeight, subdivisions, starIndices); } } Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const { IWCoGParams p; + int divisions = subdivisions; + int min = 0; + if(imageWidth > imageHeight) { + min = imageWidth; + } else { + min = imageHeight; + } + if(min / subdivisions < 10) { + divisions = min / 10; + } std::vector result; - p.cutoff = BasicThreshold(image, imageWidth, imageHeight); + p.localCutoff = LocalThresholding(image, imageWidth, imageHeight, divisions); 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) { + if (image[i] >= p.localCutoff.at(FindSubdivision(i, imageWidth, imageHeight, subdivisions)) + && p.checkedIndices.count(i) == 0) { // TODO: store longs --Mark std::vector starIndices; //indices of the current star p.maxIntensity = 0; @@ -269,7 +417,7 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im p.isValid = true; - IWCoGHelper(p, i, image, imageWidth, imageHeight, starIndices); + IWCoGHelper(p, i, image, imageWidth, imageHeight, divisions, starIndices); xDiameter = (p.xMax - p.xMin) + 1; yDiameter = (p.yMax - p.yMin) + 1; diff --git a/src/centroiders.hpp b/src/centroiders.hpp index 5b6c2477..441b3e4c 100644 --- a/src/centroiders.hpp +++ b/src/centroiders.hpp @@ -24,14 +24,18 @@ class DummyCentroidAlgorithm: public CentroidAlgorithm { class CenterOfGravityAlgorithm : public CentroidAlgorithm { public: - CenterOfGravityAlgorithm() { }; + CenterOfGravityAlgorithm(int subdivisions) : subdivisions(subdivisions) { }; Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + private: + int subdivisions; }; class IterativeWeightedCenterOfGravityAlgorithm : public CentroidAlgorithm { public: - IterativeWeightedCenterOfGravityAlgorithm() { }; + IterativeWeightedCenterOfGravityAlgorithm(int subdivisions) : subdivisions(subdivisions) { }; Stars Go(unsigned char *image, int imageWidth, int imageHeight) const override; + private: + int subdivisions; }; } diff --git a/src/io.cpp b/src/io.cpp index 7920514d..442f48a5 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -268,11 +268,13 @@ CentroidAlgorithm *DummyCentroidAlgorithmPrompt() { } CentroidAlgorithm *CoGCentroidAlgorithmPrompt() { - return new CenterOfGravityAlgorithm(); + int subdivisions = Prompt("How many subdivisions to use"); + return new CenterOfGravityAlgorithm(subdivisions); } CentroidAlgorithm *IWCoGCentroidAlgorithmPrompt() { - return new IterativeWeightedCenterOfGravityAlgorithm(); + int subdivisions = Prompt("How many subdivisions to use"); + return new IterativeWeightedCenterOfGravityAlgorithm(subdivisions); } StarIdAlgorithm *DummyStarIdAlgorithmPrompt() { diff --git a/test/centroid.cpp b/test/centroid.cpp new file mode 100644 index 00000000..277a0b7c --- /dev/null +++ b/test/centroid.cpp @@ -0,0 +1,56 @@ +#include + +#include +#include + +namespace lost { + int RowOrColumn(long i, int leftover, int div); + int FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions); + int StartOfSubdivision(int i, int leftover, int div); +} + + +TEST_CASE("Testing Start and End Rows") { + int row = 0; + int horizontalLeftover = 5; + int horizontalDiv = 10; + CHECK(lost::StartOfSubdivision(row, horizontalLeftover, horizontalDiv) == 0); + CHECK(lost::StartOfSubdivision(row + 1, horizontalLeftover, horizontalDiv) == 11); +} + +TEST_CASE("Testing Start and End Columns") { + int col = 0; + int verticalLeftover = 5; + int verticalDiv = 10; + CHECK(lost::StartOfSubdivision(col, verticalLeftover, verticalDiv) == 0); + CHECK(lost::StartOfSubdivision(col + 1, verticalLeftover, verticalDiv) == 11); +} + +TEST_CASE("Correct Row: 4th Row") { + CHECK(lost::RowOrColumn(33795 / 1024, 1024 % 100, 1024 / 100) == 3); +} + +TEST_CASE("Correct Row: 27th Row") { + CHECK(lost::RowOrColumn((24 * 1024 * 11 + 10 * 2 * 1024 + 20) / 1024, 1024 % 100, 1024 / 100) == 26); +} + +TEST_CASE("Correct Row: Last Row 1") { + CHECK(lost::RowOrColumn((1024 * 1024 - 1) / 1024, 1024 % 100, 1024 / 100) == 99); +} + +TEST_CASE("Correct Box: Last Subdivision") { + CHECK(lost::FindSubdivision(1024 * 1024 - 1, 1024, 1024, 1024) == 1024 * 1024 - 1); +} + +TEST_CASE("Testing Inverse Conjecture") { + int size = 1024; + int subdivisions = 10; + int leftover = size % subdivisions; + int div = size / subdivisions; + for(int i = 0; i < subdivisions; i++) { + int x = lost::StartOfSubdivision(i, leftover, div); + CHECK(lost::RowOrColumn(x, leftover, div) == i); + CHECK(lost::StartOfSubdivision(lost::RowOrColumn(x, leftover, div), leftover, div) == x); + } +} +