Peter Müller / NOD3 single dish reduction software

# -*- coding: utf-8 -*- """
#!/usr/bin/env python

import numpy as np

def nint(x):
    if x > 0: return int(x+0.5000001)
    else: return int(x-0.4999999)

def nanmean(a, axis=None):
    if hasattr(np, 'nanmedian'):
       return np.nanmean(a, axis)
    else:
       from scipy import stats
       return stats.nanmean(np.array(a), axis)

def nanmedian(a, axis=None):
    if hasattr(np, 'nanmedian'):
       return np.nanmedian(a, axis)
    else:
       from scipy import stats
       return stats.nanmedian(np.array(a), axis)

def nan_check(data, val, weight=False):
    Data = 1*data
    mask = np.isfinite(Data)
    Data[~mask] = val
    if weight:
       return Data, np.array(mask, dtype=float)
    else:
       return Data

def nanmean_(datalist):
    wsum = None
    for d in datalist:
        data, w = nan_check(d, 0.0, weight=True)
        if wsum == None:
           wsum = w.copy()
           wdata = data.copy()
        else:
           wsum += w
           wdata += data
    mask = np.isfinite(wdata)
    avg = np.where(mask, wdata/wsum, np.nan)
    return avg

def erf(x):
    a0 = -1.26551223
    a1 = +1.00002368
    a2 = +0.37409196
    a3 = +0.09678418
    a4 = -0.18628806
    a5 = +0.27886807
    a6 = -1.13520398
    a7 = +1.48851578
    a8 = -0.82215223
    a9 = +0.17087227
    t = 1.0/(1.0+0.5*abs(x))

    p = ((((((((a9*t + a8)*t + a7)*t + a6)*t + a5*t) + a4)*t + a3)*t + a2)*t + a1)*t + a0
    tau = t*np.exp(-x*x + p)

    erfx = np.where(x >=0.0, 1-tau, tau-1)
    
    return erfx

def erf0(x):
    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # Save the sign of x
    sign = np.where(x < 0.0, -1, 1)
    x = abs(x)

    # A & S 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*np.exp(-x*x)

    return sign*y

def erfc(x):
    return 1.0 - erf(x)

def padsin(vector, pad_width, iaxis, kwargs):
    px = pad_width[0]
    py = pad_width[1]
    a = vector[px]
    b = vector[-py-1]
    vector[:px] = a*(np.sin(np.pi/(px-0)*np.arange(float(px))-np.pi/2.) + 1)/2.
    vector[-py:] = b-b*(np.sin(np.pi/(py-0)*np.arange(float(py))-np.pi/2.) + 1)/2.
    return vector

def deg2DMS(ra='', dec='', round=False):
    RA, DEC, rs, ds = '', '', '', ''
    if dec:
       if str(dec)[0] == '-':
          ds, dec = '-', abs(dec)
       deg = int(dec)
       decM = abs(int((dec-deg)*60))
       if round:
          decS = int((abs((dec-deg)*60)-decM)*60)
       else:
          decS = (abs((dec-deg)*60)-decM)*60
       DEC = '{0}{1} {2} {3}'.format(ds, deg, decM, decS)
  
    if ra:
       if str(ra)[0] == '-':
          rs, ra = '-', abs(ra)
       raH = int(ra)
       raM = int(((ra)-raH)*60)
       if round:
          raS = int(((((ra)-raH)*60)-raM)*60)
       else:
          raS = ((((ra)-raH)*60)-raM)*60
       RA = '{0}{1} {2} {3}'.format(rs, raH, raM, raS)
  
    if ra and dec:
       return (RA, DEC)
    else:
       return RA or DEC


def deg2HMS(ra='', dec='', round=False):
    RA, DEC, rs, ds = '', '', '', ''
    if dec:
       if str(dec)[0] == '-':
          ds, dec = '-', abs(dec)
       deg = int(dec)
       decM = abs(int((dec-deg)*60))
       if round:
          decS = int((abs((dec-deg)*60)-decM)*60)
       else:
          decS = (abs((dec-deg)*60)-decM)*60
       DEC = '{0}{1} {2} {3}'.format(ds, deg, decM, decS)
  
    if ra:
       if str(ra)[0] == '-':
          rs, ra = '-', abs(ra)
       raH = int(ra/15)
       raM = int(((ra/15)-raH)*60)
       if round:
          raS = int(((((ra/15)-raH)*60)-raM)*60)
       else:
          raS = ((((ra/15)-raH)*60)-raM)*60
       RA = '{0}{1} {2} {3}'.format(rs, raH, raM, raS)
  
    if ra and dec:
       return (RA, DEC)
    else:
       return RA or DEC


def extract(Z, shape, position, fill=np.NaN):
     """ Extract a sub-array from Z using given shape and centered on  
position.
         If some part of the sub-array is out of Z bounds, result is  
padded
         with fill value.

         **Parameters**
             `Z` : array_like
                Input array.

            `shape` : tuple
                Shape of the output array

            `position` : tuple
                Position within Z

            `fill` : scalar
                Fill value

         **Returns**
             `out` : array_like
                 Z slice with given shape and center

         **Examples**

         >>> Z = np.arange(0,16).reshape((4,4))
         >>> extract(Z, shape=(3,3), position=(0,0))
         [[ NaN  NaN  NaN]
          [ NaN   0.   1.]
          [ NaN   4.   5.]]

         Schema:

             +-----------+
             | 0   0   0 | = extract (Z, shape=(3,3), position=(0,0))
             |   +---------------+
             | 0 | 0   1 | 2   3 | = Z
             |   |       |       |
             | 0 | 4   5 | 6   7 |
             +---|-------+       |
                 | 8   9  10  11 |
                 |               |
                 | 12 13  14  15 |
                 +---------------+

         >>> Z = np.arange(0,16).reshape((4,4))
         >>> extract(Z, shape=(3,3), position=(3,3))
         [[ 10.  11.  NaN]
          [ 14.  15.  NaN]
          [ NaN  NaN  NaN]]

         Schema:

             +---------------+
             | 0   1   2   3 | = Z
             |               |
             | 4   5   6   7 |
             |       +-----------+
             | 8   9 |10  11 | 0 | = extract (Z, shape=(3,3),  
position=(3,3))
             |       |       |   |
             | 12 13 |14  15 | 0 |
             +---------------+   |
                     | 0   0   0 |
                     +-----------+
     """
#    assert(len(position) == len(Z.shape))
#    if len(shape) < len(Z.shape):
#        shape = shape + Z.shape[len(Z.shape)-len(shape):]

     R = np.ones(shape, dtype=Z.dtype)*fill
     P  = np.array(list(position)).astype(int)
     Rs = np.array(list(R.shape)).astype(int)
     Zs = np.array(list(Z.shape)).astype(int)

     R_start = np.zeros((len(shape),)).astype(int)
     R_stop  = np.array(list(shape)).astype(int)
     Z_start = (P-Rs//2)
     Z_stop  = (P+Rs//2)+Rs%2

     R_start = (R_start - np.minimum(Z_start,0)).tolist()
     Z_start = (np.maximum(Z_start,0)).tolist()
     R_stop = (R_stop - np.maximum(Z_stop-Zs,0)).tolist()
     Z_stop = (np.minimum(Z_stop,Zs)).tolist()

     r = [slice(start,stop) for start,stop in zip(R_start,R_stop)]
     z = [slice(start,stop) for start,stop in zip(Z_start,Z_stop)]

     R[r] = Z[z]

     return R

def same_size(ms, ref=False, eps=1.e-2):
    cx = ms[0].header['CRVAL1']
    cy = ms[0].header['CRVAL2']
    dx = ms[0].header['CDELT1']
    dy = ms[0].header['CDELT2']
    tx = ms[0].header['CTYPE1']
    ty = ms[0].header['CTYPE2']
    Rows = 0
    Cols = 0
    medge = None
    edges = []
    for m in ms:
        rows, cols = m.data.shape
        Rows = max(Rows, rows)
        Cols = max(Cols, cols)
        xl, xr = cx + (1-m.header['CRPIX1'])*dx, cx + (cols-m.header['CRPIX1'])*dx
        yb, yt = cy + (1-m.header['CRPIX2'])*dy, cy + (rows-m.header['CRPIX2'])*dy
        edges.append([yb, yt, xl, xr])
        if medge == None:
           medge = [yb, yt, xl, xr] 
           #b0 = cy
           #l0 = cx
           b0 = (yb+yt)/2.0
           l0 = (xl+xr)/2.0
        elif not ref:
           medge[0] = min(medge[0], yb)
           medge[1] = max(medge[1], yt)
           medge[2] = max(medge[2], xl)
           medge[3] = min(medge[3], xr)
           b0 = (medge[1]+medge[0])/2.0
           l0 = (medge[3]+medge[2])/2.0
    Rows = nint((medge[1]-medge[0]) / dy) + 1
    Cols = nint((medge[3]-medge[2]) / dx) + 1

    out = 0
    for m in ms:
        if (m.header['CDELT1'] - dx)/dx > eps and (m.header['CDELT2'] - dy)/dy < eps:
           out += 1
        if m.header['CTYPE1'] != tx and m.header['CTYPE2'] != ty:
           print tx, ty, m.header['CTYPE1'], m.header['CTYPE2']
           out += 1
    if out > 0:
       return out, []

    ms_new = []
    n = 0
    for m in ms:
        rows, cols = m.data.shape
        b = (edges[n][1]+edges[n][0])/2.0
        l = (edges[n][3]+edges[n][2])/2.0
        posx = cols/2.0 + (l0-l)/dx
        posy = rows/2.0 + (b0-b)/dy
        #print (Rows, Cols), (posx,posy)
        m.data = extract(m.data, shape=(Rows, Cols), position=(nint(posy), nint(posx)))
        m.header['CRPIX1'] = (Cols+1)/2.0 + (cx-l0)/dx
        m.header['CRPIX2'] =  (Rows+1)/2.0 + (cy-b0)/dy
        #m.header['CRPIX1'] = (Cols+1)/2.0
        #m.header['CRPIX2'] = (Rows+1)/2.0 
        #m.header['CRVAL1'] = l0
        #m.header['CRVAL2'] = b0
        m.header['CRVAL1'] = cx
        m.header['CRVAL2'] = cy
        m.header['NAXIS1'] = Cols
        m.header['NAXIS2'] = Rows
        ms_new.append(m)
        n += 1
    return out, ms_new
        
def FixNaNs(x):
    for row in range(x.shape[0]):
        xrow = x[row]
        idxs = np.nonzero(xrow == xrow)[0]
        if len(idxs) > 0:
           for i in range(0, idxs[0]):
               xrow[i] = xrow[idxs[0]]
           for i in range(idxs[-1]+1, xrow.size):
               xrow[i] = xrow[idxs[-1]]
           x[row] = xrow
    return x


def nan_interpol(data):
    mask = np.isnan(data)
    mdata = data.copy()
    mdata[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask])
    return mask, mdata

from scipy.ndimage import rotate
def map_rotate(data, p):
    mask = np.isnan(data)
    dmask = rotate(np.where(mask, 0, 1), p.angle, reshape=p.reshape, order=0, mode='nearest',
                   cval=0, prefilter=False)
    data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask])
    newdata = rotate(data, p.angle, reshape=p.reshape, order=p.order, mode=p.mode,
                     cval=np.nan, prefilter=p.prefilter)
    return np.where(dmask == 0, np.nan, newdata)

from scipy.ndimage import map_coordinates
def map_interpolate(data, x, y):
    mask = np.isnan(data)
    dmask = map_coordinates(np.where(mask, 0, 1), (y, x), order=0, mode='nearest',
                            cval=0, prefilter=False, output=np.int32)
    data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask])
    newdata = map_coordinates(data, (y, x), order=4, mode='constant',
                              cval=np.nan, prefilter=True, output=np.float32)
    return np.where(dmask == 0, np.nan, newdata)
    return newdata

from scipy.ndimage import zoom
def map_zoom(data, rebin, order=4, prefilter=True):
    mask = np.isnan(data)
    dmask = zoom(np.where(mask, 0, 1), rebin, order=0, prefilter=False, output=np.int32)
    data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask])
    newdata = zoom(data, rebin, order=order, cval=np.nan, prefilter=True, output=np.float32)
    return np.where(dmask == 0, np.nan, newdata)


def griddata(inp, points, outp, method='nearest'):
    y_points, x_points = inp
    grid_y, grid_x = outp
    shape = grid_x.shape
    grid_y = grid_y.ravel()
    grid_x = grid_x.ravel()
    store = np.zeros(len(grid_x))
    for i in range(len(grid_x)):
        distance2 = (y_points-grid_y[i])**2 + (x_points-grid_x[i])**2
        j = np.where(distance2 == np.nanmin(distance2))
        store[i] = points[j]
    return store.reshape(shape)


def nan_interpolation(data):
    mask = np.isnan(data)
    data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask]) 
    return data

def _polyfill(polygon):
    return lambda x, y: point_in_polygon((x, y), polygon)

def polymask(data, polygon):
    py, px = np.indices(data.shape)
    mask = _polyfill(polygon)(*np.array([px, py]))
    return mask.reshape(data.shape)

def sector_mask(shape,centre,radius,angle_range):
    """
    Return a boolean mask for a circular sector. The start/stop angles in  
    `angle_range` should be given in clockwise order.
    """

    x,y = np.ogrid[:shape[0],:shape[1]]
    cx,cy = centre
    tmin,tmax = np.deg2rad(angle_range)

    # ensure stop angle > start angle
    if tmax < tmin:
       tmax += 2*np.pi

    # convert cartesian --> polar coordinates
    r2 = (x-cx)*(x-cx) + (y-cy)*(y-cy)
    theta = np.arctan2(x-cx, y-cy) - tmin

    # wrap angles between 0 and 2*pi
    theta %= (2*np.pi)

    # circular mask
    rad1, rad2 = radius
    circmask1 = r2 >= rad1*rad2
    circmask2 = r2 <= rad2*rad2
    circmask = circmask1*circmask2

    # angular mask
    anglemask = theta <= (tmax-tmin)

    return circmask*anglemask

#From http://www.ariel.com.au/a/python-point-int-poly.html
# Modified by Nick ODell
from collections import namedtuple
def point_in_polygon(target, poly):
    """x,y is the point to test. poly is a list of tuples comprising the polygon."""
    point = namedtuple("Point", ("x", "y"))
    line = namedtuple("Line", ("p1", "p2"))
    target = point(*target)

    inside = False
    # Build list of coordinate pairs
    # First, turn it into named tuples

    poly = map(lambda p: point(*p), poly)

    # Make two lists, with list2 shifted forward by one and wrapped around
    list1 = poly
    list2 = poly[1:] + [poly[0]]
    poly = map(line, list1, list2)

    for l in poly:
        p1 = l.p1
        p2 = l.p2

        if p1.y == p2.y:
            # This line is horizontal and thus not relevant.
            continue
        if max(p1.y, p2.y) < target.y <= min(p1.y, p2.y):
            # This line is too high or low
            continue
        if target.x < max(p1.x, p2.x):
            # Ignore this line because it's to the right of our point
            continue
        # Now, the line still might be to the right of our target point, but 
        # still to the right of one of the line endpoints.
        rise = p1.y - p2.y
        run =  p1.x - p2.x
        try:
            slope = rise/float(run)
        except ZeroDivisionError:
            slope = float('inf')

        # Find the x-intercept, that is, the place where the line we are
        # testing equals the y value of our target point.

        # Pick one of the line points, and figure out what the run between it
        # and the target point is.
        run_to_intercept = target.x - p1.x
        x_intercept = p1.x + run_to_intercept / slope
        if target.x < x_intercept:
            # We almost crossed the line.
            continue

        inside = not inside

    return inside

def toworld(header, pix):
    """ converts pixel to world coordnates """
    posx = header['CRVAL1'] + (pix[0] - header['CRPIX1']) * header['CDELT1']
    posy = header['CRVAL2'] + (pix[1] - header['CRPIX2']) * header['CDELT2']
    return (posx, posy)

def topixel(header, pos):
    """ converts world coordinates ti pixel """
    pix = (pos[0] - header['CRVAL1']) / header['CDELT1'] + header['CRPIX1']
    piy = (pos[1] - header['CRVAL2']) / header['CDELT2'] + header['CRPIX2']
    return (pix, piy)

def toworldn(header, pix):
    """ converts pixel to world coordnates """
    p = np.array(pix)
    n,m,k = p.shape
    p = p.reshape((n*k,m)).T
    p[0] = header['CRVAL1'] + (p[0] - header['CRPIX1']) * header['CDELT1']
    p[1] = header['CRVAL2'] + (p[1] - header['CRPIX2']) * header['CDELT2']
    p = p.T
    pos = p.reshape((n,m,k))
    return pos

def topixeln(header, pos):
    """ converts world coordinates ti pixel """
    p = np.array(pos)
    n,m,k = p.shape
    p = p.reshape((n*k,m)).T
    p[0] = (p[0] - header['CRVAL1']) / header['CDELT1'] + header['CRPIX1']
    p[1] = (p[1] - header['CRVAL2']) / header['CDELT2'] + header['CRPIX2']
    p = p.T
    pix = p.reshape((n,m,k))
    return pix

def smooth(x, window_len=11, window='hanning'):
    """smooth the data using a window with requested size.
    
    This method is based on the convolution of a scaled window with the signal.
    The signal is prepared by introducing reflected copies of the signal 
    (with the window size) in both ends so that transient parts are minimized
    in the begining and end part of the output signal.
    
    input:
        x: the input signal 
        window_len: the dimension of the smoothing window; should be an odd integer
        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
            flat window will produce a moving average smoothing.

    output:
        the smoothed signal
        
    example:

    t=linspace(-2,2,0.1)
    x=sin(t)+randn(len(t))*0.1
    y=smooth(x)
    
    see also: 
    
    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
    scipy.signal.lfilter
 
    TODO: the window parameter could be the window itself if an array instead of a string
    NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.
    """
    # nan interpolation
    x = nan_interpolation(x)
    # window_len must be odd
    if 1-window_len%2: window_len += 1
    if x.ndim != 1:
       raise ValueError, "smooth only accepts 1 dimension arrays."
    if x.size < window_len:
       raise ValueError, "Input vector needs to be bigger than window size."
    if window_len < 3:
       return x
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
       raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"

    s=np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
    if window == 'flat': #moving average
       w = np.ones(window_len, 'd')
    else:
       w = eval('np.'+window+'(window_len)')
    y = np.convolve(w/w.sum(), s, mode='valid')
    return y[(window_len/2):-(window_len/2)]
    ly = (len(y)-len(x))/2 
    return y[ly:-ly]

"""
Port of Manuel Guizar's code from:
http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation
"""

import numpy as np


def _upsampled_dft(data, upsampled_region_size,
                   upsample_factor=1, axis_offsets=None):
    """
    Upsampled DFT by matrix multiplication.
    This code is intended to provide the same result as if the following
    operations were performed:
        - Embed the array "data" in an array that is ``upsample_factor`` times
          larger in each dimension.  ifftshift to bring the center of the
          image to (1,1).
        - Take the FFT of the larger array.
        - Extract an ``[upsampled_region_size]`` region of the result, starting
          with the ``[axis_offsets+1]`` element.
    It achieves this result by computing the DFT in the output array without
    the need to zeropad. Much faster and memory efficient than the zero-padded
    FFT approach if ``upsampled_region_size`` is much smaller than
    ``data.size * upsample_factor``.
    Parameters
    ----------
    data : 2D ndarray
        The input data array (DFT of original data) to upsample.
    upsampled_region_size : integer or tuple of integers, optional
        The size of the region to be sampled.  If one integer is provided, it
        is duplicated up to the dimensionality of ``data``.
    upsample_factor : integer, optional
        The upsampling factor.  Defaults to 1.
    axis_offsets : tuple of integers, optional
        The offsets of the region to be sampled.  Defaults to None (uses
        image center)
    Returns
    -------
    output : 2D ndarray
            The upsampled DFT of the specified region.
    """
    # if people pass in an integer, expand it to a list of equal-sized sections
    if not hasattr(upsampled_region_size, "__iter__"):
        upsampled_region_size = [upsampled_region_size, ] * data.ndim
    else:
        if len(upsampled_region_size) != data.ndim:
            raise ValueError("shape of upsampled region sizes must be equal "
                             "to input data's number of dimensions.")

    if axis_offsets is None:
        axis_offsets = [0, ] * data.ndim
    else:
        if len(axis_offsets) != data.ndim:
            raise ValueError("number of axis offsets must be equal to input "
                             "data's number of dimensions.")

    col_kernel = np.exp(
        (-1j * 2 * np.pi / (data.shape[1] * upsample_factor)) *
        (np.fft.ifftshift(np.arange(data.shape[1]))[:, None] -
         np.floor(data.shape[1] / 2)).dot(
             np.arange(upsampled_region_size[1])[None, :] - axis_offsets[1])
    )
    row_kernel = np.exp(
        (-1j * 2 * np.pi / (data.shape[0] * upsample_factor)) *
        (np.arange(upsampled_region_size[0])[:, None] - axis_offsets[0]).dot(
            np.fft.ifftshift(np.arange(data.shape[0]))[None, :] -
            np.floor(data.shape[0] / 2))
    )

    return row_kernel.dot(data).dot(col_kernel)


def _compute_phasediff(cross_correlation_max):
    """
    Compute global phase difference between the two images (should be
        zero if images are non-negative).
    Parameters
    ----------
    cross_correlation_max : complex
        The complex value of the cross correlation at its maximum point.
    """
    return np.arctan2(cross_correlation_max.imag, cross_correlation_max.real)


def _compute_error(cross_correlation_max, src_amp, target_amp):
    """
    Compute RMS error metric between ``src_image`` and ``target_image``.
    Parameters
    ----------
    cross_correlation_max : complex
        The complex value of the cross correlation at its maximum point.
    src_amp : float
        The normalized average image intensity of the source image
    target_amp : float
        The normalized average image intensity of the target image
    """
    error = 1.0 - cross_correlation_max * cross_correlation_max.conj() /\
        (src_amp * target_amp)
    return np.sqrt(np.abs(error))


#[1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
#    "Efficient subpixel image registration algorithms," Optics Letters 33,
#    156-158 (2008).
#======================================
# Cross-Correlation (Phase Correlation)
#======================================

def register_translation(src_image, target_image, upsample_factor=1,
                         space="real"):
    """
    Efficient subpixel image translation registration by cross-correlation.
    This code gives the same precision as the FFT upsampled cross-correlation
    in a fraction of the computation time and with reduced memory requirements.
    It obtains an initial estimate of the cross-correlation peak by an FFT and
    then refines the shift estimation by upsampling the DFT only in a small
    neighborhood of that estimate by means of a matrix-multiply DFT.
    Parameters
    ----------
    src_image : ndarray
        Reference image.
    target_image : ndarray
        Image to register.  Must be same dimensionality as ``src_image``.
    upsample_factor : int, optional
        Upsampling factor. Images will be registered to within
        ``1 / upsample_factor`` of a pixel. For example
        ``upsample_factor == 20`` means the images will be registered
        within 1/20th of a pixel.  Default is 1 (no upsampling)
    space : string, one of "real" or "fourier"
        Defines how the algorithm interprets input data.  "real" means data
        will be FFT'd to compute the correlation, while "fourier" data will
        bypass FFT of input data.  Case insensitive.
    Returns
    -------
    shifts : ndarray
        Shift vector (in pixels) required to register ``target_image`` with
        ``src_image``.  Axis ordering is consistent with numpy (e.g. Z, Y, X)
    error : float
        Translation invariant normalized RMS error between ``src_image`` and
        ``target_image``.
    phasediff : float
        Global phase difference between the two images (should be
        zero if images are non-negative).
    References
    ----------
    .. [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
           "Efficient subpixel image registration algorithms,"
           Optics Letters 33, 156-158 (2008).
    """
    # images must be the same shape
    if src_image.shape != target_image.shape:
        raise ValueError("Error: images must be same size for "
                         "register_translation")

    # only 2D data makes sense right now
    if src_image.ndim != 2 and upsample_factor > 1:
        raise NotImplementedError("Error: register_translation only supports "
                                  "subpixel registration for 2D images")

    # assume complex data is already in Fourier space
    if space.lower() == 'fourier':
        src_freq = src_image
        target_freq = target_image
    # real data needs to be fft'd.
    elif space.lower() == 'real':
        src_image = np.array(src_image, dtype=np.complex128, copy=False)
        target_image = np.array(target_image, dtype=np.complex128, copy=False)
        src_freq = np.fft.fftn(src_image)
        target_freq = np.fft.fftn(target_image)
    else:
        raise ValueError("Error: register_translation only knows the \"real\" "
                         "and \"fourier\" values for the ``space`` argument.")

    # Whole-pixel shift - Compute cross-correlation by an IFFT
    shape = src_freq.shape
    image_product = src_freq * target_freq.conj()
    cross_correlation = np.fft.ifftn(image_product)

    # Locate maximum
    maxima = np.unravel_index(np.argmax(np.abs(cross_correlation)),
                              cross_correlation.shape)
    midpoints = np.array([np.fix(axis_size / 2) for axis_size in shape])

    shifts = np.array(maxima, dtype=np.float64)
    shifts[shifts > midpoints] -= np.array(shape)[shifts > midpoints]

    if upsample_factor == 1:
        src_amp = np.sum(np.abs(src_freq) ** 2) / src_freq.size
        target_amp = np.sum(np.abs(target_freq) ** 2) / target_freq.size
        CCmax = cross_correlation.max()
    # If upsampling > 1, then refine estimate with matrix multiply DFT
    else:
        # Initial shift estimate in upsampled grid
        shifts = np.round(shifts * upsample_factor) / upsample_factor
        upsampled_region_size = np.ceil(upsample_factor * 1.5)
        # Center of output array at dftshift + 1
        dftshift = np.fix(upsampled_region_size / 2.0)
        upsample_factor = np.array(upsample_factor, dtype=np.float64)
        normalization = (src_freq.size * upsample_factor ** 2)
        # Matrix multiply DFT around the current shift estimate
        sample_region_offset = dftshift - shifts*upsample_factor
        cross_correlation = _upsampled_dft(image_product.conj(),
                                           upsampled_region_size,
                                           upsample_factor,
                                           sample_region_offset).conj()
        cross_correlation /= normalization
        # Locate maximum and map back to original pixel grid
        maxima = np.array(np.unravel_index(
                              np.argmax(np.abs(cross_correlation)),
                              cross_correlation.shape),
                          dtype=np.float64)
        maxima -= dftshift
        shifts = shifts + maxima / upsample_factor
        CCmax = cross_correlation.max()
        src_amp = _upsampled_dft(src_freq * src_freq.conj(),
                                 1, upsample_factor)[0, 0]
        src_amp /= normalization
        target_amp = _upsampled_dft(target_freq * target_freq.conj(),
                                    1, upsample_factor)[0, 0]
        target_amp /= normalization

    # If its only one row or column the shift along that dimension has no
    # effect. We set to zero.
    for dim in range(src_freq.ndim):
        if shape[dim] == 1:
            shifts[dim] = 0

    return shifts, _compute_error(CCmax, src_amp, target_amp),\
        _compute_phasediff(CCmax)


def center_of_gravity(data1, data2, clip=None):
    import scipy.ndimage as nd
    mask1 = data1 > clip*np.nanmax(data1)
    mask2 = data2 > clip*np.nanmax(data2)
    d1 = np.where(mask1*mask2, data1, np.nan)
    d2 = np.where(mask1*mask2, data2, np.nan)
    b1x, b1y = nd.measurements.center_of_mass(d1, ~np.isnan(d1))
    b2x, b2y = nd.measurements.center_of_mass(d2, ~np.isnan(d2))
    return b1x-b2x, b1y-b2y

from scipy.ndimage import fourier_shift
def subpixel_shift(data1, data2, clip=0.1):
    dx, dy = center_of_gravity(data1, data2, clip=clip)
    nanz = 0
    dxs, dys = dx, dy
    d2 = data2.copy()
    while (abs(dx) > 0.1 or abs(dy) > 0.1) and nanz < 100:
          data2 = np.fft.ifft2(fourier_shift(np.fft.fft2(data2), (dx, dy))).real
          dx, dy = center_of_gravity(data1, data2, clip=0.1)
          dxs += dx
          dys += dy
          nanz += 1
    data2 = np.fft.ifft2(fourier_shift(np.fft.fft2(d2), (dxs, dys))).real
    return dxs, dys, data2


def rebin1(a, size):
    sh = 1,1,size,a.shape[0]//size
    return a.reshape(sh).mean(-1)

def rebin(a, shape):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).mean(-1).mean(1)

def hist_range_threshold(hist, bin_edges, percent):
    hist = np.concatenate((hist, [0]))
    threshold = .5*percent/100*hist.sum()
    i_bin_min = np.cumsum(hist).searchsorted(threshold)
    i_bin_max = -1-np.cumsum(np.flipud(hist)).searchsorted(threshold)
    return bin_edges[i_bin_min], bin_edges[i_bin_max]

def color_get_histogram(data, nbins):
    from guiqwt._scaler import _histogram
    if data is None:
       return [0,], [0,1]
    _min = np.nanmin(data)
    _max = np.nanmax(data)
    bins = np.unique(np.array(np.linspace(_min, _max, nbins+1), dtype=data.dtype))
    res2 = np.zeros((bins.size+1,), np.uint32)
    _histogram(data.flatten(), bins, res2)
    res = res2[1:-1], bins
    return res

def lut_range_threshold(data, nbins, percent):
    hist, bin_edges = color_get_histogram(data, nbins)
    return hist_range_threshold(hist, bin_edges, percent)

def plotPDF(data, cmap, percent, rebin=1):
    from matplotlib import pylab as plt
    import matplotlib
    import Image
    from nodmath import nan_interpolation, map_zoom, hist_range_threshold, color_get_histogram,\
                    lut_range_threshold
    xmin, xmax = lut_range_threshold(data, 256, percent)
    data = np.clip(data[3:-3, 3:-3], xmin, xmax)
    rows, cols = data.shape
    if rebin > 1:
       data = map_zoom(data, rebin, order=3, prefilter=True)
    m = matplotlib.cm.ScalarMappable(norm=None, cmap=cmap)
    colormapped = m.to_rgba(data)*255
    outputImage = Image.fromarray(np.uint8(colormapped))
    outputImage.save("pixmap.png")
    #outputImage.show() 
    #fig = plt.figure()
    #ax = fig.add_axes((0,0,1,1))
    #ax.set_axis_off()
    #ax.matshow(data, cmap=cmap)
    #plt.show()

import scipy, scipy.signal

def savitzky_golay( y, window_size, order, deriv = 0 ):
    r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
    The Savitzky-Golay filter removes high frequency noise from data.
    It has the advantage of preserving the original shape and
    features of the signal better than other types of filtering
    approaches, such as moving averages techhniques.
    Parameters
    ----------
    y : array_like, shape (N,)
        the values of the time history of the signal.
    window_size : int
        the length of the window. Must be an odd integer number.
    order : int
        the order of the polynomial used in the filtering.
        Must be less then `window_size` - 1.
    deriv: int
        the order of the derivative to compute (default = 0 means only smoothing)
    Returns
    -------
    ys : ndarray, shape (N)
        the smoothed signal (or it's n-th derivative).
    Notes
    -----
    The Savitzky-Golay is a type of low-pass filter, particularly
    suited for smoothing noisy data. The main idea behind this
    approach is to make for each point a least-square fit with a
    polynomial of high order over a odd-sized window centered at
    the point.
    Examples
    --------
    t = np.linspace(-4, 4, 500)
    y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
    ysg = savitzky_golay(y, window_size=31, order=4)
    import matplotlib.pyplot as plt
    plt.plot(t, y, label='Noisy signal')
    plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
    plt.plot(t, ysg, 'r', label='Filtered signal')
    plt.legend()
    plt.show()
    References
    ----------
    .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of
       Data by Simplified Least Squares Procedures. Analytical
       Chemistry, 1964, 36 (8), pp 1627-1639.
    .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
       W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
       Cambridge University Press ISBN-13: 9780521880688
    """
    try:
        window_size = np.abs( np.int( window_size ) )
        order = np.abs( np.int( order ) )
    except ValueError, msg:
        raise ValueError( "window_size and order have to be of type int" )
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError( "window_size size must be a positive odd number" )
    if window_size < order + 2:
        raise TypeError( "window_size is too small for the polynomials order" )
    order_range = range( order + 1 )
    half_window = ( window_size - 1 ) // 2
    # precompute coefficients
    b = np.mat( [[k ** i for i in order_range] for k in range( -half_window, half_window + 1 )] )
    m = np.linalg.pinv( b ).A[deriv]
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = y[0] - np.abs( y[1:half_window + 1][::-1] - y[0] )
    lastvals = y[-1] + np.abs( y[-half_window - 1:-1][::-1] - y[-1] )
    y = np.concatenate( ( firstvals, y, lastvals ) )
    return np.convolve( m, y, mode = 'valid' )

def savitzky_golay_piecewise( xvals, data, kernel = 11, order = 4 ):
    turnpoint = 0
    last = len( xvals )
    if xvals[1] > xvals[0] : #x is increasing?
        for i in range( 1, last ) : #yes
            if xvals[i] < xvals[i - 1] : #search where x starts to fall
                turnpoint = i
                break
    else: #no, x is decreasing
        for i in range( 1, last ) : #search where it starts to rise
            if xvals[i] > xvals[i - 1] :
                turnpoint = i
                break
    if turnpoint == 0 : #no change in direction of x
        return savitzky_golay( data, kernel, order )
    else:
        #smooth the first piece
        firstpart = savitzky_golay( data[0:turnpoint], kernel, order )
        #recursively smooth the rest
        rest = savitzky_golay_piecewise( xvals[turnpoint:], data[turnpoint:], kernel, order )
        return numpy.concatenate( ( firstpart, rest ) )

def sgolay2d ( z, window_size, order, derivative = None ):
    """
    """
    # number of terms in the polynomial expression
    n_terms = ( order + 1 ) * ( order + 2 ) / 2.0

    if  window_size % 2 == 0:
        raise ValueError( 'window_size must be odd' )

    if window_size ** 2 < n_terms:
        raise ValueError( 'order is too high for the window size' )

    half_size = window_size // 2

    # exponents of the polynomial.
    # p(x,y) = a0 + a1*x + a2*y + a3*x^2 + a4*y^2 + a5*x*y + ...
    # this line gives a list of two item tuple. Each tuple contains
    # the exponents of the k-th term. First element of tuple is for x
    # second element for y.
    # Ex. exps = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ...]
    exps = [ ( k - n, n ) for k in range( order + 1 ) for n in range( k + 1 ) ]

    # coordinates of points
    ind = np.arange( -half_size, half_size + 1, dtype = np.float64 )
    dx = np.repeat( ind, window_size )
    dy = np.tile( ind, [window_size, 1] ).reshape( window_size ** 2, )

    # build matrix of system of equation
    A = np.empty( ( window_size ** 2, len( exps ) ) )
    for i, exp in enumerate( exps ):
        A[:, i] = ( dx ** exp[0] ) * ( dy ** exp[1] )

    # pad input array with appropriate values at the four borders
    new_shape = z.shape[0] + 2 * half_size, z.shape[1] + 2 * half_size
    Z = np.zeros( ( new_shape ) )
    # top band
    band = z[0, :]
    Z[:half_size, half_size:-half_size] = band - np.abs( np.flipud( z[1:half_size + 1, :] ) - band )
    # bottom band
    band = z[-1, :]
    Z[-half_size:, half_size:-half_size] = band + np.abs( np.flipud( z[-half_size - 1:-1, :] ) - band )
    # left band
    band = np.tile( z[:, 0].reshape( -1, 1 ), [1, half_size] )
    Z[half_size:-half_size, :half_size] = band - np.abs( np.fliplr( z[:, 1:half_size + 1] ) - band )
    # right band
    band = np.tile( z[:, -1].reshape( -1, 1 ), [1, half_size] )
    Z[half_size:-half_size, -half_size:] = band + np.abs( np.fliplr( z[:, -half_size - 1:-1] ) - band )
    # central band
    Z[half_size:-half_size, half_size:-half_size] = z

    # top left corner
    band = z[0, 0]
    Z[:half_size, :half_size] = band - np.abs( np.flipud( np.fliplr( z[1:half_size + 1, 1:half_size + 1] ) ) - band )
    # bottom right corner
    band = z[-1, -1]
    Z[-half_size:, -half_size:] = band + np.abs( np.flipud( np.fliplr( z[-half_size - 1:-1, -half_size - 1:-1] ) ) - band )

    # top right corner
    band = Z[half_size, -half_size:]
    Z[:half_size, -half_size:] = band - np.abs( np.flipud( Z[half_size + 1:2 * half_size + 1, -half_size:] ) - band )
    # bottom left corner
    band = Z[-half_size:, half_size].reshape( -1, 1 )
    Z[-half_size:, :half_size] = band - np.abs( np.fliplr( Z[-half_size:, half_size + 1:2 * half_size + 1] ) - band )

    # solve system and convolve
    if derivative == None:
        m = np.linalg.pinv( A )[0].reshape( ( window_size, -1 ) )
        return scipy.signal.fftconvolve( Z, m, mode = 'valid' )
    elif derivative == 'col':
        c = np.linalg.pinv( A )[1].reshape( ( window_size, -1 ) )
        return scipy.signal.fftconvolve( Z, -c, mode = 'valid' )
    elif derivative == 'row':
        r = np.linalg.pinv( A )[2].reshape( ( window_size, -1 ) )
        return scipy.signal.fftconvolve( Z, -r, mode = 'valid' )
    elif derivative == 'both':
        c = np.linalg.pinv( A )[1].reshape( ( window_size, -1 ) )
        r = np.linalg.pinv( A )[2].reshape( ( window_size, -1 ) )
        return scipy.signal.fftconvolve( Z, -r, mode = 'valid' ), scipy.signal.fftconvolve( Z, -c, mode = 'valid' )


def main():
    Z = np.arange(0,36).reshape((6,6))
    print Z
    print
    print "extract(Z, shape=(3,3), position=(0,0))"
    print extract(Z, shape=(3,3), position=(0,0))
    print
    print "extract(Z, shape=(3,3), position=(3,3))"
    print extract(Z, shape=(3,3), position=(3,3))
    print
    print "extract(Z, shape=(10,10), position=(3,3))"
    print extract(Z, shape=(10,10), position=(3,3))

if __name__ == '__main__':
    main()