Estimation of covariance matrices is often an integral part in many signal processing algorithms. In some applications, the covariance matrices can be assumed to have certain structure. Imposing this structure in the estimation typically leads to improved accuracy and robustness (e.g., to small sample effects). In MIMO communications or in signal modelling of EEG data the full covariance matrix can sometimes be modelled as the Kronecker product of two smaller covariance matrices. These smaller matrices may also be structured, e.g., being Toeplitz or at least persymmetric. In this paper we discuss a recently proposed closed form maximum likelihood (ML) based method for the estimation of the Kronecker factor matrices. We also extend the previously presented method to be able to impose the persymmetric constraint into the estimator. Numerical examples show that the mean square errors of the new estimator attains the Cramer-Rao bound even for very small sample sizes.