def find_sources_with_sep(img):
    """Return sources (x, y) sorted by brightness. Use SEP package.
    """
    import sep
    if isinstance(img, np.ma.MaskedArray):
        image = img.filled(fill_value=np.median(img)).astype('float32')
    else:
        image = img.astype('float32')

    bkg = sep.Background(image)
    thresh = 3. * bkg.globalrms
    try:
        sources = sep.extract(image - bkg.back(), thresh)
    except Exception as e:
        buff_message = 'internal pixel buffer full'
        if e.message[0:26] == buff_message:
            sep.set_extract_pixstack(600000)
        try:
            sources = sep.extract(image - bkg.back(), thresh)
        except Exception as e:
            if e.message[0:26] == buff_message:
                sep.set_extract_pixstack(900000)
                sources = sep.extract(image - bkg.back(), thresh)

    sources.sort(order='flux')
    return np.array([[asrc['x'], asrc['y']] for asrc in sources[::-1]])

        def compare_image(the_image):
            """Return the fraction of sources found in the original image"""
            # pixel comparison is not good, doesn't work. Compare catalogs.
            if isinstance(the_image, np.ma.MaskedArray):
                full_algn = the_image.filled(fill_value=np.median(the_image))\
                    .astype('float32')
            else:
                full_algn = the_image.astype('float32')
            full_algn[the_image == 0] = np.median(the_image)
            import sep
            bkg = sep.Background(full_algn)
            thresh = 3.0 * bkg.globalrms
            allobjs = sep.extract(full_algn - bkg.back(), thresh)
            allxy = np.array([[obj['x'], obj['y']] for obj in allobjs])

            from scipy.spatial import KDTree
            ref_coordtree = KDTree(self.star_new_pos)

            # Compare here srcs list with self.star_ref_pos
            num_sources = 0
            for asrc in allxy:
                found_source = ref_coordtree.query_ball_point(asrc, 3)
                if found_source:
                    num_sources += 1
            fraction_found = float(num_sources) / float(len(allxy))
            return fraction_found

	def __sepFindFWHM(self,tries):
		from astropy.io import fits
		import math
		import traceback
		focpos=[]
		fwhm=[]
		fwhm_min=None
		fwhm_MinimumX=None
		keys = list(tries.keys())
		keys.sort()
		ln2=math.log(2)
		for k in keys:
			try:
				fwhms=[]
				ff=fits.open(tries[k])
				# loop on images..
				for i in range(1,len(ff)):
					data=ff[i].data
					bkg=sep.Background(numpy.array(data,numpy.float))
					sources=sep.extract(data-bkg, 5.0 * bkg.globalrms)
					for s in sources:
						fwhms.append(2 * math.sqrt(ln2 * (s[12]**2 + s[13]**2)))
				im_fwhm=numpy.median(fwhms)
				# find median from fwhms measurements..
				self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,im_fwhm,len(fwhms)))
				focpos.append(k)
				fwhm.append(im_fwhm)
				if (fwhm_min is None or im_fwhm < fwhm_min):
					fwhm_MinimumX = k
					fwhm_min = im_fwhm
			except Exception as ex:
				self.log('W','offset {0}: {1} {2}'.format(k,ex,traceback.format_exc()))
		return focpos,fwhm,fwhm_min,fwhm_MinimumX

def sep_phot(data, ap, th):
    """
    Preforms photometry by SEP, similar to source extractor
    """

    # Measure a spatially variable background of some image data (np array)
    try:
        bkg = sep.Background(data)  # , mask=mask, bw=64, bh=64, fw=3, fh=3) # optional parameters
    except ValueError:
        data = data.byteswap(True).newbyteorder()
        bkg = sep.Background(data)  # , mask=mask, bw=64, bh=64, fw=3, fh=3) # optional parameters

    # Directly subtract the background from the data in place
    bkg.subfrom(data)

    # for the background subtracted data, detect objects in data given some threshold
    thresh = th * bkg.globalrms  # ensure the threshold is high enough wrt background
    objs = sep.extract(data, thresh)
    # calculate the Kron radius for each object, then we perform elliptical aperture photometry within that radius
    kronrad, krflag = sep.kron_radius(data, objs['x'], objs['y'], objs['a'], objs['b'], objs['theta'], ap)
    flux, fluxerr, flag = sep.sum_ellipse(data, objs['x'], objs['y'], objs['a'], objs['b'], objs['theta'],
                                          2.5 * kronrad, subpix=1)
    flag |= krflag  # combine flags into 'flag'

    r_min = 1.75  # minimum diameter = 3.5
    use_circle = kronrad * np.sqrt(objs['a'] * objs['b']) < r_min
    x = objs['x']
    y = objs['y']
    cflux, cfluxerr, cflag = sep.sum_circle(data, x[use_circle], y[use_circle],
                                            r_min, subpix=1)
    flux[use_circle] = cflux
    fluxerr[use_circle] = cfluxerr
    flag[use_circle] = cflag

    return objs

def detect_sources(data, threshold, elipse = False, *args, **kwargs):
    '''
    Detect the sources of the image.

    Parameters:
        data : ~numpy.ndarray~
            2D image data.
        threshold : float or ~numpy.ndarray~
            The threshold of the detection.
        elipse : bool
            Tell the program if you want the elipses parameters for each object.
        **kwargs will be passed integrally to the sep functions.

    Returns:
        x, y : ~numpy.ndarray~
            The positions of the detected sources.
        a, b, theta : ~numpy.ndarray~
            The parameters of the detected elipses.
    '''
    data = _fix_data(data)
    objs = sep.extract(data, threshold, **kwargs)

    if elipse:
        return objs['x'], objs['y'], objs['a'], objs['b'], objs['theta']
    else:
        return objs['x'], objs['y']

def extract_sources(data):
	data = data.byteswap(True).newbyteorder()
	bkg = sep.Background(data)
	back = bkg.back()
	data_woback = data-back
	thresh = 1.5 * bkg.globalrms
	objects = sep.extract(data_woback, thresh)
	return objects

    def run(self, step_name, mask=None):

        try:
            img = self.image.copy()
        except AttributeError:
            logger.error('You must setup image before running pipeline!')
            exit(1)

        logger.info('running ' + step_name)

        kws = self.step_kws[step_name].copy()

        sep_extract_kws = kws.pop('sep_extract_kws', {})
        sep_extract_kws = check_kwargs_defaults(sep_extract_kws,
                                                self.SEP_EXTRACT_DEFAULTS)

        sep_back_kws = kws.pop('sep_back_kws', {})
        sep_back_kws = check_kwargs_defaults(sep_back_kws,
                                             self.SEP_BACK_DEFAULTS)

        if mask is not None:
            logger.info('applying mask') 
            sep_extract_kws['mask'] = mask
            sep_back_kws['mask'] = mask

        # define gaussian kernel for detection
        # TODO: make filter shape optional
        num_fwhm = sep_extract_kws.pop('filter_num_fwhm', 1)
        kernel_fwhm = self.psf_fwhm * num_fwhm
        logger.info('smoothing with kernel with fwhm = {:.2f} arcsec'.\
                    format(kernel_fwhm * utils.pixscale))
        kern = Gaussian2DKernel(kernel_fwhm * gaussian_fwhm_to_sigma, 
                                mode='oversample')
        kern.normalize()
        sep_extract_kws['filter_kernel'] = kern.array

        # estimate background and subtract from image
        bkg = sep.Background(img, **sep_back_kws)
        img_sub = img - bkg

        # extract sources
        logger.info('detecting with a threshold of {} x background'.\
                    format(sep_extract_kws['thresh']))
        sources, segmap =  sep.extract(
            img_sub, err=bkg.rms(), segmentation_map=True, **sep_extract_kws)

        sources = Table(sources)
        sources = sources[sources['flux'] > 0]

        logger.info('found {} sources'.format(len(sources)))

        if kws['do_measure']:
            sources = self._measure(img, sources, mask)

        sources['seg_label'] = np.arange(1, len(sources) + 1)
        self.sources[step_name] = sources

        return sources, segmap

def createMask(data, thresh=100, title="Mask", plotMask=True):
    extraction = sep.extract(data, thresh, segmentation_map=True)
    mask = extraction[1]
    if plotMask:
        plt.figure()
        plt.imshow(mask, aspect="auto", interpolation="nearest", origin="lower")
        plt.title(title)
        plt.show()
    return mask

def test_long_error_msg():
    """Ensure that the error message is created successfully when
    there is an error detail."""

    # set extract pixstack to an insanely small value; this will trigger
    # a detailed error message when running sep.extract()
    old = sep.get_extract_pixstack()
    sep.set_extract_pixstack(5)

    data = np.ones((10, 10), dtype=np.float64)
    with pytest.raises(Exception) as excinfo:
        sep.extract(data, 0.1)
    msg = excinfo.value.args[0]
    assert type(msg) == str  # check that message is the native string type
    assert msg.startswith("internal pixel buffer full: The limit")

    # restore
    sep.set_extract_pixstack(old)

        def auto_fit(image_data):
            an_width = 3
            temp_image_data = image_data
            internal_image_data = temp_image_data.byteswap(True).newbyteorder()
            bkg = sep.Background(internal_image_data)
            thresh = 1.5 * bkg.globalback
            objects = sep.extract(internal_image_data, thresh)
            center_x = internal_image_data.shape[0]/2.0
            center_y = internal_image_data.shape[1]/2.0
            center = [center_x, center_y]
            smallest_dFoc = 1000000000
            radii = 21

            count = 0
            for i, j in zip(objects['x'], objects['y']):
                pos = [i, j]

                dFoc = math.sqrt(((pos[0]-center[0])**2) + ((pos[1]-center[1])**2))
                if abs(dFoc) < smallest_dFoc:
                    smallest_dFoc = dFoc
                    minx = objects['xmin'][count]
                    miny = objects['ymin'][count]
                    maxx = objects['xmax'][count]
                    maxy = objects['ymax'][count]
                    radii = abs(math.sqrt(((maxx-minx)**2)+((maxy-miny)**2)))
                else:
                    pass
                count += 1
            i = 0
            while i < 25:
                theta = 0
                area = 0
                while theta <= 2*pi:
                    area += (((i+.1)**2)/2) - ((i**2)/2)
                    theta += .001
                flux, fluxerr, flag = sep.sum_circann(internal_image_data, internal_image_data.shape[0]/2, internal_image_data.shape[1]/2, i, i+an_width)
                metric = (flux - bkg.globalback)/bkg.globalrms
                metric /= area
                if i > 1:
                    if metric < smallest_metric:
                        smallest_metric = metric
                        annulus = [i, i+an_width]
                else:
                    smallest_metric = metric
                    annulus = [i, i+an_width]
                i += 1
            i = 1
            an = [math.ceil(x*6) for x in annulus]
            if radii > an[0]:
                radii = an[0] - 0.5
            if radii > 30:
                radii = 30
            else:
                pass
            return {'Ap': radii, 'InAn': an[0], 'OutAn': an[1]}

def find_objects(obs, segmentation_map=False):
    import sep

    noise=np.sqrt(1.0/obs.weight[0,0])
    return sep.extract(
        obs.image,
        SEP_THRESH,
        segmentation_map=segmentation_map,
        err=noise,
        **SEP_PARS
    )

def find_stars_on_data(data, verbose = 0, useDS9 = False):
	bkg = sep.Background(data)
	bkg.subfrom(data)
	thres = 1.5 * bkg.globalrms
	if verbose > 1:
		print('global average background: {0:.2f} rms: {1:.3f} threshold: {2:.3f}'.format(bkg.globalback, bkg.globalrms, thres))
	objects = sep.extract(data, thres)
	# order by flux
	if len(objects) == 0:
		return []
	return sorted(objects, cmp=lambda x,y: cmp(y['flux'],x['flux']))

def extract_sources(data):
	try:
		bkg = sep.Background(data)
	except ValueError:
		data = data.byteswap().newbyteorder()
		bkg = sep.Background(data)
	back = bkg.back()
	data_woback = data-back
	thresh = 3.0 * bkg.globalrms
	objects = sep.extract(data_woback, thresh)
	return objects

def test_extract_segmentation_map():

    # Get some background-subtracted test data:
    data = np.copy(image_data)
    bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
    bkg.subfrom(data)

    objects, segmap = sep.extract(data, 1.5*bkg.globalrms,
                                  segmentation_map=True)

    assert type(segmap) is np.ndarray
    assert segmap.shape == data.shape
    for i in range(len(objects)):
        assert objects["npix"][i] == (segmap == i+1).sum()

    def ImageObjDetection(self,filter):
        print 'SEP object detection on %s......'%filter
        sigma = 3.0
        threshold = sigma * self.bkgRMS[filter]
        objs = sep.extract(self.dataList[filter], threshold, minarea=20)
        # add fields for mag and magerr
        desc = np.dtype([('RA','float64'),('DEC','float64')])
        newObjs = np.zeros(objs.shape, dtype = objs.dtype.descr + desc.descr)
        for name in objs.dtype.names:
            newObjs[name] = objs[name]
        # add RA and DEC to objs
        newObjs['RA'], newObjs['DEC'] = self.images[filter].xy_to_rd(objs['x'],objs['y'])

        return newObjs

def find_centroid(data, t):
  '''
  filter_kernel = makeGaussian(17, 5. , 0 , np.array([8.5,8.5]))
  source = extract(data , t, filter_kernel=filter_kernel)
  '''
  #t=20
  #source = extract(data , t)
  #source = extract(data, 100)
  #source = extract(data, t)
  #source = extract(data, 500)
  source = extract(data, t)

  a = source['a']
  b = source['b']
  #print 'a: {0}'.format(a)
  #print 'b: {0}'.format(b)
  flux = source['cflux']
  #arg = np.argsort(flux)[-1]
  try:
    arg = np.argsort(flux)[-1]
  except IndexError:
    return None
  #print flux

  x = source['x'][arg]
  y = source['y'][arg]
  try:
    #fwhm = np.sqrt(np.max(a)*np.max(b))
    fwhm = np.sqrt(a[arg]*b[arg])
    print 'fwhm:{0}'.format(fwhm)
  except ValueError:
    return None
  size = data.shape[0]
  zero = size/2 + .5
  kernel = makeGaussian(17, fwhm , 0 , np.array([8.5,8.5]))
  img = signal.convolve2d(data , kernel , mode = "same")
  max_value = np.max(img)
  xi, yi = np.unravel_index(np.argmax(img), img.shape)
  
  if (xi >= 1 and xi < img.shape[0] - 1 and yi >= 1 and yi < img.shape[1] - 1):
      ox, oy = fit_3x3(img[xi-1:xi+2, yi-1:yi+2])
  else:
      ox , oy = 0. , 0.
  if (np.absolute(ox) >3) or (np.absolute(oy)>3):
      ox, oy = 0., 0.
  #return xi + ox + .5 , yi + oy + .5, max_value, flux[arg[-1]]
  print xi,yi, ox, oy
  return xi + ox + .5 , yi + oy + .5, max_value, flux[arg]

def find_stars(fn, hdu, verbose = 0, useDS9 = False, cube = None):
	"""Find stars on the image. Returns flux ordered list of stars."""
	if cube is None:
		data = np.array(hdu[0].data,np.int32)
	else:
		data = np.array(hdu[0].data[cube],np.int32)
	bkg = sep.Background(data)
	bkg.subfrom(data)
	thres = 1.5 * bkg.globalrms
	if verbose > 1:
		print 'global average background: {0:.2f} rms: {1:.3f} threshold: {2:.3f}'.format(bkg.globalback, bkg.globalrms, thres)
	objects = sep.extract(data, thres)
	# order by flux
	if len(objects) == 0:
		return None
	return sorted(objects, cmp=lambda x,y: cmp(y['flux'],x['flux']))

def test_extract_matched_filter_at_edge():
    """Exercise bug where bright star at end of image not detected
    with noise array and matched filter on."""

    data = np.zeros((20, 20))
    err = np.ones_like(data)
    kernel = np.array([[1., 2., 1.],
                       [2., 4., 2.],
                       [1., 2., 1.]])

    data[18:20, 9:12] = kernel[0:2, :]

    objects, pix = sep.extract(data, 2.0, err=err, filter_kernel=kernel,
                               filter_type="matched", segmentation_map=True)
    assert len(objects) == 1
    assert objects["npix"][0] == 6

def extract(data):
	bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
	bkg.subfrom(data)
	objs = sep.extract(data, 1.5*bkg.globalrms)
	flux, fluxerr, flag = sep.sum_circle(data, objs['x'], objs['y'], 5.,
                                         err=bkg.globalrms)
	kr, flag = sep.kron_radius(data, objs['x'], objs['y'], objs['a'],
                               			objs['b'], objs['theta'], 6.0)
	eflux, efluxerr, eflag = sep.sum_ellipse(data, objs['x'], objs['y'],
                                          objs['a'], objs['b'],
                                          objs['theta'], r=2.5 * kr,
                                          err=bkg.globalrms, subpix=1)
	retstr = ""
	for i in range(len(objs['x'])):
		retstr = retstr+(str(objs['x'][i])+"\t"+str(objs['y'][i])+"\t"+str(flux[i])+"\t"+str(fluxerr[i])+"\t"+str(kr[i])+"\t"+str(eflux[i])+"\t"+str(efluxerr[i])+"\t"+str(flag[i])+"\n")
	return retstr

def getobjects( data, thresh=5.0 ):
	if data.dtype is not 'float32':
		data = np.array(data, dtype="float32")

	bkg = sep.Background(data)
	bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)

	back = bkg.back()

	rms = bkg.rms()
	bkg.subfrom(data)
	thresh = thresh * bkg.globalrms
	objects = sep.extract(data, thresh)

	#print "there are ", len(objects), "objects"
	return objects