def smooth_mean_per_slice(data, mask_data):

    # Create a new image keeping only the spinal cord and assigning to the other voxels the mean value across the spinal cord

    # Find the voxels that belong to the spinal cord
    X, Y, Z = (mask_data > 0).nonzero()

    # Define useful variables
    N = len(X) # number of voxels in the spinal cord segmentation
    Z_min = min(Z) # min z of the segmentation
    Z_max = max(Z) # max z of the segmentation
    Z_nb = Z_max - Z_min + 1 # number of different z-slice of the segmentation
    x = len(data) # number of voxels in the X direction
    y = len(data[1]) # number of voxels in the Y direction
    z = len(data[1][1]) # number of voxels in the Z direction

    ## Count the number of voxels belonging to the spinal cord for each slice according to the segmentation
    #nb_vox_per_slice = [0 for i in range(0,Z_nb)] # initialization
    #z_index = 0
    #while z_index < Z_nb:
    #    nb_vox_per_slice[z_index]+= 1
    #    z_index+= 1


    # Sort by z-slice the values of the voxels belonging to the spinal cord
    sc_values = [[] for Z_index in range(0, Z_nb)] # initialization
    for vox_index in range(0, N):
        sc_values[Z[vox_index] - Z_min].append(data[X[vox_index]][Y[vox_index]][Z[vox_index]])

    # Compute the mean value for each slice of the spinal cord
    print '\nCompute the mean value for each slice of the spinal cord...'
    sc_mean_per_slice = [0 for Z_index in range(0, Z_nb)] # initialization
    for Z_index in range(0, Z_nb):
        sc_mean_per_slice[Z_index] = sum(sc_values[Z_index]) / len(sc_values[Z_index])

    # Define a new image assigning to all the voxels that don't belong to the mean value across their slice
    print '\nCreate a new image to smooth keeping only the spinal cord and completing with the previously computed mean values...'
    sc_data = [[[0 for k in range(0, z)] for j in range(0, y)] for i in
               range(0, x)] # initialization by the size of the original data
    for k in range(0, z):
        for j in range(0, y):
            for i in range(0, x):

                if k < Z_min:
                    sc_data[i][j][k] = sc_mean_per_slice[0]
                elif Z_min <= k <= Z_max:
                    sc_data[i][j][k] = sc_mean_per_slice[k - Z_min]
                elif k > Z_max:
                    sc_data[i][j][k] = sc_mean_per_slice[len(sc_mean_per_slice) - 1]

    # Assign the right value to the voxels that belong to the spinal cord
    for i in range(0, N):
        sc_data[X[i]][Y[i]][Z[i]] = data[X[i]][Y[i]][Z[i]]

    # Apply the filter to these new image
    smoothed_sc_data = apply_filter(sc_data)

    # Replace the original spinal cord data by the smoothed one in the original data
    for i in range(0, N):
        data[X[i]][Y[i]][Z[i]] = smoothed_sc_data[X[i]][Y[i]][Z[i]]

    # Return the corrected data
    return data
Apply the smoothing to an image of the spinal cord completed with the mean value of
dataset/ETHPy150Open neuropoly/spinalcordtoolbox/dev/sct_apply_local_filter.py/smooth_mean_per_slice