Scale invariance and $1/f^2$ power spectrum In the paper 
Occlusion Models for Natural Images : A Statistical Study of a Scale-Invariant Dead Leaves Model; Lee, A. B. Mumford, D. B. Huang, J.; International Journal of Computer Vision
I read that full scale invariance implies a power spectrum of $1/f^2$, i.e. that for a two dimensional stationary random process $Pr\{I(x,y)\}=Pr\{I(sx,sy)\}$ implies a $1/f^2$ form of the Fourier transform of the auto-covariance matrix. 
I read this in several articles by now, but I am not able to find I proof that I could follow or to prove it myself. I am also confused about the direction of implication. One proof I found showed, that a $1/f^2$ power spectrum implies a scale invariant auto-covariance matrix. However, the authors from the paper above seem to imply the opposite. Additionally, a flat power spectrum (i.e. white noise) would also fulfill the scale invariance property from above but would clearly be not $1/f^2$.
In summary, I would be grateful if someone could  


*

*hint to literature, where this is shown in detail

*provide a proof

*clarify the direction of implication or the conditions under which this is usually shown.


Thanks a lot in advance. 
 A: This is a special case of a general principle of scale invariance in random fields.

Definition A random field $(\varphi(x))_{x\in\mathbb{R}^n}$ is called scale invariant (or self-similar) of order $H$ if the law of $(\varphi(\lambda x))_{x\in\mathbb{R}^n}$ is the same as the law of $(\lambda^H\varphi(x))_{x\in\mathbb{R}^n}$.

A stationary random field is invariant under translations, so the autocovariance function $k(x,y)$ is a function of the difference only, i.e. $k(x-y)$. 
We can then relate the behavior of the power spectral density to the behavior of the finite-dimensional probability densities, which transform as
$$p_m(\varphi_1,\lambda x_1;...;\varphi_m,\lambda x_m)=\lambda^{mH}p_m(\lambda^H\varphi_1,x_1;...;\lambda^H\varphi_m,x_m). \tag1$$
First, we use substitution $y=\lambda x$ to rewrite the power spectral density of the scaled variable as
$$S(\lambda k)=\int_{\mathbb{R}^n}e^{i\lambda k\cdot x}k(x)dx=\lambda^{-n}\int_{\mathbb{R}^n}e^{ik\cdot y}k(\lambda^{-1}y)dy. \tag2$$
Now we need to figure out what happens to the autocovariance function under scaling.
Without loss of generality we may assume that the random field has zero mean.
Then, using (1) for $m=2$, we have
$$k(\lambda x)=\int_{\mathbb{R}^2}\varphi_1\varphi_2 p_2(\varphi_1,0;\varphi_2,\lambda x)d\varphi_1 d\varphi_2 \\
=\lambda^{2H}\int_{\mathbb{R}^2}\varphi_1\varphi_2 p_2(\lambda^H\varphi_1,0;\lambda^H\varphi_2,x)d\varphi_1 d\varphi_2 \\
=\lambda^{-2H}\int_{\mathbb{R}^2}\tilde\varphi_1\tilde\varphi_2 p_2(\tilde\varphi_1,0;\tilde\varphi_2,x)d\tilde\varphi_1 d\tilde\varphi_2 \\
=\lambda^{-2H} k(x),$$
where in the second-to-last line we used the substitutions $\tilde\varphi_i=\lambda^H\varphi_i$ for $i=1,2$.
This, together with (2), implies
$$S(\lambda k)=\lambda^{2H-n}S(k).$$
So for a self-similar random field in 2d ($n=2$) we have that $S(\lambda k)=\lambda^{-2}S(k)$ iff $H=0$.
This value of $H$ is referred to as full scale invariance in the paper you cited.
