Playing around with my program I found a remarkable result. I used the function f(z)=z*z*z*z+z*z+0.000002 in the iteration. Because of the leading 4th power this multiplies the complex phase by four for large z. Thus the compensation is minus three times the phase at each iteration step. To get a more uniform background color I subtracted from the final phase the amount of the phase of the initial z. This gives:
Note that the Julia set of this function has a similar outline.
Note added later (5th March 2013):
Having a cold and being very tired I made a programming mistake, giving this result. Instead of adding I subtracted the imaginary part of the z-squared. This amounts to using its complex conjugate. Actually, we can use more general mappings of the (x,y)-plane onto itself to make self-similar images. The mapping simply has to have a contracting region and does not need to be a complex function. The code of the iterated function I used here is