1 | // This file is a part of Framsticks SDK. http://www.framsticks.com/ |
---|
2 | // Copyright (C) 1999-2020 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 PI 3.1415926535 |
---|
12 | |
---|
13 | #define FIELDSTRUCT SimilMeasureDistribution |
---|
14 | |
---|
15 | static ParamEntry simil_distribution_paramtab[] = { |
---|
16 | { "Similarity: distribution", 1, 4, "SimilMeasureDistribution", "Evaluates morphological dissimilarity using distribution measure.", }, |
---|
17 | { "simil_density", 0, 0, "Density of surface sampling", "f 1 100 10", FIELD(density), "", }, |
---|
18 | { "simil_bin_num", 0, 0, "Number of bins", "d 1 1000 128", FIELD(bin_num), "", }, |
---|
19 | { "simil_samples_num", 0, 0, "Number of samples", "d 1 1048576 1048576", FIELD(samples_num), "", }, |
---|
20 | { "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.", }, |
---|
21 | { 0, }, |
---|
22 | }; |
---|
23 | |
---|
24 | #undef FIELDSTRUCT |
---|
25 | |
---|
26 | SimilMeasureDistribution::SimilMeasureDistribution() : localpar(simil_distribution_paramtab, this) |
---|
27 | { |
---|
28 | localpar.setDefault(); |
---|
29 | SimilMeasureDistribution::distribution_fun = &SimilMeasureDistribution::D2; //D1 and D2 are the best descriptors |
---|
30 | } |
---|
31 | |
---|
32 | double SimilMeasureDistribution::getDistance() |
---|
33 | { |
---|
34 | double dist = 0; |
---|
35 | for (int i = 0; i < 2; i++) |
---|
36 | { |
---|
37 | funs[i] = new std::pair<double, float>[bin_num](); |
---|
38 | for (int j = 0; j < bin_num; j++) |
---|
39 | funs[i][j] = std::make_pair(0, 0); |
---|
40 | } |
---|
41 | |
---|
42 | for (int i = 0; i < 2; i++) |
---|
43 | sst_models[i] = new SolidsShapeTypeModel((*models[i])); |
---|
44 | |
---|
45 | SimilMeasureDistribution::calculateFuns(); |
---|
46 | dist = SimilMeasureDistribution::compareFuns(); |
---|
47 | |
---|
48 | for (int i = 0; i < 2; i++) |
---|
49 | { |
---|
50 | SAFEDELETE(sst_models[i]); |
---|
51 | SAFEDELETEARRAY(funs[i]); |
---|
52 | } |
---|
53 | return dist; |
---|
54 | } |
---|
55 | |
---|
56 | int SimilMeasureDistribution::setParams(std::vector<double> params) |
---|
57 | { |
---|
58 | for (unsigned int i = 0; i < params.size(); i++) |
---|
59 | if (params.at(i) <= 0) |
---|
60 | { |
---|
61 | logPrintf("SimilDistributionMeasure", "setParams", LOG_ERROR, "Param values should be larger than 0."); |
---|
62 | return -1; |
---|
63 | } |
---|
64 | |
---|
65 | density = params.at(0); |
---|
66 | bin_num = params.at(1); |
---|
67 | samples_num = params.at(2); |
---|
68 | |
---|
69 | return 0; |
---|
70 | } |
---|
71 | |
---|
72 | void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, Model &sampled) |
---|
73 | { |
---|
74 | int size = sampled.getPartCount(); |
---|
75 | int samples_taken = samples_num; |
---|
76 | |
---|
77 | //Check if total number of points pairs is smaller than samples number |
---|
78 | //but first, prevent exceeding int limits |
---|
79 | if (size < (int) sqrt((double) std::numeric_limits<int>::max())) |
---|
80 | samples_taken = min(samples_num, size*size); |
---|
81 | |
---|
82 | //Get sampled distribution |
---|
83 | std::vector<double> dist_vect = (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, &sampled); |
---|
84 | |
---|
85 | auto result = std::minmax_element(dist_vect.begin(), dist_vect.end()); |
---|
86 | double min = *result.first; |
---|
87 | double max = *result.second; |
---|
88 | |
---|
89 | //Create histogram |
---|
90 | float hist[bin_num] = {}; |
---|
91 | int ind = 0; |
---|
92 | |
---|
93 | for (unsigned int j = 0; j < dist_vect.size(); j++) |
---|
94 | { |
---|
95 | ind = (int) std::floor((dist_vect.at(j)-min)*1/(max-min)*bin_num); |
---|
96 | if (ind <= (bin_num-1)) |
---|
97 | hist[ind] += 1; |
---|
98 | else if (ind == bin_num) |
---|
99 | hist[bin_num-1] += 1; |
---|
100 | } |
---|
101 | |
---|
102 | //Create pairs |
---|
103 | for (int j = 0; j < bin_num; j++) |
---|
104 | { |
---|
105 | fun[j] = std::make_pair(min+(max-min)/bin_num*(j+0.5), hist[j]); |
---|
106 | } |
---|
107 | |
---|
108 | //Normalize |
---|
109 | float total_mass = 0; |
---|
110 | for (int j = 0; j < bin_num; j++) |
---|
111 | { |
---|
112 | total_mass += fun[j].second; |
---|
113 | } |
---|
114 | |
---|
115 | for (int j = 0; j < bin_num; j++) |
---|
116 | fun[j].second /= total_mass; |
---|
117 | } |
---|
118 | |
---|
119 | void SimilMeasureDistribution::calculateFuns() |
---|
120 | { |
---|
121 | for (int i = 0; i < 2; i++) |
---|
122 | { |
---|
123 | Model sampled = SimilMeasureDistribution::sampleSurface(&sst_models[i]->getModel(), density); |
---|
124 | SimilMeasureDistribution::calculateFun(funs[i], sampled); |
---|
125 | } |
---|
126 | } |
---|
127 | |
---|
128 | double SimilMeasureDistribution::compareFuns() |
---|
129 | { |
---|
130 | return SimilMeasureDistribution::EMD(funs[0], funs[1]); |
---|
131 | } |
---|
132 | |
---|
133 | vector<double> SimilMeasureDistribution::D1(int samples_taken, Model *sampled) |
---|
134 | { |
---|
135 | vector<double> dist_vect; |
---|
136 | int size = sampled->getPartCount(); |
---|
137 | double x = 0; |
---|
138 | double y = 0; |
---|
139 | double z = 0; |
---|
140 | |
---|
141 | for (int i = 0; i < size; i++) |
---|
142 | { |
---|
143 | Pt3D pos = sampled->getPart(i)->p; |
---|
144 | x += pos.x; |
---|
145 | y += pos.y; |
---|
146 | z += pos.z; |
---|
147 | } |
---|
148 | |
---|
149 | x = x/size; |
---|
150 | y = y/size; |
---|
151 | z = z/size; |
---|
152 | |
---|
153 | Pt3D centroid = {x, y, z}; |
---|
154 | |
---|
155 | for (int i = 0; i < samples_taken; i++) |
---|
156 | { |
---|
157 | int p1 = rndUint(size); |
---|
158 | double dist = sampled->getPart(p1)->p.distanceTo(centroid); |
---|
159 | if (dist > 0) |
---|
160 | dist_vect.push_back(dist); |
---|
161 | } |
---|
162 | |
---|
163 | return dist_vect; |
---|
164 | } |
---|
165 | |
---|
166 | vector<double> SimilMeasureDistribution::D2(int samples_taken, Model *sampled) |
---|
167 | { |
---|
168 | vector<double> dist_vect; |
---|
169 | int size = sampled->getPartCount(); |
---|
170 | for (int i = 0; i < samples_taken; i++) |
---|
171 | { |
---|
172 | int p1 = rndUint(size); |
---|
173 | int p2 = rndUint(size); |
---|
174 | double dist = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p); |
---|
175 | if (dist > 0) |
---|
176 | dist_vect.push_back(dist); |
---|
177 | } |
---|
178 | |
---|
179 | return dist_vect; |
---|
180 | } |
---|
181 | |
---|
182 | vector<double> SimilMeasureDistribution::D3(int samples_taken, Model *sampled) |
---|
183 | { |
---|
184 | vector<double> dist_vect; |
---|
185 | int size = sampled->getPartCount(); |
---|
186 | for (int i = 0; i < samples_taken; i++) |
---|
187 | { |
---|
188 | int p1 = rndUint(size); |
---|
189 | int p2 = rndUint(size); |
---|
190 | int p3 = rndUint(size); |
---|
191 | |
---|
192 | Pt3D v(sampled->getPart(p2)->p); |
---|
193 | Pt3D w(sampled->getPart(p3)->p); |
---|
194 | v -= sampled->getPart(p1)->p; |
---|
195 | w -= sampled->getPart(p1)->p; |
---|
196 | Pt3D cross_prod(0); |
---|
197 | cross_prod.vectorProduct(v, w); |
---|
198 | |
---|
199 | double dist = 0.5 * cross_prod.length(); |
---|
200 | if (dist > 0) |
---|
201 | dist_vect.push_back(dist); |
---|
202 | } |
---|
203 | |
---|
204 | return dist_vect; |
---|
205 | } |
---|
206 | |
---|
207 | vector<double> SimilMeasureDistribution::D4(int samples_taken, Model *sampled) |
---|
208 | { |
---|
209 | vector<double> dist_vect; |
---|
210 | int size = sampled->getPartCount(); |
---|
211 | for (int i = 0; i < samples_taken; i++) |
---|
212 | { |
---|
213 | int a = rndUint(size); |
---|
214 | int b = rndUint(size); |
---|
215 | int c = rndUint(size); |
---|
216 | int d = rndUint(size); |
---|
217 | |
---|
218 | Pt3D ad(sampled->getPart(a)->p); |
---|
219 | Pt3D bd(sampled->getPart(b)->p); |
---|
220 | Pt3D cd(sampled->getPart(c)->p); |
---|
221 | |
---|
222 | ad -= sampled->getPart(d)->p; |
---|
223 | bd -= sampled->getPart(d)->p; |
---|
224 | cd -= sampled->getPart(d)->p; |
---|
225 | |
---|
226 | Pt3D cross_prod(0); |
---|
227 | cross_prod.vectorProduct(bd, cd); |
---|
228 | cross_prod.entrywiseProduct(ad); |
---|
229 | |
---|
230 | double dist = cross_prod.length()/6; |
---|
231 | if (dist > 0) |
---|
232 | dist_vect.push_back(dist); |
---|
233 | } |
---|
234 | |
---|
235 | return dist_vect; |
---|
236 | } |
---|
237 | |
---|
238 | vector<double> SimilMeasureDistribution::A3(int samples_taken, Model *sampled) |
---|
239 | { |
---|
240 | vector<double> dist_vect; |
---|
241 | int size = sampled->getPartCount(); |
---|
242 | for (int i = 0; i < samples_taken; i++) |
---|
243 | { |
---|
244 | int p1 = rndUint(size); |
---|
245 | int p2 = rndUint(size); |
---|
246 | int p3 = rndUint(size); |
---|
247 | double a = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p3)->p); |
---|
248 | double b = sampled->getPart(p3)->p.distanceTo(sampled->getPart(p2)->p); |
---|
249 | double c = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p); |
---|
250 | double beta = acos((a*a + b*b - c*c)/(2*a*b)); |
---|
251 | |
---|
252 | if (!std::isnan(beta)) |
---|
253 | dist_vect.push_back(beta); |
---|
254 | } |
---|
255 | |
---|
256 | return dist_vect; |
---|
257 | } |
---|
258 | |
---|
259 | |
---|
260 | float dist(feature_t* F1, feature_t* F2) |
---|
261 | { |
---|
262 | return abs((*F1)-(*F2)); |
---|
263 | } |
---|
264 | |
---|
265 | |
---|
266 | void SimilMeasureDistribution::fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights) |
---|
267 | { |
---|
268 | for (int j = 0; j < bin_num; j++) |
---|
269 | { |
---|
270 | points[j] = {fun[j].first}; |
---|
271 | weights[j] = fun[j].second; |
---|
272 | } |
---|
273 | } |
---|
274 | |
---|
275 | double SimilMeasureDistribution::EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2) |
---|
276 | { |
---|
277 | feature_t *points[2]; |
---|
278 | float *weights[2]; |
---|
279 | |
---|
280 | for (int i = 0; i < 2; i++) |
---|
281 | { |
---|
282 | points[i] = new feature_t[bin_num]; |
---|
283 | weights[i] = new float[bin_num](); |
---|
284 | } |
---|
285 | SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]); |
---|
286 | SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]); |
---|
287 | |
---|
288 | signature_t sig1 = {bin_num, points[0], weights[0]}, |
---|
289 | sig2 = {bin_num, points[1], weights[1]}; |
---|
290 | |
---|
291 | float e = emd(&sig1, &sig2, dist, 0, 0); |
---|
292 | |
---|
293 | for (int i = 0; i < 2; i++) |
---|
294 | { |
---|
295 | delete[] points[i]; |
---|
296 | delete[] weights[i]; |
---|
297 | } |
---|
298 | |
---|
299 | return e; |
---|
300 | } |
---|