File Exchange

image thumbnail

Radially averaged surface roughness/topograph​y power spectrum (PSD)

version 1.2.0.0 (3.49 KB) by Mona Mahboob Kanafi
Calculates radially averaged 2D power spectrum for a surface roughness/topography

23 Downloads

Updated 19 Dec 2016

View Version History

View License

In surface roughness analysis, one of the powerful tools for roughness characterization is surface roughness power spectrum. If the surface under study has isotropic roughness characteristics, then one can do a radial average on the calculated discrete Fourier transform of the surface topography and obtain its 2D power spectrum, namely, 2D PSD.
With this function you can calculate 2D PSD of a surface topography, which the topography is normally obtained by any 3D profilometry techniques, such as AFM (Atomic Force Microscopy), WLI (White Light Interferometry) and many other optical profilers.
As inputs, you first need to have a matrix (n by m) of your height values (z). Second, you need to know you PixelWidth (spatial resolution) which is obtained easily by dividing you image length to the number of pixels in length, i.e. :

PixelWidth = Lx / m; % in SI units

In order to plot the output:
loglog(q,C)

In order to check 2D FFT of your topography:
imagesc(1+log10(abs(PSD.Hm)))

An extended version of this code, calculates the surface roughness power spectrum only for a portion of the top or bottom part of the surface roughness/topography. This has been uploaded in here:

Radially averaged surface roughness power spectrum (PSD) only on top or bottom part of a topography
http://se.mathworks.com/matlabcentral/fileexchange/60771-radially-averaged-surface-roughness-power-spectrum--psd--only-on-top-or-bottom-part-of-a-topography

To calculate 1-dimensional suface roughness power spectrum, refer to:

1-Dimensional surface roughness power spectrum of a profile or topography
http://se.mathworks.com/matlabcentral/fileexchange/54315-psd-1d-z--pixelwidth--dim-

To generate an artificial randomly rough surface/topography to try these codes, refer to:

Surface generator: artificial randomly rough surfaces
http://se.mathworks.com/matlabcentral/fileexchange/60817-surface-generator--artificial-randomly-rough-surfaces

Cite As

Mona Mahboob Kanafi (2021). Radially averaged surface roughness/topography power spectrum (PSD) (https://www.mathworks.com/matlabcentral/fileexchange/54297-radially-averaged-surface-roughness-topography-power-spectrum-psd), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (41)

Elijah Uche

Hi Mona,
Please I would like to be able to obtain the underlying equations and maths on which this rough surface model was built on so I could fully appreciate and understand how it works
Thanks

Gautier Duroux

Hi,

Some tips to convert a RGB topography image (.png or .JPG) into .csv file ?
I know the color scale in height but not how to get a .csv file with each cell corresponding to a pixel with xyz coordinates from that.

Thanks in advance.

Gautier Duroux

Hi everyone,

I'm new to MatLab and I must admit it's a bit difficult for me to run this code (which seems really nice by the way).
I need to get the 2D PSD of a topographic image (RGB).

My first problem is that I don't see any line where I can enter my image file name in order to measure its PSD. It seems to only have the height topography z, the matric n*m and the piwel width as inputs but no image file name. Did I missed it or do I have to write something for this ?

My second problem is that I am not sure where to enter the inputs. It's a basic thing but very problematic for me. Any idea of the concerned lines ?

Thanks for your help !

Elijah Uche

Hi Mona,
Could you assist me with integrating some dielectric constants like conductivity, permittivity and permeability into a rough surface model
Second, any ideas on how to do a sanity check on a rough surface model to ascertain it does what it says it does...
Thanks
Elijah

Guangzhi XU

Hello Mona, could you help me with the line that deals with the 2D PSD scaling:
```
Cq = (1/U)*(a^2/((n*m)*((2*pi)^2)).*((abs((Hm))).^2));
```
The 1/U part is for the window function, (a.*abs(Hm))^2 is the squared magnitude of discrete Fourier coefficient, right? I'm puzzled at the n, m and 2*pi terms. Can you explain a bit about this normalization scheme? Many thanks.

Qian Wang

Marzieh Salehi

Michael Eller

Hi Mona,

Could you please describe how you chose your window function and why? If you want the best estimate of surface roughness, wouldn't you want to use the raw data? My thinking is that the window functions just mask the smaller components in reciprocal space. Does the Welch window have high dynamic range?

I know very little about image processing and PSD, so I'm just curious to hear why you did it this way.

Thanks!

tqf777

Hi Mona

I just deleted the flooring that you do and the problem of the scaling dissapears as expected.

Cheers!

Hi Mona
Same problem as Nathan here, the code seems to be very susceptible to sizing of the pixels. Any way we could ensure the same results can be achieved for a specific image regardless of the scaling?

Thanks!

Olivia Licata

Hello Mona, In one of your comments you said that the 'RMS roughness is the area under the power spectral density graph. To calculate it, you need to take the 2D matrix of calculated C(q) [not just radially averaged C(q)].' Is the code above only radially averaging C(q)? Should I add the additional lines that you mentioned in order to calculate RMS roughness? (in my case h(RMS) h' and h'') Thank you!

Ari Tuononen

sandra arias

Hi Mona,
Thanks for the code. It works perfectly. I would appreciate if you can help me with the following questions:
(1) My topography has some regular patterns (ripples). Can I still use this code and get a reliable PSD?
(2) I am interested in the RMS slope and RMS curvature. How can I get those if the RMS height is a value?

Thanking you in advance,
Sandra

omar Abboosh

Hi,
I am very beginner in using matlab , There is a plot that achieved it by matlab, I don't how to analysis this plot.
Please if any could help here.

Carpenter86

Hi, a very nice code! Thanks!

I am wondering, if you could also share it on a tribology web project (http://www.tribonet.org/cmdownloads/)? We are exchanging some codes there.

surajit das

Hi can you please give me the coding of surface roughness for journal bearing?

Nathan

I have some strange issues. When I report my pixel with in km, 0.1, the function runs fine. If i report the pixel width in meters, 100, all the outputs are 0. Seems odd, any suggestions?

Anyways this looks to be a pretty useful script. Thanks!

Karishma Karkera

RGB = imread('Mariela.png'); % reads the image
II = im2double(RGB); % convert to double precision
I2 = rgb2gray(II); % convert to 2d gray scale between 0 & 1
% scale gray scale image (between 0 & 1) to your height map
% you need to know your maximum and minimum z values.
z = I2 - min(I2(:));
z = ((z/range(z(:)))*(maxval-minval)) + minval;

This is regarding this code, the maxval corresponds to max value of z right?

Minkook Park

Wow, Thank you very much. I saved my time alot with your codes. If possible, I would like to ask you how to remove the tilt of the topography. tilt of 1D roughness data can be removed using polyfit polyval. But 2D... I think it's not simple. can you please add or post code about it??

Philipp Elmlinger

Hi Mona,
Thank you so much. I inserted your RMS_F2D code at the end of my psd-script in order to calculate the RMS only of spatial frequencies larger q_limit=50 µm or smaller -q_limit. For everyone who is interested, see attached code:

q_limit=50E-6; %define q_limit

%get index of negative limit on qx axis
for i=1:m./2+1;
if -q_limit>qx(i);
qx_low=i;
end;
end;
%get index of positive limit on qx axis
for i=m./2+1:m;
if q_limit>qx(i);
qx_high=i+1;
end;
end;
%get index of negative limit on qy axis
for i=1:n./2+1;
if -q_limit>qy(i);
qy_low=i;
end;
end;
%get index of positive limit on qy axis
for i=n./2+1:n;
if q_grenz>qy(i);
qy_high=i+1;
end;
end;

%Calculate undesired area S_filter under 2D-PSD between points (qx_low/qy_low), (qx_low/qx_high), (qx_high/qy_low) and (qx_high/qy_high)
S_filter=0; %define variable S_filter
for i=qx_low:qx_high;
for j=qy_low:qy_high;
S_filter=S_filter+Cq(j,i);
end;
end;

%now I use your code and subtract the undesired area S_filter from the area of the whole spectrum. Since I filter out low frequencies including the mean, I do not need to remove the mean manually. Also, I reduce the amount of pixels by the pixels of the undesired area.

RMS_F2D = sqrt(((sum(sum(Cq)))-S_filter)*(((2*pi)^2)/(((n*m)-((qx_high-qx_low)*(qy_high-qy_low)))*(PixelWidth.^2))))

If anyone wants to filter out high frequency components, change the last line to:
RMS_F2D = sqrt(S_filter*(((2*pi)^2)/((qx_high-qx_low)*(qy_high-qy_low)*(PixelWidth.^2)))

If you want to calculate the RMS between an upper and lower limit, you have to calculate one further area and subtract accordingly. A sketch of Cq vs. qx and qy including your frequency limits will help a lot.

Philipp

Mona Mahboob Kanafi

Hi Philip,

As you already know, RMS roughness is the area under the power spectral density graph. To calculate it, you need to take the 2D matrix of calculated C(q) [not just radially averaged C(q)]. For instance, below simple code gives you the RMS calculated over the whole range of wavelengths:

[n , m] = size(z);
PSD.Cq(n/2+1,m/2+1) = 0; % remove mean
RMS_F2D = sqrt((sum(sum(PSD.Cq)))*(((2*pi)^2)/(n*m*(PixelWidth.^2))));

The rests becomes a bit tricky. The RMS roughness of original surface would be:
RMS_orig = rms(PSD.z(:) - mean(PSD.z(:)))
These two must be very close to each other, but not exactly the same. The reason while they are not the same, is the fact that we have to use window function during spectral calculation to overcome non-periodic topographies. After window function is applied, the original rms of the topography drops, which then is compensated by a factor called U (normalization constant) in the code. If you are really curious to get exactly the same RMS roughness calculations from spatial and time domain, then set U in the code to be 1 (no compensation), and then compare RMS roughness of the windowed z profile (PSD.z_win) with the RMS roughness in frequency domain as explained above. These two then must be the same.

About your second question, you can calculate the area under the 2D PSD from any range you want. So, continue the same procedure as I explained above for your region of interest.

-Mona

Philipp Elmlinger

Hi, can you please add or post code how to calculate the RMS roughness from the Power Spectral Density C(q) over the whole spectrum of q or any customized spectral region q1 to q2?
Thanks

nipawan pakoktom

Mona Mahboob Kanafi

Hi Philipp,

There are some codes in MATLAB file exchange which you can use to interpolate the missing NaN elements beforehand. For example:

"Inpaint over missing data in 1-D, 2-D, 3-D,... N-D arrays" by Damien Garcia

is a nice contribution. You only need to apply this line before PSD calculations:

z = inpaintn(z);

Hope it helps.

Philipp Elmlinger

Hi,
thank you very much for this tool which works great for me if the inserted surface data z is complete. However, for some very rough surfaces I have NaNs in my z-data (topography) due to measurement issues and hence the whole fft2 results consist of NaNs. Due you have an idea how to handle such incomplete data?
Thanks Philipp

Mona Mahboob Kanafi

You can remove NaN values with any method, there is no difference. About the weird result, if you provide the code the correct inputs, it should give you the correct PSD, I'm not aware of the nature of your data, so the word weird is ambiguous for me.

Sunghoon Lee

I tried to remove nan values (isfinite)and calculate mean values (mean), but the results seems weird. Please let me know your suggestions.

Thanks

Sunghoon Lee

Hi Kanafi,
Thanks a lot! Problem is solved, but another error is came out. 'nanmean' which is in line 106. (I don't have statistics toolbox neither.) From the manual, I found out that it ignores NaN values, thus when I remove that line, it works now. Do you think it will be okay? Or if you have any suggestions, please let me know.

Many thanks,

Mona Mahboob Kanafi

Hi Lee,

Just replace "U = (winInt/recInt); % Normalization constant" with below line:

U = sum(sum(win.^2))/(n*m);

Then, you won't need rectwin at all. Let me know if it didn't fix your problem.

Sunghoon Lee

Because I don't have Signal process toolbox, I cannot use "rectwin" function. Does anyone can let me know how can I modify it? Or, what are the outputs from rectwin(n) or rectwin(m)?

Mona Mahboob Kanafi

Hi Matheu,

You probably don't provide any inputs to the code. You need to inputs: a matrix of height values (z), and your pixelwidth (lateral spacing) between your data points.

Like this:
[q , C, PSD] = psd_2D(z , PixelWidth)

If you still have problem, please explain more of the nature of your data, so that I can help you.

matheu Broom

Hi,

I cant seem to get this to work, I keep getting this error.

Not enough input arguments.
Error in psd_2D (line 36)
[n,m] = size(z);

Any help would be most appreciated

Matheu

Amos

OK, thanks for the explanation.

Mona Mahboob Kanafi

Hello Amos,

That line is just to prevent floating point error calculations in MATLAB, i.e. MATLAB calculates (1 - 2/3 - 1/3) ~=0

'rho' is later used in line 105, which there, for this part 'rho >= q(j)'
maybe MATLAB can not understand that there is a point where rho = q(j).

But since radially averaged PSD is a rough estimate in real roughness topographies, this minor effect will not have impact on the final results, so I wouldn't worry about that; and just ignore that line.

Amos

Hi, thanks for the great submission. It is really usefull for us.
One question: can anybody tell me for what line 95 is?
rho = floor(rho);
In our case, this leads to very few points in the plot of the radially averaged PSD - if we comment line 95, everything looks fine. Can we safely skip that line?

Mona Mahboob Kanafi

Hi Mariela,

Here I assume that you have a roughness color image. As an input for the code you need a double precision point cloud (height map). All you need to do is to convert your RGB image to double precision matrix. Here I will give you an example how you can convert your RGB image:

RGB = imread('Mariela.png'); % reads the image
II = im2double(RGB); % convert to double precision
I2 = rgb2gray(II); % convert to 2d gray scale between 0 & 1

% scale gray scale image (between 0 & 1) to your height map
% you need to know your maximum and minimum z values.
z = I2 - min(I2(:));
z = ((z/range(z(:)))*(maxval-minval)) + minval;

Then with this z, you can continue as the code instructed.

If this is not what you're looking for, please let me know, with more explanation of your application area.

MarielaH

Hi, I have an image (1200x1600x3 uint 8), Can you please give more details how to apply your code for images?

MarielaH

Hi, I have a image (1200x1600x3 uint 8), Can you please give more details how to apply your code for image?

Eric_sch

Thanks for the code, works well and saved me a lot of time.

Christian Wilbs

As far as I can evaluate...AWESOME!!!
Exactly what I needed. I will test the software within my Masterthesis in detail and report later on.
Thank you very much!!!

MATLAB Release Compatibility
Created with R2016a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!