New Local Thresholding Algorithm

.vscode/c_cpp_properties.json

@@ -0,0 +1,16 @@
+{
+    "configurations": [
+        {
+            "name": "Linux",
+            "includePath": [
+                "${workspaceFolder}/**"
+            ],
+            "defines": [],
+            "compilerPath": "/usr/bin/gcc",
+            "cStandard": "gnu17",
+            "cppStandard": "gnu++14",
+            "intelliSenseMode": "linux-gcc-x64"
+        }
+    ],
+    "version": 4
+}

.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"
     }
 }

Tests/Test 12/Error

@@ -0,0 +1 @@
+Segmentation fault identified at 10 subdivisions

Tests/Testing Format

@@ -0,0 +1,3 @@
+All tests conducted using form:
+./lost pipeline png <IMG Location> 0 30 centroid cog centroid_magnitude_filter 
+4 database kvector.dat starid pyramid .05 100 .0001 attitude dqm done plot_output output.png done

src/centroiders.cpp

@@ -14,7 +14,9 @@
 
 namespace lost {
 
-// DUMMY
+// DUMMYS
+
+long FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions);
 
 std::vector<Star> DummyCentroidAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const {
     std::vector<Star> result;
@@ -92,6 +94,84 @@ int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight) {
     return mean + (std * 5);
 }
 
+int Limit(bool test, int i, int leftover, int div) {
+    if(test) {
+        return i * (div + 1);
+    } else {
+        return leftover * (div + 1) + (i - leftover) * div;
+    }
+}
+// pixel * divisions/ size floor
+// subdivisionsize * size / divisions ceil
+
+// Box is which subdivision
+int Box(int box, unsigned char *image, int imageWidth, int imageHeight, int subdivisions, int horizontalLeftover, int horizontalDiv, int verticalLeftover, int verticalDiv) {
+    int row = box / subdivisions; // Finds which row on the subdivisions we're in
+    int col = box % subdivisions; // Finds the column on the subdivisions we're in
+    double average = 0;
+    double squareSum = 0;
+    long count = 0;
+    //std::cout << "Box " << box << "\n";
+    //std::cout << Limit(row < horizontalLeftover, row, horizontalLeftover, horizontalDiv) << " " << Limit(row < horizontalLeftover, row + 1, horizontalLeftover, horizontalDiv) << "\n";
+    //std::cout << Limit(col < verticalLeftover, col, verticalLeftover, verticalDiv) << " " << Limit(col < verticalLeftover, col + 1, verticalLeftover, verticalDiv) << "\n";
+    //std::cout << "\n";
+    int secondCount = 0;
+    int number;
+    for(int i = Limit(row < horizontalLeftover, row, horizontalLeftover, horizontalDiv); i < Limit(row < horizontalLeftover, row + 1, horizontalLeftover, horizontalDiv); i++) {
+        for(int j = Limit(col < verticalLeftover, col, verticalLeftover, verticalDiv); j < Limit(col < verticalLeftover, col + 1, verticalLeftover, verticalDiv); j++) {
+            average += image[i * imageWidth + j];
+            squareSum += image[i * imageWidth + j] * image[i * imageWidth + j];
+            count++;
+
+            // Checks for Error
+            if(FindSubdivision(i * imageWidth + j, imageWidth, imageHeight, subdivisions) != box) {
+                secondCount++;
+                number = FindSubdivision(i * imageWidth + j, imageWidth, imageHeight, subdivisions);
+            }
+        }
+    }
+
+    // Checks for Error
+    if(secondCount != 0) {
+        std::cout << "ERROR, Index-Subdivision Mismatching - " << "Box: " << box << " Number: " << number << " Count: "<< count << " Wrong: " <<  secondCount << "\n";
+    }
+    // double varianceNaive = 0;
+    // for(int i = Limit(row < horizontalLeftover, row, horizontalLeftover, horizontalDiv); i < Limit(row < horizontalLeftover, row + 1, horizontalLeftover, horizontalDiv); i++) {
+    //    for(int j = Limit(col < verticalLeftover, col, verticalLeftover, verticalDiv); j < Limit(col < verticalLeftover, col + 1, verticalLeftover, verticalDiv); j++) {
+    //        varianceNaive += std::pow(image[i*imageWidth+j]-(average/count),2);
+    //    }
+    // }
+    // std::cout << std::sqrt(varianceNaive/(count-1)) << " ";
+    average /= count;
+    // std::cout << average + (5 * std::sqrt((squareSum - count * average * average) / (count - 1))) << "\n";
+
+    // std::cout << squareSum << " "; 
+    // std::cout << average << " ";
+    // std::cout << sqrt((squareSum - count * average * average)/(count-1));
+    // std::cout << "\n";
+    return average + (5 * std::sqrt((squareSum - count * average * average) / (count - 1)));
+}
+
+// Local Basic Thresholding with divisions as parameter, returns a vector containing thresholds 
+// corresponding to different divisions.
+// Divisions are horizontally based
+// Uses the same statistical algorithm as BasicThresholding
+std::vector<int> LocalBasicThresholding(unsigned char *image, int imageWidth, int imageHeight, int subdivisions) {
+    // run Basic Threshold on all elements in certain subdivisions
+    int horizontalDiv = imageHeight / subdivisions; // Minimum number of lines for each division
+    int verticalDiv = imageWidth / subdivisions;
+    int horizontalLeftover = imageHeight % subdivisions; // Determines the first few lines that have 1 more line
+    int verticalLeftover = imageWidth % subdivisions;
+    std::vector<int> standardDeviations;
+    // std::cout << subdivisions << "\n";
+    // Sets threshold for the first few subdivisions that have 1 more line than the rest
+    for(int i = 0; i < subdivisions * subdivisions; i++) {
+        standardDeviations.push_back(Box(i, image, imageWidth, imageHeight, subdivisions, horizontalLeftover, horizontalDiv, verticalLeftover, verticalDiv));
+    }
+    // Return values of previous method as a vector
+    return standardDeviations;
+}
+
 // basic thresholding, but do it faster (trade off of some accuracy?)
 int BasicThresholdOnePass(unsigned char *image, int imageWidth, int imageHeight) {
     unsigned long totalMag = 0;
@@ -105,7 +185,7 @@ 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 (int) (mean + (std * 5));
 }
 
 struct CentroidParams {
@@ -117,14 +197,33 @@ struct CentroidParams {
     int yMin;
     int yMax;
     int cutoff;
+    std::vector<int> localCutoff;
     bool isValid;
     std::unordered_set<int> checkedIndices;
 };
 
+// int iToDivision
+
+// int Row(i, imageWidth, subdivisions, imageHeight / subdivisions, imageHeight % subdivisions)
+// int Column(i, imageHeight, subdivisions, imageWidth / subdivisions, imageWidth % subdivisions)
+long RowOrColumn(long i, int size, int subdivisions, int div, int leftover) {
+    if(i < (div + 1) * leftover) {
+        return i / (div + 1);
+    } else {
+        return leftover + (i - (div + 1) * leftover) / (div);
+    }
+}
+
+// For a given i and picture dimensions, determines which subdivision i is in (Zero Based)
+long FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions) {
+    return RowOrColumn(i / imageWidth, imageWidth, subdivisions, imageHeight / subdivisions, imageHeight % subdivisions) * subdivisions + 
+    RowOrColumn(i % imageWidth, imageHeight, subdivisions, imageWidth / subdivisions, imageWidth % subdivisions);
+}
+
 //recursive helper here
-void CogHelper(CentroidParams &p, long i, unsigned char *image, int imageWidth, int imageHeight) {
+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) {
             p.isValid = false;
@@ -144,25 +243,35 @@ 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);
     }
 }
 
+//CenterOfGravityAlgorithm, except the threshold changes depending on the vertical height in the image
+//Subdivisions refers to how many horizontal sections with different thresholds are present
 std::vector<Star> 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<Star> result;
-
-    p.cutoff = BasicThreshold(image, imageWidth, imageHeight);
+    p.localCutoff = LocalBasicThresholding(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 +287,7 @@ std::vector<Star> 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;
 

src/centroiders.hpp

@@ -24,8 +24,10 @@ 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 {

src/io.cpp

@@ -268,7 +268,8 @@ CentroidAlgorithm *DummyCentroidAlgorithmPrompt() {
 }
 
 CentroidAlgorithm *CoGCentroidAlgorithmPrompt() {
-    return new CenterOfGravityAlgorithm();
+    int subdivisions = Prompt<int>("How many subdivisions to use");
+    return new CenterOfGravityAlgorithm(subdivisions);
 }
 
 CentroidAlgorithm *IWCoGCentroidAlgorithmPrompt() {

test/centroid.cpp

@@ -0,0 +1,62 @@
+#include <catch.hpp>
+
+#include <math.h>
+#include <iostream>
+
+namespace lost {
+    int RowOrColumn(long i, int size, int subdivisions, int div, int leftover);
+    int FindSubdivision(long i, int imageWidth, int imageHeight, int subdivisions);
+    int Limit(bool test, int i, int leftover, int div);
+}
+
+
+TEST_CASE("Testing Start and End Rows") {
+    int row = 0;
+    int horizontalLeftover = 5;
+    int horizontalDiv = 10;
+    CHECK(lost::Limit(row < horizontalLeftover, row, horizontalLeftover, horizontalDiv) == 0);
+    CHECK(lost::Limit(row < horizontalLeftover, row + 1, horizontalLeftover, horizontalDiv) == 11);
+}
+
+TEST_CASE("Testing Start and End Columns") {
+    int col = 0;
+    int verticalLeftover = 5;
+    int verticalDiv = 10;
+    CHECK(lost::Limit(col < verticalLeftover, col, verticalLeftover, verticalDiv) == 0);
+    CHECK(lost::Limit(col < verticalLeftover, col + 1, verticalLeftover, verticalDiv) == 11);
+}
+
+// int Row(i, imageWidth, subdivisions, imageHeight / subdivisions, imageHeight % subdivisions)
+// int Column(i, imageHeight, subdivisions, imageWidth / subdivisions, imageWidth % subdivisions)
+TEST_CASE("Correct Row: 4th Row") {
+    CHECK(lost::RowOrColumn(33795, 1024, 100, 1024 / 100, 1024 % 100) == 3);
+}
+
+TEST_CASE("Correct Row: 27th Row") {
+    CHECK(lost::RowOrColumn(24 * 1024 * 11 + 10 * 2 * 1024 + 20, 1024, 100, 1024 / 100, 1024 % 100) == 26);
+}
+
+TEST_CASE("Everything is in the Box") {
+    int row = 0;
+    int horizontalLeftover = 5;
+    int horizontalDiv = 10;
+    int col = 0;
+    int verticalLeftover = 5;
+    int verticalDiv = 10;
+
+    for(int i = Limit(row < horizontalLeftover, row, horizontalLeftover, horizontalDiv); i < Limit(row < horizontalLeftover, row + 1, horizontalLeftover, horizontalDiv); i++) {
+        for(int j = Limit(col < verticalLeftover, col, verticalLeftover, verticalDiv); j < Limit(col < verticalLeftover, col + 1, verticalLeftover, verticalDiv); j++) {
+            if(lost::FindSubdivision(i * ))
+        }
+    }
+}
+
+TEST_CASE("Correct Row: Last Row 1") {
+    CHECK(lost::RowOrColumn(1024 * 1024 - 1, 1024, 100, 1024 / 100, 1024 % 100) == 99);
+}
+
+TEST_CASE("Correct Box: Last Subdivision") {
+    CHECK(lost::FindSubdivision(1024 * 1024 - 1, 1024, 1024, 1024) == 1024 * 1024 - 1);
+}
+
+