def get_3d_peaks(image, mask=None, threshold=0., nn=18, order_th=0):
    """
    returns all the peaks of image that are with the mask
    and above the provided threshold

    Parameters
    ----------
    image, (3d) test image
    mask=None, (3d) mask image
        By default no masking is performed
    threshold=0., float, threshold value above which peaks are considered
    nn=18, int, number of neighbours of the topological spatial model
    order_th=0, int, threshold on topological order to validate the peaks

    Returns
    -------
    peaks, a list of dictionaries, where each dict has the fields:
    vals, map value at the peak
    order, topological order of the peak
    ijk, array of shape (1,3) grid coordinate of the peak
    pos, array of shape (n_maxima,3) mm coordinates (mapped by affine)
        of the peaks
    """
    # Masking
    if mask is not None:
        bmask = mask.get_data().ravel()
        data = image.get_data().ravel()[bmask > 0]
        xyz = np.array(np.where(bmask > 0)).T
    else:
        shape = image.shape
        data = image.get_data().ravel()
        xyz = np.reshape(np.indices(shape), (3, np.prod(shape))).T
    affine = get_affine(image)

    if not (data > threshold).any():
        return None

    # Extract local maxima and connex components above some threshold
    ff = field_from_graph_and_data(wgraph_from_3d_grid(xyz, k=18), data)
    maxima, order = ff.get_local_maxima(th=threshold)

    # retain only the maxima greater than the specified order
    maxima = maxima[order > order_th]
    order = order[order > order_th]

    n_maxima = len(maxima)
    if n_maxima == 0:
        # should not occur ?
        return None

    # reorder the maxima to have decreasing peak value
    vals = data[maxima]
    idx = np.argsort(- vals)
    maxima = maxima[idx]
    order = order[idx]

    vals = data[maxima]
    ijk = xyz[maxima]
    pos = np.dot(np.hstack((ijk, np.ones((n_maxima, 1)))), affine.T)[:, :3]
    peaks = [{'val': vals[k], 'order': order[k], 'ijk': ijk[k], 'pos': pos[k]}
             for k in range(n_maxima)]

    return peaks