(Abridged) A variety of formation scenarios was proposed to explain the diversity of properties observed in bulges. Studying their intrinsic shape can help in constraining the dominant mechanism at the epochs of their assembly. 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 (PDF) of the intrinsic equatorial ellipticity of bulges. It is presented a new fitting algorithm (GASP2D) to perform the two-dimensional photometric decomposition of 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 PDF of the equatorial ellipticity of the bulges was derived from the distribution of the observed ellipticities of bulges and misalignments between bulges and disks. 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 = 0.85. There is not significant dependence of their PDF on morphology, light concentration, and luminosity. The interplay between bulge and disk parameters favors scenarios in which bulges assembled from mergers and/or grew over long times through disk secular evolution. But all these mechanisms have to be tested against the derived distribution of bulge intrinsic ellipticities.