def pw_affine(fromim,toim,fp,tp,tri):
    """ Warp triangular patches from an image.
        fromim = image to warp 
        toim = destination image
        fp = from points in hom. coordinates
        tp = to points in hom.  coordinates
        tri = triangulation. """
                
    im = toim.copy()
    
    # check if image is grayscale or color
    is_color = len(fromim.shape) == 3
    
    # create image to warp to (needed if iterate colors)
    im_t = zeros(im.shape, 'uint8') 
    
    for t in tri:
        # compute affine transformation
        H = homography.Haffine_from_points(tp[:,t],fp[:,t])
        
        if is_color:
            for col in range(fromim.shape[2]):
                im_t[:,:,col] = ndimage.affine_transform(
                    fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        else:
            im_t = ndimage.affine_transform(
                    fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        
        # alpha for triangle
        alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
        
        # add triangle to image
        im[alpha>0] = im_t[alpha>0]
        
    return im

def pw_affine(fromim,toim,fp,tp,tri):
    """ 画像の三角形パッチを変形する。
        fromim = 変形する画像
        toim = 画像の合成先
        fp = 基準点（同次座標系）
        tp = 対応店（同次座標系）
        tri = 三角形分割　"""
    
    im = toim.copy()
    
    # 画像がグレースケールかカラーか調べる
    is_color = len(fromim.shape) == 3
    
    #　変形する先の画像を作成する
    im_t = np.zeros(im.shape, 'uint8')
    for t in tri:
        #アフィン変換を計算する
        H = homography.Haffine_from_points(tp[:,t],fp[:,t])
        if is_color:
            for col in range(fromim.shape[2]):
                im_t[:,:,col] = ndimage.affine_transform(
                    fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        else:
            im_t = ndimage.affine_transform(
                fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        
        # 三角形の透明度マップ
        alpha = alpha_for_triangle(tp[:,t].astype('int'),im.shape[0],im.shape[1])
        
        # 三角形を画像に追加する
        im[alpha>0] = im_t[alpha>0]
    
    return im

def _resample_one_img(data, A, b, target_shape,
                      interpolation_order, out, copy=True):
    "Internal function for resample_img, do not use"
    if data.dtype.kind in ('i', 'u'):
        # Integers are always finite
        has_not_finite = False
    else:
        not_finite = np.logical_not(np.isfinite(data))
        has_not_finite = np.any(not_finite)
    if has_not_finite:
        warnings.warn("NaNs or infinite values are present in the data "
                        "passed to resample. This is a bad thing as they "
                        "make resampling ill-defined and much slower.",
                        RuntimeWarning, stacklevel=2)
        if copy:
            # We need to do a copy to avoid modifying the input
            # array
            data = data.copy()
        #data[not_finite] = 0
        from ..masking import _extrapolate_out_mask
        data = _extrapolate_out_mask(data, np.logical_not(not_finite),
                                     iterations=2)[0]

    # See https://github.com/nilearn/nilearn/issues/346 Copying the
    # array makes it C continuous and as such the int32 index in the C
    # code is a lot less likely to overflow
    if (LooseVersion(scipy.__version__) < LooseVersion('0.14.1')):
        data = data.copy()

    # Suppresses warnings in https://github.com/nilearn/nilearn/issues/1363
    with warnings.catch_warnings():
        if LooseVersion(scipy.__version__) >= LooseVersion('0.18'):
            warnings.simplefilter("ignore", UserWarning)
        # The resampling itself
        ndimage.affine_transform(data, A,
                                 offset=b,
                                 output_shape=target_shape,
                                 output=out,
                                 order=interpolation_order)

    # Bug in ndimage.affine_transform when out does not have native endianness
    # see https://github.com/nilearn/nilearn/issues/275
    # Bug was fixed in scipy 0.15
    if (LooseVersion(scipy.__version__) < LooseVersion('0.15') and
        not out.dtype.isnative):
        out.byteswap(True)

    if has_not_finite:
        # Suppresses warnings in https://github.com/nilearn/nilearn/issues/1363
        with warnings.catch_warnings():
            if LooseVersion(scipy.__version__) >= LooseVersion('0.18'):
                warnings.simplefilter("ignore", UserWarning)
            # We need to resample the mask of not_finite values
            not_finite = ndimage.affine_transform(not_finite, A,
                                                offset=b,
                                                output_shape=target_shape,
                                                order=0)
        out[not_finite] = np.nan
    return out

    def transform(self, transform, grid_coords=False, reference=None, 
                  dtype=None, interp_order=_INTERP_ORDER):
        """
        Apply a transformation to the image considered as 'floating'
        to bring it into the same grid as a given 'reference'
        image. The transformation is assumed to go from the
        'reference' to the 'floating'.
        
        transform: nd array
    
        either a 4x4 matrix describing an affine transformation
        
        or a 3xN array describing voxelwise displacements of the
        reference grid points
        
        precomputed : boolean
        True for a precomputed transformation, False for affine

        grid_coords : boolean

        True if the transform maps to grid coordinates, False if it maps
        to world coordinates
    
        reference: reference image, defaults to input. 
        """
        if reference == None: 
            reference = self
        
        if dtype == None: 
            dtype = self._get_dtype()

        # Prepare data arrays
        data = self._get_data()
        output = np.zeros(reference._shape, dtype=dtype)
        t = np.asarray(transform)

        # Case: affine transform
        if t.shape[-1] == 4: 
            if not grid_coords:
                t = np.dot(self._inv_affine, np.dot(t, reference._affine))
            ndimage.affine_transform(data, t[0:3,0:3], offset=t[0:3,3],
                                     order=interp_order, cval=self._background, 
                                     output_shape=output.shape, output=output)
    
        # Case: precomputed displacements
        else:
            if not grid_coords:
                t = apply_affine(self._inv_affine, t)
            output = ndimage.map_coordinates(data, np.rollaxis(t, 3, 0), 
                                             order=interp_order, 
                                             cval=self._background,
                                             output=dtype)
    
        return Image(output, affine=reference._affine, world=self._world)

def testRigidTransformEstimation(inImg, level, dTheta, displacement, thr):
    left=ndimage.rotate(inImg, dTheta)
    right=ndimage.rotate(inImg, -dTheta)
    left=ndimage.affine_transform(left , np.eye(2), offset=-1*displacement)
    right=ndimage.affine_transform(right, np.eye(2), offset=displacement)
    
    rightPyramid=[i for i in transform.pyramid_gaussian(right, level)]
    leftPyramid=[i for i in transform.pyramid_gaussian(left, level)]
    sel=level
    beta=estimateRigidTransformation(leftPyramid[sel], rightPyramid[sel], 2.0*dTheta, thr)
    return beta

def affine_transform2(im, rot, shift):
    '''
        Perform affine transform for 2/3D images.
    '''
    if ndim(im) == 2:
        return ndimage.affine_transform(im, rot, shift)
    else:
        imr = ndimage.affine_transform(im[:, :, 0], rot, shift)
        img = ndimage.affine_transform(im[:, :, 1], rot, shift)
        imb = ndimage.affine_transform(im[:, :, 2], rot, shift)

        return dstack((imr, img, imb))

def fast_warp_forward(img, tf, output_shape=(50, 50), mode='constant', order=1):
    """
    This wrapper function is faster than skimage.transform.warp
    """
    m = tf._inv_matrix
    res = np.zeros(shape=(output_shape[0], output_shape[1], 3), dtype=floatX)
    from scipy.ndimage import affine_transform
    trans, offset = m[:2, :2], (m[0, 2], m[1, 2])
    res[:, :, 0] = affine_transform(img[:, :, 0].T, trans, offset=offset, output_shape=output_shape, mode=mode, order=order)
    res[:, :, 1] = affine_transform(img[:, :, 1].T, trans, offset=offset, output_shape=output_shape, mode=mode, order=order)
    res[:, :, 2] = affine_transform(img[:, :, 2].T, trans, offset=offset, output_shape=output_shape, mode=mode, order=order)
    return res

    def apply_resize(self, workspace, input_image_name, output_image_name):
        image = workspace.image_set.get_image(input_image_name)
        image_pixels = image.pixel_data
        if self.size_method == R_BY_FACTOR:
            factor = self.resizing_factor.value
            shape = (np.array(image_pixels.shape[:2])*factor+.5).astype(int)
        elif self.size_method == R_TO_SIZE:
            if self.use_manual_or_image == C_MANUAL:
                shape = np.array([self.specific_height.value,
                                  self.specific_width.value])
            elif self.use_manual_or_image == C_IMAGE:
                shape = np.array(workspace.image_set.get_image(self.specific_image.value).pixel_data.shape).astype(int)
        #
        # Little bit of wierdness here. The input pixels are numbered 0 to
        # shape-1 and so are the output pixels. Therefore the affine transform
        # is the ratio of the two shapes-1
        #
        ratio = ((np.array(image_pixels.shape[:2]).astype(float)-1) /
                 (shape.astype(float)-1))
        transform = np.array([[ratio[0], 0],[0,ratio[1]]])
        if self.interpolation not in I_ALL:
            raise NotImplementedError("Unsupported interpolation method: %s" %
                                      self.interpolation.value)
        order = (0 if self.interpolation == I_NEAREST_NEIGHBOR
                 else 1 if self.interpolation == I_BILINEAR
                 else 2)
        if image_pixels.ndim == 3:
            output_pixels = np.zeros((shape[0],shape[1],image_pixels.shape[2]), 
                                     image_pixels.dtype)
            for i in range(image_pixels.shape[2]):
                affine_transform(image_pixels[:,:,i], transform,
                                 output_shape = tuple(shape),
                                 output = output_pixels[:,:,i],
                                 order = order)
        else:
            output_pixels = affine_transform(image_pixels, transform,
                                             output_shape = shape,
                                             order = order)
        output_image = cpi.Image(output_pixels, parent_image=image)
        workspace.image_set.add(output_image_name, output_image)

        if self.show_window:
            if not hasattr(workspace.display_data, 'input_images'):
                workspace.display_data.input_images = [image.pixel_data]
                workspace.display_data.output_images = [output_image.pixel_data]
                workspace.display_data.input_image_names = [input_image_name]
                workspace.display_data.output_image_names = [output_image_name]
            else:
                workspace.display_data.input_images += [image.pixel_data]
                workspace.display_data.output_images += [output_image.pixel_data]
                workspace.display_data.input_image_names += [input_image_name]
                workspace.display_data.output_image_names += [output_image_name]

def scale_images(pixel_array, mask):

    param = random.uniform(0.95,1.05)
    S = trans.zooms.zfdir2mat(param)
    pixel_array = affine_transform(pixel_array,S)
    mask = affine_transform(mask,S)
    
    scaled_data = {
        "pixel_array" : pixel_array,
        "mask" : mask, 
        "params" : param
    }
    
    return scaled_data

def shear_images(pixel_array, mask):
    # only transform x axis?
    param = random.uniform(0, 0.05)
    S = [param, 0, 0]

    pixel_array = affine_transform(pixel_array,trans.shears.striu2mat(S))
    mask = affine_transform(mask,trans.shears.striu2mat(S))

    sheared_data = {
        "pixel_array" : pixel_array,
        "mask" : mask, 
        "params" : param
    }
    
    return sheared_data

def _resample_one_img(data, A, b, target_shape,
                      interpolation_order, out, copy=True,
                      fill_value=0):
    "Internal function for resample_img, do not use"
    if data.dtype.kind in ('i', 'u'):
        # Integers are always finite
        has_not_finite = False
    else:
        not_finite = np.logical_not(np.isfinite(data))
        has_not_finite = np.any(not_finite)
    if has_not_finite:
        warnings.warn("NaNs or infinite values are present in the data "
                        "passed to resample. This is a bad thing as they "
                        "make resampling ill-defined and much slower.",
                        RuntimeWarning, stacklevel=2)
        if copy:
            # We need to do a copy to avoid modifying the input
            # array
            data = data.copy()
        #data[not_finite] = 0
        from ..masking import _extrapolate_out_mask
        data = _extrapolate_out_mask(data, np.logical_not(not_finite),
                                     iterations=2)[0]

    # Suppresses warnings in https://github.com/nilearn/nilearn/issues/1363
    with warnings.catch_warnings():
        if LooseVersion(scipy.__version__) >= LooseVersion('0.18'):
            warnings.simplefilter("ignore", UserWarning)
        # The resampling itself
        ndimage.affine_transform(data, A,
                                 offset=b,
                                 output_shape=target_shape,
                                 output=out,
                                 cval=fill_value,
                                 order=interpolation_order)

    if has_not_finite:
        # Suppresses warnings in https://github.com/nilearn/nilearn/issues/1363
        with warnings.catch_warnings():
            if LooseVersion(scipy.__version__) >= LooseVersion('0.18'):
                warnings.simplefilter("ignore", UserWarning)
            # We need to resample the mask of not_finite values
            not_finite = ndimage.affine_transform(not_finite, A,
                                                offset=b,
                                                output_shape=target_shape,
                                                order=0)
        out[not_finite] = np.nan
    return out

def transformImage(img, xfactor, angle=0, msg=False):
        """
        rotates then stretches or compresses an image only along the x-axis
        """
        if xfactor > 1.0:
                mystr = "_S"
        else:
                mystr = "_C"

        if msg is True:
                if xfactor > 1:
                        apDisplay.printMsg("stretching image by "+str(round(xfactor,3)))
                else:
                        apDisplay.printMsg("compressing image by "+str(round(xfactor,3)))
        ### image has swapped coordinates (y,x) from particles
        transMat = numpy.array([[ 1.0, 0.0 ], [ 0.0, 1.0/xfactor ]])
        #print "transMat\n",transMat
        #apImage.arrayToJpeg(img, "img"+mystr+".jpg")

        stepimg  = ndimage.rotate(img, -1.0*angle, mode='reflect')
        stepimg = apImage.frame_cut(stepimg, img.shape)
        #apImage.arrayToJpeg(stepimg, "rotate"+mystr+".jpg")

        newimg  = ndimage.affine_transform(stepimg, transMat, mode='reflect')
        #apImage.arrayToJpeg(newimg, "last_transform"+mystr+".jpg")

        return newimg

def rigid_alignment(faces,path,plotflag=False):
  """ 画像を位置合わせし、新たな画像として保存する。
      pathは、位置合わせした画像の保存先
      plotflag=Trueなら、画像を表示する """

  # 最初の画像の点を参照点とする
  refpoints = faces.values()[0]

  # 各画像を相似変換で変形する
  for face in faces:
    points = faces[face]

    R,tx,ty = compute_rigid_transform(refpoints, points)
    T = array([[R[1][1], R[1][0]], [R[0][1], R[0][0]]])

    im = array(Image.open(os.path.join(path,face)))
    im2 = zeros(im.shape, 'uint8')

    # 色チャンネルごとに変形する
    for i in range(len(im.shape)):
      im2[:,:,i] = ndimage.affine_transform(im[:,:,i],linalg.inv(T),
                                            offset=[-ty,-tx])
    if plotflag:
      imshow(im2)
      show()

    # 境界で切り抜き、位置合わせした画像を保存する
    h,w = im2.shape[:2]
    border = (w+h)/20
    imsave(os.path.join(path, 'aligned/'+face),
          im2[border:h-border,border:w-border,:])

def test_map_coordinates_dts():
    # check that ndimage accepts different data types for interpolation
    data = np.array([[4, 1, 3, 2],
                     [7, 6, 8, 5],
                     [3, 5, 3, 6]])
    shifted_data = np.array([[0, 0, 0, 0],
                             [0, 4, 1, 3],
                             [0, 7, 6, 8]])
    idx = np.indices(data.shape)
    dts = (np.uint8, np.uint16, np.uint32, np.uint64,
           np.int8, np.int16, np.int32, np.int64,
           np.intp, np.uintp, np.float32, np.float64)
    for order in range(0, 6):
        for data_dt in dts:
            these_data = data.astype(data_dt)
            for coord_dt in dts:
                # affine mapping
                mat = np.eye(2, dtype=coord_dt)
                off = np.zeros((2,), dtype=coord_dt)
                out = ndimage.affine_transform(these_data, mat, off)
                assert_array_almost_equal(these_data, out)
                # map coordinates
                coords_m1 = idx.astype(coord_dt) - 1
                coords_p10 = idx.astype(coord_dt) + 10
                out = ndimage.map_coordinates(these_data, coords_m1, order=order)
                assert_array_almost_equal(out, shifted_data)
                # check constant fill works
                out = ndimage.map_coordinates(these_data, coords_p10, order=order)
                assert_array_almost_equal(out, np.zeros((3,4)))
            # check shift and zoom
            out = ndimage.shift(these_data, 1)
            assert_array_almost_equal(out, shifted_data)
            out = ndimage.zoom(these_data, 1)
            assert_array_almost_equal(these_data, out)

def rigid_alignment(faces,path,plotflag=False):
    """    Align images rigidly and save as new images.
        path determines where the aligned images are saved
        set plotflag=True to plot the images. """
    
    # take the points in the first image as reference points
    refpoints = faces.values()[0]
    
    # warp each image using affine transform
    for face in faces:
        points = faces[face]
        
        R,tx,ty = compute_rigid_transform(refpoints, points)
        T = array([[R[1][1], R[1][0]], [R[0][1], R[0][0]]])    
        
        im = array(Image.open(os.path.join(path,face)))
        im2 = zeros(im.shape, 'uint8')
        
        # warp each color channel
        for i in range(len(im.shape)):
            im2[:,:,i] = ndimage.affine_transform(im[:,:,i],linalg.inv(T),offset=[-ty,-tx])
            
        if plotflag:
            imshow(im2)
            show()
            
        # crop away border and save aligned images
        h,w = im2.shape[:2]
        border = (w+h)/20
        
        # crop away border
        imsave(os.path.join(path, 'aligned/'+face),im2[border:h-border,border:w-border,:])

def makeDeepDreamStepped(img, classifier, end, name, iterations, octaves, octave_scale, start_sigma, end_sigma, start_jitter, end_jitter, start_step_size, end_step_size, scaleZoom):
	newImagePath = pathToOutput+'/'+name+"_i"+str(iterations)+"_o"+str(octaves)+"_os"+str(octave_scale)+"_j"+str(start_jitter)+'_'+str(end_jitter)+'.png'
	frame = deepdream_stepped(classifier, img, iterations, octaves, octave_scale, end, start_sigma, end_sigma, start_jitter, end_jitter, start_step_size, end_step_size)
	PIL.Image.fromarray(np.uint8(np.clip(frame, 0, 255))).save(newImagePath, 'png')
	h, w = frame.shape[:2]
	frame = nd.affine_transform(frame, [1-scaleZoom,1-scaleZoom,1], [h*scaleZoom/2,w*scaleZoom/2,0], order=1)
	return frame

def test_flirt2aff():
    from os.path import join as pjoin
    from nose.tools import assert_true
    import scipy.ndimage as ndi
    import nibabel as nib
    
    '''
    matfile = pjoin('fa_data',
                    '1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_affine_transf.mat')
    in_fname = pjoin('fa_data',
                     '1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_bet_FA.nii.gz')
    '''
    
    matfile=flirtaff
    in_fname = ffa
    
    ref_fname = '/usr/share/fsl/data/standard/FMRIB58_FA_1mm.nii.gz'
    res = flirt2aff_files(matfile, in_fname, ref_fname)
    mat = np.loadtxt(matfile)
    in_img = nib.load(in_fname)
    ref_img = nib.load(ref_fname)
    assert_true(np.all(res == flirt2aff(mat, in_img, ref_img)))
    # mm to mm transform
    # mm_in2mm_ref =  np.dot(ref_img.affine,
    #                        np.dot(res, npl.inv(in_img.affine)))
    # make new in image thus transformed
    in_data = in_img.get_data()
    ires = npl.inv(res)
    in_data[np.isnan(in_data)] = 0
    resliced_data = ndi.affine_transform(in_data,
                                         ires[:3,:3],
                                         ires[:3,3],
                                         ref_img.shape)
    resliced_img = nib.Nifti1Image(resliced_data, ref_img.affine)
    nib.save(resliced_img, 'test.nii')

    def transform(self,x):

#        import matplotlib.pyplot as plt
#
#        print x[:,:,0].min()
#        print x[:,:,1].min()
#        print x[:,:,2].min()
#        print x[:,:,0].max()
#        print x[:,:,1].max()
#        print x[:,:,2].max()
#
##        plt.imshow(x.astype('uint8'))
#        plt.imshow(x)
#        plt.show()

        rx = gen_mm_random(np.random.__dict__[self.random_fct],self.fct_settings,self.min,self.max)
        ry = gen_mm_random(np.random.__dict__[self.random_fct],self.fct_settings,self.min,self.max)

#        x = ndimage.affine_transform(x.astype('uint8'),self.I,offset=(rx,ry,0))
        x = ndimage.affine_transform(x,self.I,offset=(rx,ry,0))

#        plt.imshow(x.astype('uint8'))
#        plt.imshow(x)
#        plt.show()
        return x

def make_truth(img_file='./truth.png'):
    """
    Load in the truth image data, embed it into a 3d array, then rotate it in a
    weird way
    """
    pixels = 1.0 - mplimg.imread(img_file)[:, :, 0]  # b&w image, just grab red

    # Now embed in 3d array and rotate in some interesting way
    voxels = np.zeros(pixels.shape + (pixels.shape[0], ))
    voxels[:, :, voxels.shape[2] // 2] = pixels
    rot_ang_axis = np.array((2, 1, 0.5))  # Something "interesting"
    aff_matrix = angle_axis_to_matrix(rot_ang_axis)
    # Rotate about center, but affine_transform offset parameter is dumb
    center = np.array(voxels.shape) / 2  # whatever close enough
    offset = -(center - center.dot(aff_matrix)).dot(np.linalg.inv(aff_matrix))
    voxels = ndimage.affine_transform(voxels, aff_matrix, offset=offset)

    # Remake the truth figure in 3D
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    x, y, z = np.meshgrid(np.arange(voxels.shape[0]),
                          np.arange(voxels.shape[1]),
                          np.arange(voxels.shape[2]),
                          indexing='ij')
    disp_vox = voxels > 0.3
    ax.scatter(x[disp_vox], y[disp_vox], z[disp_vox])
    plt.savefig(img_file.replace('.png', '_3d.png'))
    # plt.show()

    print("truth:", voxels.shape)
    return voxels

def normalize_rotation(img,shape=(30,30),normalize=1):
    if type(img)==ListType:
        return [normalize_rotation(image,shape=shape) for image in img]
    # rescale into 0-1 range
    image = normalize_range(img)
    # compute image statistics
    cx,cy,cxx,cxy,cyy = compute_stats(image)
    # compute eigenvectors
    mat = array([[cxx,cxy],[cxy,cyy]])
    v,d = linalg.eig(mat)
    alpha0 = atan2upper(d[1,0],d[0,0]) - pi/2
    alpha1 = atan2upper(d[1,1],d[1,0]) - pi/2
    if abs(alpha0)<abs(alpha1):
        alpha = -alpha0
    else:
        alpha = -alpha1
    # perform the actual transformation
    affine = array([[cos(alpha),-sin(alpha)],[sin(alpha),cos(alpha)]])
    ccenter = array((cx,cy))
    ocenter = array((shape[0]/2,shape[1]/2))
    offset = ccenter - dot(affine,ocenter)
    if normalize: img = image
    return ndimage.affine_transform(img,affine,
                                    offset=offset,
                                    output_shape=shape)