<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def bilateral_filter(image, mask, sigma_spatial, sigma_range, sampling_spatial = None, sampling_range = None):
 """Bilateral filter of an image image - image to be bilaterally filtered mask - mask of significant points in image sigma_spatial - standard deviation of the spatial Gaussian sigma_range - standard deviation of the range Gaussian sampling_spatial - amt to reduce image array extents when sampling default is 1/2 sigma_spatial sampling_range - amt to reduce the range of values when sampling default is 1/2 sigma_range The bilateral filter is described by the following equation: sum(Fs(||p - q||)Fr(|Ip - Iq|)Iq) / sum(Fs(||p-q||)Fr(|Ip - Iq)) where the sum is over all points in the kernel p is all coordinates in the image q is the coordinates as perturbed by the mask Ip is the intensity at p Iq is the intensity at q Fs is the spatial convolution function, for us a Gaussian that falls off as the distance between falls off Fr is the "range" distance which falls off as the difference in intensity increases. 1 / sum(Fs(||p-q||)Fr(|Ip - Iq)) is the weighting for point p """

    # The algorithm is taken largely from code by Jiawen Chen which miraculously
    # extends to the masked case:
    # http://groups.csail.mit.edu/graphics/bilagrid/bilagrid_web.pdf
    #
    # Form a 3-d array whose extent is reduced in the i,j directions
    # by the spatial sampling parameter and whose extent is reduced in the
    # z (image intensity) direction by the range sampling parameter.

    # Scatter each significant pixel in the image into the nearest downsampled
    # array address where the pixel's i,j coordinate gives the corresponding
    # i and j in the matrix and the intensity value gives the corresponding z
    # in the array.

    # Count the # of values entered into each 3-d array element to form a
    # weight.

    # Similarly convolve the downsampled value and weight arrays with a 3-d
    # Gaussian kernel whose i and j Gaussian is the sigma_spatial and whose
    # z is the sigma_range.
    #
    # Divide the value by the weight to scale each z value appropriately
    #
    # Linearly interpolate using an i x j x 3 array where [:,:,0] is the
    # i coordinate in the downsampled array, [:,:,1] is the j coordinate
    # and [:,:,2] is the unrounded index of the z-slot
    #
    # One difference is that I don't pad the intermediate arrays. The
    # weights bleed off the edges of the intermediate arrays and this
    # accounts for the ring of zero values used at the border bleeding
    # back into the intermediate arrays during convolution
    #

    if sampling_spatial is None:
        sampling_spatial = sigma_spatial / 2.0
    if sampling_range is None:
        sampling_range = sigma_range / 2.0

    if np.all(np.logical_not(mask)):
        return image
    masked_image = image[mask]
    image_min = np.min(masked_image)
    image_max = np.max(masked_image)
    image_delta = image_max - image_min
    if image_delta == 0:
        return image

    #
    # ds = downsampled. Calculate the ds array sizes and sigmas.
    #
    ds_sigma_spatial = sigma_spatial / sampling_spatial
    ds_sigma_range   = sigma_range / sampling_range
    ds_i_limit       = int(image.shape[0] / sampling_spatial) + 2
    ds_j_limit       = int(image.shape[1] / sampling_spatial) + 2
    ds_z_limit       = int(image_delta / sampling_range) + 2

    grid_data    = np.zeros((ds_i_limit, ds_j_limit, ds_z_limit))
    grid_weights = np.zeros((ds_i_limit, ds_j_limit, ds_z_limit))
    #
    # Compute the downsampled i, j and z coordinates at each point
    #
    di,dj = np.mgrid[0:image.shape[0],
                     0:image.shape[1]].astype(float) / sampling_spatial
    dz = (masked_image - image_min) / sampling_range
    #
    # Treat this as a list of 3-d coordinates from now on
    #
    di = di[mask]
    dj = dj[mask]
    #
    # scatter the unmasked image points into the data array and
    # scatter a value of 1 per point into the weights
    #
    grid_data[(di + .5).astype(int),
              (dj + .5).astype(int),
              (dz + .5).astype(int)] += masked_image
    grid_weights[(di + .5).astype(int),
                 (dj + .5).astype(int),
                 (dz + .5).astype(int)] += 1
    #
    # Make a Gaussian kernel
    #
    kernel_spatial_limit = int(2 * ds_sigma_spatial) + 1
    kernel_range_limit   = int(2 * ds_sigma_range) + 1
    ki,kj,kz = np.mgrid[-kernel_spatial_limit : kernel_spatial_limit+1,
                        -kernel_spatial_limit : kernel_spatial_limit+1,
                        -kernel_range_limit : kernel_range_limit+1]
    kernel = np.exp(-.5 * ((ki**2 + kj**2) / ds_sigma_spatial ** 2 +
                           kz**2 / ds_sigma_range ** 2))

    blurred_grid_data = convolve(grid_data, kernel, mode='constant')
    blurred_weights = convolve(grid_weights, kernel, mode='constant')
    weight_mask = blurred_weights > 0
    normalized_blurred_grid = np.zeros(grid_data.shape)
    normalized_blurred_grid[weight_mask] = ( blurred_grid_data[weight_mask] /
                                             blurred_weights[weight_mask])
    #
    # Now use di, dj and dz to find the coordinate of the point within
    # the blurred grid to use. We actually interpolate between points
    # here (both in the i,j direction to get intermediate z values and in
    # the z direction to get the slot, roughly where we put our original value)
    #
    dijz = np.vstack((di, dj, dz))
    image_copy = image.copy()
    image_copy[mask] = map_coordinates(normalized_blurred_grid, dijz,
                                       order = 1)
    return image_copy