For some time now I am looking into the field of night sky image processing. It is a huge topic and there are already a lot of different approaches and solutions to many of the existing problems. I have collected many star images over the time. To me it is very interesting to try different algorithms with my own data. This article starts with a problem every night sky photographer and astronomer is aware of: noise. I did some research and came across “anisotropic diffusion”.

What is anisotropic diffusion?

In short, anisotropic diffusion is a non-linear and space-variant transformation i.e. the transformation depends on the image content. The effect is that the resulting images preserve linear structures while at the same time smoothing is made along these structures. Additional information on anisotropic diffusion is available for example on wikipedia and here. Furthermore I found this youtube explanation quite helpful.

A simple C++ 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.

I was looking for an easy way to apply this algorithm to a FITS file using C++. Finally I found one using the CImg library which already contains an implementation for anisotropic diffusion. For reading and writing the FITS file I am using the CCFits library. The listing below shows how I got it to work. Furthermore the image below shows the effect of the code.

#include <iostream>
#include <CImg.h>
#include <CCfits/CCfits>
 
using namespace cimg_library;
using namespace CCfits;
using namespace std;
 
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)]; }  
}
 
void writeFile(const CImg<float> & inImg, const string & inFilename) {
  long naxis = 2;
  long naxes[2] = { inImg.width(), inImg.height() };
  std::auto_ptr<FITS> pFits(0);
  pFits.reset(new FITS(string("!") + inFilename , USHORT_IMG , naxis , naxes) );
  
  // NOTE: At this point we assume that there is only 1 layer.
  long nelements = std::accumulate(& naxes[0], & naxes[naxis], 1, std::multiplies<long>());
  std::valarray<int> array(nelements);
  cimg_forXY(inImg, x, y) { array[inImg.offset(x, y)] = inImg(x, inImg.height() - y -1); }
 
  long fpixel(1);
  pFits->pHDU().write(fpixel, nelements, array);
}
 
int main(void) {
  CImg<float> img;
 
  try {
    readFile(img, "test.fit");
  } catch (FitsException &) {
    cerr << "Read FITS failed." << endl;
    return 1;
  }
 
  CImgDisplay imgDisp(img, "Original image");
  while (! imgDisp.is_closed()) {
    imgDisp.wait();
  }
  
  // https://cimg.sourceforge.net/reference/structcimg__library_1_1CImg.html
  CImg<float> & resImg = img.blur_anisotropic(30.0f, /*amplitude*/
                                               0.7f, /*sharpness*/
                                               0.3f, /*anisotropy*/
                                               0.6f, /*alpha*/
                                               1.1f, /*sigma*/
                                               0.8f, /*dl*/
                                               30,   /*da*/
                                               2,    /*gauss_prec*/
                                               0,    /*interpolation_type*/
                                               false /*fast_approx*/
                                               );
 
  CImgDisplay imgDisp2(resImg, "Result image");
  while (! imgDisp2.is_closed()) {
    imgDisp2.wait();
  }
  
  try {
    writeFile(resImg, "test_res.fit");
  } catch (FitsException &) {
    cerr << "Writing FITS failed." << endl;
    return 1;
  }
  
  return 0;
}

The file is also available for download here. The test image test.fits can be downloaded here (~17MB). The following command I used to build it:

g++ anisotropic_diffusion_cimg.C -o anisotropic_diffusion_cimg -lCCfits -lX11 -lpthread -lcfitsio

Thanks for reading. I hope this was helpful in a way. In Part 2 of my “Night sky image processing” series I will have a closer look onto image thresholding.

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