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 FIELDSTRUCT SimilMeasureDistribution |
---|
12 | |
---|
13 | static ParamEntry simil_distribution_paramtab[] = { |
---|
14 | { "Similarity: 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 1048576", FIELD(samples_num), "", }, |
---|
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, Model &sampled) |
---|
71 | { |
---|
72 | int size = sampled.getPartCount(); |
---|
73 | int samples_taken = samples_num; |
---|
74 | |
---|
75 | //Check if total number of points pairs is smaller than samples number |
---|
76 | //but first, prevent exceeding int limits |
---|
77 | if (size < (int) sqrt((double) std::numeric_limits<int>::max())) |
---|
78 | samples_taken = min(samples_num, size*size); |
---|
79 | |
---|
80 | //Get sampled distribution |
---|
81 | std::vector<double> dist_vect = (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, &sampled); |
---|
82 | |
---|
83 | auto result = std::minmax_element(dist_vect.begin(), dist_vect.end()); |
---|
84 | double min = *result.first; |
---|
85 | double max = *result.second; |
---|
86 | |
---|
87 | //Create histogram |
---|
88 | vector<float> hist(bin_num); |
---|
89 | int ind = 0; |
---|
90 | |
---|
91 | for (unsigned int j = 0; j < dist_vect.size(); j++) |
---|
92 | { |
---|
93 | ind = (int) std::floor((dist_vect.at(j)-min)*1/(max-min)*bin_num); |
---|
94 | if (ind <= (bin_num-1)) |
---|
95 | hist[ind] += 1; |
---|
96 | else if (ind == bin_num) |
---|
97 | hist[bin_num-1] += 1; |
---|
98 | } |
---|
99 | |
---|
100 | //Create pairs |
---|
101 | for (int j = 0; j < bin_num; j++) |
---|
102 | { |
---|
103 | fun[j] = std::make_pair(min+(max-min)/bin_num*(j+0.5), hist[j]); |
---|
104 | } |
---|
105 | |
---|
106 | //Normalize |
---|
107 | float total_mass = 0; |
---|
108 | for (int j = 0; j < bin_num; j++) |
---|
109 | { |
---|
110 | total_mass += fun[j].second; |
---|
111 | } |
---|
112 | |
---|
113 | for (int j = 0; j < bin_num; j++) |
---|
114 | fun[j].second /= total_mass; |
---|
115 | } |
---|
116 | |
---|
117 | void SimilMeasureDistribution::calculateFuns() |
---|
118 | { |
---|
119 | for (int i = 0; i < 2; i++) |
---|
120 | { |
---|
121 | Model sampled = SimilMeasureDistribution::sampleSurface(&sst_models[i]->getModel(), density); |
---|
122 | SimilMeasureDistribution::calculateFun(funs[i], sampled); |
---|
123 | } |
---|
124 | } |
---|
125 | |
---|
126 | double SimilMeasureDistribution::compareFuns() |
---|
127 | { |
---|
128 | return SimilMeasureDistribution::EMD(funs[0], funs[1]); |
---|
129 | } |
---|
130 | |
---|
131 | vector<double> SimilMeasureDistribution::D1(int samples_taken, Model *sampled) |
---|
132 | { |
---|
133 | vector<double> dist_vect; |
---|
134 | int size = sampled->getPartCount(); |
---|
135 | double x = 0; |
---|
136 | double y = 0; |
---|
137 | double z = 0; |
---|
138 | |
---|
139 | for (int i = 0; i < size; i++) |
---|
140 | { |
---|
141 | Pt3D pos = sampled->getPart(i)->p; |
---|
142 | x += pos.x; |
---|
143 | y += pos.y; |
---|
144 | z += pos.z; |
---|
145 | } |
---|
146 | |
---|
147 | x = x/size; |
---|
148 | y = y/size; |
---|
149 | z = z/size; |
---|
150 | |
---|
151 | Pt3D centroid = {x, y, z}; |
---|
152 | |
---|
153 | for (int i = 0; i < samples_taken; i++) |
---|
154 | { |
---|
155 | int p1 = rndUint(size); |
---|
156 | double dist = sampled->getPart(p1)->p.distanceTo(centroid); |
---|
157 | if (dist > 0) |
---|
158 | dist_vect.push_back(dist); |
---|
159 | } |
---|
160 | |
---|
161 | return dist_vect; |
---|
162 | } |
---|
163 | |
---|
164 | vector<double> SimilMeasureDistribution::D2(int samples_taken, Model *sampled) |
---|
165 | { |
---|
166 | vector<double> dist_vect; |
---|
167 | int size = sampled->getPartCount(); |
---|
168 | for (int i = 0; i < samples_taken; i++) |
---|
169 | { |
---|
170 | int p1 = rndUint(size); |
---|
171 | int p2 = rndUint(size); |
---|
172 | double dist = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p); |
---|
173 | if (dist > 0) |
---|
174 | dist_vect.push_back(dist); |
---|
175 | } |
---|
176 | |
---|
177 | return dist_vect; |
---|
178 | } |
---|
179 | |
---|
180 | vector<double> SimilMeasureDistribution::D3(int samples_taken, Model *sampled) |
---|
181 | { |
---|
182 | vector<double> dist_vect; |
---|
183 | int size = sampled->getPartCount(); |
---|
184 | for (int i = 0; i < samples_taken; i++) |
---|
185 | { |
---|
186 | int p1 = rndUint(size); |
---|
187 | int p2 = rndUint(size); |
---|
188 | int p3 = rndUint(size); |
---|
189 | |
---|
190 | Pt3D v(sampled->getPart(p2)->p); |
---|
191 | Pt3D w(sampled->getPart(p3)->p); |
---|
192 | v -= sampled->getPart(p1)->p; |
---|
193 | w -= sampled->getPart(p1)->p; |
---|
194 | Pt3D cross_prod(0); |
---|
195 | cross_prod.vectorProduct(v, w); |
---|
196 | |
---|
197 | double dist = 0.5 * cross_prod.length(); |
---|
198 | if (dist > 0) |
---|
199 | dist_vect.push_back(dist); |
---|
200 | } |
---|
201 | |
---|
202 | return dist_vect; |
---|
203 | } |
---|
204 | |
---|
205 | vector<double> SimilMeasureDistribution::D4(int samples_taken, Model *sampled) |
---|
206 | { |
---|
207 | vector<double> dist_vect; |
---|
208 | int size = sampled->getPartCount(); |
---|
209 | for (int i = 0; i < samples_taken; i++) |
---|
210 | { |
---|
211 | int a = rndUint(size); |
---|
212 | int b = rndUint(size); |
---|
213 | int c = rndUint(size); |
---|
214 | int d = rndUint(size); |
---|
215 | |
---|
216 | Pt3D ad(sampled->getPart(a)->p); |
---|
217 | Pt3D bd(sampled->getPart(b)->p); |
---|
218 | Pt3D cd(sampled->getPart(c)->p); |
---|
219 | |
---|
220 | ad -= sampled->getPart(d)->p; |
---|
221 | bd -= sampled->getPart(d)->p; |
---|
222 | cd -= sampled->getPart(d)->p; |
---|
223 | |
---|
224 | Pt3D cross_prod(0); |
---|
225 | cross_prod.vectorProduct(bd, cd); |
---|
226 | cross_prod.entrywiseProduct(ad); |
---|
227 | |
---|
228 | double dist = cross_prod.length()/6; |
---|
229 | if (dist > 0) |
---|
230 | dist_vect.push_back(dist); |
---|
231 | } |
---|
232 | |
---|
233 | return dist_vect; |
---|
234 | } |
---|
235 | |
---|
236 | vector<double> SimilMeasureDistribution::A3(int samples_taken, Model *sampled) |
---|
237 | { |
---|
238 | vector<double> dist_vect; |
---|
239 | int size = sampled->getPartCount(); |
---|
240 | for (int i = 0; i < samples_taken; i++) |
---|
241 | { |
---|
242 | int p1 = rndUint(size); |
---|
243 | int p2 = rndUint(size); |
---|
244 | int p3 = rndUint(size); |
---|
245 | double a = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p3)->p); |
---|
246 | double b = sampled->getPart(p3)->p.distanceTo(sampled->getPart(p2)->p); |
---|
247 | double c = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p); |
---|
248 | double beta = acos((a*a + b*b - c*c)/(2*a*b)); |
---|
249 | |
---|
250 | if (!std::isnan(beta)) |
---|
251 | dist_vect.push_back(beta); |
---|
252 | } |
---|
253 | |
---|
254 | return dist_vect; |
---|
255 | } |
---|
256 | |
---|
257 | |
---|
258 | float dist(feature_t* F1, feature_t* F2) |
---|
259 | { |
---|
260 | return abs((*F1)-(*F2)); |
---|
261 | } |
---|
262 | |
---|
263 | |
---|
264 | void SimilMeasureDistribution::fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights) |
---|
265 | { |
---|
266 | for (int j = 0; j < bin_num; j++) |
---|
267 | { |
---|
268 | points[j] = {fun[j].first}; |
---|
269 | weights[j] = fun[j].second; |
---|
270 | } |
---|
271 | } |
---|
272 | |
---|
273 | double SimilMeasureDistribution::EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2) |
---|
274 | { |
---|
275 | feature_t *points[2]; |
---|
276 | float *weights[2]; |
---|
277 | |
---|
278 | for (int i = 0; i < 2; i++) |
---|
279 | { |
---|
280 | points[i] = new feature_t[bin_num]; |
---|
281 | weights[i] = new float[bin_num](); |
---|
282 | } |
---|
283 | SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]); |
---|
284 | SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]); |
---|
285 | |
---|
286 | signature_t sig1 = {bin_num, points[0], weights[0]}, |
---|
287 | sig2 = {bin_num, points[1], weights[1]}; |
---|
288 | |
---|
289 | float e = emd(&sig1, &sig2, dist, 0, 0); |
---|
290 | |
---|
291 | for (int i = 0; i < 2; i++) |
---|
292 | { |
---|
293 | delete[] points[i]; |
---|
294 | delete[] weights[i]; |
---|
295 | } |
---|
296 | |
---|
297 | return e; |
---|
298 | } |
---|