Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65802269
surface_statistics.cpp
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, Jun 6, 08:15
Size
16 KB
Mime Type
text/x-c
Expires
Sat, Jun 8, 08:15 (2 d)
Engine
blob
Format
Raw Data
Handle
18129444
Attached To
rTAMAAS tamaas
surface_statistics.cpp
View Options
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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
Log In to Comment