We propose an efficient method for calculating level anti-crossing spectra (LAC spectra) of interacting paramagnetic defect centers in crystals. By LAC spectra we mean the magnetic field dependence of the photoluminescence intensity of paramagnetic color centers: such field dependences often exhibit sharp features, such as peaks or dips, originating from LACs in the spin system. Our approach takes into account the electronic Zeeman interaction with the external magnetic field, dipole-dipole interaction of paramagnetic centers, hyperfine coupling of paramagnetic defects to magnetic nuclei and zero-field splitting. By using this method, not only can we obtain the positions of lines in LAC spectra, but also reproduce their shapes as well as the relative amplitudes of different lines. As a striking example, we present a calculation of LAC spectra in diamond crystals containing negatively charged NV centers.