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

# ShockGraph.cpp

Go to the documentation of this file.

    00001
    00046 #include "DAG.h"
    00047 #include "ShockGraph.h"
    00048 #include "HelperFunctions.h"
    00049 #include "Exceptions.h"
    00050 #include <float.h>
    00051 #include <stdio.h>
    00052 #include <ctype.h>
    00053
    00054 #include "SGDistMeasurer.h"
    00055
    00056 #include <LEDA/graphwin.h>
    00057 #include <LEDA/graph.h>
    00058 #include <LEDA/node_array.h>
    00059
    00060 #include "VisualDAG.h"
    00061 #include "ApproxPoly.h"
    00062
    00063 #define NSAMPLES         5
    00064 #define VDIM             4
    00065 #define DIMxSAMPLES  (VDIM * NSAMPLES)
    00066
    00067 #define SA_GROWING_FACTOR 50
    00068
    00069 #ifdef WIN32
    00070 #define FILE_SEP '\\'
    00071 #else
    00072 #define FILE_SEP '/'
    00073 #endif
    00074
    00075 struct ER // Pair of edge and radius values
    00076 {
    00077         leda_edge e;
    00078         double r; // avg radius of target(e)
    00079
    00080         ER() { e = nil; r = 0.0; }
    00081         ER(leda_edge edge, const double& radius) { e = edge; r = radius; }
    00082
    00092         static int CompareEdgesDescending(const void *elem1, const void *elem2 )
    00093         {
    00094                 ER* p1 = (ER*) elem1;
    00095                 ER* p2 = (ER*) elem2;
    00096
    00097                 return (p1->r < p2->r) ? 1:(p1->r > p2->r) ? -1:0;
    00098         }
    00099 };
    00100
    00101 void SmoothFluxPoints(sg::FluxPointList& fpl, const double& dSigma, const double& dSigmaRange);
    00102 void GaussianSmooth(double *input, double *output, int length, const double& sigma, const double& sigma_range);
    00103
    00104 const char* ShockGraph::GetNextIdx()
    00105 {
    00106         static char szNum[100];
    00107
    00108         sprintf(szNum, "%d", ++nLastIndexUsed);
    00109         return szNum;
    00110 }
    00111
    00112 ShockGraph::ShockGraph()
    00113 {
    00114         area = 0;
    00115         perimeter = 0;
    00116         xmin = ymin = 0;
    00117         xmax = ymax = 0;
    00118         nLastIndexUsed = 0;
    00119         m_pSkeleton = NULL;
    00120
    00121         ASSERT(!pDistMeasurer);
    00122
    00123         if (s_matchParams.nNodeDistFunc == 1)
    00124                 pDistMeasurer = new SGDistMeasurer1;
    00125         else if (s_matchParams.nNodeDistFunc == 2)
    00126                 pDistMeasurer = new SGDistMeasurer2;
    00127         else if (s_matchParams.nNodeDistFunc == 3)
    00128                 pDistMeasurer = new SGDistMeasurer3;
    00129         else
    00130                 ASSERT(false);
    00131 }
    00132
    00133 ShockGraph::~ShockGraph()
    00134 {
    00135         delete pDistMeasurer;
    00136         //WARNING(m_pSkeleton, "Skeletton destructor isn't working at this point. Memory Leak!");
    00137         //delete m_pSkeleton;
    00138 }
    00139
    00140 void ShockGraph::Clear()
    00141 {
    00142         DAG::Clear();
    00143
    00144         area = 0.0;
    00145         perimeter = 0.0;
    00146         xmin = ymin = 0;
    00147         xmax = ymax = 0;
    00148         nLastIndexUsed = 0;
    00149
    00150         delete m_pSkeleton;
    00151         m_pSkeleton = NULL;
    00152 }
    00153
    00154 void ShockGraph::ComputeObjectViewName()
    00155 {
    00156         int i, j, k, lastFileSep;
    00157         String viewName;
    00158         const String path = GetDAGLbl();
    00159         const int nLen = path.Len();
    00160
    00161         strObjectName.Clear();
    00162
    00163         // Find first char of the last directory name in the path
    00164
    00165         for (i = nLen - 1; i >= 0 && path[i] != FILE_SEP; i--);
    00166
    00167         lastFileSep = i;
    00168
    00169         if (i > 0) // we have a string with format '/ObjName/prefixNNN.ext'
    00170         {
    00171                 // Search for next '/'
    00172                 for (j = i - 1; j >= 0 && path[j] != FILE_SEP; j--);
    00173
    00174                 // j point to the second '/'
    00175                 strObjectName.ReSize(i - j);
    00176
    00177                 for (j++, k = 0;  j < i; j++, k++)
    00178                         strObjectName[k] = path[j];
    00179
    00180                 strObjectName[k] = '\0';
    00181
    00182                 // Find beginning of view number
    00183                 for (i++; i < nLen && !isdigit(path[i]); i++);
    00184         }
    00185         else // we have a string with format 'ObjNameNNN.ext'
    00186         {
    00187                 // Find beginning of view number
    00188                 for (i = 0; i < nLen && !isdigit(path[i]); i++);
    00189
    00190                 // Copy obj name
    00191
    00192                 strObjectName.ReSize(i + 1);
    00193
    00194                 for (j = 0; j < i; j++)
    00195                         strObjectName[j] = path[j];
    00196
    00197                 strObjectName[j] = '\0';
    00198         }
    00199
    00200         // Find view number
    00201         viewName.ReSize(nLen - i + 1);
    00202
    00203         for (k = 0;  i <= nLen; i++, k++)
    00204                 viewName[k] = path[i];
    00205
    00206         nViewNumber = (viewName[0] != '\0') ? atoi(viewName):1;
    00207 }
    00208
    00210 void ShockGraph::ComputeDerivedValues()
    00211 {
    00212         DAG::ComputeDerivedValues();
    00213         ComputeObjectViewName();
    00214 }
    00215
    00217 istream& ShockGraph::Read(istream& is, bool bOnlyDataForMatching /*= false*/)
    00218 {
    00219         DAG::Read(is, bOnlyDataForMatching);
    00220
    00221         is.read((char*) &area, sizeof(area));
    00222         is.read((char*) &perimeter, sizeof(perimeter));
    00223         is.read((char*) &xmin, sizeof(xmin));
    00224         is.read((char*) &xmax, sizeof(xmax));
    00225         is.read((char*) &ymin, sizeof(ymin));
    00226         is.read((char*) &ymax, sizeof(ymax));
    00227
    00228         if (!bOnlyDataForMatching)
    00229         {
    00230                 // pSkeleton works as a boolean variable here
    00231                 is.read((char*) &m_pSkeleton, sizeof(m_pSkeleton));
    00232                 if (m_pSkeleton)
    00233                 {
    00234                 m_pSkeleton = new McGillSkeleton();
    00235                 m_pSkeleton->Read(is);
    00236                 }
    00237         }
    00238
    00239         return is;
    00240 }
    00241
    00243 ostream& ShockGraph::Write(ostream& os) const
    00244 {
    00245         DAG::Write(os);
    00246
    00247         os.write((char*) &area, sizeof(area));
    00248         os.write((char*) &perimeter, sizeof(perimeter));
    00249         os.write((char*) &xmin, sizeof(xmin));
    00250         os.write((char*) &xmax, sizeof(xmax));
    00251         os.write((char*) &ymin, sizeof(ymin));
    00252         os.write((char*) &ymax, sizeof(ymax));
    00253
    00254         // pSkeleton works as a boolean variable here
    00255         os.write((char*) &m_pSkeleton, sizeof(m_pSkeleton));
    00256         if (m_pSkeleton)
    00257           m_pSkeleton->Write(os);
    00258
    00259         return os;
    00260 }
    00261
    00263 void ShockGraph::Print(ostream& os) const
    00264 {
    00265         DAG::Print(os);
    00266
    00267         os << "\nShockGraph specific:\n";
    00268
    00269         PRINT_OPEN(os, area);
    00270         PRINT(os, perimeter);
    00271         PRINT(os, xmin);
    00272         PRINT(os, xmax);
    00273         PRINT(os, ymin);
    00274         PRINT_CLOSE(os, ymax);
    00275 }
    00276
    00278 DAG& ShockGraph::operator=(const DAG& rhs)
    00279 {
    00280         DAG::operator=(rhs);
    00281
    00282         const ShockGraph* pRhs = dynamic_cast<const ShockGraph*>(&rhs);
    00283
    00284         if (!pRhs)
    00285                 THROW_EXCEPTION("Invalid pointer type.");
    00286
    00287         area      = pRhs->area;
    00288         perimeter = pRhs->perimeter;
    00289
    00290         xmin = pRhs->xmin;
    00291         ymin = pRhs->ymin;
    00292         xmax = pRhs->xmax;
    00293         ymax = pRhs->ymax;
    00294
    00295         return *this;
    00296 }
    00297
    00301 double ShockGraph::NodeDistance(leda_node u, const DAG& from, leda_node v) const
    00302 {
    00303         ASSERT(pDistMeasurer);
    00304
    00305         return pDistMeasurer->GetDistance(*this, u, from, v);
    00306 }
    00307
    00312 bool ShockGraph::AreNodesRelated(leda_node u, const DAG& from, leda_node v) const
    00313 {
    00314         return true;
    00315 }
    00316
    00318 DAG* ShockGraph::CreateObject() const
    00319 {
    00320         return new ShockGraph();
    00321 }
    00322
    00324 DAGNodePtr ShockGraph::CreateNodeObject(NODE_LABEL lbl) const
    00325 {
    00326         return new SGNode(lbl, TEST_NODE);
    00327 }
    00328
    00333 DAGNodePtr ShockGraph::ReadNode(istream& is) const
    00334 {
    00335         DAGNodePtr node(new SGNode);
    00336
    00337         node->Read(is);
    00338
    00339         return node;
    00340 }
    00341
    00343 String ShockGraph::ClassName() const
    00344 {
    00345         return "ShockGraph";
    00346 }
    00347
    00348 NODE_LABEL ShockGraph::GetNodeLbl(leda_node v) const
    00349 {
    00350         char szLbl[100];
    00351         const char* szNodeLbl = (const char*)DAG::GetNodeLbl(v);
    00352
    00353         sprintf(szLbl, "%s:%d", szNodeLbl, NodeType(v) % 100);
    00354
    00355         return szLbl;
    00356 }
    00357
    00358 /*************************************************************************************************************************/
    00359 /* Interface with shock graph code */
    00360
    00361 bool ShockGraph::ComputeFromPPMFile(const char* szFileName,
    00362                                                                         double cutoff, double sigma, double range)
    00363 {
    00364 #ifdef COMPILE_OLD_SG_CODE
    00365         MGShockGraph sg;
    00366
    00367         if (!sg.ComputeFromPPMFile(szFileName, cutoff, sigma, range) || !Translate(sg))
    00368                 return false;
    00369
    00370         SetDAGLbl(szFileName);
    00371
    00372         ComputeDerivedValues();
    00373 #endif //COMPILE_OLD_SG_CODE
    00374
    00375         return true;
    00376 }
    00377
    00378 #ifdef COMPILE_OLD_SG_CODE
    00379
    00380 bool ShockGraph::ComputeFromMemory(Pixmap& image, double cutoff, double sigma, double range)
    00381 {
    00382         MGShockGraph sg;
    00383
    00384         if (!sg.ComputeFromMemory(image, cutoff, sigma, range) || !Translate(sg))
    00385                 return false;
    00386
    00387         SetDAGLbl(image.name);
    00388
    00389         ComputeDerivedValues();
    00390
    00391         return true;
    00392 }
    00393
    00395 bool ShockGraph::Translate(const MGShockGraph& sg)
    00396 {
    00397         leda_d_array<long, leda_node> map;
    00398         int i, j, ptCount;
    00399         SGNode* pNode;
    00400         double seg1, seg2;
    00401
    00402         Clear();
    00403
    00404         for (i = 0; i < sg.GetNodeCount(); i++)
    00405         {
    00406                 const Vertex* v = sg.GetNode(i);
    00407                 MGShockGraph::SGSegment branch = sg.GetSegment(v);
    00408
    00409                 pNode = new SGNode(branch.GetId(), branch.GetType());
    00410
    00411                 branch.GetContourLength(seg1, seg2);
    00412                 pNode->SetContourLength(seg1, seg2);
    00413
    00414                 ptCount = branch.GetPointCount();
    00415                 pNode->m_shocks.ReSize((ptCount > 0) ? ptCount:0);
    00416
    00417                 for (j = 0; j < ptCount; j++)
    00418                 {
    00419                         const MGShockGraph::SGPoint& pt = branch[j];
    00420                         ShockInfo& si = pNode->m_shocks[j];
    00421
    00422                         si.xcoord = pt.x;
    00423                         si.ycoord = pt.y;
    00424                         si.radius = pt.r;
    00425                         si.speed = pt.s;
    00426                         si.dr_ds = pt.dr_ds;
    00427                         si.color = pt.color;
    00428
    00429                         // Set bounding rectangle for drawing purposes
    00430                         if (pt.x < xmin) xmin = pt.x;
    00431                         if (pt.x > xmax) xmax = pt.x;
    00432                         if (pt.y < ymin) ymin = pt.y;
    00433                         if (pt.y > ymax) ymax = pt.y;
    00434                 }
    00435
    00436                 map[(long)v] = NewNode(pNode);
    00437         }
    00438
    00439         for (i = 0; i < sg.GetEdgeCount(); i++)
    00440         {
    00441                 const Edge* e = sg.GetEdge(i);
    00442
    00443                 ASSERT(map.defined((long)e->vFrom) && map.defined((long)e->vTo));
    00444
    00445                 // From an To nodes have to be swapped because we consider the
    00446                 // the directions are from the root to the leaves.
    00447                 NewEdge(map[(long)e->vTo], map[(long)e->vFrom]);
    00448         }
    00449
    00450         return true;
    00451 }
    00452
    00453 #endif //COMPILE_OLD_SG_CODE
    00454
    00455
    00456 // test code from McGill
    00457 char* get_arr(int *rows, int *cols, const char* shapefile)
    00458 {
    00459         char *arr;
    00460         int k, i, j;
    00461
    00462         pixval maxval;
    00463         pixel **pixels;
    00464
    00465         FILE* ppmshape = fopen(shapefile, "r");
    00466
    00467         pixels = ppm_readppm(ppmshape, cols, rows, &maxval);
    00468         fclose(ppmshape);
    00469
    00470         arr = (char*)malloc((*cols)*(*rows));
    00471
    00472         for(i = 0, k = 0; i < *rows; i++)
    00473                 for(j = 0; j < *cols; j++, k++)
    00474                 {
    00475                         char r, g, b;
    00476
    00477                         r = (0xff * PPM_GETR(pixels[i][j])) / maxval;
    00478                         g = (0xff * PPM_GETR(pixels[i][j])) / maxval;
    00479                         b = (0xff * PPM_GETR(pixels[i][j])) / maxval;
    00480
    00481                         if(r==0 && b==0 && g==0)
    00482                                 arr[k] = 0;
    00483                         else
    00484                                 arr[k] = 0xff;
    00485                 }
    00486
    00487                 ppm_freearray(pixels, *rows);
    00488
    00489                 return arr;
    00490 };
    00491
    00497 void ShockGraph::LabelEndPoints(SGNode* pNode, sg::FluxPointList& fpl, sg::DDSNode* n1, sg::DDSNode* n2)
    00498 {
    00499         sg::DDSEdgeVect *e0, *eN;
    00500
    00501         if (n1->getPoint() == fpl[0].p)
    00502         {
    00503                 ASSERT(n2->getPoint() == fpl[fpl.size() - 1].p);
    00504
    00505                 e0 = &n1->getEdges();
    00506                 eN = &n2->getEdges();
    00507         }
    00508         else
    00509         {
    00510                 ASSERT(n1->getPoint() == fpl[fpl.size() - 1].p && n2->getPoint() == fpl[0].p);
    00511
    00512                 e0 = &n2->getEdges();
    00513                 eN = &n1->getEdges();
    00514         }
    00515
    00516         pNode->m_nEndPt0 = (e0->size() == 1) ? TERMINAL_POINT:JOINT_POINT;
    00517         pNode->m_nEndPtN = (eN->size() == 1) ? TERMINAL_POINT:JOINT_POINT;
    00518 }
    00519
    00520 SGNode* ShockGraph::CopyNodeInfo(sg::DDSEdge* ddsEdge, leda_node u, SGNode* pNode)
    00521 {
    00522         sg::FluxPointList& fpl = ddsEdge->getFluxPoints();
    00523         int num_pts = fpl.size();
    00524         ASSERT(num_pts > 0);
    00525         double dx, dy;
    00526
    00527         //SmoothFluxPoints(fpl, compParams.dSigma, compParams.dSigmaRange);
    00528
    00529         pNode->m_shocks.ReSize(num_pts);
    00530
    00531         for (int j = 0; j < num_pts; j++)
    00532         {
    00533                 const sg::FluxPoint& pt = fpl[j];
    00534                 ShockInfo& si = pNode->m_shocks[j];
    00535
    00536                 si.xcoord = pt.p.x;
    00537                 si.ycoord = pt.p.y;
    00538                 si.radius = fabs(pt.dist);
    00539                 si.color = 0; //pt.col; // not used anymore. see type instead
    00540
    00541                 // Some NaN have been found in a value for radius
    00542                 if (isnan(si.radius) && !finite(si.radius))
    00543                 {
    00544                         cerr << "\nShock point has a NaN as radius." << endl;
    00545                         // Try to recover from this
    00546                         si.radius = (j > 0) ? pNode->m_shocks[j - 1].radius:0.5;
    00547                 }
    00548
    00549                 if (j > 0)
    00550                 {
    00551                         ShockInfo& prevSi = pNode->m_shocks[j - 1];
    00552
    00553                         dx = si.xcoord - prevSi.xcoord;
    00554                         dy = si.ycoord - prevSi.ycoord;
    00555
    00556                         si.dr = si.radius - prevSi.radius;
    00557                         si.speed = sqrt(dx * dx + dy * dy);
    00558                         si.dr_ds = si.dr / si.speed;
    00559
    00560                         // dir is either -1, +1 or 0
    00561                         si.dir = (fabs(si.dr_ds) <= compParams.dMinSlope) ? 0:int(si.dr / fabs(si.dr));
    00562                 }
    00563                 else
    00564                 {
    00565                         si.dr = 0;
    00566                         si.speed = 0;
    00567                         si.dr_ds = 0;
    00568                         si.dir = 0;
    00569                 }
    00570         }
    00571
    00572         LabelEndPoints(pNode, fpl, ddsEdge->getN1(), ddsEdge->getN2());
    00573
    00574         return GroupShockPoints(u, pNode);
    00575 }
    00576
    00581 SGNode* ShockGraph::GroupShockPoints(leda_node u, SGNode* pNode)
    00582 {
    00583         ShockBranch& pts = pNode->m_shocks;
    00584         //int nSize = pts.GetSize();
    00585         SGNode* pNewNode;
    00586         SGNode* pReturnNode = NULL;
    00587
    00588         int i, from;
    00589         double m, nextSlope;
    00590
    00591         int d0, dN;     // if d0 or dN == 1, means their end points were chopped off
    00592         POINTS data = pNode->GetVelocityRadiusArray(d0, dN);
    00593
    00594         double x0 = data[0].x;
    00595
    00596         // Time to fit the data with a few lines
    00597         ApproxPoly poly(data.GetSize() / compParams.dMinError, compParams.dMinSlope, 10.0);
    00598
    00599         if (DAG::IsDbgMode())
    00600         {
    00601                 cout << "disp('Fitting node: " << pNode->GetNodeLbl() << " in " << GetDAGLbl() << "');\n";
    00602                 poly.SetDbgMode(true); // outputs matlab data in dbg mode
    00603         }
    00604
    00605         poly.Fit(data);
    00606
    00607         int nSize = poly.m_knots.GetSize();
    00608         ASSERT(nSize >= 1);
    00609         int dir, ptIdx;
    00610         SmartArray<SEGMENT> approx_segs;
    00611
    00612         for (from = 0, i = 0; i < nSize; i++)
    00613         {
    00614                 approx_segs.AddTail(poly.m_knots[i].seg);
    00615
    00616                 dir = poly.m_knots[i].dir;
    00617                 ptIdx = poly.m_knots[i].nIndex + d0;
    00618
    00619                 if (i < nSize - 1)
    00620                 {
    00621                         // If the segment is too small, skip it.
    00622                         if (ptIdx - from < 2)
    00623                                 continue;
    00624
    00625                         // If there is a significant change in acceleration,
    00626                         // we create a new node, even though the direction is the same.
    00627                         if (dir == poly.m_knots[i + 1].dir)
    00628                         {
    00629                                 if (dir == 0)
    00630                                         continue;
    00631
    00632                                 //poly.CompAcuteAngle(i, i + 1) < compParams.dMaxAccelChg)
    00633                                 double m0 = fabs(poly.m_knots[i].seg.m);
    00634                                 double m1 = fabs(poly.m_knots[i + 1].seg.m);
    00635                                 double mr = fabs(m1 - m0) / MAX(m0, m1);
    00636
    00637                                 if (DAG::IsDbgMode())
    00638                 {
    00639                         cout << "disp('Acc change: " << mr << ", Max:" << compParams.dMaxAccelChg << "');\n";
    00640                         poly.SetDbgMode(true); // outputs matlab data in dbg mode
    00641                 }
    00642
    00643                                 if (mr < compParams.dMaxAccelChg)
    00644                                         continue;
    00645                         }
    00646
    00647                         pNewNode = SplitNode(u, pNode, ptIdx - from, approx_segs[approx_segs.GetSize() - 1], poly.m_knots[i + 1].seg);
    00648                 }
    00649                 else
    00650                         pNewNode = pNode;
    00651
    00652                 // Record the first node of the chain
    00653                 if (pReturnNode == NULL)
    00654                         pReturnNode = pNewNode;
    00655
    00656                 pNewNode->SetNodeType(dir == 0 ? SHOCK_3:SHOCK_1);
    00657                 pNewNode->m_shocks.SetDir(dir);
    00658
    00659                 // Update the fitting segments. Note that if <from> != 0, the x coord of points in the
    00660                 // branch will be shifted. SetSegments will correct this gived the displacement.
    00661                 pNewNode->m_shocks.SetSegments(approx_segs, from == 0 ? 0.0:approx_segs[0].p0.x - x0);
    00662
    00663                 pNewNode->ComputeDerivedValues();
    00664
    00665                 from = ptIdx;
    00666                 approx_segs.Clear();
    00667         }
    00668
    00669         /*for (int i = 1; i < nSize; i++)
    00670         {
    00671                 // Do we need to split the node?
    00672                 if (pts[i].type != pts[0].type)
    00673                 {
    00674                         pNewNode = SplitNode(u, pNode, i - 1, pts[0].type);
    00675                         pNewNode->ComputeDerivedValues();
    00676                         GroupShockPoints(u, pNode);
    00677                         break;
    00678                 }
    00679         }
    00680
    00681         pNode->ComputeDerivedValues();*/
    00682
    00683         return pReturnNode;
    00684 }
    00685
    00686 /*
    00687 @brief Splits in two the shock points in a node. It returns the new node created that
    00688 contains the first points. It also links the two nodes.
    00689 */
    00690 SGNode* ShockGraph::SplitNode(leda_node v, SGNode* pNode, int nEnd, SEGMENT& leftSeg, SEGMENT& rightSeg)
    00691 {
    00692         SGNode* pNewNode = new SGNode(GetNextIdx(), 0);
    00693         ShockBranch& newPts = pNewNode->m_shocks;
    00694         ShockBranch& pts = pNode->m_shocks;
    00695
    00696         // We must repeat the endpoint where the nodes connect. nEnd is inclusive
    00697         pNewNode->m_shocks = pts.Left(nEnd);
    00698
    00699         // Now remove the left part. nEnd is inclusive
    00700         pts = pts.Right(nEnd);
    00701
    00702         // Create new node
    00703         leda_node u = NewNode(pNewNode);
    00704
    00705         // Update the m_nEndPt0 and m_nEndPtN member variables of both nodes
    00706         pNewNode->m_nEndPt0 = pNode->m_nEndPt0;
    00707         pNewNode->m_nEndPtN = SPLIT_POINT;
    00708         pNode->m_nEndPt0 = SPLIT_POINT;
    00709
    00710         // Update the Y value (radius) at the joint pt, to avoid outliers in this crusial
    00711         // function. i.e. type 2 and 4 will depend on accurate values here
    00712         double jointRadius = (leftSeg.Y1() + rightSeg.Y0()) / 2.0;
    00713         pNewNode->m_shocks[pNewNode->m_shocks.GetSize() - 1].radius = jointRadius;
    00714         pts[0].radius = jointRadius;
    00715         leftSeg.p1.y = jointRadius;
    00716         rightSeg.p0.y = jointRadius;
    00717
    00718         // If 'v' has a parent node, move it to be a parent of 'u' instead
    00719         leda_edge e = first_in_edge(v);
    00720
    00721         if (e != nil)
    00722     {
    00723                 NewEdge(source(e), u);
    00724                 ASSERT(in_succ(e) == nil); // There shoukd be only one edge
    00725                 DelEdge(e);
    00726     }
    00727
    00728         // Link u to v
    00729         NewEdge(u, v);
    00730
    00731         return pNewNode;
    00732 }
    00733
    00737 void ShockGraph::ConnectNodes(leda_node u, SGNode* pNode, sg::DDSEdge* e, sg::DDSNode* n)
    00738 {
    00739         using namespace sg;
    00740
    00741         // Then, keep creating and connecting nodes
    00742         DDSEdgeVect& edges = n->getEdges();
    00743         DDSEdgeVect::iterator I;
    00744         DDSEdge* ee;
    00745         SGNode* pNewNode;
    00746         leda_node v;
    00747
    00748         for(I = edges.begin(); I != edges.end(); I++)
    00749         {
    00750                 ee = *I;
    00751
    00752                 if (e != ee)
    00753                 {
    00754                         pNewNode = new SGNode(GetNextIdx(), SHOCK_base);
    00755                         v = NewNode(pNewNode);
    00756                         NewEdge(u, v);
    00757                         //CopyNodeInfo(ee->getFluxPoints(), v, pNewNode);
    00758                         CopyNodeInfo(ee, v, pNewNode);
    00759                         ConnectNodes(v, pNewNode, ee, ee->getN1() != n ? ee->getN1():ee->getN2());
    00760                 }
    00761         }
    00762 }
    00763
    00773 void ShockGraph::RevertEdge(leda_edge pxEdge)
    00774 {
    00775         leda_node parent = source(pxEdge);
    00776         leda_node x = target(pxEdge);
    00777         leda_edge e;
    00778
    00779         rev_edge(pxEdge);
    00780
    00781         // Move all children of x's parent (that are also connected to x)
    00782         // to the left side of x
    00783         const ShockInfo& inPt = GetJointPoint(parent, x);
    00784
    00785         forall_adj_edges(e, parent)
    00786                 if (GetJointPoint(parent, target(e)) == inPt)
    00787                         move_edge(e, x, target(e));
    00788 }
    00789
    00794 void ShockGraph::Relabel3As2or4(leda_node v)
    00795 {
    00796         ASSERT(NodeType(v) == SHOCK_3);
    00797
    00798         int parentType = NodeType(GetFirstParent(v));
    00799
    00800         if (parentType == ROOT)
    00801                 UnsafeGetSGNode(v)->SetNodeType(SHOCK_4);
    00802         else if(outdeg(v) == 0 && indeg(v) == 2)
    00803                 UnsafeGetSGNode(v)->SetNodeType(SHOCK_2);
    00804
    00805         // If we are interested in spliting 3's into left and
    00806         // right nodes, we can do it here.
    00807         if (compParams.nSlipt3s != 0 && parentType == ROOT)
    00808         {
    00809                 const SGNode& x = *GetSGNode(v);
    00810
    00811                 const ShockInfo& leftPt = x[0];
    00812                 const ShockInfo& rightPt = x[x.GetSize() - 1];
    00813
    00814                 leda_node leftChild = NewNode(new SGNode(x));
    00815                 leda_node rightChild = NewNode(new SGNode(x));
    00816                 leda_edge e;
    00817
    00818                 forall_adj_edges(e, v)
    00819                 {
    00820                         if (GetJointPoint(v, target(e)) == leftPt)
    00821                                 move_edge(e, leftChild, target(e));
    00822                         else
    00823                         {
    00824                                 ASSERT(GetJointPoint(v, target(e)) == rightPt);
    00825                                 move_edge(e, rightChild, target(e));
    00826                         }
    00827                 }
    00828
    00829                 // Link x to its left and right sides
    00830                 NewEdge(v, rightChild);
    00831                 NewEdge(v, leftChild);
    00832         }
    00833 }
    00834
    00843 void ShockGraph::Insert4sAnd2s()
    00844 {
    00845         leda_node u, v;
    00846         leda_edge e;
    00847         bool bAllType1SameDir;
    00848         int dir;
    00849
    00850         forall_nodes(u, *this)
    00851         {
    00852                 if (NodeType(u) == SHOCK_1 && outdeg(u) > 0)
    00853                 {
    00854                         bAllType1SameDir = true;
    00855                         dir = GetBranchDir(u, GetFirstChild(u));
    00856
    00857                         // Check if all children of <u> are type 1
    00858                         forall_adj_nodes(v, u)
    00859                         {
    00860                                 if (NodeType(v) != SHOCK_1 || GetBranchDir(v, u) != dir)
    00861                                 {
    00862                                         bAllType1SameDir = false;
    00863                                         break;
    00864                                 }
    00865                         }
    00866
    00867                         if (bAllType1SameDir)
    00868                         {
    00869                                 const ShockBranch& pts = GetSGNode(u)->m_shocks;
    00870                                 int type = dir == 1 ? SHOCK_2:SHOCK_4;
    00871
    00872                                 if (type == SHOCK_2 && outdeg(u) != 1)
    00873                                         return; // we only consider necks between two branches
    00874
    00875                                 SGNode* pNewNode = new SGNode(GetNextIdx(), type);
    00876                                 pNewNode->m_shocks.AddTail(GetJointPoint(u, GetFirstChild(u)));
    00877                                 pNewNode->ComputeDerivedValues();
    00878
    00879                                 bool localMinOrMax = false;
    00880
    00881                                 //Before inserting the node, in case of a type 4,
    00882                                 // we must make sure it's really a local maximum pt
    00883                                 if (type == SHOCK_4)
    00884                                 {
    00885                                         if (pNewNode->m_shocks > GetSGNode(u)->m_shocks)
    00886                                         {
    00887                                                 localMinOrMax = true;
    00888
    00889                                                 forall_adj_nodes(v, u)
    00890                                                 {
    00891                                                         if (!(pNewNode->m_shocks > GetSGNode(v)->m_shocks))
    00892                                                         {
    00893                                                                 localMinOrMax = false;
    00894                                                                 break;
    00895                                                         }
    00896                                                 }
    00897                                         }
    00898                                 }
    00899                                 else if (type == SHOCK_2)
    00900                                 {
    00901                                         if (GetSGNode(u)->m_shocks > pNewNode->m_shocks && GetSGNode(GetFirstChild(u))->m_shocks > pNewNode->m_shocks)
    00902                                                 localMinOrMax = true;
    00903                                 }
    00904
    00905                                 if (localMinOrMax) // We have either a type 2 or a type 4
    00906                                 {
    00907                                 v = NewNode(pNewNode);
    00908
    00909                                 forall_adj_edges(e, u)
    00910                                         move_edge(e, v, target(e));
    00911
    00912                                 NewEdge(u, v); // keep the old order u -> new (v) -> u's children
    00913                         }
    00914                         else
    00915                                 delete pNewNode;
    00916                         }
    00917                 }
    00918         }
    00919 }
    00920
    00921 void ShockGraph::ComputeNodeRole(leda_node v)
    00922 {
    00923         NODE_ROLE role;
    00924         int nType = NodeType(v);
    00925
    00926         if (nType == ROOT)
    00927                 role = BIRTH;
    00928         else if (nType == SHOCK_2 || (nType == SHOCK_3 && outdeg(v) == 0))
    00929                 role = UNION;
    00930         else if (outdeg(v) == 0)
    00931                 role = DEATH;
    00932         else
    00933                 role = PROTRUSION;
    00934
    00935         UnsafeGetSGNode(v)->SetNodeRole(role);
    00936 }
    00937
    00941 sg::DDSkeleton* ShockGraph::ComputeDDSkeleton(const char* szPPMFileName)
    00942 {
    00943         using namespace sg;
    00944
    00945         cout << endl << "Reading " << szPPMFileName << "... " << flush;
    00946
    00947         int i, j, rows, cols, format;
    00948         pixval maxval;
    00949
    00950         FILE* fp = fopen(szPPMFileName, "r");
    00951
    00952         if (fp == NULL)
    00953     {
    00954                 cerr << "Can't open " << szPPMFileName << endl;
    00955                 return NULL;
    00956     }
    00957
    00958         ppm_readppminit(fp, &cols, &rows, &maxval, &format);
    00959         fclose(fp);
    00960
    00961         //  cout << "Getting array...";
    00962         char *a = get_arr(&rows, &cols, szPPMFileName);
    00963
    00964         //  cout << "done." << endl;
    00965         // first, we must make a Shape
    00966         sg::SimpleShapeMaker ssm(cols, rows);
    00967
    00968         // update the SimpleShape
    00969         for (i=0; i < rows; i++)
    00970                 for (j=0; j < cols; j++)
    00971                         ssm(j,i) = a[i*cols + j];
    00972
    00973         // free the input array
    00974         free(a);
    00975
    00976         // the shape is being made here
    00977         SimpleShape ss = *((SimpleShape*)ssm.getShape());
    00978
    00979         //double xmin,xmax,ymin,ymax; //use the member variables instead
    00980         ss.getBounds(&xmin,&xmax, &ymin,&ymax);
    00981
    00982         //std::cerr << "(xmin,ymin): ("<<xmin<<", "<<ymin<<")\n";
    00983         //std::cerr << "(xmax,ymax): ("<<xmax<<", "<<ymax<<")\n";
    00984
    00985         // we can get the contour of a SimpleShape as is consists of a
    00986         // single Curve (actually a DiscreteSegCurve)
    00987
    00988         //DiscreteSegCurve dsc = *(DiscreteSegCurve *)ssm.getContour();
    00989         //std::cerr << dsc << "\n";
    00990
    00991         // now, we can create a DistanceTransform object
    00992         DistanceTransform dt(&ss);
    00993
    00994         // and a DivergenceMap
    00995         DivergenceMap dm(dt);
    00996
    00997         cout << "making skeleton... " << flush;
    00998
    00999         // to supply to the SkeletonMaker
    01000         // cout << "Making the skeleton...";
    01001         DivergenceSkeletonMaker dsm(dm);
    01002         double thresh=-2.0, mag=1.0;
    01003
    01004         ASSERT(compParams.dSkTreshold >= 0);
    01005
    01006         return dsm.getDiscreteDivergenceSkeleton(1.0 / mag, -compParams.dSkTreshold);
    01007 }
    01008
    01013 void ShockGraph::DetectLigatureNodes(sg::DDSkeleton* sk)
    01014 {
    01015         using namespace sg;
    01016
    01017         DDSEdgeVect& edges = sk->getEdges();
    01018         DDSEdgeVect::iterator I;
    01019         DDSEdge* e;
    01020
    01021         vector<Curve*>* v = sk->getShape()->getCurves();
    01022         Curve* c = *v->begin();
    01023         DiscreteSegCurve dsc = *(DiscreteSegCurve *)c;
    01024
    01025         PRINT(cout, v->size());
    01026
    01027         for(I = edges.begin(); I != edges.end(); I++)
    01028         {
    01029                 e = *I;
    01030
    01031                 Point pRt1 = dsc.atT(e->rightT1());
    01032                 Point pRt2 = dsc.atT(e->rightT2());
    01033                 Point pLt1 = dsc.atT(e->leftT1());
    01034                 Point pLt2 = dsc.atT(e->leftT2());
    01035
    01036                 PRINT_OPEN(cout, e->hasRight());
    01037                 PRINT(cout, pRt1);
    01038                 PRINT(cout, pRt2);
    01039                 PRINTN(cout, e->hasLeft());
    01040                 PRINT(cout, pLt1);
    01041                 PRINT(cout, pLt2);
    01042                 PRINTN(cout, e->getFluxPoints().front().val);
    01043                 PRINT_CLOSE(cout, e->getFluxPoints().back().val);
    01044         }
    01045 }
    01046
    01053 bool ShockGraph::ComputeSGFromDDSkeleton(sg::DDSkeleton* sk)
    01054 {
    01055         using namespace sg;
    01056
    01057         cerr << "labeling skeleton... " << flush;
    01058
    01059         //DetectLigatureNodes(sk);
    01060
    01061         SGNode* pNode = new SGNode(GetNextIdx(), SHOCK_base);
    01062         leda_node u = NewNode(pNode);
    01063         DDSEdge* se = *sk->getEdges().begin();
    01064
    01065         SGNode* pFirstNode = CopyNodeInfo(se, u, pNode);
    01066
    01067         ConnectNodes(u, pNode, se, se->getN1());
    01068         ConnectNodes(u, pFirstNode, se, se->getN2());
    01069
    01070         Insert4sAnd2s();
    01071         //ShowGraph(*this);
    01072
    01073         bool bContinue;
    01074         leda_edge e;
    01075         int nNoInfLoop = 0, i;
    01076
    01077         // Collect the edges that need to be reverted according to time of formation.
    01078         // Edges should be directed from last-to-form to earliest-to-form.
    01079         do
    01080         {
    01081                 ASSERT(nNoInfLoop++ < 1000);
    01082                 bContinue = false;
    01083
    01084                 forall_edges(e, *this)
    01085                 {
    01086                         const ShockBranch& tgt = GetSGNode(target(e))->m_shocks;
    01087                         const ShockBranch& src = GetSGNode(source(e))->m_shocks;
    01088
    01089                         if (tgt > src)
    01090                         {
    01091                                 RevertEdge(e);
    01092                                 bContinue = true;
    01093                         }
    01094                 }
    01095         } while (bContinue);
    01096
    01097         // Nodes that do not have a parent should be children of the root node.
    01098         leda_node root = NewNode(new SGNode("#", ROOT));
    01099         leda_node v;
    01100         int nodeType;
    01101
    01102         // BTW, type 3 nodes with 2 type 1 parent/children should
    01103         // be relabeled as 2/4 respectivelly.
    01104         forall_nodes(v, *this)
    01105         {
    01106                 if (v != root && indeg(v) == 0)
    01107                         NewEdge(root, v);
    01108
    01109                 nodeType = NodeType(v);
    01110
    01111                 if (nodeType == SHOCK_3) // May need to be relabeled.
    01112                         Relabel3As2or4(v);
    01113
    01114                 // Make sure the grammar is satisfied
    01115                 ASSERT (nodeType == ROOT ||
    01116                         nodeType == SHOCK_1 ||
    01117                         nodeType == SHOCK_3 ||
    01118                         (nodeType == SHOCK_2 && indeg(v) == 2 && outdeg(v) == 0) ||
    01119                         (nodeType == SHOCK_4 && indeg(v) == 1 && NodeType(GetFirstParent(v)) == ROOT));
    01120
    01121                 ComputeNodeRole(v);
    01122         }
    01123
    01124         ComputeDerivedValues();
    01125
    01126         cerr << "(" << GetObjName() << ", " << GetViewNumber() << ") done!" << flush;
    01127
    01128         return true;
    01129 }
    01130
    01131 bool ShockGraph::ComputeFromPPMFile2(const char* szFileName, const SGCP& sgparams)
    01132 {
    01133         Clear(); //make sure the current shock graph is empty
    01134
    01135         compParams = sgparams;
    01136         SetDAGLbl(szFileName);
    01137
    01138         sg::DDSkeleton* pDDS = ComputeDDSkeleton(szFileName);
    01139
    01140         if (!pDDS)
    01141                 return false;
    01142
    01143         m_pSkeleton = new McGillSkeleton(pDDS);
    01144
    01145         return ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
    01146 }
    01147
    01148 bool ShockGraph::RecomputeFromSkeleton(const SGCP& sgparams)
    01149 {
    01150         // We must empty the graph but without deleting the skeleton, label, mins and maxs.
    01151         String strLbl = GetDAGLbl();
    01152         DAG::Clear();
    01153         nLastIndexUsed = 0;
    01154
    01155         compParams = sgparams;
    01156         SetDAGLbl(strLbl);
    01157
    01158         return ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
    01159 }
    01160
    01161 void SmoothFluxPoints(sg::FluxPointList& fpl, const double& dSigma, const double& dSigmaRange)
    01162 {
    01163         int i, nLen = fpl.size();
    01164
    01165         double* in = new double[nLen];
    01166         double* out = new double[nLen];
    01167
    01168         for (i = 0; i < nLen; i++)
    01169                 in[i] = fpl[i].dist;
    01170
    01171         GaussianSmooth(in, out, nLen, dSigma, dSigmaRange);
    01172
    01173         for (i = 0; i < nLen; i++)
    01174                 fpl[i].dist = out[i];
    01175
    01176         delete[] in;
    01177         delete[] out;
    01178 }
    01179
    01189 void GaussianSmooth(double *input, double *output, int length,
    01190                                         const double& sigma, const double& sigma_range)
    01191 {
    01192         static double* kernel = NULL;
    01193         static double prev_sigma, prev_sigma_range;
    01194         static double kernelsum;
    01195         static int w2;
    01196         int i, x;
    01197
    01198         if (kernel && (sigma != prev_sigma || sigma_range != prev_sigma_range))
    01199         {
    01200                 delete[] kernel;
    01201                 kernel = NULL;
    01202         }
    01203
    01204         // first compute the kernel
    01205         if (!kernel)
    01206         {
    01207                 w2 = (int)rint(sigma_range * sigma);
    01208                 if (w2 < 1) w2 = 1;
    01209                 int w = (2 * w2) + 1;
    01210
    01211                 kernel = new double[w];
    01212                 kernelsum = 0.0;
    01213                 prev_sigma = sigma;
    01214                 prev_sigma_range = sigma_range;
    01215
    01216                 for (double* kp = kernel, x = -w2; x <= w2; kp++,x++)
    01217                 {
    01218                         *kp = (1 / (sqrt(2 * M_PI) * sigma)) * exp((double)(x * x)/(-2 * sigma * sigma));
    01219                         kernelsum += *kp;
    01220                 }
    01221         }
    01222
    01223         PRINT(cerr, kernelsum);
    01224
    01225         double sum;
    01226
    01227         // now perform the convolution
    01228         for (x = 0; x < length; x++)
    01229         {
    01230                 sum = 0.0;
    01231
    01232                 for (i = -w2; i <= w2; i++)
    01233                 {
    01234                         if ((x + i) < 0){
    01235                                 sum += kernel[i + w2] * input[0]; /* extend left side */
    01236                         }
    01237                         else if ((x + i) >= length){
    01238                                 sum += kernel[i + w2] * input[length - 1]; /* extend right side */
    01239                         }
    01240                         else{
    01241                                 sum += kernel[i + w2] * input[x + i];
    01242                         }
    01243                 }
    01244
    01245                 output[x] = sum / kernelsum;
    01246         }
    01247 }
    01248
    01253 int ShockGraph::GetBranchDir(leda_node u, leda_node wrtV) const
    01254 {
    01255         const ShockInfo& pt = GetJointPoint(u, wrtV);
    01256
    01257         return GetSGNode(u)->GetBranchDir(pt == GetSGNode(u)->m_shocks[0] ?
    01258                 ShockBranch::FORWARD:ShockBranch::BACKWARD);
    01259 }
    01260
    01261 const ShockInfo& ShockGraph::GetJointPoint(leda_node u, leda_node v) const
    01262 {
    01263         const SGNode* pNodeU = GetSGNode(u);
    01264         const SGNode* pNodeV = GetSGNode(v);
    01265
    01266         const ShockInfo& uStart = pNodeU->m_shocks[0];
    01267         const ShockInfo& uEnd   = pNodeU->m_shocks.GetTail();
    01268
    01269         const ShockInfo& vStart = pNodeV->m_shocks[0];
    01270         const ShockInfo& vEnd   = pNodeV->m_shocks.GetTail();
    01271
    01272         if (uStart == vStart || uStart == vEnd)
    01273                 return uStart;
    01274         else if (uEnd == vStart || uEnd == vEnd)
    01275                 return uEnd;
    01276
    01277         ASSERT(false);
    01278         return vStart;
    01279 }
    01280
    01281
    01282