@@ -26,3 +26,6 @@ out.fits
 dirty.fits
 /doc
 /.ipynb_checkpoints
+psf.fits
+.DS_Store
+/data
@@ -0,0 +1,14 @@
+Crocodile contributors (sorted alphabetically)
+==============================================
+
+* **[Tim Cornwell](http://github.com/timcornwell)**
+
+  * Theory, Design, Implementation
+
+* **[Bojan Nikolic](http://github.com/bnikolic)**
+
+  * Conception, Theory, Design, Implementation
+  
+* **[Peter Wortmann](http://github.com/scpmw)**
+
+  * Theory, Design, Implementation
@@ -4,20 +4,23 @@
 # subclasses of astropy classes.
 #
 
-import matplotlib.pyplot as plt
+import os
 
+import matplotlib.pyplot as plt
 import numpy as numpy
-
-from astropy.wcs import WCS
+from astropy.coordinates import SkyCoord
 from astropy.io import fits
+from astropy.wcs import WCS
 from astropy.wcs.utils import skycoord_to_pixel, pixel_to_skycoord
-from astropy.coordinates import SkyCoord
+from reproject import reproject_interp
+
 from arl.skycomponent import SkyComponent, create_skycomponent
 
 """
 Functions that define and manipulate images. Images are just data and a World Coordinate System.
 """
 
+
 class Image:
     """Image class with Image data (as a numpy.array) and optionally the AstroPy WCS.
     
@@ -31,13 +34,29 @@ class Image:
     - in numpy, the order is (frequency, polarisation, latitude, longitude)
     
     """
-
+    
     def __init__(self):
         self.data = None
         self.wcs = None
 
 
-def image_show(im: Image, fig = None, title: str = ''):
+def create_test_image(canonical=True):
+    """Create a useful test image
+    
+    This is the test image M31 widely used in ALMA and other simulations. It is actually part of an Halpha region in
+    M31.
+    
+    :param canonical: Make the image into a 4 dimensional image
+    :returns: Image
+    """
+    chome = os.environ['CROCODILE']
+    im = image_from_fits("%s/data/models/M31.MOD" % chome)
+    if canonical:
+        im=image_replicate(im)
+    return im
+
+
+def image_show(im: Image, fig=None, title: str = ''):
     """ Show an Image with coordinates using matplotlib
     
     :param im:
@@ -52,7 +71,11 @@ def image_show(im: Image, fig = None, title: str = ''):
         fig = plt.figure()
     fig.add_subplot(111, projection=im.wcs.sub(['longitude', 'latitude']))
     plt.clf()
-    plt.imshow(im.data[0, 0, :, :], origin='lower', cmap='rainbow')
+    if len(im.data.shape) == 4:
+        print(im.data[0, 0, :, :])
+        plt.imshow(im.data[0, 0, :, :], origin='lower', cmap='rainbow')
+    elif len(im.data.shape) == 2:
+        plt.imshow(im.data[:, :], origin='lower', cmap='rainbow')
     plt.xlabel('RA---SIN')
     plt.ylabel('DEC--SIN')
     plt.title(title)
@@ -95,7 +118,7 @@ def image_to_fits(im: Image, fitsfile: str = 'imaging.fits'):
     :param fitsfile: Name of output fits file
     :type str:
     """
-    return fits.writeto(fitsfile, im.data, im.wcs.to_header(), clobber=True)
+    return fits.writeto(filename=fitsfile, data=im.data, header=im.wcs.to_header(), clobber=True)
 
 
 def image_from_fits(fitsfile: str):
@@ -109,11 +132,12 @@ def image_from_fits(fitsfile: str):
     if fitsfile[0] == '.':
         import os
         chome = os.environ['CROCODILE']
-        fitsfile="%s/%s" % (chome, fitsfile)
+        fitsfile = "%s/%s" % (chome, fitsfile)
     hdulist = fits.open(fitsfile)
     fim = Image()
     fim.data = hdulist[0].data
     fim.wcs = WCS(fitsfile)
+    hdulist.close()
     print("image_from_fits: Max, min in %s = %.6f, %.6f" % (fitsfile, fim.data.max(), fim.data.min()))
     return fim
 
@@ -131,7 +155,7 @@ def image_add_wcs(im: Image, wcs: WCS):
     return im
 
 
-def image_replicate(im: Image, shape=None):
+def image_replicate(im: Image, shape=None, frequency=1.4e9):
     """ Make a new canonical shape Image, extended along third and fourth axes by replication.
     
     The order is [chan, pol, dec, ra]
@@ -143,46 +167,47 @@ def image_replicate(im: Image, shape=None):
     :type 4-sequence:
     :returns: Image
     """
-    # TODO: Fill in extra axes in wcs
-    if shape is None:
-        shape = [1, 1, 1, 1]
+    if shape == None:
+        shape=[1,1,1,1]
+
     if len(im.data.shape) == 2:
         fim = Image()
-        image_add_wcs(fim, im.wcs)
+
+        newwcs = WCS(naxis=4)
+
+        newwcs.wcs.crpix = [im.wcs.wcs.crpix[0], im.wcs.wcs.crpix[1], 1.0, 1.0]
+        newwcs.wcs.cdelt = [im.wcs.wcs.cdelt[0], im.wcs.wcs.cdelt[1], 1.0, 1.0]
+        newwcs.wcs.crval = [im.wcs.wcs.crval[0], im.wcs.wcs.crval[1], 1.0, frequency]
+        newwcs.wcs.ctype = [im.wcs.wcs.ctype[0], im.wcs.wcs.ctype[1], 'STOKES', 'FREQ']
+
+        image_add_wcs(fim, newwcs)
         fshape = [shape[3], shape[2], im.data.shape[1], im.data.shape[0]]
         fim.data = numpy.zeros(fshape)
         print("image_replicate: replicating shape %s to %s" % (im.data.shape, fim.data.shape))
         for i3 in range(shape[3]):
             for i2 in range(shape[2]):
                 fim.data[i3, i2, :, :] = im.data[:, :]
-        print(fim.wcs)
-    elif len(im.data.shape) == 3:
-        fim = Image()
-        image_add_wcs(fim, im.wcs)
-        # TODO: fix this for non-square images!
-        fshape = [shape[3], im.data.shape[2], im.data.shape[1], im.data.shape[0]]
-        fim.data = numpy.array(fshape)
-        print("Image: replicating shape %s to %s" % (im.data.shape, fim.data.shape))
-        for i3 in range(shape[3]):
-            fim.data[i3, :, :, :] = im.data[:, :, :]
     else:
         fim = im
-
+    
     return fim
 
 
 def image_reproject(im: Image, newwcs: WCS, shape=None):
     """ Reproject an image to a new coordinate system
+    
+    Currently uses the reproject python package.
+    TODO: Write tailored reproject routine
 
-    :param im:
+    :param im: Image to be reprojected
     :type Image:
     :param newwcs: New WCS
     :type WCS:
-    :returns: Image
+    :returns: Reprojected Image, Footprint Image
     """
-    print("image_reproject: Not yet implemented")
-    
-    return im
+    rep, foot = reproject_interp((im.data, im.wcs), newwcs, shape, order='bicubic',
+                                 independent_celestial_slices=False)
+    return image_from_array(rep, newwcs), image_from_array(foot, newwcs)
 
 
 def image_FFT(im: Image, **kwargs):
@@ -211,6 +236,7 @@ def image_add(im1: Image, im2: Image, checkwcs=False):
     assert not checkwcs, "Checking WCS not yet implemented"
     return image_from_array(im1.data + im2.data, im1.wcs)
 
+
 def point_source_find(im: Image, **kwargs) -> SkyComponent:
     """ Find components in Image, return SkyComponent, just find the peak for now
     
@@ -220,14 +246,14 @@ def point_source_find(im: Image, **kwargs) -> SkyComponent:
     """
     # TODO: Implement full image fitting of components
     print("imaging.point_source_find: Finding components in Image")
-
+    
     # Beware: The index sequencing is opposite in wcs and Python!
     locpeak = numpy.array(numpy.unravel_index((numpy.abs(im.data)).argmax(), im.data.shape))
     print("imaging.point_source_find: Found peak at pixel coordinates %s" % str(locpeak))
     w = im.wcs.sub(['longitude', 'latitude'])
     sc = pixel_to_skycoord(locpeak[3], locpeak[2], im.wcs, 0, 'wcs')
     print("imaging.point_source_find: Found peak at world coordinates %s" % str(sc))
-    flux=im.data[:,:,locpeak[2],locpeak[3]]
+    flux = im.data[:, :, locpeak[2], locpeak[3]]
     print("imaging.point_source_find: Flux is %s" % flux)
     # We also need the frequency values
     w = im.wcs.sub(['spectral'])
@@ -249,18 +275,20 @@ def flux_at_direction(im: Image, sc: SkyCoord, **kwargs) -> SkyComponent:
     w = im.wcs.sub(['longitude', 'latitude'])
     pixloc = skycoord_to_pixel(sc, im.wcs, 0, 'wcs')
     print("imaging.flux_at_direction: Extracting flux at pixel coordinates %d %d" % (pixloc[0], pixloc[1]))
-    flux=im.data[:,:,int(pixloc[1]+0.5),int(pixloc[0]+0.5)]
+    flux = im.data[:, :, int(pixloc[1] + 0.5), int(pixloc[0] + 0.5)]
     print("imaging.flux_at_direction: Flux is %s" % flux)
-
+    
     # We also need the frequency values
     w = im.wcs.sub(['spectral'])
     frequency = w.wcs_pix2world(range(im.data.shape[0]), 0)
-
+    
     return create_skycomponent(direction=sc, flux=flux, frequency=frequency, shape='point')
 
+
 if __name__ == '__main__':
     import os
-    chome=os.environ['CROCODILE']
+    
+    chome = os.environ['CROCODILE']
     kwargs = {}
     m31model = image_from_fits("%s/data/models/M31.MOD" % chome)
     m31model_by_array = image_from_array(m31model.data, m31model.wcs)
 
@@ -2,9 +2,6 @@
 
 import numpy
 
-from .synthesis import sortw, doimg, wslicimg, wslicfwd
-
-
 def overlapIndices(a1, a2,
                    shiftx, shifty):
     """ Find the indices where two arrays overlapIndices
@@ -57,18 +54,12 @@ def hogbom(dirty,
     Hogbom CLEAN (1974A&AS...15..417H)
 
     :param dirty: The dirty Image, i.e., the Image to be deconvolved
-
     :param psf: The point spread-function
-
     :param window: Regions where clean components are allowed. If True, all of the dirty Image is assumed to be
     allowed for clean components
-
     :param gain: The "loop gain", i.e., the fraction of the brightest pixel that is removed in each iteration
-
     :param thresh: Cleaning stops when the maximum of the absolute deviation of the residual is less than this value
-
     :param niter: Maximum number of components to make if the threshold `thresh` is not hit
-
     :returns clean SkyComponent Image, residual Image
     """
 
@@ -93,34 +84,3 @@ def hogbom(dirty,
         if numpy.fabs(res).max() < thresh:
             break
     return comps, res
-
-
-def majorcycle(T2, L2,
-               p, v,
-               gain,
-               nmajor,
-               nminor,
-               wstep):
-    """Major cycle clean
-
-    :param v:
-    :param gain:
-    :param T2: Field of view in radians
-    :param L2: Observing wavelength (m)
-    :param p: UVWs of visiblities (m)
-    :param nmajor: Number of major cycles
-    :param nminor: Number of minor cycles
-    :param wstep: Step in w (pixels)
-    
-    NOT USED
-
-    """
-    ps, vs = sortw(p, v)
-    for i in range(nmajor):
-        dirty, psf, _ = doimg(T2, L2, ps, vs, lambda *x: wslicimg(*x, wstep=wstep, Qpx=1))
-        cc, rres = hogbom(dirty, psf, True, gain, 0,
-                          nminor)
-        xuv = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.ifftshift(cc)))
-        ps, vsp = wslicfwd(xuv, T2, L2, p, wstep=wstep)
-        vs = vs - vsp
-    return ps, vs, cc, res
@@ -6,63 +6,6 @@
 from crocodile.clean import overlapIndices, argmax
 from crocodile.synthesis import sortw, doimg, dopredict, simplepredict, simpleimg
 
-
-def majorcycle(T2, L2,
-               p, v,
-               gain,
-               nmajor,
-               nminor,
-               scales,
-               thresh=0.0,
-               fracthresh=0.1,
-               predfn=simplepredict,
-               imgfn=simpleimg):
-    """Major cycle for MultiScale Clean
-
-    :param gain:
-    :param scales:
-    :param thresh:
-    :param fracthresh:
-    :param predfn:
-    :param imgfn:
-    :param T2: Field of view in radians
-    :param L2: Maximum uv (for setting Image pixel sizes)
-    :param p: UVWs of visibilities (m)
-    :param v: Values of visibilities
-    :param nmajor: Number of major cycles
-    :param nminor: Number of minor cycles
-    :param thresh: Stopping threshold (for scale=0)
-    :param fracthresh: Minor Cycle stopping threshold (for scale=0) fraction of peak
-    
-    NOT USED
-    """
-
-    # The model is added to each major cycle and then the visibilities are
-    # calculated from the full model
-    global vsp
-    ps, vso = sortw(p, v)
-    dirty, psf, sumwt = doimg(T2, L2, ps, vso, imgfn)
-    comps = 0.0 * dirty.copy()
-    for i in range(nmajor):
-        print("Start of major cycle %d" % i)
-        cc, res = msclean(dirty, psf, True, gain, thresh, nminor, scales, fracthresh)
-        plt.clf()
-        plt.imshow(res, cmap='rainbow', origin='lower')
-        plt.colorbar()
-        plt.show()
-        comps += cc
-        # dopredict resorts the data
-        pss, vsp = dopredict(T2, L2, ps, comps, predfn)
-        vsr = vso - vsp
-        dirty, psf, sumwt = doimg(T2, L2, ps, vsr, imgfn)
-        if np.abs(dirty).max() < 1.1 * thresh:
-            print("Reached stopping threshold %.6f Jy" % thresh)
-            break
-        print("End of major cycle")
-    print("End of major cycles")
-    return ps, vsp, comps, dirty
-
-
 def msclean(dirty,
             psf,
             window,