The mean isophotal axis ratios are computed for a complete sample of 119 brightest cluster galaxies (BCGs) in Abell clusters. The cumulative distribution function F(q) of axis ratios for the BCGs differs from F(q) for a sample of ordinary elliptical galaxies at only the P = 0.14 probability level, as measured by a Kolmogorov-Smirnov test. In contrast, a sample of 30 s brightest cluster ellipticals has axis ratios which are significantly rounder (P(K-S) = 0.01) than those of the BCGs. We model the distribution of intrinsic axis ratios f(beta, gamma) as an isotropic Gaussian in the space 0 less-than-or-equal-to gamma less-than-or-equal-to beta less-than-or-equal-to 1, where beta and gamma are the scaled lengths of the intermediate and shortest axes. The best fitting model for BCGs shapes is prolate in character, peaking at beta0 = 0.80, gamma0 = 0.76, with a small dispersion sigma0 = 0.079, while the best model for ordinary elliptical galaxies is more oblate and less homogeneous, with beta0 = 0.98, gamma0 = 0.69, and sigma0 = 0.123. However, the apparent shapes of BCGs and of ordinary ellipticals are jointly well fitted by a Gaussian centered on beta0 = 0.95, gamma0 = 0.70, with sigma0 = 0.10. The second brightest cluster ellipticals are best fitted by a function with a most probable shape which is rounder (beta0 = 0.93, gamma0 = 0.82), and with a narrower spread in axis ratios (sigma0 = 0.082). The difference in shapes may be the result of different merging histories.