The J-band structural parameters of a magnitude-limited sample of 148 unbarred S0--Sb galaxies were obtained using the GASP2D algorithm and then analyzed to recover the intrinsic shape of their bulges. We developed a new method to derive the intrinsic shape of bulges based only on photometric data and on the geometrical relationships between the apparent and intrinsic shapes of bulges and disks. The method is conceived as completely independent of the studied class of objects, and it can be applied whenever a triaxial ellipsoid embedded in (or embedding) an axisymmetric component is considered. We found that the intrinsic shape is well constrained for a subsample of 115 bulges with favorable viewing angles. A large fraction of them are characterized by an elliptical section (B/A<0.9). This fraction is 33%, 55%, and 43% if using their maximum, mean, or median equatorial ellipticity, respectively. Most of them are flattened along their polar axis (C<(A+B)/2). The distribution of triaxiality is strongly bimodal. This bimodality is driven by bulges with Sérsic index n > 2 or equivalently, by the bulges of galaxies with a bulge-to-total ratio B/T > 0.3. Bulges with n ≤ 2 and with B/T ≤ 0.3 follow a similar distribution, which is different from that of bulges with n > 2 and with B/T > 0.3. In particular, bulges with n≤2 and with B/T ≤ 0.3 show a larger fraction of oblate axisymmetric (or nearly axisymmetric) bulges, a smaller fraction of triaxial bulges, and fewer prolate axisymmetric (or nearly axisymmetric) bulges with respect to bulges with n > 2 and with B/T > 0.3, respectively. According to predictions of the numerical simulations of bulge formation, bulges with n ≤ 2, which show a high fraction of oblate axisymmetric (or nearly axisymmetric) shapes and have B/T ≤ 0.3, could be the result of dissipational minor mergers. Both major dissipational and dissipationless mergers seem to be required to explain the variety of shapes found for bulges with n > 2 and B/T > 0.3.