Skip to content
2020
Introduction
Book Summaries
Projects
Other Readings/ Learning/Articles
Courses
More
Share
Explore
Detection of Pulsars using Median Image Stacking (Binapprox)

icon picker
Mean Image Stacking

Implementing mean image stacking is easy. I implemented it in two ways. One where i assumed i knew the number of images to process and second where i imagined a stream of images instead (running mean calculations)
‌

Implementation where we know the number of images to process.

‌
If the number of images are small, you can very well bring them all in memory and then sum each cell up and divide the total of each cell by number of images. If the number is large, you can load them each one by one, keep a matrix for the sum and in the end again divide by the total by the number of images.
‌

Implementation where we assume it to be a stream of images.

‌
Here we will use "Welford's online algorithm" to calculate the mean. It is an easy to implement algorithm as we will see below. This algorithm will also be used in binapprox to optimise the running mean and finding std deviation there.
‌
(Here Basic knowledge of python is assumed)
‌

Some common functions to implement first

‌

Loading the fits file

‌
We will use astropy library to load fits:
from astropy.io import fits
def load_fits (fname) :
hdulist = fits.open(fname)
data = np.array(hdulist[0].data)
hdulist.close() #closes the file
return data;
‌

Plotting the fits data

‌
We will use matplotlib library.
import matplotlib.pyplot as plt
def plotImage (data) :
plt.imshow(data, cmap=plt.cm.viridis)
plt.xlabel('x-pixels (RA)')
plt.ylabel('y-pixels (Dec)')
plt.colorbar()
plt.title
plt.show()
‌

Simple Solution

def mean_datasets (datasetArr):
final = 0
for data in datasetArr:
final += data
size = len(datasetArr)
final /= size
return final
def mean_fits (fnameArr) : #fnameArr is the array containing filenames
dataArr = []
for fname in fnameArr :
dataArr.append(load_fits(fname))
mean = mean_datasets(dataArr)
return mean
‌

Welford's online algorithm

‌
Let Mo be the mean of already processed values (n-1) , then with the new incoming value x, the new mean (Mn) would be
Mn = ( (Mo ∗ (n−1) ) + x) / n
=
Mo+(x−Mo)/nMn
=
((Mo * (n-1)) + x ) / n
=
Mo + (x-Mo)/nMn
=
((Mo∗(n−1))+x)/n
=
Mo+(x−Mo)/n

def calcMeanStdDev (imageArr):
mean = 0
count = 0
for fname in imageArr :
data = load_fits(fname);
mat = np.array(data)
count += 1
delta = mat - mean;
mean += delta/count;
return mean;
Implementation where we know the number of images to process.
Implementation where we assume it to be a stream of images.
Some common functions to implement first
Loading the fits file
Plotting the fits data
Simple Solution
Welford's online algorithm
Want to print your doc?
This is not the way.
Try clicking the ⋯ next to your doc name or using a keyboard shortcut (
CtrlP
) instead.