source: cpp/frams/model/similarity/measure-distribution.cpp @ 1044

// This file is a part of Framsticks SDK.  http://www.framsticks.com/
// Copyright (C) 1999-2020  Maciej Komosinski and Szymon Ulatowski.
// See LICENSE.txt for details.

#include "measure-distribution.h"
#include <common/nonstd_math.h>
#include <limits>
#include "EMD/emd.c"
#include <iostream>

#define PI 3.1415926535

#define FIELDSTRUCT SimilMeasureDistribution

static ParamEntry simil_distribution_paramtab[] = {
                { "Similarity: distribution", 1, 4, "SimilMeasureDistribution", "Evaluates morphological dissimilarity using distribution measure.", },
                { "simil_density", 0, 0, "Density of surface sampling", "f 1 100 10", FIELD(density), "", },
                { "simil_bin_num", 0, 0, "Number of bins", "d 1 1000 128", FIELD(bin_num), "", },
                { "simil_samples_num", 0, 0, "Number of samples", "d 1 1048576 1048576", FIELD(samples_num), "", },
                { "evaluateDistance", 0, PARAM_DONTSAVE | PARAM_USERHIDDEN, "evaluate model dissimilarity", "p f(oGeno,oGeno)", PROCEDURE(p_evaldistance), "Calculates dissimilarity between two models created from Geno objects.", },
                { 0, },
};

#undef FIELDSTRUCT

SimilMeasureDistribution::SimilMeasureDistribution() : localpar(simil_distribution_paramtab, this)
{
        localpar.setDefault();
        SimilMeasureDistribution::distribution_fun = &SimilMeasureDistribution::D2; //D1 and D2 are the best descriptors
}

double SimilMeasureDistribution::getDistance()
{
        double dist = 0;
        for (int i = 0; i < 2; i++)
        {
                funs[i] = new std::pair<double, float>[bin_num]();
                for (int j = 0; j < bin_num; j++)
                        funs[i][j] = std::make_pair(0, 0);
        }
        
        for (int i = 0; i < 2; i++)
                sst_models[i] = new SolidsShapeTypeModel((*models[i]));
        
        SimilMeasureDistribution::calculateFuns();      
        dist = SimilMeasureDistribution::compareFuns();
        
        for (int i = 0; i < 2; i++)
        {
                SAFEDELETE(sst_models[i]);
                SAFEDELETEARRAY(funs[i]);
        }
        return dist;
}

int SimilMeasureDistribution::setParams(std::vector<double> params)
{
        for (unsigned int i = 0; i < params.size(); i++)
                if (params.at(i) <= 0)
                {
                        logPrintf("SimilDistributionMeasure", "setParams", LOG_ERROR, "Param values should be larger than 0.");
                        return -1;
                }
        
        density = params.at(0);
        bin_num = params.at(1);
        samples_num = params.at(2);
        
        return 0;
}

void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, Model &sampled)
{
        int size = sampled.getPartCount();
        int samples_taken = samples_num;

        //Check if total number of points pairs is smaller than samples number
        //but first, prevent exceeding int limits
        if (size < (int) sqrt((double) std::numeric_limits<int>::max()))
                samples_taken = min(samples_num, size*size);

        //Get sampled distribution
        std::vector<double> dist_vect = (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, &sampled);

        auto result = std::minmax_element(dist_vect.begin(), dist_vect.end());
        double min = *result.first;
        double max = *result.second;

        //Create histogram
        float hist[bin_num] = {};       
        int ind = 0;

        for (unsigned int j = 0; j < dist_vect.size(); j++)
        {
                ind = (int) std::floor((dist_vect.at(j)-min)*1/(max-min)*bin_num);
                if (ind <= (bin_num-1))
                        hist[ind] += 1;
                else if (ind == bin_num)
                        hist[bin_num-1] += 1;
        }

        //Create pairs
        for (int j = 0; j < bin_num; j++)
        {
                fun[j] = std::make_pair(min+(max-min)/bin_num*(j+0.5), hist[j]);
        }

        //Normalize
        float total_mass = 0;
        for (int j = 0; j < bin_num; j++)
    {
        total_mass += fun[j].second;
    }

        for (int j = 0; j < bin_num; j++)
                fun[j].second /= total_mass;
}

void SimilMeasureDistribution::calculateFuns()
{
        for (int i = 0; i < 2; i++)
        {
                Model sampled = SimilMeasureDistribution::sampleSurface(&sst_models[i]->getModel(), density);
                SimilMeasureDistribution::calculateFun(funs[i], sampled);
        }
}

double SimilMeasureDistribution::compareFuns()
{
        return SimilMeasureDistribution::EMD(funs[0], funs[1]);
}

vector<double> SimilMeasureDistribution::D1(int samples_taken, Model *sampled)
{
        vector<double> dist_vect;
        int size = sampled->getPartCount();
        double x = 0;
        double y = 0;
        double z = 0;
        
        for (int i = 0; i < size; i++)
        {
                Pt3D pos = sampled->getPart(i)->p;
                x += pos.x;
                y += pos.y;
                z += pos.z;
        }
        
        x = x/size;
        y = y/size;
        z = z/size;
        
        Pt3D centroid = {x, y, z};
        
        for (int i = 0; i < samples_taken; i++)
        {
                int p1 = rndUint(size);
                double dist = sampled->getPart(p1)->p.distanceTo(centroid);
                if (dist > 0)
                        dist_vect.push_back(dist);
        }
        
        return dist_vect;
}

vector<double> SimilMeasureDistribution::D2(int samples_taken, Model *sampled)
{
        vector<double> dist_vect;
        int size = sampled->getPartCount();
        for (int i = 0; i < samples_taken; i++)
        {
                int p1 = rndUint(size);
                int p2 = rndUint(size);
                double dist = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p);
                if (dist > 0)
                        dist_vect.push_back(dist);
        }
        
        return dist_vect;
}

vector<double> SimilMeasureDistribution::D3(int samples_taken, Model *sampled)
{
        vector<double> dist_vect;
        int size = sampled->getPartCount();
        for (int i = 0; i < samples_taken; i++)
        {
                int p1 = rndUint(size);
                int p2 = rndUint(size);
                int p3 = rndUint(size);
                
                Pt3D v(sampled->getPart(p2)->p);
                Pt3D w(sampled->getPart(p3)->p);
                v -= sampled->getPart(p1)->p;
                w -= sampled->getPart(p1)->p;
                Pt3D cross_prod(0);
                cross_prod.vectorProduct(v, w);
                
                double dist = 0.5 * cross_prod.length();
                if (dist > 0)
                        dist_vect.push_back(dist);
        }
        
        return dist_vect;
}

vector<double> SimilMeasureDistribution::D4(int samples_taken, Model *sampled)
{
        vector<double> dist_vect;
        int size = sampled->getPartCount();
        for (int i = 0; i < samples_taken; i++)
        {
                int a = rndUint(size);
                int b = rndUint(size);
                int c = rndUint(size);
                int d = rndUint(size);
                
                Pt3D ad(sampled->getPart(a)->p);
                Pt3D bd(sampled->getPart(b)->p);
                Pt3D cd(sampled->getPart(c)->p);
                
                ad -= sampled->getPart(d)->p;
                bd -= sampled->getPart(d)->p;
                cd -= sampled->getPart(d)->p;
                
                Pt3D cross_prod(0);
                cross_prod.vectorProduct(bd, cd);
                cross_prod.entrywiseProduct(ad);
                
                double dist = cross_prod.length()/6;
                if (dist > 0)
                        dist_vect.push_back(dist);
        }
        
        return dist_vect;
}

vector<double> SimilMeasureDistribution::A3(int samples_taken, Model *sampled)
{
        vector<double> dist_vect;
        int size = sampled->getPartCount();
        for (int i = 0; i < samples_taken; i++)
        {
                int p1 = rndUint(size);
                int p2 = rndUint(size);
                int p3 = rndUint(size);
                double a = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p3)->p);
                double b = sampled->getPart(p3)->p.distanceTo(sampled->getPart(p2)->p);
                double c = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p);
                double beta = acos((a*a + b*b - c*c)/(2*a*b));
        
                if (!std::isnan(beta))
                        dist_vect.push_back(beta);
        }
        
        return dist_vect;
}


float dist(feature_t* F1, feature_t* F2)
{ 
        return abs((*F1)-(*F2)); 
}


void SimilMeasureDistribution::fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights)
{
        for (int j = 0; j < bin_num; j++)
        {
                points[j] = {fun[j].first};
                weights[j] = fun[j].second;
        }
}

double SimilMeasureDistribution::EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2)
{
        feature_t *points[2];
        float *weights[2];
        
        for (int i = 0; i < 2; i++)
        {
                points[i] = new feature_t[bin_num];
                weights[i] = new float[bin_num]();
        }
        SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]);
        SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]);
        
        signature_t sig1 = {bin_num, points[0], weights[0]},
        sig2 = {bin_num, points[1], weights[1]};

        float e = emd(&sig1, &sig2, dist, 0, 0);
        
        for (int i = 0; i < 2; i++)
        {
                delete[] points[i];
                delete[] weights[i];
        }

        return e;
}