# Night sky image processing – Part 4: Calculate the star centroid with sub-pixel accuracy

In Part 3 of my “Night sky image processing” Series I wrote about star-clustering. The result of this step was a list of stars and their pixels. The next step is to find the center of each star – the star centroid.

To keep it simple, in this example I apply the algorithm to calculate the star centroid only for one single star which I load from a FITS file. Mohammad Vali Arbabmir et. al.. proposed this algorithm in “Improving night sky star image processing algorithm for star sensors”.

Certainly it is possible to further optimize the implementation shown below. However, to me this is a good compromise between clarity and efficiency. The paper mentioned above will probably help a lot to get a better understanding of the implementation. Generally there are two steps:

### Step 1: Calculation of the intensity weighted center (IWC) (to get an idea of where star center might be)

The calculation of IWC is based on the center of gravity CoG. The CoG is the same calculation as in physics but in this context just applied to an image. Each pixel brightness value has a weight – the brighter the “heavier”. For example to get the “center” x position $x_{cog}$ (this is the x position where the image has in x-direction the “brightest” value in mean. “In mean” means that all pixel values are considered). This directly leads to the following formulas for $x_{cog}$ and $y_{cog}$:

$$x_{cog} = \frac{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1} x \cdot I(x,y) }{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1}I(x,y)}$$

$$y_{cog} = \frac{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1} y \cdot I(x,y) }{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1}I(x,y)}$$

The only thing that changes between the two formulas is $x$ and $y$. The denominator just normalized the expression so that a valid position comes out. In addition, $N$ is the image width, $Y$ is the image height and $I(x,y)$ is the pixel value at position $(x,y)$.

Now, the only difference in the calculation of IWC is that each pixel value is additionally weighted with it’s own value:

$$x_{iwc} = \frac{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1} x \cdot I^2(x,y) }{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1}I^2(x,y)}$$

$$y_{iwc} = \frac{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1} y \cdot I^2(x,y) }{\sum_{x=0}^{N-1}\sum_{y=0}^{M-1}I^2(x,y)}$$

So why should one do that? Because weighting the sum with the intensity distribution minimizes the centroid estimation error. This has some advantages and some disadvantages and it always depends on the image properties if IWC gives the best results or not. I found those two papers very helpful to get a rough idea of the different calculation methods:

### Step 2: Based on the result of step 1, iteratively apply a linear sub-pixel interpolation to (hopefully) get closer to the star centroid

We round the calculated IWC position from step 1. and use it as a base. In a few words the algorithm does the following:

• From step 1. we create a 3×3 matrix of pixels where the center pixel is the brightest one
• We divide this center pixel into 4 sub-pixels all having the same size
• Then we take for each sub-pixel the surrounding pixels (top, right, bottom, left) into account.

Then we calculate a new intensity value from those values (and of course from the sub-pixel value itself) using simple linear interpolation. This happens for each of the 4 sub-pixels. Then the brightest of the 4 sub-pixels is chosen. This pixel is then the new center and the process starts from the beginning (for example for n=10 iterations). The algorithm converges quadratically. The result is the (x,y) sub-pixel position of the star-center. A detailed description of the algorithm is presented in the paper mentioned above.

The following FITS file can be used for testing:
test3.fits – Single star in FITS format to test centroiding

### Implementation

NOTE: All the source codes of this “Night sky image processing” series presented in part 1-7 are also available on github in the star_recognizer repository.

#include <iostream>
#include <assert.h>
#include <CImg.h>
#include <CCfits/CCfits>

using namespace cimg_library;
using namespace CCfits;
using namespace std;

// See https://heasarc.gsfc.nasa.gov/fitsio/ccfits/html/cookbook.html
void readFile(CImg<float> & inImg, const string & inFilename) {

PHDU & image = pInfile->pHDU();
inImg.resize(image.axis(0) /*x*/, image.axis(1) /*y*/, 1/*z*/, 1 /*1 color*/);

// NOTE: At this point we assume that there is only 1 layer.
std::valarray<unsigned long> imgData;
cimg_forXY(inImg, x, y) { inImg(x, inImg.height() - y - 1) = imgData[inImg.offset(x, y)]; }
}

float calcIx2(const CImg<float> & img, int x) {
float Ix = 0;
cimg_forY(img, y) { Ix += pow(img(x, y), 2.0) * (float) x; }
return Ix;
}

float calcJy2(const CImg<float> & img, int y) {
float Iy = 0;
cimg_forX(img, x) { Iy += pow(img(x, y), 2.0) * (float) y; }
return Iy;
}

// Calculate Intensity Weighted Center (IWC)
void calcIntensityWeightedCenter(const CImg<float> & inImg, float * outX, float * outY) {
assert(outX && outY);
assert(inImg.width() == inImg.height());
const size_t L = inImg.width();

// Determine weighted centroid - See https://cdn.intechopen.com/pdfs-wm/26716.pdf
float Imean2 = 0, Jmean2 = 0, Ixy2 = 0;

for(size_t i = 0; i < L; ++i) {
Imean2 += calcIx2(inImg, i);
Jmean2 += calcJy2(inImg, i);
cimg_forY(inImg, y) { Ixy2 += pow(inImg(i, y), 2.0); }
}

*outX = Imean2 / Ixy2;
*outY = Jmean2 / Ixy2;
}

void calcSubPixelCenter(const CImg<float> & inImg, float * outX, float * outY,
size_t inNumIter = 10 /*num iterations*/) {
// Sub pixel interpolation
float c, a1, a2, a3, a4, b1, b2, b3, b4;
float a1n, a2n, a3n, a4n, b1n, b2n, b3n, b4n;

assert(inImg.width() == 3 && inImg.height() == 3);

b1 = inImg(0, 0); a2 = inImg(1, 0); b2 = inImg(2, 0);
a1 = inImg(0, 1);  c = inImg(1, 1); a3 = inImg(2, 1);
b4 = inImg(0, 2); a4 = inImg(1, 2); b3 = inImg(2, 2);

for (size_t i = 0; i < inNumIter; ++i) {
float c2 = 2 * c;
float sp1 = (a1 + a2 + c2) / 4;
float sp2 = (a2 + a3 + c2) / 4;
float sp3 = (a3 + a4 + c2) / 4;
float sp4 = (a4 + a1 + c2) / 4;

// New maximum is center
float newC = std::max({ sp1, sp2, sp3, sp4 });

// Calc position of new center
float ad = pow(2.0, -((float) i + 1));

if (newC == sp1) {
*outX = *outX - ad; // to the left
*outY = *outY - ad; // to the top

// Calculate new sub pixel values
b1n = (a1 + a2 + 2 * b1) / 4;
b2n = (c + b2 + 2 * a2) / 4;
b3n = sp3;
b4n = (b4 + c + 2 * a1) / 4;
a1n = (b1n + c + 2 * a1) / 4;
a2n = (b1n + c + 2 * a2) / 4;
a3n = sp2;
a4n = sp4;

} else if (newC == sp2) {
*outX = *outX + ad; // to the right
*outY = *outY - ad; // to the top

// Calculate new sub pixel values
b1n = (2 * a2 + b1 + c) / 4;
b2n = (2 * b2 + a3 + a2) / 4;
b3n = (2 * a3 + b3 + c) / 4;
b4n = sp4;
a1n = sp1;
a2n = (b2n + c + 2 * a2) / 4;
a3n = (b2n + c + 2 * a3) / 4;
a4n = sp3;
} else if (newC == sp3) {
*outX = *outX + ad; // to the right
*outY = *outY + ad; // to the bottom

// Calculate new sub pixel values
b1n = sp1;
b2n = (b2 + 2 * a3 + c) / 4;
b3n = (2 * b3 + a3 + a4) / 4;
b4n = (2 * a4 + b4 + c) / 4;
a1n = sp4;
a2n = sp2;
a3n = (b3n + 2 * a3 + c) / 4;
a4n = (b3n + 2 * a4 + c) / 4;
} else {
*outX = *outX - ad; // to the left
*outY = *outY + ad; // to the bottom

// Calculate new sub pixel values
b1n = (2 * a1 + b1 + c) / 4;
b2n = sp2;
b3n = (c + b3 + 2 * a4) / 4;
b4n = (2 * b4 + a1 + a4) / 4;
a1n = (b4n + 2 * a1 + c) / 4;
a2n = sp1;
a3n = sp3;
a4n = (b4n + 2 * a4 + c) / 4;
}

c = newC; // Oi = Oi+1

a1 = a1n;
a2 = a2n;
a3 = a3n;
a4 = a4n;

b1 = b1n;
b2 = b2n;
b3 = b3n;
b4 = b4n;
}
}

int main(void) {
CImg<float> img;
float xc, yc;

try {
} catch (FitsException &) {
cerr << "Read FITS failed." << endl;
return 1;
}

// 1. Calculate the IWC
calcIntensityWeightedCenter(img, & xc, & yc);

// 2. Round xc, yc to nearest integer and then iteratively improve.
int xi = floor(xc + 0.5);
int yi = floor(yc + 0.5);

CImg<float> img3x3 = img.get_crop(xi - 1 /*x0*/, yi - 1 /*y0*/, xi + 1 /*x1*/, yi + 1 /*y1*/);

// 3. Interpolate using sub-pixel algorithm
float xsc = xi, ysc = yi;
calcSubPixelCenter(img3x3, & xsc, & ysc, 10 /*num iterations*/);

cerr << "xc: " << xc << " --> xi: " << xi << ", yc: " << yc << " --> yi: " << yi << endl;
cerr << "xsc: " << xsc << ", ysc: " << ysc << endl;

return 0;
}

Here is the command to compile the above code:

g++ -o sub_pixel_det -std=c++0x sub_pixel_det.C -lCCfits -lcfitsio -lX11 -lpthread

Result:

@:~/sub_pixel_det\$ ./sub_pixel_det
xc: 2.56769 --> xi: 3, yc: 2.51555 --> yi: 3
xsc: 2.55957, ysc: 2.0752

In Part 5 of my “Night sky image processing” Series I am going to write about curve fitting and determination of the FWHM value of a star. As a base the centroid calculated here will be used.

Last updated: June 9, 2022 at 23:43 pm

1. Bjarne Lassen says:

hi Carsten
My name is Bjarne. I’m from Denmark.
Thank’s for a very inspirering articles.. I’m not very famiiar with C programming.. so I tried to rewrite this program in to Delphi..
My first question is: running the intensetyweightedCenter routine show that the first line(x) of the InImg and the first coloum(Y) becomes a very low number due to the fact that the upper fraction of the sum is zero as
pow(inImg[X,Y) *X / inImg[x,Y] where X := zero .. same with Y.
Is that intentionly ? second : I can not grasp the function of the
float ad = pow(2.0, -((float) i + 1)) .. calculation seems that Ad become the 2.root of 2 and 3.root of 2…and so on? is that right? and for what purpose?
Bjarne
hope You can understand my questions.

1. Hi Bjarne,

thank you for your questions. I am happy that the content on my blog is useful to you. The code you refer to I wrote ~9 years ago, but I will try to remember what I did as good as I can.

Regarding your first question, I remember that I found the calculation of the “Intensity Weighted Centroiding” (IWC) in this paper: https://cdn.intechopen.com/pdfs-wm/26716.pdf (page 5 (170), formula (5)). Looking at this formula again, I think you found a bug. By starting the calculation with index 0, the first line and the first column are indeed “ignored”. The paper does not talk about the index by which the iteration should start. This did not happen on purpose. I should correct this as soon as I find the time. Just as a hint: In the paper which I attached, I just saw that the IWC is also calculated – but the way it is done is slightly different. There, the index start indeed with 1. I did not check if in the end both methods from both papers are equal.

As a bottom line, I think you found a problem. For typical star images I guess the impact of ignoring the first column and the first line is probably small – but I have not fully thought this through, yet. Still, in the end the currently calculated value is not correct. So thanks, for the hint (y)

For your second question, the paper “Improving night sky star image processing algorithm for star sensors” on page 4 (797) provides the answer. The full algorithm to calculate the centroid based on sub-pixels is described there. I guess, unfortunately some parts of the code are not fully understandable without having this paper. To my knowledge this paper is unfortunately not freely available. Therefore, I cannot provide it here and I also cannot link to it.

Step (b) of this algorithm talks about the calculation of “ad”. As far as I understood this, for each iteration the selected pixel is divided into 4 smaller sub-pixels. “ad” is the distance by which the current centroid position (x_c, y_c) is changed to the new, more accurate centroid position (x_c_new, y_c_new) after each iteration. The algorithm always chooses the brightest of the four new sub-pixels and adds/subtracts “ad” to the centroid position accordingly.

By again dividing the selected pixel into 4 sub-pixels, the distance by which the final centroid position changes, gets smaller with each iteration. If you calculate “ad” for the first couple of iterations, you get

i=0 -> 2^-(i+1) = 1/2^(i+1) = 1/2^1 = 0.5
i=1 -> 0.25
i=2 -> 0.125
i=3 -> 0.0625

One more hint: I think there is a mistake in the paper in step (b). The way “ad” is used to calculate the next centroid position (x_c, y_c) is wrong for sp_2 and for sp_4. I once wrote a mail to the author but never got a confirmation. At least in my code, I changed the calculation of (x_c_new, y_x_new) for the two cases when sp_2 or sp_4 are the brightest pixels (compared to what is described in the paper). Let me know what you think about this!

I will probably reshape this explanation a little and post it without the hint to the paper on the blog (as reply to your questions). So hopefully it will be helpful for others in the future.

Hope this helps!

Best regards
Carsten