Page MenuHomec4science

Computing rms surface slope
Closed, ResolvedPublic



I am using tamaas to generate a self-affine fractal surface with the required rms surface slope. The tutorial on the documentation page was quite helpful in this regard. However, I have two questions related to rms surface slope calculation. (1) the value given by "spectrum.rmsSlopes()" function has an extra factor of (2*pi) when compared with the analytical relation given for instance in Persson, B. N. J. (2014). Perhaps this is probably due to the difference in definitions of q0 and q1. (2) when the surface slope is calculated in python using the gradient function, the calculated value is smaller by a factor of sqrt(2). To be clear, consider for instance the following:

q0 = 1
q1 = 1
q2 = 128
H = 0.5

spectrum = tamaas.Isopowerlaw2D()
spectrum.q0 = q0
spectrum.q1 = q1
spectrum.q2 = q2
spectrum.hurst = H

rms_heights = spectrum.rmsHeights()
rms_slopes = spectrum.rmsSlopes() #This value has an extra factor of 2*pi
#For this case, rms_slopes=177

N = 256  #Number of points (should be at least 2 * q2)

generator = tamaas.SurfaceGeneratorRandomPhase2D([N, N])
generator.spectrum = spectrum
generator.random_seed = 42

surface = generator.buildSurface()

grad = numpy.gradient(surface,1/(N-1),1/(N-1))
rms_slope_cal = numpy.sqrt(numpy.mean(grad[0]**2+grad[1]**2))  #this value is off by a factor of sqrt(2) from rms_slopes. 
#For this case, rms_slope_cal=129.1 and rms_slopes/rms_slope_cal=1.3699991200892385

Therefore, I request you to please clarify these two things.

Thanks and regards,
Yaswanth Sai

Event Timeline

yaswanth1947 created this object in space S1 c4science.
yaswanth1947 created this object with visibility "Public (No Login Required)".
frerot triaged this task as Normal priority.EditedJun 2 2021, 22:34

Hello Yaswanth,

Thanks for reporting the discrepancies. For your first point, I don't fully understand Persson's definition of kappa^2: the integral term is the second moment of the PSD, but since q in the integral is a frequency coordinate it seems to me that there's a missing 2pi factor (due to taking squaring the gradient components in the Fourier domain), as well as a factor 2 for the sum of the x and y gradient components. In Tamaas, the rms of slopes is 2pi * sqrt(2 * m02), where m02 is the second moment of the PSD defined with frequency coordinates (cf. Yastrebov et al, IJSS, 2015).

For your second point, the finite difference estimator you use is biased due to discretization error, see figure below (from chapter 1 of my thesis).

To reduce the bias you can either use tm.Statistics2D.computeSpectralRMSSlope or increase the number of points compared to q2.

Hope that answers your questions,


Thank you very much, Lucas. Increasing the number of points compared to q2 reduced the bias like you suggested.

For the first point, I looked at Yastrebov et al, (IJSS, 2015) and the rms slope is given to be sqrt(m20+m02) in page 88, which becomes sqrt(2*m20)=sqrt(2*m02) for an isotropic random field. However, you have mentioned it to be 2*pi*sqrt(2*m02). I am not sure which is the correct definition.

Also, they have given the expression for m20 in Appendix A, which when calculated making appropriate substitutions like A = C0/(q0^(-2*(1+H))), T(2) = pi, xi = 1, and zeta = (q1/q0), to compare with the expression given in Persson, B. N. J. (2014), we get 2*m20 = kappa^2.

frerot added a comment.Jun 3 2021, 21:30

It's all a matter of how the Fourier transform is normalized (the 1/2pi factors in forward/backward transforms) and which system of coordinates (frequencies or angular frequencies). Some people define the psd with a normalizing factor, some don't. However, the answer computed by Tamaas (both analytically in the spectrum object and numericaly with computeSpectralRMSSlope) agree with a finite difference calculation, which means that the normalization I use is at least self-consistent.

Thanks for the clarification. I just wanted to confirm as I was already confused with some other things in finding the rms slope.

frerot closed this task as Resolved.Jun 3 2021, 21:46