def build_rnn_dataset(steps):
    rnntrain,rnntest=makernnset(steps)
    print "trainset.shape, testset.shape =", rnntrain.shape, rnntest.shape
    X_train, y_train = np.dsplit(rnntrain,[5])
    X_valid, y_valid = np.dsplit(rnntest,[5])
                                           
    X_train = normalize(X_train)
    X_valid = normalize(X_valid)
    print 'X_train.shape, y_train.shape =', X_train.shape, y_train.shape
    return X_train, y_train, X_valid, y_valid

def generate_rgb_array(filename):
	imageframes = []

	if len(filename) == 1:
		filename = filename[0]



	# reading tif image into memory and turning into numpy array
	im = Image.open(filename)
	im.load()
	imarray = np.array(im)
	print np.unique(np.dsplit(imarray,3)[0]),np.unique(np.dsplit(imarray,3)[1]),np.unique(np.dsplit(imarray,3)[2])

	return imarray

def cat_trials(x):
    """ Concatenate trials along time axis.

    Parameters
    ----------
    x : array_like
        Segmented input data of shape [`n`,`m`,`t`], with `n` time samples, `m` signals, and `t` trials.

    Returns
    -------
    out : ndarray
        Trials are concatenated along the first (time) axis. Shape of the output is [`n``t`,`m`].

    See also
    --------
    cut_segments : Cut segments from continuous data

    Examples
    --------
    >>> x = np.random.randn(150, 4, 6)
    >>> y = cat_trials(x)
    >>> y.shape
    (900, 4)
    """
    x = np.atleast_3d(x)
    t = x.shape[2]
    return np.squeeze(np.vstack(np.dsplit(x, t)), axis=2)

def median_filter(image, selem=None):
    if selem is None:
        # default mask is 5x5 square
        selem = square(5)
    depth = image.shape[2]
    return np.dstack(median(channel[...,0], selem)
                     for channel in np.dsplit(image, depth)) / 255.

def conditional_maximum_correlation_pmf(pmf):
    """
    Compute the conditional maximum correlation from a 3-dimensional
    pmf. The maximum correlation is computed between the first two dimensions
    given the third.

    Parameters
    ----------
    pmf : np.ndarray
        The probability distribution.

    Returns
    -------
    rho_max : float
        The conditional maximum correlation.
    """
    pXYgZ = pmf / pmf.sum(axis=(0,1), keepdims=True)
    pXgZ = pXYgZ.sum(axis=1, keepdims=True)
    pYgZ = pXYgZ.sum(axis=0, keepdims=True)
    Q = np.where(pmf, pXYgZ / (np.sqrt(pXgZ)*np.sqrt(pYgZ)), 0)
    Q[np.isnan(Q)] = 0

    rho_max = max([svdvals(np.squeeze(m))[1] for m in np.dsplit(Q, Q.shape[2])])

    return rho_max

def fun_add_channels(img, weights):
    """
    A functional implementation of the same function
    """
    return (reduce(lambda a, b: a+b,
                   [weight * channel for weight, channel in list(
                    zip(weights, np.dsplit(img, img.shape[2])))]))[:, :, 0]

def IMT_find_col(imin, col):
    # This is the color distance from the reference point
    coldist = imin.colorDistance( col )
    lolout = Image( npboost(  np.squeeze(np.dsplit( coldist.getNumpy(), 3)[0] )  ) )
    modim = lolout.stretch(0,20)
    openim = m_open(modim,1).binarize()
    return IMT_calccentroid(openim)

def main():
	# 构造一个数组
	a = arange(9).reshape(3,3)
	print a
	# [[0 1 2]
	#  [3 4 5]
	#  [6 7 8]]


	# 横向拆分
	hs = hsplit(a,3)
	print hs
	# [
	# 	array([[0],[3],[6]]), 
	# 	array([[1],[4],[7]]), 
	# 	array([[2], [5], [8]])
	# ]


	# 纵向拆分
	vs = vsplit(a,3)
	print vs
	# [
	# 	array([[0, 1, 2]]), 
	# 	array([[3, 4, 5]]), 
	# 	array([[6, 7, 8]])
	# ]


	# 深向拆分,数组需要为三维数组
	b = arange(27).reshape(3,3,3)
	ds = dsplit(b,3)
	print ds

def surf2CV(surf, cvImage):
	"""
	Given a Pygame surface, convert to an OpenCv cvArray format.
	Either Ipl image or cvMat.
	surf2CV( pygame.Surface src, cv.Image dest )
	(From http://facial-expression-recognition.googlecode.com/svn/trunk/code/conversion.py)
	( Extracted 2012-Jul-16 22:37EDT by GKF)
	"""
	from numpy import dsplit, dstack
	cv.Set(cvImage, (0,0,0))
	arr = pygame.surfarray.pixels3d(surf).transpose(1,0,2) # Reshape to 320x240
	r,g,b = dsplit(arr,3)
	arr = dstack((b,g,r))
	dtype2depth = {
	'uint8': cv.IPL_DEPTH_8U,
	'int8': cv.IPL_DEPTH_8S,
	'uint16': cv.IPL_DEPTH_16U,
	'int16': cv.IPL_DEPTH_16S,
	'int32': cv.IPL_DEPTH_32S,
	'float32': cv.IPL_DEPTH_32F,
	'float64': cv.IPL_DEPTH_64F,
	}
	try:
		nChannels = arr.shape[2]
	except:
		nChannels = 3
	try:
		cv.SetData(cvImage, arr.tostring(),arr.dtype.itemsize*nChannels*arr.shape[1])
	except:
		print "Error is: ",
		print sys.exc_info()[0]

def multiply_3x3_mat(src, mat):
    """RGBの各ピクセルに対して3x3の行列演算を行う"""

    # 正規化用の係数を調査
    normalize_val = (2 ** (8 * src.itemsize)) - 1

    # 0 .. 1 に正規化して RGB分離
    b, g, r = np.dsplit(src / normalize_val, 3)

    # 行列計算
    ret_r = r * mat[0][0] + g * mat[0][1] + b * mat[0][2]
    ret_g = r * mat[1][0] + g * mat[1][1] + b * mat[1][2]
    ret_b = r * mat[2][0] + g * mat[2][1] + b * mat[2][2]

    # オーバーフロー確認(実は Matrixの係数を調整しているので不要)
    ret_r = cv2.min(ret_r, 1.0)
    ret_g = cv2.min(ret_g, 1.0)
    ret_b = cv2.min(ret_b, 1.0)

    # アンダーフロー確認(実は Matrixの係数を調整しているので不要)
    ret_r = cv2.max(ret_r, 0.0)
    ret_g = cv2.max(ret_g, 0.0)
    ret_b = cv2.max(ret_b, 0.0)

    # RGB結合
    ret_mat = np.dstack( (ret_b, ret_g, ret_r) )

    # 0 .. 255 に正規化
    ret_mat *= normalize_val

    return np.uint8(ret_mat)

def colors_peripheral_vs_central(image_roi, attrs={}, debug=False):
    image_roi, center = pad_for_rotation(image_roi)
    lesion_mask = image_roi[..., 3]

    goal = lesion_mask.sum() * 0.7
    inner = lesion_mask.copy()
    while inner.sum() > goal:
        inner = binary_erosion(inner, disk(1))
    outer = np.logical_and(lesion_mask, np.logical_not(inner))

    if debug:
        print """\
=== Colors Peripheral vs Central ===
lesion area: %d
inner goal: %d
inner area: %d
outer area: %d
""" % (lesion_mask.sum(), goal, inner.sum(), outer.sum())

    if debug:
        plt.subplot(131)
        plt.imshow(lesion_mask)
        plt.subplot(132)
        plt.imshow(inner)
        plt.subplot(133)
        plt.imshow(outer)
        plt.show()

    outer = np.nonzero(outer)
    inner = np.nonzero(inner)

    image_lab = rgb2lab(image_roi[..., :3])
    L, a, b = np.dsplit(image_lab, 3)

    delta_L = np.mean(L[outer]) - np.mean(L[inner])
    delta_a = np.mean(a[outer]) - np.mean(a[inner])
    delta_b = np.mean(b[outer]) - np.mean(b[inner])

    density_L = (
        np.histogram(L[outer], 100, (0.,100.), density=True)[0] *
        np.histogram(L[inner], 100, (0.,100.), density=True)[0]
    ).sum()
    density_a = (
        np.histogram(a[outer], 254, (-127.,127.), density=True)[0] *
        np.histogram(a[inner], 254, (-127.,127.), density=True)[0]
    ).sum()
    density_b = (
        np.histogram(b[outer], 254, (-127.,127.), density=True)[0] *
        np.histogram(b[inner], 254, (-127.,127.), density=True)[0]
    ).sum()

    attrs.update([
        ('Colors PvsC mean difference L', delta_L),
        ('Colors PvsC mean difference a', delta_a),
        ('Colors PvsC mean difference b', delta_b),
        ('Colors PvsC density baysian L', density_L),
        ('Colors PvsC density baysian a', density_a),
        ('Colors PvsC density baysian b', density_b),
    ])

def montage(vol, ncols=None):
    """Returns a 2d image monage given a 3d volume."""
    ncols = ncols if ncols else int(np.ceil(np.sqrt(vol.shape[2])))
    rows = np.array_split(vol, range(ncols,vol.shape[2],ncols), axis=2)
    # ensure the last row is the same size as the others
    rows[-1] = np.dstack((rows[-1], np.zeros(rows[-1].shape[0:2] + (rows[0].shape[2]-rows[-1].shape[2],))))
    im = np.vstack([np.squeeze(np.hstack(np.dsplit(row, ncols))) for row in rows])
    return(im)

def _swaplch(LCH):
    "Reverse the order of an LCH numpy dstack or tuple for analysis."
    try:  # Numpy array
        L, C, H = np.dsplit(LCH, 3)
        return np.dstack((H, C, L))
    except:  # Tuple
        L, C, H = LCH
        return H, C, L

def alpha_blend(image, background):
    "Une dos imagenes con opacidad, usando numpy"

    image, alpha = np.dsplit(image, np.array([3]))
    image = image
    alpha = 1 - alpha 
    resultado = image * alpha + background * (1 - alpha)
    return resultado

def rgb_to_hsv(img):
    h,w,d = img.shape
    r, g, b = np.dsplit(img,3)
    maxc = img.max(axis=2).reshape(h,w,1)
    minc = img.min(axis=2).reshape(h,w,1)
    s = (maxc-minc) / maxc
    v = maxc
    imgc = (maxc-img)/(maxc-minc)
    rc, gc, bc = np.dsplit(imgc,3)
    h = np.where(maxc==r, bc-gc,
        np.where(maxc==g, 2.0+rc-bc,
                             4.0+gc-rc))
    h = (h/6.0) % 1.0
    hsv = np.dstack([h,s,v])
    v0 = np.dstack([np.zeros_like(h),np.zeros_like(h),v])
    mask = minc == maxc
    mask = np.dstack([mask,mask,mask])
    return np.where(mask, v0, hsv)

def shift_image(im, shift):
    delta_y = shift[0]
    delta_x = shift[1]
    imOut = np.zeros(im.shape)
    for i, c in enumerate(np.dsplit(im, 3)):
        c = c[:, :, 0]
        Y = np.arange(c.shape[0])
        X = np.arange(c.shape[1])
        f = interp2d(X + delta_x, Y + delta_y, c)
        imOut[:, :, i] = f(X, Y)
    return imOut

    def write(cls, metadata, imagedata, outbase, voxel_order='LPS'):
        """
        Create png files for each image in a list of pixel data.

        Parameters
        ----------
        metadata : object
            fully loaded instance of a NIMSReader.
        imagedata : dict
            dictionary of np.darrays. label suffix as keys, with np.darrays as values.
        outbase : str
            output name prefix.
        voxel_order : str [default None]
            three character string indicating the voxel order, ex. 'LPS'.

        Returns
        -------
        results : list
            list of files written.

        Raises
        ------
        NIMSDataError
            metadata or data is None.

        """
        super(NIMSPNG, cls).write(metadata, imagedata, outbase, voxel_order)  # XXX FAIL! unexpected imagedata = None
        results = []
        for data_label, data in imagedata.iteritems():
            if data is None:
                continue
            if voxel_order and metadata.qto_xyz:  # cannot reorder if no affine
                data, qto_xyz = cls.reorder_voxels(data, metadata.qto_xyz, voxel_order)
            else:
                qto_xyz = metadata.qto_xyz
            outname = outbase + data_label
            data = np.dsplit(data, len(metadata._dcm_list))  # cut the darray
            data = [image.squeeze() for image in data]  # squeeze; remove axis with 1 val
            for i, data in enumerate(data):
                filepath = outname + '_%d' % (i + 1) + '.png'
                if data.ndim == 2:
                    data = data.astype(np.int32)
                    data = data.clip(0, (data * (data != (2**15 - 1))).max())  # -32768->0; 32767->brain.max
                    data = data * (2**8 - 1) / data.max()  # scale to full 8-bit range
                    Image.fromarray(data.astype(np.uint8), 'L').save(filepath, optimize=True)
                elif data.ndim == 3:
                    data = data.reshape((data.shape[1], data.shape[2], data.shape[0]))
                    Image.fromarray(data, 'RGB').save(filepath, optimize=True)
                log.debug('generated %s' % os.path.basename(filepath))
                results.append(filepath)
            log.debug('returning:  %s' % filepath)
        return results

def IMT_calccentroid(imin, *args):
    immod = np.squeeze(np.dsplit( imin.getNumpy().astype(float), 3)[0] )/255

    xvals = np.arange(0,imin.width)
    xstack = np.transpose( np.tile( xvals, (imin.height,1) ) )
    xmoments = np.multiply(xstack, immod).flatten()
    xbar = np.sum(xmoments)/np.sum(immod)

    yvals = np.arange(0,imin.height)
    ystack = np.tile( yvals, (imin.width,1) )
    ymoments = np.multiply(ystack, immod).flatten()
    ybar = np.sum(ymoments)/np.sum(immod)
    return [xbar,ybar]

def dir_tournant(x, y, a=0, b=0):
    """ Permet de generer un champ tournant vers la gauche
            et centre en (a,b)"""
    x, y = translate(x, y, a, b)
    if type(x) in [int, float, np.float64]:
        return normalize([-y, x])
    # vect = np.array()
    z = np.dstack((-y, x))
    U, V = np.dsplit(np.apply_along_axis(normalize, axis=2, arr=z), 2)
    U = np.apply_along_axis(lambda x: x[0], axis=2, arr=U)
    V = np.apply_along_axis(lambda x: x[0], axis=2, arr=V)
    # print vect, '-' * 100
    return np.array([U, V])

def lab_to_xyz(LAB, wp):

    L, a, b = np.dsplit(LAB, 3)
    fy = (L + 16) / 116.0
    fz = fy - b / 200.0
    fx = a / 500.0 + fy

    def finv(y):
        y = copy.copy(y)  # CEBALERT: why copy?
        eps3 = EPS ** 3
        return np.where(y > eps3, np.power(y, 3), (116 * y - 16) / KAP)

    xr, yr, zr = finv(fx), finv(fy), finv(fz)
    return np.dstack((xr * wp[0], yr * wp[1], zr * wp[2]))