Abstract We discuss the statistical distribution of galaxy shapes and viewing angles under the assumption of triaxiality by deprojecting observed surface brightness profiles of 56 brightest cluster galaxies (BCGs) coming from a recently published large deep-photometry sample. For the first time, we address this issue by directly measuring axis ratio profiles without limiting ourselves to a statistical analysis of average ellipticities. We show that these objects are strongly triaxial, with triaxiality parameters 0.39 ≤ T ≤ 0.72, they have average axis ratios 〈p(r)〉 = 0.84 and 〈q(r)〉 = 0.68, and they are more spherical in the central regions but flatten out at large radii. Measured shapes in the outskirts agree well with the shapes found for simulated massive galaxies and their dark matter halos from both the IllustrisTNG and the Magneticum simulations, possibly probing the nature of dark matter. In contrast, both simulations fail to reproduce the observed inner regions of BCGs, producing objects that are too flattened.