#include ".\Coonssdfinterp.h"

using namespace SDF;

CoonsSdfInterp::CoonsSdfInterp(void)
{
	this->reset();
}

CoonsSdfInterp::~CoonsSdfInterp(void)
{
}

void CoonsSdfInterp::setNode(octreeNode *node) {
	this->reset();
	this->node = node;
}

double CoonsSdfInterp::interpolate(const Point &from) {
	return (double)findCOONSInterpDist(from);
}

void CoonsSdfInterp::reset() {
	for(int i = 0; i < 8; ++i) {
		if(this->edges[i] != NULL) {
			delete this->edges[i];
			this->edges[i] = NULL;
		}
	}
}

/*
 * Most of the code below is copied from Young's OctreeADF.cpp
 */

float CoonsSdfInterp::findCOONSInterpDist(float u, float v, float w){
	float term1 = (1-u)*evalP0(v,w)+ u*evalP1(v,w)+ (1-v)*evalQ0(u,w)+ v*evalQ1(u,w) + (1-w)*evalR0(u,v)+ w*evalR1(u,v); 
	
	float term2 = (1-u)*(1-w)*evalP0(v,0) + u*(1-w)*evalP1(v,0)+ (1-v)*(1-w)*evalQ0(u,0)+ v*(1-w)*evalQ1(u,0)
		         +(1-u)*w*evalP0(v,1)+ u*w*evalP1(v,1)+ (1-v)*w*evalQ0(u,1)+ v*w*evalQ1(u,1)
				 +(1-u)*(1-v)*evalP0(0,w)+ +(1-u)*v*evalP0(1,w)+ u*(1-v)*evalP1(0,w) + u*v*evalP1(1,w);
	
	float term3 = 	node->m_d000 * (1 - u) * (1 - v) * (1 - w) + node->m_d100 * u * (1 - v) * (1 - w) +
					node->m_d010 * (1 - u) * v * (1 - w) + node->m_d001 * (1 - u) * (1 - v) * w +
					node->m_d101 * u * (1 - v) * w + 	node->m_d011 * (1 - u) * v * w  +	node->m_d110 * u * v * (1 - w) +
					node->m_d111 * u * v * w;
//	printf("term1 = %f term2 = %f \n",term1,term2);
	return term1- term2 + term3;
}


float CoonsSdfInterp::findCOONSInterpDist(const Point& pt){
	float u,v,w;
	u = (pt(0) - node->m_v[0]) / node->m_edgeLength;
	v = (pt(1) - node->m_v[1]) / node->m_edgeLength;
	w = (pt(2) - node->m_v[2]) / node->m_edgeLength;
	return findCOONSInterpDist(u, v, w);
}

bool CoonsSdfInterp::isEdgeExtra(octreeEdge* e, octreeNode* no, octreeNode::Direction d1, octreeNode::Direction d2 ){
	octreeNode* n1 = no->findNeighbor(d1);
	if(n1 && n1->m_edgeLength > no->m_edgeLength){
		Point sp= e->getStartPoint();
		Point ep= e->getEndPoint();
		sp[3] = 1.0;
		ep[3] = 1.0;
		if(sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,0.0,n1->m_edgeLength/2.0f,0.0f))
		 ||sp == (n1->m_v + Point(0.0,n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,n1->m_edgeLength,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f, n1->m_edgeLength, n1->m_edgeLength/2.0f,0.0f))
	     ||sp == (n1->m_v + Point(n1->m_edgeLength, n1->m_edgeLength/2.0f, n1->m_edgeLength/2.0f,0.0f))
		)
		 return true;
		
		else if(ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,0.0,n1->m_edgeLength/2.0f,0.0f))
		 ||ep == (n1->m_v + Point(0.0,n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,n1->m_edgeLength,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f, n1->m_edgeLength, n1->m_edgeLength/2.0f,0.0f))
	     ||ep == (n1->m_v + Point(n1->m_edgeLength, n1->m_edgeLength/2.0f, n1->m_edgeLength/2.0f,0.0f))
		)
		return true;
	}	
	n1 = no->findNeighbor(d2);
	if(n1 && n1->m_edgeLength > no->m_edgeLength){
		Point sp= e->getStartPoint();
		Point ep= e->getEndPoint();
		sp[3] = 1.0;
		ep[3] = 1.0;
		if(sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,0.0,n1->m_edgeLength/2.0f,0.0f))
		 ||sp == (n1->m_v + Point(0.0,n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,n1->m_edgeLength,0.0f))
		 ||sp == (n1->m_v + Point(n1->m_edgeLength/2.0f, n1->m_edgeLength, n1->m_edgeLength/2.0f,0.0f))
	     ||sp == (n1->m_v + Point(n1->m_edgeLength, n1->m_edgeLength/2.0f, n1->m_edgeLength/2.0f,0.0f))
		)
		 return true;
		
		else if(ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,0.0,n1->m_edgeLength/2.0f,0.0f))
		 ||ep == (n1->m_v + Point(0.0,n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f,n1->m_edgeLength/2.0f,n1->m_edgeLength,0.0f))
		 ||ep == (n1->m_v + Point(n1->m_edgeLength/2.0f, n1->m_edgeLength, n1->m_edgeLength/2.0f,0.0f))
	     ||ep == (n1->m_v + Point(n1->m_edgeLength, n1->m_edgeLength/2.0f, n1->m_edgeLength/2.0f,0.0f))
		)
		return true;
	}
	return false;
}


void CoonsSdfInterp::processOctreeEdges(octreeNode *no){
	//octreeEdge* e = new octreeEdge(Point(50,-50,0),Point(100,-50,0));
	//isEdgeExtra(e,m_vecNodes[12],octreeNode::D,octreeNode::F);
	//for edge 0 (shared by bottom and front plane)
	float _a,_b,_c;
	this->edges[0] = new octreeEdge(no->m_v, no->m_v+ Point(no->m_edgeLength,0,0));
	this->edges[0]->m_nDir[0] = octreeNode::D;
	this->edges[0]->m_nDir[1] = octreeNode::F;

	//if(getRootNode()->isEdgeOnBoundary(this->edges[0]) == false && 
	 if(!isEdgeExtra(this->edges[0],no,octreeNode::D,octreeNode::F)){
		
		octreeNode* nD = no->findNeighbor(octreeNode::D);
		octreeNode* nF = no->findNeighbor(octreeNode::F);
		octreeNode* nDF = NULL;
		if(nD) nDF = nD->findNeighbor(octreeNode::F);
		if(nD && nD->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nD->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nD->m_d001,no->m_d000,no->m_d100,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d100,nD->m_d101,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nD && nD->isLeafNode() == false){
			float dM = nD->subnodeArray[0][0][1]->m_d101;
			this->edges[0]->computeABCFromSamples(no->m_d000,dM,no->m_d100);
			this->edges[0]->m_edgeType = octreeEdge::LONG;
		}
		else if(nF && nF->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nF->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nF->m_d010,no->m_d000,no->m_d100,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d100,nF->m_d110,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nF && nF->isLeafNode() == false){
			float dM = nF->subnodeArray[0][1][0]->m_d110;
			this->edges[0]->computeABCFromSamples(no->m_d000,dM,no->m_d100);
			this->edges[0]->m_edgeType = octreeEdge::LONG;
		}

		else if(nDF && nDF->m_edgeLength > no->m_edgeLength){
			if(no->m_v[0] > nDF->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nDF->m_d011,no->m_d000,no->m_d100,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d100,nDF->m_d111,_a,_b,_c);
				this->edges[0]->setABC(_a,_b,_c);
				this->edges[0]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nDF && nDF->isLeafNode() == false){
			float dM = nDF->subnodeArray[0][1][1]->m_d111;
			this->edges[0]->computeABCFromSamples(no->m_d000,dM,no->m_d100);
			this->edges[0]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[0]->setABC(0.0,no->m_d100 - no->m_d000, no->m_d000);
			this->edges[0]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[0]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[0] > no->m_parent->m_v[0]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::L);
			dv1 = neighbor->m_d000;
			dv2 = no->m_d000;
			dv3 = no->m_d100;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[0]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::R);
			dv1 = no->m_d000;
			dv2 = no->m_d100;
			dv3 = neighbor->m_d100;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[0]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 0

	//********************************************************************************
	// edge 1
	//********************************************************************************
	this->edges[1] = new octreeEdge(no->m_v + Point(no->m_edgeLength,0,0), no->m_v+ Point(no->m_edgeLength,no->m_edgeLength,0));
	this->edges[1]->m_nDir[0] = octreeNode::D;
	this->edges[1]->m_nDir[1] = octreeNode::R;

	//if(getRootNode()->isEdgeOnBoundary(this->edges[0]) == false && 
	 if(!isEdgeExtra(this->edges[1],no,octreeNode::D,octreeNode::R)){
		
		octreeNode* nD = no->findNeighbor(octreeNode::D);
		octreeNode* nR = no->findNeighbor(octreeNode::R);
		octreeNode* nDR = NULL;
		if(nD) nDR = nD->findNeighbor(octreeNode::R);
		if(nD && nD->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nD->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nD->m_d101,no->m_d100,no->m_d110,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d110,nD->m_d111,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nD && nD->isLeafNode() == false){
			float dM = nD->subnodeArray[1][0][1]->m_d111;
			this->edges[1]->computeABCFromSamples(no->m_d100,dM,no->m_d110);
			this->edges[1]->m_edgeType = octreeEdge::LONG;
		}
		else if(nR && nR->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nR->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nR->m_d000,no->m_d100,no->m_d110,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d110,nR->m_d010,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nR && nR->isLeafNode() == false){
			float dM = nR->subnodeArray[0][0][0]->m_d010;
			this->edges[1]->computeABCFromSamples(no->m_d100,dM,no->m_d110);
			this->edges[1]->m_edgeType = octreeEdge::LONG;
		}

		else if(nDR && nDR->m_edgeLength > no->m_edgeLength){
			if(no->m_v[1] > nDR->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nDR->m_d001,no->m_d100,no->m_d110,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d110,nDR->m_d011,_a,_b,_c);
				this->edges[1]->setABC(_a,_b,_c);
				this->edges[1]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nDR && nDR->isLeafNode() == false){
			float dM = nDR->subnodeArray[0][0][1]->m_d011;
			this->edges[1]->computeABCFromSamples(no->m_d100,dM,no->m_d110);
			this->edges[1]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[1]->setABC(0.0,no->m_d110 - no->m_d100, no->m_d100);
			this->edges[1]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[1]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[1] > no->m_parent->m_v[1]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::F);
			dv1 = neighbor->m_d100;
			dv2 = no->m_d100;
			dv3 = no->m_d110;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[1]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::B);
			dv1 = no->m_d100;
			dv2 = no->m_d110;
			dv3 = neighbor->m_d110;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[1]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 1

	//********************************************************************************
	// edge 2
	//********************************************************************************
	this->edges[2] = new octreeEdge(no->m_v + Point(0,no->m_edgeLength,0), no->m_v+ Point(no->m_edgeLength,no->m_edgeLength,0));
	this->edges[2]->m_nDir[0] = octreeNode::D;
	this->edges[2]->m_nDir[1] = octreeNode::B;

	if(!isEdgeExtra(this->edges[2],no,octreeNode::D,octreeNode::B)){
		
		octreeNode* nD = no->findNeighbor(octreeNode::D);
		octreeNode* nB = no->findNeighbor(octreeNode::B);
		octreeNode* nDB = NULL;
		if(nD) nDB = nD->findNeighbor(octreeNode::B);
		if(nD && nD->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nD->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nD->m_d011,no->m_d010,no->m_d110,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d110,nD->m_d111,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nD && nD->isLeafNode() == false){
			float dM = nD->subnodeArray[0][1][1]->m_d111;
			this->edges[2]->computeABCFromSamples(no->m_d010,dM,no->m_d110);
			this->edges[2]->m_edgeType = octreeEdge::LONG;
		}
		else if(nB && nB->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nB->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nB->m_d000,no->m_d010,no->m_d110,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d110,nB->m_d100,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nB && nB->isLeafNode() == false){
			float dM = nB->subnodeArray[0][0][0]->m_d100;
			this->edges[2]->computeABCFromSamples(no->m_d010,dM,no->m_d110);
			this->edges[2]->m_edgeType = octreeEdge::LONG;
		}

		else if(nDB && nDB->m_edgeLength > no->m_edgeLength){
			if(no->m_v[0] > nDB->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nDB->m_d001,no->m_d010,no->m_d110,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d110,nDB->m_d101,_a,_b,_c);
				this->edges[2]->setABC(_a,_b,_c);
				this->edges[2]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nDB && nDB->isLeafNode() == false){
			float dM = nDB->subnodeArray[0][0][1]->m_d101;
			this->edges[2]->computeABCFromSamples(no->m_d010,dM,no->m_d110);
			this->edges[2]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[2]->setABC(0.0,no->m_d110 - no->m_d010, no->m_d010);
			this->edges[2]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[2]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[0] > no->m_parent->m_v[0]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::L);
			dv1 = neighbor->m_d010;
			dv2 = no->m_d010;
			dv3 = no->m_d110;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[2]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::R);
			dv1 = no->m_d010;
			dv2 = no->m_d110;
			dv3 = neighbor->m_d110;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[2]->setABC(aS,bS,cS);
		}

	}
	//end of process for edge 2
	
	//********************************************************************************
	// edge 3
	//********************************************************************************
	this->edges[3] = new octreeEdge(no->m_v, no->m_v+ Point(0,no->m_edgeLength,0));
	this->edges[3]->m_nDir[0] = octreeNode::D;
	this->edges[3]->m_nDir[1] = octreeNode::L;
	//if(getRootNode()->isEdgeOnBoundary(this->edges[0]) == false && 
	 if(!isEdgeExtra(this->edges[3],no,octreeNode::D,octreeNode::L)){
		
		octreeNode* nD = no->findNeighbor(octreeNode::D);
		octreeNode* nL = no->findNeighbor(octreeNode::L);
		octreeNode* nDL = NULL;
		if(nD) nDL = nD->findNeighbor(octreeNode::L);
		if(nD && nD->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nD->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nD->m_d001,no->m_d000,no->m_d010,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d010,nD->m_d011,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nD && nD->isLeafNode() == false){
			float dM = nD->subnodeArray[0][0][1]->m_d011;
			this->edges[3]->computeABCFromSamples(no->m_d000,dM,no->m_d010);
			this->edges[3]->m_edgeType = octreeEdge::LONG;
		}
		else if(nL && nL->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nL->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nL->m_d100,no->m_d000,no->m_d010,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d010,nL->m_d110,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nL && nL->isLeafNode() == false){
			float dM = nL->subnodeArray[1][0][0]->m_d110;
			this->edges[3]->computeABCFromSamples(no->m_d000,dM,no->m_d010);
			this->edges[3]->m_edgeType = octreeEdge::LONG;
		}

		else if(nDL && nDL->m_edgeLength > no->m_edgeLength){
			if(no->m_v[1] > nDL->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nDL->m_d101,no->m_d000,no->m_d010,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d010,nDL->m_d111,_a,_b,_c);
				this->edges[3]->setABC(_a,_b,_c);
				this->edges[3]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nDL && nDL->isLeafNode() == false){
			float dM = nDL->subnodeArray[1][0][1]->m_d111;
			this->edges[3]->computeABCFromSamples(no->m_d000,dM,no->m_d010);
			this->edges[3]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[3]->setABC(0.0,no->m_d010 - no->m_d000, no->m_d000);
			this->edges[3]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[3]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[1] > no->m_parent->m_v[1]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::F);
			dv1 = neighbor->m_d000;
			dv2 = no->m_d000;
			dv3 = no->m_d010;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[3]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::B);
			dv1 = no->m_d000;
			dv2 = no->m_d010;
			dv3 = neighbor->m_d010;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[3]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 3
	
	//********************************************************************************
	// edge 4
	//********************************************************************************
	this->edges[4] = new octreeEdge(no->m_v, no->m_v+ Point(0,0,no->m_edgeLength));
	this->edges[4]->m_nDir[0] = octreeNode::F;
	this->edges[4]->m_nDir[1] = octreeNode::L;
	if(!isEdgeExtra(this->edges[4],no,octreeNode::F,octreeNode::L)){
		
		octreeNode* nF = no->findNeighbor(octreeNode::F);
		octreeNode* nL = no->findNeighbor(octreeNode::L);
		octreeNode* nFL = NULL;
		if(nF) nFL = nF->findNeighbor(octreeNode::L);
		if(nF && nF->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nF->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nF->m_d010,no->m_d000,no->m_d001,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d001,nF->m_d011,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nF && nF->isLeafNode() == false){
			float dM = nF->subnodeArray[0][1][0]->m_d011;
			this->edges[4]->computeABCFromSamples(no->m_d000,dM,no->m_d001);
			this->edges[4]->m_edgeType = octreeEdge::LONG;
		}
		else if(nL && nL->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nL->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nL->m_d100,no->m_d000,no->m_d001,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d001,nL->m_d101,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nL && nL->isLeafNode() == false){
			float dM = nL->subnodeArray[1][0][0]->m_d101;
			this->edges[4]->computeABCFromSamples(no->m_d000,dM,no->m_d001);
			this->edges[4]->m_edgeType = octreeEdge::LONG;
		}

		else if(nFL && nFL->m_edgeLength > no->m_edgeLength){
			if(no->m_v[2] > nFL->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nFL->m_d110,no->m_d000,no->m_d001,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d000,no->m_d001,nFL->m_d111,_a,_b,_c);
				this->edges[4]->setABC(_a,_b,_c);
				this->edges[4]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nFL && nFL->isLeafNode() == false){
			float dM = nFL->subnodeArray[1][1][0]->m_d111;
			this->edges[4]->computeABCFromSamples(no->m_d000,dM,no->m_d001);
			this->edges[4]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[4]->setABC(0.0,no->m_d001 - no->m_d000, no->m_d000);
			this->edges[4]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[4]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[2] > no->m_parent->m_v[2]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::D);
			dv1 = neighbor->m_d000;
			dv2 = no->m_d000;
			dv3 = no->m_d001;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[4]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::U);
			dv1 = no->m_d000;
			dv2 = no->m_d001;
			dv3 = neighbor->m_d001;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[4]->setABC(aS,bS,cS);
		}

		
	}
	//end of process for edge 4
	
	//********************************************************************************
	// edge 5
	//********************************************************************************
	this->edges[5] = new octreeEdge(no->m_v + Point(no->m_edgeLength,0,0), no->m_v+ Point(no->m_edgeLength,0,no->m_edgeLength));
	this->edges[5]->m_nDir[0] = octreeNode::F;
	this->edges[5]->m_nDir[1] = octreeNode::R;
	if(!isEdgeExtra(this->edges[5],no,octreeNode::F,octreeNode::R)){
		
		octreeNode* nF = no->findNeighbor(octreeNode::F);
		octreeNode* nR = no->findNeighbor(octreeNode::R);
		octreeNode* nFR = NULL;
		if(nF) nFR = nF->findNeighbor(octreeNode::R);
		if(nF && nF->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nF->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nF->m_d110,no->m_d100,no->m_d101,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d101,nF->m_d111,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nF && nF->isLeafNode() == false){
			float dM = nF->subnodeArray[1][1][0]->m_d111;
			this->edges[5]->computeABCFromSamples(no->m_d100,dM,no->m_d101);
			this->edges[5]->m_edgeType = octreeEdge::LONG;
		}
		else if(nR && nR->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nR->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nR->m_d000,no->m_d100,no->m_d101,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d101,nR->m_d001,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nR && nR->isLeafNode() == false){
			float dM = nR->subnodeArray[0][0][0]->m_d001;
			this->edges[5]->computeABCFromSamples(no->m_d100,dM,no->m_d101);
			this->edges[5]->m_edgeType = octreeEdge::LONG;
		}

		else if(nFR && nFR->m_edgeLength > no->m_edgeLength){
			if(no->m_v[2] > nFR->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nFR->m_d010,no->m_d100,no->m_d101,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d100,no->m_d101,nFR->m_d011,_a,_b,_c);
				this->edges[5]->setABC(_a,_b,_c);
				this->edges[5]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nFR && nFR->isLeafNode() == false){
			float dM = nFR->subnodeArray[0][1][0]->m_d011;
			this->edges[5]->computeABCFromSamples(no->m_d100,dM,no->m_d101);
			this->edges[5]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[5]->setABC(0.0,no->m_d101 - no->m_d100, no->m_d100);
			this->edges[5]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[5]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[2] > no->m_parent->m_v[2]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::D);
			dv1 = neighbor->m_d100;
			dv2 = no->m_d100;
			dv3 = no->m_d101;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[5]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::U);
			dv1 = no->m_d100;
			dv2 = no->m_d101;
			dv3 = neighbor->m_d101;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[5]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 5

	
	//********************************************************************************
	// edge 6
	//********************************************************************************
	this->edges[6] = new octreeEdge(no->m_v + Point(no->m_edgeLength, no->m_edgeLength,0), no->m_v+ Point(no->m_edgeLength,no->m_edgeLength,no->m_edgeLength));
	this->edges[6]->m_nDir[0] = octreeNode::B;
	this->edges[6]->m_nDir[1] = octreeNode::R;
	if(!isEdgeExtra(this->edges[6],no,octreeNode::B,octreeNode::R)){
		
	
		octreeNode* nB = no->findNeighbor(octreeNode::B);
		octreeNode* nR = no->findNeighbor(octreeNode::R);
		octreeNode* nBR = NULL;
		if(nB) nBR = nB->findNeighbor(octreeNode::R);
		if(nB && nB->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nB->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nB->m_d100,no->m_d110,no->m_d111,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d110,no->m_d111,nB->m_d101,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nB && nB->isLeafNode() == false){
			float dM = nB->subnodeArray[1][0][0]->m_d101;
			this->edges[6]->computeABCFromSamples(no->m_d110,dM,no->m_d111);
			this->edges[6]->m_edgeType = octreeEdge::LONG;
		}
		else if(nR && nR->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nR->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nR->m_d010,no->m_d110,no->m_d111,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d110,no->m_d111,nR->m_d011,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nR && nR->isLeafNode() == false){
			float dM = nR->subnodeArray[0][1][0]->m_d011;
			this->edges[6]->computeABCFromSamples(no->m_d110,dM,no->m_d111);
			this->edges[6]->m_edgeType = octreeEdge::LONG;
		}

		else if(nBR && nBR->m_edgeLength > no->m_edgeLength){
			if(no->m_v[2] > nBR->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nBR->m_d000,no->m_d110,no->m_d111,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d110,no->m_d111,nBR->m_d001,_a,_b,_c);
				this->edges[6]->setABC(_a,_b,_c);
				this->edges[6]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nBR && nBR->isLeafNode() == false){
			float dM = nBR->subnodeArray[0][0][0]->m_d001;
			this->edges[6]->computeABCFromSamples(no->m_d110,dM,no->m_d111);
			this->edges[6]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[6]->setABC(0.0,no->m_d111 - no->m_d110, no->m_d110);
			this->edges[6]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{

		this->edges[6]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[2] > no->m_parent->m_v[2]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::D);
			dv1 = neighbor->m_d110;
			dv2 = no->m_d110;
			dv3 = no->m_d111;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[6]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::U);
			dv1 = no->m_d110;
			dv2 = no->m_d111;
			dv3 = neighbor->m_d111;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[6]->setABC(aS,bS,cS);
		}

	}
	//end of process for edge 6

	//********************************************************************************
	// edge 7
	//********************************************************************************
	this->edges[7] = new octreeEdge(no->m_v + Point(0,no->m_edgeLength,0), no->m_v+ Point(0,no->m_edgeLength,no->m_edgeLength));
	this->edges[7]->m_nDir[0] = octreeNode::B;
	this->edges[7]->m_nDir[1] = octreeNode::L;
	if(!isEdgeExtra(this->edges[7],no,octreeNode::B,octreeNode::L)){
		
		octreeNode* nB = no->findNeighbor(octreeNode::B);
		octreeNode* nL = no->findNeighbor(octreeNode::L);
		octreeNode* nBL = NULL;
		if(nB) nBL = nB->findNeighbor(octreeNode::L);
		if(nB && nB->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nB->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nB->m_d000,no->m_d010,no->m_d011,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d011,nB->m_d001,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nB && nB->isLeafNode() == false){
			float dM = nB->subnodeArray[0][0][0]->m_d001;
			this->edges[7]->computeABCFromSamples(no->m_d010,dM,no->m_d011);
			this->edges[7]->m_edgeType = octreeEdge::LONG;
		}
		else if(nL && nL->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[2] > nL->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nL->m_d110,no->m_d010,no->m_d011,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d011,nL->m_d111,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nL && nL->isLeafNode() == false){
			float dM = nL->subnodeArray[1][1][0]->m_d111;
			this->edges[7]->computeABCFromSamples(no->m_d010,dM,no->m_d011);
			this->edges[7]->m_edgeType = octreeEdge::LONG;
		}

		else if(nBL && nBL->m_edgeLength > no->m_edgeLength){
			if(no->m_v[2] > nBL->m_v[2]){ //edge is a back-short edge
				octreeEdge::decompose(false,nBL->m_d100,no->m_d010,no->m_d011,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d010,no->m_d011,nBL->m_d101,_a,_b,_c);
				this->edges[7]->setABC(_a,_b,_c);
				this->edges[7]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nBL && nBL->isLeafNode() == false){
			float dM = nBL->subnodeArray[1][0][0]->m_d101;
			this->edges[7]->computeABCFromSamples(no->m_d010,dM,no->m_d011);
			this->edges[7]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[7]->setABC(0.0,no->m_d011 - no->m_d010, no->m_d010);
			this->edges[7]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{

		if(no->m_id ==10){
			int br = 0;
		}

		this->edges[7]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[2] > no->m_parent->m_v[2]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::D);
			dv1 = neighbor->m_d010;
			dv2 = no->m_d010;
			dv3 = no->m_d011;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[7]->setABC(aS,bS,cS);
		}
		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::U);
			dv1 = no->m_d010;
			dv2 = no->m_d011;
			dv3 = neighbor->m_d011;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[7]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 7

	//********************************************************************************
	// edge 8
	//********************************************************************************
	this->edges[8] = new octreeEdge(no->m_v + Point(0,0,no->m_edgeLength), no->m_v+ Point(no->m_edgeLength,0,no->m_edgeLength));
	this->edges[8]->m_nDir[0] = octreeNode::U;
	this->edges[8]->m_nDir[1] = octreeNode::F;
	if(!isEdgeExtra(this->edges[8],no,octreeNode::U,octreeNode::F)){
		
		octreeNode* nU = no->findNeighbor(octreeNode::U);
		octreeNode* nF = no->findNeighbor(octreeNode::F);
		octreeNode* nUF = NULL;
		if(nU) nUF = nU->findNeighbor(octreeNode::F);
		if(nU && nU->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nU->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nU->m_d000,no->m_d001,no->m_d101,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d101,nU->m_d100,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nU && nU->isLeafNode() == false){
			float dM = nU->subnodeArray[0][0][0]->m_d100;
			this->edges[8]->computeABCFromSamples(no->m_d001,dM,no->m_d101);
			this->edges[8]->m_edgeType = octreeEdge::LONG;
		}
		else if(nF && nF->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nF->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nF->m_d011,no->m_d001,no->m_d101,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d101,nF->m_d111,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nF && nF->isLeafNode() == false){
			float dM = nF->subnodeArray[0][1][1]->m_d111;
			this->edges[8]->computeABCFromSamples(no->m_d001,dM,no->m_d101);
			this->edges[8]->m_edgeType = octreeEdge::LONG;
		}

		else if(nUF && nUF->m_edgeLength > no->m_edgeLength){
			if(no->m_v[0] > nUF->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nUF->m_d010,no->m_d001,no->m_d101,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d101,nUF->m_d110,_a,_b,_c);
				this->edges[8]->setABC(_a,_b,_c);
				this->edges[8]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nUF && nUF->isLeafNode() == false){
			float dM = nUF->subnodeArray[0][1][0]->m_d110;
			this->edges[8]->computeABCFromSamples(no->m_d001,dM,no->m_d101);
			this->edges[8]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[8]->setABC(0.0,no->m_d101 - no->m_d001, no->m_d001);
			this->edges[8]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[8]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[0] > no->m_parent->m_v[0]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::L);
			dv1 = neighbor->m_d001;
			dv2 = no->m_d001;
			dv3 = no->m_d101;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[8]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::R);
			dv1 = no->m_d001;
			dv2 = no->m_d101;
			dv3 = neighbor->m_d101;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[8]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 8

	//********************************************************************************
	// edge 9
	//********************************************************************************
	this->edges[9] = new octreeEdge(no->m_v + Point(no->m_edgeLength,0,no->m_edgeLength), no->m_v+ Point(no->m_edgeLength,no->m_edgeLength,no->m_edgeLength));
	this->edges[9]->m_nDir[0] = octreeNode::U;
	this->edges[9]->m_nDir[1] = octreeNode::R;
	if(!isEdgeExtra(this->edges[9],no,octreeNode::U,octreeNode::R)){
		
		octreeNode* nU = no->findNeighbor(octreeNode::U);
		octreeNode* nR = no->findNeighbor(octreeNode::R);
		octreeNode* nUR = NULL;
		if(nU) nUR = nU->findNeighbor(octreeNode::R);
		if(nU && nU->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nU->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nU->m_d100,no->m_d101,no->m_d111,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d101,no->m_d111,nU->m_d110,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nU && nU->isLeafNode() == false){
			float dM = nU->subnodeArray[1][0][0]->m_d110;
			this->edges[9]->computeABCFromSamples(no->m_d101,dM,no->m_d111);
			this->edges[9]->m_edgeType = octreeEdge::LONG;
		}
		else if(nR && nR->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nR->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nR->m_d001,no->m_d101,no->m_d111,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d101,no->m_d111,nR->m_d011,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nR && nR->isLeafNode() == false){
			float dM = nR->subnodeArray[0][1][1]->m_d001;
			this->edges[9]->computeABCFromSamples(no->m_d101,dM,no->m_d111);
			this->edges[9]->m_edgeType = octreeEdge::LONG;
		}

		else if(nUR && nUR->m_edgeLength > no->m_edgeLength){
			if(no->m_v[1] > nUR->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nUR->m_d000,no->m_d101,no->m_d111,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d101,no->m_d111,nUR->m_d010,_a,_b,_c);
				this->edges[9]->setABC(_a,_b,_c);
				this->edges[9]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nUR && nUR->isLeafNode() == false){
			float dM = nUR->subnodeArray[0][0][0]->m_d010;
			if(no->m_id ==2) { 
				int b=0;
			}
			this->edges[9]->computeABCFromSamples(no->m_d101,dM,no->m_d111);
			this->edges[9]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[9]->setABC(0.0,no->m_d111 - no->m_d101, no->m_d101);
			this->edges[9]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[9]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[1] > no->m_parent->m_v[1]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::F);
			dv1 = neighbor->m_d101;
			dv2 = no->m_d101;
			dv3 = no->m_d111;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[9]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::B);
			dv1 = no->m_d101;
			dv2 = no->m_d111;
			dv3 = neighbor->m_d111;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[9]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 9

	
	//********************************************************************************
	// edge 10
	//********************************************************************************
	this->edges[10] = new octreeEdge(no->m_v + Point(0,no->m_edgeLength,no->m_edgeLength), no->m_v+ Point(no->m_edgeLength,no->m_edgeLength,no->m_edgeLength));
	this->edges[10]->m_nDir[0] = octreeNode::U;
	this->edges[10]->m_nDir[1] = octreeNode::B;
	if(!isEdgeExtra(this->edges[10],no,octreeNode::U,octreeNode::B)){
		
		octreeNode* nU = no->findNeighbor(octreeNode::U);
		octreeNode* nB = no->findNeighbor(octreeNode::B);
		octreeNode* nUB = NULL;
		if(nU) nUB = nU->findNeighbor(octreeNode::B);
		if(nU && nU->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nU->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nU->m_d010,no->m_d011,no->m_d111,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d011,no->m_d111,nU->m_d110,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nU && nU->isLeafNode() == false){
			float dM = nU->subnodeArray[1][1][0]->m_d010;
			this->edges[10]->computeABCFromSamples(no->m_d011,dM,no->m_d111);
			this->edges[10]->m_edgeType = octreeEdge::LONG;
		}
		else if(nB && nB->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[0] > nB->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nB->m_d001,no->m_d011,no->m_d111,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d011,no->m_d111,nB->m_d101,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nB && nB->isLeafNode() == false){
			float dM = nB->subnodeArray[0][0][1]->m_d101;
			this->edges[10]->computeABCFromSamples(no->m_d011,dM,no->m_d111);
			this->edges[10]->m_edgeType = octreeEdge::LONG;
		}

		else if(nUB && nUB->m_edgeLength > no->m_edgeLength){
			if(no->m_v[0] > nUB->m_v[0]){ //edge is a back-short edge
				octreeEdge::decompose(false,nUB->m_d000,no->m_d011,no->m_d111,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d011,no->m_d111,nUB->m_d100,_a,_b,_c);
				this->edges[10]->setABC(_a,_b,_c);
				this->edges[10]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nUB && nUB->isLeafNode() == false){
			float dM = nUB->subnodeArray[0][0][0]->m_d100;
			this->edges[10]->computeABCFromSamples(no->m_d011,dM,no->m_d111);
			this->edges[10]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[10]->setABC(0.0,no->m_d111 - no->m_d011, no->m_d011);
			this->edges[10]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[10]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[0] > no->m_parent->m_v[0]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::L);
			dv1 = neighbor->m_d011;
			dv2 = no->m_d011;
			dv3 = no->m_d111;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[10]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::R);
			dv1 = no->m_d011;
			dv2 = no->m_d111;
			dv3 = neighbor->m_d111;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[10]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 10

	//********************************************************************************
	// edge 11
	//********************************************************************************
	this->edges[11] = new octreeEdge(no->m_v + Point(0,0,no->m_edgeLength), no->m_v+ Point(0,no->m_edgeLength,no->m_edgeLength));
	this->edges[11]->m_nDir[0] = octreeNode::U;
	this->edges[11]->m_nDir[1] = octreeNode::L;
	if(!isEdgeExtra(this->edges[11],no,octreeNode::U,octreeNode::L)){
		
		octreeNode* nU = no->findNeighbor(octreeNode::U);
		octreeNode* nL = no->findNeighbor(octreeNode::L);
		octreeNode* nUL = NULL;
		if(nU) nUL = nU->findNeighbor(octreeNode::L);
		if(nU && nU->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nU->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nU->m_d000,no->m_d001,no->m_d011,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d011,nU->m_d010,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nU && nU->isLeafNode() == false){
			float dM = nU->subnodeArray[0][0][0]->m_d010;
			this->edges[11]->computeABCFromSamples(no->m_d001,dM,no->m_d011);
			this->edges[11]->m_edgeType = octreeEdge::LONG;
		}
		else if(nL && nL->m_edgeLength > no->m_edgeLength) {
			//edge is a "short" edge
			if(no->m_v[1] > nL->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nL->m_d101,no->m_d001,no->m_d011,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d011,nL->m_d111,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}
		}
		else if(nL && nL->isLeafNode() == false){
			float dM = nL->subnodeArray[1][1][1]->m_d101;
			this->edges[11]->computeABCFromSamples(no->m_d001,dM,no->m_d011);
			this->edges[11]->m_edgeType = octreeEdge::LONG;
		}

		else if(nUL && nUL->m_edgeLength > no->m_edgeLength){
			if(no->m_v[1] > nUL->m_v[1]){ //edge is a back-short edge
				octreeEdge::decompose(false,nUL->m_d100,no->m_d001,no->m_d011,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}
			else{
				octreeEdge::decompose(true,no->m_d001,no->m_d011,nUL->m_d110,_a,_b,_c);
				this->edges[11]->setABC(_a,_b,_c);
				this->edges[11]->m_edgeType = octreeEdge::SHORT;
			}		
		}
		else if(nUL && nUL->isLeafNode() == false){
			float dM = nUL->subnodeArray[1][1][0]->m_d100;
			this->edges[11]->computeABCFromSamples(no->m_d001,dM,no->m_d011);
			this->edges[11]->m_edgeType = octreeEdge::LONG;
		}
		
		//edge is regular
		else{
			this->edges[11]->setABC(0.0,no->m_d011 - no->m_d001, no->m_d001);
			this->edges[11]->m_edgeType = octreeEdge::REGULAR;
		}
	}
     //edge is extra 
	else{
		this->edges[11]->m_edgeType = octreeEdge::TSHORT;
		float dv1,dv2,dv3;
		float aS,bS,cS;
		if(no->m_v[1] > no->m_parent->m_v[1]){
			octreeNode* neighbor = no->findNeighbor(octreeNode::F);
			dv1 = neighbor->m_d001;
			dv2 = no->m_d001;
			dv3 = no->m_d011;
			octreeEdge::decompose(false,dv1,dv2,dv3,aS,bS,cS);
			this->edges[11]->setABC(aS,bS,cS);
		}

		else{
			octreeNode* neighbor = no->findNeighbor(octreeNode::B);
			dv1 = no->m_d001;
			dv2 = no->m_d011;
			dv3 = neighbor->m_d011;
			octreeEdge::decompose(true,dv1,dv2,dv3,aS,bS,cS);
			this->edges[11]->setABC(aS,bS,cS);
		}
	}
	//end of process for edge 11
}

void CoonsSdfInterp::processOctreeFaces(octreeNode *no){
	//face 0 with edge 0,1,2 and 3, R0(u,v)
	float d00,d10,d01,d11;
	int ind0,ind1,ind2,ind3;
	/////////////////////////////////////////////////////////////////////////////
	// bottom face (face 0) with edge 0,1,2 and 3
	/////////////////////////////////////////////////////////////////////////////
	ind0=3; ind1= 1; ind2 = 0; ind3 =2;
	d00 = no->m_d000; d10 = no->m_d100; d01 = no->m_d010; d11 = no->m_d110;
	octreeNode* tempNeighbor = no->findNeighbor(octreeNode::D);
	
	//case 1: the number of samples on this face is 4~8...
	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for R0(u,v);
		m_R0[0]=0.0;
		m_R0[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_R0[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_R0[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() + 
			-d00+d01+d10-d11;
		m_R0[4]= this->edges[ind2]->getA();
		m_R0[5]= this->edges[ind0]->getA();
		m_R0[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_R0[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_R0[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the bottom face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::D);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[0][0][1]; 
		d0= no->m_d000; d2 = no->m_d100; d6 = no->m_d010; d8 = no->m_d110;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d111;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_R0[0]= face->getC(0); m_R0[1]= face->getC(1); m_R0[2]= face->getC(2); 
		m_R0[3]= face->getC(3); m_R0[4]= face->getC(4); m_R0[5]= face->getC(5); 
		m_R0[6]= face->getC(6); m_R0[7]= face->getC(7); m_R0[8]= face->getC(8);
		delete face;
	}
	else {
		//the bottom face is a small face that takes 1/4 of the neighbor
		if(no->m_id ==10){
			float br = 0;
		}
		octreeNode* node = no->findNeighbor(octreeNode::D);
		octreeNode* n00,*n10,*n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[0][0][0]; 
		n10 = no->m_parent->subnodeArray[1][0][0];
		n01 = no->m_parent->subnodeArray[0][1][0]; 
		d0= node->m_d001; d2 = node->m_d101; d6 = node->m_d011; d8 = node->m_d111;
		d1= n00->m_d100;  d3 = n00->m_d010;  d4 = n00->m_d110;  d5 = n10->m_d110;
		d7= n01->m_d110;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::LFD:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::RFD:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::LBD:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::RBD:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_R0[0]= Sface.getC(0); m_R0[1]= Sface.getC(1); m_R0[2]= Sface.getC(2); 
		m_R0[3]= Sface.getC(3); m_R0[4]= Sface.getC(4); m_R0[5]= Sface.getC(5); 
		m_R0[6]= Sface.getC(6); m_R0[7]= Sface.getC(7); m_R0[8]= Sface.getC(8);
		delete Bface;
	}
	/////////////////////////////////////////////////////////////////////////////
	// top face (face 1) with edge 0,1,2 and 3
	/////////////////////////////////////////////////////////////////////////////
	ind0=11; ind1= 9; ind2 = 8; ind3 =10;
	d00 = no->m_d001; d10 = no->m_d101; d01 = no->m_d011; d11 = no->m_d111;
	tempNeighbor = no->findNeighbor(octreeNode::U);
	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for R1(u,v);
		m_R1[0]=0.0;
		m_R1[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_R1[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_R1[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() + 
			-d00+d01+d10-d11;
		m_R1[4]= this->edges[ind2]->getA();
		m_R1[5]= this->edges[ind0]->getA();
		m_R1[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_R1[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_R1[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the top face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::U);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[0][0][0]; 
		d0= no->m_d001; d2 = no->m_d101; d6 = no->m_d011; d8 = no->m_d111;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d110;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_R1[0]= face->getC(0); m_R1[1]= face->getC(1); m_R1[2]= face->getC(2); 
		m_R1[3]= face->getC(3); m_R1[4]= face->getC(4); m_R1[5]= face->getC(5); 
		m_R1[6]= face->getC(6); m_R1[7]= face->getC(7); m_R1[8]= face->getC(8);
		delete face;
	}
	
	else {
		//the top face is a small face that takes 1/4 of the neighbor
		octreeNode* node = no->findNeighbor(octreeNode::U);
		octreeNode* n00,*n10,*n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[0][0][1]; 
		n10 = no->m_parent->subnodeArray[1][0][1];
		n01 = no->m_parent->subnodeArray[0][1][1]; 
		d0= node->m_d000; d2 = node->m_d100; d6 = node->m_d010; d8 = node->m_d110;
		d1= n00->m_d101;  d3 = n00->m_d011;  d4 = n00->m_d111;  d5 = n10->m_d111;
		d7= n01->m_d111;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::LFU:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::RFU:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::LBU:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::RBU:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_R1[0]= Sface.getC(0); m_R1[1]= Sface.getC(1); m_R1[2]= Sface.getC(2); 
		m_R1[3]= Sface.getC(3); m_R1[4]= Sface.getC(4); m_R1[5]= Sface.getC(5); 
		m_R1[6]= Sface.getC(6); m_R1[7]= Sface.getC(7); m_R1[8]= Sface.getC(8);
		delete Bface;
	}

	/////////////////////////////////////////////////////////////////////////////
	// front face (face 2) with edge 0,5,8 and 4
	/////////////////////////////////////////////////////////////////////////////

	ind0=4; ind1= 5; ind2 = 0; ind3 =8;
	d00 = no->m_d000; d10 = no->m_d100; d01 = no->m_d001; d11 = no->m_d101;
	
	tempNeighbor = no->findNeighbor(octreeNode::F);
	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for Q0(u,w);
		m_Q0[0]=0.0;
		m_Q0[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_Q0[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_Q0[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() //+ 
			-d00+d01+d10-d11;
		m_Q0[4]= this->edges[ind2]->getA();
		m_Q0[5]= this->edges[ind0]->getA();
		m_Q0[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_Q0[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_Q0[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the bottom face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::F);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[0][1][0]; 
		d0= no->m_d000; d2 = no->m_d100; d6 = no->m_d001; d8 = no->m_d101;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d111;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_Q0[0]= face->getC(0); m_Q0[1]= face->getC(1); m_Q0[2]= face->getC(2); 
		m_Q0[3]= face->getC(3); m_Q0[4]= face->getC(4); m_Q0[5]= face->getC(5); 
		m_Q0[6]= face->getC(6); m_Q0[7]= face->getC(7); m_Q0[8]= face->getC(8);
		delete face;
	}
	
	else {
		//the front face is a small face that takes 1/4 of the neighbor
		octreeNode* node = no->findNeighbor(octreeNode::F);
		octreeNode* n00,*n10,*n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[0][0][0]; 
		n10 = no->m_parent->subnodeArray[1][0][0];
		n01 = no->m_parent->subnodeArray[0][0][1]; 
		d0= node->m_d010; d2 = node->m_d110; d6 = node->m_d011; d8 = node->m_d111;
		d1= n00->m_d100;  d3 = n00->m_d001;  d4 = n00->m_d101;  d5 = n10->m_d101;
		d7= n01->m_d101;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::LFD:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::RFD:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::LFU:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::RFU:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_Q0[0]= Sface.getC(0); m_Q0[1]= Sface.getC(1); m_Q0[2]= Sface.getC(2); 
		m_Q0[3]= Sface.getC(3); m_Q0[4]= Sface.getC(4); m_Q0[5]= Sface.getC(5); 
		m_Q0[6]= Sface.getC(6); m_Q0[7]= Sface.getC(7); m_Q0[8]= Sface.getC(8);
		delete Bface;
	}

	/////////////////////////////////////////////////////////////////////////////
	// back face (face 3) with edge 2,6,10 and 7
	/////////////////////////////////////////////////////////////////////////////

	ind0=7; ind1= 6; ind2 = 2; ind3 =10;
	d00 = no->m_d010; d10 = no->m_d110; d01 = no->m_d011; d11 = no->m_d111;
	tempNeighbor = no->findNeighbor(octreeNode::B);
	if(no->m_id ==10){
		int br = 0;
	}

	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for Q1(u,w);
		m_Q1[0]=0.0;
		m_Q1[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_Q1[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_Q1[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() + 
			-d00+d01+d10-d11;
		m_Q1[4]= this->edges[ind2]->getA();
		m_Q1[5]= this->edges[ind0]->getA();
		m_Q1[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_Q1[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_Q1[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the bottom face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::B);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[0][0][0]; 
		d0= no->m_d010; d2 = no->m_d110; d6 = no->m_d011; d8 = no->m_d111;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d101;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_Q1[0]= face->getC(0); m_Q1[1]= face->getC(1); m_Q1[2]= face->getC(2); 
		m_Q1[3]= face->getC(3); m_Q1[4]= face->getC(4); m_Q1[5]= face->getC(5); 
		m_Q1[6]= face->getC(6); m_Q1[7]= face->getC(7); m_Q1[8]= face->getC(8);
		delete face;
	}
	
	else {
		//the bottom face is a small face that takes 1/4 of the neighbor
		octreeNode* node = no->findNeighbor(octreeNode::B);
		octreeNode* n00,*n10,*n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[0][1][0]; 
		n10 = no->m_parent->subnodeArray[1][1][0];
		n01 = no->m_parent->subnodeArray[0][1][1]; 
		d0= node->m_d000; d2 = node->m_d100; d6 = node->m_d001; d8 = node->m_d101;
		d1= n00->m_d110;  d3 = n00->m_d011;  d4 = n00->m_d111;  d5 = n10->m_d111;
		d7= n01->m_d111;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::LBD:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::RBD:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::LBU:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::RBU:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_Q1[0]= Sface.getC(0); m_Q1[1]= Sface.getC(1); m_Q1[2]= Sface.getC(2); 
		m_Q1[3]= Sface.getC(3); m_Q1[4]= Sface.getC(4); m_Q1[5]= Sface.getC(5); 
		m_Q1[6]= Sface.getC(6); m_Q1[7]= Sface.getC(7); m_Q1[8]= Sface.getC(8);
		delete Bface;
	}

	
	/////////////////////////////////////////////////////////////////////////////
	// left face (face 3) with edge 3,7,11 and 4
	/////////////////////////////////////////////////////////////////////////////
	ind0=4; ind1= 7; ind2 = 3; ind3 =11;
	d00 = no->m_d000; d10 = no->m_d010; d01 = no->m_d001; d11 = no->m_d011;
	tempNeighbor = no->findNeighbor(octreeNode::L);
	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for P0(v,w);
		m_P0[0]=0.0;
		m_P0[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_P0[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_P0[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() + 
			-d00+d01+d10-d11;
		m_P0[4]= this->edges[ind2]->getA();
		m_P0[5]= this->edges[ind0]->getA();
		m_P0[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_P0[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_P0[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the back face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::L);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[1][0][0]; 
		d0= no->m_d000; d2 = no->m_d010; d6 = no->m_d001; d8 = no->m_d011;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d111;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_P0[0]= face->getC(0); m_P0[1]= face->getC(1); m_P0[2]= face->getC(2); 
		m_P0[3]= face->getC(3); m_P0[4]= face->getC(4); m_P0[5]= face->getC(5); 
		m_P0[6]= face->getC(6); m_P0[7]= face->getC(7); m_P0[8]= face->getC(8);
		delete face;
	}
	
	else {
		//the left face is a small face that takes 1/4 of the neighbor
		if(no->m_id == 10){
			int br = 0;
		}
		octreeNode* node = no->findNeighbor(octreeNode::L);
		octreeNode* n00,*n10,*n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[0][0][0]; 
		n10 = no->m_parent->subnodeArray[0][1][0];
		n01 = no->m_parent->subnodeArray[0][0][1]; 
		d0= node->m_d100; d2 = node->m_d110; d6 = node->m_d101; d8 = node->m_d111;
		d1= n00->m_d010;  d3 = n00->m_d001;  d4 = n00->m_d011;  d5 = n10->m_d011;
		d7= n01->m_d011;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::LFD:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::LBD:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::LFU:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::LBU:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_P0[0]= Sface.getC(0); m_P0[1]= Sface.getC(1); m_P0[2]= Sface.getC(2); 
		m_P0[3]= Sface.getC(3); m_P0[4]= Sface.getC(4); m_P0[5]= Sface.getC(5); 
		m_P0[6]= Sface.getC(6); m_P0[7]= Sface.getC(7); m_P0[8]= Sface.getC(8);
		delete Bface;
	}
	
	/////////////////////////////////////////////////////////////////////////////
	// right face (face 3) with edge 1,6,9 and 5
	/////////////////////////////////////////////////////////////////////////////
	ind0=5; ind1= 6; ind2 = 1; ind3 =9;
	d00 = no->m_d100; d10 = no->m_d110; d01 = no->m_d101; d11 = no->m_d111;
	tempNeighbor = no->findNeighbor(octreeNode::R);
	if(tempNeighbor == NULL || (tempNeighbor->isLeafNode() && (tempNeighbor->m_edgeLength == no->m_edgeLength)  ))
	{
		//setting coefficient for P1(v,w);
		m_P1[0]=0.0;
		m_P1[1]= -this->edges[ind2]->getA() + this->edges[ind3]->getA();
		m_P1[2]= -this->edges[ind0]->getA() + this->edges[ind1]->getA();
		m_P1[3]= -this->edges[ind0]->getB() + this->edges[ind1]->getB() - this->edges[ind2]->getB() + this->edges[ind3]->getB() + 
			-d00+d01+d10-d11;
		m_P1[4]= this->edges[ind2]->getA();
		m_P1[5]= this->edges[ind0]->getA();
		m_P1[6]= this->edges[ind2]->getB() - this->edges[ind0]->getC() + this->edges[ind1]->getC() + d00 - d10;
		m_P1[7]= this->edges[ind0]->getB() - this->edges[ind2]->getC() + this->edges[ind3]->getC() + d00 - d01;
		m_P1[8]= this->edges[ind0]->getC() + this->edges[ind2]->getC() - d00;
	}
	else if (this->edges[ind0]->m_edgeType == octreeEdge::LONG && this->edges[ind1]->m_edgeType == octreeEdge::LONG &&
	         this->edges[ind2]->m_edgeType == octreeEdge::LONG && this->edges[ind3]->m_edgeType == octreeEdge::LONG) {
		//the bottom face is a big face with 9 samples.
		octreeNode* node = no->findNeighbor(octreeNode::R);
		octreeNode* n001;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n001 = node->subnodeArray[0][0][0]; 
		d0= no->m_d100; d2 = no->m_d110; d6 = no->m_d101; d8 = no->m_d111;
		//d1= n001->m_d101; d3 = n001->m_d011; d4 = n001->m_d111; d5 = n101->m_d111;
		//d7= n011->m_d111;
			
		d1 = this->edges[ind2]->evaluateEdge(0.5);
		d3 = this->edges[ind0]->evaluateEdge(0.5);
		d5 = this->edges[ind1]->evaluateEdge(0.5);
		d7 = this->edges[ind3]->evaluateEdge(0.5);
		d4 = n001->m_d011;
		
		octreeFace* face = new octreeFace();
		face->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		m_P1[0]= face->getC(0); m_P1[1]= face->getC(1); m_P1[2]= face->getC(2); 
		m_P1[3]= face->getC(3); m_P1[4]= face->getC(4); m_P1[5]= face->getC(5); 
		m_P1[6]= face->getC(6); m_P1[7]= face->getC(7); m_P1[8]= face->getC(8);
		delete face;
	}
	
	else {
		//the right face is a small face that takes 1/4 of the neighbor
		octreeNode* node = no->findNeighbor(octreeNode::R);
		octreeNode* n00,*n10, *n01;
		float d0,d1,d2,d3,d4,d5,d6,d7,d8;
		n00 = no->m_parent->subnodeArray[1][0][0]; 
		n10 = no->m_parent->subnodeArray[1][1][0];
		n01 = no->m_parent->subnodeArray[1][0][1]; 
		d0= node->m_d000; d2 = node->m_d010; d6 = node->m_d001; d8 = node->m_d011;
		d1= n00->m_d110;  d3 = n00->m_d101;  d4 = n00->m_d111;  d5 = n10->m_d111;
		d7= n01->m_d111;
		octreeFace* Bface = new octreeFace();
		octreeFace Sface;
		Bface->computeABCFromSamples(d0,d1,d2,d3,d4,d5,d6,d7,d8);
		switch(no->m_octant){
		case octreeNode::RFD:
				Bface->decompose(0,(*Bface),Sface);
				break;
		case octreeNode::RBD:
				Bface->decompose(1,(*Bface),Sface);
				break;
		case octreeNode::RFU:
				Bface->decompose(2,(*Bface),Sface);
				break;
		case octreeNode::RBU:
				Bface->decompose(3,(*Bface),Sface);
				break;
		}
		m_P1[0]= Sface.getC(0); m_P1[1]= Sface.getC(1); m_P1[2]= Sface.getC(2); 
		m_P1[3]= Sface.getC(3); m_P1[4]= Sface.getC(4); m_P1[5]= Sface.getC(5); 
		m_P1[6]= Sface.getC(6); m_P1[7]= Sface.getC(7); m_P1[8]= Sface.getC(8);
		delete Bface;
	}
}

float evalInterpolationSurface(float* d,float u,float v){
	float eval = d[0]*u*u*v*v + d[1]*u*u*v + d[2]*u*v*v + d[3]*u*v + d[4]*u*u + d[5]*v*v + d[6]*u + d[7]*v + d[8]; 
	return eval;
}


///////////////////////////////////////////////////////////////////
// auxiliary evaluation functions
float CoonsSdfInterp::evalR0(float u, float v){
	return evalInterpolationSurface(m_R0,u,v);
}

float CoonsSdfInterp::evalR1(float u, float v){
	return evalInterpolationSurface(m_R1,u,v);
}

float CoonsSdfInterp::evalP0(float v, float w){
	return evalInterpolationSurface(m_P0,v,w);
}

float CoonsSdfInterp::evalP1(float v, float w){
	return evalInterpolationSurface(m_P1,v,w);
}

float CoonsSdfInterp::evalQ0(float u, float w){
	return evalInterpolationSurface(m_Q0,u,w);
}

float CoonsSdfInterp::evalQ1(float u, float w){
	return evalInterpolationSurface(m_Q1,u,w);
}