typedef leda_list<leda_edge> EdgeList;
typedef leda_list<leda_node> NodeList;
typedef leda_d_array<leda_node, leda_node> NodeMap;
typedef leda_d_array<leda_edge, leda_node> EdgeNodeMap;
typedef leda_d_array<leda_node, leda_list_item> NodeListItemMap;
typedef leda_d_array<leda_node, int> NodeIndexMap;

typedef leda_list<DAGNodePtr> DAGNodeList;
typedef leda_d_array<DAGNodePtr, DAGNodePtr> DAGNodeMap;


/*
 * For a given node, the	 following method counts the number of its siblings,
 * parents and children, if any.
 * This number is necassary to create an input file to
 * feed the earth mover's distance code.
 */
int GestureGraph::Get_number_of_relations(int node){

 int counter = 0; //The first number accounts for the node itself.
 int nVertices = hierarchyLevels.GetSize();

//cout << "no_vertices from size :" << hierarchyLevels.GetSize() << endl;

 for(int i=0; i < nVertices; i++){

        if (hierarchicalMatrix[node][i]!=0) //count all the relations.
                 counter++;

 }

 return counter;

}


int GestureGraph::Get_nth_relation(int n, int node){

  int no = 0;
  int nVertices = hierarchyLevels.GetSize();

        for (int i=0; i<nVertices; i++){
                if (hierarchicalMatrix[node][i] != 0)
                        no++;

                if (no==n)
                        return i;
        }

 return 0;
}

/*
 * Computing the Euclidean distance between 2 nodes.
 */

float GestureGraph::Get_distance(int node1, int node2){

 const DAGNode *dag_ptr1 = GetNode(GetNode(node1));
 const GGNode *p1 = dynamic_cast<const GGNode*>(dag_ptr1);

 const DAGNode *dag_ptr2 = GetNode(GetNode(node2));
 const GGNode *p2 = dynamic_cast<const GGNode*>(dag_ptr2);

 int x1 = p1->m_nXPos;
 int y1 = p1->m_nYPos;
 int x2 = p2->m_nXPos;
 int y2 = p2->m_nYPos;
 float temp = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2);

 return sqrt(temp);
}


int GestureGraph::Get_xpos(int node){

 const DAGNode *dag_ptr = GetNode(GetNode(node));
 const GGNode *p = dynamic_cast<const GGNode*>(dag_ptr);

 return p->m_nXPos;

}

int GestureGraph::Get_ypos(int node){

 const DAGNode *dag_ptr = GetNode(GetNode(node));
 const GGNode *p = dynamic_cast<const GGNode*>(dag_ptr);

 return p->m_nYPos;

}

void GetPath(int i, int last, int tgt)
{
	COctahedronPoints pts(3);
	int* pNeigPts = pts.GetNeighbours(i - 1, 2);
	int n1 = pNeigPts[0] + 1;
	int n2 = pNeigPts[1] + 1;
	free(pNeigPts);

	int n = n1 == last ? n2:n1;

	if (n == tgt)
		return;
	else
	{
		cout << endl << n;
		GetPath(n, i, tgt);
	}
}

SmartArray<int> GetClosestViewsByRange(int i, int res, const double& range)
{
	COctahedronPoints pts(res);
	int n;

	//int* pNeigPts = pts.GetNeighbours(i - 1, range, n);
	n = 126;
	int* pNeigPts = pts.GetNeighbours(i - 1, n);
	SmartArray<int> views(n);
	
	for (int j = 0; j < n; j++)
		views[j] = pNeigPts[j] + 1;

	free(pNeigPts);

	return views;
}

SmartArray<int> GetSymmetricViewMap()
{
	const int nViewsPerObj = 128;
	SmartArray<int> symmap, views;
	int viewmap[nViewsPerObj], i, j, k;

	symmap.ReSize(nViewsPerObj + 1, true); // view 0 isn't used

	for (i = 1; i <= nViewsPerObj; i++)
	{
		views = GetClosestViewsByRange(i, 3, 178); // 178 = a bit less than 180.

		// Turn off all views
		memset(viewmap, 0, nViewsPerObj * sizeof(int));

		// Turn on views that are within 178 deg.
		for (j = 0; j < views.GetSize(); j++)
			viewmap[views[j] - 1] = 1;

		// Turn on the target view
		viewmap[i - 1] = 1;

		// There should be only one view off = 0
		for (k = 0; k < nViewsPerObj; k++)
			if (viewmap[k] == 0 && symmap[i] == 0)
				symmap[i] = k + 1;
			else if (viewmap[k] == 0 && symmap[i] != 0)
			{
				ShowError("Multiple symmetric views.");
				ASSERT(false);
			}
			else if (k == nViewsPerObj - 1 && symmap[i] == 0 && viewmap[k] != 0)
			{
				ShowError("One view is unmapped.");
				ASSERT(false);
			}
	}

	return symmap;
}



	/*//Check if it is type 4
	if (outdeg(v) == 2 && indeg(v) == 1 && parentType == ROOT)
	{
		int d1 = GetBranchDir(GetFirstChild(v), v);
		int d2 = GetBranchDir(GetSecondChild(v), v);
		
		// If children are decreasing must relabel
		if (d1 == -1 && d2 == -1)
			UnsafeGetSGNode(v)->SetNodeType(SHOCK_4);
	}
	//Check if it is type 2
	else if(outdeg(v) == 0 && indeg(v) == 2)
	{
		int d1 = GetBranchDir(GetFirstParent(v), v);
		int d2 = GetBranchDir(GetSecondParent(v), v);
		
		// If parents are increasing must relabel
		if (d1 == 1 && d2 == 1)
			UnsafeGetSGNode(v)->SetNodeType(SHOCK_2);
		else // it's a 3 with 1's decresing. change that!!!
		{
			WARNING(d1 == -1 || d2 == -1, "Should not occurr in this context");
			forall_in_edges(e, v)
				if (GetBranchDir(source(e), v) == -1)
					rev_edge(e);
		}
	}*/


	/*DBG_VAL(pNode1->GetSegments().GetHead().p0)
	DBG_VAL(g1Pts.GetHead())
	DBG_VAL(pNode1->GetSegments().GetTail().p1)
	DBG_VAL(g1Pts.GetTail())
	
	DBG_VAL(pNode2->GetSegments().GetHead().p0)
	DBG_VAL(g2Pts.GetHead())
	DBG_VAL(pNode2->GetSegments().GetTail().p1)
	DBG_VAL(g2Pts.GetTail())*/

int GetLeaves(TREE_NODE root, leda_graph& g, leda_list<TREE_NODE>& leaves)
{
	int nMaxHeight = root.height;

	if (g.outdeg(root.node) > 0)
	{
		leda_node v, w = root.node;
		int h;

		root.height++;

		forall_adj_nodes(v, w)
		{
			// If loopy node, hide it so as to process it last
			if (indeg(v) > 1)
				g.hide_node(v);
			else
			{
				root.node = v;
    			h = GetLeaves(root, g, leaves);

    			if (h > nMaxHeight)
    				nMaxHeight = h;
			}
		}
		
		// Let's now process loopy nodes, which are hidden
		list_item it;
		leda_list<leda_node> hnlist = g.hidden_nodes();
		leda_edge e;
		
		forall_items(it, hnlist)
		{
			// Add the loopy node. It will be last
			// among its siblings
			v = hnlist[it];	
			root.node = v;	
			h = GetLeaves(root, g, leaves);

    		if (h > nMaxHeight)
    			nMaxHeight = h;
    			
    		// Now that we porcess all the children of w,
    		// we continue with the other parents of the
    		// loopy nodes. We process these parents and

    		// the hide them, so that they don't get process twice.
    			
    		forall_in_edges(e, v)
    		{
				if (g.source(e) != w)
				{
					root.node = g.source(e);
					root.height = GetNodeHeight(g, root.node);
					
					h = GetLeaves(root, g, leaves);
					g.hide_node(root.node);
					
					if (h > nMaxHeight)
    					nMaxHeight = h;
				}
			}
					
			g.restore_node(v);
		}
		
		// If we are dealing with the root, it is time to unhide
		// all th eparents of loopy nodes, becaus ewe are done.
		if (g.indeg(w) == 0)
			g.restore_all_edges();
	}
	else
		leaves.append(root);

	return nMaxHeight;
}



void ShockGraph::RevertEdge(leda_edge pxEdge)
{
	leda_node parent = source(pxEdge);
	leda_node x = target(pxEdge);
	leda_edge e;
	SmartArray<leda_edge> edgesToDelete(1, SA_GROWING_FACTOR);
	
	edgesToDelete.AddTail(pxEdge); // pxEdge will be del at the end
	
	// if x not type 3 or has no children
	if (NodeType(x) != SHOCK_3 || outdeg(x) == 0)
	{
		// 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 (e != pxEdge && GetJointPoint(parent, target(e)) == inPt)
			{
				NewEdge(x, target(e));
				edgesToDelete.AddTail(e);
			}
		}
		
		NewEdge(x, parent); // creates x->p edge
	}
	else // There are nodes coming from different sides
	{
		ASSERT(NodeType(x) == SHOCK_3);
		
		leda_node leftChild = NewNode(new SGNode(*GetSGNode(x)));
		leda_node rightChild = NewNode(new SGNode(*GetSGNode(x)));
		
		// 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 (e != pxEdge && GetJointPoint(parent, target(e)) == inPt)
			{
				NewEdge(leftChild, target(e));
				edgesToDelete.AddTail(e);
			}
		}
		
		NewEdge(leftChild, parent); // creates x's left->p edge
		
		// Move all children of x to the right side of x
		forall_adj_edges(e, x)
		{
			NewEdge(rightChild, target(e));
			edgesToDelete.AddTail(e);
		}
		
		// Link x to its left and right sides
		NewEdge(x, rightChild);
		NewEdge(x, leftChild);
	}
	
	for (int i = 0; i < edgesToDelete.GetSize(); i++)
		//DelEdge(edgesToDelete[i]);
		hide_edge(edgesToDelete[i]);
}
	
	
	
	
	double tau = pts.AvgRadius() * compParams.dSigma; //compParams.dSigma
	int winSize = int(ceil(pow(pts.Length(), 3.0) / (compParams.dSigmaRange * pts.Area())));

	if (winSize < 1) winSize = 1;

	//ApproxPoly poly(tau, compParams.dDerivCutoff);
	//poly.TwoPassFit(data, nSize);
	//poly.Fit(data, nSize);
	//poly.Fit(data, nSize, winSize, compParams.dDerivCutoff, compParams.dSigma);
	//poly.Fit(data, nSize, int(ceil(nSize / 12)));


DDSNode* McGillSkeleton::ReadNodesAndEdges(istream& is, DDSkeleton& dds)
{
	int fp_count, n1_adj_edges_count, i, j, tgtIdx;
	double lt1, lt2, rt1, rt2;
	FluxPoint fp;
	DDSNode* n2;

	// read n1
	is.read((char*) &fp, sizeof(fp));
	is.read((char*) &n1_adj_edges_count, sizeof(n1_adj_edges_count));

	DDSNode* n1 =  new DDSNode(fp);
	dds.addNode(n1);

	for (i = 0; i < n1_adj_edges_count; i++)
	{
	        is.read((char*) &tgtIdx, sizeof(tgtIdx));

		// Read target node
		n2 = ReadNodesAndEdges(is, dds);

		//Read edge
		is.read((char*) &lt1, sizeof(lt1));
		is.read((char*) &lt2, sizeof(lt2));
		is.read((char*) &rt1, sizeof(rt1));
		is.read((char*) &rt2, sizeof(rt2));
		is.read((char*) &fp_count, sizeof(fp_count));
		ASSERT(fp_count > 0);
		FluxPointList fpl;

		for (j = 0; j < fp_count; j++) 
		{
			is.read((char*) &fp, sizeof(fp));
			fpl.push_back(fp);
		}

		//Add edge
		DDSEdge* e = new DDSEdge(&dds, fpl, tgtIdx == 2 ? n1:n2, tgtIdx == 2 ? n2:n1);
		if (lt1 != -1 || lt2 != -1) e->setLeftSegment(lt1, lt2);
		if (rt1 != -1 || rt2 != -1) e->setRightSegment(rt1, rt2);
		dds.addEdge(e);
	}

	return n1;

  leda_d_array<long, leda_node> map;


}

void McGillSkeleton::WriteNodesAndEdges(ostream& os, DDSNode* src, DDSEdge* pInEdge)
{
	FluxPoint fp;
	DDSEdgeVect::iterator I;
	DDSEdgeVect edges = src->getEdges();
	int src_adj_edges_count = (pInEdge) ? edges.size() - 1:edges.size(); // Do not save in edge
	DDSEdge* e;
	int tgtIdx;

	// Write node 1
	os.write((char*) &src->fp, sizeof(src->fp));
	os.write((char*) &src_adj_edges_count, sizeof(src_adj_edges_count));
    
	for(I = edges.begin(); I != edges.end(); I++)
	{
		e = (*I);

		// Write only adj edges
		if (pInEdge != e)
		{
		        ASSERT(src == e->getN1() || src == e->getN2());

		        // Keep the adjacency indices in the edge
		        tgtIdx = (src == e->getN1()) ? 2:1;
			os.write((char*) &tgtIdx, sizeof(tgtIdx));

			//Write target node's adj edges
			WriteNodesAndEdges(os, (tgtIdx == 1) ? e->getN1():e->getN2(), e);

			//Write edge
			FluxPointList& fpl = e->getFluxPoints();
			int fp_count = fpl.size();
			ASSERT(fp_count > 0);

			os.write((char*) &e->left_t1, sizeof(e->left_t1));
			os.write((char*) &e->left_t2, sizeof(e->left_t2));
			os.write((char*) &e->right_t1, sizeof(e->right_t1));
			os.write((char*) &e->right_t2, sizeof(e->right_t2));
			os.write((char*) &fp_count, sizeof(fp_count));

			FluxPointList::iterator II;

			for(II = fpl.begin(); II != fpl.end(); II++)
			{
				fp = *II;
				os.write((char*) &fp, sizeof(fp));
			}		
		}
	}
}








/*!
@brief Finds the branch with the largest average radius.
*/
sg::DDSEdge* ShockGraph::FindLatestNodeToForm(sg::DDSkeleton* sk) const
{
	using namespace sg;
	
	DDSEdgeVect& edges = sk->getEdges();
	DDSEdgeVect::iterator I;
	DDSEdge* e = *edges.begin();
	double dMaxAvgRadius = -1, dAvgRadius;
	int i;
	
	for(I = edges.begin(); I != edges.end(); I++)
	{
		const sg::FluxPointList& fpl = (*I)->getFluxPoints();
		
		for (i = 0, dAvgRadius = 0.0; i < fpl.size(); i++)
			dAvgRadius += fabs(fpl[i].dist);
		
		if (dAvgRadius > 0.0)
			dAvgRadius = dAvgRadius / fpl.size();
		
		// Save the gratest min point among all branches
		if (dAvgRadius > dMaxAvgRadius)
		{
			dMaxAvgRadius = dAvgRadius;
			e = *I;
		}
	}
	
	return e;
}

double SweepContrast(const SGNode *node)
{
	
	double min_rad = 999999;
	double max_rad = 0;
	double curr_rad = 0;
	double sweep_contrast;
	
	int shockCount = node->GetShockCount();  
	
	for (int j = 0; j < shockCount; j++)
    {
		double curr_rad = node->m_shocks[j].radius;	    
		if (curr_rad > max_rad)
			max_rad = curr_rad;
		if (curr_rad < min_rad)
			min_rad = curr_rad;
    }
	sweep_contrast = ((double) max_rad - (double) min_rad)/(double) max_rad;
	return(sweep_contrast);
}


/* This function labels the node as being a 1,2,3, or 4 based on first assuming
that the sweep function is second order.  We will find the max and minimum radius
along the skeleton, and look at their deviation and their relative positions,
and based on simple ordering rules, we'll classify the skeleton. 
*/

int ShockGraph::ComputeNodeLabel(const SGNode *node) const
{
	double sweep_contrast_thresh = .20;
	double min_rad = 999999;
	double max_rad = 0;
	double curr_rad = 0;
	double sweep_contrast = SweepContrast(node);
	double avg_radius = 0;
	
	int shockCount = node->GetShockCount();  
	int node_type;
	
	double first_rad = node->m_shocks[0].radius;
	double last_rad = node->m_shocks[shockCount - 1].radius;
	
	for (int j = 0; j < shockCount; j++)
	{
		double curr_rad = node->m_shocks[j].radius;	    
		if (curr_rad > max_rad)
			max_rad = curr_rad;
		if (curr_rad < min_rad)
			min_rad = curr_rad;
		avg_radius = avg_radius + curr_rad;
	}
	avg_radius = avg_radius / (double) shockCount;
	
	// If both end points have significantly larger radius than min, then type 2
	if ((min_rad < first_rad) && (min_rad < last_rad) &&
		((first_rad - min_rad)/first_rad >= sweep_contrast_thresh) &&
		((last_rad - min_rad)/last_rad >= sweep_contrast_thresh))
	{
		node_type = SHOCK_2;
	}
	// Else if both endpoints have significantly less radius than max, then type 4
	else if ((max_rad > first_rad) && (max_rad > last_rad) &&
		((max_rad - first_rad)/max_rad >= sweep_contrast_thresh) &&
		((max_rad - last_rad)/max_rad >= sweep_contrast_thresh))
	{
		node_type = SHOCK_4;
	}
	// If average radius is near the max or the min, then type 3.
	else if (((max_rad - avg_radius)/max_rad <= sweep_contrast_thresh) ||
		((avg_radius - min_rad)/max_rad <= sweep_contrast_thresh))
	{
		node_type = SHOCK_3;
	}
	else
	{
		node_type = SHOCK_1;
	}
	return node_type;
}

NODE_ROLE ShockGraph::ComputeNodeRole(leda_node v) const
{
	int nType = NodeType(v);
	
	// Is birth. Rule 1
	if (nType == ROOT)
		return BIRTH;
	
	if (GetSGNode(v)->GetShockCount() == 1 && nType != SHOCK_2 && nType != SHOCK_4)
		return DUMMY;
	
	int nParentType = NodeType(source(first_in_edge(v)));
	
	// Is union
	if (indeg(v) == 2 && (nType == SHOCK_2 || nType == SHOCK_3))
	{
		// Rule 5 and 6a
		//ASSERT(nType == SHOCK_2 || nType == SHOCK_3);
		return UNION;
	}
	else if (nType == SHOCK_3 && nParentType == SHOCK_1)
		// Rule 6b
		return UNION;
	
	// Is death. Rules 7, 8, 9 and 10
	if (outdeg(v) == 0)
		return DEATH;
	
	

	int nFirstOrderSeg = 0;
	leda_node w;
	
	forall_adj_nodes(w, v)

		if (NodeType(w) == SHOCK_1)
			nFirstOrderSeg++;
		
		bool bAllOnes = nFirstOrderSeg == outdeg(v);
		
		// Is protrusion
		
		// Rule 4
		if (nType == SHOCK_1 && nFirstOrderSeg >= 2 && bAllOnes)
			return PROTRUSION;
		
		// Rule 2
		else if (nType == SHOCK_4 && bAllOnes)
			return PROTRUSION;
		
		// Rule 3b
		else if (nType == SHOCK_3 && nParentType == ROOT && bAllOnes)
			return PROTRUSION;
		
		else if (nType == SHOCK_3) // and more
			return PROTRUSION;
		
			/*cerr << endl << GetDAGLbl() << " in node " << GetNodeLbl(v) << 
			"\nUnknown role: " << 
			"\n\tNode type: " << nType <<
			"\n\tParent type: " << nParentType <<
			"\n\tnFirstOrderSeg: " << nFirstOrderSeg <<
			"\n\tbAllOnes: " << bAllOnes << endl << flush;
			
			 
		DBG_MSG("Unknown SG node role!!!");*/
		
		return UNK_ROLE;
}
void ShockGraph::EliminateSpuriousNodes()
{
	leda_node v;
	bool flag = false;
	
	forall_nodes(v, *this)
	{

		if (NodeType(v) == SHOCK_3 && 
			indeg(v) == 1 && outdeg(v) == 1 && 
			GetSGNode(v)->GetShockCount() == 1)
		{
			leda_edge outEdge = first_adj_edge(v);
			
			if (NodeType(target(outEdge)) == SHOCK_1)
			{
				leda_edge inEdge = first_in_edge(v);
				
				NewEdge(source(inEdge), target(outEdge));
				
				del_edge(inEdge);
				del_edge(outEdge);
				del_node(v);
				
				flag = true;
			}
		}
	}
	
	if (flag)
		cerr << endl << GetDAGLbl() << " has been cleaned up!";
}


bool ShockGraph::ComputeFromPPMFile2(const char* szFileName, 

		                     double cutoff, double sigma, double range)
{
  using namespace sg;
  
  // do the skeleton work, Alex's code
  //  cout << endl << "Generating the skeleton:" << endl;
  
  int i, j, rows, cols, format;
  pixval maxval;
  
  FILE* fp = fopen(szFileName, "r");
  
  if (fp == NULL)
    {
      cerr << "Can't open " << szFileName << endl;
      return false;
    }
  
  ppm_readppminit(fp, &cols, &rows, &maxval, &format);
  fclose(fp);
  
  //  cout << "Getting array...";
  char *a = get_arr(&rows, &cols, szFileName);  
  //  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;
  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);
  
  // to supply to the SkeletonMaker
  //  cout << "Making the skeleton...";
  DivergenceSkeletonMaker dsm(dm);
  double thresh=-2.0, mag=1.0;
  DDSkeleton *skeleton = dsm.getDiscreteDivergenceSkeleton(1.0/mag, thresh); 
  //  cout << "done." << endl;  

  /*************************************NEW CODE********************************/
  
  // analyse the points, create the LEDA graph (SG)
  //  cout << endl << "Generating the shock graph:" << endl;
  
  Clear(); //make sure the graph is empty
  
  // stores the shock points
  std::vector<DiscreteDivergenceSkeletonNode*> nodes;
  std::vector<DiscreteDivergenceSkeletonEdge*> edges;
  
  nodes = skeleton->getNodes();	// get all nodes
  edges = skeleton->getEdges();
  
  //DiscreteDivergenceSkeletonEdge* SKparent;
  SGNode* pNode = NULL;
  //	SGNode* SGparent = NULL;
  leda_d_array<long, leda_node> map;
  char szNodeLbl[10];
  
  // make some LEDA nodes
  for ( int i=0; i<edges.size(); i++ )
    {
      long pointerID = (long)(edges[i]);
      
      sprintf(szNodeLbl, "%d", i);
      pNode = new SGNode(szNodeLbl, 0);
      
      FluxPointList fpl = (edges[i])->getFluxPoints();
      int num_pts = fpl.size();
      
      pNode->m_shocks.ReSize((num_pts>0) ? num_pts:0);
      for ( int j=0; j<num_pts; j++ )
	{
	  FluxPoint pt = fpl[j];
	  ShockInfo& si = pNode->m_shocks[j];
	  si.xcoord = pt.p.x;
	  si.ycoord = pt.p.y;
	  si.radius = -1.0 * pt.dist; // guessing here
	  si.speed  = pt.val;  // guessing here too
	  si.color = 0;
	  
	  // Set bounding rectangle for drawing purposes
	  /*if (pt.p.x - si.radius < this->xmin) this->xmin = pt.p.x - si.radius;
	  if (pt.p.x + si.radius > this->xmax) this->xmax = pt.p.x + si.radius;
	  if (pt.p.y - si.radius < this->ymin) this->ymin = pt.p.y - si.radius;
	  if (pt.p.y + si.radius > this->ymax) this->ymax = pt.p.y + si.radius;*/

	  if (pt.p.x < this->xmin) this->xmin = pt.p.x;
	  if (pt.p.x > this->xmax) this->xmax = pt.p.x;
	  if (pt.p.y < this->ymin) this->ymin = pt.p.y;
	  if (pt.p.y > this->ymax) this->ymax = pt.p.y;
	}
		
      pNode->ComputeDerivedValues();
      pNode->SetNodeRole(NODE_ROLE(0));

      map[pointerID] = NewNode(pNode); // make a map entry
      

      //      cout << "Node #" << i << " has " << pNode->GetShockCount() << " points and average radius " 
      //	   << pNode->m_shocks.AvgRadius() << ". Type is " << pNode->GetType() << endl;
    }
  
  //  cout << "Drawing the edges (v6)..." << endl;
  DiscreteDivergenceSkeletonEdge* root_edge;  
  sprintf(szNodeLbl, "%c", '#');
  SGNode* hash = new SGNode(szNodeLbl, 0);
  hash->SetNodeRole(NODE_ROLE(BIRTH));
  hash->SetNodeType(ROOT);
  map[(long)hash] = NewNode(hash);
  
  while (( root_edge = GetGreatestUntouched(map,edges) ) != NULL )
    {
      // split if necessary, otherwise returns itput node
      pNode = (SGNode*)GetSGNode(map[(long)root_edge]);
      SGNode* pNode2 = SplitSGNode(pNode, map);

      // draw edges from # to each root edge
      NewEdge( map[(long)hash], map[(long)root_edge], 1.0 );
      //      cout << "Edge drawn from # to avg rad " << GetSGNode(map[(long)root_edge])->m_shocks.AvgRadius();
      if ( pNode != pNode2 ) // if there's a split we need a second root
	{
	  NewEdge( map[(long)hash], map[(long)pNode2], 1.0 );
	  //	  cout << " and to it's mirror side.";
	}
      //      cout << endl;
      // call DrawEdges to recursively traverse the skeleton and create edges
      if ( root_edge->getN1() != NULL ) DrawEdges(root_edge, root_edge->getN1(), map, pNode);
      if ( root_edge->getN2() != NULL ) DrawEdges(root_edge, root_edge->getN2(), map, pNode2);
    }
  ComputeDerivedValues();
  SetDAGLbl(szFileName);
  return true;
}

SGNode* ShockGraph::SplitSGNode(SGNode* pNode, leda_d_array<long, leda_node>& map)
{
  int type = ComputeNodeLabel(pNode);
  if ( type == SHOCK_3 || type == SHOCK_2 )
    {
      SGNode* newNode = new SGNode(*pNode);
      map[(long)newNode] = NewNode(newNode);
      // should add labels 'a' and 'b' to each side
      return newNode;
    }
  else
    return pNode;
}

sg::DiscreteDivergenceSkeletonEdge* ShockGraph::GetGreatestUntouched(leda_d_array<long, leda_node>& map, 
	std::vector<sg::DiscreteDivergenceSkeletonEdge*> edges)
{
  using namespace sg;
  DiscreteDivergenceSkeletonEdge* greatest_edge = NULL;
  for ( int i=0; i<edges.size(); i++ )
    {
      const SGNode* node = GetSGNode(map[(long)edges[i]]);
      if ( node->GetType() == 0 ) 
	{
	  if ( greatest_edge == NULL || node->m_shocks.AvgRadius() > GetSGNode(map[(long)greatest_edge])->m_shocks.AvgRadius() )
	    {
	      // select new greatest edge
	      greatest_edge = edges[i];
	    }
	}
    }
  return greatest_edge;
}

/*!
	@brief Recursive edge maker
*/
void ShockGraph::DrawEdges(sg::DiscreteDivergenceSkeletonEdge* SKparent,
	sg::DiscreteDivergenceSkeletonNode* node, leda_d_array<long, leda_node>& map, SGNode* SGparent)
{
  using namespace sg;
  std::vector<DiscreteDivergenceSkeletonEdge*> l_edges;	// stores local edges
  long pointerID;
  
  // in the case of 3's and 2's the SKparent will return a different pointer than
  // what is expected.  This is because of the duplicate entries for 3's and 2's
  // we want to use the correct key, in this case SGparent itself, not SKparent
  SGparent->SetNodeType(ComputeNodeLabel(SGparent)); // set the type
  if ( SGparent == GetSGNode(map[(long)SKparent]) )
    pointerID = (long)SKparent;
  else
    pointerID = (long)SGparent;
  
  if ( node->getEdges().size() > 1 )
    {
      l_edges = node->getEdges();
      for ( int j=0; j<l_edges.size(); j++ )
	{ 
	  // compare radii (do not recurse if a larger radius is found
	  SGNode* SGchild = (SGNode*)GetSGNode(map[(long)l_edges[j]]);
	  if ( SKparent != l_edges[j] && SGparent->m_shocks.AvgRadius() > SGchild->m_shocks.AvgRadius() )
	    {
	      // draw edges from parent to (each) child, check for splits
	      SGNode* SGchild2 = SplitSGNode(SGchild, map);
	      NewEdge( map[pointerID], map[(long)l_edges[j]], 1.0); // draw the edge
	      //	      cout << "Edge drawn from avg rad " << SGparent->m_shocks.AvgRadius()
		  // << " to avg rad " << SGchild->m_shocks.AvgRadius();	      
	      if ( SGchild != SGchild2 ) // if there's a split we need a second root
		{
		  NewEdge( map[pointerID], map[(long)SGchild2], 1.0 );
		  //		  cout << " and to its mirror side";
		}
	      //	      cout << endl;
	      // recursive call
	      if ( SGchild->GetType() == 0 && l_edges[j]->getN1() != NULL && l_edges[j]->getN1() != node )
		DrawEdges(l_edges[j], l_edges[j]->getN1(), map, SGchild);
	      else 
		SGchild->SetNodeType(ComputeNodeLabel(SGchild));
	      // more recursion...
	      if ( SGchild2->GetType() == 0 && l_edges[j]->getN2() != NULL && l_edges[j]->getN2() != node )
		DrawEdges(l_edges[j], l_edges[j]->getN2(), map, SGchild2);
	      else
		SGchild2->SetNodeType(ComputeNodeLabel(SGchild2));
	    }
	}
    }
}


##############################################################################

int IsAncestor(leda_node ancestor, leda_node child, int nMaxDist) const;

int ComputeMatchDeg(leda_node root, const DAG& g, const NODE_MAP& matchMap, 
						 leda_node& bestRoot, int& nMaxDeg) const;

leda_node GetBestMatchedSubtree(const DAG& g, const NODE_MAP& matchMap) const;

int DAG::IsAncestor(leda_node ancestor, leda_node child, int nMaxDist) const
{
	leda_edge e;

	if (nMaxDist == 0)
		return 0;

	forall_in_edges(e, child)
	{
		parent = source(e);

		if (parent == ancestor)
			return nMaxDist;
		else if (n = IsAncestor(ancestor, parent, nMaxDist - 1))
			return n;
	}

	return 0;
}

int DAG::ComputeMatchDeg(leda_node root, const DAG& g, const NODE_MAP& matchMap, 
						 leda_node& bestRoot, int& nMaxDeg) const
{
	int n = 0, childDeg;
	int nMaxDist = 2;

	forall_adj_nodes(v, root)
	{
		childDeg = ComputeMatchDeg(v);

		if (g.IsAncestor(matchMap[root], matchMap[v], nMaxDist))
			n += 1 + childDeg;
	}

	if (n > nMaxDeg)
	{
		nMaxDeg = n;
		bestRoot = root;
	}

	return n;
}

leda_node DAG::GetBestMatchedSubtree(const DAG& g, const NODE_MAP& matchMap) const
{
	leda_node root = GetFirstRoot();
	leda_node bestRoot;
	int nMaxDeg = 0;

	ComputeMatchDeg(root, g, matchMap, bestRoot, nMaxDeg);

	return bestRoot;
}

double TreeSimilarity(r1, r2)
{
	if (r1&r2 are leaves)
		return NodeSimilarity();
	else
		...
}

double Match()
{
	Construct bipartite graph with links on related nodes.

	BFS()
	{

	}
}