We estimate the distribution of intrinsic shapes of APM galaxy clusters from the distribution of their apparent shapes. We measure the projected cluster ellipticities using two alternative methods. The first method is based on moments of the discrete galaxy distribution while the second is based on moments of the smoothed galaxy distribution. We study the performance of both methods using Monte Carlo cluster simulations covering the range of APM cluster distances and including a random distribution of background galaxies. We find that the first method suffers from severe systematic biases, whereas the second is more reliable. After excluding clusters dominated by substructure and quantifying the systematic biases in our estimated shape parameters, we recover a corrected distribution of projected ellipticities. We use the non-parametric kernel method to estimate the smooth apparent ellipticity distribution, and numerically invert a set of integral equations to recover the corresponding distribution of intrinsic ellipticities under the assumption that the clusters are either oblate or prolate spheroids. The prolate spheroidal model fits the APM cluster data best.