Page MenuHomec4science
Files F61026097

surface_statistics.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 01:48

surface_statistics.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "surface_statistics.hh"
#include "fftw_engine.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace {
FFTWEngine engine; ///< local to this translation unit
}
/* -------------------------------------------------------------------------- */
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;
#pragma omp parallel for
for (int k = 0; k < N; ++k)
distrib_fft[k] = 0.;
#pragma omp parallel for
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<Real>(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<Real>(1., angle);
distrib_fft[N / 2] += phase;
}
#pragma omp parallel for
for (int k = 1; k < N; ++k)
distrib_fft[k] /= n * n;
distrib_fft[0] = std::polar(1., 0.);
#pragma omp parallel for
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
#pragma omp parallel for
for (int k = 0; k < N; ++k) {
distrib_fft[k] = conj(distrib_fft[k]);
}
auto* 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);
#pragma omp parallel for
for (int k = 0; k < N; ++k) {
distrib[k] = conj(distrib_fft[k]);
}
return distrib;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeAverageAmplitude(const Surface<Real>& surface) {
return surface.mean();
}
/* -------------------------------------------------------------------------- */
Complex SurfaceStatistics::computeAverageAmplitude(
const SurfaceComplex<Real>& surface) {
return surface.mean();
}
/* -------------------------------------------------------------------------- */
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();
#pragma omp parallel for collapse(2), reduction(+ : area)
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);
Surface<Real> acf(correlation_heights.size(), correlation_heights.getL());
engine.backward(acf, correlation_heights);
return acf;
}
/* -------------------------------------------------------------------------- */
SurfaceComplex<Real>
SurfaceStatistics::computePowerSpectrum(Surface<Real>& surface) {
SurfaceComplex<Real> power_spectrum(surface.size(), surface.getL());
engine.forward(surface, power_spectrum);
UInt n = power_spectrum.dataSize();
UInt factor = power_spectrum.size() * power_spectrum.size();
#pragma omp parallel for
for (UInt i = 0; i < n; ++i)
power_spectrum(i) /= factor;
power_spectrum.makeItRealBySquare();
return power_spectrum;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSkewness(const Surface<Real>& surface, Real avg,
Real std) {
Real skw = 0.0;
UInt n = surface.size();
#pragma omp parallel for reduction(+ : skw)
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();
#pragma omp parallel for reduction(+ : kur)
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();
#pragma omp parallel for reduction(max : _max)
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();
#pragma omp parallel for reduction(min : _min)
for (UInt i = 0; i < n * n; i++) {
_min = std::min(surface(i), _min);
}
return _min;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeRMSSlope(const Surface<Real>& surface) {
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));
UInt xh = ((i == (n - 1)) ? (k - n + 1) : (k + 1));
// yl = ((j==0)?(k+n*n-n):(k-n));
UInt yh = ((j == (n - 1)) ? (k - n * n + n) : (k + n));
Real gx = surface(xh) - surface(k);
Real gy = surface(yh) - surface(k);
grad(m++) = gx * gx + gy * gy;
}
}
Real avg = SurfaceStatistics::computeAverage(grad);
Real 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;
// /!\ commented code was causing race conditions with dummy variable in
// GridHermitian
#pragma omp parallel for collapse(2), reduction(+ : rms_slopes)
for (UInt i = 1; i < n / 2; ++i)
for (UInt j = 1; j < n / 2; ++j) {
rms_slopes += 2. * ((i * i) + (j * j)) * power(i, j).real();
// rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
rms_slopes += 2. * ((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)
#pragma omp parallel for reduction(+ : rms_slopes)
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 += 2. * (i * i) * power(0, i).real();
// rms_slopes += 1.*(i*i)*power(0,n-i).real();
}
// internal line (i or j == n/2)
#pragma omp parallel for reduction(+ : rms_slopes)
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 += 1. * (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 Grid<Real, 2>& surface,
Real epsilon) {
UInt n = surface.dataSize();
UInt cpt = 0;
Real std = 0.0;
Real avg = 0.0;
#pragma omp parallel for reduction(+ : avg, std, cpt)
for (UInt i = 0; i < 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();
#pragma omp parallel for collapse(2), reduction(+ : rms_heights)
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;
#pragma omp parallel for collapse(2), reduction(+ : rms_curvatures)
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)
#pragma omp parallel for reduction(+ : rms_curvatures)
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)
#pragma omp parallel for reduction(+ : rms_curvatures)
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();
#pragma omp parallel for collapse(2), reduction(+ : perimeter)
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 _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");
}
Real exp = order - 2. * hurst;
Real res = A * Cq * pow(_k1, exp) / exp * (pow(xi, exp) - 1);
return res;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas

Event Timeline

c4science · Help