This paper is concerned with the direct and inverse scattering of time-harmonic electromagnetic waves from bi-anisotropic media. For the direct problem, we establish an integro-differential equation formulation, its Fredholm property, and the uniqueness of a weak solution. Using this integro-differential formulation we study a fast spectral Galerkin method for the numerical solution to the direct problem. Numerical examples are presented and convergence of the spectral method is proved via G;rding estimates for (strongly singular) integro-differential equations. We solve the inverse problem of recovering bi-anisotropic scatterers from far-field data using orthogonality sampling methods. These methods aim to construct imaging functionals which are robust to noise, computationally cheap, and require data for only one or a few incident fields. We provide some theoretical analysis as well as numerical simulations for the proposed imaging functionals.