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 | |
---|
13 | static 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 | |
---|
24 | SimilMeasureDistribution::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 | |
---|
30 | double 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 | |
---|
54 | int 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 | |
---|
70 | void 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 points). |
---|
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 = std::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 | |
---|
125 | void 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 | |
---|
134 | double SimilMeasureDistribution::compareFuns() |
---|
135 | { |
---|
136 | return SimilMeasureDistribution::EMD(funs[0], funs[1]); |
---|
137 | } |
---|
138 | |
---|
139 | void 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 | |
---|
169 | void 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 | |
---|
181 | void 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 | |
---|
202 | void 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 | |
---|
229 | void 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 | |
---|
247 | float dist(feature_t* F1, feature_t* F2) |
---|
248 | { |
---|
249 | return abs((*F1) - (*F2)); |
---|
250 | } |
---|
251 | |
---|
252 | |
---|
253 | void 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 | |
---|
262 | double 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 | } |
---|