/**------------------------------------------------------------------------
 * @file
 *
 * Copyright 2001 by the University of Toronto (UT).
 *
 *                       All Rights Reserved
 * 
 * Permission to copy and modify this software and its documentation 
 * only for internal use in your organization is hereby granted, 
 * provided that this notice is retained thereon and on all copies.  
 * UT makes no representations as to the suitability and operability 
 * of this software for any purpose.  It is provided "as is" without 
 * express or implied warranty.
 *
 * UT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING
 * ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT
 * SHALL UT BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
 * OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
 * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 * PERFORMANCE OF THIS SOFTWARE.
 *
 * No other rights, including for example, the right to redistribute 
 * this software and its documentation or the right to prepare 
 * derivative works, are granted unless specifically provided in a 
 * separate license agreement.
 *
 * Author: Diego Macrini based on Ali Shokoufandeh's original code
 * 
 * $Log: ShockGraph.cpp,v $
 * Revision 1.4  2002/06/13 21:46:37  dmac
 * Matching algorithm updated.
 *
 * Revision 1.3  2002/06/10 18:37:33  dmac
 * New class that handles the node distance function added.
 *
 * Revision 1.2  2002/05/16 22:50:11  dmac
 * Several bugs fixed and a new DAG-derived class added
 *
 * Revision 1.1.1.1  2002/03/10 07:38:07  dmac
 * Imported sources
 *
 *
 *-----------------------------------------------------------------------*/

//#include <LEDA/graphwin.h>
//#include <LEDA/graph.h>
//#include <LEDA/node_array.h>
#include <float.h>
#include <stdio.h>
#include <ctype.h>
#include <vector>

// VC++6 -> Include CImg before DAG.h
#include "../CImg/CImg.h"

#include "DAG.h"
#include "ShockGraph.h"
//#include "VisualDAG.h"

#include "ImageProcessing.h"

#include "HelperFunctions.h"
#include "Exceptions.h"
#include "SGDistMeasurer.h"
#include "ApproxPoly.h"

#include "skeleton.h"
#include "BasicUtils.h"

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

#define SA_GROWING_FACTOR 50

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

using namespace std;
using namespace dml;
using namespace leda;

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; }
	
	/*!
	@brief Helper function to compare the edges that will be reverted.
	We want to consider the nodes with greatest radius first.
	
	  The compare function ussuallry performs elem1 - elem2, and elemnts get
	  sorted from smallest to largest. i.e., if elem1 < elem 2, we should return
	  -1 to get an ascending ordering.
	  We want descending sort, so elem1 < elem2 should return 1.
	*/
	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;
	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;
}

//! Computes the object name and view number from the DAG's label
void ShockGraph::ComputeDerivedValues()
{
	DAG::ComputeDerivedValues();
	ComputeObjectViewName();
}

//! Extends the behaviour of the same function in the base class.
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);
	  	}
		
		m_boundaryPts.Read(is);
	}
	
	return is;
}

//! Extends the behaviour of the same function in the base class.
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);
	
	m_boundaryPts.Write(os);
		
	return os;
}

//! Extends the behaviour of the same function in the base class.
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);
}

/*! 
	We MUST have this operator, otherwise C++ would create one for us
	with this exact parameters (and highest precedence over every other
	operator=), and that's no good since m_pSkeleton is a pointer!!!!
*/
ShockGraph& ShockGraph::operator=(const ShockGraph& rhs)
{
	DAG::operator=(rhs);

	area      = rhs.area;
	perimeter = rhs.perimeter;
	
	xmin = rhs.xmin;
	ymin = rhs.ymin;
	xmax = rhs.xmax;
	ymax = rhs.ymax;
	
	// TODO: get a copy of the skel ptr by 
	// m_pSkeleton = new McGillSkeleton(*m_pSkeleton->GetSkeleton());
	// right know, m_pSkeleton is set to NULL. Ie, isn't copied.
	
	m_boundaryPts = rhs.m_boundaryPts;

	return *this;
}

/*! 
	Extends the behaviour of the same function in the base class. This function may
	never be called, but we leave it to make sure that if it is ever called, it does the 
	right thing.
*/
DAG& ShockGraph::operator=(const DAG& rhs)
{
	// Make sure we have a ShockGraph on the rhs
	const ShockGraph* pRhs = dynamic_cast<const ShockGraph*>(&rhs);
	
	if (!pRhs)
		THROW_EXCEPTION("Invalid pointer type.");

	// Now, can can proceed as usual with the right operator=
	return operator=(*pRhs);
}

/*! 
Calls the appropiate node distance function.
*/ 
double ShockGraph::NodeDistance(leda_node u, const DAG& from, leda_node v) const
{
	ASSERT(pDistMeasurer);
	
	return pDistMeasurer->GetDistance(*this, u, from, v);
}

/*! 
Relation function for two Shock Graph nodes.
For now, all nodes of the same type are related. 
*/
bool ShockGraph::AreNodesRelated(leda_node u, const DAG& from, leda_node v) const
{
	return true;
}

//! Creates a new shock graph object.
DAG* ShockGraph::CreateObject() const
{
	return new ShockGraph();
}

//! Creates a new shock graph node object.
DAGNodePtr ShockGraph::CreateNodeObject(NODE_LABEL lbl) const
{
	return new SGNode(lbl, TEST_NODE);
}

/*! 
Creates a new SGNode and fill its contents with the current 
information in the stream.
*/
DAGNodePtr ShockGraph::ReadNode(istream& is) const
{
	DAGNodePtr node(new SGNode);
	
	node->Read(is);
	
	return node;
}

//! Returns the classs name
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)
{
	return true;
}

/*!
	@brief Determines whether the end shock points in the given branch are joints or terminal points.
	
	In general, joint points will be outliers wrt the points in the branch.
*/
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);
}

/*!
@brief Groups the points according to their type. It also computes the
derived values right now, since the points will not be subdevided further.
*/
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;
	
	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, compParams.dMaxYDiff);
	
	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 = poly.m_knots[i].seg.m;
				double m1 = poly.m_knots[i + 1].seg.m;
// 				if(m0 < 0 && m1 < 0)
// 				{
// 					m0 = fabs(m0);
// 					m1 = fabs(m1);
// 				}
				double mr = fabs(m1 - m0) / fabs((fabs(m0) >= fabs(m1))? 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;
}

/*!
@brief Traverses the skeleton and creates a "pre" Shock Graph
*/
/* static */ void ShockGraph::ShowBranchInfo(sg::DDSEdge* e, sg::DDSNode* n, int level)
{
	using namespace sg;
	
	// Then, keep creating and connecting nodes
	DDSEdgeVect& edges = n->getEdges();
	DDSEdgeVect::iterator I;
	DDSEdge* ee;
	
	for (int i = 0; i < level; i++)
		cerr << "\t";
				
	cerr << e->getN1()->fp.p << ", " << e->getN2()->fp.p << endl;
	
	// Check that the flow is fine wrt end points
#ifdef _DEBUG
	sg::FluxPointList& fpl = e->getFluxPoints();
	int last = fpl.size() - 1;
	
	if (last > 2 && fpl[0].p.distanceToPt(fpl[1].p) >=  fpl[0].p.distanceToPt(fpl[last - 1].p))
	{
		DBG_MSG("Wrong flow direction")
		cerr << "dist from " << fpl[0].p << " to " << fpl[1].p << " >= dist from "
			<< fpl[0].p << " to " << fpl[last - 1].p;
	}
	
	ASSERT(last <= 2 || fpl[last].p.distanceToPt(fpl[last - 1].p) <  fpl[last].p.distanceToPt(fpl[1].p));
#endif //_DEBUG

	for(I = edges.begin(); I != edges.end(); I++)
	{
		ee = *I;
		
		if (e != ee)
			ShowBranchInfo(ee, ee->getN1() != n ? ee->getN1():ee->getN2(), level + 1);
	}
}

/*!
@brief Traverses the skeleton and creates a "pre" Shock Graph
*/
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, v, pNewNode);
			ConnectNodes(v, pNewNode, ee, ee->getN1() != n ? ee->getN1():ee->getN2());
		}
	}
}

/*!
@brief Revers the direction of an edge and moves all the edges that meet
at a junction.

  All the children of a node are linked with each other at a juntion, so
  when we move a node up, we need to redistribute the relations in the juntion as well.
  
	we start with parent->x and check for to cases: parent->x or parent->x->children
*/
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));
}

/*!
	@brief "Tries" to relabel type 3 nodes acting as a 2s or 4s. 
	i.e., a type 3 branch that is a sort of neck between two type 1s 
	of opposite directions; or a type 3 that is an implicit root. 
*/
void ShockGraph::TryRelabel3As2or4(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);
	}
}

/*!
  @brief Inserts type 4 and type 2 nodes, which are local maxima or minima
  SINGLE shock points.

  If there are a bunch of type 1 nodes conectd at a point and they have
  the same direction with respect to the join point, then the joint is either a
  type 2 or a type 4.
*/
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);
}

/*!
@brief Computes the discrete divergence skeleton of a closed curve in a ppm file.
*/
sg::DDSkeleton* ShockGraph::ComputeDDSkeleton(const char* szPPMFileName)
{
	using namespace sg;
	
	cout << endl << "Reading " << szPPMFileName << "... " << flush;
	
	cimg_library::CImg<unsigned char> image(szPPMFileName);

	if (image.dimx() == 0 && image.dimy() == 0)
	{
		cerr << "Can't open " << szPPMFileName << endl;
		return NULL;
	}
	
	sg::SimpleShapeMaker ssm(image.dimx(), image.dimy());

	int x, y;
	
	// update the SimpleShape
	for (y=0; y < image.dimy(); y++)
		for (x=0; x < image.dimx(); x++)
			ssm(x,y) = (image(x, y) != 0);
	
	cout << "Computing the flux-based skeleton..." << flush;

	// the shape is being made here
	SimpleShape ss = *((SimpleShape*)ssm.getShape());
	
	ss.getBounds(&xmin,&xmax, &ymin,&ymax);
	
	// 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
	DivergenceSkeletonMaker dsm(dm);
	
	double thresh=-2.0, mag=1.0;
	
	ASSERT(compParams.dSkTreshold >= 0);
	
	return dsm.getDiscreteDivergenceSkeleton(1.0 / mag, -compParams.dSkTreshold);
}

/*!
	right_t1 and right_t2 parameter values between which the contour segment on the right is.
	left_t1 and left_t1 parameter values between which the contour segment on the left is.
*/
void ShockGraph::DetectLigatureNodes(sg::DDSkeleton* sk)
{
	using namespace sg;

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

	std::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);
	}
}

/*!
	@brief Collect the edges that need to be reverted according to time of formation.
	Edges must be directed from last-to-form to earliest-to-form.
*/
void ShockGraph::DirectEdgesAcordingToTime()
{
	bool bContinue;
	leda_edge e;
	int nNoInfLoop = 0;
	
	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);
}

/*!
	@brief Add common root node, generalize type 2 and 4, compute node roles.
	
	common root node: add a dummy root node so that there is a single root.
	generalize 2 and 4: Type 3 that are children of root are type 4 and type 3s 
		that are terminal nodes with two type 1s as parents are type 2.
	node roles: assigns roles according to the sg grammar.
*/
void ShockGraph::FinalizeSGConstruction()
{
	// Nodes that do not have a parent should be children of the root node.
	leda_node root = NewNode(new SGNode("#", ROOT));
	leda_node v, p1;
	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);
		
		/* OPTIONAL STEP */
		if (NodeType(v) == SHOCK_3) // MAY need to be relabeled.
			TryRelabel3As2or4(v);	
			
		nodeType = NodeType(v);
		
		// At this point some basic grammatical rules MUST be 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));
		
		// Relabel node using grammatical knowledge whenever there are errors...
		// Rule 4 errors: A 1 with many 3s as children => all 3's are 1's
		if (nodeType == SHOCK_3)
		{
			p1 = GetFirstParent(v);
			
			if (NodeType(p1) == SHOCK_1 && indeg(v) == 1 && outdeg(p1) > 1)
				UnsafeGetSGNode(v)->SetNodeType(SHOCK_1);
			else if (NodeType(p1) == SHOCK_3 && indeg(v) == 1 && outdeg(p1) > 1)
				UnsafeGetSGNode(v)->SetNodeType(SHOCK_1);
		}

		ComputeNodeRole(v);
	}
}

/*!
@brief Maps the edges in the DDSkeleton to nodes in a shock graph and copies the shocks' info.

  The DDSkeleton has the shock point info stored in the edges of the graph. The endpoint shocks are the nodes
  in this graph and the edges conect these endpoints.
*/
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);

#ifdef _DEBUG
//	cerr << "Side 1" << endl;
//	ShowBranchInfo(se, se->getN1(), 1);
//	cerr << "Side 2" << endl;
//	ShowBranchInfo(se, se->getN2(), 1);
#endif //_DEBUG

	ConnectNodes(u, pNode, se, se->getN1());

	ConnectNodes(u, pFirstNode, se, se->getN2());

	//ShowGraph(*this);
	Insert4sAnd2s();

	DirectEdgesAcordingToTime();

	FinalizeSGConstruction();

	ComputeDerivedValues();
	
	// New functionality (6/Jun/2006)
	if (m_pSkeleton)
		m_pSkeleton->ApproxSkeletonWithLines(compParams.dMinError, 
			compParams.dMinSlope, 10, compParams.dMaxYDiff);

#ifdef _DEBUG	
	// Validate the grammatical consistence of the graph
//	if ((u = CheckGrammar()) != nil)
//	{
//		cerr << "\nThe SG grammar of this graph is invalid."
//			" However, this is not a serious problem."
//			"\n The first grammatical error was found at node: " 
//			<< GetSGNode(u)->GetLblForGraph() << endl;
//	}
#endif //_DEBUG

	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;
}

/*!
    @brief 1D Gaussian smoothing of input into output, where length
	is the number of valid entries.

	Author(s): Carlos Phillips. McGill Univeristy
	Created: 20 Jul 1999
	Modified: 20 Jan 2003 by Diego Macrini
 
*/
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(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;
	}
}

/*!
	@brief The direction of branch in u with respect to the joint
	point with v.
*/
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;
	
	pNodeU->m_shocks.Print(cerr);
	pNodeV->m_shocks.Print(cerr);
	ASSERT(false);
	return vStart;
}

/*!
	@brief Checks that there is a grammatical rule for each node and its children.
	
	It traverses the tree (rooted DAG) in order. It stops atthe first rule violation.
	
	Note: subgraphs rooted at a node with multiple parents or with an ancestor
	with multiple parents whould be checked multiple times. While this is easy to fix,
	it turns out that in the SG grammar, those cases are trivial (subgraphs with a single node),
	and so it cheaper to checke them multiple times (at most two) than to keep track of visits.
	
	@param root is the root node of the shock graph
	@return leda_node nil if all rules are valid. o.w., first node (in depth-first order) with invalid rule.
*/
leda_node ShockGraph::CheckGrammar(leda_node root /* = nil */) const
{
	if (root == nil)
		root = GetFirstRootNode();
		
	if (GetRewriteRule(root) == INV_RULE)
		return root;
	
	leda_node child, invNode;
	
	forall_adj_nodes(child, root)
	{
		invNode = CheckGrammar(child);
		
		if (invNode != nil)
			return invNode;
	}
	
	return nil;
}

/*!
	@brief Checks that all adjacent nodes are of valid types and have exactly 1 in-edge
	
	@param r root node
	@param validType arrat of valid node tipes, ie, SHOCK_1, SHOCK_2, ...
	@param pInvJointPt (optinal) pointer to invalid joint point between r and it children.
	@return true if valid. false otherwise.
*/
bool ShockGraph::CheckChildNodeTypes(leda_node r, SmartArray<int> validTypes, 
	const ShockInfo* pInvJointPt /* = NULL */) const
{
	int i, n = validTypes.GetTailIdx() + 1;
	leda_node c;
	bool righLbl;
	
	forall_adj_nodes(c, r)
	{
		if (indeg(c) > 1)
			return false;
			
		if (pInvJointPt != NULL && (*pInvJointPt) == GetJointPoint(r, c))
			return false;
		
		righLbl = false;
		
		for (i = 0; !righLbl && i < n; i++)
			if (NodeType(c) == validTypes[i])
				righLbl = true;
		
		if (!righLbl)
			return false;
	}
	
	return true;
}

/*!
	@brief Gets the generating rule for a node and its children, ie, the rewrite rule.
	
	Note: subgraphs rooted at a node with multiple parents or with an ancestor
	with multiple parents whould be checked multiple times. While this is easy to fix,
	it turns out that in the SG grammar, those cases are trivial (subgraphs with a single node),
	and so it cheaper to checke them multiple times (at most two) than to keep track of visits.
*/
GRAMMAR_RULE ShockGraph::GetRewriteRule(leda_node r) const
{	
	SmartArray<int> validTypes(2); // don't need more than 2 elements
	
	if (outdeg(r) == 0)
	{
		switch (NodeType(r))
		{
			case SHOCK_4: return RULE_7;
			case SHOCK_3: return RULE_8;
			case SHOCK_2: return RULE_9;
			case SHOCK_1: return RULE_10;
		}
	}
	else
	{
		switch (NodeType(r))
		{
			case ROOT:
				validTypes.AddTail(SHOCK_3);
				validTypes.AddTail(SHOCK_4);
				return CheckChildNodeTypes(r, validTypes) ? RULE_1:INV_RULE;
			case SHOCK_4:
				validTypes.AddTail(SHOCK_1);
				return CheckChildNodeTypes(r, validTypes) ? RULE_2:INV_RULE;
			case SHOCK_3:
				{
					leda_node p1 = GetFirstParent(r);
					int p1Type = NodeType(p1);
					
					if (indeg(r) == 1 && (p1Type == SHOCK_1 || p1Type == ROOT))
					{
						const ShockInfo* pPt = NULL;
						
						if (p1Type == SHOCK_1)
							pPt = &GetJointPoint(r, p1);
						
						validTypes.AddTail(SHOCK_1);
						return CheckChildNodeTypes(r, validTypes, pPt) ? RULE_3:INV_RULE;
					}
				}
				return INV_RULE;
			case SHOCK_1:
				if (outdeg(r) > 1)
				{
					validTypes.AddTail(SHOCK_1);
					if (CheckChildNodeTypes(r, validTypes))
						return RULE_4;
				}
				else /* outdeg(r) == 1 */
				{
					leda_node c = GetFirstChild(r);
					leda_node p2 = GetFirstParent(c); // there is at least 1 parent
					int childInDeg = indeg(c);
					int childType = NodeType(c);
					
					if (childInDeg >= 2 && p2 == r)
						p2 = GetSecondParent(c); // ensures p2 is second parent
					
					if (childType == SHOCK_1 && childInDeg == 1) 
						return RULE_4; // NEW RULE 4, allow chain of type 1s
					else if (childType == SHOCK_2 && childInDeg == 2 && NodeType(p2) == SHOCK_1)
						return RULE_5;
					else if (childType == SHOCK_3 && childInDeg == 2 && NodeType(p2) == SHOCK_1 &&
						 GetJointPoint(r, c) != GetJointPoint(p2, c))
						return RULE_6;
					else if (childType == SHOCK_3 && childInDeg == 1)
						return RULE_6;
					
				}
				return INV_RULE;
		}
	}	
	return INV_RULE;
}

bool ShockGraph::ComputeFromPPMFile3(const char* szPPMFileName, const SGCP& sgparams)
{
	cimg_library::CImg<float> img(szPPMFileName);	
	
	// If color image, we only keep the "red" chanel
	if (img.dimv() > 1)
		img = img.channel(0);

	// Removes "white holes" and "black satellites" from the image
	RemoveNonMaxComponents(img);

	IMAGEINFO imgInfo;

	imgInfo.pField = CopyImageToField(img);
	imgInfo.strFileName = szPPMFileName;
	
	return ComputeFromPPMFile3(imgInfo, sgparams);
}

/*!
	@brief Gets a set of "TreeSkel" skeletons given by pointers to the skeletons in a vector
	
	The function determine_skeletons() is called which will return
	a vector with pointers to skeletons of the type "TreeSkel"
	In order for this function to work it asks a couple of parameters:
	1 char* input file name
	2 float& start threshold boundary simplification
	3 float& start threshold structural simplification
	4 int& (0 or 1) indicating not to do a 1D boundary search, or to actually do it
	5 int& (0 or 1) indicating not to do a 1D structural search, or to actually do it
	6 float& upperbound for threshold boundary simplification
	7 float& upperbound for threshold structural simplification
	8 float& stepsize for the threshold for boundary simplification
	9 float& stepsize for the threshold for structural simplification
	10 bool give a unsimpified skeleton or not
	Params in sgparams: {dSkTreshold}
	sk_upper, sk_step, sh_upper, sh_step, compDerivSkel
*/
bool ShockGraph::ComputeFromPPMFile3(const IMAGEINFO& imgInfo, const SGCP& sgparams)
{
	Clear(); //make sure the current shock graph is empty

	compParams = sgparams;
	SetDAGLbl(imgInfo.strFileName);
//!!! boundary points from ImgPreProcess
//	m_boundaryPts = imgInfo.boundaryPts;

// 	for (int k = 0; k < m_boundaryPts.GetSize(); k++)
// 	{
// 		if (xmin == 0 || xmin >  m_boundaryPts[k].x)
// 			xmin = m_boundaryPts[k].x;
// 			
// 		if (xmax <  m_boundaryPts[k].x)	
// 			xmax = m_boundaryPts[k].x;
// 			
// 		if (xmin == 0 || ymin >  m_boundaryPts[k].y)
// 			ymin = m_boundaryPts[k].y;
// 			
// 		if (ymax <  m_boundaryPts[k].y)	
// 			ymax = m_boundaryPts[k].y;
// 	}
	//xmin = ymin = 0;
	//xmax = ymax = 300;

	bool compDerivSkel = true;

	SkelCompParams parameters;
	
	//some hard coded parameters... these are suitable for a certain file
	//called anim7.pgm
	
	parameters.pInputfile = imgInfo.strFileName;
	parameters.fThresSimp1 = sgparams.tsTau1;
	parameters.fThresSimp2 = sgparams.tsTau2;
	parameters.bOnlyBound = false;
	parameters.bOnlyStruc = false;
	parameters.bDo2DField = false;
	parameters.fSimp1UpperBound = 71;
	parameters.fSimp1StepSize = 5;
	parameters.fSimp2UpperBound = 101;
	parameters.fSimp2StepSize = 5;
	parameters.bGiveSimpleSkel = false;
	parameters.dMinSlope = sgparams.dMinSlope;
	parameters.dMinError = sgparams.dMinError;
	parameters.dMaxAccelChg = sgparams.dMaxAccelChg;
	parameters.nRecErrorWeightBnd = sgparams.nBndryRecErrorWeight;
	parameters.nRecErrorWeightStr = sgparams.nStrucRecErrorWeight;
	parameters.pInputField = imgInfo.pField;

	/*parameters.pInputField = NULL;

	if (!imgInfo.ptrImg.IsNull())
	{
		const cimg_library::CImg<float>& img = *imgInfo.ptrImg;
		parameters.pInputField = new FIELD<float>(img.dimx(), img.dimy());
		float* const p = parameters.pInputField->data();
		
		for (int y = 0; y < img.dimy(); y++)
			for (int x = 0; x < img.dimx(); x++)
				*(p + y * img.dimx() + x) = img(x, y);
	}*/
	
	if(true)
	{

		AFMMSkeleton* pDDS = AFMMSkeleton::MakeSkeleton(parameters);
		
		if(sgparams.nSimplifyComplete == 1)
			pDDS->Simplify();
		else if(sgparams.nSimplifyExternal == 1)
			pDDS->SimplifyExternal();
		else if(sgparams.nSimplifyInternal == 1)
			pDDS->SimplifyInternal();

		pDDS->GetDimensions(&xmin,&xmax,&ymin,&ymax);	

		if(!pDDS)
			return false;
		
		m_pSkeleton = new McGillSkeleton(pDDS);

		bool thingsOK = ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());

		return thingsOK;
	}
	
	
	
	
	std::vector<TreeSkel*> set_of_skeletons; //vector to store the pointer to the TreeSkel skeletons
	Skeleton get_skeletons;

	set_of_skeletons = get_skeletons.determine_skeletons(parameters);

	TreeSkel* get_one_skel;

	for(unsigned int j = 0; j < set_of_skeletons.size(); j++)  // for each skeleton, test whether it exists. do this by counting the number of branches
	{
		get_one_skel = set_of_skeletons[j];
		int number_branch = 0;
		for(int i = 0 ; i < get_one_skel->branches.Count(); i++) //count the number of branches
		{
			if(get_one_skel->branches[i]->geomLength()) { number_branch++; }
		}
		cerr << "The number of branches in skeleton number " << j << " is " << number_branch << endl;
	}

	if (set_of_skeletons.size() == 0)
		return false;

	m_pSkeleton = TranslateTreeSkelIntoDDS(set_of_skeletons[0], sgparams); // This returns a McGillSkeleton
	delete set_of_skeletons[0];

	bool thingsOK = ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
	
	/*if (thingsOK)
	{
		// Add the SGs from alternative tree skeletons...
		m_derivSGs.ReSize(set_of_skeletons.size() - 1);
			
		for (unsigned int i = 1; thingsOK && i < set_of_skeletons.size(); i++)
		{
			ShockGraph* pSG = new ShockGraph;
			pSG->m_pSkeleton = TranslateTreeSkelIntoDDS(set_of_skeletons[i], sgparams); // This returns a McGillSkeleton
			delete set_of_skeletons[i];
			
			thingsOK = pSG->ComputeSGFromDDSkeleton(((McGillSkeleton*)m_pSkeleton)->GetSkeleton());
			
			if (thingsOK)
				m_derivSGs[i] = pSG;
		}
	}*/

	return thingsOK;
}


/*!
	@brief Text
	
	More text
	Params in sgparams: {dSkTreshold}
*/
/*static*/ McGillSkeleton* ShockGraph::TranslateTreeSkelIntoDDS(TreeSkel* pTreeSkel, const SGCP& sgparams)
{
	sg::DiscreteSegCurve c;
	sg::SimpleShape s(&c,0,300,0,300);
	sg::DistanceTransform dt(&s);
	sg::DivergenceMap dm(dt);
	sg::DDSkeleton* pDDS = new sg::DDSkeleton(dm);
	leda_d_array<long, sg::DDSNode*> map;
	sg::DDSNode *pEndPt1, *pEndPt2;
	sg::DDSEdge* pEdge;
	int j;

//################################ ADDED JANUARI 13 ##############################

// for(int i = 0; i < pTreeSkel->branches.Count(); i++)
// {
// 	cerr << *pTreeSkel->branches[i] << endl;
// }

//################################################################################
	
	
	
	
	//We take the first branch out of the list of branches
	//and add this one to the DDSkeleton
	int nFirstNonZeroLength = 0;
	TreeSkel::TreeBranch& br = *pTreeSkel->branches[nFirstNonZeroLength];
	if(pTreeSkel->branches.Count() > 1)
	{
		while( ((br.end1->branches.Count() != 1) && (br.end2->branches.Count() != 1)) ||
			((br.end1->branches.Count() == 1) && (br.end2->branches.Count() == 1)))
		{	
 			cerr << "Going to next branch in search for an endbranch\n";
			nFirstNonZeroLength++;
			br = *pTreeSkel->branches[nFirstNonZeroLength];
		}
	}
	//We now found the first branch that has a non zero length
	//and will add this branch to pDDS
	
// 	cerr << "Number of branches at endpoint 1: " << br.end1->branches.Count() << " and at end2 : " << br.end2->branches.Count() << endl;
	
	//Add endpoints to DDSNodes
	if(br.end1->branches.Count() == 1)
	{
		pEndPt1 = new sg::DDSNode(TranslateTreePtIntoFluxPt(br.end1));
		//cerr << "Endpoint 1: (" << br.end1->x << "," << br.end1->y << ") number of neighbor branches at this point " << br.end1->branches.Count() << endl;
		pDDS->addNode(pEndPt1);
		pEndPt2 = new sg::DDSNode(TranslateTreePtIntoFluxPt(br.end2));
		//cerr << "Endpoint 2: (" << br.end2->x << "," << br.end2->y << ") number of neighbor branches at this point " << br.end2->branches.Count() << endl;
		pDDS->addNode(pEndPt2);
	}
	else
	{
		pEndPt1 = new sg::DDSNode(TranslateTreePtIntoFluxPt(br.end2));
		//cerr << "Endpoint 1: (" << br.end1->x << "," << br.end1->y << ") number of neighbor branches at this point " << br.end1->branches.Count() << endl;
		pDDS->addNode(pEndPt1);
		pEndPt2 = new sg::DDSNode(TranslateTreePtIntoFluxPt(br.end1));
		//cerr << "Endpoint 2: (" << br.end2->x << "," << br.end2->y << ") number of neighbor branches at this point " << br.end2->branches.Count() << endl;
		pDDS->addNode(pEndPt2);
	}	
	
	
	sg::FluxPointList fluxPts;
	TreeSkel::TreePoint* pTrSkelEnd1 = NULL;
	TreeSkel::TreePoint* pTrSkelEnd2 = NULL;
	//Add the edge to the sDDS, in "normal" order
	//so from end1 -> end2
	if(br.end1->branches.Count() == 1)
	{
		//sg::FluxPointList fluxPts;
		cout << endl << endl << "%new branch" << endl;
		fluxPts.push_back(pEndPt1->fp);
		for (j = 0; j < br.points.Count(); j++)
			fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
		fluxPts.push_back(pEndPt2->fp);
		pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
		pDDS->addEdge(pEdge);
		
		/*TreeSkel::TreePoint* */ pTrSkelEnd1 = br.end1;
		/*TreeSkel::TreePoint* */ pTrSkelEnd2 = br.end2;
	}
	else
	{
		//sg::FluxPointList fluxPts;
		cout << endl << endl << "%new branch" << endl;
		fluxPts.push_back(pEndPt1->fp);
		for (j = br.points.Count() -1; j >= 0; j--)
			fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
		fluxPts.push_back(pEndPt2->fp);
		pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
		pDDS->addEdge(pEdge);
		
		/*TreeSkel::TreePoint* */ pTrSkelEnd1 = br.end2;
		/*TreeSkel::TreePoint* */ pTrSkelEnd2 = br.end1;
	}	
	
	//Recursively add a new edge and a node to the existing graph
	//given: pDDS, current branch, current endpoint
	//Node that belongs to the endpoint
	cerr << "Start adding edges and nodes to the graph\n";
	bool bFlowAway = false;
	ConnectTreeSkelNodes(pDDS, br, pTrSkelEnd1, pEndPt1, bFlowAway);
	bFlowAway = true;
	ConnectTreeSkelNodes(pDDS, br, pTrSkelEnd2, pEndPt2, bFlowAway);
	cerr << "Done with all the adding of branches\n";
	
	return (new McGillSkeleton(pDDS));
}

/*static*/ void ShockGraph::ConnectTreeSkelNodes(sg::DDSkeleton* pDDS, const TreeSkel::TreeBranch& existingBr, 
	const TreeSkel::TreePoint* pExistingEndPt, sg::DDSNode* pExistingNode, bool flowAway)
{
	using namespace sg;
	
	DDSNode* pNewEndPt;
	DDSNode *pEndPt1, *pEndPt2;
	DDSEdge* pEdge;
	int i,j;
	//for all branches that are connected to the existing endpoint
	for (i = 0; i < pExistingEndPt->branches.Count(); i++)
	{
		//cerr << "\nLooking at 1 branch..." << endl;
		//for all branches different from the one that already
		//existed for this endpoint
		const TreeSkel::TreeBranch& br = *pExistingEndPt->branches[i];

// 		if (br.geomLength() == 0)
// 			ASSERT(true)

// 		cerr << "\nNumber of innerpoints in the branch is: " << br.points.Count() << endl;

		//ASSERT(br.geomLength() != 0)

		if (br != existingBr)
		{
			//cerr << "\nWe actually want to add it to the structure..." << endl;
			sg::FluxPointList fluxPts;//(br.points.Count() + 2); // plus endpoints
			
			if (br.end1 != pExistingEndPt)
			//Endpoint 1 is not yet a node
			{
				pNewEndPt = new sg::DDSNode(TranslateTreePtIntoFluxPt(br.end1));
				pDDS->addNode(pNewEndPt);
				//Add edge in from end2 -> end1, thus away from the given one
				if(flowAway == true) //start from pExistingNode which is end2
				{
					pEndPt1 = pExistingNode;
					pEndPt2 = pNewEndPt;
					//cout << endl << endl << "%new branch" << endl;
					fluxPts.push_back(pEndPt1->fp);
					for (j = br.points.Count() -1; j > -1; j--)
						fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
					fluxPts.push_back(pEndPt2->fp);
					pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
					pDDS->addEdge(pEdge);
					ConnectTreeSkelNodes(pDDS, br, br.end1, pEndPt2, flowAway);
				}
				else //flowAway == false, so start from pNewEndPt which is end1
				{
					pEndPt1 = pNewEndPt;
					pEndPt2 = pExistingNode;
					//cout << endl << endl << "%new branch" << endl;
					fluxPts.push_back(pEndPt1->fp);
					for (j = 0; j < br.points.Count(); j++)
						fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
					fluxPts.push_back(pEndPt2->fp);
					pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
					pDDS->addEdge(pEdge);
					ConnectTreeSkelNodes(pDDS, br, br.end1, pEndPt2, flowAway);
				}
			}
			else if (br.end2 != pExistingEndPt)
			//Endpoint 2 is not yet a node
			{
				pNewEndPt = new DDSNode(TranslateTreePtIntoFluxPt(br.end2));
				pDDS->addNode(pNewEndPt);
				
				if(flowAway == true) //start from pExistingNode which is end1
				{
					pEndPt1 = pExistingNode;
					pEndPt2 = pNewEndPt;
					//cout << endl << endl << "%new branch" << endl;
					fluxPts.push_back(pEndPt1->fp);
					for (j = 0; j < br.points.Count(); j++)
						fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
					fluxPts.push_back(pEndPt2->fp);
					pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
					pDDS->addEdge(pEdge);
					ConnectTreeSkelNodes(pDDS, br, br.end2, pEndPt2, flowAway);
				}
				else //flowAway == false, so start from pNewEndPt which is end2
				{	
					pEndPt1 = pNewEndPt;
					pEndPt2 = pExistingNode;
					//cout << endl << endl << "%new branch" << endl;
					fluxPts.push_back(pEndPt1->fp);
					for (j = br.points.Count() -1; j > -1; j--)
						fluxPts.push_back(TranslateTreePtIntoFluxPt(br.points[j]));
					fluxPts.push_back(pEndPt2->fp);
					pEdge = new sg::DDSEdge(pDDS, fluxPts, pEndPt1, pEndPt2);
					pDDS->addEdge(pEdge);
					ConnectTreeSkelNodes(pDDS, br, br.end2, pEndPt2, flowAway);
				}
			}
			else
				cerr << "We don't deal with circular branches right now\n";
		}
	}
}

/*!
	@brief Text
	
	More text
*/
/*static*/ sg::FluxPoint ShockGraph::TranslateTreePtIntoFluxPt(const TreeSkel::TreePoint* pTreePt)
{
	sg::FluxPoint fp;
	
	fp.p.x = pTreePt->x;
	fp.p.y = pTreePt->y;
	fp.dist = pTreePt->dRadius; //pTreePt.radius;
	//cout << pTreePt->dRadius << endl;
	fp.col = 0; // not used
	fp.val = 0; // not used
	
	return fp;
}


