def demean_image(in_file, in_mask=None, out_file=None):
    """
    Demean image data inside mask
    """
    import numpy as np
    import nibabel as nb
    import os.path as op
    from nipype.utils import NUMPY_MMAP
    if out_file is None:
        fname, fext = op.splitext(op.basename(in_file))
        if fext == '.gz':
            fname, _ = op.splitext(fname)
        out_file = op.abspath('./%s_demean.nii.gz' % fname)
    im = nb.load(in_file, mmap=NUMPY_MMAP)
    data = im.get_data().astype(np.float32)
    msk = np.ones_like(data)
    if in_mask is not None:
        msk = nb.load(in_mask, mmap=NUMPY_MMAP).get_data().astype(np.float32)
        msk[msk > 0] = 1.0
        msk[msk < 1] = 0.0
    mean = np.median(data[msk == 1].reshape(-1))
    data[msk == 1] = data[msk == 1] - mean
    nb.Nifti1Image(data, im.affine, im.header).to_filename(out_file)
    return out_file
