def _create_cifti_image(bold_file, label_file, annotation_files, gii_files,
                            volume_target, surface_target, tr):
        """
        Generate CIFTI image in target space

        Parameters
            bold_file : 4D BOLD timeseries
            label_file : label atlas
            annotation_files : FreeSurfer annotations
            gii_files : 4D BOLD surface timeseries in GIFTI format
            volume_target : label atlas space
            surface_target : gii_files space
            tr : repetition timeseries

        Returns
            out_file : BOLD data as CIFTI dtseries
        """

        label_img = nb.load(label_file)
        bold_img = resample_to_img(bold_file, label_img)

        bold_data = bold_img.get_data()
        timepoints = bold_img.shape[3]
        label_data = label_img.get_data()

        # set up CIFTI information
        series_map = ci.Cifti2MatrixIndicesMap((0, ),
                                               'CIFTI_INDEX_TYPE_SERIES',
                                               number_of_series_points=timepoints,
                                               series_exponent=0,
                                               series_start=0.0,
                                               series_step=tr,
                                               series_unit='SECOND')
        # Create CIFTI brain models
        idx_offset = 0
        brainmodels = []
        bm_ts = np.empty((timepoints, 0))

        for structure, labels in CIFTI_STRUCT_WITH_LABELS.items():
            if labels is None:  # surface model
                model_type = "CIFTI_MODEL_TYPE_SURFACE"
                # use the corresponding annotation
                hemi = structure.split('_')[-1]
                annot = nb.freesurfer.read_annot(annotation_files[hemi == "RIGHT"])
                # currently only supports L/R cortex
                gii = nb.load(gii_files[hemi == "RIGHT"])
                # calculate total number of vertices
                surf_verts = len(annot[0])
                # remove medial wall for CIFTI format
                vert_idx = np.nonzero(annot[0] != annot[2].index(b'unknown'))[0]
                # extract values across volumes
                ts = np.array([tsarr.data[vert_idx] for tsarr in gii.darrays])

                vert_idx = ci.Cifti2VertexIndices(vert_idx)
                bm = ci.Cifti2BrainModel(index_offset=idx_offset,
                                         index_count=len(vert_idx),
                                         model_type=model_type,
                                         brain_structure=structure,
                                         vertex_indices=vert_idx,
                                         n_surface_vertices=surf_verts)
                bm_ts = np.column_stack((bm_ts, ts))
                idx_offset += len(vert_idx)
                brainmodels.append(bm)
            else:
                model_type = "CIFTI_MODEL_TYPE_VOXELS"
                vox = []
                ts = None
                for label in labels:
                    ijk = np.nonzero(label_data == label)
                    ts = (bold_data[ijk] if ts is None
                          else np.concatenate((ts, bold_data[ijk])))
                    vox += [[ijk[0][ix], ijk[1][ix], ijk[2][ix]]
                            for ix, row in enumerate(ts)]

                bm_ts = np.column_stack((bm_ts, ts.T))

                vox = ci.Cifti2VoxelIndicesIJK(vox)
                bm = ci.Cifti2BrainModel(index_offset=idx_offset,
                                         index_count=len(vox),
                                         model_type=model_type,
                                         brain_structure=structure,
                                         voxel_indices_ijk=vox)
                idx_offset += len(vox)
                brainmodels.append(bm)

        volume = ci.Cifti2Volume(
            bold_img.shape[:3],
            ci.Cifti2TransformationMatrixVoxelIndicesIJKtoXYZ(-3, bold_img.affine))
        brainmodels.append(volume)

        # create CIFTI geometry based on brainmodels
        geometry_map = ci.Cifti2MatrixIndicesMap((1, ),
                                                 'CIFTI_INDEX_TYPE_BRAIN_MODELS',
                                                 maps=brainmodels)
        # provide some metadata to CIFTI matrix
        meta = {
            "target_surface": surface_target,
            "target_volume": volume_target,
        }
        # generate and save CIFTI image
        matrix = ci.Cifti2Matrix()
        matrix.append(series_map)
        matrix.append(geometry_map)
        matrix.metadata = ci.Cifti2MetaData(meta)
        hdr = ci.Cifti2Header(matrix)
        img = ci.Cifti2Image(bm_ts, hdr)
        img.nifti_header.set_intent('NIFTI_INTENT_CONNECTIVITY_DENSE_SERIES')

        _, out_base, _ = split_filename(bold_file)
        out_file = "{}.dtseries.nii".format(out_base)
        ci.save(img, out_file)
        return os.path.join(os.getcwd(), out_file)
Generate CIFTI image in target space

        Parameters
            bold_file : 4D BOLD timeseries
            label_file : label atlas
            annotation_files : FreeSurfer annotations
            gii_files : 4D BOLD surface timeseries in GIFTI format
            volume_target : label atlas space
            surface_target : gii_files space
            tr : repetition timeseries

        Returns
            out_file : BOLD data as CIFTI dtseries