Why can't the surfaces of stars be observed? If I'm correctly informed, only three stars: Sol, Betelgeuse and Altair have surfaces which have been resolved by telescopes. All other stars are only point sources of light, even in the greatest of telescopes. Is it just the huge distance to the stars which reduce them to dots? Or is there some other optical explanation? I mean, the Sun is about as large as the Moon on the sky, but I cannot see any details on its surface because it blinds me. However, the next nearest star is invisible to the naked eye.

Some exoplanets have been directly imaged by removing the light from their star. It seems funny that we cannot detect any features on the photosphere of such a star, such as spots or flares or its shape. Plenty of photons are available.

Actually, there are several stars which have been directly imaged, maybe my sopposition is quite wrong:
https://en.wikipedia.org/wiki/List_of_stars_with_resolved_images
There are many different ways to get spatial information about the surface of a star besides direct imaging.

Direct imaging is difficult because the angular resolution available goes as $\lambda/D$.
For a 8-m telescope and light at 500 nm, one can resolve $6\times10^{-8}$ radians (assuming the blurring of the atmosphere can be overcome by adaptive optics or similar).

The nearest stars are a couple of parsecs away, so the smallest spatial scales that could be resolved are $\sim 2\times 3.1\times10^{16} \times 6\times10^{-8} =  3.7\times10^{9}$ m, or about 500 solar radii. Hence no surface features or even a disk could be resolved.

Of course you could use interferometric techniques to effectively increase the size of $D$ and measurements of the angular radii are now possible for many nearby stars or giant stars at larger distances.

Surface imaging is harder. Indirect techniques are much more common. These include doppler imaging and eclipse mapping. The former uses the fact that there is a relationship between the position of a bright/dark feature on a rotating star and the doppler shift of the light from that feature. By observing a time-series of spectra, the lumps and bumps in spectral lines can be inverted to produce a "doppler map" of the surface. The technique is usually limited to stars that are considerably more rapidly rotating than the Sun. There is a lot of ambiguity in the image reconstruction process - many surfaces could lead to the same observable sigature and clever statistical techniques (and even philosphies) have to be deployed to choose between them. Many stars have published "doppler maps" of their surfaces. Here is a typical example of such a study and below I show an example of a "doppler image" for the star II Peg (a K-type subgiant), from Gu et al. (2003), showing dark spotted regions. A typical resolution for such an image is about 10 degrees on the star.



Eclipse mapping, for which I can't easily locate a good link, uses the fact that a star/accretion disk is orbited by another star or planet that periodically eclipses it. What happens to the light from the system in and out of eclipse can be used to probe the surface of the eclipsed object. There are of course limitations to the spatial resolution that can be obtained, depending on the size of the eclipsing object, how long the eclipse takes and how wide the orbit is. But useful constraints can be made on the structure of accretion disks, sizes of starspots etc, although "maps" are not usually produced. A recent example using the transits of a planet to probe starspot structure is Roettenbacher et al. 2013.

Another possibility is rotational modulation. Features on the surface that rotate around are self-eclipsed by the star and produce modulation of the observed light. This can be used to try and estimate e.g. the size and location of starspots. Again there are many degeneracies and ambiguities, but this has become a growth industry since the delivery of thousands of extremely high quality light curves from the Kepler satellite.