the code is in the regionprops function, you can access it via "edit regionprops", and checking the "ComputeEllipseParams" function.
The principle is to compute the inertia ellipse of a shape. It is obtained by computing normalised centered moments. The matrix of Inertia is obtained from the normalised moments of order 2 : Ixx, Ixy, and Iyy. By extracting eigen values and eigen vectors of the inertia matrix, it is possible to obtain the direction of the first eigen vector, that is converted into angle.