Context: A variety of formation scenarios have been proposed to explain the diversity of properties observed in bulges. Studying their intrinsic shape can help to constrain the dominant mechanisms at the epochs of their assembly. Aims: The structural parameters of a magnitude-limited sample of 148 unbarred S0-Sb galaxies were derived in order to study the correlations between bulges and disks, as well as the probability distribution function of the intrinsic equatorial ellipticity of bulges. Methods: We present a new fitting algorithm (GASP2D) to perform two-dimensional photometric decomposition of the galaxy surface-brightness distribution. This was assumed to be the sum of the contribution of a bulge and disk component characterized by elliptical and concentric isophotes with constant (but possibly different) ellipticity and position angles. Bulge and disk parameters of the sample galaxies were derived from the J-band images, which were available in the Two Micron All Sky Survey. The probability distribution function of the equatorial ellipticity of the bulges was derived from the distribution of the observed ellipticities of bulges and misalignments between bulges and disks. Results: Strong correlations between the bulge and disk parameters were found. About 80% of bulges in unbarred lenticular and early-to-intermediate spiral galaxies are not oblate but triaxial ellipsoids. Their mean axial ratio in the equatorial plane is < B/A > = 0.85. Their probability distribution function is not significantly dependent on morphology, light concentration or luminosity. The possible presence of nuclear bars does not influence our results. Conclusions: The interplay between bulge and disk parameters favors scenarios in which bulges have assembled from mergers and/or have grown over long times through disk secular evolution. However, all these mechanisms have to be tested against the derived distribution of bulge intrinsic ellipticities. Table 3 is only available in electronic form at http://www.aanda.org