def main():
    '''Input a set of R G B images and return a color image'''

    # Setup command line options
    parser = argparse.ArgumentParser(description='Coadd R,G,B fits images and return color png')
    parser.add_argument("rFile", help='Input R image', default='r.fits')
    parser.add_argument("gFile", help='Input G image', default='g.fits')
    parser.add_argument("bFile", help='Input B image', default='b.fits')
    parser.add_argument('--addNoise', dest="noise", nargs='+', type=float, help='Gaussian noise to add to each image', default=[0,0,0])
    parser.add_argument('--outputFile', dest="outputFile", help='Output PNG file name', default='rgb.png')
    args = parser.parse_args()

    images = []
    images.append(pyfits.getdata(args.rFile,0))
    images.append(pyfits.getdata(args.gFile,0))
    images.append(pyfits.getdata(args.bFile,0))

    # add noise
    if (args.noise != [0,0,0]):
        addNoise(images,args.noise)

    # scale image
    scaledImages = scaleAsinh(images, scales =[1.,1.,1.])

    # create RGB
    mode='RGB'
    pngImage = Image.new(mode, scaledImages[0].shape)
    pngImage.paste(createRGB(scaledImages),(0,0))

    pngImage.save(args.outputFile)

def plot_fits_reg_vs_out(fits_dir, regular, outliers, objids2name):
    reg_fits = [os.path.join(fits_dir, objids2name[p]) for p in regular]
    outl_fits = [os.path.join(fits_dir, objids2name[p]) for p in outliers]
    reg_fits_data = [pyfits.getdata(fit) for fit in reg_fits]
    outl_fits_data = [pyfits.getdata(fit) for fit in outl_fits]
    plot_fits_by_size(reg_fits_data, name='-reg')
    plot_fits_by_size(outl_fits_data, name='-outl')

def make_voronoi_intens(targetSN, w1, w2):
    """ Make image"""
    image = "collapsed_w{0}_{1}.fits".format(w1, w2)
    intens = pf.getdata(image)
    extent = calc_extent(image)
    vordata = pf.getdata("voronoi_sn{0}_w{1}_{2}.fits".format(targetSN, w1,
                                                              w2))
    vordata = np.ma.array(vordata, mask=np.isnan(vordata))
    bins = np.unique(vordata)[:-1]
    combined = np.zeros_like(intens)
    combined[:] = np.nan
    for j, bin in enumerate(bins):
        idx, idy = np.where(vordata == bin)
        flux = intens[idx,idy]
        combined[idx,idy] = np.nanmean(flux)
    vmax = np.nanmedian(intens) + 4 * np.nanstd(intens)
    fig = plt.figure(1)
    plt.minorticks_on()
    make_contours()
    plt.imshow(combined, cmap="cubehelix_r", origin="bottom", vmax=vmax,
                    extent=extent, vmin=0)
    plt.xlabel("X [kpc]")
    plt.ylabel("Y [kpc]")
    cbar = plt.colorbar()
    cbar.set_label("Flux [$10^{-20}$ erg s$^{-1}$ cm$^{-2}$]")
    plt.savefig("figs/intens_sn{0}.png".format(targetSN), dpi=300)
    pf.writeto("figs/intens_sn{0}.fits".format(targetSN), combined,
               clobber=True)
    return

def fitsread(imgname, header = False):
    """
    Read CSUSB telescope FITS image cube.
    
    Parameters
    ----------
    image : string
        FITS image name
        
    header : boolean
        Return FITS image header?
        
    Returns
    -------
    img_data : numpy array
        2D or 3D numpy array
    """
    try:
        if header:
            img_data, header = pyfits.getdata(imgname, ignore_missing_end = True, header = True)
            return img_data, header
        else:
            img_data = pyfits.getdata(imgname, ignore_missing_end = True)
            return img_data
    except IOError:
        print "FITSREAD: Unable to open FITS image %s" %imgname
    
    return

def mk_image(galaxy):
    base = './../../images_v5/GS_2.5as_matched/gs_all_'

    i_img = pyf.getdata(base+str(galaxy)+'_I.fits')
    j_img = pyf.getdata(base+str(galaxy)+'_J.fits')
    h_img = pyf.getdata(base+str(galaxy)+'_H.fits')

    #include 90% of pixels
    x = pyl.hstack(i_img)
    i_lim = scoreatpercentile(x,99)
    x = pyl.hstack(j_img)
    j_lim = scoreatpercentile(x,99)
    x = pyl.hstack(h_img)
    h_lim = scoreatpercentile(x,99)

    print galaxy, i_lim, j_lim, h_lim

    img = pyl.zeros((h_img.shape[0], h_img.shape[1], 3), dtype=float)
    img[:,:,0] = img_scale.asinh(h_img, scale_min=-0.1*h_lim, scale_max=h_lim,
            non_linear=0.5)
    img[:,:,1] = img_scale.asinh(j_img, scale_min=-0.1*j_lim, scale_max=j_lim,
            non_linear=0.5)
    img[:,:,2] = img_scale.asinh(i_img, scale_min=-0.1*i_lim, scale_max=i_lim,
            non_linear=0.5)

    return img

def irstack(myfiles):
    data = []
    tempfiles = glob("templist")+glob("*b.fits")+glob("*r.fits")+glob("flat.fits")
    for myfile in tempfiles:
        os.remove(myfile)
    for myfile in myfiles:
        data.append(pyfits.getdata(myfile))
    #sky subtract and replace with median
    mediansky = mean([median(data[i]) for i in range(len(data))])
    for i in range(len(myfiles)):
        fred = pyfits.open(myfiles[i])
        im = fred[0].data
        im2 = (im.transpose() - median(im,axis=0)).transpose() + mediansky
        fred[0].data = im2
        fred.writeto("%ibb.fits"%i,'ignore',True)
    iraf.imcomb("*bb.fits","flat",combine="median",reject="sigclip",lsigma=3,hsigma=2)
    flat = pyfits.getdata('flat.fits')
    flat /= median(flat)
    fred.writeto('flat.fits','ignore',True)
    for i in range(len(myfiles)):
        fred = pyfits.open(myfiles[i])
        im = fred[0].data
        im2 = ((im/flat).transpose() - median(im,axis=0)).transpose()
        fred[0].data = im2
        fred.writeto("%irb.fits"%i,'ignore',True)
    iraf.files("*rb.fits",Stdout="templist")
    iraf.stack("templist","output",1,'none')

def reduceData(data):
    for filter in data:
        for file in data[filter]:
            fh = pyfits.open(file)
            images = []
            for ext in [1,2,3,4]:
                bias = pyfits.getdata('BIAS%i.fits' % ext, ext=0)
                flat = pyfits.getdata('FLAT_%s_%i.fits' % (filter, ext), ext=0)

                image = fh[ext].data
                hdr = fh[ext].header

                biassec = hdr['BIASSEC'].strip().replace('[', '').replace(']','').replace(',',':').split(':')

                overscan = numpy.median(image[int(biassec[2])+1:int(biassec[3])-1,
                                              int(biassec[0])+1:int(biassec[1])-1].copy().ravel())

                print 'subtracting bias of about ', numpy.mean(bias)

                if overscan > 5000:
                    img = (1.*image) - bias
                else:
                    img = (1.*image) - (bias/overscan) - overscan

                img /= flat
                images.append(img)

            fh.close()

            fh = pyfits.open(file)
            for ext in [1,2,3,4]:
                fh[ext].data = images[ext-1]

            fh.writeto('RED%s' % (file))

def combine_bins(targetSN, bins=None, outfile="a.fits"):
    """ Combine spectra of given bins. """
    data = pf.getdata("binned_sn{0}.fits".format(targetSN), 0)
    error = pf.getdata("binned_sn{0}.fits".format(targetSN), 1)
    binimg = pf.getdata("voronoi_sn{0}.fits".format(targetSN)).flatten()
    for i,bin in enumerate(bins):
        N = np.where(binimg==bin)[0].size
        spec = data[bin-1,:]
        specerr = error[bin-1,:]
        if i == 0:
            combined = N * spec
            comberr = N * specerr
            Ntot = N
        else:
            combined += N * spec
            Ntot += N * spec
            comberr += N * specerr
    h = pf.getheader("binned_sn{0}.fits".format(targetSN))
    h["NAXIS"] = 1
    del h["NAXIS2"]
    hdu0 = pf.PrimaryHDU(combined, h)
    hdu1 = pf.ImageHDU(comberr, h)
    hdulist = pf.HDUList([hdu0, hdu1])
    hdulist.writeto(outfile, clobber=True)
    return

def collapse_cube(w1, w2):
    """ Collapse a MUSE data cube.

    Arguments

    cube : MUSE data cube name containing both data and stat extensions.
    iext : Initial extension to be used. Default is one for combined cubes.

    """
    fits = "slice_w{0}_{1}.fits".format(w1, w2)
    outfits = "collapsed_w{0}_{1}.fits".format(w1, w2)
    data = pf.getdata(fits, 0)
    error = pf.getdata(fits, 1)
    h = pf.getheader(fits, 0)
    h2 = pf.getheader(fits, 1)
    h["NAXIS"] = 2
    del h["NAXIS3"]
    h2["NAXIS"] = 2
    del h2["NAXIS3"]
    print "Starting collapsing process..."
    start = time.time()
    w = wavelength_array(fits)
    # newdata = np.trapz(data, dx=np.diff(w)[0], axis=0)
    # newdata = np.nansum(data, axis=0) * np.diff(w)[0]
    newdata = np.nanmedian(data, axis=0)
    noise = 1.482602 / np.sqrt(6.) * np.nanmedian(np.abs(2.* data - \
           np.roll(data, 2, axis=0) - np.roll(data, -2, axis=0)), \
           axis=0)
    end = time.time()
    print "Collapsing lasted {0} minutes.".format((end - start)/60.)
    hdu = pf.PrimaryHDU(newdata, h)
    hdu2 = pf.ImageHDU(noise, h2)
    hdulist = pf.HDUList([hdu, hdu2])
    hdulist.writeto(outfits, clobber=True)
    return

def spt_mapping_images_to_jpeg():

    spti = pyfits.getdata("swarp.SPT-CLJ0307-5042.2x2.i.fits")
    sptv = pyfits.getdata("swarp.SPT-CLJ0307-5042.2x2.r.fits")
    sptb = pyfits.getdata("swarp.SPT-CLJ0307-5042.2x2.g.fits")

    # sptb=sptb[5066-587:5066+587,5262-587:5262+587]
    ##sptv=sptv[5441-587:5441+587,5372-587:5372+587]
    # sptv=sptv[5066-587:5066+587,5262-587:5262+587]
    # spti=spti[5066-587:5066+587,5262-587:5262+587]

    spti = spti[5262 - 587 : 5262 + 587, 5066 - 587 : 5066 + 587]
    sptv = sptv[5372 - 587 : 5372 + 587, 5441 - 587 : 5441 + 587]
    sptb = sptb[5262 - 587 : 5262 + 587, 5066 - 587 : 5066 + 587]

    skypi, sigpi = plotim(spti)
    skypv, sigpv = plotim(sptv)
    skypb, sigpb = plotim(sptb)

    print skypi, sigpi, skypv, sigpv, skypb, sigpb

    sptii = spti - skypi
    sptii = sptii / sigpi
    sptiv = sptv - skypv
    sptiv = sptiv / sigpv
    sptib = sptb - skypb
    sptib = sptib / sigpb

    acut = 0.3
    icut = 50.0 * acut
    vcut = 60.0 * acut
    bcut = 30.0 * acut

    sptii = cut_out_tails(sptii, 0, icut)
    sptiv = cut_out_tails(sptiv, 0, vcut)
    sptib = cut_out_tails(sptib, 0, bcut)

    sptii = (sptii * 255 / icut).astype(int)
    sptiv = (sptiv * 255 / vcut).astype(int)
    sptib = (sptib * 255 / bcut).astype(int)

    im = np.zeros((3, 1174, 1174), dtype=np.uint8)

    im[0, :, :] = sptii
    im[1, :, :] = sptiv
    im[2, :, :] = sptib

    # figure()
    # contour(sptii)
    # figure()
    # contour(sptiv)
    # figure()
    # contour(sptib)
    # show()

    im = np.uint8(im)

    scipy.misc.imsave("spt_comp.jpg", im)

    return 0

def twilightFlatMaker(flatImagesPath,flatDarkImagesPath,masterFlatSavePath,plots=False):
    '''
    Make a master flat using a series of images taken at twilight
    by fitting the individual pixel intensities over time using least-squares
    and use the intercept as the normalizing factor in the master flat.
    
    INPUTS: flatImagesPath - Path to the flat field exposures
    
            flatDarkImagesPath - Path to the flat field darks
            
            masterFlatSavePath - Where to save the master flat that is created
            
            plots - Plot the master flat on completion when plots=True
    '''
    ## Create zero array with the dimensions of the first image for the flat field
    [dim1, dim2] = np.shape(pyfits.getdata(flatImagesPath[0]))
    flatSum = np.zeros([dim1, dim2])

    ## Create N-dimensional array for N dark frames, where the first 
    ##    two dimensions are the dimensions of the first image
    darks = np.zeros([len(flatDarkImagesPath),dim1,dim2])

    ## Take mean of all darks
    for i in range(0,len(flatDarkImagesPath)):
        darks[i,:,:] = pyfits.getdata(flatDarkImagesPath[i])
    dark = np.mean(darks,axis=0)

    ## Create N-dimensional array for N flat frames, where the first 
    ##    two dimensions are the dimensions of the first image
    flats = np.zeros([len(flatImagesPath),dim1,dim2])

    ## Assemble data cube of flats
    for i in range(0,len(flatImagesPath)):
        flats[i,:,:] = pyfits.getdata(flatImagesPath[i]) - dark

    def linearFitIntercept(x,y):
        '''Use least-squares to find the best-fit y-intercept '''
        return np.linalg.lstsq(np.vstack([x,np.ones(len(x))]).T,y)[0][1] ## Returns intercept

    flat = np.zeros([dim1,dim2])
    for i in range(0,dim1):
        print 'Master flat computing step:',i+1,'of',dim1
        for j in range(0,dim2):
            flat[i,j] = linearFitIntercept(range(len(flats[:,i,j])),flats[:,i,j])

    masterFlat = flat/np.mean(flat)

    if plots:
        ## If plots == True, plot the resulting master flat
        fig = plt.figure()
        a = plt.imshow(masterFlat,interpolation='nearest')
        a.set_cmap('gray')
        plt.title('Normalized Master Flat Field')
        fig.colorbar(a)
        fig.canvas.set_window_title('oscaar2.0 - Master Flat') 
        plt.show()

    ## Write out both a Numpy pickle (.NPY) and a FITS file
    np.save(masterFlatSavePath+'.npy',masterFlat)
    pyfits.writeto(masterFlatSavePath+'.fits',masterFlat)

	def __init__(self, nfiles, mode='obs'):
		if mode=='obs': self.files = glob.glob("Buzzard*.fit")
		self.nfil = nfiles
		self.cat = pyfits.getdata(self.files[0])
		for i in range(1,nfiles): self.cat = np.concatenate((self.cat,pyfits.getdata(self.files[i]))) 
		self.ngal = len(self.cat)
		print 'Got %2.2f M galaxies.' %(self.ngal/1.0e6)

def shiftRGB(redF,greenF,blueF,blueshiftr=0,blueshiftc=0,greenshiftr=0,greenshiftc=0,redshiftr=0,redshiftc=0,ext=None):
    """
    this code shift the pixels of three r, g, b images. Using g image as reference and shift the other two images. It will return the shifted r,g,b images.
    each row goes along ra direction
    each col goes along dec direction
    CRVAL1 ; ra direction
    CRVAL2: dec direction
    """
    blueHdr = pf.getheader(blueF,ext)
    greenHdr = pf.getheader(greenF,ext)
    redHdr = pf.getheader(redF,ext)
    bluerow = blueHdr['crval1']*3600./0.27
    bluecol = blueHdr['crval2']*3600./0.27
    greenrow = greenHdr['crval1']*3600./0.27
    greencol = greenHdr['crval2']*3600./0.27
    redrow = redHdr['crval1']*3600./0.27
    redcol = redHdr['crval2']*3600./0.27
    """
    col0=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[0].split(':')[0])-1
    col1=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[0].split(':')[1]) 
    row0=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[1].split(':')[0])-1
    row1=int(blueHdr['datasec'].split('[')[1].split(']')[0].split(',')[1].split(':')[1]) 
    """
    blue = pf.getdata(blueF,ext)
    green = pf.getdata(greenF,ext)
    red = pf.getdata(redF,ext)

    ctgreenrow = (bluerow+greenrow+redrow)/3.
    ctgreencol = (bluecol+greencol+redcol)/3.
    blue = nd.shift(blue,[bluerow - ctgreenrow+blueshiftr,bluecol-ctgreencol+blueshiftc],mode='nearest',order=1)
    green = nd.shift(green,[greenrow - ctgreenrow+greenshiftr,greencol-ctgreencol+greenshiftc],mode='nearest',order=1)
    red = nd.shift(red,[redrow - ctgreenrow+redshiftr, redcol-ctgreencol+redshiftc],mode='nearest',order=1)
    return red,green,blue

def calc_sky_from_seg(infile,segfile):
   """
   Description: Calculates the sky level in an image using only 
      those regions in which SExtractor's segmentation file has
      a value of 0.

   Inputs:
    infile:   input fits file
    segfile:  SExtractor segmentation file associated with infile.
   """

   """ Load data """
   indat = pf.getdata(infile)
   segdat = pf.getdata(segfile)

   """ Set mask region and select associated regions of indat """
   mask = segdat == 0  
   sky = indat[mask]  
   # NB: These preceding 2 lines could have been combined as
   #   sky = indat[segdat==0]

   """ Calculate statistics """
   print "Statistics of sky outside masked regions"
   print "----------------------------------------"
   print "  N_pix  = %d" % sky.size
   print "  Median = %f" % n.median(sky)
   print "  Mean   = %f" % n.mean(sky)

   return

def do_nods(filelist):
    """
    """

    headers = [pyfits.getheader(fn) for fn in filelist]

    for ii,fn in (enumerate(filelist)):
        if ii==len(filelist)-1:
            break

        try:
            if headers[ii]['BOREOFFX'] == 0 and headers[ii+1]['BOREOFFX'] == -5.9220000000000E-03:
                fitsfile = pyfits.open(fn)
                fitsfile[0].data -= pyfits.getdata(filelist[ii+1])

                matches = difflib.SequenceMatcher(None,fn,filelist[ii+1]).get_matching_blocks()
                outfilename = fn[matches[0].a:matches[0].size] + fn[matches[0].size:matches[1].a] + "-" + filelist[ii+1][matches[0].size:matches[1].b] + fn[matches[1].a:matches[1].size+matches[1].a]
                fitsfile.writeto(outfilename)
                print matches, outfilename

            elif headers[ii+1]['BOREOFFX'] == 0 and headers[ii]['BOREOFFX'] == -5.9220000000000E-03:
                fitsfile = pyfits.open(fn)
                fitsfile[0].data = pyfits.getdata(filelist[ii+1]) - fitsfile[0].data

                matches = difflib.SequenceMatcher(None,fn,filelist[ii+1]).get_matching_blocks()
                outfilename = fn[matches[0].a:matches[0].size] + filelist[ii+1][matches[0].size:matches[1].b] + "-" + fn[matches[0].size:matches[1].a] + fn[matches[1].a:matches[1].size+matches[1].a]
                fitsfile.writeto(outfilename)
                print matches, outfilename
        except IOError:
            pass

def imWeightedAve( image1, image2, weight1, weight2, outfile, clobber=False, verbose=False):
    """
     construct a weighted average of image1 and image2:

     (weight1*image1 + weight2*image2) / (weight1+weight2)

     Mean image is written to outfile.
    """
    import os
    import pyfits
    from numpy import ndarray, nan_to_num
    import exceptions

    if os.path.isfile(outfile)  : 
        if clobber : 
            os.unlink( outfile )
        else : 
            print( "%s exists. Not clobbering."%outfile )
            return( outfile )
        
    # read in the sci and wht images 
    im1hdr = pyfits.getheader( image1 )
    im1 = pyfits.getdata( image1 )
    im2 = pyfits.getdata( image2 )
    wht1 = pyfits.getdata( weight1 )
    wht2 = pyfits.getdata( weight2 )

    meanim = nan_to_num( (wht1*im1 + wht2*im2)/(wht1+wht2) )
    
    # TODO : make a useful header
    outdir = os.path.dirname( outfile )
    if not os.path.isdir(outdir): 
        os.makedirs( outdir )
    pyfits.writeto( outfile, meanim, header=im1hdr )
    return( outfile )

def plot_image(fits,ax, rms=np.nan, F=np.nan, cont=None, contcol='r', stretch='sqrt'):
  ax.set_frame_color('k')
  ax.set_tick_color('k')
  ax.tick_labels.set_font(size='xx-small')
  ax.tick_labels.set_xformat('ddd.dd')
  ax.tick_labels.set_yformat('ddd.dd')
  
  
  head = pf.getheader(fits)
  try:
    image = pf.getdata(fits)
  except:
    print 'problem with image: ',fits
    return
  ximgsize = head.get('NAXIS1')
  yimgsize = head.get('NAXIS2')
  pixscale = head.get('CDELT1')   
   
    
  if np.isnan(rms):
    rms = np.std(image)
    av = np.median(image)
    n_one = np.ones(image.shape)
    for i in range(5):
      #print rms, av
      within1sigma = np.where((image - av*n_one) < 3.*rms*n_one)
      av = np.median(image[within1sigma])
      rms = np.std(image[within1sigma])
  if np.isnan(F):
    F = image.max()
  
  
  dat = pf.getdata(fits)
  mean = np.median(dat)
  sig = np.std(dat)
  for i in range(10):
    dat = np.ma.masked_where(abs(dat - mean) > 5.*sig,dat).compressed()
    mean = np.median(dat)
    sig = np.std(dat)
  #, vmid = mean+20*sig
  ax.show_grayscale(vmin = mean-0.5*sig, vmax = mean+15*sig, stretch=stretch, invert=True)
  

  
  if cont!=None:
    image = pf.getdata(cont)
    rms = np.std(image)
    av = np.median(image)
    n_one = np.ones(image.shape)
    for i in range(5):
      within1sigma = np.where((image - av*n_one) < 3.*rms*n_one)
      av = np.median(image[within1sigma])
      rms = np.std(image[within1sigma])
    cont_levels = 3.*rms*np.array([-2.**(1 * j ) for j in range(0, 2 )] + [ 2.**(1. * j ) for j in range(0, 10 ) ] )
    ax.show_contour(cont, hdu=0, layer='contour', levels=cont_levels, filled=False, cmap=None, colors=contcol, returnlevels=False, convention=None, slices=[0,1], smooth=None, kernel='gauss')
  #else:
    #cont_levels = 3.*rms*np.array([-2.**(1 * j ) for j in range(0, 2 )] + [ 2.**(1. * j ) for j in range(0, 10 ) ] )
    #ax.show_contour(fits, hdu=0, layer='contour', levels=cont_levels, filled=False, cmap=None, colors=contcol, returnlevels=False, convention=None, slices=[0,1], smooth=None, kernel='gauss')
    
  return

def make_color_ao_image():
    h_img = pyfits.getdata(workdir + 'mag05jullgs_h_rot.fits')
    h_imgScl = img_scale.sqrt(h_img, scale_min=50, scale_max=5500)
    kp_img = pyfits.getdata(workdir + 'mag05jullgs_kp_rot.fits')
    kp_imgScl = img_scale.sqrt(kp_img, scale_min=10, scale_max=17000)
    lp_img = pyfits.getdata(workdir + 'mag05jullgs_lp_rot.fits')
    lp_imgScl = img_scale.sqrt(lp_img, scale_min=-100, scale_max=60000)

    sgra = np.array([540, 410])
    scale = 0.00995
    h_xextent = np.array([0, h_img.shape[0]])
    h_yextent = np.array([0, h_img.shape[0]])
    h_xextent = (h_xextent - sgra[0]) * -scale
    h_yextent = (h_yextent - sgra[1]) * scale
    h_extent = [h_xextent[0], h_xextent[-1], h_yextent[0], h_yextent[-1]]

    img = np.zeros((h_img.shape[0], h_img.shape[1], 3), dtype=float)
    img[:,:,0] = lp_imgScl
    img[:,:,1] = kp_imgScl
    img[:,:,2] = h_imgScl

    py.figure(4, figsize=(10, 10))
    py.clf()
    py.subplots_adjust(left=0, right=1, bottom=0, top=1)
    py.imshow(img, extent=h_extent)
    
    ax = py.gca()
    py.setp(ax.get_xticklabels(), visible=False)
    py.setp(ax.get_yticklabels(), visible=False)

    py.xlim(5, -5)
    py.ylim(-4, 6)

    py.savefig(workdir + 'img_lgsao_color.png')

def get_1d_2d_spectra(IDtuple):
    """ Given a tuple of mask id, quadrant number, slit_id and object
    number, return the 1d extracted flux, error, the 2d image of the
    slit in which the object falls, and the y-coordinate of the objet
    in the slit.

    Note each slit can have several objects (up to 9?), each one has a
    different object number. The object y centres and edges are in
    object_?, start_?  and end_?, where '?' is the object number in
    the slit.

    Returns
    -------
    wa, fl, er, sky, image, (ystart, yobj, yend)
    """
    quad, iext, slit_id, obj = IDtuple
    image = pyfits.getdata('mos_science_flux_extracted.fits', iext)
    fluxes = pyfits.getdata('mos_science_flux_reduced.fits', iext)
    skies = pyfits.getdata('mos_sci_sky_reduced.fits', iext)
    errors = pyfits.getdata('mos_sci_error_flux_reduced.fits', iext)
    hd = pyfits.getheader('mos_science_flux_extracted.fits', iext)

    wa = hd['CRVAL1'] + np.arange(hd['NAXIS1']) * hd['CD1_1']

    i0, i1, pos, ind = get_pos_ind(IDtuple)

    fl = fluxes[ind]
    er = errors[ind]
    sky = skies[ind]

    return wa, fl, er, sky, image[i0:i1, :], pos-i0

def bfixpix(image_file, mask_file, outsuffix='_f', msksuffix='_s'):
    """
    Inputs
    ---------
    image_file : string
        input image file to fix bad pixels on

    mask_file : string
        mask file (0 == good pixels, >0 == bad pixels

    outsuffix : string
        suffix for fixed image. default = '_f'

    msksuffix : string
        suffix for bad pixels significance mask. default = '_s'
    """
    outf = image_file.replace('.fits', outsuffix + '.fits')
    outm = image_file.replace('.fits', msksuffix + '.fits')
    
    util.rmall([outf, outm])
    print("bfixpix: {0} -> {1}".format(image_file, outf))

    # fetch the image, fetch the mask
    img, hdr = pyfits.getdata(image_file, header=True)
    msk = pyfits.getdata(mask_file)

    # median the image
    medimg = ndimage.median_filter(img, 3, mode='nearest')

    # generate the pixel files
    outf_img = np.where(msk == 0, img, medimg)
    outm_img = np.where(msk == 1, (img - medimg), 0)

    pyfits.writeto(outf, outf_img, hdr)
    pyfits.writeto(outm, outm_img, hdr)