Python scipy.ndimage() Examples

def predict_multiscale(net, image, tile_size, scales, classes, flip_evaluation, recurrence):
    """
    Predict an image by looking at it with different scales.
        We choose the "predict_whole_img" for the image with less than the original input size,
        for the input of larger size, we would choose the cropping method to ensure that GPU memory is enough.
    """
    image = image.data
    N_, C_, H_, W_ = image.shape
    full_probs = np.zeros((H_, W_, classes))  
    for scale in scales:
        scale = float(scale)
        print("Predicting image scaled by %f" % scale)
        scale_image = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
        scaled_probs = predict_whole(net, scale_image, tile_size, recurrence)
        if flip_evaluation == True:
            flip_scaled_probs = predict_whole(net, scale_image[:,:,:,::-1].copy(), tile_size, recurrence)
            scaled_probs = 0.5 * (scaled_probs + flip_scaled_probs[:,::-1,:])
        full_probs += scaled_probs
    full_probs /= len(scales)
    return full_probs 

def scipy_conv_c01b(self, images, filters):
        """
        Emulate c01b convolution with scipy
        """
        assert images.ndim == 4
        assert filters.ndim == 4
        in_chans, rows, cols, bs = images.shape
        in_chans_, rows_, cols_, out_chans = filters.shape
        assert in_chans_ == in_chans

        out_bc01 = [
            [sum(scipy.ndimage.filters.convolve(images[c, :, :, b],
                                                filters[c, ::-1, ::-1, i])
                 for c in xrange(in_chans))
             for i in xrange(out_chans)]
            for b in xrange(bs)]
        out_c01b = numpy.array(out_bc01).transpose(1, 2, 3, 0)
        return out_c01b 

def predict_multiscale(net, image, tile_size, scales, classes, flip_evaluation, recurrence):
    """
    Predict an image by looking at it with different scales.
        We choose the "predict_whole_img" for the image with less than the original input size,
        for the input of larger size, we would choose the cropping method to ensure that GPU memory is enough.
    """
    image = image.data
    N_, C_, H_, W_ = image.shape
    full_probs = np.zeros((H_, W_, classes))  
    for scale in scales:
        scale = float(scale)
        #print("Predicting image scaled by %f" % scale)
        scale_image = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
        scaled_probs = predict_whole(net, scale_image, tile_size, recurrence)
        if flip_evaluation == True:
            flip_scaled_probs = predict_whole(net, scale_image[:,:,:,::-1].copy(), tile_size, recurrence)
            scaled_probs = 0.5 * (scaled_probs + flip_scaled_probs[:,::-1,:])
        full_probs += scaled_probs
    full_probs /= len(scales)
    return full_probs 

def test_binary_ops(funcname,
                    input,
                    structure,
                    origin):
    da_func = getattr(da_ndm, funcname)
    sp_func = getattr(spnd, funcname)

    da_result = da_func(
        input,
        structure=structure,
        origin=origin
    )

    sp_result = sp_func(
        input,
        structure=structure,
        origin=origin
    )

    dau.assert_eq(sp_result, da_result) 

def predict_multiscale(net, image, tile_size, scales, classes, flip_evaluation, recurrence):
    """
    Predict an image by looking at it with different scales.
        We choose the "predict_whole_img" for the image with less than the original input size,
        for the input of larger size, we would choose the cropping method to ensure that GPU memory is enough.
    """
    image = image.data
    N_, C_, H_, W_ = image.shape
    full_probs = np.zeros((H_, W_, classes))  
    for scale in scales:
        scale = float(scale)
        print("Predicting image scaled by %f" % scale)
        scale_image = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
        scaled_probs = predict_whole(net, scale_image, tile_size, recurrence)
        if flip_evaluation == True:
            flip_scaled_probs = predict_whole(net, scale_image[:,:,:,::-1].copy(), tile_size, recurrence)
            scaled_probs = 0.5 * (scaled_probs + flip_scaled_probs[:,::-1,:])
        full_probs += scaled_probs
    full_probs /= len(scales)
    return full_probs 

def test_binary_ops_iter(funcname,
                         input,
                         structure,
                         iterations,
                         origin):
    da_func = getattr(da_ndm, funcname)
    sp_func = getattr(spnd, funcname)

    da_result = da_func(
        input,
        structure=structure,
        iterations=iterations,
        origin=origin
    )

    sp_result = sp_func(
        input,
        structure=structure,
        iterations=iterations,
        origin=origin
    )

    dau.assert_eq(sp_result, da_result) 

def applyFilter():
    global imageRaw

    # Repeat array four times for Lepton2
    if leptonVersion == 0:
        array2d = leptonValues.reshape(60, 80)
        array2dBig = array2d.repeat(4, axis=0).repeat(4, axis=1)
        imageRaw = numpy.transpose(array2dBig)

    # Repeat array two times for Lepton3
    else:
        array2d = leptonValues.reshape(120, 160)
        array2dBig = array2d.repeat(2, axis=0).repeat(2, axis=1)
        imageRaw = numpy.transpose(array2dBig)

    # Apply the gaussian blur filter
    if filterType == 1:
        imageRaw = scipy.ndimage.filters.gaussian_filter(imageRaw, 1.33, mode='nearest')
    # Apply box blur filter
    if filterType == 2:
        imageRaw = scipy.ndimage.filters.uniform_filter(imageRaw, 3)


# Converts the Lepton raw values to RGB colors 

def bwdist(bwvol):
    """
    positive distance transform from positive entries in logical image

    Parameters
    ----------
    bwvol : nd array
        The logical volume

    Returns
    -------
    possdtrf : nd array
        the positive distance transform

    See Also
    --------
    bw2sdtrf
    """

    # reverse volume to run scipy function
    revbwvol = np.logical_not(bwvol)

    # get distance
    return scipy.ndimage.morphology.distance_transform_edt(revbwvol) 

def shear(x, intensity=0.1, is_random=False, row_index=0, col_index=1, channel_index=2,
                 fill_mode='nearest', cval=0.):
    """Shear an image randomly or non-randomly.

    Parameters
    -----------
    x : numpy array
        An image with dimension of [row, col, channel] (default).
    intensity : float
        Percentage of shear, usually -0.5 ~ 0.5 (is_random==True), 0 ~ 0.5 (is_random==False),
        you can have a quick try by shear(X, 1).
    is_random : boolean, default False
        If True, randomly shear.
    row_index, col_index, channel_index : int
        Index of row, col and channel, default (0, 1, 2), for theano (1, 2, 0).
    fill_mode : string
        Method to fill missing pixel, default ‘nearest’, more options ‘constant’, ‘reflect’ or ‘wrap’.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    cval : scalar, optional
        Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    """
    if is_random:
        shear = np.random.uniform(-intensity, intensity)
    else:
        shear = intensity
    shear_matrix = np.array([[1, -np.sin(shear), 0],
                             [0, np.cos(shear), 0],
                             [0, 0, 1]])

    h, w = x.shape[row_index], x.shape[col_index]
    transform_matrix = transform_matrix_offset_center(shear_matrix, h, w)
    x = apply_transform(x, transform_matrix, channel_index, fill_mode, cval)
    return x 

def __call__(self, image):
        angle = self.random_state.uniform(
            self.angle_range[0], self.angle_range[1])
        if isinstance(image, np.ndarray):
            mi, ma = image.min(), image.max()
            image = scipy.ndimage.interpolation.rotate(
                image, angle, reshape=False, axes=self.axes, mode=self.mode)
            return np.clip(image, mi, ma)
        elif isinstance(image, Image.Image):
            return image.rotate(angle)
        else:
            raise Exception('unsupported type') 

def derive(self, samples):
		samples = samples[0]
		uuids = [s['uuid'] for s in samples]
		image_paths = [os.path.join(self.path, 'images', '{}.png'.format(u)) for u in uuids]

		derived = numpy.array([scipy.ndimage.imread(p).transpose() for p in image_paths])
		derived = derived.astype(numpy.float32)
		derived /= 255
		return derived 

def __call__(self, image):
        angle = self.random_state.uniform(
            self.angle_range[0], self.angle_range[1])
        if isinstance(image, np.ndarray):
            mi, ma = image.min(), image.max()
            image = scipy.ndimage.interpolation.rotate(
                image, angle, reshape=False, axes=self.axes, mode=self.mode)
            return np.clip(image, mi, ma)
        elif isinstance(image, Image.Image):
            return image.rotate(angle)
        else:
            raise Exception('unsupported type') 

def __call__(self, image):
        if isinstance(image, np.ndarray):
            return scipy.ndimage.interpolation.zoom(image, self.zoom, order=1)
        elif isinstance(image, Image.Image):
            return image.resize(self.size, Image.BILINEAR)
        else:
            raise Exception('unsupported type') 

def test_binary_ops_expanded(funcname,
                             input,
                             structure,
                             iterations,
                             mask,
                             border_value,
                             origin,
                             brute_force):
    da_func = getattr(da_ndm, funcname)
    sp_func = getattr(spnd, funcname)

    da_result = da_func(
        input,
        structure=structure,
        iterations=iterations,
        mask=mask,
        border_value=border_value,
        origin=origin,
        brute_force=brute_force
    )

    sp_result = sp_func(
        input,
        structure=structure,
        iterations=iterations,
        mask=mask,
        border_value=border_value,
        origin=origin,
        brute_force=brute_force
    )

    dau.assert_eq(sp_result, da_result) 

def rotation(x, rg=20, is_random=False, row_index=0, col_index=1, channel_index=2,
                    fill_mode='nearest', cval=0.):
    """Rotate an image randomly or non-randomly.

    Parameters
    -----------
    x : numpy array
        An image with dimension of [row, col, channel] (default).
    rg : int or float
        Degree to rotate, usually 0 ~ 180.
    is_random : boolean, default False
        If True, randomly rotate.
    row_index, col_index, channel_index : int
        Index of row, col and channel, default (0, 1, 2), for theano (1, 2, 0).
    fill_mode : string
        Method to fill missing pixel, default ‘nearest’, more options ‘constant’, ‘reflect’ or ‘wrap’

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    cval : scalar, optional
        Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_

    Examples
    ---------
    >>> x --> [row, col, 1] greyscale
    >>> x = rotation(x, rg=40, is_random=False)
    >>> tl.visualize.frame(x[:,:,0], second=0.01, saveable=True, name='temp',cmap='gray')
    """
    if is_random:
        theta = np.pi / 180 * np.random.uniform(-rg, rg)
    else:
        theta = np.pi /180 * rg
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta), 0],
                                [np.sin(theta), np.cos(theta), 0],
                                [0, 0, 1]])

    h, w = x.shape[row_index], x.shape[col_index]
    transform_matrix = transform_matrix_offset_center(rotation_matrix, h, w)
    x = apply_transform(x, transform_matrix, channel_index, fill_mode, cval)
    return x 

def shift(x, wrg=0.1, hrg=0.1, is_random=False, row_index=0, col_index=1, channel_index=2,
                 fill_mode='nearest', cval=0.):
    """Shift an image randomly or non-randomly.

    Parameters
    -----------
    x : numpy array
        An image with dimension of [row, col, channel] (default).
    wrg : float
        Percentage of shift in axis x, usually -0.25 ~ 0.25.
    hrg : float
        Percentage of shift in axis y, usually -0.25 ~ 0.25.
    is_random : boolean, default False
        If True, randomly shift.
    row_index, col_index, channel_index : int
        Index of row, col and channel, default (0, 1, 2), for theano (1, 2, 0).
    fill_mode : string
        Method to fill missing pixel, default ‘nearest’, more options ‘constant’, ‘reflect’ or ‘wrap’.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    cval : scalar, optional
        Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    """
    h, w = x.shape[row_index], x.shape[col_index]
    if is_random:
        tx = np.random.uniform(-hrg, hrg) * h
        ty = np.random.uniform(-wrg, wrg) * w
    else:
        tx, ty = hrg * h, wrg * w
    translation_matrix = np.array([[1, 0, tx],
                                   [0, 1, ty],
                                   [0, 0, 1]])

    transform_matrix = translation_matrix  # no need to do offset
    x = apply_transform(x, transform_matrix, channel_index, fill_mode, cval)
    return x 

def shear(x, intensity=0.1, is_random=False, row_index=0, col_index=1, channel_index=2,
                 fill_mode='nearest', cval=0.):
    """Shear an image randomly or non-randomly.

    Parameters
    -----------
    x : numpy array
        An image with dimension of [row, col, channel] (default).
    intensity : float
        Percentage of shear, usually -0.5 ~ 0.5 (is_random==True), 0 ~ 0.5 (is_random==False),
        you can have a quick try by shear(X, 1).
    is_random : boolean, default False
        If True, randomly shear.
    row_index, col_index, channel_index : int
        Index of row, col and channel, default (0, 1, 2), for theano (1, 2, 0).
    fill_mode : string
        Method to fill missing pixel, default ‘nearest’, more options ‘constant’, ‘reflect’ or ‘wrap’.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    cval : scalar, optional
        Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0.

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    """
    if is_random:
        shear = np.random.uniform(-intensity, intensity)
    else:
        shear = intensity
    shear_matrix = np.array([[1, -np.sin(shear), 0],
                             [0, np.cos(shear), 0],
                             [0, 0, 1]])

    h, w = x.shape[row_index], x.shape[col_index]
    transform_matrix = transform_matrix_offset_center(shear_matrix, h, w)
    x = apply_transform(x, transform_matrix, channel_index, fill_mode, cval)
    return x 

def zoom(self, zoom_factor, inplace=False, order=3):
        """
        Resize the image to a new pixel size (but keep the real size) by pixel interpolation.

        Parameters
        ----------
        zoom_factor : float
            > 1: up sampling
            < 1: down sampling
        order : int
            The spline interpolation order to use. (default: 3). Use 0 for binary or very sharp images.
        inplace : bool
            create a new image?
        """
        from scipy.ndimage.interpolation import zoom
        if not inplace:        
            new = copy.deepcopy(self)
            new.pixels = zoom(new.pixels, zoom_factor, order=order)
            new.size['pixels']['x'] = new.pixels.shape[1]
            new.size['pixels']['y'] = new.pixels.shape[0]
            return new
        else:
            self.pixels = zoom(self.pixels, zoom_factor, order=order)
            self.size['pixels']['x'] = self.pixels.shape[1]
            self.size['pixels']['y'] = self.pixels.shape[0]
            return self

# Note: The following functions are not part of the SPM_image class.
# All following functions are performed on numpy arrays 

def load_images_torcs_4():
    import glob
    filepath = 'datasets/torcs_4/training_set/*.png'
    filenames = glob.glob(filepath)

    sample = scipy.ndimage.imread(filenames[0])
    num_images = len(filenames)
    images = np.zeros((num_images, sample.shape[0], sample.shape[1]), dtype=np.uint8)

    for i in range(num_images):
        images[i] = scipy.ndimage.imread(filenames[i])

    images = images.reshape(len(images), 1, 64, 64)

    return images 

def predict_whole_img(net, image, classes, method, scale):
    """
         Predict the whole image w/o using multiple crops.
         The scale specify whether rescale the input image before predicting the results.
    """
    N_, C_, H_, W_ = image.shape
    if torch_ver == '0.4':
        interp = nn.Upsample(size=(H_, W_), mode='bilinear', align_corners=True)
    else:
        interp = nn.Upsample(size=(H_, W_), mode='bilinear')
    if scale != 1:
        scaled_img = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
    else:
        scaled_img = image
    
    full_prediction_ = net(Variable(torch.from_numpy(scaled_img), volatile=True).cuda(), )
    if 'dsn' in method or 'center' in method or 'fuse' in method:
        full_prediction = full_prediction_[-1]
    else:
        full_prediction = full_prediction_

    if torch_ver == '0.4':
        full_prediction = F.upsample(input=full_prediction, size=(H_, W_), mode='bilinear', align_corners=True)
    else:
        full_prediction = F.upsample(input=full_prediction, size=(H_, W_), mode='bilinear')
    result = full_prediction.cpu().data.numpy().transpose(0,2,3,1)
    return result 

def predict_whole_img_w_label(net, image, classes, method, scale, label):
    """
         Predict the whole image w/o using multiple crops.
         The scale specify whether rescale the input image before predicting the results.
    """
    N_, C_, H_, W_ = image.shape
    if torch_ver == '0.4':
        interp = nn.Upsample(size=(H_, W_), mode='bilinear', align_corners=True)
    else:
        interp = nn.Upsample(size=(H_, W_), mode='bilinear')

#     bug
#     if scale > 1:
#         scaled_img = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
#     else:
#         scaled_img = image

    scaled_img = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
    
    full_prediction_ = net(Variable(torch.from_numpy(scaled_img), volatile=True).cuda(), label)
    if 'dsn' in method or 'center' in method or 'fuse' in method:
        full_prediction = full_prediction_[-1]
    else:
        full_prediction = full_prediction_

    full_prediction = F.upsample(input=full_prediction, size=(H_, W_), mode='bilinear', align_corners=True)
    result = full_prediction.cpu().data.numpy().transpose(0,2,3,1)
    return result 

def predict_whole_img(net, image, classes, method, scale):
    """
         Predict the whole image w/o using multiple crops.
         The scale specify whether rescale the input image before predicting the results.
    """
    N_, C_, H_, W_ = image.shape

    if torch_ver == '0.4':
        interp = nn.Upsample(size=(H_, W_), mode='bilinear', align_corners=True)
    else:
        interp = nn.Upsample(size=(H_, W_), mode='bilinear')

#     bug
#     if scale > 1:
#         scaled_img = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
#     else:
#         scaled_img = image

    scaled_img = ndimage.zoom(image, (1.0, 1.0, scale, scale), order=1, prefilter=False)
    
    full_prediction_ = net(Variable(torch.from_numpy(scaled_img), volatile=True).cuda(), )
    if 'dsn' in method:
        full_prediction = full_prediction_[-1]
    else:
        full_prediction = full_prediction_

    if torch_ver == '0.4':
        full_prediction = F.upsample(input=full_prediction, size=(H_, W_), mode='bilinear', align_corners=True)
    else:
        full_prediction = F.upsample(input=full_prediction, size=(H_, W_), mode='bilinear')
    
    result = full_prediction.cpu().data.numpy().transpose(0,2,3,1)

    # result = interp(full_prediction).cpu().data.numpy().transpose(0,2,3,1)
    return result 

def lineProfile(zValue,originalPoint,endPoint): 
    num=equalDistance+1
    
    x0=originalPoint[:,0].reshape(-1,1)
    x1=endPoint[:,:,0]
    # print(x0.dtype)
    # print(x1.dtype)
    x=np.linspace(x0, x1, num,dtype=np.int) #可以不用修改数组类型。出于内存优化考虑,会加快后续np.stack计算速度
    del x0,x1
    # fp[:]=np.linspace(x0, x1, num,dtype=np.int)[:]
    
    
    y0=originalPoint[:,1].reshape(-1,1)  
    y1=endPoint[:,:,1]
    y=np.linspace(y0, y1, num,dtype=np.int)
    del y0,y1
    
    xStack=np.stack(x,axis=-1)
    # xStackSplit=np.array_split(xStack,saveN) #为减缓单次数组计算量，切分数组为saveN个单个数组，切分量越大，单个数组越小，计算所占用内存越小
    # del xStack #清空该变量，释放内存.关于numpy array占据内存的大小，查找确认
    # print("xStack finished!")
    
    yStack=np.stack(y,axis=-1)
    # yStackSplit=np.array_split(yStack,saveN)
    # del yStack
    del x,y

    zi = scipy.ndimage.map_coordinates(zValue,[xStack,yStack],cval=0,mode="nearest",order=0) #根据数组索引值，提取实际值
    del xStack,yStack
    

    # #记录每次存储的数据路径为.txt文件，最后通过该文件读取所有数据，组合与写入raster格式文件
    # with open(os.path.join(saveFp,"SVFValue_filenames.txt"), "a") as text_file:
    #     text_file.write("%s\n"%SVFValueFn)
    return zi

#read raster blocks 分批读取较大栅格数据，并计算 

def lineProfile(zValue,originalPoint,endPoint): 
    num=equalDistance+1
    
    x0=originalPoint[:,0].reshape(-1,1)
    x1=endPoint[:,:,0]
    # print(x0.dtype)
    # print(x1.dtype)
    x=np.linspace(x0, x1, num,dtype=np.int) #可以不用修改数组类型。出于内存优化考虑,会加快后续np.stack计算速度
    del x0,x1
    # fp[:]=np.linspace(x0, x1, num,dtype=np.int)[:]
    
    
    y0=originalPoint[:,1].reshape(-1,1)  
    y1=endPoint[:,:,1]
    y=np.linspace(y0, y1, num,dtype=np.int)
    del y0,y1
    
    xStack=np.stack(x,axis=-1)
    # xStackSplit=np.array_split(xStack,saveN) #为减缓单次数组计算量，切分数组为saveN个单个数组，切分量越大，单个数组越小，计算所占用内存越小
    # del xStack #清空该变量，释放内存.关于numpy array占据内存的大小，查找确认
    # print("xStack finished!")
    
    yStack=np.stack(y,axis=-1)
    # yStackSplit=np.array_split(yStack,saveN)
    # del yStack
    del x,y

    zi = scipy.ndimage.map_coordinates(zValue,[xStack,yStack],cval=0,mode="nearest",order=0) #根据数组索引值，提取实际值
    del xStack,yStack
    

    # #记录每次存储的数据路径为.txt文件，最后通过该文件读取所有数据，组合与写入raster格式文件
    # with open(os.path.join(saveFp,"SVFValue_filenames.txt"), "a") as text_file:
    #     text_file.write("%s\n"%SVFValueFn)
    return zi

#read raster blocks 分批读取较大栅格数据，并计算 

def make_mask_array(drz_image):
    """
    Creates _msk.fits mask file that contains pixel values of 1 outside the drizzled image footprint and pixel values
    of 0 inside the footprint. This file is used by subroutine hla_nexp_flags().

    Parameters
    ----------
    drz_image : string
        drizzled image filename

    Returns
    -------
    mask : numpy.ndarray object
        mask array
    """
    mask = fits.open(drz_image)[1].data != 0
    dilate = scipy.ndimage.morphology.binary_dilation
    erode = scipy.ndimage.morphology.binary_erosion
    kernel1 = numpy.ones((25, 25), dtype=int)
    kernel2 = numpy.ones((31, 31), dtype=int)
    # add padding around the edge so pixels close to image boundary are correct
    padding = 13
    bigmask = numpy.pad(mask, padding, 'constant')
    # strip the padding back off after creating mask
    mask = (erode(dilate(bigmask, kernel1), kernel2) == 0)[padding:-padding, padding:-padding]
    mask = mask.astype(numpy.int16)
    return mask


# ====================================================================================================================== 

def findBrightestSpot(image, minpct=99.99):
    """
    Find the sun position (in pixels, in the current projection) using the image.
    """
    if isinstance(image, envmap.EnvironmentMap):
        image = image.data

    # Gaussian filter
    filteredimg = scipy.ndimage.filters.gaussian_filter(image, (5, 5, 0))

    # Intensity image
    if filteredimg.ndim == 2 or filteredimg.shape[2] > 1:
        intensityimg = np.dot( filteredimg[:,:,:3], [.299, .587, .114] )
    else:
        intensityimg = filteredimg
    intensityimg[~np.isfinite(intensityimg)] = 0

    # Look for the value a the *minpct* percentage and threshold at this value
    # We do not take into account the pixels with a value of 0
    minval = np.percentile(intensityimg[intensityimg > 0], minpct)
    thresholdmap = intensityimg >= minval

    # Label the regions in the thresholded image
    labelarray, n = scipy.ndimage.measurements.label(thresholdmap, np.ones((3, 3), dtype="bool8"))

    # Find the size of each of them
    funcsize = lambda x: x.size
    patchsizes = scipy.ndimage.measurements.labeled_comprehension(intensityimg,
                                                                  labelarray,
                                                                  index=np.arange(1, n+1),
                                                                  func=funcsize,
                                                                  out_dtype=np.uint32,
                                                                  default=0.0)

    # Find the biggest one (we must add 1 because the label 0 is the background)
    biggestPatchIdx = np.argmax(patchsizes) + 1

    # Obtain the center of mass of the said biggest patch (we suppose that it is the sun)
    centerpos = scipy.ndimage.measurements.center_of_mass(intensityimg, labelarray, biggestPatchIdx)

    return centerpos 

def apply_transform(x,
                    transform_matrix,
                    channel_axis=0,
                    fill_mode='nearest',
                    cval=0.):
  """Apply the image transformation specified by a matrix.

  Arguments:
      x: 2D numpy array, single image.
      transform_matrix: Numpy array specifying the geometric transformation.
      channel_axis: Index of axis for channels in the input tensor.
      fill_mode: Points outside the boundaries of the input
          are filled according to the given mode
          (one of `{'constant', 'nearest', 'reflect', 'wrap'}`).
      cval: Value used for points outside the boundaries
          of the input if `mode='constant'`.

  Returns:
      The transformed version of the input.
  """
  x = np.rollaxis(x, channel_axis, 0)
  final_affine_matrix = transform_matrix[:2, :2]
  final_offset = transform_matrix[:2, 2]
  channel_images = [
      ndi.interpolation.affine_transform(
          x_channel,
          final_affine_matrix,
          final_offset,
          order=0,
          mode=fill_mode,
          cval=cval) for x_channel in x
  ]
  x = np.stack(channel_images, axis=0)
  x = np.rollaxis(x, 0, channel_axis + 1)
  return x 

def rotation(x, rg=20, is_random=False, row_index=0, col_index=1, channel_index=2,
                    fill_mode='nearest', cval=0.):
    """Rotate an image randomly or non-randomly.

    Parameters
    -----------
    x : numpy array
        An image with dimension of [row, col, channel] (default).
    rg : int or float
        Degree to rotate, usually 0 ~ 180.
    is_random : boolean, default False
        If True, randomly rotate.
    row_index, col_index, channel_index : int
        Index of row, col and channel, default (0, 1, 2), for theano (1, 2, 0).
    fill_mode : string
        Method to fill missing pixel, default ‘nearest’, more options ‘constant’, ‘reflect’ or ‘wrap’

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_
    cval : scalar, optional
        Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0

        - `scipy ndimage affine_transform <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.affine_transform.html>`_

    Examples
    ---------
    >>> x --> [row, col, 1] greyscale
    >>> x = rotation(x, rg=40, is_random=False)
    >>> tl.visualize.frame(x[:,:,0], second=0.01, saveable=True, name='temp',cmap='gray')
    """
    if is_random:
        theta = np.pi / 180 * np.random.uniform(-rg, rg)
    else:
        theta = np.pi /180 * rg
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta), 0],
                                [np.sin(theta), np.cos(theta), 0],
                                [0, 0, 1]])

    h, w = x.shape[row_index], x.shape[col_index]
    transform_matrix = transform_matrix_offset_center(rotation_matrix, h, w)
    x = apply_transform(x, transform_matrix, channel_index, fill_mode, cval)
    return x 

def __call__(self, image):
        if isinstance(image, np.ndarray):
            return scipy.ndimage.interpolation.zoom(image, self.zoom)
        elif isinstance(image, Image.Image):
            return image.resize(self.size, Image.BILINEAR)
        else:
            raise Exception('unsupported type') 

def transform_array(self, X, y, w):
    """Transform the data in a set of (X, y, w) arrays."""
    from PIL import Image
    images = [scipy.ndimage.imread(x, mode='RGB') for x in X]
    images = [Image.fromarray(x).resize(self.size) for x in images]
    return np.array(images), y, w