  We describe a fast, optimal estimator for measuring angular bispectra between
two correlated weakly non-Gaussian fields ($Y$ and $Z$) from observational
datasets, based on a separable modal bispectrum expansion. Our methodology is
applicable to (1) any shape of the input theoretical bispectrum templates
(factorizable or not), (2) both even and odd $\ell_1 + \ell_2 + \ell_3$
multipole domains and (3) both amplitude ($f_{\rm NL})$ bispectrum estimation
and full bispectrum reconstruction, considering either joint estimation of
($YYY$, $ZZZ$, $YYZ$ and $ZZY$) shapes, or independent estimation of
auto-bispectra ($YYY$ or $ZZZ$) and cross-bispectra ($YYZ$ or $ZZY$); hence, it
has quite high versatility. The methodology described here was implemented and
used for the official analysis of temperature and polarization cosmic microwave
background maps from the $Planck$ satellite.
