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

Last change on this file since 1121 was 1121, checked in by Maciej Komosinski, 9 days ago

Added comments, removed unused variables

File size: 8.6 KB
Line 
1// This file is a part of Framsticks SDK.  http://www.framsticks.com/
2// Copyright (C) 1999-2021  Maciej Komosinski and Szymon Ulatowski.
3// See LICENSE.txt for details.
4
5#include "measure-distribution.h"
6#include <common/nonstd_math.h>
7#include <limits>
8#include "EMD/emd.c"
9#include <iostream>
10
11#define FIELDSTRUCT SimilMeasureDistribution
12
13static ParamEntry simil_distribution_paramtab[] = {
14                { "Creature: Similarity: Descriptor distribution", 1, 4, "SimilMeasureDistribution", "Evaluates morphological dissimilarity using distribution measure.", },
15                { "simil_density", 0, 0, "Density of surface sampling", "f 1 100 10", FIELD(density), "", },
16                { "simil_bin_num", 0, 0, "Number of bins", "d 1 1000 128", FIELD(bin_num), "", },
17                { "simil_samples_num", 0, 0, "Number of samples", "d 1 1048576 10000", FIELD(samples_num), "", }, //based on experiments, not much accuracy to gain when this is increased above 1000
18                { "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.", },
19                { 0, },
20};
21
22#undef FIELDSTRUCT
23
24SimilMeasureDistribution::SimilMeasureDistribution() : localpar(simil_distribution_paramtab, this)
25{
26        localpar.setDefault();
27        SimilMeasureDistribution::distribution_fun = &SimilMeasureDistribution::D2; //D1 and D2 are the best descriptors
28}
29
30double SimilMeasureDistribution::getDistance()
31{
32        double dist = 0;
33        for (int i = 0; i < 2; i++)
34        {
35                funs[i] = new std::pair<double, float>[bin_num]();
36                for (int j = 0; j < bin_num; j++)
37                        funs[i][j] = std::make_pair(0, 0);
38        }
39
40        for (int i = 0; i < 2; i++)
41                sst_models[i] = new SolidsShapeTypeModel((*models[i]));
42
43        SimilMeasureDistribution::calculateFuns();
44        dist = SimilMeasureDistribution::compareFuns();
45
46        for (int i = 0; i < 2; i++)
47        {
48                SAFEDELETE(sst_models[i]);
49                SAFEDELETEARRAY(funs[i]);
50        }
51        return dist;
52}
53
54int SimilMeasureDistribution::setParams(std::vector<double> params)
55{
56        for (unsigned int i = 0; i < params.size(); i++)
57                if (params.at(i) <= 0)
58                {
59                        logPrintf("SimilDistributionMeasure", "setParams", LOG_ERROR, "Param values should be larger than 0.");
60                        return -1;
61                }
62
63        density = params.at(0);
64        bin_num = params.at(1);
65        samples_num = params.at(2);
66
67        return 0;
68}
69
70void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, const Model &sampled)
71{
72        int samples_taken = samples_num;
73
74        //Check if total number of point pairs is smaller than samples number (just to avoid the calculation of the same stats for the same pairs of Parts)
75        //This optimization turned out to have a minor effect, only present for very high simil_samples_num, and was not perfect anyway:
76        //- samples are selected randomly so there are no guarantees that they will be repeated,
77        //- even if they do, it has the benefit of averaging the result that becomes more stable,
78        //- the concept of "point pairs" is not relevant when we randomly select more than two points, as some descriptors do.
79        //int size = sampled.getPartCount();
80        //if (size < (int) sqrt((double) std::numeric_limits<int>::max())) //prevent exceeding int limits
81        //      samples_taken = min(samples_num, size*size);
82
83        rndgen.seed(55); //For determinism. Otherwise the descriptors (that choose samples pseudo-randomly) for the same Model can yield different values and so the dissimilarity between the object and its copy will not be 0.
84        std::uniform_int_distribution<> uniform_distrib(0, sampled.getPartCount() - 1);
85
86        //Get sampled distribution
87        std::vector<double> dist_vect;
88        dist_vect.reserve(samples_taken); //we will add up to samples_taken elements to this vector
89        (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, uniform_distrib, sampled, dist_vect);
90
91        auto result = std::minmax_element(dist_vect.begin(), dist_vect.end());
92        double min = *result.first;
93        double max = *result.second;
94
95        //Create histogram
96        vector<int> hist(bin_num);
97        int ind = 0;
98
99        for (unsigned int j = 0; j < dist_vect.size(); j++)
100        {
101                ind = (int)std::floor((dist_vect.at(j) - min) * 1 / (max - min) * bin_num);
102                if (ind <= (bin_num - 1))
103                        hist[ind]++;
104                else if (ind == bin_num)
105                        hist[bin_num - 1]++;
106        }
107
108        //Create pairs
109        for (int j = 0; j < bin_num; j++)
110        {
111                fun[j] = std::make_pair(min + (max - min) / bin_num * (j + 0.5), hist[j]);
112        }
113
114        //Normalize
115        float total_mass = 0;
116        for (int j = 0; j < bin_num; j++)
117        {
118                total_mass += fun[j].second;
119        }
120
121        for (int j = 0; j < bin_num; j++)
122                fun[j].second /= total_mass;
123}
124
125void SimilMeasureDistribution::calculateFuns()
126{
127        for (int i = 0; i < 2; i++)
128        {
129                Model sampled = SimilMeasureDistribution::sampleSurface(&sst_models[i]->getModel(), density);
130                SimilMeasureDistribution::calculateFun(funs[i], sampled);
131        }
132}
133
134double SimilMeasureDistribution::compareFuns()
135{
136        return SimilMeasureDistribution::EMD(funs[0], funs[1]);
137}
138
139void SimilMeasureDistribution::D1(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
140{
141        int size = sampled.getPartCount();
142        double x = 0;
143        double y = 0;
144        double z = 0;
145
146        for (int i = 0; i < size; i++)
147        {
148                Pt3D pos = sampled.getPart(i)->p;
149                x += pos.x;
150                y += pos.y;
151                z += pos.z;
152        }
153
154        x = x / size;
155        y = y / size;
156        z = z / size;
157
158        Pt3D centroid = { x, y, z };
159
160        for (int i = 0; i < samples_taken; i++)
161        {
162                int p1 = distribution(rndgen);
163                double dist = sampled.getPart(p1)->p.distanceTo(centroid);
164                if (dist > 0)
165                        dist_vect.push_back(dist);
166        }
167}
168
169void SimilMeasureDistribution::D2(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
170{
171        for (int i = 0; i < samples_taken; i++)
172        {
173                int p1 = distribution(rndgen);
174                int p2 = distribution(rndgen);
175                double dist = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p);
176                if (dist > 0)
177                        dist_vect.push_back(dist);
178        }
179}
180
181void SimilMeasureDistribution::D3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
182{
183        for (int i = 0; i < samples_taken; i++)
184        {
185                int p1 = distribution(rndgen);
186                int p2 = distribution(rndgen);
187                int p3 = distribution(rndgen);
188
189                Pt3D v(sampled.getPart(p2)->p);
190                Pt3D w(sampled.getPart(p3)->p);
191                v -= sampled.getPart(p1)->p;
192                w -= sampled.getPart(p1)->p;
193                Pt3D cross_prod(0);
194                cross_prod.vectorProduct(v, w);
195
196                double dist = 0.5 * cross_prod.length();
197                if (dist > 0)
198                        dist_vect.push_back(dist);
199        }
200}
201
202void SimilMeasureDistribution::D4(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
203{
204        for (int i = 0; i < samples_taken; i++)
205        {
206                int a = distribution(rndgen);
207                int b = distribution(rndgen);
208                int c = distribution(rndgen);
209                int d = distribution(rndgen);
210
211                Pt3D ad(sampled.getPart(a)->p);
212                Pt3D bd(sampled.getPart(b)->p);
213                Pt3D cd(sampled.getPart(c)->p);
214
215                ad -= sampled.getPart(d)->p;
216                bd -= sampled.getPart(d)->p;
217                cd -= sampled.getPart(d)->p;
218
219                Pt3D cross_prod(0);
220                cross_prod.vectorProduct(bd, cd);
221                cross_prod.entrywiseProduct(ad);
222
223                double dist = cross_prod.length() / 6;
224                if (dist > 0)
225                        dist_vect.push_back(dist);
226        }
227}
228
229void SimilMeasureDistribution::A3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
230{
231        for (int i = 0; i < samples_taken; i++)
232        {
233                int p1 = distribution(rndgen);
234                int p2 = distribution(rndgen);
235                int p3 = distribution(rndgen);
236                double a = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p3)->p);
237                double b = sampled.getPart(p3)->p.distanceTo(sampled.getPart(p2)->p);
238                double c = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p);
239                double beta = acos((a * a + b * b - c * c) / (2 * a * b));
240
241                if (!std::isnan(beta))
242                        dist_vect.push_back(beta);
243        }
244}
245
246
247float dist(feature_t* F1, feature_t* F2)
248{
249        return abs((*F1) - (*F2));
250}
251
252
253void SimilMeasureDistribution::fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights)
254{
255        for (int j = 0; j < bin_num; j++)
256        {
257                points[j] = { fun[j].first };
258                weights[j] = fun[j].second;
259        }
260}
261
262double SimilMeasureDistribution::EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2)
263{
264        feature_t *points[2];
265        float *weights[2];
266
267        for (int i = 0; i < 2; i++)
268        {
269                points[i] = new feature_t[bin_num];
270                weights[i] = new float[bin_num]();
271        }
272        SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]);
273        SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]);
274
275        signature_t sig1 = { bin_num, points[0], weights[0] },
276                sig2 = { bin_num, points[1], weights[1] };
277
278        float e = emd(&sig1, &sig2, dist, 0, 0, bin_num, bin_num);
279
280        for (int i = 0; i < 2; i++)
281        {
282                delete[] points[i];
283                delete[] weights[i];
284        }
285
286        return e;
287}
Note: See TracBrowser for help on using the repository browser.