Ignore:
Timestamp:
05/02/19 23:50:27 (5 years ago)
Author:
Maciej Komosinski
Message:

Added another, improved way of calculating dissimilarity of two creatures/models. Details and comparisons in ​https://doi.org/10.1007/978-3-030-16692-2_8

File:
1 edited

Legend:

Unmodified
Added
Removed

  • cpp/frams/model/similarity/simil_model.cpp

    r818 r869  
    11// This file is a part of Framsticks SDK.  http://www.framsticks.com/
    2 // Copyright (C) 1999-2017  Maciej Komosinski and Szymon Ulatowski.
     2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
    33// See LICENSE.txt for details.
    44
    5 // simil_model.cpp: implementation of the ModelSimil class.
    6 //
    7 //////////////////////////////////////////////////////////////////////
     5// implementation of the ModelSimil class.
     6
    87#include "SVD/matrix_tools.h"
     8#include "hungarian/hungarian.h"
    99#include "simil_model.h"
    1010#include "simil_match.h"
    … …  
    3535
    3636static ParamEntry MSparam_tab[] = {
    37                 { "Creature: Similarity", 1, 7, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016", },
     37                { "Creature: Similarity", 1, 8, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016\nhttps://doi.org/10.1007/978-3-030-16692-2_8", },
     38                { "simil_method", 0, 0, "Similarity algorithm", "d 0 1 0 ~New (flexible criteria order, optimal matching)~Old (vertex degree order, greedy matching)", FIELD(matching_method), "",},
    3839                { "simil_parts", 0, 0, "Weight of parts count", "f 0 100 0", FIELD(m_adFactors[0]), "Differing number of parts is also handled by the 'part degree' similarity component.", },
    3940                { "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
    … …  
    7677        isFuzzy = 0;
    7778        fuzzyDepth = 10;
    78        
     79
    7980        //Determines whether weighted MDS should be used.
    8081        wMDS = 0;
    81 }
    82 
    83 /**     Evaluates distance between two given genotypes. The distance depends strongly
    84         on weights set.
    85         @param G0 Pointer to the first of compared genotypes
    86         @param G1 Pointer to the second of compared genotypes.
    87         @return Distance between two genotypes.
    88         @sa m_adFactors, matching_method
    89         */
    90 double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
     82        //Determines whether best matching should be saved using hungarian similarity measure.
     83        saveMatching = 0;
     84}
     85
     86double ModelSimil::EvaluateDistanceGreedy(const Geno *G0, const Geno *G1)
    9187{
    9288        double dResult = 0;
    … …  
    9894        if (m_Gen[0] == NULL || m_Gen[1] == NULL)
    9995        {
    100                 DB(printf("ModelSimil::EvaluateDistance - invalid genotypes pointers\n");)
    101                         return 0.0;
     96                DB(printf("ModelSimil::EvaluateDistanceGreedy - invalid genotype(s) pointers\n");) //TODO convert all such error printfs to legacy error messages, since if's are not in DB(). Example below.
     97                        logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "NULL genotype pointer(s)");
     98                return 0.0;
    10299        }
    103100        // create models of objects to compare
    … …  
    108105        if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
    109106        {
    110                 DB(printf("ModelSimil::EvaluateDistance - invalid models pointers\n");)
     107                DB(printf("ModelSimil::EvaluateDistanceGreedy - invalid model(s) pointers\n");)
    111108                        return 0.0;
    112109        }
    … …  
    143140        if ((this->*pfMatchingFunction)() == 0)
    144141        {
    145                 DB(printf("ModelSimil::EvaluateDistance - MatchParts() error\n");)
     142                DB(printf("ModelSimil::EvaluateDistanceGreedy - MatchParts() error\n");)
    146143                        return 0.0;
    147144        }
    … …  
    155152                if (CountPartsDistance() == 0)
    156153                {
    157                 DB(printf("ModelSimil::EvaluateDistance - CountPartDistance() error\n");)
    158                         return 0.0;
     154                        DB(printf("ModelSimil::EvaluateDistanceGreedy - CountPartDistance() error\n");)
     155                                return 0.0;
    159156                }
    160157
    … …  
    340337                if (tdn1[1] < tdn2[1])
    341338                {
    342                 // when degree - tdn2[1] - is BIGGER
    343                 return 1;
     339                        // when degree - tdn2[1] - is BIGGER
     340                        return 1;
     341                }
     342                else
     343                {
     344                        return 0;
     345                }
     346}
     347
     348/** Comparison function required for qsort() call. Used while sorting groups of
     349        Parts with respect to fuzzy vertex degree. Compares two TDN structures
     350        with respect to their [4] field ( fuzzy vertex degree). Highest degree goes first.
     351        @param pElem1 Pointer to the TDN structure of some Part.
     352        @param pElem2 Pointer to the TDN structure of some Part.
     353        @return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
     354        */
     355int ModelSimil::CompareFuzzyDegrees(const void *pElem1, const void *pElem2)
     356{
     357        int *tdn1 = (int *)pElem1;
     358        int *tdn2 = (int *)pElem2;
     359
     360        if (tdn1[4] > tdn2[4])
     361        {
     362                // when degree - tdn1[4] - is BIGGER
     363                return -1;
     364        }
     365        else
     366                if (tdn1[4] < tdn2[4])
     367                {
     368                        // when degree - tdn2[4] - is BIGGER
     369                        return 1;
    344370                }
    345371                else
    … …  
    374400                if (tdn1[NEURO_CONNS] < tdn2[NEURO_CONNS])
    375401                {
    376                 // when number of NConn Elem1 is SMALLER
    377                 return 1;
     402                        // when number of NConn Elem1 is SMALLER
     403                        return 1;
    378404                }
    379405                else
    … …  
    389415                                if (tdn1[NEURONS] < tdn2[NEURONS])
    390416                                {
    391                                 // when number of Neu is SMALLER for Elem1
    392                                 return 1;
     417                                        // when number of Neu is SMALLER for Elem1
     418                                        return 1;
    393419                                }
    394420                                else
    … …  
    406432                                                if (tdn1[ORIG_IND] < tdn2[ORIG_IND])
    407433                                                {
    408                                                 // when the number of NIt Deg1 id SMALLER
    409                                                 return 1;
     434                                                        // when the number of NIt Deg1 id SMALLER
     435                                                        return 1;
    410436                                                }
    411437                                                else
    … …  
    494520}
    495521
     522/** Prints one array of parts fuzzy neighbourhood.
     523        @param index of organism
     524        */
     525void ModelSimil::_PrintFuzzyNeighbourhood(int o)
     526{
     527        assert(m_fuzzyNeighb[o] != NULL);
     528        printf("Fuzzy neighbourhhod of organism %i\n", o);
     529        int size = m_Mod[o]->getPartCount();
     530        for (int i = 0; i < size; i++)
     531        {
     532                for (int j = 0; j < fuzzyDepth; j++)
     533                {
     534                        printf("%f ", m_fuzzyNeighb[o][i][j]);
     535                }
     536                printf("\n");
     537        }
     538}
     539
    496540/** Creates arrays holding information about organisms' Parts (m_aDegrees) andm_Neigh
    497541        fills them with initial data (original indices and zeros).
    … …  
    526570                                for (j = 0; j < partCount; j++)
    527571                                {
    528                                 m_aDegrees[i][j][0] = j;
    529                                 m_aDegrees[i][j][1] = 0;
    530                                 m_aDegrees[i][j][2] = 0;
    531                                 m_aDegrees[i][j][3] = 0;
    532                                 m_aDegrees[i][j][4] = 0;
    533 
    534                                 // sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
    535                                 assert(m_aDegrees[i][j] != NULL);
    536 
    537                                 if (isFuzzy)
    538                                 {
    539                                         m_Neighbours[i][j] = new int[partCount];
    540                                         for (int k = 0; k < partCount; k++)
     572                                        m_aDegrees[i][j][0] = j;
     573                                        m_aDegrees[i][j][1] = 0;
     574                                        m_aDegrees[i][j][2] = 0;
     575                                        m_aDegrees[i][j][3] = 0;
     576                                        m_aDegrees[i][j][4] = 0;
     577
     578                                        // sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
     579                                        assert(m_aDegrees[i][j] != NULL);
     580
     581                                        if (isFuzzy)
    541582                                        {
    542                                                 m_Neighbours[i][j][k] = 0;
     583                                                m_Neighbours[i][j] = new int[partCount];
     584                                                for (int k = 0; k < partCount; k++)
     585                                                {
     586                                                        m_Neighbours[i][j][k] = 0;
     587                                                }
     588
     589                                                m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
     590                                                for (int k = 0; k < fuzzyDepth; k++)
     591                                                {
     592                                                        m_fuzzyNeighb[i][j][k] = 0;
     593                                                }
     594
     595                                                assert(m_Neighbours[i][j] != NULL);
     596                                                assert(m_fuzzyNeighb[i][j] != NULL);
    543597                                        }
    544 
    545                                         m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
    546                                         for (int k = 0; k < fuzzyDepth; k++)
    547                                         {
    548                                                 m_fuzzyNeighb[i][j][k] = 0;
    549                                         }
    550 
    551                                         assert(m_Neighbours[i][j] != NULL);
    552                                         assert(m_fuzzyNeighb[i][j] != NULL);
    553                                 }
    554598
    555599                                }
    … …  
    671715}
    672716
    673 //store information about identity of parts "fuzzy degrees" in the m_aDegrees[4]
    674 void ModelSimil::CheckFuzzyIdentity()
    675 {
    676         int partCount = 0;
     717void ModelSimil::FuzzyOrder()
     718{
     719        int i, depth, partInd, prevPartInd, partCount;
    677720        for (int mod = 0; mod < 2; mod++)
    678721        {
    679                 //for subsequent pairs of parts
    680722                partCount = m_Mod[mod]->getPartCount();
    681                 m_aDegrees[mod][partCount - 1][FUZZ_DEG] = 0;
    682                 for (int partInd = (partCount - 2); partInd >= 0; partInd--)
    683                 {
    684                         m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][partInd + 1][FUZZ_DEG];
    685                         for (int depth = 1; depth < fuzzyDepth; depth++)
    686                         {
    687                                 if (m_fuzzyNeighb[mod][partInd][depth] != m_fuzzyNeighb[mod][partInd + 1][depth])
    688                                 {
    689                                         m_aDegrees[mod][partInd][FUZZ_DEG] += 1;
     723                partInd = m_fuzzyNeighb[mod][partCount - 1][0];
     724                m_aDegrees[mod][partInd][FUZZ_DEG] = 0;
     725
     726                for (i = (partCount - 2); i >= 0; i--)
     727                {
     728                        prevPartInd = partInd;
     729                        partInd = m_fuzzyNeighb[mod][i][0];
     730                        m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][prevPartInd][FUZZ_DEG];
     731                        for (depth = 1; depth < fuzzyDepth; depth++)
     732                        {
     733                                if (m_fuzzyNeighb[mod][i][depth] != m_fuzzyNeighb[mod][i + 1][depth])
     734                                {
     735                                        m_aDegrees[mod][partInd][FUZZ_DEG]++;
    690736                                        break;
    691737                                }
    … …  
    761807
    762808        SortFuzzyNeighb();
     809        FuzzyOrder();
    763810}
    764811
    … …  
    884931
    885932        int i;
     933        int(*pfDegreeFunction) (const void*, const void*) = NULL;
     934        pfDegreeFunction = (isFuzzy == 1) ? &CompareFuzzyDegrees : &CompareDegrees;
    886935        // sortowanie obu tablic wg stopni punktów - TDN[1]
    887         if (isFuzzy != 1)
    888         {
    889                 for (i = 0; i < 2; i++)
    890                 {
    891                         DB(_PrintDegrees(i));
    892                         std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
    893                                 sizeof(TDN), ModelSimil::CompareDegrees);
    894                         DB(_PrintDegrees(i));
    895                 }
    896         }//sortowanie wg romzmytego stopnia wierzcholka
    897 
    898         else
    899         {
    900                 SortPartInfoFuzzy();
    901         }
    902 
     936        for (i = 0; i < 2; i++)
     937        {
     938                DB(_PrintDegrees(i));
     939                std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
     940                        sizeof(TDN), pfDegreeFunction);
     941                DB(_PrintDegrees(i));
     942        }
    903943
    904944        // sprawdzenie wartosci parametru
    … …  
    935975
    936976                                        std::qsort(m_aDegrees[i][iPocz], (size_t)(j - iPocz),
    937                                         sizeof(TDN), ModelSimil::CompareConnsNo);
     977                                                sizeof(TDN), ModelSimil::CompareConnsNo);
    938978                                DB(_PrintDegrees(i));
    939979                                // wyswietlamy z jedna komorka po zakonczeniu tablicy
    … …  
    948988}
    949989
    950 int ModelSimil::SortPartInfoFuzzy()
    951 {
    952 
    953         // sprawdz zalozenie - o modelach
    954         assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
    955         assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
    956 
    957         // sprawdz zalozenie - o tablicach
    958         assert(m_aDegrees[0] != NULL);
    959         assert(m_aDegrees[1] != NULL);
    960         // sprawdz zalozenie - o tablicach
    961         assert(m_fuzzyNeighb[0] != NULL);
    962         assert(m_fuzzyNeighb[1] != NULL);
    963 
    964 
    965         TDN * m_aDegreesTmp[2];
    966 
    967         for (int i = 0; i < 2; i++)
    968         {
    969                 int partCount = m_Mod[i]->getPartCount();
    970                 m_aDegreesTmp[i] = new TDN[partCount];
    971 
    972                 for (int j = 0; j < partCount; j++)
    973                 {
    974                         for (int k = 0; k < TDN_SIZE; k++)
    975                         {
    976                                 m_aDegreesTmp[i][j][k] = m_aDegrees[i][j][k];
    977                         }
    978                 }
    979         }
    980 
    981         int newInd = 0;
    982         for (int i = 0; i < 2; i++)
    983         {
    984                 int partCount = m_Mod[i]->getPartCount();
    985                 for (int j = 0; j < partCount; j++)
    986                 {
    987                         newInd = (int)m_fuzzyNeighb[i][j][0];
    988                         for (int k = 0; k < TDN_SIZE; k++)
    989                         {
    990                                 m_aDegrees[i][j][k] = m_aDegreesTmp[i][newInd][k];
    991                         }
    992                 }
    993         }
    994 
    995         SAFEDELETEARRAY(m_aDegreesTmp[0]);
    996         SAFEDELETEARRAY(m_aDegreesTmp[1]);
    997 
    998         CheckFuzzyIdentity();
    999 
    1000         return 1;
    1001 }
    1002 
    1003 /** Checks if given Parts have identical physical and biological properties
    1004         (except for geometry that might differ).
    1005         @param P1 Pointer to first Part.
    1006         @param P2 Pointer to second Part.
    1007         @return 1 - identical properties, 0 - there are differences.
    1008         */
    1009 int ModelSimil::CheckPartsIdentity(Part *P1, Part *P2)
    1010 {
    1011         // sprawdz, czy te Parts chociaz sa w sensownym miejscu pamieci
    1012         assert((P1 != NULL) && (P2 != NULL));
    1013 
    1014         if ((P1->assim != P2->assim) ||
    1015                 (P1->friction != P2->friction) ||
    1016                 (P1->ingest != P2->ingest) ||
    1017                 (P1->mass != P2->mass) ||
    1018                 (P1->size != P2->size) ||
    1019                 (P1->density != P2->density))
    1020                 // gdy znaleziono jakas roznice w parametrach fizycznych i
    1021                 // biologicznych
    1022                 return 0;
    1023         else
    1024                 // gdy nie ma roznic
    1025                 return 1;
    1026 }
    1027990
    1028991/** Prints the state of the matching object. Debug method.
    … …  
    11581121                        aiKoniecGrupyDopasowania[0]);)
    11591122                        DB(printf("Organizm 1: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[1],
    1160                         aiKoniecGrupyDopasowania[1]);)
     1123                                aiKoniecGrupyDopasowania[1]);)
    11611124
    11621125                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
    11631126                        {
    11641127
    1165                         // iterujemy i - Parts organizmu 0
    1166                         // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
    1167 
    1168                         if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
    1169                         {
    1170                                 // interesuja nas tylko te niedopasowane jeszcze (i)
    1171                                 for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
    1172                                 {
    1173 
    1174                                         // iterujemy j - Parts organizmu 1
    1175                                         // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
    1176 
    1177                                         if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
     1128                                // iterujemy i - Parts organizmu 0
     1129                                // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
     1130
     1131                                if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
     1132                                {
     1133                                        // interesuja nas tylko te niedopasowane jeszcze (i)
     1134                                        for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
    11781135                                        {
    1179                                                 // interesuja nas tylko te niedopasowane jeszcze (j)
    1180                                                 // teraz obliczymy lokalne różnice pomiędzy Parts
    1181                                                 iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
    1182                                                         - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
    1183 
    1184                                                 //iNit currently is not a component of distance measure           
    1185                                                 //iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
    1186                                                 //           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
    1187 
    1188                                                 iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
    1189                                                         - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
    1190 
    1191                                                 // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
    1192                                                 // find original indices of Parts for both of the models
    1193                                                 iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
    1194                                                 iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
    1195                                                 // now compute the geometrical distance of these Parts (use m_aPositions
    1196                                                 // which should be computed by SVD)
    1197                                                 Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
    1198                                                 Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
    1199                                                 dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
    1200 
    1201                                                 // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
    1202                                                 // identyczności pozostałych własności (oprócz geometrii)
    1203                                                 // - żeby móc stwierdzić w ogóle identyczność Parts
    1204 
    1205                                                 // i ostateczna odleglosc indukowana przez te roznice
    1206                                                 // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
    1207                                                 aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
    1208                                                         m_adFactors[2] * double(iNeu) +
    1209                                                         m_adFactors[3] * dGeo;
    1210                                                 DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
    1211                                                         aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
    1212                                                         DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
    1213                                                         DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
    1214                                                         DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
     1136
     1137                                                // iterujemy j - Parts organizmu 1
     1138                                                // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
     1139
     1140                                                if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
     1141                                                {
     1142                                                        // interesuja nas tylko te niedopasowane jeszcze (j)
     1143                                                        // teraz obliczymy lokalne różnice pomiędzy Parts
     1144                                                        iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
     1145                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
     1146
     1147                                                        //iNit currently is not a component of distance measure           
     1148                                                        //iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
     1149                                                        //           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
     1150
     1151                                                        iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
     1152                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
     1153
     1154                                                        // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
     1155                                                        // find original indices of Parts for both of the models
     1156                                                        iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
     1157                                                        iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
     1158                                                        // now compute the geometrical distance of these Parts (use m_aPositions
     1159                                                        // which should be computed by SVD)
     1160                                                        Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
     1161                                                        Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
     1162                                                        dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
     1163
     1164                                                        // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
     1165                                                        // identyczności pozostałych własności (oprócz geometrii)
     1166                                                        // - żeby móc stwierdzić w ogóle identyczność Parts
     1167
     1168                                                        // i ostateczna odleglosc indukowana przez te roznice
     1169                                                        // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
     1170                                                        aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
     1171                                                                m_adFactors[2] * double(iNeu) +
     1172                                                                m_adFactors[3] * dGeo;
     1173                                                        DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
     1174                                                                aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
     1175                                                                DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
     1176                                                                DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
     1177                                                                DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
     1178                                                }
    12151179                                        }
    12161180                                }
    1217                         }
    12181181                        }
    12191182
    … …  
    13881351                                        if (PartZ1NaLiscie0 || PartZ0NaLiscie1)
    13891352                                        {
    1390                                         // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
    1391                                         // mozliwych dopasowan
    1392                                         // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
    1393                                         // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
    1394                                         // duza czesc obliczen (znalezc mu nowa mozliwa pare)
    1395 
    1396                                         // indeks organizmu, ktorego Part nie ma dopasowywanego Part
    1397                                         // z drugiego organizmu na swojej liscie
    1398                                         int iBezDrugiego;
    1399 
    1400                                         // okreslenie indeksu organizmu z dopasowywanym Part
    1401                                         if (!PartZ1NaLiscie0)
    1402                                         {
    1403                                                 iBezDrugiego = 0;
    1404                                         }
    1405                                         else
    1406                                         {
    1407                                                 iBezDrugiego = 1;
    1408                                         }
    1409                                         // sprawdz indeks organizmu
    1410                                         assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
    1411 
    1412                                         int iDopasowywany = -1;
    1413                                         // poszukujemy pierwszego z listy dopasowania
    1414                                         for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
    1415                                         {
    1416                                                 if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
     1353                                                // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
     1354                                                // mozliwych dopasowan
     1355                                                // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
     1356                                                // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
     1357                                                // duza czesc obliczen (znalezc mu nowa mozliwa pare)
     1358
     1359                                                // indeks organizmu, ktorego Part nie ma dopasowywanego Part
     1360                                                // z drugiego organizmu na swojej liscie
     1361                                                int iBezDrugiego;
     1362
     1363                                                // okreslenie indeksu organizmu z dopasowywanym Part
     1364                                                if (!PartZ1NaLiscie0)
    14171365                                                {
    1418                                                         iDopasowywany = i;
    1419                                                         break;
     1366                                                        iBezDrugiego = 0;
    14201367                                                }
    1421                                         }
    1422                                         // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
    1423                                         // nieujemny i w odpowiedniej grupie!
    1424                                         assert((iDopasowywany >= 0) &&
    1425                                                 (iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
    1426 
    1427                                         // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
    1428                                         m_pMatching->Match(
    1429                                                 iBezDrugiego,
    1430                                                 aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
    1431                                                 1 - iBezDrugiego,
    1432                                                 aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
    1433                                         DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
    1434                                                 aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
    1435 
    1436                                                 // zmniejsz liczby niedopasowanych elementow w grupach
    1437                                                 aiRozmiarGrupy[0]--;
    1438                                         aiRozmiarGrupy[1]--;
    1439 
    1440                                         // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
    1441                                         // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
    1442                                         if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
    1443                                         {
    1444                                                 // jesli grupy sie jeszcze nie wyczrpaly
    1445                                                 // to jest mozliwosc dopasowania w organizmie
    1446 
    1447                                                 int iZDrugim = 1 - iBezDrugiego;
     1368                                                else
     1369                                                {
     1370                                                        iBezDrugiego = 1;
     1371                                                }
    14481372                                                // sprawdz indeks organizmu
    1449                                                 assert((iZDrugim == 0) || (iZDrugim == 1));
    1450 
    1451                                                 // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
    1452                                                 // poprzednie obliczenia minimum
    1453                                                 // dla organizmu 1 (o rozmiarze grupy z 0)
    1454                                                 for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
     1373                                                assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
     1374
     1375                                                int iDopasowywany = -1;
     1376                                                // poszukujemy pierwszego z listy dopasowania
     1377                                                for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
    14551378                                                {
    1456                                                         apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
    1457                                                 }
    1458                                                 // dla organizmu 0 (o rozmiarze grup z 1)
    1459                                                 for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
    1460                                                 {
    1461                                                         apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
    1462                                                 }
    1463 
    1464                                                 // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
    1465                                                 // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
    1466                                                 dMinimum = HUGE_VAL;
    1467                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1468                                                 {
    1469                                                         if (!(m_pMatching->IsMatched(
    1470                                                                 1 - iZDrugim,
    1471                                                                 aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
    1472                                                         {
    1473                                                                 // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
    1474                                                                 // niedopasowanych jeszcze Parts
    1475                                                                 if (iZDrugim == 0)
    1476                                                                 {
    1477                                                                         // teraz niestety musimy rozpoznac odpowiedni organizm
    1478                                                                         // zeby moc indeksowac niesymetryczna tablice
    1479                                                                         if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
    1480                                                                         {
    1481                                                                                 dMinimum = aadOdleglosciParts[iIndex[0]][i];
    1482                                                                         }
    1483                                                                         // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
    1484                                                                         assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
    1485 
    1486                                                                 }
    1487                                                                 else
    1488                                                                 {
    1489                                                                         if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
    1490                                                                         {
    1491                                                                                 dMinimum = aadOdleglosciParts[i][iIndex[1]];
    1492                                                                         }
    1493                                                                         // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
    1494                                                                         assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
    1495                                                                 }
    1496                                                         }
    1497                                                 }
    1498                                                 // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
    1499                                                 assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
    1500 
    1501                                                 // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
    1502                                                 // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
    1503                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1504                                                 {
    1505                                                         if (!(m_pMatching->IsMatched(
    1506                                                                 1 - iZDrugim,
    1507                                                                 aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
    1508                                                         {
    1509                                                                 if (iZDrugim == 0)
    1510                                                                 {
    1511                                                                         // teraz niestety musimy rozpoznac odpowiedni organizm
    1512                                                                         // zeby moc indeksowac niesymetryczna tablice
    1513                                                                         if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
    1514                                                                         {
    1515                                                                                 // jesli w danym miejscu tablicy odleglosci jest faktyczne
    1516                                                                                 // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
    1517                                                                                 apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
    1518                                                                         }
    1519                                                                 }
    1520                                                                 else
    1521                                                                 {
    1522                                                                         if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
    1523                                                                         {
    1524                                                                                 apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
    1525                                                                         }
    1526                                                                 }
    1527                                                         }
    1528                                                 }
    1529 
    1530                                                 // a teraz - po znalezieniu potencjalnych elementow do dopasowania
    1531                                                 // dopasowujemy pierwszy z potencjalnych
    1532                                                 iDopasowywany = -1;
    1533                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1534                                                 {
    1535                                                         if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
     1379                                                        if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
    15361380                                                        {
    15371381                                                                iDopasowywany = i;
    … …  
    15391383                                                        }
    15401384                                                }
    1541                                                 // musielismy znalezc jakiegos dopasowywanego
     1385                                                // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
     1386                                                // nieujemny i w odpowiedniej grupie!
    15421387                                                assert((iDopasowywany >= 0) &&
    1543                                                         (iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
    1544 
    1545                                                 // no to juz mozemy dopasowac
     1388                                                        (iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
     1389
     1390                                                // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
    15461391                                                m_pMatching->Match(
    1547                                                         iZDrugim,
    1548                                                         aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
    1549                                                         1 - iZDrugim,
    1550                                                         aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
    1551                                                 DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
    1552 
    1553                                                         //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
    1554                                                         //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1392                                                        iBezDrugiego,
     1393                                                        aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
     1394                                                        1 - iBezDrugiego,
     1395                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
     1396                                                DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
     1397                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
     1398
    15551399                                                        // zmniejsz liczby niedopasowanych elementow w grupach
    15561400                                                        aiRozmiarGrupy[0]--;
    15571401                                                aiRozmiarGrupy[1]--;
    15581402
    1559                                         }
    1560                                         else
    1561                                         {
    1562                                                 // jedna z grup sie juz wyczerpala
    1563                                                 // wiec nie ma mozliwosci dopasowania tego drugiego Partu
    1564                                                 /// i trzeba poczekac na zmiane grupy
    1565                                         }
    1566 
    1567                                         DB(printf("Przypadek 2.\n");)
     1403                                                // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
     1404                                                // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
     1405                                                if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
     1406                                                {
     1407                                                        // jesli grupy sie jeszcze nie wyczrpaly
     1408                                                        // to jest mozliwosc dopasowania w organizmie
     1409
     1410                                                        int iZDrugim = 1 - iBezDrugiego;
     1411                                                        // sprawdz indeks organizmu
     1412                                                        assert((iZDrugim == 0) || (iZDrugim == 1));
     1413
     1414                                                        // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
     1415                                                        // poprzednie obliczenia minimum
     1416                                                        // dla organizmu 1 (o rozmiarze grupy z 0)
     1417                                                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
     1418                                                        {
     1419                                                                apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
     1420                                                        }
     1421                                                        // dla organizmu 0 (o rozmiarze grup z 1)
     1422                                                        for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
     1423                                                        {
     1424                                                                apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
     1425                                                        }
     1426
     1427                                                        // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
     1428                                                        // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
     1429                                                        dMinimum = HUGE_VAL;
     1430                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1431                                                        {
     1432                                                                if (!(m_pMatching->IsMatched(
     1433                                                                        1 - iZDrugim,
     1434                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
     1435                                                                {
     1436                                                                        // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
     1437                                                                        // niedopasowanych jeszcze Parts
     1438                                                                        if (iZDrugim == 0)
     1439                                                                        {
     1440                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
     1441                                                                                // zeby moc indeksowac niesymetryczna tablice
     1442                                                                                if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
     1443                                                                                {
     1444                                                                                        dMinimum = aadOdleglosciParts[iIndex[0]][i];
     1445                                                                                }
     1446                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
     1447                                                                                assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
     1448
     1449                                                                        }
     1450                                                                        else
     1451                                                                        {
     1452                                                                                if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
     1453                                                                                {
     1454                                                                                        dMinimum = aadOdleglosciParts[i][iIndex[1]];
     1455                                                                                }
     1456                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
     1457                                                                                assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
     1458                                                                        }
     1459                                                                }
     1460                                                        }
     1461                                                        // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
     1462                                                        assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
     1463
     1464                                                        // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
     1465                                                        // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
     1466                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1467                                                        {
     1468                                                                if (!(m_pMatching->IsMatched(
     1469                                                                        1 - iZDrugim,
     1470                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
     1471                                                                {
     1472                                                                        if (iZDrugim == 0)
     1473                                                                        {
     1474                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
     1475                                                                                // zeby moc indeksowac niesymetryczna tablice
     1476                                                                                if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
     1477                                                                                {
     1478                                                                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
     1479                                                                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
     1480                                                                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
     1481                                                                                }
     1482                                                                        }
     1483                                                                        else
     1484                                                                        {
     1485                                                                                if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
     1486                                                                                {
     1487                                                                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
     1488                                                                                }
     1489                                                                        }
     1490                                                                }
     1491                                                        }
     1492
     1493                                                        // a teraz - po znalezieniu potencjalnych elementow do dopasowania
     1494                                                        // dopasowujemy pierwszy z potencjalnych
     1495                                                        iDopasowywany = -1;
     1496                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1497                                                        {
     1498                                                                if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
     1499                                                                {
     1500                                                                        iDopasowywany = i;
     1501                                                                        break;
     1502                                                                }
     1503                                                        }
     1504                                                        // musielismy znalezc jakiegos dopasowywanego
     1505                                                        assert((iDopasowywany >= 0) &&
     1506                                                                (iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
     1507
     1508                                                        // no to juz mozemy dopasowac
     1509                                                        m_pMatching->Match(
     1510                                                                iZDrugim,
     1511                                                                aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
     1512                                                                1 - iZDrugim,
     1513                                                                aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
     1514                                                        DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
     1515
     1516                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1517                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1518                                                                // zmniejsz liczby niedopasowanych elementow w grupach
     1519                                                                aiRozmiarGrupy[0]--;
     1520                                                        aiRozmiarGrupy[1]--;
     1521
     1522                                                }
     1523                                                else
     1524                                                {
     1525                                                        // jedna z grup sie juz wyczerpala
     1526                                                        // wiec nie ma mozliwosci dopasowania tego drugiego Partu
     1527                                                        /// i trzeba poczekac na zmiane grupy
     1528                                                }
     1529
     1530                                                DB(printf("Przypadek 2.\n");)
    15681531
    15691532                                        }// PRZYPADEK 2.
    … …  
    17241687                if (m_adFactors[3] > 0)
    17251688                {
    1726                 // tutaj zacznij pętlę po przekształceniach  geometrycznych
    1727                 const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
    1728                 // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
    1729                 // pomiędzy transformacjami;
    1730                 // wartości orginalne transformacji dOrig uzyskuje się przez:
    1731                 // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
    1732                 //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
    1733                 const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
    1734                 const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
    1735                 const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
     1689                        // tutaj zacznij pętlę po przekształceniach  geometrycznych
     1690                        const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
     1691                        // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
     1692                        // pomiędzy transformacjami;
     1693                        // wartości orginalne transformacji dOrig uzyskuje się przez:
     1694                        // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
     1695                        //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
     1696                        const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
     1697                        const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
     1698                        const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
    17361699
    17371700#ifdef max
    17381701#undef max //this macro would conflict with line below
    17391702#endif
    1740                 double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
    1741                 int iMinSimTransform = -1; // index of the transformation with the minimum similarity
    1742                 SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
    1743 
    1744                 // remember the original positions of model 0 as computed by SVD in order to restore them later, after
    1745                 // all transformations have been computed
    1746                 Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
    1747                 // create the stored positions
    1748                 StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
    1749                 assert(StoredPositions != NULL);
    1750                 // copy the original positions of model 0 (store them)
    1751                 int iPart; // a counter of Parts
    1752                 for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1753                 {
    1754                         StoredPositions[iPart].x = m_aPositions[0][iPart].x;
    1755                         StoredPositions[iPart].y = m_aPositions[0][iPart].y;
    1756                         StoredPositions[iPart].z = m_aPositions[0][iPart].z;
    1757                 }
    1758                 // now the original positions of model 0 are stored
    1759 
    1760 
    1761                 int iTransform; // a counter of geometric transformations
    1762                 for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
    1763                 {
    1764                         // for each geometric transformation to be done
    1765                         // entry conditions:
    1766                         // - models (m_Mod) exist and are available
    1767                         // - matching (m_pMatching) exists and is empty
    1768                         // - all properties are created and available (m_aDegrees and m_aPositions)
    1769 
    1770                         // recompute geometric properties according to the transformation iTransform
    1771                         // but only for model 0
     1703                        double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
     1704                        int iMinSimTransform = -1; // index of the transformation with the minimum similarity
     1705                        SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
     1706
     1707                        // remember the original positions of model 0 as computed by SVD in order to restore them later, after
     1708                        // all transformations have been computed
     1709                        Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
     1710                        // create the stored positions
     1711                        StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
     1712                        assert(StoredPositions != NULL);
     1713                        // copy the original positions of model 0 (store them)
     1714                        int iPart; // a counter of Parts
    17721715                        for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    17731716                        {
    1774                                 // for each iPart, a part of the model iMod
    1775                                 m_aPositions[0][iPart].x *= dMulX[iTransform];
    1776                                 m_aPositions[0][iPart].y *= dMulY[iTransform];
    1777                                 m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1717                                StoredPositions[iPart].x = m_aPositions[0][iPart].x;
     1718                                StoredPositions[iPart].y = m_aPositions[0][iPart].y;
     1719                                StoredPositions[iPart].z = m_aPositions[0][iPart].z;
    17781720                        }
    1779                         // now the positions are recomputed
    1780                         ComputeMatching();
    1781 
    1782                         // teraz m_pMatching istnieje i jest pełne
    1783                         assert(m_pMatching != NULL);
    1784                         assert(m_pMatching->IsFull() == true);
    1785 
    1786                         // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
    1787                         int iTempRes = CountPartsDistance();
    1788                         // załóż sukces
    1789                         assert(iTempRes == 1);
    1790                         double dCurrentSim = m_adFactors[0] * double(m_iDV) +
    1791                                 m_adFactors[1] * double(m_iDD) +
    1792                                 m_adFactors[2] * double(m_iDN) +
    1793                                 m_adFactors[3] * double(m_dDG);
    1794                         // załóż poprawną wartość podobieństwa
    1795                         assert(dCurrentSim >= 0.0);
    1796 
    1797                         // porównaj wartość obliczoną z dotychczasowym minimum
    1798                         if (dCurrentSim < dMinSimValue)
    1799                         {
    1800                                 // jeśli uzyskano mniejszą wartość dopasowania,
    1801                                 // to zapamiętaj to przekształcenie geometryczne i dopasowanie
    1802                                 dMinSimValue = dCurrentSim;
    1803                                 iMinSimTransform = iTransform;
    1804                                 SAFEDELETE(pMinSimMatching);
    1805                                 pMinSimMatching = new SimilMatching(*m_pMatching);
    1806                                 assert(pMinSimMatching != NULL);
     1721                        // now the original positions of model 0 are stored
     1722
     1723
     1724                        int iTransform; // a counter of geometric transformations
     1725                        for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
     1726                        {
     1727                                // for each geometric transformation to be done
     1728                                // entry conditions:
     1729                                // - models (m_Mod) exist and are available
     1730                                // - matching (m_pMatching) exists and is empty
     1731                                // - all properties are created and available (m_aDegrees and m_aPositions)
     1732
     1733                                // recompute geometric properties according to the transformation iTransform
     1734                                // but only for model 0
     1735                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1736                                {
     1737                                        // for each iPart, a part of the model iMod
     1738                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
     1739                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
     1740                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1741                                }
     1742                                // now the positions are recomputed
     1743                                ComputeMatching();
     1744
     1745                                // teraz m_pMatching istnieje i jest pełne
     1746                                assert(m_pMatching != NULL);
     1747                                assert(m_pMatching->IsFull() == true);
     1748
     1749                                // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
     1750                                int iTempRes = CountPartsDistance();
     1751                                // załóż sukces
     1752                                assert(iTempRes == 1);
     1753                                double dCurrentSim = m_adFactors[0] * double(m_iDV) +
     1754                                        m_adFactors[1] * double(m_iDD) +
     1755                                        m_adFactors[2] * double(m_iDN) +
     1756                                        m_adFactors[3] * double(m_dDG);
     1757                                // załóż poprawną wartość podobieństwa
     1758                                assert(dCurrentSim >= 0.0);
     1759
     1760                                // porównaj wartość obliczoną z dotychczasowym minimum
     1761                                if (dCurrentSim < dMinSimValue)
     1762                                {
     1763                                        // jeśli uzyskano mniejszą wartość dopasowania,
     1764                                        // to zapamiętaj to przekształcenie geometryczne i dopasowanie
     1765                                        dMinSimValue = dCurrentSim;
     1766                                        iMinSimTransform = iTransform;
     1767                                        SAFEDELETE(pMinSimMatching);
     1768                                        pMinSimMatching = new SimilMatching(*m_pMatching);
     1769                                        assert(pMinSimMatching != NULL);
     1770                                }
     1771
     1772                                // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
     1773                                m_pMatching->Empty();
     1774                        } // for ( iTransform )
     1775
     1776                        // po pętli przywróć najlepsze dopasowanie
     1777                        delete m_pMatching;
     1778                        m_pMatching = pMinSimMatching;
     1779
     1780                        DB(printf("Matching has been chosen!\n");)
     1781                                // print the name of the chosen transformation:
     1782                                // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
     1783
     1784                                // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
     1785                                // - najpierw przywroc Parts pozycje orginalne, po SVD
     1786                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1787                                {
     1788                                        m_aPositions[0][iPart].x = StoredPositions[iPart].x;
     1789                                        m_aPositions[0][iPart].y = StoredPositions[iPart].y;
     1790                                        m_aPositions[0][iPart].z = StoredPositions[iPart].z;
     1791                                }
     1792                        // - usun teraz zapamietane pozycje Parts
     1793                        delete[] StoredPositions;
     1794                        // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
     1795                        for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
     1796                        {
     1797                                // for each transformation before and INCLUDING iMinTransform
     1798                                // do the transformation (only model 0)
     1799                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1800                                {
     1801                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
     1802                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
     1803                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1804                                }
    18071805                        }
    1808 
    1809                         // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
    1810                         m_pMatching->Empty();
    1811                 } // for ( iTransform )
    1812 
    1813                 // po pętli przywróć najlepsze dopasowanie
    1814                 delete m_pMatching;
    1815                 m_pMatching = pMinSimMatching;
    1816 
    1817                 DB(printf("Matching has been chosen!\n");)
    1818                         // print the name of the chosen transformation:
    1819                         // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
    1820 
    1821                         // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
    1822                         // - najpierw przywroc Parts pozycje orginalne, po SVD
    1823                         for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1824                         {
    1825                         m_aPositions[0][iPart].x = StoredPositions[iPart].x;
    1826                         m_aPositions[0][iPart].y = StoredPositions[iPart].y;
    1827                         m_aPositions[0][iPart].z = StoredPositions[iPart].z;
    1828                         }
    1829                 // - usun teraz zapamietane pozycje Parts
    1830                 delete[] StoredPositions;
    1831                 // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
    1832                 for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
    1833                 {
    1834                         // for each transformation before and INCLUDING iMinTransform
    1835                         // do the transformation (only model 0)
    1836                         for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1837                         {
    1838                                 m_aPositions[0][iPart].x *= dMulX[iTransform];
    1839                                 m_aPositions[0][iPart].y *= dMulY[iTransform];
    1840                                 m_aPositions[0][iPart].z *= dMulZ[iTransform];
    1841                         }
    1842                 }
    18431806
    18441807                }
    … …  
    18551818
    18561819        DB(_PrintPartsMatching();)
    1857 
    18581820
    18591821                return 1;
    … …  
    19851947                int nSize = m_Mod[iMod]->getPartCount();
    19861948
    1987                 double *pDistances =  new double[nSize * nSize];
     1949                double *pDistances = new double[nSize * nSize];
    19881950
    19891951                for (int i = 0; i < nSize; i++)
    … …  
    19971959                Pt3D P1Pos, P2Pos; // positions of parts
    19981960                double dDistance; // the distance between Parts
    1999                
     1961
    20001962                double *weights = new double[nSize];
    20011963                for (int i = 0; i < nSize; i++)
    20021964                {
    2003                         if (wMDS==1)
     1965                        if (wMDS == 1)
    20041966                                weights[i] = 0;
    20051967                        else
    20061968                                weights[i] = 1;
    20071969                }
    2008                                
    2009                 if (wMDS==1)
     1970
     1971                if (wMDS == 1)
    20101972                        for (int i = 0; i < pModel->getJointCount(); i++)
    20111973                        {
    20121974                                weights[pModel->getJoint(i)->p1_refno]++;
    2013                                 weights[pModel->getJoint(i)->p2_refno]++; 
     1975                                weights[pModel->getJoint(i)->p2_refno]++;
    20141976                        }
    2015                
     1977
    20161978                for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
    20171979                {
    … …  
    20572019}
    20582020
     2021/**     Evaluates distance between two given genotypes. The distance depends strongly
     2022        on weights set and the matching algorithm used.
     2023        @param G0 Pointer to the first of compared genotypes
     2024        @param G1 Pointer to the second of compared genotypes.
     2025        @return Distance between two genotypes.
     2026        @sa m_adFactors, matching_method
     2027        */
     2028double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
     2029{
     2030        return matching_method == 0 ? EvaluateDistanceHungarian(G0, G1) : EvaluateDistanceGreedy(G0, G1);
     2031}
     2032
    20592033void ModelSimil::p_evaldistance(ExtValue *args, ExtValue *ret)
    20602034{
    … …  
    20662040                ret->setDouble(EvaluateDistance(g1, g2));
    20672041}
     2042
     2043void ModelSimil::FillPartsDistances(double*& dist, int bigger, int smaller, bool geo)
     2044{
     2045        for (int i = 0; i < bigger; i++)
     2046        {
     2047                for (int j = 0; j < bigger; j++)
     2048                {
     2049                        // assign penalty for unassignment for vertex from bigger model
     2050                        if (j >= smaller)
     2051                        {
     2052                                if (geo)
     2053                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].length();
     2054                                else
     2055                                        dist[i*bigger + j] = m_adFactors[1] * m_aDegrees[1 - m_iSmaller][i][DEGREE] + m_adFactors[2] * m_aDegrees[1 - m_iSmaller][i][NEURONS];
     2056                        }
     2057                        // compute distance between parts
     2058                        else
     2059                        {
     2060                                if (geo)
     2061                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].distanceTo(m_aPositions[m_iSmaller][j]);
     2062                                else
     2063                                        dist[i*bigger + j] = m_adFactors[1] * abs(m_aDegrees[1 - m_iSmaller][i][DEGREE] - m_aDegrees[m_iSmaller][j][DEGREE])
     2064                                        + m_adFactors[2] * abs(m_aDegrees[1 - m_iSmaller][i][NEURONS] - m_aDegrees[m_iSmaller][j][NEURONS]);
     2065                        }
     2066
     2067                }
     2068        }
     2069}
     2070
     2071double ModelSimil::EvaluateDistanceHungarian(const Geno *G0, const Geno *G1)
     2072{
     2073        double dResult = 0;
     2074
     2075        m_Gen[0] = G0;
     2076        m_Gen[1] = G1;
     2077
     2078        // check whether pointers are not NULL
     2079        if (m_Gen[0] == NULL || m_Gen[1] == NULL)
     2080        {
     2081                DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid genotypes pointers\n");)
     2082                        return 0.0;
     2083        }
     2084        // create models of objects to compare
     2085        m_Mod[0] = new Model(*(m_Gen[0]));
     2086        m_Mod[1] = new Model(*(m_Gen[1]));
     2087
     2088        // validate models
     2089        if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
     2090        {
     2091                DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid models pointers\n");)
     2092                        return 0.0;
     2093        }
     2094
     2095        //Get information about vertex degrees, neurons and 3D location
     2096        if (!CreatePartInfoTables())
     2097                return 0;
     2098        if (!CountPartDegrees())
     2099                return 0;
     2100        if (!GetPartPositions())
     2101                return 0;
     2102        if (!CountPartNeurons())
     2103                return 0;
     2104
     2105        m_iSmaller = m_Mod[0]->getPartCount() <= m_Mod[1]->getPartCount() ? 0 : 1;
     2106        int nSmaller = m_Mod[m_iSmaller]->getPartCount();
     2107        int nBigger = m_Mod[1 - m_iSmaller]->getPartCount();
     2108
     2109        double* partsDistances = new double[nBigger*nBigger]();
     2110        FillPartsDistances(partsDistances, nBigger, nSmaller, false);
     2111        int *assignment = new int[nBigger]();
     2112
     2113        HungarianAlgorithm hungarian;
     2114
     2115        if (m_adFactors[3] > 0)
     2116        {
     2117                if (!ComputePartsPositionsBySVD())
     2118                {
     2119                        return 0;
     2120                }
     2121
     2122                // tutaj zacznij pętlę po przekształceniach  geometrycznych
     2123                const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
     2124                // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
     2125                // pomiędzy transformacjami;
     2126                const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
     2127                const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
     2128                const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
     2129
     2130                std::vector<int> minAssignment(nBigger);
     2131#ifdef max
     2132#undef max //this macro would conflict with line below
     2133#endif
     2134                double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
     2135
     2136                int iTransform; // a counter of geometric transformations
     2137                for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
     2138                {
     2139                        // for each geometric transformation to be done
     2140                        // entry conditions:
     2141                        // - models (m_Mod) exist and are available
     2142                        // - all properties are created and available (m_aDegrees and m_aPositions)
     2143                        double* tmpPartsDistances = new double[nBigger*nBigger]();
     2144                        std::copy(partsDistances, partsDistances + nBigger * nBigger, tmpPartsDistances);
     2145                        // recompute geometric properties according to the transformation iTransform
     2146                        // but only for model 0
     2147                        for (int iPart = 0; iPart < m_Mod[m_iSmaller]->getPartCount(); iPart++)
     2148                        {
     2149                                // for each iPart, a part of the model iMod
     2150                                m_aPositions[m_iSmaller][iPart].x *= dMulX[iTransform];
     2151                                m_aPositions[m_iSmaller][iPart].y *= dMulY[iTransform];
     2152                                m_aPositions[m_iSmaller][iPart].z *= dMulZ[iTransform];
     2153                        }
     2154                        // now the positions are recomputed
     2155
     2156                        FillPartsDistances(tmpPartsDistances, nBigger, nSmaller, true);
     2157                        std::fill_n(assignment, nBigger, 0);
     2158                        double dCurrentSim = hungarian.Solve(tmpPartsDistances, assignment, nBigger, nBigger);
     2159
     2160                        delete[] tmpPartsDistances;
     2161                        // załóż poprawną wartość podobieństwa
     2162                        assert(dCurrentSim >= 0.0);
     2163
     2164                        // porównaj wartość obliczoną z dotychczasowym minimum
     2165                        if (dCurrentSim < dMinSimValue)
     2166                        {
     2167                                dMinSimValue = dCurrentSim;
     2168                                if (saveMatching == 1)
     2169                                {
     2170                                        minAssignment.clear();
     2171                                        minAssignment.insert(minAssignment.begin(), assignment, assignment + nBigger);
     2172                                }
     2173                        }
     2174                }
     2175
     2176                dResult = dMinSimValue;
     2177                if (saveMatching == 1)
     2178                        std::copy(minAssignment.begin(), minAssignment.end(), assignment);
     2179        }
     2180
     2181        else
     2182        {
     2183                dResult = hungarian.Solve(partsDistances, assignment, nBigger, nBigger);
     2184        }
     2185
     2186        //add difference in anywhere and onJoint neurons
     2187        dResult += m_adFactors[2] * (abs(m_aOnJoint[0][3] - m_aOnJoint[1][3]) + abs(m_aAnywhere[0][3] - m_aAnywhere[1][3]));
     2188        //add difference in part numbers
     2189        dResult += (nBigger - nSmaller) * m_adFactors[0];
     2190
     2191        // delete degree arrays created in CreatePartInfoTables
     2192        SAFEDELETEARRAY(m_aDegrees[0]);
     2193        SAFEDELETEARRAY(m_aDegrees[1]);
     2194
     2195        // and position arrays
     2196        SAFEDELETEARRAY(m_aPositions[0]);
     2197        SAFEDELETEARRAY(m_aPositions[1]);
     2198
     2199        // delete created models
     2200        SAFEDELETE(m_Mod[0]);
     2201        SAFEDELETE(m_Mod[1]);
     2202
     2203        delete[] assignment;
     2204        delete[] partsDistances;
     2205
     2206        return dResult;
     2207}
Note: See TracChangeset for help on using the changeset viewer.