One of the main problems with using supershapes for my purpose is normal calculation. Since the equation is periodic in polar coordinates, most of the generated shapes cannot be seen as single axis rotational surface silhouettes. Even for ones that could, the shader cannot sample all needed information to calculate the normal. So my previous normal generation approach wouldn't work.
I tried to adapt it to pola shapes, treating each sampling angle as a rotational surface slice. As expected, this created harsh streaky patterns radiating from the center:
To somewhat remedy this, I mixed it with spherical normals, where sphere radius is matching shape radius at each radial slice. The mix factor falls off towards the edges, resulting in spherical normal being more prominent near the center, masking out the most distracting streaks area:
After some experimentation I realized that the superformula can be mapped onto linear axis (instead of angle). In this form it yields less varying shapes but they all can be seen as rotational surface silhouettes. This is better suited for tree stylization I'm after. The mapping also allows for more straightforward normal calculation: