#!/usr/bin/env python3


import os
import sys
import pyfits
import numpy
import scipy
import scipy.ndimage


if __name__ == "__main__":

    fn = sys.argv[1]
    median_file = sys.argv[2]
    diff_file = sys.argv[3]

    hdulist = pyfits.open(fn)
    img_data = hdulist[0].data
    print(img_data.shape)

    # now smooth the image to subtract large(ish)-scale features
    print("Applying median-filter for background removal")
    if (os.path.isfile(median_file)):
        h = pyfits.open(median_file)
        large_scale = h[0].data
    else:
        large_scale = scipy.ndimage.median_filter(img_data, cval=0, size=5)
        pyfits.PrimaryHDU(data=large_scale, header=hdulist[0].header).writeto(median_file, clobber=True)

    print(large_scale.shape)

    diff = img_data - large_scale
    pyfits.PrimaryHDU(data=diff, header=hdulist[0].header).writeto(
        diff_file, clobber=True)