<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _calc_snr(volume, mask, dilation=5, reference_tr=None, ):
 """ Calculate the the SNR of a volume Calculates the Signal to Noise Ratio, the mean of brain voxels divided by the standard deviation across non-brain voxels. Specify a TR value to calculate the mean and standard deviation for that TR. To calculate the standard deviation of non-brain voxels we can subtract any baseline structure away first, hence getting at deviations due to the system noise and not something like high baseline values in non-brain parts of the body. Parameters volume : 4d array, float Take a volume time series mask : 3d array, binary A binary mask the same size as the volume dilation : int How many binary dilations do you want to perform on the mask to determine the non-brain voxels. If you increase this the SNR increases and the non-brain voxels (after baseline subtraction) more closely resemble a gaussian reference_tr : int or list Specifies the TR to calculate the SNR for. If multiple are supplied then it will use the average of them. Returns ------- snr : float  The SNR of the volume """


    # If no TR is specified then take all of them
    if reference_tr is None:
        reference_tr = list(range(volume.shape[3]))

    # Dilate the mask in order to ensure that non-brain voxels are far from
    # the brain
    if dilation > 0:
        mask_dilated = ndimage.morphology.binary_dilation(mask,
                                                          iterations=dilation)
    else:
        mask_dilated = mask

    # Make a matrix of brain and non_brain voxels, selecting the timepoint/s
    brain_voxels = volume[mask > 0][:, reference_tr]
    nonbrain_voxels = (volume[:, :, :, reference_tr]).astype('float64')

    # If you have multiple TRs
    if len(brain_voxels.shape) > 1:
        brain_voxels = np.mean(brain_voxels, 1)
        nonbrain_voxels = np.mean(nonbrain_voxels, 3)

    nonbrain_voxels = nonbrain_voxels[mask_dilated == 0]

    # Take the means of each voxel over time
    mean_voxels = np.nanmean(brain_voxels)

    # Find the standard deviation of the voxels
    std_voxels = np.nanstd(nonbrain_voxels)

    # Return the snr
    return mean_voxels / std_voxels