def setupPatches(self, dims):
     self.dims = dims
     self.im = mpimg.imread(self.inputfile)
     self.block_shape = (self.dims[0], self.dims[1], self.im.shape[2]) #height, width
     margin=np.mod(self.im.shape,self.block_shape)
     self.im_crop = self.im[:(self.im.shape-margin)[0],:(self.im.shape-margin)[1],:(self.im.shape-margin)[2]]
     self.view = view_as_blocks(self.im_crop, self.block_shape)

def test_view_as_blocks_2D_array():

    A = np.arange(4 * 4).reshape(4, 4)
    B = view_as_blocks(A, (2, 2))
    assert_equal(B[0, 1], np.array([[2, 3],
                                   [6, 7]]))
    assert_equal(B[1, 0, 1, 1], 13)

def test_view_as_blocks_3D_array():
    A = np.arange(4 * 4 * 6).reshape(4, 4, 6)
    B = view_as_blocks(A, (1, 2, 2))
    assert_equal(B.shape, (4, 2, 3, 1, 2, 2))
    assert_equal(B[2:, 0, 2], np.array([[[[52, 53],
                                          [58, 59]]],
                                        [[[76, 77],
                                          [82, 83]]]]))

def small_spectrogram_max_pooling(spec):
    """
    Using code adapated from:
    http://scikit-image.org/docs/dev/auto_examples/plot_view_as_blocks.html
    """

    spec = force_spectrogram_length(spec, 384)
    im_norm = (spec - spec.mean()) / spec.var()

    view = view_as_blocks(im_norm, (32, 32))
    flatten_view = view.reshape(view.shape[0], view.shape[1], -1)

    return np.max(flatten_view, axis=2).flatten()

def _downsample(array, factors, sum=True):
    """Performs downsampling with integer factors.

    Parameters
    ----------
    array : ndarray
        Input n-dimensional array.
    factors: tuple
        Tuple containing downsampling factor along each axis.
    sum : bool
        If True, downsampled element is the sum of its corresponding
        constituent elements in the input array. Default is True.

    Returns
    -------
    array : ndarray
        Downsampled array with same number of dimensions as that of input
        array.

    """

    pad_size = []
    if len(factors) != array.ndim:
        raise ValueError("'factors' must have the same length "
                         "as 'array.shape'")
    else:
        for i in range(len(factors)):
            if array.shape[i] % factors[i] != 0:
                pad_size.append(factors[i] - (array.shape[i] % factors[i]))
            else:
                pad_size.append(0)

    for i in range(len(pad_size)):
        array = _pad_asymmetric_zeros(array, pad_size[i], i)

    out = view_as_blocks(array, factors)
    block_shape = out.shape

    if sum:
        for i in range(len(block_shape) // 2):
            out = out.sum(-1)
    else:
        for i in range(len(block_shape) // 2):
            out = out.mean(-1)
    return out

def frequency_max_pooling(spec, normalise=True):
    """
    Using code adapated from:
    http://scikit-image.org/docs/dev/auto_examples/plot_view_as_blocks.html
    """

    if normalise:
        im_norm = (spec - spec.mean()) / spec.var()
    else:
        im_norm = spec

    view = view_as_blocks(im_norm, (8, spec.shape[1]))
    flatten_view = view.reshape(view.shape[0], view.shape[1], -1)

    A = np.max(flatten_view, axis=2).flatten()
    B = np.var(flatten_view, axis=2).flatten()
    C = np.mean(flatten_view, axis=2).flatten()

    return np.hstack((A, B, C))

    def get_subvolume(self, box_zyx, scale=0):
        """
        Extract the subvolume, specified in new (scaled) coordinates from the
        original volume service, then scale result accordingly before returning it.
        """
        true_scale = scale + self.scale_delta
        
        if true_scale in self.original_volume_service.available_scales:
            # The original source already has the data at the necessary scale.
            return self.original_volume_service.get_subvolume( box_zyx, true_scale )

        # Start with the closest scale we've got
        base_scales = np.array(self.original_volume_service.available_scales)
        i_best = np.abs(base_scales - true_scale).argmin()
        best_base_scale = base_scales[i_best]
        
        delta_from_best = true_scale - best_base_scale

        if delta_from_best > 0:
            orig_box_zyx = box_zyx * 2**delta_from_best
            orig_data = self.original_volume_service.get_subvolume(orig_box_zyx, best_base_scale)

            if self.dtype == np.uint64:
                # Assume that uint64 means labels.
                downsampled_data, _ = downsample_labels_3d( orig_data, 2**self.scale_delta )
            else:
                downsampled_data = downsample_raw( orig_data, self.scale_delta )[-1]
            return downsampled_data
        else:
            upsample_factor = int(2**-delta_from_best)
            orig_box_zyx = downsample_box(box_zyx, np.array(3*(upsample_factor,)))
            orig_data = self.original_volume_service.get_subvolume(orig_box_zyx, best_base_scale)

            orig_shape = np.array(orig_data.shape)
            upsampled_data = np.empty( orig_shape * upsample_factor, dtype=self.dtype )
            v = view_as_blocks(upsampled_data, 3*(upsample_factor,))
            v[:] = orig_data[:,:,:,None, None, None]

            relative_box = box_zyx - upsample_factor*orig_box_zyx[0]
            requested_data = upsampled_data[box_to_slicing(*relative_box)]

            # Force contiguous so caller doesn't have to worry about it.
            return np.asarray(requested_data, order='C')

def demo_upsample_nearest():
    a = sp.misc.lena() / 1.
    a = sp.misc.lena() / 1.
    a.shape = a.shape[:2] + (1,)
    print a.shape
    #print np.tile(a, 2).shape
    #a = np.dstack((a, -a))
    N = 96
    a = np.tile(a, N)
    a[:,:,95] = -a[:,:,95]

    #r = np.tile(a, (2, 2, 1))
    #np.kron(a, np.ones((2,2,1))).shape

    # -- loop
    #a = a[:, :, 0].reshape(256, 256, 1)
    #r = np.empty((1024, 1024, 1))
    #r[0::2, 0::2] = a
    #r[0::2, 1::2] = a
    #r[1::2, 0::2] = a
    #r[1::2, 1::2] = a

    # -- block view
    r = np.empty((1024, 1024, N))
    b = view_as_blocks(r, (2, 2, 1))
    print b.shape

    a2 = a.reshape(a.shape + (1, 1, 1))
    #a[:, :, :, np.newaxis, np.newaxis, np.newaxis]
    b[:] = a2
    #import IPython; ipshell = IPython.embed; ipshell(banner1='ipshell')


    #b2 = b.swapaxes(1, 3).reshape(r.shape)
    b2 = b.transpose((0, 3, 1, 4, 2, 5)).reshape(r.shape)

    pl.matshow(b2[:,:,0])
    pl.matshow(b2[:,:,1])
    pl.matshow(b2[:,:,95])
    pl.show()

def extract_blocks(img_gray):
  patch_sizes = [128, 256]
  result = []
  for size in patch_sizes:
    blocks = view_as_blocks(img_gray, (size, size))
    for row in range(blocks.shape[0]):
      for col in range(blocks.shape[1]):
        block = blocks[row, col]
        pred = rc.predict(block)
        if pred == None:
          continue
        pred_prob = rc.predict_prob(block)[0]
        top1 = numpy.argsort(pred_prob)[-1:][0]
        top1_prob = pred_prob[top1]
        tops = numpy.argsort(pred_prob)[-5:]
        tops = tops[::-1]
        result.append((top1_prob, pred[0], row, col, size, block))
        #print "Size", size, "Prediction:", pred, "Argmax:", numpy.argmax(pred_prob), "Class:", classes[numpy.argmax(pred_prob)]
        #for idx, top in enumerate(tops):
        #  print "", idx+1, ": ", classes[top], " : ", pred_prob[top]
        #print "="*80
  return result

def test_view_as_blocks_1D_array():

    A = np.arange(10)
    B = view_as_blocks(A, (5,))
    assert_equal(B, np.array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]]))

def test_view_as_blocks_1D_array_wrong_block_shape():

    A = np.arange(10)
    view_as_blocks(A, (3,))

def test_view_as_blocks_wrong_block_dimension():

    A = np.arange(10)
    view_as_blocks(A, (2, 2))

def test_view_as_blocks_block_too_large():

    A = np.arange(10)
    view_as_blocks(A, (11,))

def test_view_as_blocks_negative_shape():

    A = np.arange(10)
    view_as_blocks(A, (-2,))

def main():
    print "Load data"

    sourceHeight = 'height_masked_final.asc'
    sourceBiomass = 'ndvi_masked_final.asc'

    sourceCoverageModel = 'testCov.txt'
    sourceHeightModel = 'heightModel.txt'
    sourceLCC = "LCC.asc"

    heightGrid = numpy.loadtxt(sourceHeight, skiprows=6)
    heightGrid = heightGrid[500:1050, 300:750]
    #    heightGrid = heightGrid[640:685,455:500]
    #    heightGrid = heightGrid[860:890,600:700]
    rgb = (heightGrid - numpy.min(heightGrid)) / (numpy.max(heightGrid) - numpy.min(heightGrid))
    rgb *= 255
    heightRGBA = numpy.zeros((heightGrid.shape[0], heightGrid.shape[1], 3), dtype=numpy.uint8)
    heightRGBA[:, :, 0:3] = rgb[:, :, numpy.newaxis]
    # misc.imsave('heightMap_Paulinapolder.png',heightRGBA)
    ndviGrid = numpy.loadtxt(sourceBiomass, skiprows=6)
    ndviGrid = ndviGrid[500:1050, 300:750]
    #    ndviGrid = ndviGrid[640:685,455:500]
    #    ndviGrid = ndviGrid[860:890,600:700]
    rgb = (ndviGrid - numpy.min(ndviGrid)) / (numpy.max(ndviGrid) - numpy.min(ndviGrid))
    rgb *= 255
    ndviRGBA = numpy.zeros((ndviGrid.shape[0], ndviGrid.shape[1], 3), dtype=numpy.uint8)
    ndviRGBA[:, :, 0:3] = rgb[:, :, numpy.newaxis]
    misc.imsave('ndviMap_Paulinapolder.png', ndviRGBA)
    lccGrid = numpy.loadtxt(sourceLCC, skiprows=6)
    lccGrid = lccGrid[500:1050, 300:750]
    heightModelGrid = numpy.loadtxt(sourceHeightModel)
    coverageModelGrid = numpy.loadtxt(sourceCoverageModel)
    PaulinaPolder = False
    NDVI = True
    LCC = False
    if PaulinaPolder:
        if NDVI and LCC:
            vegetationMask = ndviGrid > 0
        elif NDVI:
            vegetationMask = ndviGrid > 0.08  # 0.02 demo
        #             figure()
        #             tempveg = numpy.zeros((ndviGrid.shape[0],ndviGrid.shape[1],3))
        #             tempveg[:,:,:] = (ndviGrid[:,:,numpy.newaxis]+1) / 2.0;
        #             imshow(tempveg)
        #             show()
        #             figure()
        #             tempveg = numpy.zeros((vegetationMask.shape[0],vegetationMask.shape[1],3))
        #             tempveg[:,:,:] = vegetationMask[:,:,numpy.newaxis];
        #             imshow(tempveg)
        #             show()
        #             figure()
        #             tempveg = numpy.zeros((vegetationMask.shape[0],vegetationMask.shape[1],3))
        #             tempveg[:,:,:] = heightGrid[:,:,numpy.newaxis] / numpy.max(heightGrid);
        #             imshow(tempveg)
        #             show()
        elif LCC:
            vegetationMask = lccGrid > 0
        heightValues = heightGrid[vegetationMask]
        baseValues = ndviGrid[vegetationMask]
        lccValues = lccGrid[vegetationMask]
        lengthX, lengthY = heightGrid.shape
        nTypes = 7
        area = "NDVI"

        lXTemp = lengthX
        lYTemp = lengthY
        if lengthX % 2 == 1:
            lXTemp += 1
        if lengthY % 2 == 1:
            lYTemp += 1
        vegetationMaskExtended = np.zeros((lXTemp, lYTemp), dtype=bool)
        vegetationMaskExtended[0:lengthX, 0:lengthY] = vegetationMask
        res = 2.0  # 2.0
        wangGridLengthX = np.ceil(lengthX / res)
        wangGridLengthY = np.ceil(lengthY / res)
        xWangIndices, yWangIndices = numpy.indices((wangGridLengthX, wangGridLengthY))
        blocks = view_as_blocks(vegetationMaskExtended, block_shape=(int(res), int(res)))
        blocks = blocks.reshape(wangGridLengthX, wangGridLengthY, res * res)
        blocks_summed = np.sum(blocks, axis=2)
        wangVegetationMask = blocks_summed > 0
        print "wvm", wangVegetationMask.shape
        print "vm", vegetationMask.shape
    else:
        vegetationMask = coverageModelGrid > 0
        heightValues = heightModelGrid[vegetationMask] * 1
        baseValues = coverageModelGrid[vegetationMask] * .01
        lengthX, lengthY = heightModelGrid.shape
        nTypes = 3
        heightValues += numpy.fabs(numpy.min(heightValues))
        minHeight = numpy.min(heightValues)
        maxHeight = numpy.max(heightValues)
        heightValues = (heightValues - minHeight) / (maxHeight - minHeight)
        rgb = (heightModelGrid - numpy.min(heightModelGrid)) / (numpy.max(heightModelGrid) - numpy.min(heightModelGrid))
        rgb *= 255
        heightRGBA = numpy.zeros((heightModelGrid.shape[0], heightModelGrid.shape[1], 3), dtype=numpy.uint8)
        heightRGBA[:, :, 0:3] = rgb[:, :, numpy.newaxis]
        # misc.imsave('heightMap_Ecomodel.png',heightRGBA)

        rgb = (coverageModelGrid - numpy.min(coverageModelGrid)) / (
        numpy.max(coverageModelGrid) - numpy.min(coverageModelGrid))
#.........这里部分代码省略.........

import matplotlib.cm as cm

from skimage import data
from skimage import color
from skimage.util.shape import view_as_blocks


# -- get `astronaut` from skimage.data in grayscale
l = color.rgb2gray(data.astronaut())

# -- size of blocks
block_shape = (4, 4)

# -- see `astronaut` as a matrix of blocks (of shape
#    `block_shape`)
view = view_as_blocks(l, block_shape)

# -- collapse the last two dimensions in one
flatten_view = view.reshape(view.shape[0], view.shape[1], -1)

# -- resampling `astronaut` by taking either the `mean`,
#    the `max` or the `median` value of each blocks.
mean_view = np.mean(flatten_view, axis=2)
max_view = np.max(flatten_view, axis=2)
median_view = np.median(flatten_view, axis=2)

# -- display resampled images
fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharex=True, sharey=True)
ax0, ax1, ax2, ax3 = axes.ravel()

ax0.set_title("Original rescaled with\n spline interpolation (order=3)")

def test_view_as_blocks_wrong_block_dimension():
    A = np.arange(10)
    with testing.raises(ValueError):
        view_as_blocks(A, (2, 2))

def test_view_as_blocks_1D_array_wrong_block_shape():
    A = np.arange(10)
    with testing.raises(ValueError):
        view_as_blocks(A, (3,))

def test_view_as_blocks_block_too_large():
    A = np.arange(10)
    with testing.raises(ValueError):
        view_as_blocks(A, (11,))

def test_view_as_blocks_block_not_a_tuple():

    A = np.arange(10)
    view_as_blocks(A, [5])