def twosample_test(data_images, vardata_images, mask_images, labels, stat_id,
                   permutations=0, cluster_forming_th=0.01):
    """
    Helper function for permutation-based mass univariate twosample group
    analysis. Labels is a binary vector (1-2). Regions more active for group
    1 than group 2 are inferred.
    """
    # Prepare arrays
    data, vardata, xyz, mask = prepare_arrays(data_images, vardata_images,
                                              mask_images)

    # Create two-sample permutation test instance
    if vardata_images is None:
        ptest = permutation_test_twosample(
            data[labels == 1], data[labels == 2], xyz, stat_id=stat_id)
    else:
        ptest = permutation_test_twosample(
            data[labels == 1], data[labels == 2], xyz,
            vardata1=vardata[labels == 1], vardata2=vardata[labels == 2],
            stat_id=stat_id)

    # Compute z-map image
    zmap = np.zeros(data_images[0].shape).squeeze()
    zmap[list(xyz)] = ptest.zscore()
    zimg = Image(zmap, get_affine(data_images[0]))

    # Compute mask image
    maskimg = Image(mask, get_affine(data_images[0]))

    # Multiple comparisons
    if permutations <= 0:
        return zimg, maskimg
    else:
        # Cluster definition: (threshold, diameter)
        cluster_def = (ptest.height_threshold(cluster_forming_th), None)

        # Calibration
        voxel_res, cluster_res, region_res = \
            ptest.calibrate(nperms=permutations, clusters=[cluster_def])
        nulls = {}
        nulls['zmax'] = ptest.zscore(voxel_res['perm_maxT_values'])
        nulls['s'] = cluster_res[0]['perm_size_values']
        nulls['smax'] = cluster_res[0]['perm_maxsize_values']

        # Return z-map image, mask image and dictionary of null
        # distribution for cluster sizes (s), max cluster size (smax)
        # and max z-score (zmax)
        return zimg, maskimg, nulls