[1044] | 1 | // This file is a part of Framsticks SDK. http://www.framsticks.com/ |
---|
[1120] | 2 | // Copyright (C) 1999-2021 Maciej Komosinski and Szymon Ulatowski. |
---|
[1044] | 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[] = { |
---|
[1054] | 14 | { "Creature: Similarity: Descriptor distribution", 1, 4, "SimilMeasureDistribution", "Evaluates morphological dissimilarity using distribution measure.", }, |
---|
[1044] | 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), "", }, |
---|
[1120] | 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 |
---|
[1046] | 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.", }, |
---|
[1044] | 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 | } |
---|
[1120] | 39 | |
---|
[1044] | 40 | for (int i = 0; i < 2; i++) |
---|
| 41 | sst_models[i] = new SolidsShapeTypeModel((*models[i])); |
---|
[1120] | 42 | |
---|
| 43 | SimilMeasureDistribution::calculateFuns(); |
---|
[1044] | 44 | dist = SimilMeasureDistribution::compareFuns(); |
---|
[1120] | 45 | |
---|
[1044] | 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 | } |
---|
[1120] | 62 | |
---|
[1044] | 63 | density = params.at(0); |
---|
| 64 | bin_num = params.at(1); |
---|
| 65 | samples_num = params.at(2); |
---|
[1120] | 66 | |
---|
[1044] | 67 | return 0; |
---|
| 68 | } |
---|
| 69 | |
---|
[1120] | 70 | void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, const Model &sampled) |
---|
[1044] | 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 |
---|
[1120] | 77 | //if (size < (int) sqrt((double) std::numeric_limits<int>::max())) |
---|
| 78 | // samples_taken = min(samples_num, size*size); |
---|
[1044] | 79 | |
---|
[1120] | 80 | 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. |
---|
| 81 | std::uniform_int_distribution<> uniform_distrib(0, sampled.getPartCount() - 1); |
---|
| 82 | |
---|
[1044] | 83 | //Get sampled distribution |
---|
[1120] | 84 | std::vector<double> dist_vect; |
---|
| 85 | dist_vect.reserve(samples_taken); //we will add up to samples_taken elements to this vector |
---|
| 86 | (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, uniform_distrib, sampled, dist_vect); |
---|
[1044] | 87 | |
---|
| 88 | auto result = std::minmax_element(dist_vect.begin(), dist_vect.end()); |
---|
| 89 | double min = *result.first; |
---|
| 90 | double max = *result.second; |
---|
| 91 | |
---|
| 92 | //Create histogram |
---|
[1120] | 93 | vector<int> hist(bin_num); |
---|
[1044] | 94 | int ind = 0; |
---|
| 95 | |
---|
| 96 | for (unsigned int j = 0; j < dist_vect.size(); j++) |
---|
| 97 | { |
---|
[1120] | 98 | ind = (int)std::floor((dist_vect.at(j) - min) * 1 / (max - min) * bin_num); |
---|
| 99 | if (ind <= (bin_num - 1)) |
---|
| 100 | hist[ind]++; |
---|
[1044] | 101 | else if (ind == bin_num) |
---|
[1120] | 102 | hist[bin_num - 1]++; |
---|
[1044] | 103 | } |
---|
| 104 | |
---|
| 105 | //Create pairs |
---|
| 106 | for (int j = 0; j < bin_num; j++) |
---|
| 107 | { |
---|
[1120] | 108 | fun[j] = std::make_pair(min + (max - min) / bin_num * (j + 0.5), hist[j]); |
---|
[1044] | 109 | } |
---|
| 110 | |
---|
| 111 | //Normalize |
---|
| 112 | float total_mass = 0; |
---|
| 113 | for (int j = 0; j < bin_num; j++) |
---|
[1120] | 114 | { |
---|
| 115 | total_mass += fun[j].second; |
---|
| 116 | } |
---|
[1044] | 117 | |
---|
| 118 | for (int j = 0; j < bin_num; j++) |
---|
| 119 | fun[j].second /= total_mass; |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | void SimilMeasureDistribution::calculateFuns() |
---|
| 123 | { |
---|
| 124 | for (int i = 0; i < 2; i++) |
---|
| 125 | { |
---|
| 126 | Model sampled = SimilMeasureDistribution::sampleSurface(&sst_models[i]->getModel(), density); |
---|
| 127 | SimilMeasureDistribution::calculateFun(funs[i], sampled); |
---|
| 128 | } |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | double SimilMeasureDistribution::compareFuns() |
---|
| 132 | { |
---|
| 133 | return SimilMeasureDistribution::EMD(funs[0], funs[1]); |
---|
| 134 | } |
---|
| 135 | |
---|
[1120] | 136 | void SimilMeasureDistribution::D1(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect) |
---|
[1044] | 137 | { |
---|
[1120] | 138 | int size = sampled.getPartCount(); |
---|
[1044] | 139 | double x = 0; |
---|
| 140 | double y = 0; |
---|
| 141 | double z = 0; |
---|
[1120] | 142 | |
---|
[1044] | 143 | for (int i = 0; i < size; i++) |
---|
| 144 | { |
---|
[1120] | 145 | Pt3D pos = sampled.getPart(i)->p; |
---|
[1044] | 146 | x += pos.x; |
---|
| 147 | y += pos.y; |
---|
| 148 | z += pos.z; |
---|
| 149 | } |
---|
[1120] | 150 | |
---|
| 151 | x = x / size; |
---|
| 152 | y = y / size; |
---|
| 153 | z = z / size; |
---|
| 154 | |
---|
| 155 | Pt3D centroid = { x, y, z }; |
---|
| 156 | |
---|
[1044] | 157 | for (int i = 0; i < samples_taken; i++) |
---|
| 158 | { |
---|
[1120] | 159 | int p1 = distribution(rndgen); |
---|
| 160 | double dist = sampled.getPart(p1)->p.distanceTo(centroid); |
---|
[1044] | 161 | if (dist > 0) |
---|
| 162 | dist_vect.push_back(dist); |
---|
| 163 | } |
---|
| 164 | } |
---|
| 165 | |
---|
[1120] | 166 | void SimilMeasureDistribution::D2(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect) |
---|
[1044] | 167 | { |
---|
[1120] | 168 | int size = sampled.getPartCount(); |
---|
[1044] | 169 | for (int i = 0; i < samples_taken; i++) |
---|
| 170 | { |
---|
[1120] | 171 | int p1 = distribution(rndgen); |
---|
| 172 | int p2 = distribution(rndgen); |
---|
| 173 | double dist = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p); |
---|
[1044] | 174 | if (dist > 0) |
---|
| 175 | dist_vect.push_back(dist); |
---|
| 176 | } |
---|
| 177 | } |
---|
| 178 | |
---|
[1120] | 179 | void SimilMeasureDistribution::D3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect) |
---|
[1044] | 180 | { |
---|
| 181 | for (int i = 0; i < samples_taken; i++) |
---|
| 182 | { |
---|
[1120] | 183 | int p1 = distribution(rndgen); |
---|
| 184 | int p2 = distribution(rndgen); |
---|
| 185 | int p3 = distribution(rndgen); |
---|
| 186 | |
---|
| 187 | Pt3D v(sampled.getPart(p2)->p); |
---|
| 188 | Pt3D w(sampled.getPart(p3)->p); |
---|
| 189 | v -= sampled.getPart(p1)->p; |
---|
| 190 | w -= sampled.getPart(p1)->p; |
---|
[1044] | 191 | Pt3D cross_prod(0); |
---|
| 192 | cross_prod.vectorProduct(v, w); |
---|
[1120] | 193 | |
---|
[1044] | 194 | double dist = 0.5 * cross_prod.length(); |
---|
| 195 | if (dist > 0) |
---|
| 196 | dist_vect.push_back(dist); |
---|
| 197 | } |
---|
| 198 | } |
---|
| 199 | |
---|
[1120] | 200 | void SimilMeasureDistribution::D4(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect) |
---|
[1044] | 201 | { |
---|
| 202 | for (int i = 0; i < samples_taken; i++) |
---|
| 203 | { |
---|
[1120] | 204 | int a = distribution(rndgen); |
---|
| 205 | int b = distribution(rndgen); |
---|
| 206 | int c = distribution(rndgen); |
---|
| 207 | int d = distribution(rndgen); |
---|
| 208 | |
---|
| 209 | Pt3D ad(sampled.getPart(a)->p); |
---|
| 210 | Pt3D bd(sampled.getPart(b)->p); |
---|
| 211 | Pt3D cd(sampled.getPart(c)->p); |
---|
| 212 | |
---|
| 213 | ad -= sampled.getPart(d)->p; |
---|
| 214 | bd -= sampled.getPart(d)->p; |
---|
| 215 | cd -= sampled.getPart(d)->p; |
---|
| 216 | |
---|
[1044] | 217 | Pt3D cross_prod(0); |
---|
| 218 | cross_prod.vectorProduct(bd, cd); |
---|
| 219 | cross_prod.entrywiseProduct(ad); |
---|
[1120] | 220 | |
---|
| 221 | double dist = cross_prod.length() / 6; |
---|
[1044] | 222 | if (dist > 0) |
---|
| 223 | dist_vect.push_back(dist); |
---|
| 224 | } |
---|
| 225 | } |
---|
| 226 | |
---|
[1120] | 227 | void SimilMeasureDistribution::A3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect) |
---|
[1044] | 228 | { |
---|
[1120] | 229 | int size = sampled.getPartCount(); |
---|
[1044] | 230 | for (int i = 0; i < samples_taken; i++) |
---|
| 231 | { |
---|
[1120] | 232 | int p1 = distribution(rndgen); |
---|
| 233 | int p2 = distribution(rndgen); |
---|
| 234 | int p3 = distribution(rndgen); |
---|
| 235 | double a = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p3)->p); |
---|
| 236 | double b = sampled.getPart(p3)->p.distanceTo(sampled.getPart(p2)->p); |
---|
| 237 | double c = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p); |
---|
| 238 | double beta = acos((a * a + b * b - c * c) / (2 * a * b)); |
---|
| 239 | |
---|
[1044] | 240 | if (!std::isnan(beta)) |
---|
| 241 | dist_vect.push_back(beta); |
---|
| 242 | } |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | float dist(feature_t* F1, feature_t* F2) |
---|
[1120] | 247 | { |
---|
| 248 | return abs((*F1) - (*F2)); |
---|
[1044] | 249 | } |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | void SimilMeasureDistribution::fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights) |
---|
| 253 | { |
---|
| 254 | for (int j = 0; j < bin_num; j++) |
---|
| 255 | { |
---|
[1120] | 256 | points[j] = { fun[j].first }; |
---|
[1044] | 257 | weights[j] = fun[j].second; |
---|
| 258 | } |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | double SimilMeasureDistribution::EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2) |
---|
| 262 | { |
---|
| 263 | feature_t *points[2]; |
---|
| 264 | float *weights[2]; |
---|
[1120] | 265 | |
---|
[1044] | 266 | for (int i = 0; i < 2; i++) |
---|
| 267 | { |
---|
| 268 | points[i] = new feature_t[bin_num]; |
---|
| 269 | weights[i] = new float[bin_num](); |
---|
| 270 | } |
---|
| 271 | SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]); |
---|
| 272 | SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]); |
---|
| 273 | |
---|
[1120] | 274 | signature_t sig1 = { bin_num, points[0], weights[0] }, |
---|
| 275 | sig2 = { bin_num, points[1], weights[1] }; |
---|
| 276 | |
---|
[1062] | 277 | float e = emd(&sig1, &sig2, dist, 0, 0, bin_num, bin_num); |
---|
[1120] | 278 | |
---|
[1044] | 279 | for (int i = 0; i < 2; i++) |
---|
| 280 | { |
---|
| 281 | delete[] points[i]; |
---|
| 282 | delete[] weights[i]; |
---|
| 283 | } |
---|
| 284 | |
---|
| 285 | return e; |
---|
| 286 | } |
---|