Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/surface_statistics.hh
e96d4bebeabcfeatures/netcdf-mpi

Last Change
Raw File

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__ */
c4science · Help