Skip to content

Add files via upload #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
simonpi429 wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Add files via upload #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
271cd7f
Add files via upload
simonpi429 Nov 23, 2025
19d2999
centroiders.cpp is now the faster (?) version. Changed decimal class …
simonpi429 Nov 23, 2025
0fb3486
cleaned up calculations
simonpi429 Nov 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 89 additions & 80 deletions src/centroiders.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,17 +1,13 @@
#include "centroiders.hpp"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <cmath>
#include <vector>
#include <iostream>
#include <unordered_map>
#include <unordered_set>

#include "decimal.hpp"

namespace lost {

// DUMMY
Expand All @@ -21,7 +17,7 @@ std::vector<Star> DummyCentroidAlgorithm::Go(unsigned char *, int imageWidth, in

unsigned int randomSeed = 123456;
for (int i = 0; i < numStars; i++) {
result.push_back(Star(rand_r(&randomSeed) % imageWidth, rand_r(&randomSeed) % imageHeight, DECIMAL(10.0)));
result.push_back(Star(rand_r(&randomSeed) % imageWidth, rand_r(&randomSeed) % imageHeight, 10.0));
}

return result;
Expand All @@ -41,11 +37,11 @@ int BadThreshold(unsigned char *image, int imageWidth, int imageHeight) {
int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) {
// code here, duh
long total = imageWidth * imageHeight;
//decimal top = 255;
decimal sumB = 0;
decimal sum1 = 0;
decimal wB = 0;
decimal maximum = 0;
//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];
Expand All @@ -59,12 +55,12 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) {
sum1 += i * histogram[i];
}
for (int i = 0; i < 256; i ++) {
decimal wF = total - wB;
float wF = total - wB;
//std::cout << "wF\n" << wB << "\n";
//std::cout << "wB\n" << wF << "\n";
if (wB > 0 && wF > 0) {
decimal mF = (sum1 - sumB) / wF;
decimal val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF);
float mF = (sum1 - sumB) / wF;
float val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF);
//std::cout << val << "\n";
if (val >= maximum) {
level = i;
Expand All @@ -80,74 +76,78 @@ int OtsusThreshold(unsigned char *image, int imageWidth, int imageHeight) {
// a simple, but well tested thresholding algorithm that works well with star images
int BasicThreshold(unsigned char *image, int imageWidth, int imageHeight) {
unsigned long totalMag = 0;
decimal std = 0;
float std = 0;
long totalPixels = imageHeight * imageWidth;
for (long i = 0; i < totalPixels; i++) {
totalMag += image[i];
}
decimal mean = totalMag / totalPixels;
float mean = totalMag / totalPixels;
for (long i = 0; i < totalPixels; i++) {
std += DECIMAL_POW(image[i] - mean, 2);
std += std::pow(image[i] - mean, 2);
}
std = DECIMAL_SQRT(std / totalPixels);
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;
decimal std = 0;
decimal sq_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];
}
decimal mean = totalMag / totalPixels;
decimal variance = (sq_totalMag / totalPixels) - (mean * mean);
std = DECIMAL_SQRT(variance);
float mean = totalMag / totalPixels;
float variance = (sq_totalMag / totalPixels) - (mean * mean);
std = std::sqrt(variance);
return mean + (std * 5);
}

struct CentroidParams {
decimal yCoordMagSum;
decimal xCoordMagSum;
float yCoordMagSum;
float xCoordMagSum;
long magSum;
int xMin;
int xMax;
int yMin;
int yMax;
int cutoff;
bool isValid;
std::unordered_set<int> checkedIndices;
std::vector<bool> checkedIndices;
int pixelCount;
};

//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) {
if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && !p->checkedIndices[i]) {
int x = i % imageWidth;
int y = i / imageWidth;
//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 (x == 0 || x == imageWidth - 1 || y == 0 || y == 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;
p->checkedIndices[i] = true;
p.pixelCount++;
if (x > p->xMax) {
p->xMax = x;
} else if (x < p->xMin) {
p->xMin = x;
}
if (i / imageWidth > p->yMax) {
p->yMax = i / imageWidth;
} else if (i / imageWidth < p->yMin) {
p->yMin = i / imageWidth;
if (y > p->yMax) {
p->yMax = y;
} else if (y < p->yMin) {
p->yMin = y;
}
p->magSum += image[i];
p->xCoordMagSum += ((i % imageWidth)) * image[i];
p->yCoordMagSum += ((i / imageWidth)) * image[i];
if (i % imageWidth != imageWidth - 1) {
p->xCoordMagSum += ((x)) * image[i];
p->yCoordMagSum += ((y)) * image[i];
if (x != imageWidth - 1) {
CogHelper(p, i + 1, image, imageWidth, imageHeight);
}
if (i % imageWidth != 0) {
if (x != 0) {
CogHelper(p, i - 1, image, imageWidth, imageHeight);
}
CogHelper(p, i + imageWidth, image, imageWidth, imageHeight);
Expand All @@ -157,47 +157,54 @@ void CogHelper(CentroidParams *p, long i, unsigned char *image, int imageWidth,

std::vector<Star> CenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const {
CentroidParams p;
p.checkedIndices.assign(imageWidth * imageHeight, false);

std::vector<Star> 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) {
if (image[i] >= p.cutoff && !p.checkedIndices[i]) {

int x = i % imageWidth;
int y = i / imageWidth;
//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.xMax = x;
p.xMin = x;
p.yMax = y;
p.yMin = y;
p.isValid = true;

int sizeBefore = p.checkedIndices.size();
p.pixelCount = 0;

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
decimal xCoord = (p.xCoordMagSum / (p.magSum * DECIMAL(1.0)));
decimal yCoord = (p.yCoordMagSum / (p.magSum * DECIMAL(1.0)));
float xCoord = p.xCoordMagSum / p.magSum;
float yCoord = p.yCoordMagSum / p.magSum;

if (p.isValid) {
result.push_back(Star(xCoord + DECIMAL(0.5), yCoord + DECIMAL(0.5), (xDiameter)/DECIMAL(2.0), (yDiameter)/DECIMAL(2.0), p.checkedIndices.size() - sizeBefore));
result.push_back(Star(xCoord + 0.5f, yCoord + 0.5f,
((float)(xDiameter))/2.0f,
((float)(yDiameter))/2.0f,
p.pixelCount));
}
}
}
std::cout << "centroiding: " << result.size() << std::endl;
return result;
}

//Determines how accurate and how much iteration is done by the IWCoG algorithm,
//smaller means more accurate and more iterations.
decimal iWCoGMinChange = DECIMAL(0.0002);
float iWCoGMinChange = 0.0002;

struct IWCoGParams {
int xMin;
Expand All @@ -208,16 +215,17 @@ struct IWCoGParams {
int maxIntensity;
int guess;
bool isValid;
std::unordered_set<int> checkedIndices;
std::vector<bool> checked;

};

void IWCoGHelper(IWCoGParams *p, long i, unsigned char *image, int imageWidth, int imageHeight, std::vector<int> *starIndices) {
if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && p->checkedIndices.count(i) == 0) {
if (i >= 0 && i < imageWidth * imageHeight && image[i] >= p->cutoff && !p->checked[i]) {
//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);
p->checked[i] = true;
starIndices->push_back(i);
if (image[i] > p->maxIntensity) {
p->maxIntensity = image[i];
Expand Down Expand Up @@ -246,22 +254,23 @@ void IWCoGHelper(IWCoGParams *p, long i, unsigned char *image, int imageWidth, i

Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int imageWidth, int imageHeight) const {
IWCoGParams p;
p.checked.assign(imageWidth * imageHeight, false);
int imageSize = imageWidth * imageHeight;
std::vector<Star> result;
p.cutoff = BasicThreshold(image, imageWidth, imageHeight);
for (long i = 0; i < imageHeight * imageWidth; i++) {
for (long i = 0; i < imageSize; 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
if (image[i] >= p.cutoff && !p.checked[i]) {
std::vector<int> starIndices; //indices of the current star
p.maxIntensity = 0;
int xDiameter = 0;
int yDiameter = 0;
decimal yWeightedCoordMagSum = 0;
decimal xWeightedCoordMagSum = 0;
decimal weightedMagSum = 0;
decimal fwhm; //fwhm variable
decimal standardDeviation;
decimal w; //weight value
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;
Expand All @@ -276,20 +285,20 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im
yDiameter = (p.yMax - p.yMin) + 1;

//calculate fwhm
decimal count = 0;
float count = 0;
for (int j = 0; j < (int) starIndices.size(); j++) {
if (image[starIndices.at(j)] > p.maxIntensity / 2) {
count++;
}
}
fwhm = DECIMAL_SQRT(count);
standardDeviation = fwhm / (DECIMAL(2.0) * DECIMAL_SQRT(DECIMAL(2.0) * DECIMAL_LOG(2.0)));
decimal modifiedStdDev = DECIMAL(2.0) * DECIMAL_POW(standardDeviation, 2);
// TODO: Why are these decimals? --Mark
decimal guessXCoord = (p.guess % imageWidth);
decimal guessYCoord = (p.guess / imageWidth);
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
decimal change = INFINITY;
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) {
Expand All @@ -300,24 +309,24 @@ Stars IterativeWeightedCenterOfGravityAlgorithm::Go(unsigned char *image, int im
stop++;
for (long j = 0; j < (long)starIndices.size(); j++) {
//calculate w
decimal currXCoord = starIndices.at(j) % imageWidth;
decimal currYCoord = starIndices.at(j) / imageWidth;
w = p.maxIntensity * DECIMAL_EXP(DECIMAL(-1.0) * ((DECIMAL_POW(currXCoord - guessXCoord, 2) / modifiedStdDev) + (DECIMAL_POW(currYCoord - guessYCoord, 2) / modifiedStdDev)));
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 * DECIMAL(image[starIndices.at(j)]);
yWeightedCoordMagSum += w * currYCoord * DECIMAL(image[starIndices.at(j)]);
weightedMagSum += w * DECIMAL(image[starIndices.at(j)]);
xWeightedCoordMagSum += w * currXCoord * ((float) image[starIndices.at(j)]);
yWeightedCoordMagSum += w * currYCoord * ((float) image[starIndices.at(j)]);
weightedMagSum += w * ((float) image[starIndices.at(j)]);
}
decimal xTemp = xWeightedCoordMagSum / weightedMagSum;
decimal yTemp = yWeightedCoordMagSum / weightedMagSum;
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 + DECIMAL(0.5), guessYCoord + DECIMAL(0.5), xDiameter/DECIMAL(2.0), yDiameter/DECIMAL(2.0), starIndices.size()));
result.push_back(Star(guessXCoord + 0.5f, guessYCoord + 0.5f, ((float)(xDiameter))/2.0f, ((float)(yDiameter))/2.0f, starIndices.size()));
}
}
}
Expand Down
Loading