def resc(patch):
    """
    :param patch:  [image,mask]
    :return: random rescaling of the pair [image,mask]

    --- Rescaling reinforces axons size diversity ---
    """


    s = random.choice([0.5, 0.75, 1.0, 1.5, 2.0])
    data_rescale=[]
    for scale in s:

        image_rescale = rescale(patch[0], scale)
        mask_rescale = rescale(patch[1], scale)
        s_r = mask_rescale.shape[0]
        q_h, r_h = divmod(256-s_r,2)

        if q_h > 0 :
            image_rescale = np.pad(image_rescale,(q_h, q_h+r_h), mode = "reflect")
            mask_rescale = np.pad(mask_rescale,(q_h, q_h+r_h), mode = "reflect")
        else :
            patches = extract_patch(image_rescale,mask_rescale, 256)
            i = np.random.randint(len(patches), size=1)
            image_rescale,mask_rescale = patches[i]

        mask_rescale = preprocessing.binarize(np.array(mask_rescale), threshold=0.001)
        data_rescale = [image_rescale, mask_rescale]

    return data_rescale

def read_image(name, size=None, debug=False):
    """ read image and segmentation, returns RGB + alpha composite """
    image = imread(name) / 255.

    if image.shape[2] == 4:
        alpha = image[...,3]
        image = image[...,:3]
    else:
        segmentation_name = os.path.splitext(name)[0][:-6] + '-label.png'
        segmentation = imread(segmentation_name)
        alpha = np.logical_or(segmentation[...,0], segmentation[...,1]) * 1.

    if size is not None:
        scale_x = float(size[0]) / image.shape[1]
        scale_y = float(size[1]) / image.shape[0]
        scale = min(scale_x, scale_y)

        if debug:
            print name, size[0], size[1], image.shape[1], image.shape[0], scale, image.shape[1]*scale, image.shape[0]*scale

        if scale > 1.0:
            print 'Image %s smaller than requested size' % name

        if scale != 1.0:
            image = rescale(image, scale, order=3)
            alpha = rescale(alpha, scale, order=0)

    return np.dstack((image, alpha))

def test_iradon_angles():
    """
    Test with different number of projections
    """
    size = 100
    # Synthetic data
    image = np.tri(size) + np.tri(size)[::-1]
    # Large number of projections: a good quality is expected
    nb_angles = 200
    radon_image_200 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                     endpoint=False))
    reconstructed = iradon(radon_image_200)
    delta_200 = np.mean(abs(rescale(image) - rescale(reconstructed)))
    assert delta_200 < 0.03
    # Lower number of projections
    nb_angles = 80
    radon_image_80 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                    endpoint=False))
    # Test whether the sum of all projections is approximately the same
    s = radon_image_80.sum(axis=0)
    assert np.allclose(s, s[0], rtol=0.01)
    reconstructed = iradon(radon_image_80)
    delta_80 = np.mean(abs(image / np.max(image) -
                           reconstructed / np.max(reconstructed)))
    # Loss of quality when the number of projections is reduced
    assert delta_80 > delta_200

def standardize_roi_height(roi, standard_height=20, rel_height=0.666):
    h = standard_height * 5

    y = roi.sum(1)
    yp = y[y >= y.max() * rel_height]

    pw = yp.shape[0]
    sf = standard_height / pw

    sh = int(np.ceil(float(roi.shape[0]) * sf))
    if sh <= h:                 # if the scale factor estimation isn't off try to rescale according to the central part
        res = rescale(roi, sf)
    else:
        if h < roi.shape[0]:    # if the thing is too big, squeez it down
            sf = h / roi.shape[0]
            res = rescale(roi, sf)
        else:                   # if the scale factor estimation is off,
            res = roi           # but the image is still smaller than the standard, just center.

    # TODO: the centering should depend on the symmetry of the word (4 cases: are, gone, to, for)
    # w = res.shape[1]
    # c = int(h / 2)
    # os_p = int(np.floor(res.shape[0] / 2))
    # os_m = int(np.ceil(res.shape[0] / 2))
    # uni = np.zeros((h, w))
    # uni[c - os_m: c + os_p, :] = res

    # Pad
    zer = np.zeros((1, res.shape[1]))
    uni = np.append(zer, res, axis=0)
    uni = np.append(uni, zer, axis=0)
    uni = uni / uni.max()

    return uni

 def scale_images(self, data, scaleRatio=-1):
     if scaleRatio <= 0:
         scaleRatio = self.scaleRatio
     if data.ndim == 3:
         for i in xrange(data.shape[0]):
             scale_dim = self.rng.random_integers(1)
             if scale_dim:
                 #Scale in first dimension
                 img_scaled = st.rescale(data[i], scale=(scaleRatio, 1))
                 scaled_shape = img_scaled.T.shape
                 I = numpy.eye(scaled_shape[0], scaled_shape[1])
                 data[i] = numpy.dot(I, img_scaled)
             else:
                 #Scale in the second dimension
                 img_scaled = st.rescale(data[i], scale=(1, scaleRatio))
                 scaled_shape = img_scaled.T.shape
                 I = numpy.eye(scaled_shape[0], scaled_shape[1])
                 data[i] = numpy.dot(img_scaled, I)
     else:
         scale_dim = self.rng.random_integers(1)
         if scale_dim:
             #Scale in first dimension
             img_scaled = st.rescale(data, scale=(scaleRatio, 1))
             scaled_shape = img_scaled.T.shape
             I = numpy.eye(scaled_shape[0], scaled_shape[1])
             data = numpy.dot(I, img_scaled)
         else:
             #Scale in the second dimension
             img_scaled = st.rescale(data, scale=(1, scaleRatio))
             scaled_shape = img_scaled.T.shape
             I = numpy.eye(scaled_shape[0], scaled_shape[1])
             data = numpy.dot(img_scaled, I)
     return data

def test_rescale_invalid_scale():
    x = np.zeros((10, 10, 3))
    with testing.raises(ValueError):
        rescale(x, (2, 2),
                multichannel=False, anti_aliasing=False, mode='constant')
    with testing.raises(ValueError):
        rescale(x, (2, 2, 2),
                multichannel=True, anti_aliasing=False, mode='constant')

def test_rescale_multichannel_defaults():
    # multichannel should always default to False as of 0.16
    x = np.zeros((8, 3), dtype=np.double)
    scaled = rescale(x, 2, order=0, anti_aliasing=False, mode='constant')
    assert_equal(scaled.shape, (16, 6))

    x = np.zeros((8, 8, 3), dtype=np.double)
    scaled = rescale(x, 2, order=0, anti_aliasing=False, mode='constant')
    assert_equal(scaled.shape, (16, 16, 6))

def pyramid(I1, I2):
	if (I1.shape[0]<=500 and I1.shape[1]<=500):
		return findBestRoll(I1, I2, range(-15, 15), range(-15, 15))
	else:
		I1R = rescale(I1, 0.5)
		I2R = rescale(I2, 0.5)
		rangeValues = pyramid(I1R, I2R)
		print(rangeValues)
		return findBestRoll(I1, I2, range(rangeValues[0]*2-2, rangeValues[0]*2+3), range(rangeValues[1]*2-2, rangeValues[1]*2+3))

def rescale_images(im1, im2, pts):
    p1, p2, p3, p4 = pts
    len1 = np.sqrt((p2[1] - p1[1])**2 + (p2[0] - p1[0])**2)
    len2 = np.sqrt((p4[1] - p3[1])**2 + (p4[0] - p3[0])**2)
    dscale = len2/len1
    if dscale < 1:
        im1 = sktr.rescale(im1, dscale)
    else:
        im2 = sktr.rescale(im2, 1./dscale)
    return im1, im2

def test_warp_clip():
    x = np.zeros((5, 5), dtype=np.double)
    x[2, 2] = 1

    outx = rescale(x, 3, order=3, clip=False)
    assert outx.min() < 0

    outx = rescale(x, 3, order=3, clip=True)
    assert_almost_equal(outx.min(), 0)
    assert_almost_equal(outx.max(), 1)

def scale_lesion(lesion, size):
    """ scale segmented lesion to uniform size """
    image = lesion[...,:3]
    alpha = lesion[...,3]

    scale = float(size) / max(*alpha.shape)
    if scale != 1.0:
        image = rescale(image, scale, order=3)
        alpha = rescale(alpha, scale, order=0)

    return np.dstack((image, alpha))

def fit_in_box(s, s_mask, s_max, t_max):
    # Resize foreground and mask so area fits in box
    y_ratio = float(t_max[0])/s_max[0]
    x_ratio = float(t_max[1])/s_max[1]
    if y_ratio > x_ratio:
        s = rescale(s, x_ratio)
        s_mask = rescale(s_mask, x_ratio)
    else:
        s = rescale(s, y_ratio)
        s_mask = rescale(s_mask, y_ratio)
    return s, s_mask

def test_warp_clip():
    x = np.zeros((5, 5), dtype=np.double)
    x[2, 2] = 1

    outx = rescale(x, 3, order=3, clip=False,
                   multichannel=False, anti_aliasing=False, mode='constant')
    assert outx.min() < 0

    outx = rescale(x, 3, order=3, clip=True,
                   multichannel=False, anti_aliasing=False, mode='constant')
    assert_almost_equal(outx.min(), 0)
    assert_almost_equal(outx.max(), 1)

def test_warp_clip():
    x = np.zeros((5, 5), dtype=np.double)
    x[2, 2] = 1

    with expected_warnings(['The default mode', 'The default multichannel']):
        outx = rescale(x, 3, order=3, clip=False)
    assert outx.min() < 0

    with expected_warnings(['The default mode', 'The default multichannel']):
        outx = rescale(x, 3, order=3, clip=True)
    assert_almost_equal(outx.min(), 0)
    assert_almost_equal(outx.max(), 1)

def random_backgrond(weights=[2, 8, 4, 3, 0.2, 0.2, 0.5], start_level=1, end_level=None):
    img = np.random.normal(0, 1, (start_level, start_level)) * weights[0]
    i = start_level - 1
    for i, w in enumerate(weights[1:], start_level):
        r = np.random.normal(0, 1, (2**i, 2**i))
        img = rescale(img, 2) + w*r

    if end_level and end_level != i:
        img = rescale(img, 2**(end_level - i))

    img = (img - img.min())
    img /= img.max()
    return img

def test_rescale_multichannel_defaults():
    # ensure multichannel=None matches the previous default behaviour

    # 2D: multichannel should default to False
    x = np.zeros((8, 3), dtype=np.double)
    with expected_warnings(['The default mode', 'The default multichannel']):
        scaled = rescale(x, 2, order=0)
    assert_equal(scaled.shape, (16, 6))

    # 3D: multichannel should default to True
    x = np.zeros((8, 8, 3), dtype=np.double)
    with expected_warnings(['The default mode', 'The default multichannel']):
        scaled = rescale(x, 2, order=0,)
    assert_equal(scaled.shape, (16, 16, 3))

def agrandissement(img):
    img_bilinear = transform.rescale(img, scale=2, mode='wrap', preserve_range=True, order=1)
    # order=1 <=> bi-linear (skimage.transform.wrap)
    img_nearNeighbor = transform.rescale(image=img, scale=2, mode='wrap', preserve_range=True, order=0)
    # order=0 <=> nearest neighbor (skimage.transform.wrap)

    print(quadratique(img_nearNeighbor, img_bilinear))  # debug

    fig = plt.figure()
    fig.add_subplot(1, 2, 1)
    plt.imshow(img_bilinear, cmap='gray')
    fig.add_subplot(1, 2, 2)
    plt.imshow(img_nearNeighbor, cmap='gray')
    plt.show()

def test_keep_range():
    image = np.linspace(0, 2, 25).reshape(5, 5)

    out = rescale(image, 2, preserve_range=False, clip=True, order=0)
    assert out.min() == 0
    assert out.max() == 2

    out = rescale(image, 2, preserve_range=True, clip=True, order=0)
    assert out.min() == 0
    assert out.max() == 2

    out = rescale(image.astype(np.uint8), 2, preserve_range=False,
                  clip=True, order=0)
    assert out.min() == 0
    assert out.max() == 2 / 255.0

def random_scaling(image, normal, size, a, b):
    scale = random.uniform(a, b);
    #print scale
    #resize images
    img_r = transform.rescale(image, scale);
    norm_r = transform.rescale(normal, scale);
    img_r, norm_r = random_crop(img_r, norm_r, size);
    #TODO modify depth : divide by scale
    #modify normals
    #for line in range(norm_r.shape[0]):
    #    for col in range(norm_r.shape[1]):
    #        norm_r[line,col,2] = norm_r[line,col,2] * scale;
    #        norm = np.linalg.norm(norm_r[line,col]);
    #        norm_r[line,col] = norm_r[line,col]/norm;
    return img_r, norm_r;

def test_iradon_sart():
    debug = False

    image = rescale(PHANTOM, 0.8)
    theta_ordered = np.linspace(0., 180., image.shape[0], endpoint=False)
    theta_missing_wedge = np.linspace(0., 150., image.shape[0], endpoint=True)
    for theta, error_factor in ((theta_ordered, 1.),
                                (theta_missing_wedge, 2.)):
        sinogram = radon(image, theta, circle=True)
        reconstructed = iradon_sart(sinogram, theta)

        if debug:
            from matplotlib import pyplot as plt
            plt.figure()
            plt.subplot(221)
            plt.imshow(image, interpolation='nearest')
            plt.subplot(222)
            plt.imshow(sinogram, interpolation='nearest')
            plt.subplot(223)
            plt.imshow(reconstructed, interpolation='nearest')
            plt.subplot(224)
            plt.imshow(reconstructed - image, interpolation='nearest')
            plt.show()

        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration) =', delta)
        assert delta < 0.02 * error_factor
        reconstructed = iradon_sart(sinogram, theta, reconstructed)
        delta = np.mean(np.abs(reconstructed - image))
        print('delta (2 iterations) =', delta)
        assert delta < 0.014 * error_factor
        reconstructed = iradon_sart(sinogram, theta, clip=(0, 1))
        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration, clip) =', delta)
        assert delta < 0.018 * error_factor

        np.random.seed(1239867)
        shifts = np.random.uniform(-3, 3, sinogram.shape[1])
        x = np.arange(sinogram.shape[0])
        sinogram_shifted = np.vstack(np.interp(x + shifts[i], x,
                                               sinogram[:, i])
                                     for i in range(sinogram.shape[1])).T
        reconstructed = iradon_sart(sinogram_shifted, theta,
                                    projection_shifts=shifts)
        if debug:
            from matplotlib import pyplot as plt
            plt.figure()
            plt.subplot(221)
            plt.imshow(image, interpolation='nearest')
            plt.subplot(222)
            plt.imshow(sinogram_shifted, interpolation='nearest')
            plt.subplot(223)
            plt.imshow(reconstructed, interpolation='nearest')
            plt.subplot(224)
            plt.imshow(reconstructed - image, interpolation='nearest')
            plt.show()

        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration, shifted sinogram) =', delta)
        assert delta < 0.022 * error_factor