Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87047170
surface_statistics.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Oct 10, 05:52
Size
21 KB
Mime Type
text/x-c++
Expires
Sat, Oct 12, 05:52 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21507035
Attached To
rTAMAAS tamaas
surface_statistics.hh
View Options
/**
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __SURFACE_STATISTICS_HH__
#define __SURFACE_STATISTICS_HH__
/* -------------------------------------------------------------------------- */
#include "surface.hh"
#include <fftw3.h>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class
SurfaceStatistics
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public
:
SurfaceStatistics
(){};
virtual
~
SurfaceStatistics
(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public
:
inline
static
Real
computeMaximum
(
const
Surface
<
Real
>
&
s
);
inline
static
Real
computeMinimum
(
const
Surface
<
Real
>
&
s
);
inline
static
UInt
computeAverage
(
const
Surface
<
UInt
>
&
s
,
UInt
epsilon
=
0
){
SURFACE_FATAL
(
"TO IMPLEMENT"
);};
inline
static
unsigned
long
computeAverage
(
const
Surface
<
unsigned
long
>
&
s
,
unsigned
long
epsilon
=
0
){
SURFACE_FATAL
(
"TO IMPLEMENT"
);};
inline
static
Real
computeAverage
(
const
Surface
<
Real
>
&
s
,
Real
epsilon
=
-
1.
*
std
::
numeric_limits
<
Real
>::
max
()){
return
computeAverage
<
Real
>
(
s
,
epsilon
);
}
template
<
typename
T
>
inline
static
T
computeAverage
(
const
Surface
<
T
>
&
s
,
T
epsilon
=
-
1.
*
std
::
numeric_limits
<
T
>::
max
());
inline
static
Real
computeStdev
(
const
Surface
<
Real
>
&
s
,
Real
epsilon
=-
1.
*
std
::
numeric_limits
<
Real
>::
max
());
inline
static
Real
computeSkewness
(
const
Surface
<
Real
>
&
surface
,
Real
avg
,
Real
std
);
inline
static
Real
computeKurtosis
(
const
Surface
<
Real
>
&
surface
,
Real
avg
,
Real
std
);
inline
static
Real
computeRMSSlope
(
const
Surface
<
Real
>
&
surface
);
inline
static
Real
computeSpectralRMSSlope
(
Surface
<
Real
>
&
surface
);
inline
static
Real
computeSpectralMeanCurvature
(
Surface
<
Real
>
&
surface
);
inline
static
Real
computeSpectralStdev
(
Surface
<
Real
>
&
surface
);
// inline static void computeStatistics(Surface<Real> & surface);
inline
static
std
::
vector
<
Real
>
computeMoments
(
Surface
<
Real
>
&
surface
);
inline
static
Real
computeAnalyticFractalMoment
(
UInt
order
,
UInt
k1
,
UInt
k2
,
Real
hurst
,
Real
rms
,
Real
L
);
//! compute perimeter
inline
static
Real
computePerimeter
(
const
Surface
<
Real
>
&
surface
,
Real
eps
=
1e-10
);
//! get real contact area
inline
static
Real
computeContactArea
(
const
Surface
<
Real
>
&
surface
);
//! get real contact area
inline
static
Real
computeContactAreaRatio
(
const
Surface
<
Real
>
&
surface
);
template
<
bool
consider_zeros
=
false
>
inline
static
std
::
vector
<
int
>
computeDistribution
(
const
Surface
<
Real
>
&
surface
,
unsigned
int
n_bins
,
Real
min_value
=
std
::
numeric_limits
<
Real
>::
max
(),
Real
max_value
=-
1
*
std
::
numeric_limits
<
Real
>::
max
());
inline
static
std
::
vector
<
complex
>
computeSpectralDistribution
(
const
Surface
<
Real
>
&
surface
,
unsigned
int
n_bins
,
UInt
cutoff
=
std
::
numeric_limits
<
UInt
>::
max
());
//! compute average amplitude
inline
static
Real
computeAverageAmplitude
(
const
Surface
<
Real
>
&
surface
);
inline
static
Real
computeAverageAmplitude
(
const
SurfaceComplex
<
Real
>
&
surface
);
//! sum all heights
inline
static
Real
computeSum
(
const
Surface
<
Real
>
&
surface
);
//! compute auto correlation
inline
static
Surface
<
Real
>
computeAutocorrelation
(
Surface
<
Real
>
&
surface
);
//! compute powerspectrum
inline
static
SurfaceComplex
<
Real
>
computePowerSpectrum
(
Surface
<
Real
>
&
surface
);
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
//#include "surface_statistics_inline_impl.cc"
/* -------------------------------------------------------------------------- */
template
<
bool
consider_zeros
>
std
::
vector
<
int
>
SurfaceStatistics
::
computeDistribution
(
const
Surface
<
Real
>
&
surface
,
unsigned
int
n_bins
,
Real
min_value
,
Real
max_value
){
if
(
n_bins
==
0
)
SURFACE_FATAL
(
"badly set number of bins"
);
Real
min_heights
;
Real
max_heights
;
std
::
vector
<
int
>
bins
;
bins
.
resize
(
n_bins
);
bins
.
assign
(
n_bins
,
0
);
min_heights
=
min_value
;
max_heights
=
max_value
;
if
(
max_value
==
-
1.
*
std
::
numeric_limits
<
Real
>::
max
()
)
max_heights
=
computeMaximum
(
surface
);
if
(
min_value
==
std
::
numeric_limits
<
Real
>::
max
()
)
min_heights
=
computeMinimum
(
surface
);
int
counter
=
0
;
const
int
n
=
surface
.
size
();
for
(
int
i
=
0
;
i
<
n
*
n
;
++
i
)
{
Real
value
=
surface
(
i
);
if
(
!
consider_zeros
&&
value
<
zero_threshold
)
continue
;
if
(
value
>
max_heights
||
value
<
min_heights
)
SURFACE_FATAL
(
"out of bound value: "
<<
min_heights
<<
" <? "
<<
value
<<
" <? "
<<
max_heights
<<
" : bounds probably badly chosen"
);
Real
ind
=
(
value
-
min_heights
)
/
(
max_heights
-
min_heights
);
unsigned
int
index
=
ind
*
n_bins
;
if
(
index
==
n_bins
)
--
index
;
++
bins
[
index
];
++
counter
;
}
if
(
consider_zeros
&&
counter
!=
n
*
n
)
SURFACE_FATAL
(
"problem in counting points"
);
// Real offset = (max_heights - min_heights)/n_bins;
//Real x = 0;
// Real f1 = 0, f2 = 0;
// for (unsigned int i = 0 ; i < n_bins ; ++i)
// {
// x = i*offset+min_heights;
//f1 = bins[i]/Real(n*n);
//f2 = bins[i]/Real(counter);
//}
Real
min_non_zero
=
std
::
numeric_limits
<
Real
>::
max
();
for
(
int
i
=
0
;
i
<
n
*
n
;
++
i
)
{
Real
value
=
surface
(
i
);
if
(
value
<
min_non_zero
&&
value
>
1e-14
)
min_non_zero
=
value
;
}
return
bins
;
}
/* -------------------------------------------------------------------------- */
std
::
vector
<
complex
>
SurfaceStatistics
::
computeSpectralDistribution
(
const
Surface
<
Real
>
&
surface
,
unsigned
int
n_bins
,
UInt
cutoff
){
if
(
n_bins
==
0
)
SURFACE_FATAL
(
"badly set number of bins"
);
const
int
n
=
surface
.
size
();
Real
min_value
=
SurfaceStatistics
::
computeMinimum
(
surface
);
Real
max_value
=
SurfaceStatistics
::
computeMaximum
(
surface
);
std
::
vector
<
complex
>
distrib_fft
(
n_bins
);
const
int
N
=
n_bins
;
Real
L
=
max_value
-
min_value
;
for
(
int
k
=
0
;
k
<
N
;
++
k
)
distrib_fft
[
k
]
=
0.
;
for
(
int
k
=
1
;
k
<
N
/
2
;
++
k
)
{
for
(
int
i
=
0
;
i
<
n
*
n
;
++
i
)
{
Real
value
=
surface
(
i
);
Real
angle
=
-
2.
*
M_PI
*
k
*
(
value
-
min_value
)
/
L
;
complex
phase
=
std
::
polar
(
1.
,
angle
);
distrib_fft
[
k
]
+=
phase
;
distrib_fft
[
N
-
k
]
+=
phase
;
}
}
for
(
int
i
=
0
;
i
<
n
*
n
;
++
i
)
{
Real
value
=
surface
(
i
);
Real
angle
=
-
2.
*
M_PI
*
N
/
2
*
(
value
-
min_value
)
/
L
;
complex
phase
=
std
::
polar
(
1.
,
angle
);
distrib_fft
[
N
/
2
]
+=
phase
;
}
for
(
int
k
=
1
;
k
<
N
;
++
k
)
distrib_fft
[
k
]
/=
n
*
n
;
distrib_fft
[
0
]
=
std
::
polar
(
1.
,
0.
);
for
(
int
k
=
1
;
k
<
N
/
2
;
++
k
)
{
if
(
k
>=
(
int
)
cutoff
)
{
distrib_fft
[
k
]
=
0.
;
distrib_fft
[
N
-
k
]
=
0.
;
}
}
//proceed to inverse fourier transform
for
(
int
k
=
0
;
k
<
N
;
++
k
)
{
distrib_fft
[
k
]
=
std
::
conj
(
distrib_fft
[
k
]);
}
fftw_complex
*
in
=
(
fftw_complex
*
)
&
distrib_fft
[
0
];
fftw_plan
p
=
fftw_plan_dft_1d
(
N
,
in
,
in
,
FFTW_FORWARD
,
FFTW_ESTIMATE
);
fftw_execute
(
p
);
fftw_destroy_plan
(
p
);
std
::
vector
<
complex
>
distrib
;
distrib
.
resize
(
N
);
for
(
int
k
=
0
;
k
<
N
;
++
k
)
{
distrib
[
k
]
=
std
::
conj
(
distrib_fft
[
k
]);
}
return
distrib
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeAverageAmplitude
(
const
Surface
<
Real
>
&
surface
){
Real
average_amplitude
=
0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
for
(
UInt
j
=
0
;
j
<
n
;
++
j
)
average_amplitude
+=
std
::
abs
(
surface
(
i
,
j
));
average_amplitude
*=
1.
/
n
/
n
;
return
average_amplitude
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeAverageAmplitude
(
const
SurfaceComplex
<
Real
>
&
surface
){
Real
average_amplitude
=
0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
for
(
UInt
j
=
0
;
j
<
n
;
++
j
)
average_amplitude
+=
std
::
abs
(
surface
(
i
,
j
));
average_amplitude
*=
1.
/
n
/
n
;
return
average_amplitude
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeContactArea
(
const
Surface
<
Real
>
&
surface
)
{
Real
L
=
surface
.
getL
();
return
SurfaceStatistics
::
computeContactAreaRatio
(
surface
)
*
L
*
L
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeContactAreaRatio
(
const
Surface
<
Real
>
&
surface
)
{
UInt
area
=
0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
;
i
++
)
{
for
(
UInt
j
=
0
;
j
<
n
;
j
++
)
{
if
(
surface
(
i
,
j
)
>
1e-10
)
area
++
;
}
}
return
area
/
Real
(
n
*
n
);
}
/* -------------------------------------------------------------------------- */
Surface
<
Real
>
SurfaceStatistics
::
computeAutocorrelation
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
correlation_heights
=
computePowerSpectrum
(
surface
);
correlation_heights
.
FFTITransform
();
return
correlation_heights
.
real
();
}
/* -------------------------------------------------------------------------- */
SurfaceComplex
<
Real
>
SurfaceStatistics
::
computePowerSpectrum
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
power_spectrum
(
surface
);
power_spectrum
.
FFTTransform
();
UInt
n
=
power_spectrum
.
size
();
Real
factor
=
1.
/
(
n
*
n
);
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
for
(
UInt
j
=
0
;
j
<
n
;
++
j
)
power_spectrum
(
i
,
j
)
*=
factor
;
power_spectrum
.
makeItRealBySquare
();
return
power_spectrum
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
T
SurfaceStatistics
::
computeAverage
(
const
Surface
<
T
>
&
surface
,
T
epsilon
){
Real
avg
=
0.0
;
UInt
cpt
=
0
;
UInt
n
=
surface
.
size
();
#pragma omp parallel for reduction(+: avg, cpt)
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
Real
val
=
surface
(
i
);
if
(
std
::
abs
(
val
)
>
epsilon
){
avg
+=
val
;
++
cpt
;
}
}
avg
/=
cpt
;
return
avg
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeSkewness
(
const
Surface
<
Real
>
&
surface
,
Real
avg
,
Real
std
){
Real
skw
=
0.0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
Real
tmp
=
surface
(
i
)
-
avg
;
skw
+=
(
tmp
*
tmp
*
tmp
);
}
skw
/=
(
n
*
std
*
std
*
std
);
return
skw
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeKurtosis
(
const
Surface
<
Real
>
&
surface
,
Real
avg
,
Real
std
){
Real
kur
=
0.0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
Real
tmp
=
surface
.
size
()
-
avg
;
kur
+=
(
tmp
*
tmp
*
tmp
*
tmp
);
}
kur
/=
(
n
*
std
*
std
*
std
*
std
);
return
kur
-
3
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeMaximum
(
const
Surface
<
Real
>
&
surface
){
Real
_max
=
-
1.
*
std
::
numeric_limits
<
Real
>::
max
();
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
_max
=
std
::
max
(
surface
(
i
),
_max
);
}
return
_max
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeMinimum
(
const
Surface
<
Real
>
&
surface
){
Real
_min
=
std
::
numeric_limits
<
Real
>::
max
();
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
_min
=
std
::
min
(
surface
(
i
),
_min
);
}
return
_min
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeRMSSlope
(
const
Surface
<
Real
>
&
surface
){
// UInt xl, yl;
UInt
xh
,
yh
;
Real
avg
,
slp
;
Real
gx
,
gy
;
UInt
n
=
surface
.
size
();
Surface
<
Real
>
grad
(
n
,
1.
);
for
(
UInt
i
=
0
,
m
=
0
;
i
<
n
;
i
++
)
{
for
(
UInt
j
=
0
;
j
<
n
;
j
++
)
{
UInt
k
=
i
+
n
*
j
;
// xl = ((i==0)?(k+n-1):(k-1));
xh
=
((
i
==
(
n
-
1
))
?
(
k
-
n
+
1
)
:
(
k
+
1
));
// yl = ((j==0)?(k+n*n-n):(k-n));
yh
=
((
j
==
(
n
-
1
))
?
(
k
-
n
*
n
+
n
)
:
(
k
+
n
));
gx
=
surface
(
xh
)
-
surface
(
k
);
gy
=
surface
(
yh
)
-
surface
(
k
);
grad
(
m
++
)
=
gx
*
gx
+
gy
*
gy
;
}
}
avg
=
SurfaceStatistics
::
computeAverage
(
grad
);
slp
=
sqrt
(
avg
);
/* According to paper Hyun PRE 70:026117 (2004) */
Real
scale_factor
=
surface
.
getL
()
/
n
;
return
slp
/
scale_factor
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeSpectralRMSSlope
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
power
=
computePowerSpectrum
(
surface
);
UInt
n
=
surface
.
size
();
Real
rms_slopes
=
0
;
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
)
for
(
UInt
j
=
1
;
j
<
n
/
2
;
++
j
){
rms_slopes
+=
1.
*
((
i
*
i
)
+
(
j
*
j
))
*
power
(
i
,
j
).
real
();
rms_slopes
+=
1.
*
((
i
*
i
)
+
(
j
*
j
))
*
power
(
n
-
i
,
n
-
j
).
real
();
rms_slopes
+=
1.
*
((
i
*
i
)
+
(
j
*
j
))
*
power
(
n
-
i
,
j
).
real
();
rms_slopes
+=
1.
*
((
i
*
i
)
+
(
j
*
j
))
*
power
(
i
,
n
-
j
).
real
();
}
// external line (i or j == 0)
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
){
rms_slopes
+=
1.
*
(
i
*
i
)
*
power
(
i
,
0
).
real
();
rms_slopes
+=
1.
*
(
i
*
i
)
*
power
(
n
-
i
,
0
).
real
();
rms_slopes
+=
1.
*
(
i
*
i
)
*
power
(
0
,
i
).
real
();
rms_slopes
+=
1.
*
(
i
*
i
)
*
power
(
0
,
n
-
i
).
real
();
}
//internal line (i or j == n/2)
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
)
{
rms_slopes
+=
0.5
*
(
i
*
i
+
n
*
n
/
4
)
*
power
(
i
,
n
/
2
).
real
();
rms_slopes
+=
0.5
*
(
i
*
i
+
n
*
n
/
4
)
*
power
(
n
-
i
,
n
/
2
).
real
();
rms_slopes
+=
0.5
*
(
i
*
i
+
n
*
n
/
4
)
*
power
(
n
/
2
,
i
).
real
();
rms_slopes
+=
0.5
*
(
i
*
i
+
n
*
n
/
4
)
*
power
(
n
/
2
,
n
-
i
).
real
();
}
{
//i == n/2 j == n/2
rms_slopes
+=
0.25
*
(
n
*
n
/
2
)
*
power
(
n
/
2
,
n
/
2
).
real
();
//i == 0 j == n/2
rms_slopes
+=
1.
*
(
n
*
n
/
4
)
*
power
(
0
,
n
/
2
).
real
();
//i == n/2 j == 0
rms_slopes
+=
1.
*
(
n
*
n
/
4
)
*
power
(
n
/
2
,
0
).
real
();
}
rms_slopes
=
2.
*
M_PI
/
n
*
sqrt
(
rms_slopes
);
Real
scale_factor
=
surface
.
getL
()
/
n
;
return
rms_slopes
/
scale_factor
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeStdev
(
const
Surface
<
Real
>
&
surface
,
Real
epsilon
){
UInt
n
=
surface
.
size
();
UInt
cpt
=
0
;
Real
std
=
0.0
;
Real
avg
=
0.0
;
for
(
UInt
i
=
0
;
i
<
n
*
n
;
i
++
)
{
Real
val
=
surface
(
i
);
if
(
val
>
epsilon
){
avg
+=
val
;
std
+=
val
*
val
;
++
cpt
;
}
}
avg
/=
cpt
;
return
sqrt
(
std
/
cpt
-
avg
*
avg
);
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeSpectralStdev
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
power
=
computePowerSpectrum
(
surface
);
Real
rms_heights
=
0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
for
(
UInt
j
=
0
;
j
<
n
;
++
j
){
rms_heights
+=
power
(
i
,
j
).
real
();
}
rms_heights
=
sqrt
(
rms_heights
);
return
rms_heights
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeSpectralMeanCurvature
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
power
=
computePowerSpectrum
(
surface
);
UInt
n
=
surface
.
size
();
Real
rms_curvatures
=
0
;
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
)
for
(
UInt
j
=
1
;
j
<
n
/
2
;
++
j
){
Real
coeff
=
((
i
*
i
)
+
(
j
*
j
));
coeff
*=
coeff
;
rms_curvatures
+=
1.
*
coeff
*
power
(
i
,
j
).
real
();
rms_curvatures
+=
1.
*
coeff
*
power
(
n
-
i
,
n
-
j
).
real
();
rms_curvatures
+=
1.
*
coeff
*
power
(
n
-
i
,
j
).
real
();
rms_curvatures
+=
1.
*
coeff
*
power
(
i
,
n
-
j
).
real
();
}
// external line (i or j == 0)
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
){
Real
coeff
=
i
*
i
;
coeff
*=
coeff
;
rms_curvatures
+=
4.
/
3.
*
coeff
*
power
(
i
,
0
).
real
();
rms_curvatures
+=
4.
/
3.
*
coeff
*
power
(
n
-
i
,
0
).
real
();
rms_curvatures
+=
4.
/
3.
*
coeff
*
power
(
0
,
i
).
real
();
rms_curvatures
+=
4.
/
3.
*
coeff
*
power
(
0
,
n
-
i
).
real
();
}
//internal line (i or j == n/2)
for
(
UInt
i
=
1
;
i
<
n
/
2
;
++
i
)
{
Real
coeff
=
(
i
*
i
+
n
*
n
/
4
);
coeff
*=
coeff
;
rms_curvatures
+=
0.5
*
coeff
*
power
(
i
,
n
/
2
).
real
();
rms_curvatures
+=
0.5
*
coeff
*
power
(
n
-
i
,
n
/
2
).
real
();
rms_curvatures
+=
0.5
*
coeff
*
power
(
n
/
2
,
i
).
real
();
rms_curvatures
+=
0.5
*
coeff
*
power
(
n
/
2
,
n
-
i
).
real
();
}
{
//i == n/2 j == n/2
rms_curvatures
+=
0.25
*
(
n
*
n
/
2
)
*
(
n
*
n
/
2
)
*
power
(
n
/
2
,
n
/
2
).
real
();
//i == 0 j == n/2
rms_curvatures
+=
2.
/
3.
*
(
n
*
n
/
4
)
*
(
n
*
n
/
4
)
*
power
(
0
,
n
/
2
).
real
();
//i == n/2 j == 0
rms_curvatures
+=
2.
/
3.
*
(
n
*
n
/
4
)
*
(
n
*
n
/
4
)
*
power
(
n
/
2
,
0
).
real
();
}
rms_curvatures
*=
9.
/
4.
;
rms_curvatures
=
M_PI
/
n
/
n
*
sqrt
(
rms_curvatures
);
return
rms_curvatures
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computePerimeter
(
const
Surface
<
Real
>
&
surface
,
Real
eps
)
{
int
perimeter
=
0
;
UInt
n
=
surface
.
size
();
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
n
;
++
j
)
{
if
(
j
>
0
&&
(
(
surface
(
i
,(
j
-
1
))
<
eps
&&
surface
(
i
,
j
)
>
eps
)
||
(
surface
(
i
,(
j
-
1
))
>
eps
&&
surface
(
i
,
j
)
<
eps
)
)
)
perimeter
++
;
if
(
i
>
0
&&
(
(
surface
(
j
,(
i
-
1
))
<
eps
&&
surface
(
j
,
i
)
>
eps
)
||
(
surface
(
j
,(
i
-
1
))
>
eps
&&
surface
(
j
,
i
)
<
eps
))
)
perimeter
++
;
}
}
return
perimeter
/
Real
(
n
);
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeSum
(
const
Surface
<
Real
>
&
surface
){
Real
sum
=
0.
;
UInt
n
=
surface
.
size
();
#pragma omp parallel for reduction(+: sum)
for
(
UInt
i
=
0
;
i
<
n
*
n
;
++
i
)
sum
+=
surface
(
i
);
return
sum
;
}
/* -------------------------------------------------------------------------- */
std
::
vector
<
Real
>
SurfaceStatistics
::
computeMoments
(
Surface
<
Real
>
&
surface
){
SurfaceComplex
<
Real
>
power
=
computePowerSpectrum
(
surface
);
Real
L
=
surface
.
getL
();
UInt
n
=
surface
.
size
();
Real
k_factor
=
2.
*
M_PI
/
L
;
std
::
vector
<
std
::
vector
<
Real
>
>
moment
;
moment
.
resize
(
6
);
for
(
int
i
=
0
;
i
<
5
;
i
++
)
moment
[
i
].
resize
(
5
);
for
(
UInt
i
=
0
;
i
<
n
/
2
;
i
++
)
{
Real
i0
=
pow
(
Real
(
i
)
*
k_factor
,
0
);
Real
i2
=
pow
(
Real
(
i
)
*
k_factor
,
2
);
Real
i4
=
pow
(
Real
(
i
)
*
k_factor
,
4
);
Real
im0
=
pow
(
-
Real
(
i
)
*
k_factor
,
0
);
Real
im2
=
pow
(
-
Real
(
i
)
*
k_factor
,
2
);
Real
im4
=
pow
(
-
Real
(
i
)
*
k_factor
,
4
);
for
(
UInt
j
=
0
;
j
<
n
/
2
;
j
++
)
{
Real
Phi
=
power
(
i
,
j
).
real
();
Real
Phi_10
;
Real
Phi_01
;
Real
Phi_11
;
if
(
i
==
0
&&
j
==
0
)
{
Phi_10
=
0
;
Phi_11
=
0
;
Phi_01
=
0
;
}
else
if
(
i
==
0
)
{
Phi_10
=
0
;
Phi_11
=
0
;
Phi_01
=
power
(
i
,(
n
-
j
)).
real
();
}
else
if
(
j
==
0
)
{
Phi_10
=
power
(
n
-
i
,
j
).
real
();
Phi_11
=
0
;
Phi_01
=
0
;
}
else
{
Phi_10
=
power
(
n
-
i
,
j
).
real
();
Phi_11
=
power
(
n
-
i
,(
n
-
j
)).
real
();
Phi_01
=
power
(
i
,(
n
-
j
)).
real
();
}
Real
j0
=
pow
(
Real
(
j
)
*
k_factor
,
0
);
Real
j2
=
pow
(
Real
(
j
)
*
k_factor
,
2
);
Real
j4
=
pow
(
Real
(
j
)
*
k_factor
,
4
);
Real
jm0
=
pow
(
-
Real
(
j
)
*
k_factor
,
0
);
Real
jm2
=
pow
(
-
Real
(
j
)
*
k_factor
,
2
);
Real
jm4
=
pow
(
-
Real
(
j
)
*
k_factor
,
4
);
moment
[
0
][
0
]
+=
i0
*
j0
*
Phi
+
im0
*
j0
*
Phi_10
+
i0
*
jm0
*
Phi_01
+
im0
*
jm0
*
Phi_11
;
moment
[
0
][
2
]
+=
i0
*
j2
*
Phi
+
im0
*
j2
*
Phi_10
+
i0
*
jm2
*
Phi_01
+
im0
*
jm2
*
Phi_11
;
moment
[
2
][
0
]
+=
i2
*
j0
*
Phi
+
im2
*
j0
*
Phi_10
+
i2
*
jm0
*
Phi_01
+
im2
*
jm0
*
Phi_11
;
moment
[
2
][
2
]
+=
i2
*
j2
*
Phi
+
im2
*
j2
*
Phi_10
+
i2
*
jm2
*
Phi_01
+
im2
*
jm2
*
Phi_11
;
moment
[
0
][
4
]
+=
i0
*
j4
*
Phi
+
im0
*
j4
*
Phi_10
+
i0
*
jm4
*
Phi_01
+
im0
*
jm4
*
Phi_11
;
moment
[
4
][
0
]
+=
i4
*
j0
*
Phi
+
im4
*
j0
*
Phi_10
+
i4
*
jm0
*
Phi_01
+
im4
*
jm0
*
Phi_11
;
}
}
std
::
vector
<
Real
>
mean_moments
;
//m2
Real
mean_m2
=
0.5
*
(
moment
[
0
][
2
]
+
moment
[
2
][
0
]);
mean_moments
.
push_back
(
mean_m2
);
//m4
Real
mean_m4
=
(
3
*
moment
[
2
][
2
]
+
moment
[
0
][
4
]
+
moment
[
4
][
0
])
/
3.
;
mean_moments
.
push_back
(
mean_m4
);
return
mean_moments
;
}
/* -------------------------------------------------------------------------- */
Real
SurfaceStatistics
::
computeAnalyticFractalMoment
(
UInt
order
,
UInt
k1
,
UInt
k2
,
Real
hurst
,
Real
A
,
Real
L
){
Real
xi
=
k2
/
k1
;
Real
exp
=
-
2.
*
hurst
;
Real
_k1
=
2.
*
M_PI
/
L
*
k1
;
Real
Cq
;
switch
(
order
){
case
0
:
Cq
=
2.
*
M_PI
;
break
;
case
2
:
Cq
=
M_PI
;
break
;
case
4
:
Cq
=
3.
/
4.
*
M_PI
;
break
;
default
:
SURFACE_FATAL
(
"invalid order"
);
}
exp
=
order
-
2.
*
hurst
;
Real
res
=
A
*
Cq
*
pow
(
_k1
,
exp
)
/
exp
*
(
pow
(
xi
,
exp
)
-
1
);
return
res
;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif
/* __SURFACE_STATISTICS_HH__ */
Event Timeline
Log In to Comment