Quad tree search pm1quadtree.cc

From Minor Miracle Software
Jump to: navigation, search
/*
 * Jesse B. Dooley
 * Date: March 8, 2000
 * For:  CMSC 420-0201 Spring 2000
 *       Prof. Michelle Hugue
 * Part2: Minimum Spanning Tree and PM Quadtree
 * Due: 18-Apr-2000
 * Module: pm1quadtree.cc
 */

/*
 * isolated points can legally be included as a self
 * terminating route with zero length.
 *
 */

#ifndef PM1QUADTREE_CPP
#define PM1QUADTREE_CPP

#include <algorithm>
#include <iostream.h>
#include <float.h>

#include "TemasBase.h"
#include "List.cc"
#include "TwoDPointFloat.cc" // for centers
#include "site.cc" // yes, it holds just sites
#include "route.cc"

using namespace std;

class pm1quadtree {

	enum Direction{ NORTH, EAST, WEST, SOUTH, NORTHEAST, NORTHWEST, SOUTHEAST, SOUTHWEST, CENTER, ROOT, NO_DIRECTION };

	enum NodeType{ GRAY, BLACK, WHITE, NO_TYPE };

	struct quadrant {
		// center for the quadrant
		TwoDPointFloat Center_TwoDPointFloat;

		// corners
		TwoDPointFloat NW_TwoDPointFloat;
		TwoDPointFloat SE_TwoDPointFloat;
		TwoDPointFloat NE_TwoDPointFloat;
		TwoDPointFloat SW_TwoDPointFloat;

		// length for one quadrant side
		float width_float;

		// direction to parent quadrant
		Direction DIRECTION;

		NodeType TYPE; // BLACK data point; acording to Samet
					   // WHITE empty leaf
		               // GRAY internal

		quadrant* Parent_p;
		// data pointers
		// List<site*> route_List;
		// Sites are routes with zero weight
	    List<route*> route_List;

		// quadrant pointers
		quadrant* NORTHEAST_quadrant_p;
		quadrant* NORTHWEST_quadrant_p;
		quadrant* SOUTHWEST_quadrant_p;
		quadrant* SOUTHEAST_quadrant_p;
	};

	// Circle storage
	site* Circle_Site_p;
	double Radius_Circle_double;

	// Nearest Route storage
	route* NearestRoute_route_p;
	bool Intercept_bool; // true if two edges intercept
	quadrant *root_quadrant_p; // once we start quadranting

////////////METHODS////////////////////////

	void SetCorners( quadrant* p )
	{
		p->NW_TwoDPointFloat.xaxis( p->Center_TwoDPointFloat.xaxis() - p->width_float );
		p->NW_TwoDPointFloat.yaxis( p->Center_TwoDPointFloat.yaxis() + p->width_float );

		p->SE_TwoDPointFloat.xaxis( p->Center_TwoDPointFloat.xaxis() + p->width_float );
		p->SE_TwoDPointFloat.yaxis( p->Center_TwoDPointFloat.yaxis() - p->width_float );

		p->NE_TwoDPointFloat.xaxis( p->Center_TwoDPointFloat.xaxis() + p->width_float );
		p->NE_TwoDPointFloat.yaxis( p->Center_TwoDPointFloat.yaxis() + p->width_float );

		p->SW_TwoDPointFloat.xaxis( p->Center_TwoDPointFloat.xaxis() - p->width_float );
		p->SW_TwoDPointFloat.yaxis( p->Center_TwoDPointFloat.yaxis() - p->width_float );
	} // SetCorners( quadrant* p )


	void PRINT_DIRECTION( Direction D )
	{
		switch( D )
		{
			case NORTHEAST: cout << "NE";//"NORTHEAST";
				break;
			case NORTHWEST: cout << "NW";//"NORTHWEST";
				break;
			case SOUTHWEST: cout << "SW";//"SOUTHWEST";
				break;
			case SOUTHEAST: cout << "SE";//"SOUTHEAST";
				break;
//			case NORTH:     cout << "NORTH";//"NORTH";
				break;
//			case SOUTH:     cout << "SOUTH";//"SOUTH";
				break;
//			case EAST:      cout << "EAST";//"EAST";
				break;
//			case WEST:      cout << "WEST";//"WEST";
				break;
//			case CENTER:    cout << "CENTER";//"CENTER";
				break;
//			case ROOT:      cout << "ROOT";//"ROOT";
				break;
//			case NO_DIRECTION: cout << "NO_DIRECTION";
			default:
			;
		};
	} // PRINT_DIRECTION

	void PrintNodeType( NodeType N )
	{
		switch( N )
		{
			case GRAY: cout << "GRAY";
				break;
			case BLACK: cout << "BLACK";
				break;
			case WHITE: cout << "WHITE";
				break;
			case NO_TYPE: cout << "NO_TYPE";
				break;
			default:
				cout << "pm1quadtree:PrintNodeType: No NodeType found." << endl;
		};
	}

	void copynode( quadrant *ptrn, quadrant *ptr )
	{
	   if(ptr)
	   {
		 ptrn = new quadrant;

	     ptrn->Center_TwoDPointFloat = ptr->Center_TwoDPointFloat;
	     ptrn->width_float = ptr->width_float;

		ptrn->route_List = ptr->route_List;
		ptrn->TYPE = ptr->TYPE;
		ptrn->DIRECTION = ptr->DIRECTION;

	     copynode( ptrn->NORTHEAST_quadrant_p, ptr->NORTHEAST_quadrant_p );
	     copynode( ptrn->NORTHWEST_quadrant_p, ptr->NORTHWEST_quadrant_p );
	     copynode( ptrn->SOUTHWEST_quadrant_p, ptr->SOUTHWEST_quadrant_p );
	     copynode( ptrn->SOUTHEAST_quadrant_p, ptr->SOUTHEAST_quadrant_p );
	   }
	} // copynode

	void deletenode( quadrant* ptn )
	{
		if( ptn )
		{
			ptn->Center_TwoDPointFloat.~TwoDPointFloat();
			ptn->NW_TwoDPointFloat.~TwoDPointFloat();
			ptn->SE_TwoDPointFloat.~TwoDPointFloat();
			ptn->NE_TwoDPointFloat.~TwoDPointFloat();
			ptn->SW_TwoDPointFloat.~TwoDPointFloat();

			ptn->Parent_p = NULL;
			ptn->route_List.erase();
			deletenode( ptn->NORTHEAST_quadrant_p );
			deletenode( ptn->NORTHWEST_quadrant_p );
			deletenode( ptn->SOUTHEAST_quadrant_p );
			deletenode( ptn->SOUTHWEST_quadrant_p );

		    delete ptn;
		}
	} // deletenode( quadrant* ptn )

// quadrant methods
	void printQuadrant( quadrant *Q_p )
	{
		if( Q_p == NULL )
			return;

		printSquare( Q_p );
		printQuadrant( Q_p->NORTHWEST_quadrant_p );
		printQuadrant( Q_p->NORTHEAST_quadrant_p );
		printQuadrant( Q_p->SOUTHWEST_quadrant_p );
		printQuadrant( Q_p->SOUTHEAST_quadrant_p );
	}

	void printSquare( quadrant *Q_p )
	{
		if( Q_p == NULL )
			return;
//Q_p->Center_TwoDPointFloat.PRINT();

		PRINT_DIRECTION( Q_p->DIRECTION );

		List<site*> R_list;

		Q_p->route_List.begin();
	    for( int itr = 0; itr < Q_p->route_List.size(); itr++ )
	    {
			route *R = *(Q_p->route_List).getcurrent();
// if vertex, print V, else, if Qedge print Q, print Name
			if( ISOLATED_ELEMENT( R ) )
			{
				cout << "VQ" << R->StartName();
			}
			else if( PT_IN_QUADRANT( R->StartSite(), Q_p ) )
			{
				if( R_list.find( R->StartSite() ) == false )
				{
					cout << "VQ";
					R_list.AddToList( R->StartSite() );
					cout << R->StartName();
				}
			}
			else if( PT_IN_QUADRANT( R->FinishSite(), Q_p ) )
			{
				if( R_list.find( R->FinishSite() ) == false )
				{
					cout << "VQ";
					R_list.AddToList( R->FinishSite() );
					cout << R->FinishName();
				}
			}
			else if( CLIP_SQUARE( R, Q_p ) )
			{
				cout << "Q";
			}

			Q_p->route_List.next();
	    } // for

		if( Q_p->TYPE == WHITE )
		{
			cout << 'N';
		}

		if( LEAF(Q_p) )
		{
			cout << endl;
		}
	} // printSquare

	Direction calcDirection( TwoDPointFloat S_p,   // site point
							 TwoDPointFloat Q_p ) // quadrant point
	{
	// calculates which direction I'm going.

        // Center
		if( S_p == Q_p ) // quadrant center
		    return CENTER;

		// north
		if( S_p.xaxis() == Q_p.xaxis() && // X does not change
			S_p.yaxis() > Q_p.yaxis() )//  &&    // Y does change S_p->yaxis() < Q_p->yaxis() + Q_p->width_float ) // so we don't leave the quadrant
		    return NORTH;

		// south
		if( S_p.xaxis() == Q_p.xaxis() && // X does not change
			S_p.yaxis() < Q_p.yaxis() )//  &&    // Y does change S_p->yaxis() > Q_p->Center_TwoDPointFloat.yaxis() - Q_p->width_float ) // so we don't leave the quadrant
		    return SOUTH;

		// east
		if( S_p.xaxis() > Q_p.xaxis()  && // X does change
			S_p.yaxis() == Q_p.yaxis() ) // &&    // Y does not change	S_p->xaxis() < Q_p->Center_TwoDPointFloat.xaxis()  + Q_p->width_float ) // so we don't leave the quadrant
		    return EAST;

		// west
		if( S_p.xaxis() < Q_p.xaxis() && // X does change
			S_p.yaxis() == Q_p.yaxis() ) // &&    // Y does not change	S_p->xaxis() > Q_p->Center_TwoDPointFloat.xaxis()  - Q_p->width_float ) // so we don't leave the quadrant
		    return WEST;

 		// NORTHEAST is true if Center X < Point X && Center Y < Point Y
		if( S_p.xaxis() > Q_p.xaxis() &&
			S_p.yaxis() > Q_p.yaxis() )
			return NORTHEAST;

		// NORTHWEST is true if Center X > Point x && Center Y < Point Y
		if( S_p.xaxis() < Q_p.xaxis() &&
		    S_p.yaxis() > Q_p.yaxis() )
			return NORTHWEST;

		// SW is true if Center X > Point X && Center Y > Point Y
		if( S_p.xaxis() < Q_p.xaxis() &&
			S_p.yaxis() < Q_p.yaxis() )
			return SOUTHWEST;

		// SOUTHEAST is true if Center X < TwoDPointFloat X && Center Y > TwoDPointFloat Y
		if( S_p.xaxis() > Q_p.xaxis() &&
			S_p.yaxis() < Q_p.yaxis() )
			return SOUTHEAST;

		return NO_DIRECTION;
	} // calcDirection

	// quadrant methods
	void initquadrant( quadrant *P )
	{
		P->width_float = FLT_MAX;

		P->TYPE = NO_TYPE;

		P->DIRECTION = NO_DIRECTION;

		P->Parent_p = NULL;
	    P->NORTHEAST_quadrant_p = NULL;
	    P->NORTHWEST_quadrant_p = NULL;
	    P->SOUTHWEST_quadrant_p = NULL;
	    P->SOUTHEAST_quadrant_p = NULL;
	}

	bool RECTANGLE_CHECK( site* S_p, quadrant* Rectangle_p )
	{
		if( S_p == NULL || Rectangle_p == NULL )
			return false;

		return ( Rectangle_p->NW_TwoDPointFloat.xaxis() <= S_p->xaxis() &&
				 Rectangle_p->NW_TwoDPointFloat.yaxis() >= S_p->yaxis() &&
				 Rectangle_p->SE_TwoDPointFloat.xaxis() >= S_p->xaxis() &&
				 Rectangle_p->SE_TwoDPointFloat.yaxis() <= S_p->yaxis() );
	} // RECTANGLE_CHECK

// Rectangle Search Methods
	bool  RECTANGLE_CLIP_QUADRANT( quadrant* Q_p, quadrant* Rectangle_p )
	{
	// Case 0: quadrant dosen't exist
		if( Q_p == NULL || Rectangle_p == NULL )
			return false;

	// Case 1: quadrant is inside the rectangle
		if( Rectangle_p->NW_TwoDPointFloat.xaxis() <= Q_p->NW_TwoDPointFloat.xaxis() &&
			Rectangle_p->NW_TwoDPointFloat.yaxis() >= Q_p->NW_TwoDPointFloat.yaxis() &&
			Rectangle_p->NE_TwoDPointFloat.xaxis() >= Q_p->NE_TwoDPointFloat.xaxis() &&
		 	Rectangle_p->NE_TwoDPointFloat.yaxis() >= Q_p->NE_TwoDPointFloat.yaxis() &&
			Rectangle_p->SW_TwoDPointFloat.xaxis() <= Q_p->SW_TwoDPointFloat.xaxis() &&
			Rectangle_p->SW_TwoDPointFloat.yaxis() <= Q_p->SW_TwoDPointFloat.yaxis() &&
			Rectangle_p->SE_TwoDPointFloat.xaxis() >= Q_p->SE_TwoDPointFloat.xaxis() &&
			Rectangle_p->SE_TwoDPointFloat.yaxis() <= Q_p->SE_TwoDPointFloat.yaxis() )
			return true;

	// Case 2: rectangle is inside the quadrant
		if( Rectangle_p->NW_TwoDPointFloat.xaxis() >= Q_p->NW_TwoDPointFloat.xaxis() &&
			Rectangle_p->NW_TwoDPointFloat.yaxis() <= Q_p->NW_TwoDPointFloat.yaxis() &&
			Rectangle_p->NE_TwoDPointFloat.xaxis() <= Q_p->NE_TwoDPointFloat.xaxis() &&
		 	Rectangle_p->NE_TwoDPointFloat.yaxis() <= Q_p->NE_TwoDPointFloat.yaxis() &&
			Rectangle_p->SW_TwoDPointFloat.xaxis() >= Q_p->SW_TwoDPointFloat.xaxis() &&
			Rectangle_p->SW_TwoDPointFloat.yaxis() >= Q_p->SW_TwoDPointFloat.yaxis() &&
			Rectangle_p->SE_TwoDPointFloat.xaxis() <= Q_p->SE_TwoDPointFloat.xaxis() &&
			Rectangle_p->SE_TwoDPointFloat.yaxis() >= Q_p->SE_TwoDPointFloat.yaxis() )
			return true;
	// Case 3: quadrant and rectangle maybe co-terminus
	// This case covered above.

	// Case 4: quadrant and rectangle clip each other
	// Note: If one line clips then another line clips so TRUE can be returned
	// after one success. This also covers the one-point corner clip.

	// Case 4.1: Rectangle's North side clips the quadrant's East or West sides.
			if( intersect( Rectangle_p->NW_TwoDPointFloat, // North side
						   Rectangle_p->NE_TwoDPointFloat,
						   Q_p->NE_TwoDPointFloat,               // East side
						   Q_p->SE_TwoDPointFloat ) ||
				intersect( Rectangle_p->NW_TwoDPointFloat,
						   Rectangle_p->NE_TwoDPointFloat,
						   Q_p->NW_TwoDPointFloat,               // West
						   Q_p->SW_TwoDPointFloat ) )
				return true;

	// Case 4.2: Rectangle's East side clips the quadrant's North or South sides.
			if( intersect( Rectangle_p->NE_TwoDPointFloat, // East side
						   Rectangle_p->SE_TwoDPointFloat,
						   Q_p->NE_TwoDPointFloat,               // North side
						   Q_p->NW_TwoDPointFloat ) ||
				intersect( Rectangle_p->NE_TwoDPointFloat,
						   Rectangle_p->SE_TwoDPointFloat,
						   Q_p->SW_TwoDPointFloat,               // South
						   Q_p->SE_TwoDPointFloat ) )
				return true;

	// Case 4.3: Rectangle's South side clips the quadrant's East or West sides.
			if( intersect( Rectangle_p->SE_TwoDPointFloat, // South side
						   Rectangle_p->SW_TwoDPointFloat,
						   Q_p->NE_TwoDPointFloat,               // East side
						   Q_p->SE_TwoDPointFloat ) ||
				intersect( Rectangle_p->SE_TwoDPointFloat,
						   Rectangle_p->SW_TwoDPointFloat,
						   Q_p->SW_TwoDPointFloat,               // West
						   Q_p->NW_TwoDPointFloat ) )
				return true;

	// Case 4.4: Rectangle's West side clips the quadrant's North or South sides.
			if( intersect( Rectangle_p->NW_TwoDPointFloat, // West side
						   Rectangle_p->SW_TwoDPointFloat,
						   Q_p->NE_TwoDPointFloat,               // North side
						   Q_p->NW_TwoDPointFloat ) ||
				intersect( Rectangle_p->NW_TwoDPointFloat,
						   Rectangle_p->SW_TwoDPointFloat,
						   Q_p->SW_TwoDPointFloat,               // South
						   Q_p->SE_TwoDPointFloat ) )
				return true;

		return false;
	} // QUADRANT_IN_RECTANGLE( quadrant* P )

	List<string> RECTANGLE_SEARCH( quadrant* Q_p, quadrant* Rectangle_p )
	{
		List<string> L;
		if( Q_p == NULL || Rectangle_p == NULL)
			return L;

		printSquare( Q_p );
		if( RECTANGLE_CLIP_QUADRANT( Q_p, Rectangle_p ) == false )
			return L;

		Q_p->route_List.begin();
		for( int i = 0; i < Q_p->route_List.size(); i++ )
		{// test if any point in list is inside the rectangle
			if( RECTANGLE_CHECK( (*Q_p->route_List.getcurrent())->StartSite(), Rectangle_p ) )
				L += (*Q_p->route_List.getcurrent())->StartSite()->name();

			if( RECTANGLE_CHECK( (*Q_p->route_List.getcurrent())->FinishSite(), Rectangle_p ) )
				L += (*Q_p->route_List.getcurrent())->FinishSite()->name();

			Q_p->route_List.next();
		}

		L += RECTANGLE_SEARCH( Q_p->NORTHWEST_quadrant_p, Rectangle_p );
		L += RECTANGLE_SEARCH( Q_p->NORTHEAST_quadrant_p, Rectangle_p );
		L += RECTANGLE_SEARCH( Q_p->SOUTHWEST_quadrant_p, Rectangle_p );
		L += RECTANGLE_SEARCH( Q_p->SOUTHEAST_quadrant_p, Rectangle_p );

		return L;
	} // RECTANGLE_SEARCH

// end Rectangle Search Methods

/// PM methods
	bool BASE_SITE( site* S_p )
	{
		if( S_p == NULL )
			return false;

		return ( S_p->type() == 'B' || S_p->type() == 'S' || S_p->type() == 'L' );
	}

	void PM_INSERT( List<route*> P, quadrant* Q_p )
	{
		if( empty() )
		{
			repartitionRoot();
			Q_p = root_quadrant_p;
		}

		if( Q_p == NULL || P.empty() )
			return;

		if( Q_p->width_float <= MinRouteDistance )
			return;

		List<route*> Q_List = QUADRANT_CHECK( P, Q_p );

		if( Q_List.empty() == false ) // valid quadrant edges, continue
		{
			if( Intercept_bool )
				return;

			if( LEAF( Q_p ) ) // we have a leaf
			{
				Q_List += Q_p->route_List; // add the new edges to it.
				if( PM1_CHECK( Q_List, Q_p ) ) // change this to PM2, PM3 etc.
				{ // all edges valid for this
					Q_p->route_List = Q_List;
					Q_p->TYPE = BLACK; // it is now
					return;
				}
				else // something doesn't belong
				{
					SPLIT_PM_NODE( Q_p );
					Q_p->route_List.erase(); // clean up
				}
			}
			// we are going here anyway
			PM_INSERT( Q_List, Q_p->NORTHWEST_quadrant_p );
			PM_INSERT( Q_List, Q_p->NORTHEAST_quadrant_p );
			PM_INSERT( Q_List, Q_p->SOUTHWEST_quadrant_p );
			PM_INSERT( Q_List, Q_p->SOUTHEAST_quadrant_p );
		}
	} // PM_INSERT List<route*>, quadrant )


	// Does any route in L clip or terminate, for either point, in R?
	List<route*> QUADRANT_CHECK( List<route*> L, quadrant* R )
	{
		List<route*> M_List;

		if( R == NULL || L.empty() )
			return M_List;

		// Brute force route intersection test done for PM1_CHECK
		if( QEDGE_INTERCEPT_CHECK( L ) )
			return M_List;

		L.begin();
		for( int i = 0; i < L.size(); i++ )
		{
			if( CLIP_SQUARE( (*L.getcurrent()), R ) ||
				PT_IN_QUADRANT( (*L.getcurrent())->StartSite(), R ) ||
				PT_IN_QUADRANT( (*L.getcurrent())->FinishSite(), R ) )
			{
				M_List += *L.getcurrent();
			}
			L.next();
		}

		return M_List;
	} // QUADRANT_CHECK

	bool PT_IN_QUADRANT( site* S, quadrant* R )
	{// Is site inside quadrant, including border
		if( R == NULL || S == NULL )
			return false;

		if( PT_IN_SQUARE( S, R ) ||
			PT_IN_BORDER( S, R ) )
			return true;

		return false;
	} // PT_IN_QUADRANT

	List<route*> CLIP_LINES( List<route*> L, quadrant* R )
	{
		List<route*> P;

		if( L.empty() || R == NULL )
			return P;
		else
		{
			L.begin();
			for( int i = 0; i < L.size(); i++ )
			{
				if( CLIP_SQUARE( *L.getcurrent(), R ) )
					P += *L.getcurrent();
				L.next();
			}
		}
		return P;
	} // CLIP_LINES


	// Is point P inside the border of quadrant Q?
	bool PT_IN_SQUARE( site* P, quadrant* Q )
	{
		if( P == NULL || Q == NULL )
			return false;

		return ((P->point().xaxis() > Q->NW_TwoDPointFloat.xaxis() ) &&
				(P->point().yaxis() < Q->NW_TwoDPointFloat.yaxis() ) &&
				(P->point().xaxis() < Q->SE_TwoDPointFloat.xaxis() ) &&
				(P->point().yaxis() > Q->SE_TwoDPointFloat.yaxis() ) );
	}

	// Determine if edge list L  and Quadrant Q_p
	// form a valid pm1 quadtree node
	// 1. One point, i.e., Start == Finish alone
	// 2. One edge clipping the quadrant
	// 3. One or more points, all sides closed,
	//    with many edges ending in Q_p but none
	//    ending in a different point within Q_p.
	// 4. edges can not intersect.
	bool PM1_CHECK( List<route*> L, quadrant* Q_p )
	{// proofed 19-MAY-2000
		site* Site_p = NULL; // site in quadrant
		route* Qedge_p = NULL; // Qedge that clips quadrant

		if( Q_p == NULL || L.empty() )
			return false;

		if( Intercept_bool )
			return false;

		// Brute force route intersection test done for PM1_CHECK
		if( QEDGE_INTERCEPT_CHECK( L ) )
			return false;

		L.begin();
		for( int i = 0; i < L.size(); i++ )
		{
		// Case 1: Isolated route
			if( ISOLATED_ELEMENT( *L.getcurrent() ) ) // isolated check
			{// Is it inside the square?
				if( Qedge_p || Site_p )
					return false;

				if( PT_IN_SQUARE( (*L.getcurrent())->StartSite(), Q_p ) )
					Site_p = (*L.getcurrent())->StartSite();
				else
					return false;

				L.next();
				continue;
			} // end Case 1
			else if( CLIP_SQUARE( (*L.getcurrent()), Q_p ) )
			{
				if( PT_IN_SQUARE( (*L.getcurrent())->StartSite(), Q_p ) == false &&
					PT_IN_SQUARE( (*L.getcurrent())->FinishSite(), Q_p ) == false )
				{
					if( Qedge_p || Site_p )
						return false;

					Qedge_p = *L.getcurrent();

					L.next();
					continue;
				}

				if( PT_IN_SQUARE( (*L.getcurrent())->StartSite(), Q_p ) )
				{
					if( Qedge_p )
						return false;

					if( Site_p == NULL )
						Site_p = (*L.getcurrent())->StartSite();

					if( Site_p != (*L.getcurrent())->StartSite() )
						return false;

					L.next();
					continue;
				}

				if( PT_IN_SQUARE( (*L.getcurrent())->FinishSite(), Q_p ) )
				{
					if( Qedge_p )
						return false;

					if( Site_p == NULL )
						Site_p = (*L.getcurrent())->FinishSite();

					if( Site_p != (*L.getcurrent())->FinishSite() )
						return false;

					L.next();
					continue;
				}
			} // border cases
			else if( PT_IN_BORDER( (*L.getcurrent())->FinishSite(), Q_p ) == false &&
					 PT_IN_BORDER( (*L.getcurrent())->StartSite(), Q_p ) == false )
				return false;

			L.next();
		} // for

		return true;
	} // PM1_CHECK

	void SPLIT_PM_NODE( quadrant* Q_p )
	{
		Q_p->TYPE = GRAY;
		repartitionNORTHWEST( Q_p );
		repartitionNORTHEAST( Q_p );
		repartitionSOUTHWEST( Q_p );
		repartitionSOUTHEAST( Q_p );
	}

	void PM_DELETE( List<route*> P, quadrant* R )
	{
		if( R == NULL )
			return;

		List<route*> L = CLIP_LINES( P, R );

		if( L.empty() )
			return;

		if( R->TYPE == GRAY )
		{
			PM_DELETE( L, R->NORTHEAST_quadrant_p );
			PM_DELETE( L, R->NORTHWEST_quadrant_p );
			PM_DELETE( L, R->SOUTHWEST_quadrant_p );
			PM_DELETE( L, R->SOUTHEAST_quadrant_p );

			L.erase();
			if( TRY_TO_MERGE_PM1( R, R, L ) )
			{
				deletenode( R ); // condense all childeren
				if( L.empty() )
					R->TYPE = WHITE; //
				else
					R->TYPE = BLACK; //
				R->route_List = L;
			}
		}
	}

	bool POSSIBLE_PM1_MERGE( quadrant* P )
	{
		if( P == NULL )
			return false;

		return ( LEAF( P->NORTHEAST_quadrant_p ) ||
				 LEAF( P->NORTHWEST_quadrant_p ) ||
				 LEAF( P->SOUTHWEST_quadrant_p ) ||
				 LEAF( P->SOUTHEAST_quadrant_p ) );
	}

	bool LEAF( quadrant* Q_p )
	{
		if( Q_p == NULL )
			return false;

		return ( Q_p->TYPE == BLACK || Q_p->TYPE == WHITE );
	} // LEAF


	bool TRY_TO_MERGE_PM1( quadrant* P, quadrant* R, List<route*>& L )
	{
		if( LEAF( P ) )
		{
			L += P->route_List;
			return true;
		}
		return false;
	}


	// Does point P clip a border of Q_p?
	bool PT_IN_BORDER( site* P, quadrant* Q_p )
	{
		if( Q_p == NULL || P == NULL )
			return false;

	// North side
		if( PointOnLine( P->point(), Q_p->NW_TwoDPointFloat, Q_p->NE_TwoDPointFloat ) )
			return true;

	// South side
		if( PointOnLine( P->point(), Q_p->SW_TwoDPointFloat, Q_p->SE_TwoDPointFloat ) )
			return true;

	// East side
		if( PointOnLine( P->point(), Q_p->NE_TwoDPointFloat, Q_p->SE_TwoDPointFloat ) )
			return true;

	// West side
		if( PointOnLine( P->point(), Q_p->NW_TwoDPointFloat, Q_p->SW_TwoDPointFloat ) )
			return true;

		return false;
	}

	// Does the route R_p cross quadrant Q_p anywhere?
	bool CLIP_SQUARE( route *R_p, quadrant* Q_p )
	{
		if( R_p == NULL || Q_p == NULL )
			return false;

	// return true if the route crosses any quadrant border
	// start passing it sides
	// NORTH side
		if( intersect( R_p->StartSite()->point(),
					   R_p->FinishSite()->point(),
					   Q_p->NW_TwoDPointFloat,
					   Q_p->NE_TwoDPointFloat ) )
			return true;

	// SOUTH side
		if( intersect( R_p->StartSite()->point(),
					   R_p->FinishSite()->point(),
					   Q_p->SW_TwoDPointFloat,
					   Q_p->SE_TwoDPointFloat ) )
			return true;

	// EAST side
		if( intersect( R_p->StartSite()->point(),
					   R_p->FinishSite()->point(),
					   Q_p->NE_TwoDPointFloat,
					   Q_p->SE_TwoDPointFloat ) )
			return true;

	// WEST side
		if( intersect( R_p->StartSite()->point(),
					   R_p->FinishSite()->point(),
					   Q_p->NW_TwoDPointFloat,
					   Q_p->SW_TwoDPointFloat ) )
			return true;

		return false;
	} // CLIP_SQUARE

	bool ISOLATED_ELEMENT( route* R_p )
	{
		return ( R_p->StartSite() == R_p->FinishSite() );
	}

	// Brute force check for intercepts in route list L
	// returns true if two edges intersect, false otherwise
	bool QEDGE_INTERCEPT_CHECK( List<route*> L )
	{
		if( L.empty() )
			return false;

		List<route*> M_list = L; // copy the list
		L.begin();
		for(int i = 0; i < L.size(); i++ )
		{
			M_list.begin();
			for(int j = 0; j < M_list.size(); j++ )
			{
				if( intersect( (*L.getcurrent())->StartSite()->point(),
							   (*L.getcurrent())->FinishSite()->point(),
							   (*M_list.getcurrent())->StartSite()->point(),
							   (*M_list.getcurrent())->FinishSite()->point()) &&
							   ( (*L.getcurrent())->StartSite()->point() !=
								 (*M_list.getcurrent())->StartSite()->point() ) &&
							   ( (*L.getcurrent())->StartSite()->point() !=
								 (*M_list.getcurrent())->FinishSite()->point() ) &&
							   ( (*L.getcurrent())->FinishSite()->point() !=
								 (*M_list.getcurrent())->StartSite()->point() ) &&
							   ( (*L.getcurrent())->FinishSite()->point() !=
								 (*M_list.getcurrent())->FinishSite()->point() ) )
							   {
							   Intercept_bool = true;
								return true;
								}
				M_list.next();
			} // for
			L.next();
		} // for
		return false;
	} // QEDGE_INTERCEPT_CHECK

////////////////////////////////
	// Does point A lie on line BC?
	bool PointOnLine( TwoDPointFloat A, TwoDPointFloat B,
					  TwoDPointFloat C )
	{
	// Case 1: Line BC is Vertical.
		if( B.xaxis() == C.xaxis() )
		{
			if( A.xaxis() == B.xaxis() &&
				A.yaxis() >= B.yaxis() &&
				A.yaxis() <= C.yaxis() )
				return true;
			else
				return false;
		}
	// Case 2: Line BC is Horizontal
		if( B.yaxis() == C.yaxis() )
		{
			if( A.yaxis() == B.yaxis() &&
				A.xaxis() >= B.xaxis() &&
				A.xaxis() <= C.xaxis() )
				return true;
			else
				return false;
		}
	// Case 3: Line BC has a slope
	/*Given points P1 = (x1,y1) and P2 = (x2,y2), the parametric form for the line is:

	    x = x1 + t(x2-x1)
	    y = y1 + t(y2-y1)
	    0 <= t <= 1

	Note that t is the same in both equations. t is called the parameter.
	When t = 0 we get P1 and when t = 1 we get P2. As t varies between 0 and 1,
	we get all the other points on the line segment between P1 and P2.

	Because the slope of a vertical line is infinite, the slope-intercept form
	of a line can be troublesome to work with. Since the parametric form does not
	involve slope, vertical lines are no longer a special case.
	So if the point is (x3,y3)
		t1 = x3-x1 / x2-x1
		t2 = y3-y1 / y2-y1
	if 0 <= t1 and t2 <= 1 the point is inside the line.*/

		float t1 = 0.0;
		float t1_D = A.xaxis() - B.xaxis();
		float t1_N = C.xaxis() - B.xaxis();
		float t2 = 0.0;
		float t2_D = A.yaxis() - B.yaxis();
		float t2_N = C.yaxis() - B.yaxis();

		t1 = t1_D / t1_N;

		t2 = t2_D / t2_N;

		return ( 0 <= t1 &&
				 t1 <= 1 &&
				 0 <= t2 &&
				 t2 <= 1 );
	}

	// Does the line between A and B intersect with the line from C to D?
	bool intersect( TwoDPointFloat A, TwoDPointFloat B,
					TwoDPointFloat C, TwoDPointFloat D )
	{
		if( A.xaxis() == 618 && A.yaxis() ==640 &&
			C.xaxis() == 512 && C.yaxis() == 640 )
			int Too = 0;


		float k = (B.xaxis()-A.xaxis())*(D.yaxis() -C.yaxis())-(B.yaxis()-A.yaxis())*(D.xaxis()-C.xaxis());
		if( k == 0.0 )
			return false;

		float r = ((A.yaxis()-C.yaxis())*(D.xaxis()-C.xaxis())-(A.xaxis()-C.xaxis())*(D.yaxis()-C.yaxis()))/k;
		float s = ((A.yaxis()-C.yaxis())*(B.xaxis()-A.xaxis())-(A.xaxis()-C.xaxis())*(B.yaxis()-A.yaxis()))/k;


		if( r<0 || r>1 || s<0 || s>1 )
			return false;

		float XI = A.xaxis() + r *(B.xaxis() - A.xaxis());
		float YI = A.yaxis() + r * (B.yaxis() - A.yaxis());

		if ((XI == A.xaxis() && YI == A.yaxis()) || (XI == B.xaxis() && YI == B.yaxis()))
			return false;

		return true;
	} // intersect

///// end PM methods
	void repartitionRoot()
	{
		root_quadrant_p = new quadrant;

		initquadrant( root_quadrant_p );
		root_quadrant_p->width_float = MaxAxis/2;

		// center
		root_quadrant_p->Center_TwoDPointFloat.xaxis( MaxAxis/2 );
		root_quadrant_p->Center_TwoDPointFloat.yaxis( MaxAxis/2 );

		// corners
		SetCorners( root_quadrant_p );

		root_quadrant_p->TYPE = WHITE;
		root_quadrant_p->DIRECTION = ROOT;
	} // repartitionRoot()

	void repartitionNORTHEAST( quadrant* Q_p )
	{
		Q_p->NORTHEAST_quadrant_p = new quadrant;

		initquadrant( Q_p->NORTHEAST_quadrant_p );

		Q_p->NORTHEAST_quadrant_p->TYPE = WHITE;
		Q_p->NORTHEAST_quadrant_p->DIRECTION = NORTHEAST;
		Q_p->NORTHEAST_quadrant_p->Parent_p = Q_p;

		// set the new X axis
		// half the distance to the parent square border
		Q_p->NORTHEAST_quadrant_p->Center_TwoDPointFloat.xaxis( Q_p->Center_TwoDPointFloat.xaxis() + Q_p->width_float/2 );

		// set the new Y axis
		Q_p->NORTHEAST_quadrant_p->Center_TwoDPointFloat.yaxis( Q_p->Center_TwoDPointFloat.yaxis() + Q_p->width_float/2 );

		// set width
		Q_p->NORTHEAST_quadrant_p->width_float = Q_p->width_float/2;

		// corners
		SetCorners( Q_p->NORTHEAST_quadrant_p );
	} // repartitionNORTHEAST( quadrant *S )

	void repartitionNORTHWEST( quadrant* Q_p )
	{
		Q_p->NORTHWEST_quadrant_p = new quadrant;

		initquadrant( Q_p->NORTHWEST_quadrant_p );

		Q_p->NORTHWEST_quadrant_p->TYPE = WHITE;
		Q_p->NORTHWEST_quadrant_p->DIRECTION = NORTHWEST;
		Q_p->NORTHWEST_quadrant_p->Parent_p = Q_p;

		// set the new X axis
		Q_p->NORTHWEST_quadrant_p->Center_TwoDPointFloat.xaxis( Q_p->Center_TwoDPointFloat.xaxis() - Q_p->width_float/2 );

		// set the new Y axis
		Q_p->NORTHWEST_quadrant_p->Center_TwoDPointFloat.yaxis( Q_p->Center_TwoDPointFloat.yaxis()  + Q_p->width_float/2 );

		// set width
		Q_p->NORTHWEST_quadrant_p->width_float = Q_p->width_float/2;

		// corners
		SetCorners( Q_p->NORTHWEST_quadrant_p );
	} // repartitionNORTHWEST

	void repartitionSOUTHWEST( quadrant* Q_p )
	{
		Q_p->SOUTHWEST_quadrant_p = new quadrant;

		initquadrant( Q_p->SOUTHWEST_quadrant_p );

		Q_p->SOUTHWEST_quadrant_p->TYPE = WHITE;
		Q_p->SOUTHWEST_quadrant_p->DIRECTION = SOUTHWEST;
		Q_p->SOUTHWEST_quadrant_p->Parent_p = Q_p;

		// set the new X axis
		Q_p->SOUTHWEST_quadrant_p->Center_TwoDPointFloat.xaxis( Q_p->Center_TwoDPointFloat.xaxis() - Q_p->width_float/2 );

		// set the new Y axis
		Q_p->SOUTHWEST_quadrant_p->Center_TwoDPointFloat.yaxis( Q_p->Center_TwoDPointFloat.yaxis() - Q_p->width_float/2 );

		// set width
		Q_p->SOUTHWEST_quadrant_p->width_float = Q_p->width_float/2;

		// corners
		SetCorners( Q_p->SOUTHWEST_quadrant_p );

	} // repartitionSW

	void repartitionSOUTHEAST( quadrant* Q_p )
	{
		Q_p->SOUTHEAST_quadrant_p = new quadrant;

		initquadrant( Q_p->SOUTHEAST_quadrant_p );

		Q_p->SOUTHEAST_quadrant_p->TYPE = WHITE;
		Q_p->SOUTHEAST_quadrant_p->DIRECTION = SOUTHEAST;
		Q_p->SOUTHEAST_quadrant_p->Parent_p = Q_p;

		// set the new X axis
		Q_p->SOUTHEAST_quadrant_p->Center_TwoDPointFloat.xaxis( Q_p->Center_TwoDPointFloat.xaxis() +  Q_p->width_float/2 );

		// set the new Y axis
		Q_p->SOUTHEAST_quadrant_p->Center_TwoDPointFloat.yaxis( Q_p->Center_TwoDPointFloat.yaxis() - Q_p->width_float/2 );

		// set width
		Q_p->SOUTHEAST_quadrant_p->width_float = Q_p->width_float/2;

		// corners
		SetCorners( Q_p->SOUTHEAST_quadrant_p );

	} // repartitionSOUTHEAST

	bool searchtree( route *S_p, quadrant* C_p )
	{
	// if NULL we KNOW it's not here.
		if( C_p == NULL )
			return false;
	// just in case
		if( S_p == NULL )
			return false;

		if( C_p->TYPE == WHITE )
			return false;

	// this is a data node
		if( C_p->TYPE == BLACK )
		{
			if( C_p->route_List.find( S_p ) )
				return true;

			return false;
		}

	// skip down the tree again
		switch( calcDirection( S_p->StartSite()->point(), C_p->Center_TwoDPointFloat ) )
		{
			case CENTER:
			case NORTH:
			case EAST:
			case NORTHEAST: if( C_p->NORTHEAST_quadrant_p == NULL )
								return false;
							C_p = C_p->NORTHEAST_quadrant_p;
							return searchtree( S_p, C_p );
							break;
			case NORTHWEST: if( C_p->NORTHWEST_quadrant_p == NULL )
								return false;// check the TwoDPointFloat against NORTHWEST
							C_p = C_p->NORTHWEST_quadrant_p;
							return searchtree( S_p, C_p );
							break;
			case WEST:
			case SOUTHWEST: if( C_p->SOUTHWEST_quadrant_p == NULL )
								return false;// check the TwoDPointFloat against SW
							C_p = C_p->SOUTHWEST_quadrant_p;
							return searchtree( S_p, C_p );
							break;
			case SOUTH:
			case SOUTHEAST: if( C_p->SOUTHEAST_quadrant_p == NULL )
								return false; // check the TwoDPointFloat against SOUTHEAST
							C_p = C_p->SOUTHEAST_quadrant_p;
							return searchtree( S_p, C_p );
							break;
			default:
				cout << "pm1quadtree:searchtree calcDirection failure " << endl;
		};

		return false;
	} // searchtree

	void PRINT_PATH( site *S_p, quadrant* C_p )
	{
		if( C_p == NULL || S_p == NULL )
			return;

		if( C_p->TYPE == BLACK )
		{
			cout << S_p->name() << endl;
			return;
		}

		switch( calcDirection( S_p->point(), C_p->Center_TwoDPointFloat ) )
		{
			case CENTER: cout << "NE";
				break;
			case NORTH: cout << "NE";
				break;
			case EAST:
				break;
			case NORTHEAST: cout << "NE";
				C_p = C_p->NORTHEAST_quadrant_p;
				break;
			case NORTHWEST: cout << "NW";
				C_p = C_p->NORTHWEST_quadrant_p;
				break;
			case WEST:;
				break;
			case SOUTHWEST: cout << "SW";
				C_p = C_p->SOUTHWEST_quadrant_p;
				break;
			case SOUTH:;
				break;
			case SOUTHEAST: cout << "SE";
				C_p = C_p->SOUTHEAST_quadrant_p;
				break;
			default:
				return;
		};

		if( C_p->TYPE == GRAY )
			cout << endl;

		PRINT_PATH( S_p, C_p );
	} // PRINT_PATH

// NEAREST_BASE_SEARCH
	route NEAREST_BASE_SEARCH( quadrant* Q_p )
	{
		route R; // default weight is DBL_MAX
		if( Q_p == NULL || CIRCLE_CHECK( Q_p ) == false )
			return R;

		printSquare( Q_p );

		if( Q_p->TYPE == WHITE )
			return R;

		route RouteTemp;
		Q_p->route_List.begin();
		for( int i = 0; i < Q_p->route_List.size(); i++ )
		{
			if( ISOLATED_ELEMENT(((*Q_p->route_List.getcurrent() ))) == false )
			{// test for material site
				RouteTemp.make( Circle_Site_p, (*Q_p->route_List.getcurrent())->StartSite() );
				if( RouteTemp.getweight() < Radius_Circle_double )
				{
					R = RouteTemp;
					Radius_Circle_double = RouteTemp.getweight();
				}

				RouteTemp.make( Circle_Site_p, (*Q_p->route_List.getcurrent())->FinishSite() );
				if( RouteTemp.getweight() < Radius_Circle_double )
				{
					R = RouteTemp;
					Radius_Circle_double = RouteTemp.getweight();
				}
			}

			Q_p->route_List.next();
		} // for

		if( LEAF( Q_p ) )
			return R;

		RouteTemp = NEAREST_BASE_SEARCH( Q_p->NORTHWEST_quadrant_p );
		route NE_route = NEAREST_BASE_SEARCH( Q_p->NORTHEAST_quadrant_p );
		route SW_route = NEAREST_BASE_SEARCH( Q_p->SOUTHWEST_quadrant_p );
		route SE_route = NEAREST_BASE_SEARCH( Q_p->SOUTHEAST_quadrant_p );

		if( RouteTemp.getweight() < R.getweight() )
			R = RouteTemp;

		if( NE_route.getweight() < R.getweight() )
			R = NE_route;

		if( SW_route.getweight() < R.getweight() )
			R = SW_route;

		if( SE_route.getweight() < R.getweight() )
			R = SE_route;

		if( R.getweight() < Radius_Circle_double )
			Radius_Circle_double = R.getweight();

		return R;
	} // NEAREST_BASE_SEARCH

// NEAREST_ROUTE_SEARCH
	void NEAREST_ROUTE_SEARCH( quadrant* Q_p )
	{
		if( Q_p == NULL || CIRCLE_CHECK( Q_p ) == false )
			return;

		printSquare( Q_p );

		if( Q_p->TYPE == WHITE )
			return;

		route* RouteTemp = NULL;
		double Temp_double = DBL_MAX;
		Q_p->route_List.begin();
		for( int i = 0; i < Q_p->route_List.size(); i++ )
		{
			if( ISOLATED_ELEMENT(((*Q_p->route_List.getcurrent() ))) == false )
			{// test for material site
				Temp_double = PT_TO_EDGE( (*Q_p->route_List.getcurrent()) );
				if( Temp_double < Radius_Circle_double )
				{
					Radius_Circle_double = Temp_double;
					NearestRoute_route_p = (*Q_p->route_List.getcurrent());
				}
			}
			Q_p->route_List.next();
		} // for

		if( LEAF( Q_p ) )
			return;

		NEAREST_ROUTE_SEARCH( Q_p->NORTHWEST_quadrant_p );
		NEAREST_ROUTE_SEARCH( Q_p->NORTHEAST_quadrant_p );
		NEAREST_ROUTE_SEARCH( Q_p->SOUTHWEST_quadrant_p );
		NEAREST_ROUTE_SEARCH( Q_p->SOUTHEAST_quadrant_p );

	} // NEAREST_ROUTE_SEARCH

	List<string> RADIUS_SITE_SEARCH( quadrant* Q_p )
	{
		List<string> R; // default weight is DBL_MAX
		if( Q_p == NULL || CIRCLE_CHECK( Q_p ) == false )
			return R;

		printSquare( Q_p );
		route RouteTemp;
		Q_p->route_List.begin();
		for( int i = 0; i < Q_p->route_List.size(); i++ )
		{
			if( (*Q_p->route_List.getcurrent())->StartSite()->material() == false )
			{
				RouteTemp.make( Circle_Site_p, (*Q_p->route_List.getcurrent())->StartSite() );
				if( RouteTemp.getweight() <= Radius_Circle_double )
				{
					R += (*Q_p->route_List.getcurrent())->StartName();
				}

				RouteTemp.make( Circle_Site_p, (*Q_p->route_List.getcurrent())->FinishSite() );
				if( RouteTemp.getweight() <= Radius_Circle_double )
				{
					R += (*Q_p->route_List.getcurrent())->FinishName();
				}
			}
			Q_p->route_List.next();
		} // for

		R += RADIUS_SITE_SEARCH( Q_p->NORTHWEST_quadrant_p );
		R += RADIUS_SITE_SEARCH( Q_p->NORTHEAST_quadrant_p );
		R += RADIUS_SITE_SEARCH( Q_p->SOUTHWEST_quadrant_p );
		R += RADIUS_SITE_SEARCH( Q_p->SOUTHEAST_quadrant_p );

		return R;
	} // RADIUS_SITE_SEARCH

	bool CIRCLE_CHECK( quadrant* Q_p )
	{// add CLIP_CIRCLE later
		if( Q_p == NULL )
			return false;

		return (QUADRANT_IN_CIRCLE( Q_p ) || QUADRANT_CLIP_CIRCLE( Q_p ));
	} // CIRCLE_CHECK

	bool QUADRANT_CLIP_CIRCLE( quadrant* Q_p )
	{
		if( Q_p == NULL || Circle_Site_p == NULL)
			return false;

		if( Radius_Circle_double == DBL_MAX )
			return true;

		return ( CLIP_CIRCLE( Q_p->NW_TwoDPointFloat, Q_p->NE_TwoDPointFloat ) || // north face
				 CLIP_CIRCLE( Q_p->SW_TwoDPointFloat, Q_p->SE_TwoDPointFloat ) || // south face
				 CLIP_CIRCLE( Q_p->SE_TwoDPointFloat, Q_p->NE_TwoDPointFloat ) || // east face
				 CLIP_CIRCLE( Q_p->NW_TwoDPointFloat, Q_p->SW_TwoDPointFloat ) ); // west face
	} // QUADRANT_CLIP_CIRCLE

	bool CLIP_CIRCLE( TwoDPointFloat One, TwoDPointFloat Two )
	{
		if( Circle_Site_p == NULL)
			return false;

		if( Radius_Circle_double == DBL_MAX )
			return true;

		double M = SLOPE( Two, One );

		double b = (Two.yaxis() + One.yaxis() - M * Two.xaxis() - M * One.xaxis()) / 2;

		double h =  Circle_Site_p->xaxis();

		double k =  Circle_Site_p->yaxis();

		double A = 1 + M * M;

		double B = -2*(h - b*M + k * M );

		double C = -2*b*k + b*b + k*k - Radius_Circle_double * Radius_Circle_double + h * h;

		double SA = B*B - 4*(A * C );

		if( SA < 0 )
			return false;
		else
			return true;

//		double SA_sqrt = sqrt( SA );
//		double PositiveSide = 0.0;
//		PositiveSide = (-B + SA_sqrt ) / 2*A;
//		double NegativeSide = 0.0;
//		NegativeSide = (-B - SA_sqrt ) / 2*A;
//
//		return true;
	} // CLIP_CIRCLE

	bool QUADRANT_IN_CIRCLE( quadrant* Q_p )
	{// Is quadrant Q_p within the circle defined
	 // by Circle_Site_p and Radius_Circle_double?

		if( Q_p == NULL || Circle_Site_p == NULL)
			return false;

		if( Radius_Circle_double == DBL_MAX )
			return true;

		route NW_route;
		route NE_route;
		route SW_route;
		route SE_route;

		site NW_site;
		site NE_site;
		site SW_site;
		site SE_site;

		NW_site.xaxis( Q_p->NW_TwoDPointFloat.xaxis() );
		NW_site.yaxis( Q_p->NW_TwoDPointFloat.yaxis() );
		NE_site.xaxis( Q_p->NE_TwoDPointFloat.xaxis() );
		NE_site.yaxis( Q_p->NE_TwoDPointFloat.yaxis() );
		SW_site.xaxis( Q_p->SW_TwoDPointFloat.xaxis() );
		SW_site.yaxis( Q_p->SW_TwoDPointFloat.yaxis() );
		SE_site.xaxis( Q_p->SE_TwoDPointFloat.xaxis() );
		SE_site.yaxis( Q_p->SE_TwoDPointFloat.yaxis() );

		NW_route.make( Circle_Site_p, &NW_site );
		NE_route.make( Circle_Site_p, &NE_site );
		SW_route.make( Circle_Site_p, &SW_site );
		SE_route.make( Circle_Site_p, &SE_site );

		return ( NW_route.getweight() <= Radius_Circle_double &&
				 NE_route.getweight() <= Radius_Circle_double &&
				 SW_route.getweight() <= Radius_Circle_double &&
				 SE_route.getweight() <= Radius_Circle_double );
	} // QUADRANT_IN_CIRCLE

	double SLOPE( TwoDPointFloat One, TwoDPointFloat Two )
	{
		double M = One.xaxis() - Two.xaxis();
		if( M == 0.0 )
			return 0.0;
		else
			return (One.yaxis() - Two.yaxis()) / M;
	} // SLOPE

	double PT_TO_EDGE( route* R )
	{
		if( R == NULL || Circle_Site_p == NULL )
			return DBL_MAX;

		double R0 = PT_TO_EDGE_RIGHTANGLE( R );
		double R1 = DBL_MAX;
		double R2 = DBL_MAX;
		route R1_route;
		route R2_route;

		// that didn't work so try each end point
		R1_route.make( Circle_Site_p, R->FinishSite() );
		R2_route.make( Circle_Site_p, R->StartSite() );

		R1 = R1_route.getweight();
		R2 = R2_route.getweight();
		if( R1 < R2 &&
			R1 < R0 )
			return R1;

		if( R2 < R1 &&
			R2 < R0 )
			return R2;

		return R0;
	} // PT_TO_EDGE

	double PT_TO_EDGE_RIGHTANGLE( route* R )
	{
		if( R == NULL || Circle_Site_p == NULL )
			return DBL_MAX;

		double M = SLOPE( R->StartSite()->point(), R->FinishSite()->point() );
		double M_inv = 0.0;
		if( M == 0.0 )
			M_inv = 0.0;
		else
			M_inv = -(1 / M);

		double b = R->StartSite()->yaxis() - M*R->StartSite()->xaxis();
		double c = Circle_Site_p->yaxis() - M_inv*Circle_Site_p->xaxis();
		double x1 = ( c - b ) / ( M - M_inv );
		double y1 = M_inv * x1 + c;

		site S_new;

		try
		{
			S_new.xaxis( x1 );
			S_new.yaxis( y1 );
		}
		catch( BadAxisException e )
		{// no closest approach
			return DBL_MAX;
		}


		float t = (x1 - R->StartSite()->xaxis()) / (R->FinishSite()->xaxis() - R->StartSite()->xaxis() );
		route R_new;
		R_new.make( Circle_Site_p, &S_new );

		if( t > 0.0 && t < 1.0 )
			return R_new.getweight();
		else
			return DBL_MAX;
	} // PT_TO_EDGE_RIGHTANGLE

public:

	pm1quadtree()
	{
		root_quadrant_p = NULL;
		Radius_Circle_double = DBL_MIN;
		Circle_Site_p = NULL;
		Intercept_bool = false;
		NearestRoute_route_p = NULL;
	} // pm1quadtree

	pm1quadtree( const pm1quadtree &P )
	{
		if( P.root_quadrant_p )
		{
			root_quadrant_p = new quadrant;
			root_quadrant_p = P.root_quadrant_p;
			// recursivly fill the tree
			copynode( root_quadrant_p, P.root_quadrant_p);
		}
		else
		{
			root_quadrant_p = NULL;
		}
		Intercept_bool = P.Intercept_bool;
		Radius_Circle_double = P.Radius_Circle_double;
		Circle_Site_p = P.Circle_Site_p;
		NearestRoute_route_p = P.NearestRoute_route_p;
	}

	~pm1quadtree()
	{
		deletenode( root_quadrant_p );
		root_quadrant_p = NULL;
		Circle_Site_p = NULL;
		Intercept_bool = false;
		NearestRoute_route_p = NULL;
	} // ~pm1quadtree

	bool empty()
	{
		return ( root_quadrant_p == NULL );
	} // empty

	void erase()
	{
		deletenode( root_quadrant_p );
		root_quadrant_p = NULL;
		Radius_Circle_double = FLT_MIN;
		Intercept_bool = false;
	} // erase

	bool find( site* S_p )
	{
		route R_route;
		R_route.make( S_p, S_p );
		quadrant* C_p = root_quadrant_p;

		return searchtree( &R_route, C_p );
	} // find( TwoDPointFloat S )

	bool ADD_QEDGE( route *R )
	{
		List<route*> R_List;
		if( R == NULL )
		{
			cout << "ERROR: Given qedge is void." << endl;
			return false;
		}

		R_List.AddToList( R );
		quadrant* Q_p = root_quadrant_p;
		PM_INSERT( R_List, Q_p );

		if( Intercept_bool )
		{
			cout << "ERROR: intersecting routes found. Clearing PM1 Quadtree." << endl;
			erase();
			return false;
		}
		return true;
	} // insert

	void LIST_PATH( site *S_p )
	{
		if( find( S_p ) )
			PRINT_PATH( S_p, root_quadrant_p );
	} // LIST_PATH( site * )

	void PRINT_PMTREE() // print all nodes in the tree
	{
	// Print in order NW, NE, SW, SE corners, then quadrant
		if( empty() )
		{
			cout << "ERROR: The PM1Quadtree is empty." << endl;
			return;
		}

		cout << "PM Quadtree listing:" << endl;
		printQuadrant( root_quadrant_p );
	}

	void RECTANGLE_SITES( TwoDPointFloat NW_point,
						  TwoDPointFloat NE_point,
						  TwoDPointFloat SW_point,
						  TwoDPointFloat SE_point )
	{
	// Set reference corners
		if( empty() )
		{
			cout << "ERROR: The PM1Quadtree is empty." << endl;
			return;
		}

		quadrant RectangleHolder;
		initquadrant( &RectangleHolder );

		RectangleHolder.NW_TwoDPointFloat = NW_point;
		RectangleHolder.NE_TwoDPointFloat = NE_point;
		RectangleHolder.SE_TwoDPointFloat = SE_point;
		RectangleHolder.SW_TwoDPointFloat = SW_point;

		cout << "Traversing PM Quadtree..." << endl;

		List<string> L_list;
		L_list = RECTANGLE_SEARCH( root_quadrant_p, &RectangleHolder );
		cout << endl;

		if( L_list.empty() == false )
		{
			cout << "List of sites within the region:" << endl;
			L_list.begin();
			for( int i = 0; i < L_list.size(); i++ )
			{
				cout << (*L_list.getcurrent()) << " ";
				L_list.next();
			}
		}
		else
			cout << "ERROR: No sites in the region" << endl;

	} // RECTANGLE_SITES

	void NEAREST_BASE( site* S )
	{
		if( S == NULL )
			return;

		if( empty() )
		{
			cout << "ERROR: The PM1Quadtree is empty." << endl;
			return;
		}

		route R;

		Radius_Circle_double = DBL_MAX;
		Circle_Site_p = S;

		cout << "Traversing PM Quadtree..." << endl;

		R = NEAREST_BASE_SEARCH( root_quadrant_p );
		Radius_Circle_double = DBL_MAX;
		if( R.getweight() == DBL_MAX )
		{
			cout << "ERROR: No base site was found." << endl;
			return;
		}

		cout << "Nearest base site to raw material site " << S->name();
		cout << " is ";

		cout << R.FinishName();
		cout << " of type ";
		cout << R.FinishSite()->type();
		cout << " at a distance ";
		R.RouteWeight().OutPut();
		cout << "." << endl;
	} // NEAREST_BASE

	List<string> RADIUS_SITES( site* Target_Site_p, double Site_Radius_double )
	{
		List<string> S_List;
		Radius_Circle_double = Site_Radius_double;
		Circle_Site_p = Target_Site_p;

		if( empty() )
		{
			cout << "ERROR: The PM1Quadtree is empty." << endl;
			return S_List;
		}

		if( Target_Site_p == NULL )
		{
			cout << "ERROR: No valid site given." << endl;
			return S_List;
		}

	   	cout << "Traversing PM Quadtree ..." << endl;
		S_List = RADIUS_SITE_SEARCH( root_quadrant_p );

		// clean up
		Radius_Circle_double = DBL_MAX;
		Circle_Site_p = NULL;
		return S_List;
	} // RADIUS_SITES

	void NEAREST_ROUTE( site* S )
	{
		if( S == NULL )
		{
			cout << "ERROR: An invalid site was given." << endl;
			return;
		}

		if( empty() )
		{
			cout << "ERROR: The PM1Quadtree is empty." << endl;
			return;
		}

		NearestRoute_route_p = NULL;

		Radius_Circle_double = DBL_MAX;
		Circle_Site_p = S;

		cout << "Traversing PM Quadtree..." << endl;

		NEAREST_ROUTE_SEARCH( root_quadrant_p );
		if( NearestRoute_route_p == NULL )
		{
			cout << "ERROR: No base site was found." << endl;
			return;
		}

		cout << "Nearest route to raw material site " << S->name();
		cout << " is ";
		cout << NearestRoute_route_p->StartName();
		cout << " <--> ";
		cout << NearestRoute_route_p->FinishName();
		cout << " at a distance " << Radius_Circle_double << "." << endl;

		// cleanup
		Radius_Circle_double = DBL_MAX;
		NearestRoute_route_p = NULL;
		cout << endl;
	} // NEAREST_ROUTE
};
#endif // PM1QUADTREE_CC