# -*- coding: utf-8 -*-
"""
Created on Nov 2 2014
@author: Erlend Ronnekleiv
This code is written for Pythonxy v2.7.
The Spyder GUI can be used to run the program:
https://code.google.com/p/spyderlib/ (download version for python 2.7)
If you googe python you will find a lot of tutorials and Q&As
about how to use this language.
You probably have to download some missing libraries to get the code working.
For this you can use pip. In your system command line interface
(cmd.exe under windows) run pipto install missing packages.
Example: pip install pyfits
Suggested workflow:
1) Set file directory and name: dir0 and fname_in
2) Set parameters for starFilter to match the recorded star shape.
3) Adjust parameters for nebFilter to have a radius > detected star radius
for proper interpolation.
4) Adjust starThr to a few times the image st.dev. after applying starFilter
Adjust starThrBgrFactor to represent the increase in st.dev with increased
nebula bacground.
5) Set initial unsharp mask parameters usmr and usmw. Unsharpening allows the
stars to be shapened to take up a smaller diameter before thay are removed.
6) Run the program, review results, adjust parameters, and run again until
results are acceptable.
7) Open _StarsRemoved.fit and _usm.fit in an image editor
(depending on the editor you may first have to convert to 16 bit tiff).
Review the difference. Sharp nebula details may have been removed in the
process. This can be repaired by "cloning" back the original usm image and
manually retouching stars in those regions.
Note: I wrote this program to use it myself. It is not very user friendly,
allthough I have added a number of comments throughout the code. I would be
pleased to help if someone wants to make the program more user friendly.
Please contact erlend.web@eronn.net
"""
from __future__ import division
import numpy as np
import pyfits
from numpy.fft import rfft2, irfft2, ifftshift
import time
import scipy.ndimage as ndimg
#import matplotlib.pyplot as pl
#import PIL
#import scipy.special
#import matplotlib.cm as cm # colormaps
#import scipy.misc.imsave as imsave
nx=np.newaxis
def cosap2d(r1,r2):
"""
Cosine apodised circular filter response. Flat inside ratius r1.
Apodized down to zero at radius r2.
"""
a = int(np.ceil(r2))-1
x = np.arange(-a,a+1);
r = np.sqrt(x**2 + x[:,nx]**2)
ap = np.zeros(np.shape(r))
ap[r<=r1]=1;
ii = (r>r1)&(rr3] = 0
return filt/np.sum(filt)
def star(r1,r2):
a = int(np.ceil(r2))-1
x = np.arange(-a,a+1);
ap = np.zeros(np.shape(x))
ap[x<=r1]=1;
ii = (x>r1)&(x