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:

• “Optimization of Existing Centroiding Algorithms for Shack Hartmann Sensor” by Akondi Vyas, M B Roopashree and B R Prasad
• “Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor” by Akondi Vyas, M B Roopashree and B R Prasad 

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) {
 
  std::auto_ptr<FITS> pInfile(new FITS(inFilename, Read, true));
  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;
  image.read(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 {
    readFile(img, "star3.fits");
  } 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