---

ShockGraph.cpp

Go to the documentation of this file.

#include "DAG.h"
#include "ShockGraph.h"
#include "HelperFunctions.h"
#include "Exceptions.h"
#include <float.h>
#include <stdio.h>
#include <ctype.h>

#include "SGDistMeasurer.h"

#include <LEDA/graphwin.h>
#include <LEDA/graph.h>
#include <LEDA/node_array.h>

#include "VisualDAG.h"
#include "ApproxPoly.h"

#define NSAMPLES         5
#define VDIM             4
#define DIMxSAMPLES  (VDIM * NSAMPLES)

#define SA_GROWING_FACTOR 50

#ifdef WIN32
#define FILE_SEP '\\'
#else
#define FILE_SEP '/'
#endif

struct ER // Pair of edge and radius values
{
        leda_edge e;
        double r; // avg radius of target(e)
        
        ER() { e = nil; r = 0.0; }
        ER(leda_edge edge, const double& radius) { e = edge; r = radius; }
        
        static int CompareEdgesDescending(const void *elem1, const void *elem2 ) 
        {
                ER* p1 = (ER*) elem1;
                ER* p2 = (ER*) elem2;
                
                return (p1->r < p2->r) ? 1:(p1->r > p2->r) ? -1:0;
        }
};

void SmoothFluxPoints(sg::FluxPointList& fpl, const double& dSigma, const double& dSigmaRange);
void GaussianSmooth(double *input, double *output, int length, const double& sigma, const double& sigma_range);

const char* ShockGraph::GetNextIdx()
{ 
        static char szNum[100];
        
        sprintf(szNum, "%d", ++nLastIndexUsed);
        return szNum;
}

ShockGraph::ShockGraph()
{
        area = 0;
        perimeter = 0;
        xmin = ymin = 0;
        xmax = ymax = 0;
        nLastIndexUsed = 0;
        m_pSkeleton = NULL;
        
        ASSERT(!pDistMeasurer);
        
        if (s_matchParams.nNodeDistFunc == 1)
                pDistMeasurer = new SGDistMeasurer1;
        else if (s_matchParams.nNodeDistFunc == 2)
                pDistMeasurer = new SGDistMeasurer2;
        else if (s_matchParams.nNodeDistFunc == 3)
                pDistMeasurer = new SGDistMeasurer3;
        else
                ASSERT(false);
}

ShockGraph::~ShockGraph()
{
        delete pDistMeasurer;
        //WARNING(m_pSkeleton, "Skeletton destructor isn't working at this point. Memory Leak!");
        //delete m_pSkeleton;
}

void ShockGraph::Clear()
{ 
        DAG::Clear();
        
        area = 0.0;
        perimeter = 0.0;
        xmin = ymin = 0;
        xmax = ymax = 0;
        nLastIndexUsed = 0;

        delete m_pSkeleton;
        m_pSkeleton = NULL;
}

void ShockGraph::ComputeObjectViewName()
{
        int i, j, k, lastFileSep;
        String viewName;
        const String path = GetDAGLbl();
        const int nLen = path.Len();
        
        strObjectName.Clear();
        
        // Find first char of the last directory name in the path
        
        for (i = nLen - 1; i >= 0 && path[i] != FILE_SEP; i--);
        
        lastFileSep = i;
        
        if (i > 0) // we have a string with format '/ObjName/prefixNNN.ext'
        {
                // Search for next '/'
                for (j = i - 1; j >= 0 && path[j] != FILE_SEP; j--);
                
                // j point to the second '/'
                strObjectName.ReSize(i - j);
                
                for (j++, k = 0;  j < i; j++, k++)
                        strObjectName[k] = path[j];
                
                strObjectName[k] = '\0';
                
                // Find beginning of view number
                for (i++; i < nLen && !isdigit(path[i]); i++);
        }
        else // we have a string with format 'ObjNameNNN.ext'
        {
                // Find beginning of view number
                for (i = 0; i < nLen && !isdigit(path[i]); i++);
                
                // Copy obj name
                
                strObjectName.ReSize(i + 1);
                
                for (j = 0; j < i; j++)
                        strObjectName[j] = path[j];
                
                strObjectName[j] = '\0';
        }
        
        // Find view number
        viewName.ReSize(nLen - i + 1);
        
        for (k = 0;  i <= nLen; i++, k++)
                viewName[k] = path[i];
        
        nViewNumber = (viewName[0] != '\0') ? atoi(viewName):1;
}

void ShockGraph::ComputeDerivedValues()
{
        DAG::ComputeDerivedValues();
        ComputeObjectViewName();
}

istream& ShockGraph::Read(istream& is, bool bOnlyDataForMatching /*= false*/)
{
        DAG::Read(is, bOnlyDataForMatching);
        
        is.read((char*) &area, sizeof(area));
        is.read((char*) &perimeter, sizeof(perimeter));
        is.read((char*) &xmin, sizeof(xmin));
        is.read((char*) &xmax, sizeof(xmax));
        is.read((char*) &ymin, sizeof(ymin));
        is.read((char*) &ymax, sizeof(ymax));

        if (!bOnlyDataForMatching)
        {
                // pSkeleton works as a boolean variable here
                is.read((char*) &m_pSkeleton, sizeof(m_pSkeleton));
                if (m_pSkeleton)
                {
                m_pSkeleton = new McGillSkeleton();
                m_pSkeleton->Read(is);
                }
        }
        
        return is;
}

ostream& ShockGraph::Write(ostream& os) const
{
        DAG::Write(os);
        
        os.write((char*) &area, sizeof(area));
        os.write((char*) &perimeter, sizeof(perimeter));
        os.write((char*) &xmin, sizeof(xmin));
        os.write((char*) &xmax, sizeof(xmax));
        os.write((char*) &ymin, sizeof(ymin));
        os.write((char*) &ymax, sizeof(ymax));

        // pSkeleton works as a boolean variable here
        os.write((char*) &m_pSkeleton, sizeof(m_pSkeleton));
        if (m_pSkeleton)
          m_pSkeleton->Write(os);
        
        return os;
}

void ShockGraph::Print(ostream& os) const
{
        DAG::Print(os);
        
        os << "\nShockGraph specific:\n";
        
        PRINT_OPEN(os, area);
        PRINT(os, perimeter);
        PRINT(os, xmin);
        PRINT(os, xmax);
        PRINT(os, ymin);
        PRINT_CLOSE(os, ymax);
}

DAG& ShockGraph::operator=(const DAG& rhs)
{
        DAG::operator=(rhs);
        
        const ShockGraph* pRhs = dynamic_cast<const ShockGraph*>(&rhs);
        
        if (!pRhs)
                THROW_EXCEPTION("Invalid pointer type.");
        
        area      = pRhs->area;
        perimeter = pRhs->perimeter;
        
        xmin = pRhs->xmin;
        ymin = pRhs->ymin;
        xmax = pRhs->xmax;
        ymax = pRhs->ymax;
        
        return *this;
}

double ShockGraph::NodeDistance(leda_node u, const DAG& from, leda_node v) const
{
        ASSERT(pDistMeasurer);
        
        return pDistMeasurer->GetDistance(*this, u, from, v);
}

bool ShockGraph::AreNodesRelated(leda_node u, const DAG& from, leda_node v) const
{
        return true;
}

DAG* ShockGraph::CreateObject() const
{
        return new ShockGraph();
}

DAGNodePtr ShockGraph::CreateNodeObject(NODE_LABEL lbl) const
{
        return new SGNode(lbl, TEST_NODE);
}

DAGNodePtr ShockGraph::ReadNode(istream& is) const
{
        DAGNodePtr node(new SGNode);
        
        node->Read(is);
        
        return node;
}

String ShockGraph::ClassName() const
{
        return "ShockGraph";
}

NODE_LABEL ShockGraph::GetNodeLbl(leda_node v) const
{
        char szLbl[100];
        const char* szNodeLbl = (const char*)DAG::GetNodeLbl(v);
        
        sprintf(szLbl, "%s:%d", szNodeLbl, NodeType(v) % 100);
        
        return szLbl;
}

/*************************************************************************************************************************/
/* Interface with shock graph code */

bool ShockGraph::ComputeFromPPMFile(const char* szFileName, 
                                                                        double cutoff, double sigma, double range)
{
#ifdef COMPILE_OLD_SG_CODE
        MGShockGraph sg;
        
        if (!sg.ComputeFromPPMFile(szFileName, cutoff, sigma, range) || !Translate(sg))
                return false;
        
        SetDAGLbl(szFileName);
        
        ComputeDerivedValues();
#endif //COMPILE_OLD_SG_CODE
        
        return true;
}

#ifdef COMPILE_OLD_SG_CODE

bool ShockGraph::ComputeFromMemory(Pixmap& image, double cutoff, double sigma, double range)
{
        MGShockGraph sg;
        
        if (!sg.ComputeFromMemory(image, cutoff, sigma, range) || !Translate(sg))
                return false;
        
        SetDAGLbl(image.name);
        
        ComputeDerivedValues();
        
        return true;
}

bool ShockGraph::Translate(const MGShockGraph& sg)
{
        leda_d_array<long, leda_node> map;
        int i, j, ptCount;
        SGNode* pNode;
        double seg1, seg2;
        
        Clear();
        
        for (i = 0; i < sg.GetNodeCount(); i++)
        {
                const Vertex* v = sg.GetNode(i);
                MGShockGraph::SGSegment branch = sg.GetSegment(v);
                
                pNode = new SGNode(branch.GetId(), branch.GetType());
                
                branch.GetContourLength(seg1, seg2);
                pNode->SetContourLength(seg1, seg2);
                
                ptCount = branch.GetPointCount();
                pNode->m_shocks.ReSize((ptCount > 0) ? ptCount:0);
                
                for (j = 0; j < ptCount; j++)
                {
                        const MGShockGraph::SGPoint& pt = branch[j];
                        ShockInfo& si = pNode->m_shocks[j];
                        
                        si.xcoord = pt.x;
                        si.ycoord = pt.y;
                        si.radius = pt.r;
                        si.speed = pt.s;
                        si.dr_ds = pt.dr_ds;
                        si.color = pt.color;
                        
                        // Set bounding rectangle for drawing purposes
                        if (pt.x < xmin) xmin = pt.x;
                        if (pt.x > xmax) xmax = pt.x;
                        if (pt.y < ymin) ymin = pt.y;
                        if (pt.y > ymax) ymax = pt.y;
                }
                
                map[(long)v] = NewNode(pNode);
        }
        
        for (i = 0; i < sg.GetEdgeCount(); i++)
        {
                const Edge* e = sg.GetEdge(i);
                
                ASSERT(map.defined((long)e->vFrom) && map.defined((long)e->vTo));
                
                // From an To nodes have to be swapped because we consider the
                // the directions are from the root to the leaves.
                NewEdge(map[(long)e->vTo], map[(long)e->vFrom]);
        }
        
        return true;
}

#endif //COMPILE_OLD_SG_CODE


// test code from McGill
char* get_arr(int *rows, int *cols, const char* shapefile)
{
        char *arr;
        int k, i, j;
        
        pixval maxval;
        pixel **pixels;
        
        FILE* ppmshape = fopen(shapefile, "r");  
        
        pixels = ppm_readppm(ppmshape, cols, rows, &maxval);
        fclose(ppmshape);
        
        arr = (char*)malloc((*cols)*(*rows));
        
        for(i = 0, k = 0; i < *rows; i++)
                for(j = 0; j < *cols; j++, k++)
                {
                        char r, g, b;
                        
                        r = (0xff * PPM_GETR(pixels[i][j])) / maxval;
                        g = (0xff * PPM_GETR(pixels[i][j])) / maxval;
                        b = (0xff * PPM_GETR(pixels[i][j])) / maxval;
                        
                        if(r==0 && b==0 && g==0)
                                arr[k] = 0;
                        else
                                arr[k] = 0xff;
                }
                
                ppm_freearray(pixels, *rows);
                
                return arr;
};

void ShockGraph::LabelEndPoints(SGNode* pNode, sg::FluxPointList& fpl, sg::DDSNode* n1, sg::DDSNode* n2)
{
        sg::DDSEdgeVect *e0, *eN;
        
        if (n1->getPoint() == fpl[0].p)
        {
                ASSERT(n2->getPoint() == fpl[fpl.size() - 1].p);
                
                e0 = &n1->getEdges();
                eN = &n2->getEdges();
        }
        else
        {
                ASSERT(n1->getPoint() == fpl[fpl.size() - 1].p && n2->getPoint() == fpl[0].p);
                
                e0 = &n2->getEdges();
                eN = &n1->getEdges();
        }

        pNode->m_nEndPt0 = (e0->size() == 1) ? TERMINAL_POINT:JOINT_POINT;
        pNode->m_nEndPtN = (eN->size() == 1) ? TERMINAL_POINT:JOINT_POINT;
}

SGNode* ShockGraph::CopyNodeInfo(sg::DDSEdge* ddsEdge, leda_node u, SGNode* pNode)
{
        sg::FluxPointList& fpl = ddsEdge->getFluxPoints();
        int num_pts = fpl.size();
        ASSERT(num_pts > 0);
        double dx, dy;

        //SmoothFluxPoints(fpl, compParams.dSigma, compParams.dSigmaRange);
        
        pNode->m_shocks.ReSize(num_pts);
        
        for (int j = 0; j < num_pts; j++)
        {
                const sg::FluxPoint& pt = fpl[j];
                ShockInfo& si = pNode->m_shocks[j];
                
                si.xcoord = pt.p.x;
                si.ycoord = pt.p.y;
                si.radius = fabs(pt.dist);
                si.color = 0; //pt.col; // not used anymore. see type instead
                
                // Some NaN have been found in a value for radius
                if (isnan(si.radius) && !finite(si.radius))
                {
                        cerr << "\nShock point has a NaN as radius." << endl;
                        // Try to recover from this
                        si.radius = (j > 0) ? pNode->m_shocks[j - 1].radius:0.5;
                }
                
                if (j > 0)
                {
                        ShockInfo& prevSi = pNode->m_shocks[j - 1];
                        
                        dx = si.xcoord - prevSi.xcoord;
                        dy = si.ycoord - prevSi.ycoord;
                        
                        si.dr = si.radius - prevSi.radius;
                        si.speed = sqrt(dx * dx + dy * dy);
                        si.dr_ds = si.dr / si.speed;

                        // dir is either -1, +1 or 0
                        si.dir = (fabs(si.dr_ds) <= compParams.dMinSlope) ? 0:int(si.dr / fabs(si.dr));
                }
                else
                {
                        si.dr = 0;
                        si.speed = 0;
                        si.dr_ds = 0;
                        si.dir = 0;
                }
        }

        LabelEndPoints(pNode, fpl, ddsEdge->getN1(), ddsEdge->getN2());

        return GroupShockPoints(u, pNode);
}

SGNode* ShockGraph::GroupShockPoints(leda_node u, SGNode* pNode)
{
        ShockBranch& pts = pNode->m_shocks;
        //int nSize = pts.GetSize();
        SGNode* pNewNode;
        SGNode* pReturnNode = NULL;

        int i, from;
        double m, nextSlope;
        
        int d0, dN;     // if d0 or dN == 1, means their end points were chopped off
        POINTS data = pNode->GetVelocityRadiusArray(d0, dN);
        
        double x0 = data[0].x;
                
        // Time to fit the data with a few lines
        ApproxPoly poly(data.GetSize() / compParams.dMinError, compParams.dMinSlope, 10.0);
        
        if (DAG::IsDbgMode())
        {
                cout << "disp('Fitting node: " << pNode->GetNodeLbl() << " in " << GetDAGLbl() << "');\n";
                poly.SetDbgMode(true); // outputs matlab data in dbg mode
        }
        
        poly.Fit(data);
        
        int nSize = poly.m_knots.GetSize();
        ASSERT(nSize >= 1);
        int dir, ptIdx;
        SmartArray<SEGMENT> approx_segs;

        for (from = 0, i = 0; i < nSize; i++)
        {
                approx_segs.AddTail(poly.m_knots[i].seg);
                
                dir = poly.m_knots[i].dir;
                ptIdx = poly.m_knots[i].nIndex + d0;
                
                if (i < nSize - 1)
                {
                        // If the segment is too small, skip it.
                        if (ptIdx - from < 2)
                                continue;
                                
                        // If there is a significant change in acceleration,
                        // we create a new node, even though the direction is the same.
                        if (dir == poly.m_knots[i + 1].dir)
                        {
                                if (dir == 0)
                                        continue;
                                
                                //poly.CompAcuteAngle(i, i + 1) < compParams.dMaxAccelChg)
                                double m0 = fabs(poly.m_knots[i].seg.m);
                                double m1 = fabs(poly.m_knots[i + 1].seg.m);
                                double mr = fabs(m1 - m0) / MAX(m0, m1);
                                
                                if (DAG::IsDbgMode())
                {
                        cout << "disp('Acc change: " << mr << ", Max:" << compParams.dMaxAccelChg << "');\n";
                        poly.SetDbgMode(true); // outputs matlab data in dbg mode
                }
                                
                                if (mr < compParams.dMaxAccelChg)
                                        continue;
                        }

                        pNewNode = SplitNode(u, pNode, ptIdx - from, approx_segs[approx_segs.GetSize() - 1], poly.m_knots[i + 1].seg);
                }
                else
                        pNewNode = pNode;

                // Record the first node of the chain
                if (pReturnNode == NULL)
                        pReturnNode = pNewNode;

                pNewNode->SetNodeType(dir == 0 ? SHOCK_3:SHOCK_1);
                pNewNode->m_shocks.SetDir(dir);
                
                // Update the fitting segments. Note that if <from> != 0, the x coord of points in the
                // branch will be shifted. SetSegments will correct this gived the displacement.
                pNewNode->m_shocks.SetSegments(approx_segs, from == 0 ? 0.0:approx_segs[0].p0.x - x0);
                
                pNewNode->ComputeDerivedValues();

                from = ptIdx;
                approx_segs.Clear();
        }

        /*for (int i = 1; i < nSize; i++) 
        {
                // Do we need to split the node?
                if (pts[i].type != pts[0].type) 
                {
                        pNewNode = SplitNode(u, pNode, i - 1, pts[0].type);
                        pNewNode->ComputeDerivedValues();
                        GroupShockPoints(u, pNode);
                        break;
                }
        }
        
        pNode->ComputeDerivedValues();*/

        return pReturnNode;
}

/*
@brief Splits in two the shock points in a node. It returns the new node created that
contains the first points. It also links the two nodes.
*/
SGNode* ShockGraph::SplitNode(leda_node v, SGNode* pNode, int nEnd, SEGMENT& leftSeg, SEGMENT& rightSeg)
{
        SGNode* pNewNode = new SGNode(GetNextIdx(), 0);
        ShockBranch& newPts = pNewNode->m_shocks;
        ShockBranch& pts = pNode->m_shocks;

        // We must repeat the endpoint where the nodes connect. nEnd is inclusive
        pNewNode->m_shocks = pts.Left(nEnd);

        // Now remove the left part. nEnd is inclusive
        pts = pts.Right(nEnd);
        
        // Create new node
        leda_node u = NewNode(pNewNode);
        
        // Update the m_nEndPt0 and m_nEndPtN member variables of both nodes
        pNewNode->m_nEndPt0 = pNode->m_nEndPt0;
        pNewNode->m_nEndPtN = SPLIT_POINT;
        pNode->m_nEndPt0 = SPLIT_POINT;
        
        // Update the Y value (radius) at the joint pt, to avoid outliers in this crusial
        // function. i.e. type 2 and 4 will depend on accurate values here
        double jointRadius = (leftSeg.Y1() + rightSeg.Y0()) / 2.0;
        pNewNode->m_shocks[pNewNode->m_shocks.GetSize() - 1].radius = jointRadius;
        pts[0].radius = jointRadius;
        leftSeg.p1.y = jointRadius;
        rightSeg.p0.y = jointRadius;
        
        // If 'v' has a parent node, move it to be a parent of 'u' instead
        leda_edge e = first_in_edge(v);
        
        if (e != nil)
    {
                NewEdge(source(e), u);
                ASSERT(in_succ(e) == nil); // There shoukd be only one edge
                DelEdge(e);
    }
        
        // Link u to v
        NewEdge(u, v);
        
        return pNewNode;
}

void ShockGraph::ConnectNodes(leda_node u, SGNode* pNode, sg::DDSEdge* e, sg::DDSNode* n)
{
        using namespace sg;
        
        // Then, keep creating and connecting nodes
        DDSEdgeVect& edges = n->getEdges();
        DDSEdgeVect::iterator I;
        DDSEdge* ee;
        SGNode* pNewNode;
        leda_node v;

        for(I = edges.begin(); I != edges.end(); I++)
        {
                ee = *I;
                
                if (e != ee)
                {
                        pNewNode = new SGNode(GetNextIdx(), SHOCK_base);
                        v = NewNode(pNewNode);
                        NewEdge(u, v);
                        //CopyNodeInfo(ee->getFluxPoints(), v, pNewNode);
                        CopyNodeInfo(ee, v, pNewNode);
                        ConnectNodes(v, pNewNode, ee, ee->getN1() != n ? ee->getN1():ee->getN2());
                }
        }
}

void ShockGraph::RevertEdge(leda_edge pxEdge)
{
        leda_node parent = source(pxEdge);
        leda_node x = target(pxEdge);
        leda_edge e;

        rev_edge(pxEdge);
        
        // Move all children of x's parent (that are also connected to x)
        // to the left side of x
        const ShockInfo& inPt = GetJointPoint(parent, x);
                
        forall_adj_edges(e, parent)
                if (GetJointPoint(parent, target(e)) == inPt)
                        move_edge(e, x, target(e));
}

void ShockGraph::Relabel3As2or4(leda_node v)
{
        ASSERT(NodeType(v) == SHOCK_3);
        
        int parentType = NodeType(GetFirstParent(v));
        
        if (parentType == ROOT)
                UnsafeGetSGNode(v)->SetNodeType(SHOCK_4);
        else if(outdeg(v) == 0 && indeg(v) == 2)
                UnsafeGetSGNode(v)->SetNodeType(SHOCK_2);

        // If we are interested in spliting 3's into left and
        // right nodes, we can do it here.
        if (compParams.nSlipt3s != 0 && parentType == ROOT)
        {
                const SGNode& x = *GetSGNode(v);
                
                const ShockInfo& leftPt = x[0];
                const ShockInfo& rightPt = x[x.GetSize() - 1];
                
                leda_node leftChild = NewNode(new SGNode(x));
                leda_node rightChild = NewNode(new SGNode(x));
                leda_edge e;
                
                forall_adj_edges(e, v)
                {
                        if (GetJointPoint(v, target(e)) == leftPt)
                                move_edge(e, leftChild, target(e));
                        else
                        {
                                ASSERT(GetJointPoint(v, target(e)) == rightPt);
                                move_edge(e, rightChild, target(e));
                        }
                }
                                
                // Link x to its left and right sides
                NewEdge(v, rightChild);
                NewEdge(v, leftChild);
        }
}

void ShockGraph::Insert4sAnd2s()
{
        leda_node u, v;
        leda_edge e;
        bool bAllType1SameDir;
        int dir;
        
        forall_nodes(u, *this)
        {
                if (NodeType(u) == SHOCK_1 && outdeg(u) > 0)
                {
                        bAllType1SameDir = true;
                        dir = GetBranchDir(u, GetFirstChild(u));
                        
                        // Check if all children of <u> are type 1
                        forall_adj_nodes(v, u)
                        {
                                if (NodeType(v) != SHOCK_1 || GetBranchDir(v, u) != dir)
                                {
                                        bAllType1SameDir = false; 
                                        break;
                                }
                        }
                                
                        if (bAllType1SameDir)
                        {
                                const ShockBranch& pts = GetSGNode(u)->m_shocks;
                                int type = dir == 1 ? SHOCK_2:SHOCK_4;
                                
                                if (type == SHOCK_2 && outdeg(u) != 1)
                                        return; // we only consider necks between two branches

                                SGNode* pNewNode = new SGNode(GetNextIdx(), type);
                                pNewNode->m_shocks.AddTail(GetJointPoint(u, GetFirstChild(u)));
                                pNewNode->ComputeDerivedValues();
                                
                                bool localMinOrMax = false;
                                
                                //Before inserting the node, in case of a type 4,
                                // we must make sure it's really a local maximum pt
                                if (type == SHOCK_4)
                                {                               
                                        if (pNewNode->m_shocks > GetSGNode(u)->m_shocks)
                                        {
                                                localMinOrMax = true;
                                                
                                                forall_adj_nodes(v, u)
                                                {
                                                        if (!(pNewNode->m_shocks > GetSGNode(v)->m_shocks))
                                                        {
                                                                localMinOrMax = false;          
                                                                break;
                                                        }
                                                }
                                        }       
                                }
                                else if (type == SHOCK_2)
                                {
                                        if (GetSGNode(u)->m_shocks > pNewNode->m_shocks && GetSGNode(GetFirstChild(u))->m_shocks > pNewNode->m_shocks)
                                                localMinOrMax = true;
                                }
                                
                                if (localMinOrMax) // We have either a type 2 or a type 4
                                {                                                               
                                v = NewNode(pNewNode);
                                
                                forall_adj_edges(e, u)
                                        move_edge(e, v, target(e));
                                                                                        
                                NewEdge(u, v); // keep the old order u -> new (v) -> u's children                                       
                        }
                        else
                                delete pNewNode;
                        }
                } 
        }
}

void ShockGraph::ComputeNodeRole(leda_node v)
{
        NODE_ROLE role;
        int nType = NodeType(v);
        
        if (nType == ROOT)
                role = BIRTH;
        else if (nType == SHOCK_2 || (nType == SHOCK_3 && outdeg(v) == 0))
                role = UNION;
        else if (outdeg(v) == 0)
                role = DEATH;
        else
                role = PROTRUSION;
        
        UnsafeGetSGNode(v)->SetNodeRole(role);
}

sg::DDSkeleton* ShockGraph::ComputeDDSkeleton(const char* szPPMFileName)
{
        using namespace sg;
        
        cout << endl << "Reading " << szPPMFileName << "... " << flush;
        
        int i, j, rows, cols, format;
        pixval maxval;
        
        FILE* fp = fopen(szPPMFileName, "r");
        
        if (fp == NULL)
    {
                cerr << "Can't open " << szPPMFileName << endl;
                return NULL;
    }
        
        ppm_readppminit(fp, &cols, &rows, &maxval, &format);
        fclose(fp);
        
        //  cout << "Getting array...";
        char *a = get_arr(&rows, &cols, szPPMFileName);  
        
        //  cout << "done." << endl;
        // first, we must make a Shape
        sg::SimpleShapeMaker ssm(cols, rows);
        
        // update the SimpleShape
        for (i=0; i < rows; i++)
                for (j=0; j < cols; j++)
                        ssm(j,i) = a[i*cols + j];
                
        // free the input array 
        free(a);
        
        // the shape is being made here
        SimpleShape ss = *((SimpleShape*)ssm.getShape());
        
        //double xmin,xmax,ymin,ymax; //use the member variables instead
        ss.getBounds(&xmin,&xmax, &ymin,&ymax);
        
        //std::cerr << "(xmin,ymin): ("<<xmin<<", "<<ymin<<")\n";
        //std::cerr << "(xmax,ymax): ("<<xmax<<", "<<ymax<<")\n";
        
        // we can get the contour of a SimpleShape as is consists of a
        // single Curve (actually a DiscreteSegCurve) 
        
        //DiscreteSegCurve dsc = *(DiscreteSegCurve *)ssm.getContour();
        //std::cerr << dsc << "\n";
        
        // now, we can create a DistanceTransform object
        DistanceTransform dt(&ss);
        
        // and a DivergenceMap
        DivergenceMap dm(dt);

        cout << "making skeleton... " << flush;
        
        // to supply to the SkeletonMaker
        // cout << "Making the skeleton...";
        DivergenceSkeletonMaker dsm(dm);
        double thresh=-2.0, mag=1.0;
        
        ASSERT(compParams.dSkTreshold >= 0);
        
        return dsm.getDiscreteDivergenceSkeleton(1.0 / mag, -compParams.dSkTreshold);
}

void ShockGraph::DetectLigatureNodes(sg::DDSkeleton* sk)
{
        using namespace sg;

        DDSEdgeVect& edges = sk->getEdges();
        DDSEdgeVect::iterator I;
        DDSEdge* e;

        vector<Curve*>* v = sk->getShape()->getCurves();
        Curve* c = *v->begin();
        DiscreteSegCurve dsc = *(DiscreteSegCurve *)c;

        PRINT(cout, v->size());
        
        for(I = edges.begin(); I != edges.end(); I++)
        {
                e = *I;
                
                Point pRt1 = dsc.atT(e->rightT1());
                Point pRt2 = dsc.atT(e->rightT2());
                Point pLt1 = dsc.atT(e->leftT1());
                Point pLt2 = dsc.atT(e->leftT2());

                PRINT_OPEN(cout, e->hasRight());
                PRINT(cout, pRt1);
                PRINT(cout, pRt2);
                PRINTN(cout, e->hasLeft());
                PRINT(cout, pLt1);
                PRINT(cout, pLt2);
                PRINTN(cout, e->getFluxPoints().front().val);
                PRINT_CLOSE(cout, e->getFluxPoints().back().val);
        }
}

bool ShockGraph::ComputeSGFromDDSkeleton(sg::DDSkeleton* sk)
{
        using namespace sg;

        cerr << "labeling skeleton... " << flush;       

        //DetectLigatureNodes(sk);
        
        SGNode* pNode = new SGNode(GetNextIdx(), SHOCK_base);
        leda_node u = NewNode(pNode);
        DDSEdge* se = *sk->getEdges().begin();
        
        SGNode* pFirstNode = CopyNodeInfo(se, u, pNode);
        
        ConnectNodes(u, pNode, se, se->getN1());
        ConnectNodes(u, pFirstNode, se, se->getN2());
        
        Insert4sAnd2s();
        //ShowGraph(*this);
        
        bool bContinue;
        leda_edge e;
        int nNoInfLoop = 0, i;
        
        // Collect the edges that need to be reverted according to time of formation.
        // Edges should be directed from last-to-form to earliest-to-form.
        do
        {
                ASSERT(nNoInfLoop++ < 1000);
                bContinue = false;
                
                forall_edges(e, *this)
                {
                        const ShockBranch& tgt = GetSGNode(target(e))->m_shocks;
                        const ShockBranch& src = GetSGNode(source(e))->m_shocks;
                                
                        if (tgt > src)
                        {
                                RevertEdge(e);
                                bContinue = true;
                        }
                }                       
        } while (bContinue);
        
        // Nodes that do not have a parent should be children of the root node.
        leda_node root = NewNode(new SGNode("#", ROOT));
        leda_node v;
        int nodeType;
        
        // BTW, type 3 nodes with 2 type 1 parent/children should
        // be relabeled as 2/4 respectivelly.
        forall_nodes(v, *this)
        {
                if (v != root && indeg(v) == 0)
                        NewEdge(root, v);
                        
                nodeType = NodeType(v);
                
                if (nodeType == SHOCK_3) // May need to be relabeled.
                        Relabel3As2or4(v);
                
                // Make sure the grammar is satisfied
                ASSERT (nodeType == ROOT ||
                        nodeType == SHOCK_1 ||
                        nodeType == SHOCK_3 ||
                        (nodeType == SHOCK_2 && indeg(v) == 2 && outdeg(v) == 0) ||
                        (nodeType == SHOCK_4 && indeg(v) == 1 && NodeType(GetFirstParent(v)) == ROOT));

                ComputeNodeRole(v);
        }

        ComputeDerivedValues();

        cerr << "(" << GetObjName() << ", " << GetViewNumber() << ") done!" << flush;
        
        return true;
}

bool ShockGraph::ComputeFromPPMFile2(const char* szFileName, const SGCP& sgparams)
{
        Clear(); //make sure the current shock graph is empty

        compParams = sgparams;
        SetDAGLbl(szFileName);
        
        sg::DDSkeleton* pDDS = ComputeDDSkeleton(szFileName);
        
        if (!pDDS)
                return false;

        m_pSkeleton = new McGillSkeleton(pDDS);

        return ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
}

bool ShockGraph::RecomputeFromSkeleton(const SGCP& sgparams)
{
        // We must empty the graph but without deleting the skeleton, label, mins and maxs.
        String strLbl = GetDAGLbl();
        DAG::Clear();
        nLastIndexUsed = 0;

        compParams = sgparams;
        SetDAGLbl(strLbl);

        return ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
}

void SmoothFluxPoints(sg::FluxPointList& fpl, const double& dSigma, const double& dSigmaRange)
{
        int i, nLen = fpl.size();

        double* in = new double[nLen];
        double* out = new double[nLen];

        for (i = 0; i < nLen; i++)
                in[i] = fpl[i].dist;

        GaussianSmooth(in, out, nLen, dSigma, dSigmaRange);

        for (i = 0; i < nLen; i++)
                fpl[i].dist = out[i];

        delete[] in;
        delete[] out;
}

void GaussianSmooth(double *input, double *output, int length, 
                                        const double& sigma, const double& sigma_range)
{
        static double* kernel = NULL;
        static double prev_sigma, prev_sigma_range;
        static double kernelsum;
        static int w2;
        int i, x;

        if (kernel && (sigma != prev_sigma || sigma_range != prev_sigma_range))
        {
                delete[] kernel;
                kernel = NULL;
        }

        // first compute the kernel
        if (!kernel)
        {
                w2 = (int)rint(sigma_range * sigma);
                if (w2 < 1) w2 = 1;
                int w = (2 * w2) + 1;
                
                kernel = new double[w];
                kernelsum = 0.0;
                prev_sigma = sigma;
                prev_sigma_range = sigma_range;
                
                for (double* kp = kernel, x = -w2; x <= w2; kp++,x++)
                {
                        *kp = (1 / (sqrt(2 * M_PI) * sigma)) * exp((double)(x * x)/(-2 * sigma * sigma));
                        kernelsum += *kp;
                }
        }
        
        PRINT(cerr, kernelsum);

        double sum;

        // now perform the convolution
        for (x = 0; x < length; x++)
        {
                sum = 0.0;

                for (i = -w2; i <= w2; i++)
                {
                        if ((x + i) < 0){
                                sum += kernel[i + w2] * input[0]; /* extend left side */
                        }
                        else if ((x + i) >= length){
                                sum += kernel[i + w2] * input[length - 1]; /* extend right side */
                        }
                        else{
                                sum += kernel[i + w2] * input[x + i];
                        }
                }

                output[x] = sum / kernelsum;
        }
}

int ShockGraph::GetBranchDir(leda_node u, leda_node wrtV) const
{
        const ShockInfo& pt = GetJointPoint(u, wrtV);
        
        return GetSGNode(u)->GetBranchDir(pt == GetSGNode(u)->m_shocks[0] ?
                ShockBranch::FORWARD:ShockBranch::BACKWARD);
}

const ShockInfo& ShockGraph::GetJointPoint(leda_node u, leda_node v) const
{
        const SGNode* pNodeU = GetSGNode(u);
        const SGNode* pNodeV = GetSGNode(v);
        
        const ShockInfo& uStart = pNodeU->m_shocks[0];
        const ShockInfo& uEnd   = pNodeU->m_shocks.GetTail();

        const ShockInfo& vStart = pNodeV->m_shocks[0];
        const ShockInfo& vEnd   = pNodeV->m_shocks.GetTail();
        
        if (uStart == vStart || uStart == vEnd)
                return uStart;
        else if (uEnd == vStart || uEnd == vEnd)
                return uEnd;
        
        ASSERT(false);
        return vStart;
}

---