misc.h

Go to the documentation of this file.

/* misc.cpp :  various  numerical functions
 * Copyright (C) 2020  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU 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 General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef ATOMPROBE_MATHS_MISC_H
#define ATOMPROBE_MATHS_MISC_H

#include "gsl/gsl_matrix.h"

#include <numeric>
#include <vector>
#include <limits>
#include <cmath>

#include "atomprobe/primitives/point3D.h"

namespace AtomProbe
{

enum PointDir{  POINTDIR_TOGETHER =0,
                POINTDIR_IN_COMMON,
                POINTDIR_APART
             };

//Multiply two matrices naively. If alloc=true, then the input matrix is overwritten, and resultant matrix must be deallocated with gsl_matrix_free().
void gsl_matrix_mult(const gsl_matrix *A, const gsl_matrix *B, gsl_matrix* &res,
    bool alloc=false);

// Uses SVD internally. Returns <0 on error.
unsigned int estimateRank(const gsl_matrix *m, 
        float tolerance=sqrt(std::numeric_limits<float>::epsilon()));

// returns false if constrained. Note that input matrix M is overwritten
bool solveLeastSquares(gsl_matrix *m, gsl_vector *b, gsl_vector* &x);

template<class T>
T weightedMean(const std::vector<T> &values, const std::vector<T> &weight)
{
    //ASSERT(weight.size() == values.size());
    if(values.empty())
        return 0.0f;
    
    //sum of (value[]*weight[])
    T sumWiXi=(T)0;
    // sum of(weight[])
    T wTotal=(T)0;
    wTotal= std::accumulate(weight.begin(),weight.end(),(T)0);
    for(size_t ui=0;ui<values.size();ui++)
    {
        //ASSERT(weight[ui] >=0.0f);
        sumWiXi+=values[ui]*weight[ui];
    }

    if(wTotal == 0.0f)
        return 0.0f;

    return sumWiXi/wTotal;
}

template<typename T>
void meanAndStdev(const std::vector<T> &f,float &meanVal, 
            float &stdevVal,bool normalCorrection=true)
{
    double meanDbl=0;
    for(size_t ui=0;ui<f.size();ui++)
        meanDbl+=f[ui];
    meanVal =meanDbl/f.size();

    
    double stdevDbl=0;
    for(size_t ui=0;ui<f.size();ui++)
    {
        float delta;
        delta=f[ui]-meanVal;
        stdevDbl+=delta*delta;
    }
    stdevVal = sqrtf( stdevDbl/float(f.size()-1));

    //Perform bias correction, assuming the input data is normally distributed
    if(normalCorrection)
    {
        float n=f.size();
        //Approximation to C4 = sqrt(2/(n-1))*gamma(n/2)/gamma((n-1)/2)
        // multiplier must be applied to 1/(n-1) normalised standard deviation
        // Citation: 
        //  http://support.sas.com/documentation/cdl/en/qcug/63922/HTML/default/qcug_functions_sect007.htm
        //  https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation
        stdevVal*=(1.0 - 1.0/(4.0*n) - 7.0/(32.0*n*n) - 19.0/(128.0*n*n*n));
    }
}


unsigned int vectorPointDir(const Point3D &pA, const Point3D &pB, 
                const Point3D &vC, const Point3D &vD);

//Distance between a line segment and a point in 3D space
float distanceToSegment(const Point3D &fA, const Point3D &fB, const Point3D &p);


//TODO: Convert to vector input?
float signedDistanceToFacet(const Point3D &fA, const Point3D &fB, 
            const Point3D &fC,  const Point3D &normal,const Point3D &p);

float distanceToFacet(const Point3D &fA, const Point3D &fB, 
            const Point3D &fC, const Point3D &normal,const Point3D &p);


inline float dotProduct(float a1, float a2, float a3, 
            float b1, float b2, float b3)
{
    return a1*b1 + a2*b2 + a3* b3;
}


//Calculate the determinant of a 3x3 matrix, C-order layout.
double det3by3(const double *ptArray);

//Determines the volume of a quadrilateral pyramid
//input points "planarpts" must be adjacent (connected) by 
//0 <-> 1 <-> 2 <-> 0, all points connected to apex
double pyramidVol(const Point3D *planarPts, const Point3D &apex);

//Return strue if a triangle is degenerate (has colinear points)
bool triIsDegenerate(const Point3D &fA, const Point3D &fB, 
            const Point3D &fC);

}

#endif