Main Page Namespace List Class Hierarchy Compound List File List Namespace Members Compound Members File Members Related Pages

# SGDistMeasurer.cpp

Go to the documentation of this file.

    00001
    00032 #include "SGDistMeasurer.h"
    00033 #include "ShockGraph.h"
    00034 #include "DAG.h"
    00035 #include "HelperFunctions.h"
    00036 #include "SGNode.h"
    00037 #include <float.h>
    00038 #include <stdio.h>
    00039 #include <ctype.h>
    00040 #include <assert.h>
    00041
    00042 #include "ApproxPoly.h"
    00043
    00044 #ifdef WIN32
    00045 #define M_PI 3.14159265359
    00046 #endif
    00047
    00048 #define MINDIST 0
    00049 #define MAXDIST 1
    00050
    00052 // SGDistMeasurer3
    00053
    00054 void SGDistMeasurer3::LogResults(const ModelFit& m1,
    00055         const ModelFit& m2, POINTS g1Pts, POINTS g2Pts, bool bReverseOrder) const
    00056 {
    00057         static bool bFirstTime = true;
    00058         static fstream log;
    00059
    00060         //Log the results if debug mode. In any case, empty the log file
    00061         if (bFirstTime)
    00062         {
    00063                 bFirstTime = false;
    00064                 log.open("nodedist.m", ios::out | ios::trunc);
    00065         }
    00066         //else
    00067         //      log.open("nodedist.m", ios::out | ios::app);
    00068
    00069         int d0, dN;
    00070
    00071         // If reverted, put them back in place before displaying
    00072         if (bReverseOrder)
    00073         {
    00074                 g1Pts = pNode1->GetVelocityRadiusArray(d0, dN, false);
    00075                 g2Pts = pNode2->GetVelocityRadiusArray(d0, dN, false);
    00076         }
    00077
    00078         log << "disp('" << pG1->GetDAGLbl() << "');\n";
    00079         log << "disp('" << pG2->GetDAGLbl() << "');\n";
    00080
    00081         log << "disp('" << pNode1->GetNodeLbl()
    00082                 << " vs " << pNode2->GetNodeLbl()
    00083                 << " [DFS indices: " << pG1->GetNodeDFSIndex(v1) << " vs " << pG2->GetNodeDFSIndex(v2) << "]');\n";
    00084         log << "% SEARCH CODE: " << pG1->GetViewNumber() << pNode1->GetNodeLbl()
    00085                 << pG2->GetViewNumber() << pNode2->GetNodeLbl() << "\n";
    00086
    00087         m1.Plot(log, g2Pts, pNode2->GetSegments());
    00088
    00089         log << "disp('" << pG1->GetDAGLbl() << "');\n";
    00090         log << "disp('" << pG2->GetDAGLbl() << "');\n";
    00091
    00092         log << "disp('" << pNode2->GetNodeLbl()
    00093                 << " vs " << pNode1->GetNodeLbl()
    00094                 << " [DFS indices: " << pG2->GetNodeDFSIndex(v2) << " vs " << pG1->GetNodeDFSIndex(v1) << "]');\n";
    00095
    00096         m2.Plot(log, g1Pts, pNode1->GetSegments());
    00097
    00098         //log << "disp('Avg error: " << error << "');\n";
    00099 }
    00100
    00104 double SGDistMeasurer3::CompDistance() const
    00105 {
    00106         /*if (pG1->GetViewNumber() == 79 && pG2->GetViewNumber() == 6 &&
    00107                 pG1->GetNodeDFSIndex(v1) == 1 && pG2->GetNodeDFSIndex(v2) == 1)
    00108         {
    00109                 cout << "\n" << pNode1->GetNodeLbl() << "\n" << pNode2->GetNodeLbl() << "\n";
    00110         }*/
    00111
    00112         int n1Type = pG1->NodeType(v1);
    00113         int n2Type = pG2->NodeType(v2);
    00114
    00115         // If label is different, do not compare nodes
    00116         if (n1Type != n2Type)
    00117                 return MAXDIST;
    00118         else if (n1Type == ROOT || n2Type == ROOT)
    00119                 return n1Type == n2Type ? MINDIST:MAXDIST;
    00120         else if (pNode1->m_shocks.GetSize() == 1 && pNode2->m_shocks.GetSize() == 1)
    00121         {
    00122                 double r1 = pNode1->m_shocks.GetHead().radius;
    00123                 double r2 = pNode2->m_shocks.GetHead().radius;
    00124
    00125                 return 1 - MIN(r1, r2) / MAX(r1, r2);
    00126         }
    00127         else if (pNode1->m_shocks.GetSize() == 1 || pNode2->m_shocks.GetSize() == 1)
    00128         {
    00129                 if (pNode1->m_shocks.GetSize() > 4 || pNode2->m_shocks.GetSize() > 4)
    00130                         return MAXDIST; // if it is too long, it's unlikely that they are the same
    00131
    00132                 double r1 = pNode1->m_shocks.AvgRadius();
    00133                 double r2 = pNode2->m_shocks.AvgRadius();
    00134                 double l1 = pNode1->m_shocks.Length();
    00135                 double l2 = pNode2->m_shocks.Length();
    00136
    00137                 double rdiff = 1 - MIN(r1, r2) / MAX(r1, r2);
    00138                 double lendiff = 1 - MIN(l1, l2) / MAX(l1, l2);
    00139
    00140                 return 0.85 * rdiff + 0.15 * lendiff;
    00141         }
    00142
    00143         ASSERT(pNode1->GetSegments().GetSize() > 0 && pNode2->GetSegments().GetSize() > 0);
    00144
    00145         // We are dealing with same label, non-root nodes, with more that 1 shock point
    00146         // so...
    00147         leda_node par1 = pG1->GetFirstParent(v1);
    00148         leda_node par2 = pG2->GetFirstParent(v2);
    00149         int n1Dir, n2Dir;
    00150
    00151         if (pG1->NodeType(par1) != ROOT && pG2->NodeType(par2) != ROOT)
    00152         {
    00153                 n1Dir = pG1->GetBranchDir(v1, par1);
    00154                 n2Dir = pG2->GetBranchDir(v2, par2);
    00155
    00156         // If the nodes' directions wrt their parents differ,
    00157         // we can say they directions don't match
    00158         if (n1Dir != n2Dir)
    00159                 return MAXDIST;
    00160         }
    00161
    00162         // Now, given that their directions match, we need to make
    00163         // sure that we get them right. The may go from 0 to N or form N to 0.
    00164         n1Dir = pNode1->GetBranchDir(ShockBranch::FORWARD);
    00165         n2Dir = pNode2->GetBranchDir(ShockBranch::FORWARD);
    00166         bool bReverseOrder = n1Dir != n2Dir;
    00167
    00168         // what if dir == 0 but there is a slope? We should still reverse it.
    00169         // THIS IS A TO DO
    00170
    00171         int d0, dN;
    00172         POINTS g1Pts = pNode1->GetVelocityRadiusArray(d0, dN, bReverseOrder);
    00173         POINTS g2Pts = pNode2->GetVelocityRadiusArray(d0, dN, bReverseOrder);
    00174
    00175         ASSERT( pNode1->GetSegments().GetHead().p0.y == (bReverseOrder ? g1Pts.GetTail().y:g1Pts.GetHead().y) );
    00176         ASSERT( pNode1->GetSegments().GetTail().p1.y == (bReverseOrder ? g1Pts.GetHead().y:g1Pts.GetTail().y) );
    00177         ASSERT( pNode2->GetSegments().GetHead().p0.y == (bReverseOrder ? g2Pts.GetTail().y:g2Pts.GetHead().y) );
    00178         ASSERT( pNode2->GetSegments().GetTail().p1.y == (bReverseOrder ? g2Pts.GetHead().y:g2Pts.GetTail().y) );
    00179
    00180         // Let's fit the model of node one to the data of node 2 and viceversa
    00181         ModelFit m1(2), m2(2); // 2 is max diff in lenght
    00182
    00183         double e1 = m1.Fit(g1Pts, pNode2->GetSegments());
    00184
    00185         if (e1 >= INFINITY)
    00186                 return MAXDIST;
    00187
    00188         double e2 = m2.Fit(g2Pts, pNode1->GetSegments());
    00189
    00190         if (e2 >= INFINITY)
    00191                 return MAXDIST;
    00192
    00193         if (DAG::IsDbgMode())
    00194                 LogResults(m1, m2, g1Pts, g2Pts, bReverseOrder);
    00195
    00196         /*double e, s;
    00197
    00198         for (int i = 0; i < pNode2->GetSegments().GetSize(); i++)
    00199         {
    00200                 m1.GetSegmentError(i, &e, &s);
    00201                 DBG_SHOW(e)
    00202                 DBG_SHOW(s)
    00203         }*/
    00204
    00205         // We have non-infinite errors, so we must normalize them
    00206         double avge1 = (e1 / g1Pts.GetSize()) / pNode1->m_shocks.AvgRadius();
    00207         double avge2 = (e2 / g2Pts.GetSize()) / pNode2->m_shocks.AvgRadius();
    00208         double error = (avge1 + avge2) / 2.0;
    00209
    00210         if (error >= MAXDIST)
    00211                 return MAXDIST;
    00212
    00213         double minL = MIN(pNode1->m_shocks.Length(), pNode2->m_shocks.Length());
    00214         double maxL = MAX(pNode1->m_shocks.Length(), pNode2->m_shocks.Length());
    00215
    00216         return 0.85 * error + 0.15 * (1 - minL / maxL);
    00217 }
    00218
    00220 // SGDistMeasurer1
    00221
    00222 double SGDistMeasurer1::CompDistance() const
    00223 {
    00224         double level1 = NodeDistanceLevel1();
    00225         double level2 = 0;
    00226
    00227         //      double level2 = NodeDistanceLevel2(v1,v2);
    00228
    00229         /*if(pG1->GetNodeLbl(v1) == pG2->GetNodeLbl(v2)){
    00230 //      cout << pG1->GetNodeLbl(v1) << "<>" << pG2->GetNodeLbl(v2) << endl;
    00231 //      cout << level1 << ":" <<level2 << endl;
    00232 }*/
    00233
    00234         //double value = w1 * (level1) + w2 * (level2);
    00235         double value = level1;
    00236
    00237         //cout << " s = " << value;
    00238
    00239         ASSERT(value >= 0 && value <= 1);
    00240
    00241         return value;
    00242 }
    00243
    00244 /*
    00245 Takes two numbers which represent the relative weights
    00246 of level1 vs level2 in the node distance calculation.
    00247 */
    00248 void SGDistMeasurer1::SetWeights(double level1Weight, double level2Weight)
    00249 {
    00250         ASSERT(level1Weight + level2Weight == 1);
    00251
    00252         w1 = level1Weight;
    00253         w2 = level2Weight;
    00254 }
    00255
    00256 double SGDistMeasurer1::AverageRadius(const SGNode *sb) const
    00257 {
    00258         double min_rad = 999999;
    00259         double max_rad = 0;
    00260         double avg_rad = 0;
    00261         int shockCount;
    00262
    00263         shockCount = sb->GetShockCount();
    00264
    00265         for (int j = 0; j < shockCount; j++)
    00266     {
    00267                 double curr_rad = sb->m_shocks[j].radius;
    00268                 if (curr_rad > max_rad)
    00269                         max_rad = curr_rad;
    00270                 if (curr_rad < min_rad)
    00271                         min_rad = curr_rad;
    00272
    00273                 avg_rad += curr_rad;
    00274     }
    00275
    00276         if (shockCount != 0)
    00277                 avg_rad /= shockCount;
    00278         return avg_rad;
    00279 }
    00280
    00281 double SGDistMeasurer1::Length(const SGNode *sb) const
    00282 {
    00283         int shockCount = sb->GetShockCount();
    00284
    00285         double distance = 0;
    00286         double x1,x2,y1,y2;
    00287
    00288         for(int i = 1; i < shockCount; i++)
    00289     {
    00290                 x1 =  sb->m_shocks[i-1].xcoord;
    00291                 y1 =  sb->m_shocks[i-1].ycoord;
    00292
    00293                 x2 =  sb->m_shocks[i].xcoord;
    00294                 y2 =  sb->m_shocks[i].ycoord;
    00295
    00296                 distance += sqrt(pow(x2-x1,2.0) + pow(y2-y1,2.0));
    00297     }
    00298
    00299         return distance;
    00300 }
    00301
    00302 double SGDistMeasurer1::Elongation(const SGNode *sb) const
    00303 {
    00304         double avg_radius_sb = AverageRadius(sb);
    00305         double length = Length(sb);
    00306         double elongation;
    00307
    00308
    00309         if(length != 0)
    00310           {
    00311             // This just normalizes the value.
    00312             elongation = (avg_radius_sb/length)/(avg_radius_sb + length);
    00313           }
    00314         else
    00315           elongation = avg_radius_sb;
    00316         //      cout << "elongation = " << elongation << endl;
    00317         return elongation;
    00318
    00319 }
    00320
    00321 double SGDistMeasurer1::Bend(const SGNode *u) const
    00322 {
    00323         double x1,x2,y1,y2;
    00324         int numSinModel = u->GetShockCount();
    00325
    00326         if(numSinModel <= 3)
    00327                 return M_PI;
    00328
    00329         x1 =  u->m_shocks[numSinModel/2].xcoord - u->m_shocks[0].xcoord;
    00330         x2 =  u->m_shocks[numSinModel-1].xcoord - u->m_shocks[numSinModel/2].xcoord;
    00331
    00332         y1 =  u->m_shocks[numSinModel/2].ycoord - u->m_shocks[0].ycoord;
    00333         y2 =  u->m_shocks[numSinModel-1].ycoord - u->m_shocks[numSinModel/2].ycoord;
    00334
    00335
    00336         double value = (x1*x2+y1*y2)/(sqrt(pow(x1,2.0)+pow(y1,2.0))*sqrt(pow(x2,2.0)+pow(y2,2.0)));
    00337         double theta;
    00338
    00339         if(fabs(value) <= 1)
    00340                 theta = acos(value);
    00341         else
    00342                 theta = 0;
    00343
    00344         return theta;
    00345 }
    00346
    00347
    00348 double SGDistMeasurer1::SweepContrast(const SGNode *node) const
    00349 {
    00350
    00351         double min_rad = 999999;
    00352         double max_rad = 0;
    00353         double curr_rad = 0;
    00354         double sweep_contrast = 0.0;
    00355
    00356         int shockCount = node->GetShockCount();
    00357
    00358         if (shockCount <= 0)
    00359                 return 0;
    00360
    00361         double LengthNode = Length(node);
    00362
    00363         for (int j = 0; j < shockCount; j++)
    00364     {
    00365                 double curr_rad = node->m_shocks[j].radius;
    00366
    00367                 if (curr_rad > max_rad)
    00368                         max_rad = curr_rad;
    00369
    00370                 if (curr_rad < min_rad)
    00371                         min_rad = curr_rad;
    00372     }
    00373
    00374         if (max_rad > 0.0)
    00375                 sweep_contrast = ((double) max_rad - (double) min_rad) / (double) max_rad;
    00376
    00377         return(sweep_contrast);
    00378 }
    00379
    00380
    00381 double SGDistMeasurer1::Taper(const SGNode *u) const
    00382 {
    00383
    00384         int numSinNode = u->GetShockCount();
    00385         double length = Length(u);
    00386         double slope;
    00387         double del_radius;
    00388
    00389         if (numSinNode < 2)
    00390                 return 0.0;
    00391
    00392         del_radius = fabs(u->m_shocks[0].radius - u->m_shocks[numSinNode - 1].radius);
    00393
    00394         if (del_radius == 0.0)
    00395         {
    00396                 slope = 0.0;
    00397         }
    00398         else
    00399         {
    00400                 slope = del_radius / length;
    00401         }
    00402         //      cout << "del_radius = " << del_radius << " length = " << length << " slope = " << slope << endl;
    00403
    00404         return(slope);
    00405 }
    00406
    00407 double SGDistMeasurer1::Acceleration(const SGNode *u) const
    00408 {
    00409         int numShocks = u->GetShockCount();
    00410         double accelerationSum = 0;
    00411
    00412         for(int i = 2; i < numShocks; i++)
    00413     {
    00414                 // These are the X and Y coordinates
    00415                 // in the shocks that will be used
    00416                 // to approximate the derivative.
    00417                 double x1 =  u->m_shocks[i-2].xcoord;
    00418                 double x2 =  u->m_shocks[i-1].xcoord;
    00419                 double x3 =  u->m_shocks[i].xcoord;
    00420
    00421                 double y1 =  u->m_shocks[i-2].ycoord;
    00422                 double y2 =  u->m_shocks[i-1].ycoord;
    00423                 double y3 =  u->m_shocks[i].ycoord;
    00424
    00425                 double radius1 =  u->m_shocks[i-2].radius;
    00426                 double radius2 =  u->m_shocks[i-1].radius;
    00427                 double radius3 =  u->m_shocks[i].radius;
    00428
    00429                 double D1 = sqrt(pow(x2-x1,2.0) + pow(y2-y1,2.0));
    00430                 double D2 = sqrt(pow(x3-x2,2.0) + pow(y3-y2,2.0));
    00431
    00432                 double dT1 = radius2 - radius1;
    00433                 double dT2 = radius3 - radius2;
    00434
    00435                 double V1 = D1/dT1;
    00436                 double V2 = D2/dT2;
    00437
    00438                 double A = (V2 - V1)/dT2 - dT2;
    00439
    00440                 accelerationSum += A;
    00441     }
    00442         return accelerationSum;
    00443 }
    00444
    00445 double SGDistMeasurer1::ElongationDiff() const
    00446 {
    00447         double u_elongation = Elongation(pNode1);
    00448         double v_elongation = Elongation(pNode2);
    00449
    00450         double elongation_diff;
    00451
    00452         if(fabs(v_elongation + u_elongation) != 0)
    00453                 elongation_diff = fabs(v_elongation - u_elongation)/fabs(v_elongation + u_elongation);
    00454         else{
    00455                 elongation_diff =  fabs(v_elongation - u_elongation);
    00456         }
    00457
    00458         return elongation_diff;
    00459 }
    00460
    00461 double SGDistMeasurer1::TaperDiff() const
    00462 {
    00463         double taper_u = Taper(pNode1);
    00464         double taper_v = Taper(pNode2);
    00465
    00466         double taper_max, taper_diff;
    00467
    00468         if (taper_u >= taper_v)
    00469         {
    00470                 taper_max = taper_u;
    00471         }
    00472         else
    00473         {
    00474                 taper_max = taper_v;
    00475         }
    00476         if (taper_max == 0)
    00477         {
    00478                 taper_diff = 0;
    00479         }
    00480         else
    00481         {
    00482                 taper_diff = fabs(taper_u - taper_v)/taper_max;
    00483         }
    00484         return(taper_diff);
    00485 }
    00486
    00487 double SGDistMeasurer1::SweepContrastDiff() const
    00488 {
    00489         double u_sweep_contrast = SweepContrast(pNode1);
    00490         double v_sweep_contrast = SweepContrast(pNode2);
    00491         double max_sweep_contrast = 0.0;
    00492
    00493
    00494         double sweep_contrast_diff = 0.0;
    00495
    00496         if (u_sweep_contrast >= v_sweep_contrast)
    00497         {
    00498                 max_sweep_contrast = u_sweep_contrast;
    00499         }
    00500         else
    00501         {
    00502                 max_sweep_contrast = v_sweep_contrast;
    00503         }
    00504
    00505         if((u_sweep_contrast == 0.0) && (v_sweep_contrast == 0.0))
    00506         {
    00507                 sweep_contrast_diff == 0;
    00508         }
    00509         else
    00510         {
    00511                 sweep_contrast_diff = fabs(u_sweep_contrast - v_sweep_contrast)/max_sweep_contrast;
    00512         }
    00513
    00514         return (sweep_contrast_diff);
    00515 }
    00516
    00517 /*
    00518 This takes two nodes and returns a number representing
    00519 the distance between these two nodes in terms of the
    00520 2nd derivative of the radius vs time "graph".
    00521 */
    00522
    00523 double SGDistMeasurer1::AccelerationDiff() const
    00524 {
    00525         double u_acc = Acceleration(pNode1);
    00526         double v_acc = Acceleration(pNode2);
    00527
    00528         double acc_diff;
    00529
    00530         if(fabs(v_acc) + fabs(u_acc) != 0)
    00531                 acc_diff = fabs(v_acc -u_acc)/(fabs(v_acc)+fabs(u_acc));
    00532         else
    00533                 acc_diff =  fabs(v_acc - u_acc);
    00534
    00535         return acc_diff;
    00536 }
    00537
    00538 /*
    00539 This takes two nodes and returns a number representing
    00540 the distance between these two nodes in terms of an
    00541 approximation of an angle that represents it's medial
    00542 axis.
    00543 */
    00544
    00545 double SGDistMeasurer1::BendDiff() const
    00546 {
    00547         double u_bend = Bend(pNode1);
    00548         double v_bend = Bend(pNode2);
    00549
    00550         double bend_diff = fabs(u_bend - v_bend)/M_PI;
    00551
    00552         return bend_diff;
    00553 }
    00554
    00555 /*
    00556 This function is the control function for the first level of
    00557 matching.
    00558 */
    00559 double SGDistMeasurer1::NodeDistanceLevel1() const
    00560 {
    00561         int uType = pG1->NodeType(v1);
    00562         int vType = pG2->NodeType(v2);
    00563
    00564         if (uType == ROOT && vType == ROOT)
    00565                 return 0;
    00566         else if(uType == ROOT || vType == ROOT)
    00567                 return 1;
    00568
    00569         //cout << endl << pG1->GetNodeLbl(v1) << " <> " << pG2->GetNodeLbl(v2);
    00570
    00571         if(uType == SHOCK_1 && vType == SHOCK_1)
    00572     {
    00573                 if (pNode1->GetShockCount() <= 1 || pNode2->GetShockCount() <= 1)
    00574                 {
    00575                         DBG_MSG("Warning: Type 1 node with 1 or 0 shocks!")
    00576                         return 1;
    00577                 }
    00578
    00579                 double taper = TaperDiff();
    00580                 //              double acc = AccelerationDiff();
    00581                 double bend = BendDiff();
    00582                 double elongation = ElongationDiff();
    00583
    00584                 //              cout << "taper diff = " << taper << " bend diff = " << bend << " elongation diff = " << elongation << endl;
    00585
    00586
    00587                 double value = 0.5*taper + 0.2*bend + 0.3*elongation; // + 0.1*acc;
    00588
    00589                 if(value <= 1 && value >= 0)
    00590                         return value;
    00591                 else if(value < 0)
    00592                 {
    00593                         cerr << "Error in the distance calculation type 1 SHOCK!!" << value << endl;
    00594                         assert(false);
    00595                         return -1;
    00596                 }
    00597                 else
    00598                 {
    00599                         cerr << "Error in the normalization type 1 SHOCK!!" << value << endl;
    00600                         assert(false);
    00601                         return -1;
    00602                 }
    00603     }
    00604         else if(uType == SHOCK_3 && vType == SHOCK_3)
    00605     {
    00606                 if (pNode1->GetShockCount() <= 1 || pNode2->GetShockCount() <= 1)
    00607                 {
    00608                         DBG_MSG("Warning: Type 3 node with 1 or 0 shocks!")
    00609                         return 1;
    00610                 }
    00611
    00612                 double bend = BendDiff();
    00613                 double elongation = ElongationDiff();
    00614
    00615                 if(bend <= 1 && bend >= 0)
    00616                 {
    00617                         double value = 0.8*bend + 0.2*elongation;
    00618
    00619                         if(value >= 0 && value <= 1)
    00620                                 return value;
    00621                         else
    00622                         {
    00623                                 cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
    00624                                 assert(false);
    00625                                 return -1;
    00626                         }
    00627                 }
    00628                 else if(bend < 0 || elongation < 0)
    00629                 {
    00630                         cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
    00631                         assert(false);
    00632                         return -1;
    00633                 }
    00634                 else
    00635                 {
    00636                         cerr << "Error in the normalization type 3 SHOCK!!" << bend << endl;
    00637                         assert(false);
    00638                         return -11;
    00639                 }
    00640     }
    00641         else if(uType == SHOCK_2 && vType == SHOCK_2)
    00642     {
    00643                 if (pNode1->GetShockCount() <= 0 || pNode2->GetShockCount() <= 0)
    00644                 {
    00645                         DBG_MSG("Warning: Type 2 node with 0 shocks!")
    00646                         return 1;
    00647                 }
    00648
    00649                 double bend = BendDiff();
    00650                 double elongation = ElongationDiff();
    00651                 double contrast_sweep = SweepContrastDiff();
    00652
    00653                 //              cout << "cont sweep diff = " << contrast_sweep << " bend diff = " << bend << " elongation diff = " << elongation << endl;
    00654
    00655                 if(bend <= 1 && bend >= 0)
    00656                 {
    00657                         double value = 0.3 * bend + 0.2 * elongation + 0.5 * contrast_sweep;
    00658
    00659                         if(value >= 0 && value <= 1)
    00660                                 return value;
    00661                         else
    00662                         {
    00663                                 cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
    00664                                 assert(false);
    00665                                 return -1;
    00666                         }
    00667                 }
    00668                 else if(bend < 0 || elongation < 0)
    00669                 {
    00670                         cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
    00671                         assert(false);
    00672                         return -1;
    00673                 }
    00674                 else{
    00675                         cerr << "Error in the normalization type 2 SHOCK!!" << bend << endl;
    00676                         assert(false);
    00677                         return -1;
    00678                 }
    00679     }
    00680         else if(uType == SHOCK_4 && vType == SHOCK_4)  //updated to compare non-perfect blobs, using same criteria as type 2.
    00681     {
    00682                 double bend = BendDiff();
    00683                 double elongation = ElongationDiff();
    00684
    00685                 if(bend <= 1 && bend >= 0)
    00686                 {
    00687                         double value = 0.8*bend + 0.2*elongation;
    00688
    00689                         if(value >= 0 && value <= 1)
    00690                                 return value;
    00691                         else
    00692                         {
    00693                                 cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
    00694                                 assert(false);
    00695                                 return -1;
    00696                         }
    00697                 }
    00698                 else if(bend < 0 || elongation < 0)
    00699                 {
    00700                         cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
    00701                         assert(false);
    00702                         return -1;
    00703                 }
    00704                 else{
    00705                         cerr << "Error in the normalization type 2 SHOCK!!" << bend << endl;
    00706                         assert(false);
    00707                         return -1;
    00708                 }
    00709     }
    00710         else if((vType == SHOCK_1 && uType == SHOCK_3) || (vType == SHOCK_3 && uType == SHOCK_1))
    00711     {
    00712                 double taper = TaperDiff();
    00713                 double bend = BendDiff();
    00714                 double elongation = ElongationDiff();
    00715
    00716                 double value = 0.5*taper + 0.3*elongation +  0.2*bend;
    00717
    00718                 if(value <= 1)
    00719                         return value;
    00720                 else if(value < 0){
    00721                         cerr << "Error in comparison of TYPE 3 and 1!! " << value << endl;
    00722                         assert(false);
    00723                         return -1;
    00724                 }
    00725                 else
    00726                 {
    00727                         cerr << "Error in normalization of TYPE 3 and 1!! " << value << endl;
    00728                         assert(false);
    00729                         return -1;
    00730                 }
    00731     }
    00732         /*      else if((uType == SHOCK_3 && vType == SHOCK_2) || (uType == SHOCK_2 && vType == SHOCK_3))
    00733     {
    00734         double bend = BendDiff();
    00735         double elongation = ElongationDiff();
    00736
    00737          if(bend <= 1 && bend >= 0)
    00738          {
    00739          double value = 0.8*bend + 0.2*elongation;
    00740
    00741           if(value >= 0 && value <= 1)
    00742           return value;
    00743           else
    00744           {
    00745           cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
    00746           assert(false);
    00747           return -1;
    00748           }
    00749           }
    00750           else if(bend < 0 || elongation < 0)
    00751           {
    00752           cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
    00753           assert(false);
    00754           return -1;
    00755           }
    00756           else
    00757           {
    00758           cerr << "Error in the normalization type 3 SHOCK!!" << bend << endl;
    00759           assert(false);
    00760           return 1;
    00761           }
    00762           }
    00763         */
    00764         else{
    00765
    00766                 return 1;
    00767         }
    00768 }
    00769
    00770 /*
    00771 This function implements rule 1 of the grammar. It tries to quantify how
    00772 close the roots (a and b) are to each other in terms of the sub tree they
    00773 generated.
    00774 */
    00775
    00776 double SGDistMeasurer1::rule1(leda_node a, leda_node b) const
    00777 {
    00778         int aType = pG1->NodeType(a);
    00779         int bType = pG2->NodeType(b);
    00780
    00781         double num_3A = 0, num_3B = 0, num_4A = 0,num_4B = 0;
    00782
    00783         leda_node w,s;
    00784
    00785         if(aType == ROOT &&  bType == ROOT)
    00786     {
    00787                 forall_adj_nodes(w, a)
    00788                 {
    00789                         if (pG1->NodeType(w) == SHOCK_3)
    00790                                 num_3A++;
    00791                         else if(pG1->NodeType(w) == SHOCK_4)
    00792                                 num_4A++;
    00793                 }
    00794
    00795                 forall_adj_nodes(s, b)
    00796                 {
    00797                         if (pG2->NodeType(s) == SHOCK_3)
    00798                                 num_3B++;
    00799                         else if(pG2->NodeType(s) == SHOCK_4)
    00800                                 num_4B++;
    00801                 }
    00802
    00803                 if(num_3A + num_4A == 0){
    00804                         cerr << "Tree has no children of the root" << endl;
    00805                         assert(false);
    00806                         return -1;
    00807                 }
    00808
    00809                 if(num_3B + num_4B == 0){
    00810                         cerr << "Tree has no children of the root" << endl;
    00811                         assert(false);
    00812                         return -1;
    00813                 }
    00814
    00815                 double score;
    00816
    00817                 if(num_3A + num_3B > 0)
    00818                         score = fabs(num_3A - num_3B)/(num_3A + num_3B);
    00819                 else if(num_3A + num_3B == 0)
    00820                 {
    00821                         score = 0;
    00822                 }
    00823                 else
    00824                 {
    00825                         cerr << "Tree has a negative number of children of the root node" << endl;
    00826                         assert(false);
    00827                         return -1;
    00828                 }
    00829
    00830                 if(num_4A + num_4B > 0)
    00831                         score += fabs(num_4A - num_4B)/(num_4A + num_4B);
    00832                 else if(num_4A + num_4B == 0)
    00833                 {
    00834                         score += 0;
    00835                 }
    00836                 else
    00837                 {
    00838                         cerr << "Tree has a negative number of children of the root node" << endl;
    00839                         assert(false);
    00840                         return -1;
    00841                 }
    00842
    00843                 return score/2;
    00844     }
    00845         else
    00846     {
    00847                 return 1;
    00848     }
    00849 }
    00850
    00851 /*
    00852 This function tries to check how closely related both nodes
    00853 that were generated from rule 2 match.
    00854 */
    00855
    00856 double SGDistMeasurer1::rule2(leda_node a, leda_node b)  const
    00857 {
    00858         leda_node w,s;
    00859         int aType = pG1->NodeType(a);
    00860         int bType = pG2->NodeType(b);
    00861
    00862         int num_child_a = 0, num_child_b = 0;
    00863
    00864         if(aType == SHOCK_4 && bType == SHOCK_4)
    00865     {
    00866                 forall_adj_nodes(w, a)
    00867                 {
    00868                         if (pG1->NodeType(w) == SHOCK_1)
    00869                                 num_child_a++;
    00870                 }
    00871                 forall_adj_nodes(s, b)
    00872                 {
    00873                         if (pG2->NodeType(s) == SHOCK_1)
    00874                                 num_child_b++;
    00875                 }
    00876
    00877                 if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    00878                         return 0;
    00879                 else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
    00880                         return 1;
    00881     }
    00882         else
    00883     {
    00884                 return 1;
    00885     }
    00886
    00887         double sim = abs(num_child_a - num_child_b);
    00888
    00889
    00890         // If this isn't true then it must go into death
    00891         // and it's outdeg should be zero!! I check for
    00892         // this above.
    00893         assert(num_child_a + num_child_b > 0);
    00894
    00895         // normalizing
    00896         sim = sim/(num_child_a +num_child_b);
    00897
    00898         return sim;
    00899 }
    00900
    00901 /*
    00902 This function tries to check how closely related both nodes
    00903 that were generated from rule 3 match. NOTE this function
    00904 assumes that the nodes being passed as parameters are of
    00905 type 3.
    00906 */
    00907
    00908 double SGDistMeasurer1::rule3(leda_node a, leda_node b)  const
    00909 {
    00910
    00911         leda_edge e = pG1->first_in_edge(a);
    00912         leda_edge f = pG2->first_in_edge(b);
    00913
    00914         int aParentType, bParentType;
    00915         int num_child_a = 0, num_child_b = 0;
    00916
    00917         int aType = pG1->NodeType(a);
    00918         int bType = pG2->NodeType(b);
    00919
    00920
    00921         // If a node has been removed such that it
    00922         //  no longer has any parent then we will
    00923         // assume that it's parent is the ROOT.
    00924         if(e == nil && f == nil)
    00925     {
    00926                 aParentType = ROOT;
    00927                 bParentType = ROOT;
    00928     }
    00929         else if (e == nil || f == nil)
    00930     {
    00931                 return 1;
    00932     }
    00933         else
    00934     {
    00935                 aParentType = pG1->NodeType(pG1->source(e));
    00936                 bParentType = pG2->NodeType(pG2->source(f));
    00937     }
    00938
    00939         leda_node w,s;
    00940
    00941         // rule 3a
    00942         if(aParentType == SHOCK_1 && bParentType == SHOCK_1
    00943                 && pG1->indeg(a) == 1 && pG2->indeg(b) == 1)
    00944     {
    00945
    00946                 forall_adj_nodes(w, a)
    00947                 {
    00948                         if (pG1->NodeType(w) == SHOCK_1)
    00949                                 num_child_a++;
    00950                 }
    00951                 forall_adj_nodes(s, b)
    00952                 {
    00953                         if (pG2->NodeType(s) == SHOCK_1)
    00954                                 num_child_b++;
    00955                 }
    00956
    00957                 // If both go to death...
    00958                 if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    00959                         return 0;
    00960                 else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
    00961                         return 1;
    00962     }
    00963         else if(aParentType == ROOT && bParentType == ROOT)
    00964     {
    00965                 // Rule 3 b
    00966                 forall_adj_nodes(w, a)
    00967                 {
    00968                         if (pG1->NodeType(w) == SHOCK_1)
    00969                                 num_child_a++;
    00970                 }
    00971                 forall_adj_nodes(s, b)
    00972                 {
    00973                         if (pG2->NodeType(s) == SHOCK_1)
    00974                                 num_child_b++;
    00975                 }
    00976
    00977                 if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    00978                         return 0;
    00979                 if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
    00980                         return 1;
    00981     }
    00982         else
    00983     {
    00984                 return 1;
    00985     }
    00986
    00987         double sim = abs(num_child_a - num_child_b);
    00988
    00989         // If this isn't true then it must go into death
    00990         // and it's outdeg should be zero!! I check for
    00991         // this above.
    00992         assert(num_child_a + num_child_b > 0);
    00993
    00994         // normalizing
    00995         sim = sim/(num_child_a + num_child_b);
    00996
    00997         return sim;
    00998 }
    00999
    01000 /*
    01001 This function tries to check how closely related both nodes
    01002 that were generated from rule 4 match.
    01003 */
    01004
    01005 double SGDistMeasurer1::rule4(leda_node a, leda_node b)  const
    01006 {
    01007         int aType = pG1->NodeType(a);
    01008         int bType = pG2->NodeType(b);
    01009         int num_child_a = 0, num_child_b = 0;
    01010         leda_node w,s;
    01011
    01012         if(aType == SHOCK_1 && bType == SHOCK_1)
    01013     {
    01014                 // If they go to death
    01015                 if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    01016                         return 0;
    01017                 else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
    01018                         return 1;
    01019
    01020                 forall_adj_nodes(w, a)
    01021                 {
    01022                         if (pG1->NodeType(w) == SHOCK_1)
    01023                                 num_child_a++;
    01024                 }
    01025                 forall_adj_nodes(s, b)
    01026                 {
    01027                         if (pG2->NodeType(s) == SHOCK_1)
    01028                                 num_child_b++;
    01029                 }
    01030     }
    01031         else
    01032     {
    01033                 return 1;
    01034     }
    01035
    01036         double sim = abs(num_child_a - num_child_b);
    01037
    01038         // If this isn't true then it must go into death
    01039         // and it's outdeg should be zero!! I check for
    01040         // this above.
    01041         assert(num_child_a + num_child_b > 0);
    01042
    01043
    01044         sim = sim/(num_child_a +num_child_b);
    01045
    01046         return sim;
    01047 }
    01048
    01049
    01050 /*
    01051 This function tries to check how closely related both nodes
    01052 that were generated from rule 5 match.
    01053 */
    01054
    01055 double SGDistMeasurer1::rule5_6(leda_node a, leda_node b)  const
    01056 {
    01057
    01058         int aType = pG1->NodeType(a);
    01059         int bType = pG2->NodeType(b);
    01060         int num_2child_a = 0;
    01061         int num_2child_b = 0;
    01062         int num_3child_a = 0;
    01063         int num_3child_b = 0;
    01064         leda_node w,s;
    01065
    01066         if(aType == SHOCK_1 && bType == SHOCK_1)
    01067     {
    01068                 forall_adj_nodes(w, a)
    01069                 {
    01070                         if (pG1->NodeType(w) == SHOCK_2)
    01071                                 num_2child_a++;
    01072                         else if (pG1->NodeType(w) == SHOCK_3)
    01073                                 num_3child_a++;
    01074                 }
    01075                 forall_adj_nodes(s, b)
    01076                 {
    01077                         if (pG2->NodeType(s) == SHOCK_2)
    01078                                 num_2child_b++;
    01079                         else if (pG2->NodeType(s) == SHOCK_3)
    01080                                 num_3child_b++;
    01081                 }
    01082
    01083                 /* This checks that a type 1 doesn't have more than
    01084                 1 type 2 child. */
    01085                 assert(num_2child_a <= 1 && num_2child_b <= 1);
    01086
    01087                 /* This checks that a type 1 doesn't have more than
    01088                 1 type 3 child. */
    01089                 assert(num_3child_a <= 1 && num_3child_b <= 1);
    01090
    01091                 if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    01092                         return 0;
    01093                 else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
    01094                         return 1;
    01095
    01096                 // rule 5
    01097                 if(num_2child_a == 1 && num_2child_b == 1){
    01098                         return 0;
    01099                 }
    01100                 else if(num_3child_a == 1 && num_3child_b == 1)
    01101                 {
    01102                         // rule 6a
    01103                         leda_node t,r;
    01104
    01105                         forall_adj_nodes(t, a)
    01106                         {
    01107                                 forall_adj_nodes(r, b)
    01108                                 {
    01109                                         if(pG1->indeg(t) == 2 && pG2->indeg(r) == 2)
    01110                                                 return 0;
    01111                                         else if(pG1->indeg(t) == 1 && pG2->indeg(r) == 1)
    01112                                                 return 0;
    01113                                         else{
    01114                                                 return 1;
    01115                                         }
    01116                                 }
    01117                         }
    01118                 }
    01119                 else{
    01120                         // Nodes have different types of children
    01121                         return 1;
    01122                 }
    01123     }
    01124         else
    01125         {
    01126                 return 1;
    01127         }
    01128 }
    01129
    01130 /*
    01131 This function is the control function for the second level of
    01132 matching.
    01133 */
    01134 double SGDistMeasurer1::NodeDistanceLevel2(leda_node a, leda_node b) const
    01135 {
    01136
    01137         int aType = pG1->NodeType(a);
    01138         int bType = pG2->NodeType(b);
    01139
    01140         if (aType == ROOT && bType == ROOT)
    01141                 return rule1(a,b);
    01142         else if(aType == ROOT || bType == ROOT)
    01143                 return 1;
    01144
    01145         if(pG1->NodeRole(a) == DUMMY && pG2->NodeRole(b) == DUMMY)
    01146     {
    01147                 return 0;
    01148     }
    01149         else if(pG1->NodeRole(a) == DUMMY || pG2->NodeRole(b) == DUMMY)
    01150     {
    01151                 return 1;
    01152     }
    01153
    01154         int aParentType = pG1->NodeType(pG1->source(pG1->first_in_edge(a)));
    01155         int bParentType = pG2->NodeType(pG2->source(pG2->first_in_edge(b)));
    01156
    01157         leda_node parent_a = pG1->source(pG1->first_in_edge(a));
    01158         leda_node parent_b = pG2->source(pG2->first_in_edge(b));
    01159
    01160         if(aType == SHOCK_1 && bType == SHOCK_1)
    01161     {
    01162                 double parent_score;
    01163
    01164                 // IF BOTH EMERGED FROM RULE 2
    01165                 if(aParentType == SHOCK_4 && bParentType == SHOCK_4)
    01166                 {
    01167                         // This will return a normalized number which
    01168                         // is the difference in siblings for each of
    01169                         // a and b
    01170
    01171                         parent_score = rule2(parent_a, parent_b);
    01172                 }
    01173                 // IF BOTH EMERGED FROM RULE 3
    01174                 else if(aParentType == SHOCK_3 && bParentType == SHOCK_3)
    01175                 {
    01176                         parent_score = rule3(parent_a, parent_b);
    01177                 }
    01178                 // IF BOTH EMERGED FROM RULE 4
    01179                 else if(aParentType == SHOCK_1 && bParentType == SHOCK_1)
    01180                 {
    01181                         parent_score = rule4(parent_a, parent_b);
    01182                 }
    01183                 else
    01184                 {
    01185                         // Nodes didn't emerge from the same production rule.
    01186                         parent_score = 1;
    01187                 }
    01188
    01189
    01190
    01191                 leda_node w,s;
    01192                 bool rule4a = false, rule4b = false;
    01193
    01194                 forall_adj_nodes(w, a)
    01195                 {
    01196                         if (pG1->NodeType(w) == SHOCK_1){
    01197                                 rule4a = true;
    01198                                 break;
    01199                         }
    01200                 }
    01201
    01202                 forall_adj_nodes(s, b)
    01203                 {
    01204                         if (pG2->NodeType(s) == SHOCK_1){
    01205                                 rule4b = true;
    01206                                 break;
    01207                         }
    01208                 }
    01209
    01210                 // DESCENDANTS SIMILARITY
    01211
    01212                 double child_score;
    01213
    01214                 if(rule4a && rule4b)
    01215                 {
    01216                         child_score = rule4(a, b);
    01217                 }
    01218                 else if(!rule4a && !rule4b)
    01219                 {
    01220                         child_score = rule5_6(a, b);
    01221                 }
    01222                 else if(pG1->indeg(a) == 0 && pG2->indeg(b) == 0){
    01223                         // If both go to death
    01224                         child_score = 0;
    01225                 }
    01226                 else
    01227                 {
    01228                         child_score = 1;
    01229                 }
    01230                 return (parent_score + child_score)/2;
    01231     }
    01232         else if(aType == SHOCK_2 && bType == SHOCK_2)
    01233     {
    01234                 // ANCESTOR SIMILARITY
    01235                 if(aParentType == SHOCK_1 && bParentType == SHOCK_1 && pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
    01236                         return 0;
    01237                 else
    01238                 {
    01239                         cerr << "A type 2 shock doesn't go to death or have type 1 parents!!" << endl;
    01240                         assert(false);
    01241                         return -1;
    01242                 }
    01243     }
    01244         else if(aType == SHOCK_3 && bType == SHOCK_3)
    01245     {
    01246                 double child_score = rule3(a, b);
    01247                 double ancestor_score;
    01248
    01249                 if(aParentType == ROOT && bParentType == ROOT)
    01250                         ancestor_score = 0;
    01251                 else if(pG1->indeg(a) == 2 && pG2->indeg(b) == 2)
    01252                         ancestor_score = 0;
    01253                 else if(pG1->indeg(a) == 1 && pG2->indeg(b) == 1)
    01254                 {
    01255                         ancestor_score = 0;
    01256                 }
    01257                 else
    01258                         ancestor_score = 1;
    01259
    01260                 return (child_score + ancestor_score)/2;
    01261     }
    01262         else if(aType == SHOCK_4 && bType == SHOCK_4)
    01263     {
    01264                 // Descendants
    01265                 double child_score = rule2(a, b);
    01266                 // parents
    01267                 double parent_score = rule1(parent_a,parent_b);
    01268
    01269                 return (child_score + parent_score)/2;
    01270     }
    01271         else{
    01272                 return 1;
    01273         }
    01274 }
    01275
    01277 // SGDistMeasurer2
    01278
    01279 #define MAX_DIST       100
    01280 #define DIFF_ROLE_DIST 100
    01281 #define DIFF_TYPE_DIST 1
    01282
    01283 /*void DbgNodeSimilarity(const ShockGraph& sg1, leda_node u, const ShockGraph& sg2, leda_node v,
    01284 double scale, double dr, double dv,
    01285 double branchSim, double tsvSim, double nodeSim)
    01286 {
    01287 static bool bFirstTime = true;
    01288 static fstream dbg;
    01289
    01290  if (bFirstTime)
    01291  {
    01292  bFirstTime = false;
    01293  dbg.open("nodesim.log", ios::trunc | ios::out);
    01294  }
    01295
    01296   dbg << sg1.GetDAGLbl() << " vs " << sg2.GetDAGLbl() << endl
    01297   << sg1.GetNodeLbl(u) << " vs " << sg2.GetNodeLbl(v) << endl
    01298   << "role " << sg1.NodeRole(u) << " vs role " << sg2.NodeRole(v) << endl
    01299   << "scale = " << scale << ", dr = " << dr << ", dv = " << dv << endl
    01300   << "Branch similarity = " << branchSim << endl
    01301   << "TSV similarity = " << tsvSim << endl
    01302   << "Node similarity = " << nodeSim << endl
    01303   << endl << flush;
    01304 }*/
    01305
    01306 /*bool SGDistMeasurer1::AreNodesRelated() const
    01307 {
    01308 const ShockGraph& sgFrom = dynamic_cast<const ShockGraph&>(from);
    01309
    01310  NODE_ROLE ru = NodeRole(u), rv = sgFrom.NodeRole(v);
    01311
    01312   if (ru == UNK_ROLE || rv == UNK_ROLE)
    01313   return NodeType(u) == sgFrom.NodeType(v);
    01314   else if (ru == DUMMY || rv == DUMMY)
    01315   return false;
    01316   else
    01317   return ru == rv;
    01318 }*/
    01319
    01320 double SGDistMeasurer2::CompDistance() const
    01321 {
    01322         ASSERT(false);
    01323         return 0;
    01324 }
    01325
    01326 /*double SGDistMeasurer2::CompDistance()
    01327 {
    01328 double da = 0, dv = 0, scale = 0, branchSim;
    01329 const ShockBranch& b1 = pG1->GetSGNode(v1)->m_shocks;
    01330 const ShockBranch& b2 = pG2->GetSGNode(v2)->m_shocks;
    01331
    01332  int uType = pG1->NodeType(v1);
    01333
    01334   if (pG1->NodeRole(v1) != BIRTH)
    01335   {
    01336   if (uType == SHOCK_4)
    01337   {
    01338   branchSim = 0;
    01339   }
    01340   else
    01341   {
    01342   //scale = ComputeScale(u, sgFrom, v);
    01343
    01344           da = b1.CompareAreas(b2, 1); //use scale!!!
    01345           dv = b1.CompareVelocities(b2, 1); //use scale!!!
    01346
    01347            const double w = .7;
    01348            branchSim = w * da + (1 - w) * dv;
    01349            }
    01350            }
    01351            else
    01352            {
    01353            double n1, n2;
    01354
    01355                 n1 = pG1->NeighborsCount(v1, SHOCK_3);
    01356                 n2 = pG2->NeighborsCount(v2, SHOCK_3);
    01357
    01358                  branchSim = SIMILARITY(n1, n2, 0.5);
    01359
    01360                   n1 = pG1->NeighborsCount(v1, SHOCK_4);
    01361                   n2 = pG2->NeighborsCount(v2, SHOCK_4);
    01362
    01363                    branchSim *= SIMILARITY(n1, n2, 0.5);
    01364                    }
    01365
    01366                         double tsvSim = pG1->NodeTSVSimilarity(v1, *pG2, v2);
    01367                         double nodeSim = (1 - s_dTSVSimWeight) * branchSim + s_dTSVSimWeight * tsvSim;
    01368
    01369                          return nodeSim;
    01370                          }
    01371
    01372                           int SGDistMeasurer2::NeighborsCount(leda_node v, int nNodeType) const
    01373                           {
    01374                           int nCount = 0;
    01375                           leda_edge e;
    01376                           leda_node w;
    01377
    01378                            forall_in_edges(e, v)
    01379                            if (pG1->NodeType(pG1->source(e)) == nNodeType)
    01380                            nCount++;
    01381
    01382                                 forall_adj_nodes(w, v)
    01383                                 if (NodeType(w) == nNodeType)
    01384                                 nCount++;
    01385
    01386                                  return nCount;
    01387 }*/
    01388
    01389 // Computes the difference in scale between two branches by using their parents info.
    01390 /*double SGDistMeasurer2::ComputeScale(leda_node u, const ShockGraph& from, leda_node v) const
    01391 {
    01392 double r1, r2;
    01393 leda_node uparent = source(first_in_edge(u));
    01394 leda_node vparent = from.source(from.first_in_edge(v));
    01395
    01396  // Root node shouldn't be analysed!
    01397  ASSERT(uparent != nil && vparent != nil);
    01398
    01399   if (NodeType(uparent) == ROOT)
    01400   r1 = GetSGNode(u)->m_shocks.AvgRadius();
    01401   else
    01402   r1 = GetSGNode(uparent)->m_shocks.AvgRadius();
    01403
    01404           if (from.NodeType(vparent) == ROOT)
    01405           r2 = from.GetSGNode(v)->m_shocks.AvgRadius();
    01406           else
    01407           r2 = from.GetSGNode(vparent)->m_shocks.AvgRadius();
    01408
    01409            return  (r1 == 0 || r2 == 0) ? 1:r1 / r2;
    01410 }*/
    01411
    01412
    01413
    01414
    01415