public enum OCT_EDGE implements OCT_ENUM {

	/**
	 * BACK DOWN, Direction X,Y,Z: 0, -1, -1
	 */
	BD(0, -1, -1, OCT_VERTEX.LBD, OCT_VERTEX.RBD), // 0
	/**
	 * RIGHT DOWN, Direction X,Y,Z: 1, 0, -1
	 */
	RD(1, 0, -1, OCT_VERTEX.RBD, OCT_VERTEX.RFD), // 1
	/**
	 * FRONT DOWN, Direction X,Y,Z: 0, 1, -1
	 */
	FD(0, 1, -1, OCT_VERTEX.LFD, OCT_VERTEX.RFD), // 2
	/**
	 * LEFT DOWN, Direction X,Y,Z: -1, 0, -1
	 */
	LD(-1, 0, -1, OCT_VERTEX.LBD, OCT_VERTEX.LFD), // 3
	/**
	 * BACK UP, Direction X,Y,Z: 0, -1, -1
	 */
	BU(0, -1, 1, OCT_VERTEX.LBU, OCT_VERTEX.RBU), // 4
	/**
	 * RIGHT UP, Direction X,Y,Z: 1, 0, 1
	 */
	RU(1, 0, 1, OCT_VERTEX.RBU, OCT_VERTEX.RFU), // 5
	/**
	 * FRONT UP, Direction X,Y,Z: 0, 1, 1
	 */
	FU(0, 1, 1, OCT_VERTEX.LFU, OCT_VERTEX.RFU), // 6
	/**
	 * LEFT UP, Direction X,Y,Z: -1, 0, 1
	 */
	LU(-1, 0, 1, OCT_VERTEX.LBU, OCT_VERTEX.LFU), // 7
	/**
	 * LEFT BACK, Direction X,Y,Z: -1, -1, 0
	 */
	LB(-1, -1, 0, OCT_VERTEX.LBD, OCT_VERTEX.LBU), // 8
	/**
	 * RIGHT BACK, Direction X,Y,Z: 1, -1, 0
	 */
	RB(1, -1, 0, OCT_VERTEX.RBD, OCT_VERTEX.RBU), // 9
	/**
	 * RIGHT FRONT, Direction X,Y,Z: 1, 1, 0
	 */
	RF(1, 1, 0, OCT_VERTEX.RFD, OCT_VERTEX.RFU), // 10
	/**
	 * LEFT FRONT Direction X,Y,Z: -1, 1, 0
	 */
	LF(-1, 1, 0, OCT_VERTEX.LFD, OCT_VERTEX.LFU); // 11

	public final int r; // vector from the node center to the edge center
	public final int s;
	public final int t;

	final OCT_VERTEX v0;
	final OCT_VERTEX v1;

	OCT_EDGE(int _r, int _s, int _t, OCT_VERTEX _v0, OCT_VERTEX _v1) {
		r = _r;
		s = _s;
		t = _t;
		v0 = _v0;
		v1 = _v1;
	}

	/**
	 * Returns the direction from the node center to the edge center in the X
	 * axis.
	 */
	public final int getR() {
		return r;
	}

	/**
	 * Returns the direction from the node center to the edge center in the Y
	 * axis.
	 */
	public final int getS() {
		return s;
	}

	/**
	 * Returns the direction from the node center to the edge center in the Z
	 * axis.
	 */
	public final int getT() {
		return t;
	}

	/**
	 * Returns the position of the edge in the list of edges. Same result as the
	 * java.lang.Enum method ordinal().
	 */
	public final int getOrdinal() {
		return this.ordinal();
	}

	/**
	 * Returns the list of all edges as an ArrayList.
	 */
	public final static ArrayList<OCT_EDGE> getAll() {
		ArrayList<OCT_EDGE> t = new ArrayList<OCT_EDGE>();
		for (OCT_EDGE e : OCT_EDGE.values()) {
			t.add(e);
		}
		return t;
	}

	/**
	 * Returns the edge at the 'i' position in the list of edges.
	 */
	public final static OCT_EDGE get(int _i) {
		if (_i > -1 && _i < OCT_EDGE.values().length) {
			return OCT_EDGE.values()[_i];
		} else
			return null;
	}

	/**
	 * Returns the edge at the specified direction, if it exists. Otherwise returns null.
	 */
	public final static OCT_EDGE get(int _dirR, int _dirS, int _dirT) {
		for (OCT_EDGE e : OCT_EDGE.values()) {
			if (e.r == _dirR && e.s == _dirS && e.t == _dirT) {
				return e;
			}
		}
		return null;
	}

	/**
	 * Returns the first vertex of an edge.
	 */
	public final OCT_VERTEX getV0() {
		return v0;
	}

	/**
	 * Returns the second vertex of an edge.
	 */
	public final OCT_VERTEX getV1() {
		return v1;
	}

	// *************************************************************************************
	// GET OTHER PARTS OF THE NODE
	// *************************************************************************************

	/**
	 * Returns a list with the two vertices of the edge.
	 */
	public final ArrayList<OCT_VERTEX> getVertices() {
		ArrayList<OCT_VERTEX> list = new ArrayList<OCT_VERTEX>();
		list.add(v0);
		list.add(v1);
		return list;
	}

	/**
	 * Returns a list with the two faces of the edge.
	 */
	public final ArrayList<OCT_FACE> getFaces() {
		ArrayList<OCT_FACE> list = new ArrayList<OCT_FACE>();
		switch (this) {
		case BD:
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.D);
			break;
		case RD:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.D);
			break;
		case FD:
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.D);
			break;
		case LD:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.D);
			break;
		case BU:
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.U);
			break;
		case RU:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.U);
			break;
		case FU:
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.U);
			break;
		case LU:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.U);
			break;
		case LB:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.B);
			break;
		case RB:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.B);
			break;
		case RF:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.F);
			break;
		case LF:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.F);
			break;
		}
		return list;
	}

	/**
	 * Returns itself, a list with only one element.
	 */
	public final ArrayList<OCT_EDGE> getEdges() {
		ArrayList<OCT_EDGE> list = new ArrayList<OCT_EDGE>();
		list.add(this);
		return list;
	}

	/**
	 * Returns a list with the two octants touching the edge.
	 */
	public final ArrayList<OCT_OCTANT> getOctants() {
		ArrayList<OCT_OCTANT> list = new ArrayList<OCT_OCTANT>();
		switch (this) {
		case BD:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.RBD);
			break;
		case RD:
			list.add(OCT_OCTANT.RBD);
			list.add(OCT_OCTANT.RFD);
			break;
		case FD:
			list.add(OCT_OCTANT.LFD);
			list.add(OCT_OCTANT.RFD);
			break;
		case LD:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.LFD);
			break;
		case BU:
			list.add(OCT_OCTANT.LBU);
			list.add(OCT_OCTANT.RBU);
			break;
		case RU:
			list.add(OCT_OCTANT.RBU);
			list.add(OCT_OCTANT.RFU);
			break;
		case FU:
			list.add(OCT_OCTANT.LFU);
			list.add(OCT_OCTANT.RFU);
			break;
		case LU:
			list.add(OCT_OCTANT.LBU);
			list.add(OCT_OCTANT.LFU);
			break;
		case LB:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.LBU);
			break;
		case RB:
			list.add(OCT_OCTANT.RBD);
			list.add(OCT_OCTANT.RBU);
			break;
		case RF:
			list.add(OCT_OCTANT.RFD);
			list.add(OCT_OCTANT.RFU);
			break;
		case LF:
			list.add(OCT_OCTANT.LFD);
			list.add(OCT_OCTANT.LFU);
			break;
		}
		return list;
	}
} public interface OCT_ENUM {

	public int getR(); // direction in the X axis

	public int getS(); // direction in the Y axis

	public int getT(); // direction in the Z axis

	public int getOrdinal();

	public ArrayList<OCT_VERTEX> getVertices();

	public ArrayList<OCT_EDGE> getEdges();

	public ArrayList<OCT_FACE> getFaces();

	public ArrayList<OCT_OCTANT> getOctants();
	
	// public ArrayList<? extends OCT_ENUM> getAll();

} public enum OCT_FACE implements OCT_ENUM {

	/**
	 * LEFT, Direction X,Y,Z: -1, 0, 0
	 */
	L(-1, 0, 0, OCT_VERTEX.LFD, OCT_VERTEX.LFU, OCT_VERTEX.LBU, OCT_VERTEX.LBD,
			OCT_EDGE.LF, OCT_EDGE.LU, OCT_EDGE.LB, OCT_EDGE.LD),
	/**
	 * RIGHT, Direction X,Y,Z: 1, 0, 0
	 */
	R(1, 0, 0, OCT_VERTEX.RBD, OCT_VERTEX.RBU, OCT_VERTEX.RFU, OCT_VERTEX.RFD,
			OCT_EDGE.RB, OCT_EDGE.RU, OCT_EDGE.RF, OCT_EDGE.RD),
	/**
	 * BACK, Direction X,Y,Z: 0, -1, 0
	 */
	B(0, -1, 0, OCT_VERTEX.LBD, OCT_VERTEX.LBU, OCT_VERTEX.RBU, OCT_VERTEX.RBD,
			OCT_EDGE.LB, OCT_EDGE.BU, OCT_EDGE.RB, OCT_EDGE.BD),
	/**
	 * FRONT, Direction X,Y,Z: 0, 1, 0
	 */
	F(0, 1, 0, OCT_VERTEX.RFD, OCT_VERTEX.RFU, OCT_VERTEX.LFU, OCT_VERTEX.LFD,
			OCT_EDGE.RF, OCT_EDGE.FU, OCT_EDGE.LF, OCT_EDGE.FD),
	/**
	 * DOWN, Direction X,Y,Z: 0, 0, -1
	 */
	D(0, 0, -1, OCT_VERTEX.LBD, OCT_VERTEX.RBD, OCT_VERTEX.RFD, OCT_VERTEX.LFD,
			OCT_EDGE.BD, OCT_EDGE.RD, OCT_EDGE.FD, OCT_EDGE.LD),
	/**
	 * UP, Direction X,Y,Z: 0, 0, 1
	 */
	U(0, 0, 1, OCT_VERTEX.LBU, OCT_VERTEX.LFU, OCT_VERTEX.RFU, OCT_VERTEX.RBU,
			OCT_EDGE.LU, OCT_EDGE.FU, OCT_EDGE.RU, OCT_EDGE.BU);

	public final int r; // vector from the node center to the EDGE center
	public final int s;
	public final int t;

	final OCT_VERTEX v0; // the FACE contains these 4 vertices
	final OCT_VERTEX v1;
	final OCT_VERTEX v2;
	final OCT_VERTEX v3;

	final OCT_EDGE e0; // the FACE contains these 4 edges
	final OCT_EDGE e1;
	final OCT_EDGE e2;
	final OCT_EDGE e3;

	OCT_FACE(int _r, int _s, int _t, OCT_VERTEX _v0, OCT_VERTEX _v1,
			OCT_VERTEX _v2, OCT_VERTEX _v3, OCT_EDGE _e0, OCT_EDGE _e1,
			OCT_EDGE _e2, OCT_EDGE _e3) {
		r = _r;
		s = _s;
		t = _t;
		v0 = _v0;
		v1 = _v1;
		v2 = _v2;
		v3 = _v3;
		e0 = _e0;
		e1 = _e1;
		e2 = _e2;
		e3 = _e3;
	}

	/**
	 * Returns the direction from the node center to the face center in the X
	 * axis.
	 */
	public final int getR() {
		return r;
	}

	/**
	 * Returns the direction from the node center to the face center in the Y
	 * axis.
	 */
	public final int getS() {
		return s;
	}

	/**
	 * Returns the direction from the node center to the face center in the Z
	 * axis.
	 */
	public final int getT() {
		return t;
	}

	/**
	 * Returns the position of the face in the list of faces. Same result as the
	 * java.lang.Enum method ordinal().
	 */
	public final int getOrdinal() {
		return this.ordinal();
	}

	/**
	 * Returns the list of all faces as an ArrayList.
	 */
	public final static ArrayList<OCT_FACE> getAll() {
		ArrayList<OCT_FACE> t = new ArrayList<OCT_FACE>();
		for (OCT_FACE e : OCT_FACE.values()) {
			t.add(e);
		}
		return t;
	}

	/**
	 * Returns the face at the 'i' position in the list of faces.
	 */
	public final static OCT_FACE get(int _i) {
		if (_i > -1 && _i < OCT_FACE.values().length) {
			return OCT_FACE.values()[_i];
		} else
			return null;
	}

	/**
	 * Returns the face at the specified direction, if it exists. Otherwise returns null.
	 */
	public final static OCT_FACE get(int _dirR, int _dirS, int _dirT) {
		for (OCT_FACE f : OCT_FACE.values()) {
			if (f.r == _dirR && f.s == _dirS && f.t == _dirT) {
				return f;
			}
		}
		return null;
	}

	/**
	 * Returns the first vertex of a face.
	 */
	public final OCT_VERTEX getV0() {
		return v0;
	}

	/**
	 * Returns the second vertex of a face.
	 */
	public final OCT_VERTEX getV1() {
		return v1;
	}

	/**
	 * Returns the third vertex of a face.
	 */
	public final OCT_VERTEX getV2() {
		return v2;
	}

	/**
	 * Returns the fourth vertex of a face.
	 */
	public final OCT_VERTEX getV3() {
		return v3;
	}

	/**
	 * Returns the first edge of a face.
	 */
	public final OCT_EDGE getE0() {
		return e0;
	}

	/**
	 * Returns the second edge of a face.
	 */
	public final OCT_EDGE getE1() {
		return e1;
	}

	/**
	 * Returns the third edge of a face.
	 */
	public final OCT_EDGE getE2() {
		return e2;
	}

	/**
	 * Returns the fourth edge of a face.
	 */
	public final OCT_EDGE getE3() {
		return e3;
	}

	// *************************************************************************************
	// GET OTHER PARTS OF THE NODE
	// *************************************************************************************

	/**
	 * Returns a list with the four vertices of the face.
	 */
	public final ArrayList<OCT_VERTEX> getVertices() {
		ArrayList<OCT_VERTEX> list = new ArrayList<OCT_VERTEX>();
		list.add(v0);
		list.add(v1);
		list.add(v2);
		list.add(v3);
		return list;
	}

	/**
	 * Returns a list with the four edges of the face.
	 */
	public final ArrayList<OCT_EDGE> getEdges() {
		ArrayList<OCT_EDGE> list = new ArrayList<OCT_EDGE>();
		list.add(e0);
		list.add(e1);
		list.add(e2);
		list.add(e3);
		return list;
	}

	/**
	 * Returns itself, a list with only one element.
	 */
	public final ArrayList<OCT_FACE> getFaces() {
		ArrayList<OCT_FACE> list = new ArrayList<OCT_FACE>();
		list.add(this);
		return list;
	}

	/**
	 * Returns a list with the four octants touching the face.
	 */
	public final ArrayList<OCT_OCTANT> getOctants() {
		ArrayList<OCT_OCTANT> list = new ArrayList<OCT_OCTANT>();
		switch (this) {
		case L:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.LFD);
			list.add(OCT_OCTANT.LBU);
			list.add(OCT_OCTANT.LFU);
			break;
		case R:
			list.add(OCT_OCTANT.RBD);
			list.add(OCT_OCTANT.RFD);
			list.add(OCT_OCTANT.RBU);
			list.add(OCT_OCTANT.RFU);
			break;
		case B:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.RBD);
			list.add(OCT_OCTANT.LBU);
			list.add(OCT_OCTANT.RBU);
			break;
		case F:
			list.add(OCT_OCTANT.RFD);
			list.add(OCT_OCTANT.LFD);
			list.add(OCT_OCTANT.RFU);
			list.add(OCT_OCTANT.LFU);
			break;
		case D:
			list.add(OCT_OCTANT.LBD);
			list.add(OCT_OCTANT.RBD);
			list.add(OCT_OCTANT.RFD);
			list.add(OCT_OCTANT.LFD);
			break;
		case U:
			list.add(OCT_OCTANT.LBU);
			list.add(OCT_OCTANT.RBU);
			list.add(OCT_OCTANT.RFU);
			list.add(OCT_OCTANT.LFU);
			break;
		}
		return list;
	}
} public enum OCT_VERTEX implements OCT_ENUM {

	/**
	 * LBD, Direction X,Y,Z: -1, -1, -1
	 */
	LBD(-1, -1, -1), // 0
	/**
	 * RBD, Direction X,Y,Z: 1, -1, -1
	 */
	RBD(1, -1, -1), // 1
	/**
	 * RFD, Direction X,Y,Z: 1, 1, -1
	 */
	RFD(1, 1, -1), // 2
	/**
	 * LFD, Direction X,Y,Z: -1, 1, -1
	 */
	LFD(-1, 1, -1), // 3
	/**
	 * LBU, Direction X,Y,Z: -1, -1, 1
	 */
	LBU(-1, -1, 1), // 4
	/**
	 * RBU, Direction X,Y,Z: 1, -1, 1
	 */
	RBU(1, -1, 1), // 5
	/**
	 * RFU, Direction X,Y,Z: 1, 1, 1
	 */
	RFU(1, 1, 1), // 6
	/**
	 * LFU, Direction X,Y,Z: -1, 1, 1
	 */
	LFU(-1, 1, 1); // 7

	public final int r; // vector from the node center to the vertex
	public final int s;
	public final int t;

	OCT_VERTEX(int _r, int _s, int _t) {
		r = _r;
		s = _s;
		t = _t;
	}

	/**
	 * Returns the direction from the node center to the vertex in the Z axis.
	 */
	public final int getR() {
		return r;
	}

	/**
	 * Returns the direction from the node center to the vertex in the Z axis.
	 */
	public final int getS() {
		return s;
	}

	/**
	 * Returns the direction from the node center to the vertex in the Z axis.
	 */
	public final int getT() {
		return t;
	}

	/**
	 * Returns the position of the vertex in the list of vertices. Same result
	 * as the java.lang.Enum method ordinal().
	 */
	public final int getOrdinal() {
		return this.ordinal();
	}

	/**
	 * Returns the list of all vertices as an ArrayList.
	 */
	public final static ArrayList<OCT_VERTEX> getAll() {
		ArrayList<OCT_VERTEX> t = new ArrayList<OCT_VERTEX>();
		for (OCT_VERTEX e : OCT_VERTEX.values()) {
			t.add(e);
		}
		return t;
	}
	
	/**
	 * Returns the vertex at the 'i' position in the list of vertices.
	 */
	public final static OCT_VERTEX get(int _i) {
		if (_i > -1 && _i < OCT_VERTEX.values().length) {
			return OCT_VERTEX.values()[_i];
		} else
			return null;
	}

	/**
	 * Returns the vertex at the specified direction, if it exists. Otherwise returns null.
	 */
	public final static OCT_VERTEX get(int _dirR, int _dirS, int _dirT) {
		for (OCT_VERTEX v : OCT_VERTEX.values()) {
			if (v.r == _dirR && v.s == _dirS && v.t == _dirT) {
				return v;
			}
		}
		;
		return null;
	}

	// *************************************************************************************
	// GET OTHER PARTS OF THE NODE
	// *************************************************************************************

	/**
	 * Returns itself, a list with only one element.
	 */
	public final ArrayList<OCT_VERTEX> getVertices() {
		ArrayList<OCT_VERTEX> list = new ArrayList<OCT_VERTEX>();
		list.add(this);
		return list;
	}

	/**
	 * Returns a list with the three faces touching the vertex.
	 */
	public final ArrayList<OCT_FACE> getFaces() {
		ArrayList<OCT_FACE> list = new ArrayList<OCT_FACE>();
		switch (this) {
		case LBD:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.D);
			break;
		case RBD:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.D);
			break;
		case RFD:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.D);
			break;
		case LFD:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.D);
			break;
		case LBU:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.U);
			break;
		case RBU:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.B);
			list.add(OCT_FACE.U);
			break;
		case RFU:
			list.add(OCT_FACE.R);
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.U);
			break;
		case LFU:
			list.add(OCT_FACE.L);
			list.add(OCT_FACE.F);
			list.add(OCT_FACE.U);
			break;
		}
		return list;
	}

	/**
	 * Returns a list with the three edges touching the vertex.
	 */
	public final ArrayList<OCT_EDGE> getEdges() {
		ArrayList<OCT_EDGE> list = new ArrayList<OCT_EDGE>();
		switch (this) {
		case LBD:
			list.add(OCT_EDGE.BD);
			list.add(OCT_EDGE.LD);
			list.add(OCT_EDGE.LB);
			break;
		case RBD:
			list.add(OCT_EDGE.BD);
			list.add(OCT_EDGE.RD);
			list.add(OCT_EDGE.RB);
			break;
		case RFD:
			list.add(OCT_EDGE.FD);
			list.add(OCT_EDGE.RD);
			list.add(OCT_EDGE.RF);
			break;
		case LFD:
			list.add(OCT_EDGE.FD);
			list.add(OCT_EDGE.LD);
			list.add(OCT_EDGE.LF);
			break;
		case LBU:
			list.add(OCT_EDGE.BU);
			list.add(OCT_EDGE.LU);
			list.add(OCT_EDGE.LB);
			break;
		case RBU:
			list.add(OCT_EDGE.BU);
			list.add(OCT_EDGE.RU);
			list.add(OCT_EDGE.RB);
			break;
		case RFU:
			list.add(OCT_EDGE.FU);
			list.add(OCT_EDGE.RU);
			list.add(OCT_EDGE.RF);
			break;
		case LFU:
			list.add(OCT_EDGE.FU);
			list.add(OCT_EDGE.LU);
			list.add(OCT_EDGE.LF);
			break;
		}
		return list;
	}

	/**
	 * Returns a list with the octant corresponding to the vertex.
	 */
	public final ArrayList<OCT_OCTANT> getOctants() {
		ArrayList<OCT_OCTANT> list = new ArrayList<OCT_OCTANT>();
		switch (this) {
		case LBD:
			list.add(OCT_OCTANT.LBD);
			break;
		case RBD:
			list.add(OCT_OCTANT.RBD);
			break;
		case RFD:
			list.add(OCT_OCTANT.RFD);
			break;
		case LFD:
			list.add(OCT_OCTANT.LFD);
			break;
		case LBU:
			list.add(OCT_OCTANT.LBU);
			break;
		case RBU:
			list.add(OCT_OCTANT.RBU);
			break;
		case RFU:
			list.add(OCT_OCTANT.RFU);
			break;
		case LFU:
			list.add(OCT_OCTANT.LFU);
			break;
		}
		return list;
	}
} public interface OctFunction {

	/**
	 * To be implemented by the user, custom function.
	 */
	float compute(float x, float y, float z);

} public class Matrix implements Cloneable, java.io.Serializable {

/* ------------------------
   Class variables
 * ------------------------ */

   /** Array for internal storage of elements.
   @serial internal array storage.
   */
   private double[][] A;

   /** Row and column dimensions.
   @serial row dimension.
   @serial column dimension.
   */
   private int m, n;

/* ------------------------
   Constructors
 * ------------------------ */

   /** Construct an m-by-n matrix of zeros. 
   @param m    Number of rows.
   @param n    Number of colums.
   */

   public Matrix (int m, int n) {
      this.m = m;
      this.n = n;
      A = new double[m][n];
   }

   /** Construct an m-by-n constant matrix.
   @param m    Number of rows.
   @param n    Number of colums.
   @param s    Fill the matrix with this scalar value.
   */

   public Matrix (int m, int n, double s) {
      this.m = m;
      this.n = n;
      A = new double[m][n];
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = s;
         }
      }
   }

   /** Construct a matrix from a 2-D array.
   @param A    Two-dimensional array of doubles.
   @exception  IllegalArgumentException All rows must have the same length
   @see        #constructWithCopy
   */

   public Matrix (double[][] A) {
      m = A.length;
      n = A[0].length;
      for (int i = 0; i < m; i++) {
         if (A[i].length != n) {
            throw new IllegalArgumentException("All rows must have the same length.");
         }
      }
      this.A = A;
   }

   /** Construct a matrix quickly without checking arguments.
   @param A    Two-dimensional array of doubles.
   @param m    Number of rows.
   @param n    Number of colums.
   */

   public Matrix (double[][] A, int m, int n) {
      this.A = A;
      this.m = m;
      this.n = n;
   }

   /** Construct a matrix from a one-dimensional packed array
   @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
   @param m    Number of rows.
   @exception  IllegalArgumentException Array length must be a multiple of m.
   */

   public Matrix (double vals[], int m) {
      this.m = m;
      n = (m != 0 ? vals.length/m : 0);
      if (m*n != vals.length) {
         throw new IllegalArgumentException("Array length must be a multiple of m.");
      }
      A = new double[m][n];
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = vals[i+j*m];
         }
      }
   }

/* ------------------------
   Public Methods
 * ------------------------ */

   /** Construct a matrix from a copy of a 2-D array.
   @param A    Two-dimensional array of doubles.
   @exception  IllegalArgumentException All rows must have the same length
   */

   public static Matrix constructWithCopy(double[][] A) {
      int m = A.length;
      int n = A[0].length;
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         if (A[i].length != n) {
            throw new IllegalArgumentException
               ("All rows must have the same length.");
         }
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j];
         }
      }
      return X;
   }

   /** Make a deep copy of a matrix
   */

   public Matrix copy () {
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j];
         }
      }
      return X;
   }

   /** Clone the Matrix object.
   */

   public Object clone () {
      return this.copy();
   }

   /** Access the internal two-dimensional array.
   @return     Pointer to the two-dimensional array of matrix elements.
   */

   public double[][] getArray () {
      return A;
   }

   /** Copy the internal two-dimensional array.
   @return     Two-dimensional array copy of matrix elements.
   */

   public double[][] getArrayCopy () {
      double[][] C = new double[m][n];
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j];
         }
      }
      return C;
   }

   /** Make a one-dimensional column packed copy of the internal array.
   @return     Matrix elements packed in a one-dimensional array by columns.
   */

   public double[] getColumnPackedCopy () {
      double[] vals = new double[m*n];
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            vals[i+j*m] = A[i][j];
         }
      }
      return vals;
   }

   /** Make a one-dimensional row packed copy of the internal array.
   @return     Matrix elements packed in a one-dimensional array by rows.
   */

   public double[] getRowPackedCopy () {
      double[] vals = new double[m*n];
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            vals[i*n+j] = A[i][j];
         }
      }
      return vals;
   }

   /** Get row dimension.
   @return     m, the number of rows.
   */

   public int getRowDimension () {
      return m;
   }

   /** Get column dimension.
   @return     n, the number of columns.
   */

   public int getColumnDimension () {
      return n;
   }

   /** Get a single element.
   @param i    Row index.
   @param j    Column index.
   @return     A(i,j)
   @exception  ArrayIndexOutOfBoundsException
   */

   public double get (int i, int j) {
      return A[i][j];
   }

   /** Get a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param j0   Initial column index
   @param j1   Final column index
   @return     A(i0:i1,j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int i0, int i1, int j0, int j1) {
      Matrix X = new Matrix(i1-i0+1,j1-j0+1);
      double[][] B = X.getArray();
      try {
         for (int i = i0; i <= i1; i++) {
            for (int j = j0; j <= j1; j++) {
               B[i-i0][j-j0] = A[i][j];
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Get a submatrix.
   @param r    Array of row indices.
   @param c    Array of column indices.
   @return     A(r(:),c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int[] r, int[] c) {
      Matrix X = new Matrix(r.length,c.length);
      double[][] B = X.getArray();
      try {
         for (int i = 0; i < r.length; i++) {
            for (int j = 0; j < c.length; j++) {
               B[i][j] = A[r[i]][c[j]];
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Get a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param c    Array of column indices.
   @return     A(i0:i1,c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int i0, int i1, int[] c) {
      Matrix X = new Matrix(i1-i0+1,c.length);
      double[][] B = X.getArray();
      try {
         for (int i = i0; i <= i1; i++) {
            for (int j = 0; j < c.length; j++) {
               B[i-i0][j] = A[i][c[j]];
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Get a submatrix.
   @param r    Array of row indices.
   @param j0   Initial column index
   @param j1   Final column index
   @return     A(r(:),j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int[] r, int j0, int j1) {
      Matrix X = new Matrix(r.length,j1-j0+1);
      double[][] B = X.getArray();
      try {
         for (int i = 0; i < r.length; i++) {
            for (int j = j0; j <= j1; j++) {
               B[i][j-j0] = A[r[i]][j];
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Set a single element.
   @param i    Row index.
   @param j    Column index.
   @param s    A(i,j).
   @exception  ArrayIndexOutOfBoundsException
   */

   public void set (int i, int j, double s) {
      A[i][j] = s;
   }

   /** Set a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param j0   Initial column index
   @param j1   Final column index
   @param X    A(i0:i1,j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int i0, int i1, int j0, int j1, Matrix X) {
      try {
         for (int i = i0; i <= i1; i++) {
            for (int j = j0; j <= j1; j++) {
               A[i][j] = X.get(i-i0,j-j0);
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param r    Array of row indices.
   @param c    Array of column indices.
   @param X    A(r(:),c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int[] r, int[] c, Matrix X) {
      try {
         for (int i = 0; i < r.length; i++) {
            for (int j = 0; j < c.length; j++) {
               A[r[i]][c[j]] = X.get(i,j);
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param r    Array of row indices.
   @param j0   Initial column index
   @param j1   Final column index
   @param X    A(r(:),j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int[] r, int j0, int j1, Matrix X) {
      try {
         for (int i = 0; i < r.length; i++) {
            for (int j = j0; j <= j1; j++) {
               A[r[i]][j] = X.get(i,j-j0);
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param c    Array of column indices.
   @param X    A(i0:i1,c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int i0, int i1, int[] c, Matrix X) {
      try {
         for (int i = i0; i <= i1; i++) {
            for (int j = 0; j < c.length; j++) {
               A[i][c[j]] = X.get(i-i0,j);
            }
         }
      } catch(ArrayIndexOutOfBoundsException e) {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Matrix transpose.
   @return    A'
   */

   public Matrix transpose () {
      Matrix X = new Matrix(n,m);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[j][i] = A[i][j];
         }
      }
      return X;
   }

   /** One norm
   @return    maximum column sum.
   */

   public double norm1 () {
      double f = 0;
      for (int j = 0; j < n; j++) {
         double s = 0;
         for (int i = 0; i < m; i++) {
            s += Math.abs(A[i][j]);
         }
         f = Math.max(f,s);
      }
      return f;
   }

   /** Two norm
   @return    maximum singular value.
   */

   public double norm2 () {
      return (new SingularValueDecomposition(this).norm2());
   }

   /** Infinity norm
   @return    maximum row sum.
   */

   public double normInf () {
      double f = 0;
      for (int i = 0; i < m; i++) {
         double s = 0;
         for (int j = 0; j < n; j++) {
            s += Math.abs(A[i][j]);
         }
         f = Math.max(f,s);
      }
      return f;
   }

   /** Frobenius norm
   @return    sqrt of sum of squares of all elements.
   */

   public double normF () {
      double f = 0;
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            f = Maths.hypot(f,A[i][j]);
         }
      }
      return f;
   }

   /**  Unary minus
   @return    -A
   */

   public Matrix uminus () {
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = -A[i][j];
         }
      }
      return X;
   }

   /** C = A + B
   @param B    another matrix
   @return     A + B
   */

   public Matrix plus (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] + B.A[i][j];
         }
      }
      return X;
   }

   /** A = A + B
   @param B    another matrix
   @return     A + B
   */

   public Matrix plusEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] + B.A[i][j];
         }
      }
      return this;
   }

   /** C = A - B
   @param B    another matrix
   @return     A - B
   */

   public Matrix minus (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] - B.A[i][j];
         }
      }
      return X;
   }

   /** A = A - B
   @param B    another matrix
   @return     A - B
   */

   public Matrix minusEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] - B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element multiplication, C = A.*B
   @param B    another matrix
   @return     A.*B
   */

   public Matrix arrayTimes (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] * B.A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element multiplication in place, A = A.*B
   @param B    another matrix
   @return     A.*B
   */

   public Matrix arrayTimesEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] * B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element right division, C = A./B
   @param B    another matrix
   @return     A./B
   */

   public Matrix arrayRightDivide (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] / B.A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element right division in place, A = A./B
   @param B    another matrix
   @return     A./B
   */

   public Matrix arrayRightDivideEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] / B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element left division, C = A.\B
   @param B    another matrix
   @return     A.\B
   */

   public Matrix arrayLeftDivide (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = B.A[i][j] / A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element left division in place, A = A.\B
   @param B    another matrix
   @return     A.\B
   */

   public Matrix arrayLeftDivideEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = B.A[i][j] / A[i][j];
         }
      }
      return this;
   }

   /** Multiply a matrix by a scalar, C = s*A
   @param s    scalar
   @return     s*A
   */

   public Matrix times (double s) {
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = s*A[i][j];
         }
      }
      return X;
   }

   /** Multiply a matrix by a scalar in place, A = s*A
   @param s    scalar
   @return     replace A by s*A
   */

   public Matrix timesEquals (double s) {
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = s*A[i][j];
         }
      }
      return this;
   }

   /** Linear algebraic matrix multiplication, A * B
   @param B    another matrix
   @return     Matrix product, A * B
   @exception  IllegalArgumentException Matrix inner dimensions must agree.
   */

   public Matrix times (Matrix B) {
      if (B.m != n) {
         throw new IllegalArgumentException("Matrix inner dimensions must agree.");
      }
      Matrix X = new Matrix(m,B.n);
      double[][] C = X.getArray();
      double[] Bcolj = new double[n];
      for (int j = 0; j < B.n; j++) {
         for (int k = 0; k < n; k++) {
            Bcolj[k] = B.A[k][j];
         }
         for (int i = 0; i < m; i++) {
            double[] Arowi = A[i];
            double s = 0;
            for (int k = 0; k < n; k++) {
               s += Arowi[k]*Bcolj[k];
            }
            C[i][j] = s;
         }
      }
      return X;
   }

   /** LU Decomposition
   @return     LUDecomposition
   @see LUDecomposition
   */

   public LUDecomposition lu () {
      return new LUDecomposition(this);
   }

   /** QR Decomposition
   @return     QRDecomposition
   @see QRDecomposition
   */

   public QRDecomposition qr () {
      return new QRDecomposition(this);
   }

   /** Cholesky Decomposition
   @return     CholeskyDecomposition
   @see CholeskyDecomposition
   */

   public CholeskyDecomposition chol () {
      return new CholeskyDecomposition(this);
   }

   /** Singular Value Decomposition
   @return     SingularValueDecomposition
   @see SingularValueDecomposition
   */

   public SingularValueDecomposition svd () {
      return new SingularValueDecomposition(this);
   }

   /** Eigenvalue Decomposition
   @return     EigenvalueDecomposition
   @see EigenvalueDecomposition
   */

   public EigenvalueDecomposition eig () {
      return new EigenvalueDecomposition(this);
   }

   /** Solve A*X = B
   @param B    right hand side
   @return     solution if A is square, least squares solution otherwise
   */

   public Matrix solve (Matrix B) {
      return (m == n ? (new LUDecomposition(this)).solve(B) :
                       (new QRDecomposition(this)).solve(B));
   }

   /** Solve X*A = B, which is also A'*X' = B'
   @param B    right hand side
   @return     solution if A is square, least squares solution otherwise.
   */

   public Matrix solveTranspose (Matrix B) {
      return transpose().solve(B.transpose());
   }

   /** Matrix inverse or pseudoinverse
   @return     inverse(A) if A is square, pseudoinverse otherwise.
   */

   public Matrix inverse () {
      return solve(identity(m,m));
   }

   /** Matrix determinant
   @return     determinant
   */

   public double det () {
      return new LUDecomposition(this).det();
   }

   /** Matrix rank
   @return     effective numerical rank, obtained from SVD.
   */

   public int rank () {
      return new SingularValueDecomposition(this).rank();
   }

   /** Matrix condition (2 norm)
   @return     ratio of largest to smallest singular value.
   */

   public double cond () {
      return new SingularValueDecomposition(this).cond();
   }

   /** Matrix trace.
   @return     sum of the diagonal elements.
   */

   public double trace () {
      double t = 0;
      for (int i = 0; i < Math.min(m,n); i++) {
         t += A[i][i];
      }
      return t;
   }

   /** Generate matrix with random elements
   @param m    Number of rows.
   @param n    Number of colums.
   @return     An m-by-n matrix with uniformly distributed random elements.
   */

   public static Matrix random (int m, int n) {
      Matrix A = new Matrix(m,n);
      double[][] X = A.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            X[i][j] = Math.random();
         }
      }
      return A;
   }

   /** Generate identity matrix
   @param m    Number of rows.
   @param n    Number of colums.
   @return     An m-by-n matrix with ones on the diagonal and zeros elsewhere.
   */

   public static Matrix identity (int m, int n) {
      Matrix A = new Matrix(m,n);
      double[][] X = A.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            X[i][j] = (i == j ? 1.0 : 0.0);
         }
      }
      return A;
   }


   /** Print the matrix to stdout.   Line the elements up in columns
     * with a Fortran-like 'Fw.d' style format.
   @param w    Column width.
   @param d    Number of digits after the decimal.
   */

   public void print (int w, int d) {
      print(new PrintWriter(System.out,true),w,d); }

   /** Print the matrix to the output stream.   Line the elements up in
     * columns with a Fortran-like 'Fw.d' style format.
   @param output Output stream.
   @param w      Column width.
   @param d      Number of digits after the decimal.
   */

   public void print (PrintWriter output, int w, int d) {
      DecimalFormat format = new DecimalFormat();
      format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
      format.setMinimumIntegerDigits(1);
      format.setMaximumFractionDigits(d);
      format.setMinimumFractionDigits(d);
      format.setGroupingUsed(false);
      print(output,format,w+2);
   }

   /** Print the matrix to stdout.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
   @param format A  Formatting object for individual elements.
   @param width     Field width for each column.
   @see java.text.DecimalFormat#setDecimalFormatSymbols
   */

   public void print (NumberFormat format, int width) {
      print(new PrintWriter(System.out,true),format,width); }

   // DecimalFormat is a little disappointing coming from Fortran or C's printf.
   // Since it doesn't pad on the left, the elements will come out different
   // widths.  Consequently, we'll pass the desired column width in as an
   // argument and do the extra padding ourselves.

   /** Print the matrix to the output stream.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
   @param output the output stream.
   @param format A formatting object to format the matrix elements 
   @param width  Column width.
   @see java.text.DecimalFormat#setDecimalFormatSymbols
   */

   public void print (PrintWriter output, NumberFormat format, int width) {
      output.println();  // start on new line.
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            String s = format.format(A[i][j]); // format the number
            int padding = Math.max(1,width-s.length()); // At _least_ 1 space
            for (int k = 0; k < padding; k++)
               output.print(' ');
            output.print(s);
         }
         output.println();
      }
      output.println();   // end with blank line.
   }

   /** Read a matrix from a stream.  The format is the same the print method,
     * so printed matrices can be read back in (provided they were printed using
     * US Locale).  Elements are separated by
     * whitespace, all the elements for each row appear on a single line,
     * the last row is followed by a blank line.
   @param input the input stream.
   */

   public static Matrix read (BufferedReader input) throws java.io.IOException {
      StreamTokenizer tokenizer= new StreamTokenizer(input);

      // Although StreamTokenizer will parse numbers, it doesn't recognize
      // scientific notation (E or D); however, Double.valueOf does.
      // The strategy here is to disable StreamTokenizer's number parsing.
      // We'll only get whitespace delimited words, EOL's and EOF's.
      // These words should all be numbers, for Double.valueOf to parse.

      tokenizer.resetSyntax();
      tokenizer.wordChars(0,255);
      tokenizer.whitespaceChars(0, ' ');
      tokenizer.eolIsSignificant(true);
      java.util.Vector<Double> vD = new java.util.Vector<Double>();

      // Ignore initial empty lines
      while (tokenizer.nextToken() == StreamTokenizer.TT_EOL);
      if (tokenizer.ttype == StreamTokenizer.TT_EOF)
	throw new java.io.IOException("Unexpected EOF on matrix read.");
      do {
         vD.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
      } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);

      int n = vD.size();  // Now we've got the number of columns!
      double row[] = new double[n];
      for (int j=0; j<n; j++)  // extract the elements of the 1st row.
         row[j]=vD.elementAt(j).doubleValue();
      java.util.Vector<double[]> v = new java.util.Vector<double[]>();
      v.addElement(row);  // Start storing rows instead of columns.
      while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
         // While non-empty lines
         v.addElement(row = new double[n]);
         int j = 0;
         do {
            if (j >= n) throw new java.io.IOException
               ("Row " + v.size() + " is too long.");
            row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
         } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
         if (j < n) throw new java.io.IOException
            ("Row " + v.size() + " is too short.");
      }
      int m = v.size();  // Now we've got the number of rows.
      double[][] A = new double[m][];
      v.copyInto(A);  // copy the rows out of the vector
      return new Matrix(A);
   }


/* ------------------------
   Private Methods
 * ------------------------ */

   /** Check if size(A) == size(B) **/

   private void checkMatrixDimensions (Matrix B) {
      if (B.m != m || B.n != n) {
         throw new IllegalArgumentException("Matrix dimensions must agree.");
      }
   }

  private static final long serialVersionUID = 1;
} public class SingularValueDecomposition implements java.io.Serializable {

/* ------------------------
   Class variables
 * ------------------------ */

   /** Arrays for internal storage of U and V.
   @serial internal storage of U.
   @serial internal storage of V.
   */
   private double[][] U, V;

   /** Array for internal storage of singular values.
   @serial internal storage of singular values.
   */
   private double[] s;

   /** Row and column dimensions.
   @serial row dimension.
   @serial column dimension.
   */
   private int m, n;

/* ------------------------
   Constructor
 * ------------------------ */

   /** Construct the singular value decomposition
       Structure to access U, S and V.
   @param Arg    Rectangular matrix
   */

   public SingularValueDecomposition (Matrix Arg) {

      // Derived from LINPACK code.
      // Initialize.
      double[][] A = Arg.getArrayCopy();
      m = Arg.getRowDimension();
      n = Arg.getColumnDimension();

      /* Apparently the failing cases are only a proper subset of (m<n), 
	 so let's not throw error.  Correct fix to come later?
      if (m<n) {
	  throw new IllegalArgumentException("Jama SVD only works for m >= n"); }
      */
      int nu = Math.min(m,n);
      s = new double [Math.min(m+1,n)];
      U = new double [m][nu];
      V = new double [n][n];
      double[] e = new double [n];
      double[] work = new double [m];
      boolean wantu = true;
      boolean wantv = true;

      // Reduce A to bidiagonal form, storing the diagonal elements
      // in s and the super-diagonal elements in e.

      int nct = Math.min(m-1,n);
      int nrt = Math.max(0,Math.min(n-2,m));
      for (int k = 0; k < Math.max(nct,nrt); k++) {
         if (k < nct) {

            // Compute the transformation for the k-th column and
            // place the k-th diagonal in s[k].
            // Compute 2-norm of k-th column without under/overflow.
            s[k] = 0;
            for (int i = k; i < m; i++) {
               s[k] = Maths.hypot(s[k],A[i][k]);
            }
            if (s[k] != 0.0) {
               if (A[k][k] < 0.0) {
                  s[k] = -s[k];
               }
               for (int i = k; i < m; i++) {
                  A[i][k] /= s[k];
               }
               A[k][k] += 1.0;
            }
            s[k] = -s[k];
         }
         for (int j = k+1; j < n; j++) {
            if ((k < nct) & (s[k] != 0.0))  {

            // Apply the transformation.

               double t = 0;
               for (int i = k; i < m; i++) {
                  t += A[i][k]*A[i][j];
               }
               t = -t/A[k][k];
               for (int i = k; i < m; i++) {
                  A[i][j] += t*A[i][k];
               }
            }

            // Place the k-th row of A into e for the
            // subsequent calculation of the row transformation.

            e[j] = A[k][j];
         }
         if (wantu & (k < nct)) {

            // Place the transformation in U for subsequent back
            // multiplication.

            for (int i = k; i < m; i++) {
               U[i][k] = A[i][k];
            }
         }
         if (k < nrt) {

            // Compute the k-th row transformation and place the
            // k-th super-diagonal in e[k].
            // Compute 2-norm without under/overflow.
            e[k] = 0;
            for (int i = k+1; i < n; i++) {
               e[k] = Maths.hypot(e[k],e[i]);
            }
            if (e[k] != 0.0) {
               if (e[k+1] < 0.0) {
                  e[k] = -e[k];
               }
               for (int i = k+1; i < n; i++) {
                  e[i] /= e[k];
               }
               e[k+1] += 1.0;
            }
            e[k] = -e[k];
            if ((k+1 < m) & (e[k] != 0.0)) {

            // Apply the transformation.

               for (int i = k+1; i < m; i++) {
                  work[i] = 0.0;
               }
               for (int j = k+1; j < n; j++) {
                  for (int i = k+1; i < m; i++) {
                     work[i] += e[j]*A[i][j];
                  }
               }
               for (int j = k+1; j < n; j++) {
                  double t = -e[j]/e[k+1];
                  for (int i = k+1; i < m; i++) {
                     A[i][j] += t*work[i];
                  }
               }
            }
            if (wantv) {

            // Place the transformation in V for subsequent
            // back multiplication.

               for (int i = k+1; i < n; i++) {
                  V[i][k] = e[i];
               }
            }
         }
      }

      // Set up the final bidiagonal matrix or order p.

      int p = Math.min(n,m+1);
      if (nct < n) {
         s[nct] = A[nct][nct];
      }
      if (m < p) {
         s[p-1] = 0.0;
      }
      if (nrt+1 < p) {
         e[nrt] = A[nrt][p-1];
      }
      e[p-1] = 0.0;

      // If required, generate U.

      if (wantu) {
         for (int j = nct; j < nu; j++) {
            for (int i = 0; i < m; i++) {
               U[i][j] = 0.0;
            }
            U[j][j] = 1.0;
         }
         for (int k = nct-1; k >= 0; k--) {
            if (s[k] != 0.0) {
               for (int j = k+1; j < nu; j++) {
                  double t = 0;
                  for (int i = k; i < m; i++) {
                     t += U[i][k]*U[i][j];
                  }
                  t = -t/U[k][k];
                  for (int i = k; i < m; i++) {
                     U[i][j] += t*U[i][k];
                  }
               }
               for (int i = k; i < m; i++ ) {
                  U[i][k] = -U[i][k];
               }
               U[k][k] = 1.0 + U[k][k];
               for (int i = 0; i < k-1; i++) {
                  U[i][k] = 0.0;
               }
            } else {
               for (int i = 0; i < m; i++) {
                  U[i][k] = 0.0;
               }
               U[k][k] = 1.0;
            }
         }
      }

      // If required, generate V.

      if (wantv) {
         for (int k = n-1; k >= 0; k--) {
            if ((k < nrt) & (e[k] != 0.0)) {
               for (int j = k+1; j < nu; j++) {
                  double t = 0;
                  for (int i = k+1; i < n; i++) {
                     t += V[i][k]*V[i][j];
                  }
                  t = -t/V[k+1][k];
                  for (int i = k+1; i < n; i++) {
                     V[i][j] += t*V[i][k];
                  }
               }
            }
            for (int i = 0; i < n; i++) {
               V[i][k] = 0.0;
            }
            V[k][k] = 1.0;
         }
      }

      // Main iteration loop for the singular values.

      int pp = p-1;
      int iter = 0;
      double eps = Math.pow(2.0,-52.0);
      double tiny = Math.pow(2.0,-966.0);
      while (p > 0) {
         int k,kase;

         // Here is where a test for too many iterations would go.

         // This section of the program inspects for
         // negligible elements in the s and e arrays.  On
         // completion the variables kase and k are set as follows.

         // kase = 1     if s(p) and e[k-1] are negligible and k<p
         // kase = 2     if s(k) is negligible and k<p
         // kase = 3     if e[k-1] is negligible, k<p, and
         //              s(k), ..., s(p) are not negligible (qr step).
         // kase = 4     if e(p-1) is negligible (convergence).

         for (k = p-2; k >= -1; k--) {
            if (k == -1) {
               break;
            }
            if (Math.abs(e[k]) <=
                  tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
               e[k] = 0.0;
               break;
            }
         }
         if (k == p-2) {
            kase = 4;
         } else {
            int ks;
            for (ks = p-1; ks >= k; ks--) {
               if (ks == k) {
                  break;
               }
               double t = (ks != p ? Math.abs(e[ks]) : 0.) + 
                          (ks != k+1 ? Math.abs(e[ks-1]) : 0.);
               if (Math.abs(s[ks]) <= tiny + eps*t)  {
                  s[ks] = 0.0;
                  break;
               }
            }
            if (ks == k) {
               kase = 3;
            } else if (ks == p-1) {
               kase = 1;
            } else {
               kase = 2;
               k = ks;
            }
         }
         k++;

         // Perform the task indicated by kase.

         switch (kase) {

            // Deflate negligible s(p).

            case 1: {
               double f = e[p-2];
               e[p-2] = 0.0;
               for (int j = p-2; j >= k; j--) {
                  double t = Maths.hypot(s[j],f);
                  double cs = s[j]/t;
                  double sn = f/t;
                  s[j] = t;
                  if (j != k) {
                     f = -sn*e[j-1];
                     e[j-1] = cs*e[j-1];
                  }
                  if (wantv) {
                     for (int i = 0; i < n; i++) {
                        t = cs*V[i][j] + sn*V[i][p-1];
                        V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
                        V[i][j] = t;
                     }
                  }
               }
            }
            break;

            // Split at negligible s(k).

            case 2: {
               double f = e[k-1];
               e[k-1] = 0.0;
               for (int j = k; j < p; j++) {
                  double t = Maths.hypot(s[j],f);
                  double cs = s[j]/t;
                  double sn = f/t;
                  s[j] = t;
                  f = -sn*e[j];
                  e[j] = cs*e[j];
                  if (wantu) {
                     for (int i = 0; i < m; i++) {
                        t = cs*U[i][j] + sn*U[i][k-1];
                        U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
                        U[i][j] = t;
                     }
                  }
               }
            }
            break;

            // Perform one qr step.

            case 3: {

               // Calculate the shift.
   
               double scale = Math.max(Math.max(Math.max(Math.max(
                       Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])), 
                       Math.abs(s[k])),Math.abs(e[k]));
               double sp = s[p-1]/scale;
               double spm1 = s[p-2]/scale;
               double epm1 = e[p-2]/scale;
               double sk = s[k]/scale;
               double ek = e[k]/scale;
               double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
               double c = (sp*epm1)*(sp*epm1);
               double shift = 0.0;
               if ((b != 0.0) | (c != 0.0)) {
                  shift = Math.sqrt(b*b + c);
                  if (b < 0.0) {
                     shift = -shift;
                  }
                  shift = c/(b + shift);
               }
               double f = (sk + sp)*(sk - sp) + shift;
               double g = sk*ek;
   
               // Chase zeros.
   
               for (int j = k; j < p-1; j++) {
                  double t = Maths.hypot(f,g);
                  double cs = f/t;
                  double sn = g/t;
                  if (j != k) {
                     e[j-1] = t;
                  }
                  f = cs*s[j] + sn*e[j];
                  e[j] = cs*e[j] - sn*s[j];
                  g = sn*s[j+1];
                  s[j+1] = cs*s[j+1];
                  if (wantv) {
                     for (int i = 0; i < n; i++) {
                        t = cs*V[i][j] + sn*V[i][j+1];
                        V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
                        V[i][j] = t;
                     }
                  }
                  t = Maths.hypot(f,g);
                  cs = f/t;
                  sn = g/t;
                  s[j] = t;
                  f = cs*e[j] + sn*s[j+1];
                  s[j+1] = -sn*e[j] + cs*s[j+1];
                  g = sn*e[j+1];
                  e[j+1] = cs*e[j+1];
                  if (wantu && (j < m-1)) {
                     for (int i = 0; i < m; i++) {
                        t = cs*U[i][j] + sn*U[i][j+1];
                        U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
                        U[i][j] = t;
                     }
                  }
               }
               e[p-2] = f;
               iter = iter + 1;
            }
            break;

            // Convergence.

            case 4: {

               // Make the singular values positive.
   
               if (s[k] <= 0.0) {
                  s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
                  if (wantv) {
                     for (int i = 0; i <= pp; i++) {
                        V[i][k] = -V[i][k];
                     }
                  }
               }
   
               // Order the singular values.
   
               while (k < pp) {
                  if (s[k] >= s[k+1]) {
                     break;
                  }
                  double t = s[k];
                  s[k] = s[k+1];
                  s[k+1] = t;
                  if (wantv && (k < n-1)) {
                     for (int i = 0; i < n; i++) {
                        t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
                     }
                  }
                  if (wantu && (k < m-1)) {
                     for (int i = 0; i < m; i++) {
                        t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
                     }
                  }
                  k++;
               }
               iter = 0;
               p--;
            }
            break;
         }
      }
   }

/* ------------------------
   Public Methods
 * ------------------------ */

   /** Return the left singular vectors
   @return     U
   */

   public Matrix getU () {
      return new Matrix(U,m,Math.min(m+1,n));
   }

   /** Return the right singular vectors
   @return     V
   */

   public Matrix getV () {
      return new Matrix(V,n,n);
   }

   /** Return the one-dimensional array of singular values
   @return     diagonal of S.
   */

   public double[] getSingularValues () {
      return s;
   }

   /** Return the diagonal matrix of singular values
   @return     S
   */

   public Matrix getS () {
      Matrix X = new Matrix(n,n);
      double[][] S = X.getArray();
      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            S[i][j] = 0.0;
         }
         S[i][i] = this.s[i];
      }
      return X;
   }

   /** Two norm
   @return     max(S)
   */

   public double norm2 () {
      return s[0];
   }

   /** Two norm condition number
   @return     max(S)/min(S)
   */

   public double cond () {
      return s[0]/s[Math.min(m,n)-1];
   }

   /** Effective numerical matrix rank
   @return     Number of nonnegligible singular values.
   */

   public int rank () {
      double eps = Math.pow(2.0,-52.0);
      double tol = Math.max(m,n)*s[0]*eps;
      int r = 0;
      for (int i = 0; i < s.length; i++) {
         if (s[i] > tol) {
            r++;
         }
      }
      return r;
   }
  private static final long serialVersionUID = 1;
} public class OctNode implements java.io.Serializable {

	// *************************************************************************************
	// VARIABLES
	// *************************************************************************************

	protected short codeR;
	protected short codeS;
	protected short codeT;
	protected byte level;
	protected boolean isSelected = false;

	// *************************************************************************************
	// CONSTRUCTORS
	// *************************************************************************************

	/**
	 * The position of the node in the tree is defined by its code (R, S, T
	 * coordinates) and its level. This is not the position is space (X, Y, Z
	 * coordinates) but the position in the data structure. At the level 3,
	 * there are 8x8x8 possible nodes, so the R, S, T coordinates can be an
	 * integer or a short between 0 and 7. Higher or lower integers are valid
	 * possible but will be out of bounds. Level can be between 0 and 127.
	 */
	public OctNode(short _codeR, short _codeS, short _codeT, byte _level) {
		codeR = _codeR;
		codeS = _codeS;
		codeT = _codeT;
		if (_level < 0 || _level > 127) {
			PApplet.println("ERROR in OctNode: level not in range. Level 0 was used.");
			_level = 0;
		}
		level = _level;
	}

	public OctNode(int _codeR, int _codeS, int _codeT, int _level) {
		codeR = (short) _codeR;
		codeS = (short) _codeS;
		codeT = (short) _codeT;
		if (_level < 0 || _level > 127) {
			PApplet.println("ERROR in OctNode: level not in range. Level 0 was used.");
			_level = 0;
		}
		level = (byte) _level;
	}

	public OctNode(PVector _code, int _level) {
		this.codeR = (short) _code.x;
		this.codeS = (short) _code.y;
		this.codeT = (short) _code.z;
		if (_level < 0 || _level > 127) {
			PApplet.println("ERROR in OctNode: level not in range. Level 0 was used.");
			_level = 0;
		}
		this.level = (byte) _level;
	}

	public OctNode(OctNode _tempNode) {
		codeR = _tempNode.codeR;
		codeS = _tempNode.codeS;
		codeT = _tempNode.codeT;
		level = _tempNode.level;
	}

	public OctNode() {
		codeR = (short) 0;
		codeS = (short) 0;
		codeT = (short) 0;
		level = (byte) 0;
	}
	
	// *************************************************************************************
	// SECTION GETTERS AND SETTERS
	// *************************************************************************************

	public short getLevel() {
		return level;
	}

	public void setLevel(byte _level) {
		if (_level < 0 || _level > 127) {
			PApplet.println("ERROR in setLevel(): level not in range. Level 0 was used.");
			_level = 0;
		}
		this.level = _level;
	}

	public short getCodeR() {
		return codeR;
	}

	public void setCodeR(short _codeR) {
		this.codeR = _codeR;
	}

	public short getCodeS() {
		return codeS;
	}

	public void setCodeS(short _codeS) {
		this.codeS = _codeS;
	}

	public short getCodeT() {
		return codeT;
	}

	public void setCodeT(short _codeT) {
		this.codeT = _codeT;
	}

	public void setCode(short _codeR, short _codeS, short _codeT) {
		this.codeR = _codeR;
		this.codeS = _codeS;
		this.codeT = _codeT;
	}

	public void setCode(int _codeR, int _codeS, int _codeT) {
		this.codeR = (short) _codeR;
		this.codeS = (short) _codeS;
		this.codeT = (short) _codeT;
	}
	
	public void setCode(PVector _code) {
		this.codeR = (short) _code.x;
		this.codeS = (short) _code.y;
		this.codeT = (short) _code.z;
	}

	public PVector getCode() {
		return new PVector(codeR, codeS, codeT);
	}

	// *************************************************************************************
	// RST GEOMETRY
	// *************************************************************************************

	/**
	 * Gets a specific vertex of the node. Can be converted to XYZ coordinated
	 * with '.toXYZ' function.
	 */
	public OctRST getVertex(OCT_VERTEX _v) {
		OctRST vLocation = new OctRST();
		float grid = 1f / (1 << level);
		vLocation.r = (float) (codeR + 0.5 + (_v.r / 2f)) * grid;
		vLocation.s = (float) (codeS + 0.5 + (_v.s / 2f)) * grid;
		vLocation.t = (float) (codeT + 0.5 + (_v.t / 2f)) * grid;
		return vLocation;
	}

	/**
	 * Gets all vertices of the node. Can be converted to XYZ coordinated with
	 * '.toXYZ' function.
	 */
	public ArrayList<OctRST> getVertices() {
		ArrayList<OctRST> tempVertexList = new ArrayList<OctRST>();
		for (OCT_VERTEX v : OCT_VERTEX.values()) {
			tempVertexList.add(getVertex(v));
		}
		return tempVertexList;
	}

	/**
	 * Gets a list of vertices of the node. Can be converted to XYZ coordinated
	 * with '.toXYZ' function.
	 */
	public ArrayList<OctRST> getVertices(ArrayList<OCT_VERTEX> _v) {
		ArrayList<OctRST> tempVertexList = new ArrayList<OctRST>();
		for (OCT_VERTEX v : _v) {
			tempVertexList.add(getVertex(v));
		}
		return tempVertexList;
	}

	/**
	 * Gets the center of the node. Can be converted to XYZ coordinates with
	 * '.toXYZ' function.
	 */
	public OctRST getCenter() {
		OctRST cLocation = new OctRST();
		float grid = 1f / (1 << level);
		// add 0.5 to take the node size into account
		cLocation.r = (float) (codeR + 0.5) * grid;
		cLocation.s = (float) (codeS + 0.5) * grid;
		cLocation.t = (float) (codeT + 0.5) * grid;
		return cLocation;
	}

	public OctRST getFaceCenter(OCT_FACE _f) {
		OctRST faceC = new OctRST();
		float grid = 1f / (1 << level);
		// add 0.5 to take the node size into account
		faceC.r = (float) (codeR + 0.5 + (_f.r / 2f)) * grid;
		faceC.s = (float) (codeS + 0.5 + (_f.s / 2f)) * grid;
		faceC.t = (float) (codeT + 0.5 + (_f.t / 2f)) * grid;
		return faceC;
	}

	public ArrayList<OctRST> getFaceCenters(ArrayList<OCT_FACE> _fList) {
		ArrayList<OctRST> faceCenters = new ArrayList<OctRST>();
		for (OCT_FACE f : _fList) {
			faceCenters.add(getFaceCenter(f));
		}
		return faceCenters;
	}

	public ArrayList<OctRST> getFaceCenters() {
		ArrayList<OctRST> faceCenters = new ArrayList<OctRST>();
		for (OCT_FACE f : OCT_FACE.getAll()) {
			faceCenters.add(getFaceCenter(f));
		}
		return faceCenters;
	}

	public OctRST getEdgeCenter(OCT_EDGE _e) {
		OctRST faceC = new OctRST();
		float grid = 1f / (1 << level);
		// add 0.5 to take the node size into account
		faceC.r = (float) (codeR + 0.5 + _e.r / 2f) * grid;
		faceC.s = (float) (codeS + 0.5 + _e.s / 2f) * grid;
		faceC.t = (float) (codeT + 0.5 + _e.t / 2f) * grid;
		return faceC;
	}

	public ArrayList<OctRST> getEdgeCenters(ArrayList<OCT_EDGE> _eList) {
		ArrayList<OctRST> edgeCenters = new ArrayList<OctRST>();
		for (OCT_EDGE e : _eList) {
			edgeCenters.add(getEdgeCenter(e));
		}
		return edgeCenters;
	}

	public ArrayList<OctRST> getEdgeCenters() {
		ArrayList<OctRST> edgeCenters = new ArrayList<OctRST>();
		for (OCT_EDGE e : OCT_EDGE.getAll()) {
			edgeCenters.add(getEdgeCenter(e));
		}
		return edgeCenters;
	}

	// *************************************************************************************
	// OCTANT AND SIBLINGS
	// *************************************************************************************

	/**
	 * Gets the octant value of the node. See 'oct.enum.OCT_OCTANT' for the list
	 * and naming.
	 */
	public OCT_OCTANT getOctant() {
		if ((codeR + 1) >> 1 == codeR >> 1) {
			if ((codeS + 1) >> 1 == codeS >> 1) {
				if ((codeT + 1) >> 1 == codeT >> 1) {
					return OCT_OCTANT.LBD;
				} else {
					return OCT_OCTANT.LBU;
				}
			} else {
				if ((codeT + 1) >> 1 == codeT >> 1) {
					return OCT_OCTANT.LFD;
				} else {
					return OCT_OCTANT.LFU;
				}
			}
		} else {
			if ((codeS + 1) >> 1 == codeS >> 1) {
				if ((codeT + 1) >> 1 == codeT >> 1) {
					return OCT_OCTANT.RBD;
				} else {
					return OCT_OCTANT.RBU;
				}
			} else {
				if ((codeT + 1) >> 1 == codeT >> 1) {
					return OCT_OCTANT.RFD;
				} else {
					return OCT_OCTANT.RFU;
				}
			}
		}
	}

	/**
	 * Gets a sibling of the node at the specified position in the parent.
	 */
	public OctNode getSibling(OCT_OCTANT _o) {
		int r = 1;
		if (_o.getR() == -1)
			r = 0;
		int s = 1;
		if (_o.getS() == -1)
			s = 0;
		int t = 1;
		if (_o.getT() == -1)
			t = 0;
		OctNode tempNode = new OctNode(((this.codeR >> 1) << 1) + r,
				((this.codeS >> 1) << 1) + s, ((this.codeT >> 1) << 1) + t,
				this.level);
		return tempNode;
	}

	/**
	 * Gets the siblings of the node adjacent to the specified position of the parent.
	 */
	public ArrayList<OctNode> getSiblings(OCT_ENUM _o) {
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		for (OctNode n : this.getParent().getChildren(_o)) {
			if (this.getOctant() != n.getOctant()) {
				tempNodeList.add(n);
			}
		}
		return tempNodeList;
	}
	
	/**
	 * Gets a list of siblings of the node in the specified directions, can be
	 * an edge, face or vertex/octant. Does not check if the siblings are
	 * actually in the octree or not.
	 */
	public ArrayList<OctNode> getSiblings(ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		for (OctNode n : this.getParent().getChildren(_o)) {
			if (this.getOctant() != n.getOctant()) {
				tempNodeList.add(n);
			}
		}
		return tempNodeList;
	}
	
	/**
	 * Gets all siblings of the node, returns 7 values. Does not check if the
	 * siblings are actually in the octree or not.
	 */
	public ArrayList<OctNode> getSiblings() {
		OCT_OCTANT thisOctant = this.getOctant();
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		for (OCT_OCTANT o : OCT_OCTANT.values()) {
			if (thisOctant != o) {
				tempNodeList.add(getSibling(o));
			}
		}
		return tempNodeList;
	}

	// *************************************************************************************
	// PARENT
	// *************************************************************************************

	/**
	 * Gets the single parent of the node, one level up. Does not check if the
	 * parent is actually in the octree or not.
	 */
	public OctNode getParent() {
		return getParent(1);
	}

	/**
	 * Gets the single parent of the node '_l' levels up. Does not check if the
	 * parent is actually in the octree or not.
	 */
	public OctNode getParent(int _l) {
		if (level < _l || _l <= 0) {
			PApplet.println("ERROR in getParent(): parent level cannot be 0 or negative. 'Null' returned.");
			return null;
		}
		OctNode tempNode = new OctNode(codeR >> _l, codeS >> _l, codeT >> _l,
				level - _l);
		return tempNode;
	}

	/**
	 * Gets all parents of the node from the node level (not included) up to
	 * '_l' levels up. Does not check if the parents are actually in the octree
	 * or not.
	 */
	public ArrayList<OctNode> getAllParents(int _l) {
		int up = _l;
		/*
		 * if (up <= 0) { PApplet.println(
		 * "ERROR in getAllParents(): level not in range. 'Null' returned.");
		 * return null; }
		 */
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		for (int i = 1; i <= up; i++) {
			tempNodeList.add(getParent(i));
		}
		return tempNodeList;
	}

	/**
	 * Gets all parents of the node from the node level (not included) up to the
	 * level 0 included. Does not check if the parents are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllParents() {
		int up = this.level;
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		for (int i = 1; i <= up; i++) {
			tempNodeList.add(getParent(i));
		}
		return tempNodeList;
	}

	/**
	 * Gets all parents of the node from the node level (not included) up to the
	 * octree minimum depth. Does not check if the parents are actually in the
	 * octree or not. The octree is there just to specify the level bounds.
	 */
	public ArrayList<OctNode> getAllParents(OctOctree _octree) {
		ArrayList<OctNode> tempNodeList = new ArrayList<OctNode>();
		if (this.level < _octree.minD) {
			PApplet.println("ERROR in getAllParents(octree): node outside of octree.");
			return tempNodeList;
		}
		for (int i = 1; i <= this.level - _octree.minD; i++) {
			tempNodeList.add(getParent(i));
		}
		return tempNodeList;
	}

	// *************************************************************************************
	// CHILDREN
	// *************************************************************************************

	/**
	 * Gets one specific child of the node, one level down. Does not check if
	 * the child is actually in the octree or not. If the level is lower than
	 * 127 will cause an error, but such a small level should probably never
	 * happen.
	 */
	public OctNode getChild(OCT_OCTANT _o) {
		int r = 1;
		if (_o.r == -1)
			r = 0;
		int s = 1;
		if (_o.s == -1)
			s = 0;
		int t = 1;
		if (_o.t == -1)
			t = 0;
		OctNode tempNode = new OctNode((codeR << 1) + r, (codeS << 1) + s,
				(codeT << 1) + t, level + 1);
		return tempNode;
	}

	/**
	 * Gets the children of the node, one level down, according to a list of
	 * octants. Does not check if the children are actually in the octree or
	 * not.
	 */
	public ArrayList<OctNode> getChildren(ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		for (OCT_ENUM o : _o) {
			tempChildList.addAll(getChildren(1, o));
		}
		return tempChildList;
	}

	/**
	 * Gets the children of the node, one level down, that are adjacent to one
	 * vertex, face or edge of the parent. Does not check if the children are
	 * actually in the octree or not. A face returns 4 nodes, an edge 2.
	 */
	public ArrayList<OctNode> getChildren(OCT_ENUM _o) {
		return getChildren(1, _o);
	}

	/**
	 * Gets the children of the node, '_l' level down, that are adjacent to a
	 * list if vertex, face or edge of the parent. Does not check if the
	 * children are actually in the octree or not.
	 */
	public ArrayList<OctNode> getChildren(int _l, ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		for (OCT_ENUM o : _o) {
			tempChildList.addAll(getChildren(_l, o));
		}
		return tempChildList;
	}

	/**
	 * Gets the children of the node, '_l' level down, that are adjacent to one
	 * vertex, face or edge of the parent. Does not check if the children are
	 * actually in the octree or not.
	 */
	public ArrayList<OctNode> getChildren(int _l, OCT_ENUM _e) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (_l <= 0) {
			return tempChildList;
		}
		tempChildList.add(this);
		return getChildrenRecursion(_l, _e, tempChildList);
	}

	private ArrayList<OctNode> getChildrenRecursion(int _l, OCT_ENUM _e,
			ArrayList<OctNode> _nodeList) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		for (OctNode n : _nodeList) {
			for (OCT_OCTANT o : _e.getOctants()) {
				tempChildList.add(n.getChild(o));
			}
		}
		if (_l == 1) {
			return tempChildList;
		} else {
			_l -= 1;
			return getChildrenRecursion(_l, _e, tempChildList);
		}
	}

	/**
	 * Gets the children of the node, '_l' level down. Does not check if the
	 * children are actually in the octree or not.
	 */
	public ArrayList<OctNode> getChildren(int _l) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (_l <= 0) {
			return tempChildList;
		}
		for (int m = (int) 1 << _l; --m >= 0;) {
			for (int n = (int) 1 << _l; --n >= 0;) {
				for (int p = (int) 1 << _l; --p >= 0;) {
					OctNode tempNode = new OctNode((codeR << _l) + m,
							(codeS << _l) + n, (codeT << _l) + p, level + _l);
					tempChildList.add(tempNode);
				}
			}
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node, one level down. Does not check if the
	 * children are actually in the octree or not.
	 */
	public ArrayList<OctNode> getChildren() {
		return getChildren(1);
	}

	/**
	 * Gets all children of the node from the node level (not included) up to
	 * '_l' levels down. Does not check if the children are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllChildren(int _l) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (_l <= 0) {
			return tempChildList;
		}
		for (int i = 1; i <= _l; i++) {
			tempChildList.addAll(getChildren(i));
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node up from the node level (not included) up to
	 * '_l' levels down. Does not check if the children are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllChildren(int _l, ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (_l <= 0) {
			return tempChildList;
		}
		for (int i = 1; i <= _l; i++) {
			for (OCT_ENUM o : _o) {
				tempChildList.addAll(getChildren(i, o));
			}
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node from the node level (not included) up to
	 * '_l' levels down. Does not check if the children are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllChildren(int _l, OCT_ENUM _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (_l <= 0) {
			// PApplet.println("ERROR in getAllChildren(): level not in range. 'Null' returned.");
			return tempChildList;
		}
		for (int i = 1; i <= _l; i++) {
			tempChildList.addAll(getChildren(i, _o));
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node up from the node level (not included) up to
	 * the max depth or the octree. Does not check if the children are actually
	 * in the octree or not. If the node is at the maximum level of the octree,
	 * returns null because no child is possible.
	 */
	public ArrayList<OctNode> getAllChildren(OctOctree _octree,
			ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllChildren(octree): node outside of octree.");
			return tempChildList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			for (OCT_ENUM o : _o) {
				tempChildList.addAll(getChildren(i, o));
			}
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node from the node level (not included) up to
	 * the max depth or the octree. Does not check if the children are actually
	 * in the octree or not.
	 */
	public ArrayList<OctNode> getAllChildren(OctOctree _octree) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllChildren(octree): node outside of octree.");
			return tempChildList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			tempChildList.addAll(getChildren(i));
		}
		return tempChildList;
	}

	/**
	 * Gets all children of the node from the node level (not included) up to
	 * the max depth or the octree. Does not check if the children are actually
	 * in the octree or not.
	 */
	public ArrayList<OctNode> getAllChildren(OctOctree _octree, OCT_ENUM _o) {
		ArrayList<OctNode> tempChildList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllChildren(octree): node outside of octree.");
			return tempChildList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			tempChildList.addAll(getChildren(i, _o));
		}
		return tempChildList;
	}

	// *************************************************************************************
	// NEIGHBORS (SAME LEVEL)
	// *************************************************************************************

	/**
	 * Gets the single neighbor of the node, moving in the specified direction
	 * and at the same level. If _dirR, _dirS, _dirT are different than +1, 0,
	 * or -1, the node will not be a neighbor! But it will still be returned.
	 * Does not check if the neighbor is in the Octree or not!
	 */
	public OctNode getNbr(OCT_ENUM _e) {
		OctNode tempNode = new OctNode(codeR, codeS, codeT, level);
		tempNode.codeR += _e.getR();
		tempNode.codeS += _e.getS();
		tempNode.codeT += _e.getT();
		return tempNode;
	}

	/**
	 * Gets the neighbors of the node at the same level and in the specified
	 * directions. Does not check if the neighbors are in the Octree or not!
	 */
	public ArrayList<OctNode> getNbrs(ArrayList<? extends OCT_ENUM> _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : _n) {
			t.add(getNbr(e));
		}
		return t;
	}

	/**
	 * Gets the neighbors of the node at the same level. Does not check if the
	 * neighbors are in the Octree or not! Returns the 26 nodes surrounding the
	 * node.
	 */
	public ArrayList<OctNode> getNbrs() {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : neighborEnums()) {
			t.add(getNbr(e));
		}
		return t;
	}

	private final static ArrayList<OCT_ENUM> neighborEnums() {
		ArrayList<OCT_ENUM> eList = new ArrayList<OCT_ENUM>();
		eList.addAll(OCT_VERTEX.getAll());
		eList.addAll(OCT_FACE.getAll());
		eList.addAll(OCT_EDGE.getAll());
		return eList;
	}

	// *************************************************************************************
	// SMALLER NEIGHBORS (HIGHER IN THE OCTREE)
	// *************************************************************************************

	/**
	 * Gets the neighbors of the node 1 level down and in the specified
	 * directions. Does not check if the neighbors are in the Octree or not!
	 */
	public ArrayList<OctNode> getSmallerNbrs(OCT_ENUM _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		t.addAll(getSmallerNbrs(1, _n));
		return t;
	}
	
	/**
	 * Gets the neighbors of the node 1 level down and in the specified
	 * directions. Does not check if the neighbors are in the Octree or not!
	 */
	public ArrayList<OctNode> getSmallerNbrs(ArrayList<? extends OCT_ENUM> _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		t.addAll(getSmallerNbrs(1, _n));
		return t;
	}

	/**
	 * Gets the neighbors of the node at the same level. Does not check if the
	 * neighbors are in the Octree or not! Returns the 26 nodes surrounding the
	 * node.
	 */
	public ArrayList<OctNode> getSmallerNbrs() {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : neighborEnums()) {
			t.addAll(getSmallerNbrs(1,e));
		}
		return t;
	}

	/**
	 * Gets the neighbors of the node for the specified level and OCT_ENUM. Does
	 * not check if the neighbors are in the Octree or not! Too get all
	 * intermediate levels, use getAllSmallerNbrs!
	 */
	public ArrayList<OctNode> getSmallerNbrs(int _l, OCT_ENUM _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_l > 0) {
			ArrayList<OctNode> tempNodeList = this.getChildren(_l, _n);
			for (OctNode tempNode2 : tempNodeList) {
				t.add(tempNode2.getNbr(_n));
			}
			return t;
		}
		// PApplet.println("ERROR in getSmallerNbrs(): level not in range. 'Null' returned.");
		return t;
	}

	/**
	 * Gets the all the neighbors of the node at the specified level. Does not
	 * check if the neighbors are in the Octree or not! Too get all intermediate
	 * levels, use getAllSmallerNbrs!
	 */
	public ArrayList<OctNode> getSmallerNbrs(int _l) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_l > 0) {
			for (OCT_ENUM e : neighborEnums()) {
				t.addAll(getSmallerNbrs(_l, e));
			}
			return t;
		}
		// PApplet.println("ERROR in getSmallerNbrs(): level not in range. 'Null' returned.");
		return t;
	}

	/**
	 * Gets the neighbors of the node for the specified OCT_ENUM list. Does not
	 * check if the neighbors are in the Octree or not! Too get all intermediate
	 * levels, use getAllSmallerNbrs!
	 */
	public ArrayList<OctNode> getSmallerNbrs(int _l, ArrayList<? extends OCT_ENUM> _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_l > 0) {
			for (OCT_ENUM e : _n) {
				t.addAll(getSmallerNbrs(_l, e));
			}
			return t;
		}
		// PApplet.println("ERROR in getSmallerNbrs(): level not in range. 'Null' returned.");
		return t;
	}
	
	/**
	 * Gets all neighbors of the node from the node level (not included) up to
	 * '_l' levels down. Does not check if the neighbors are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(int _l) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (_l <= 0) {
			// PApplet.println("ERROR in getAllSmallerNbrs(): level not in range. 'Null' returned.");
			return tempList;
		}
		for (int i = 1; i <= _l; i++) {
			tempList.addAll(getSmallerNbrs(i));
		}
		return tempList;
	}

	/**
	 * Gets all neighbors of the node up from the node level (not included) up
	 * to '_l' levels down. Does not check if the neighbors are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(int _l, ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (_l <= 0) {
			// PApplet.println("ERROR in getAllSmallerNbrs(): level not in range. 'Null' returned.");
			return tempList;
		}
		for (int i = 1; i <= _l; i++) {
			for (OCT_ENUM o : _o) {
				tempList.addAll(getSmallerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets all neighbors of the node from the node level (not included) up to
	 * '_l' levels down. Does not check if the neighbors are actually in the
	 * octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(int _l, OCT_ENUM _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (_l <= 0) {
			// PApplet.println("ERROR in getAllSmallerNbrs(): level not in range. 'Null' returned.");
			return tempList;
		}
		for (int i = 1; i <= _l; i++) {
			tempList.addAll(getSmallerNbrs(i, _o));
		}
		return tempList;
	}

	/**
	 * Gets all neighbors of the node up from the node level (not included) up
	 * to the max depth or the octree. Does not check if the neighbors are
	 * actually in the octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(OctOctree _octree,
			ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllSmallerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			for (OCT_ENUM o : _o) {
				tempList.addAll(getSmallerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets all neighbors of the node from the node level (not included) up to
	 * the max depth or the octree. Does not check if the neighbors are actually
	 * in the octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(OctOctree _octree) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllSmallerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			tempList.addAll(getSmallerNbrs(i));
		}
		return tempList;
	}

	/**
	 * Gets all neighbors of the node from the node level (not included) up to
	 * the max depth or the octree. Does not check if the neighbors are actually
	 * in the octree or not.
	 */
	public ArrayList<OctNode> getAllSmallerNbrs(OctOctree _octree, OCT_ENUM _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level > _octree.maxD) {
			PApplet.println("ERROR in getAllSmallerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= _octree.maxD - this.level; i++) {
			tempList.addAll(getSmallerNbrs(i, _o));
		}
		return tempList;
	}

	// *************************************************************************************
	// BIGGER NEIGHBORS (LOWER IN THE OCTREE)
	// *************************************************************************************

	/**
	 * Gets the node that is bigger (lower level) and touching the node. Be
	 * careful with this method as it can return a null. Does not check if the
	 * node is in the octree or not!
	 */
	public OctNode getBiggerNbr(OCT_ENUM _enum) {
		if (level == 0) {
			PApplet.println("ERROR in getBiggerNbr(): node level cannot be 0 or negative. 'Null' returned.");
			return null;
		}
		OctNode n = (this.getNbr(_enum)).getParent(1);
		if (!n.isParentOf(this))
			return n;
		else
			PApplet.println("ERROR in getBiggerNbr(): no possible node. 'Null' returned.");
		return null;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getBiggerNbrs() {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : neighborEnums()) {
			t.add(this.getBiggerNbr(e));
		}
		return t;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getBiggerNbrs(ArrayList<OCT_ENUM> _o) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : _o) {
			t.add(this.getBiggerNbr(e));
		}
		return t;
	}

	/**
	 * Gets the node that is bigger (lower level) and touching the node. Can be
	 * a parent of the node. Does not check if the node is in the octree or not!
	 */
	public ArrayList<OctNode> getBiggerNbrs(int _l, OCT_ENUM _enum) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (level < _l || _l <= 0) {
			// PApplet.println("ERROR in getBiggerNbrs(): level not in range. 'Null' returned.");
			return t;
		}
		OctNode tempParentNode = (this.getNbr(_enum)).getParent(_l);
		if (!tempParentNode.isParentOf(this))
			t.add(tempParentNode);
		return t;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getBiggerNbrs(int _l) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : neighborEnums()) {
			t.addAll(getBiggerNbrs(_l, e));
		}
		return t;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getBiggerNbrs(int _l, ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_ENUM e : _o) {
			t.addAll(getBiggerNbrs(_l, e));
		}
		return t;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(int _l, ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		for (int i = 1; i <= _l; i++) {
			for (OCT_ENUM o : _o) {
				tempList.addAll(getBiggerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(int _l, OCT_ENUM _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		for (int i = 1; i <= _l; i++) {
			tempList.addAll(getBiggerNbrs(i, _o));
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(int _l) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		for (int i = 1; i <= _l; i++) {
			for (OCT_ENUM o : neighborEnums()) {
				tempList.addAll(getBiggerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs() {
		int _l = this.level;
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		for (int i = 1; i <= _l; i++) {
			for (OCT_ENUM o : neighborEnums()) {
				tempList.addAll(getBiggerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(OctOctree _octree,
			ArrayList<? extends OCT_ENUM> _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level < _octree.minD) {
			PApplet.println("ERROR in getAllBiggerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= this.level - _octree.minD; i++) {
			for (OCT_ENUM o : _o) {
				tempList.addAll(getBiggerNbrs(i, o));
			}
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(OctOctree _octree, OCT_ENUM _o) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level < _octree.minD) {
			PApplet.println("ERROR in getAllBiggerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= this.level - _octree.minD; i++) {
			tempList.addAll(getBiggerNbrs(i, _o));
		}
		return tempList;
	}

	/**
	 * Gets the nodes that are bigger (lower level) and touching the node. Can
	 * be parents of the nodes. Does not check if the nodes are in the octree or
	 * not!
	 */
	public ArrayList<OctNode> getAllBiggerNbrs(OctOctree _octree) {
		ArrayList<OctNode> tempList = new ArrayList<OctNode>();
		if (this.level < _octree.minD) {
			PApplet.println("ERROR in getAllBiggerNbrs(octree): node outside of octree.");
			return tempList;
		}
		for (int i = 1; i <= this.level - _octree.minD; i++) {
			for (OCT_ENUM o : neighborEnums()) {
				tempList.addAll(getBiggerNbrs(i, o));
			}
		}
		return tempList;
	}

	// *************************************************************************************
	// SWEEP
	// *************************************************************************************
	
	/**
	 * Generates all nodes that might be touching the node, at the same level,
	 * higher level or lower level. The octree is provided to give bounds to the
	 * list.
	 */
	public ArrayList<OctNode> sweep(OctOctree _octree, boolean _doSame,
			boolean _doUp, boolean _doDown) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_doSame)
			t.addAll(getNbrs());
		if (_doUp)
			t.addAll(getAllBiggerNbrs(this.level - _octree.minD));
		if (_doDown)
			t.addAll(getAllSmallerNbrs(_octree.maxD - this.level));
		return t;
	}

	/**
	 * Generates all nodes that might be touching the node, at the same level,
	 * higher level or lower level. The octree is provided to give bounds to the
	 * list.
	 */
	public ArrayList<OctNode> sweep(OctOctree _octree, boolean _doSame,
			boolean _doUp, boolean _doDown, OCT_ENUM _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_doSame)
			t.add(getNbr(_n));
		if (_doUp)
			t.addAll(getAllBiggerNbrs(this.level - _octree.minD, _n));
		if (_doDown)
			t.addAll(getAllSmallerNbrs(_octree.maxD - this.level, _n));
		return t;
	}

	/**
	 * Generates all nodes that might be touching the node, at the same level,
	 * higher level or lower level. The octree is provided to give bounds to the
	 * list.
	 */
	public ArrayList<OctNode> sweep(OctOctree _octree, boolean _doSame,
			boolean _doUp, boolean _doDown, ArrayList<OCT_ENUM> _n) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		if (_doSame)
			t.addAll(getNbrs(_n));
		if (_doUp)
			t.addAll(getAllBiggerNbrs(this.level - _octree.minD, _n));
		if (_doDown)
			t.addAll(getAllSmallerNbrs(_octree.maxD - this.level, _n));
		return t;
	}

	// *************************************************************************************
	// BOX AND HEXA
	// *************************************************************************************

	/**
	 * Computes the volume of the node based on its 8 vertices without any
	 * deformations. In XYZ units. Octree needs to be provided!
	 */
	public float volumeBox(OctOctree _octree) {
		float volume = _octree.getDimension().x / (1 << level)
				* _octree.getDimension().y / (1 << level)
				* _octree.getDimension().z / (1 << level);
		return volume;
	}

	public OctRST getMin() {
		OctRST m = this.getVertex(OCT_VERTEX.LBD);
		OctRST n = this.getVertex(OCT_VERTEX.RFU);
		OctRST min = new OctRST(PApplet.min(m.r,n.r),PApplet.min(m.s,n.s),PApplet.min(m.t,n.t));
		return min;
	}
	
	public OctRST getMax() {
		OctRST m = this.getVertex(OCT_VERTEX.LBD);
		OctRST n = this.getVertex(OCT_VERTEX.RFU);
		OctRST max = new OctRST(PApplet.max(m.r,n.r),PApplet.max(m.s,n.s),PApplet.max(m.t,n.t));
		return max;
	}
	
	/**
     * ADAPTED FROM https://gist.github.com/aadnk/7123926.
     * <p>
     * Calculates intersection with the given ray between a certain distance
     * interval.
     * <p>
     * Ray-box intersection is using IEEE numerical properties to ensure the
     * test is both robust and efficient, as described in:
     * <br>
     * <code>Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley: "An
     * Efficient and Robust Ray-Box Intersection Algorithm" Journal of graphics
     * tools, 10(1):49-54, 2005</code>
     */
    public OctRST intersectPoint(OctRST start, OctRST direction, float minDist, float maxDist) {
    	direction = direction.normalize();
    	OctRST invDir = new OctRST(1f / direction.r, 1f / direction.s, 1f / direction.t);
   	
    	OctRST max = this.getMax();
    	OctRST min = this.getMin();
   		
        boolean signDirX = invDir.r < 0;
        boolean signDirY = invDir.s < 0;
        boolean signDirZ = invDir.t < 0;

        OctRST bbox = signDirX ? max : min;
        double tmin = (bbox.r-start.r) * invDir.r;
        bbox = signDirX ? min : max;
        double tmax = (bbox.r-start.r) * invDir.r;
        bbox = signDirY ? max : min;
        double tymin = (bbox.s-start.s) * invDir.s;
        bbox = signDirY ? min : max;
        double tymax = (bbox.s-start.s) * invDir.s;
 
      //  PApplet.println(tmin);
      //  PApplet.println(tmax);
        
        if ((tmin > tymax) || (tymin > tmax)) {
            return null;
        }
        if (tymin > tmin) {
            tmin = tymin;
        }
        if (tymax < tmax) {
            tmax = tymax;
        }

        bbox = signDirZ ? max : min;
        double tzmin = (bbox.t-start.t) * invDir.t;
        bbox = signDirZ ? min : max;
        double tzmax = (bbox.t-start.t) * invDir.t;

        if ((tmin > tzmax) || (tzmin > tmax)) {
            return null;
        }
        if (tzmin > tmin) {
            tmin = tzmin;
        }
        if (tzmax < tmax) {
            tmax = tzmax;
        }
       // PApplet.println(tmin);
       // PApplet.println(tmax);
        if ((tmin < maxDist) && (tmax > minDist)) {
            return new OctRST(direction.normalize().scale((float) tmin).add(start));
        }
        return null;
    }

	/**
	 * Computes the volume of the node based on its 8 vertices, also valid if
	 * deformed. Based on "Efficient computation of volume of hexahedral cells",
	 * J. Grandy. http://www.osti.gov/scitech/servlets/purl/632793
	 */
	public float volumeHexa(OctOctree _octree) {
		float volume = 0;
		OctXYZ v1 = this.getVertex(OCT_VERTEX.get(6)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(0)).toXYZ(_octree));
		OctXYZ v2 = this.getVertex(OCT_VERTEX.get(3)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(0)).toXYZ(_octree));
		OctXYZ v3 = this.getVertex(OCT_VERTEX.get(2)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(7)).toXYZ(_octree));
		OctXYZ v4 = this.getVertex(OCT_VERTEX.get(4)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(0)).toXYZ(_octree));
		OctXYZ v5 = this.getVertex(OCT_VERTEX.get(7)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(5)).toXYZ(_octree));
		OctXYZ v6 = this.getVertex(OCT_VERTEX.get(1)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(0)).toXYZ(_octree));
		OctXYZ v7 = this.getVertex(OCT_VERTEX.get(5)).toXYZ(_octree)
				.sub(this.getVertex(OCT_VERTEX.get(2)).toXYZ(_octree));
		volume = PApplet.abs((v1.dot(v2.cross(v3)) + v1.dot(v4.cross(v5)) + v1
				.dot(v6.cross(v7))) / 6);
		return volume;
	}

	// *************************************************************************************
	// BOOLEAN CHECKS
	// *************************************************************************************

	/**
	 * Returns true if the node code is within the bounds of the octree. Returns
	 * false for nodes that have a negative or too large R, S or T code.
	 */
	public boolean isCodeWB(OctOctree _octree) {
		if (this.codeR >= 0 && this.codeS >= 0 && this.codeT >= 0
				&& this.codeR < 1 << level && this.codeS < 1 << level
				&& this.codeT < 1 << level)
			return true;
		else
			return false;
	}

	/**
	 * Returns true if the node level is within the current bounds of the octree
	 * (between min. and max. level).
	 */
	public boolean isLevelWB(OctOctree _octree) {
		if (this.level <= _octree.getMaxDepth()
				&& this.level >= _octree.getMinDepth())
			return true;
		else
			return false;
	}

	/**
	 * Returns true if 'ifLevelWB' AND 'ifCodeWB' are true
	 */
	public boolean isWB(OctOctree _octree) {
		if (this.isCodeWB(_octree) && this.isLevelWB(_octree))
			return true;
		else
			return false;
	}
	
	/**
	 * Returns true if node is selected.
	 */	
	public boolean isSelected() {
		if (this.isSelected)
			return true;
		else
			return false;
	}
	

	/**
	 * Returns true if the node is included in the '_t' node (check if the node
	 * is a child of _t). Is false if the two nodes have the same level or are
	 * equals.
	 */
	public boolean isChildOf(OctNode _t) {
		if (this.level == _t.level) {
			if (this.equals(_t))
				return false;
		} else if (this.level > _t.level) {
			if (this.getParent(this.level - _t.level).equals(_t))
				return true;
			else
				return false;
		} else if (this.level < _t.level)
			return false;
		return false;
	}

	/**
	 * Returns true if the node is a parent of '_t' (check if the node is a
	 * parent of _t). Is false if the two nodes have the same level.
	 */
	public boolean isParentOf(OctNode _t) {
		if (this.level == _t.level) {
			return false;
		} else if (this.level > _t.level) {
			return false;
		} else if (this.level < _t.level) {
			if (_t.getParent(_t.level - this.level).equals(this))
				return true;
			else
				return false;
		}
		return false;
	}

	/**
	 * Returns true if the node is a sibling of '_t' (check if the nodes have
	 * the same parent). Is false if the two nodes don't have the same level or
	 * are equals.
	 */
	public boolean isSiblingOf(OctNode _t) {
		if (this.getParent(1).equals(_t.getParent(1))
				&& this.equals(_t) == false) {
			return true;
		} else
			return false;
	}

	/**
	 * Returns true if the node is a neighbor of '_t'. Returns false if the node
	 * '_t' equals the node.
	 */
	public boolean isNeighborOf(OctNode _t) {
		OctNode t = new OctNode(_t);
		if (this.level == t.level) {
			if (this.getNbrs().contains(t) && this != _t)
				return true;
		}
		if (this.level > t.level) {
			for (OctNode temp : this.getNbrs()) {
				if (temp.isChildOf(t) && (this.isChildOf(t) == false))
					return true;
			}
		}
		if (this.level < t.level) {
			for (OctNode temp : t.getNbrs()) {
				if (temp.isChildOf(this) && (t.isChildOf(this) == false))
					return true;
			}
		}
		return false;
	}

	// *************************************************************************************
	// COMMON VERTICES, EDGES, FACES
	// *************************************************************************************

	/**
	 * Returns the other sibling node adjacent to the specified edge. Null if
	 * the edge is not touching the node.
	 * 
	 * @param _e
	 *            Edge to provide
	 */
	public OctNode getOtherNode(OCT_EDGE _e) {
		for (OCT_OCTANT e : _e.getOctants()) {
			if (this.getOctant() != e) {
				return this.getParent().getChild(e);
			}
		}
		return null;
	}

	/**
	 * Returns the three other sibling nodes adjacent to the specified face.
	 * Null if the edge is not touching the node.
	 * 
	 * @param _e
	 *            Edge to provide
	 */
	public ArrayList<OctNode> getOtherNodes(OCT_FACE _f) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OCT_OCTANT e : _f.getOctants()) {
			if (this.getOctant() != e) {
				t.add(this.getParent().getChild(e));
			}
		}
		return t;
	}

	// TODO generalize for nodes of different levels and not siblings?
	/**
	 * Returns the parent's edge that is shared with the sibling. Returns null
	 * if the nodes are not siblings.
	 */
	public OCT_EDGE getSharedEdge(OctNode _withThisSibling) {
		if (this.isSiblingOf(_withThisSibling)) {
			int thisOctant = this.getOctant().getOrdinal();
			int siblingOctant = _withThisSibling.getOctant().getOrdinal();
			return OctTables.commonEdgeLUT[siblingOctant][thisOctant];
		} else
			return null;
	}

	// TODO generalize for nodes of different levels and not siblings!
	/**
	 * Returns the parent's face that is shared with the sibling. Returns null
	 * if the nodes are not siblings.
	 */
	public OCT_FACE getSharedFace(OctNode _withThisSibling) {
		if (this.isSiblingOf(_withThisSibling)) {
			int tableShift = 12;
			int thisOctant = this.getOctant().getOrdinal();
			int siblingOctant = _withThisSibling.getOctant().getOrdinal();
			return OctTables.commonFaceLUT[siblingOctant + tableShift][thisOctant];
		} else
			return null;
	}

	// *************************************************************************************
	// DRAW
	// *************************************************************************************

	public void drawCenter(OctOctree _octree) {
		PVector center = getCenter().toXYZ(_octree).toPVector();
		_octree.p5.beginShape(PConstants.POINTS);
		_octree.p5.vertex(center.x, center.y, center.z);
		_octree.p5.endShape();
	}

	public void drawVertices(OctOctree _octree) {
		_octree.p5.beginShape(PConstants.POINTS);
		for (OctRST t : getVertices()) {
			_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
					t.toXYZ(_octree).z);
		}
		_octree.p5.endShape();
	}

	public void drawEdges(OctOctree _octree) {
		_octree.p5.beginShape(PConstants.LINES);
		for (OCT_EDGE e : OCT_EDGE.values()) {
			for (OctRST t : getVertices(e.getVertices())) {
				_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
						t.toXYZ(_octree).z);
			}
		}
		_octree.p5.endShape();
	}

	public void drawFaces(OctOctree _octree) {
		_octree.p5.beginShape(PConstants.QUADS);
		for (OCT_FACE f : OCT_FACE.values()) {
			for (OctRST t : getVertices(f.getVertices())) {
				_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
						t.toXYZ(_octree).z);
			}
		}
		_octree.p5.endShape();
	}

	public void drawVertex(OctOctree _octree, OCT_VERTEX _v) {
		_octree.p5.beginShape(PConstants.POINTS);
		OctRST t = this.getVertex(_v);
		_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
				t.toXYZ(_octree).z);
		_octree.p5.endShape();
	}

	public void drawEdge(OctOctree _octree, OCT_EDGE _e) {
		_octree.p5.beginShape(PConstants.LINES);
		for (OctRST t : getVertices(_e.getVertices())) {
			_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
					t.toXYZ(_octree).z);
		}
		_octree.p5.endShape();
	}

	public void drawFace(OctOctree _octree, OCT_FACE _f) {
		_octree.p5.beginShape(PConstants.QUADS);
		for (OctRST t : getVertices(_f.getVertices())) {
			_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
					t.toXYZ(_octree).z);
		}
		_octree.p5.endShape();
	}

	// *************************************************************************************
	// HASH CODE
	// *************************************************************************************

	public int hashCode() {
		return new String("r" + codeR + "s" + codeS + "t" + codeT + "l" + level)
				.hashCode();
	}

	public boolean equals(Object _t) {
		OctNode t = (OctNode) _t;
		if (t.level == this.level && t.codeR == this.codeR
				&& t.codeS == this.codeS && t.codeT == this.codeT)
			return true;
		else
			return false;
	}
} public class OctOctree implements java.io.Serializable {

	// *************************************************************************************
	// VARIABLES
	// *************************************************************************************

	protected transient PApplet p5;

	// xyz dimensions
	protected OctXYZ dimension = new OctXYZ();
	protected OctXYZ origin = new OctXYZ();

	// octree structure

	protected int minD = 0; // minimum depth, level 1 equals to 8 nodes
	protected int maxD = 8; // max depth, limits the scope of some calculations

	// list
	public HashSet<OctNode> nodeList;
	public HashSet<OctNode> selectedNodes = new HashSet<OctNode>();

	// *************************************************************************************
	// CONSTRUCTORS
	// *************************************************************************************

	/**
	 * Creates a new Octree. A default origin of (0,0,0) is assumed.
	 */
	public OctOctree(PApplet _p5, float _dimX, float _dimY, float _dimZ) {
		p5 = _p5;
		dimension.x = _dimX;
		dimension.y = _dimY;
		dimension.z = _dimZ;
		nodeList = new HashSet<OctNode>();
	}
	
	/**
	 * Creates a new Octree.
	 */
	public OctOctree(PApplet _p5, OctOctree _octree) {
		p5 = _p5;
		dimension = _octree.dimension;
		origin = _octree.origin;
		minD = _octree.minD;
		maxD = _octree.maxD;
		nodeList = _octree.nodeList;
		selectedNodes = _octree.selectedNodes;
	}

	/**
	 * Creates a new Octree with a user-defined dimension and origin.
	 */
	public OctOctree(PApplet _p5, float _dimX, float _dimY, float _dimZ,
			float _originX, float _originY, float _originZ) {
		p5 = _p5;
		dimension.x = _dimX;
		dimension.y = _dimY;
		dimension.z = _dimZ;
		origin.x = _originX;
		origin.y = _originY;
		origin.z = _originZ;
		nodeList = new HashSet<OctNode>();
	}

	public OctOctree(PApplet _p5, PVector _dimension, PVector _origin) {
		new OctOctree(_p5, _dimension.x, _dimension.y, _dimension.z, _origin.x,
				_origin.y, _origin.z);
	}

	public OctOctree(PApplet _p5, OctXYZ _dimension, OctXYZ _origin) {
		new OctOctree(_p5, _dimension.x, _dimension.y, _dimension.z, _origin.x,
				_origin.y, _origin.z); 
	}

	public OctOctree(PApplet _p5, PVector _dimension) {
		new OctOctree(_p5, _dimension.x, _dimension.y, _dimension.z);
	}

	public OctOctree(PApplet _p5, OctXYZ _dimension) {
		new OctOctree(_p5, _dimension.x, _dimension.y, _dimension.z);
	}

	// *************************************************************************************
	// GETTERS AND SETTERS
	// *************************************************************************************

	/**
	 * This resizes the octree. It is a scaling of the whole existing structure.
	 */
	public void setDimensions(OctXYZ _dimension) {
		dimension = _dimension;
	}

	public void setDimensions(PVector _dimension) {
		dimension.x = _dimension.x;
		dimension.y = _dimension.y;
		dimension.z = _dimension.z;
	}

	public void setDimensions(float _x, float _y, float _z) {
		dimension.x = _x;
		dimension.y = _y;
		dimension.z = _z;
	}

	/**
	 * This changes the origin of the octree. It is a translation of the whole
	 * existing structure.
	 */
	public void setOrigin(OctXYZ _origin) {
		origin = _origin;
	}

	public void setOrigin(float _x, float _y, float _z) {
		origin.x = _x;
		origin.y = _y;
		origin.z = _z;
	}

	public void setOrigin(PVector _origin) {
		origin.x = _origin.x;
		origin.y = _origin.y;
		origin.z = _origin.z;
	}

	public OctXYZ getMin() {
		OctXYZ t = origin.add(dimension);
		OctXYZ min = new OctXYZ(PApplet.min(origin.x,t.x),PApplet.min(origin.y,t.y),PApplet.min(origin.z,t.z));
		return min;
	}

	public OctXYZ getMax() {
		OctXYZ t = origin.add(dimension);
		OctXYZ max = new OctXYZ(PApplet.max(origin.x,t.x),PApplet.max(origin.y,t.y),PApplet.max(origin.z,t.z));
		return max;
	}	

	/**
	 * This changes the center of the octree. It is a translation of the whole
	 * existing structure. 
	 */
	public void setCenter(OctXYZ _center) {
		origin = _center.sub(dimension.scale(0.5f));
	}

	public void setCenter(float _x, float _y, float _z) {
		origin.x = _x - dimension.x / 2;
		origin.y = _y - dimension.y / 2;
		origin.z = _z - dimension.z / 2;
	}

	public void setCenter(PVector _center) {
		origin.x = _center.x - dimension.x / 2;
		origin.y = _center.y - dimension.y / 2;
		origin.z = _center.z - dimension.z / 2;
	}

	// TODO a method to create a cropping instead of a scaling

	/**
	 * Changes the minimum depth of the octree. It might be necessary to clean
	 * the octree afterwards if you use this command with nodes already in the
	 * tree. If you increase the minimum level, some of the existing nodes might
	 * be too large. Use algCleanLevels() to remove those.
	 */
	public void setMinDepth(int _minD) {
		if (_minD <= maxD && _minD >= 0) {
			minD = _minD;
		} else {
			PApplet.println("ERROR in setMinDepth(): level not in range. Level not changed.");
		}

	}

	/**
	 * Changes the maximum depth of the octree. It might be necessary to clean
	 * the octree afterwards if you use this command with nodes already in the
	 * tree. If you decrease the maximum level, some of the existing nodes might
	 * be too small and create errors. Use algCleanLevels() to remove those.
	 */
	public void setMaxDepth(int _maxD) {
		if (_maxD >= minD && _maxD <= 127) {
			maxD = _maxD;
		} else {
			PApplet.println("ERROR in setMaxDepth(): level not in range. Level not changed.");
		}
	}

	public OctXYZ getDimension() {
		return dimension;
	}

	public OctXYZ getOrigin() {
		return origin;
	}

	public OctXYZ getCenter() {
		return origin.add(dimension.scale(0.5f));
	}

	public int getMinDepth() {
		return minD;
	}

	public int getMaxDepth() {
		return maxD;
	}

	public HashSet<OctNode> getNodes() {
		return nodeList;
	}

	// *************************************************************************************
	// RESET THE TREE
	// *************************************************************************************

	/**
	 * Clears the octree. Sets minD @ 0 and maxD @ 8, origin and dimension to
	 * (0,0,0).
	 */
	public void clear() {
		nodeList.clear();
		dimension = new OctXYZ();
		origin = new OctXYZ();
		minD = 0;
		maxD = 8;
	}

	/**
	 * Clears only the list of nodes.
	 */
	public void clearNodes() {
		nodeList.clear();
	}


	// *************************************************************************************
	// ADD AND REMOVE POINTS
	// *************************************************************************************

	/**
	 * Adds the node that contains the point at the specified level.
	 */
	public void addPoint(OctXYZ _tempPoint,int _level) {
		OctRST coord = _tempPoint.toRST(this);
		nodeList.add(new OctNode(PApplet.floor(coord.r*(1<<_level)),PApplet.floor(coord.s*(1<<_level)),PApplet.floor(coord.t*(1<<_level)),_level));
	}

	// *************************************************************************************
	// ADD AND REMOVE NODES
	// *************************************************************************************

	/**
	 * Adds a node. If the node is already in the octree, does nothing. Accepts
	 * node that are out of bounds or have a level outside the min/max depth.
	 */
	public void addNode(OctNode _tempNode) {
		nodeList.add(_tempNode);
	}

	/**
	 * Removes a node. If the node is not in the octree, does nothing.
	 */
	public void deleteNode(OctNode _tempNode) {
		nodeList.remove(_tempNode);
	}

	/**
	 * Adds a node. If the node is already in the octree, prints error.
	 */
	public void addNodeVerbose(OctNode _tempNode) {
		if (_tempNode != null) {
			if (nodeList.contains(_tempNode)) {
				PApplet.println("ERROR in nodeAdd: node already in octree.");
			} else
				nodeList.add(_tempNode);
		}
	}

	/**
	 * Removes a node. If the node is not in the octree, prints error.
	 */
	public void deleteNodeVerbose(OctNode _tempNode) {
		if (_tempNode != null) {
			if (nodeList.contains(_tempNode)) {
				nodeList.remove(_tempNode);
			} else
				PApplet.println("ERROR in nodeRemove: node not in octree.");
		}
	}

	/**
	 * Adds a list of nodes.
	 */
	public void addNode(ArrayList<OctNode> _tempNodeList) {
		for (OctNode _tempNode : _tempNodeList) {
			this.addNode(_tempNode);
		}
	}

	/**
	 * Removes a list of nodes.
	 */
	public void deleteNode(ArrayList<OctNode> _tempNodeList) {
		for (OctNode _tempNode : _tempNodeList) {
			this.deleteNode(_tempNode);
		}
	}

	/**
	 * Adds a set of nodes.
	 */
	public void addNode(HashSet<OctNode> _tempNodeList) {
		for (OctNode _tempNode : _tempNodeList) {
			this.addNode(_tempNode);
		}
	}

	/**
	 * Removes a set of nodes.
	 */
	public void deleteNode(HashSet<OctNode> _tempNodeList) {
		for (OctNode _tempNode : _tempNodeList) {
			this.deleteNode(_tempNode);
		}
	}

	public void selectNode(OctXYZ start, OctXYZ end) {
		//selectedNodes.addAll( this.getRayNodes(start,end).keySet());
		OctNode closest = this.getClosestNode( start,  end);
		for (OctNode _tempNode : this.nodeList) {
			if (_tempNode.equals(closest))
			{_tempNode.isSelected=true;
			selectedNodes.add(_tempNode);
			//PApplet.println("one node selected");
			}
		}
	}	

	private HashMap<OctNode,Float> getRayNodes(OctXYZ start, OctXYZ end) {
		//PApplet.println("getraytnode");
		OctOctree tempTree = new OctOctree(p5,1000,1000,1000);
		tempTree.dimension=this.dimension;
		tempTree.origin=this.origin;
		tempTree.maxD=this.maxD;
		tempTree.minD=this.minD;
		
		HashMap<OctNode,Float> _tempNodes = new HashMap<OctNode, Float>();
		ArrayList<OctNode> listOfNodes = new ArrayList<OctNode>();
		ArrayList<OctNode> listOfNodes2 = new ArrayList<OctNode>();

		tempTree.algGenerate(this.minD);
		for (OctNode _tempNode : tempTree.nodeList) {
			//PApplet.println('1');
			OctRST point = _tempNode.intersectPoint(start.toRST(this), end.toRST(this),-999999999,999999999);
			if (point != null) {
				//PApplet.println("point found");
				_tempNodes.put(_tempNode, point.sub(start.toRST(this)).getNorm());
				listOfNodes.add(_tempNode);
				//selectedNodes.add(_tempNode);
			}
		}
		tempTree.clearNodes();
		listOfNodes2.addAll(listOfNodes);
		listOfNodes.clear();
		for (int i = 0; i<this.maxD-this.minD; i++) {
			for (OctNode _tempNode : listOfNodes2) {
				for (OctNode _tempNode2 : _tempNode.getChildren()) {

					OctRST point = _tempNode2.intersectPoint(start.toRST(this), end.toRST(this),-999999999,999999999);
					if (point != null) {
						_tempNodes.put(_tempNode2, point.sub(start.toRST(this)).getNorm());
						listOfNodes.add(_tempNode2);
						//selectedNodes.add(_tempNode2);
					}
				}
			}
			listOfNodes2.clear();
			listOfNodes2.addAll(listOfNodes);
			listOfNodes.clear();
		}
		return _tempNodes;
	}

	private OctNode getClosestNode(OctXYZ start, OctXYZ direction) {
		//PApplet.println("getclosestnode");
		OctNode _tempNode = new OctNode();
		float closest=999999999;
		//PApplet.println(start.x +"y" + start.y +"y" +start.z );
		//PApplet.println(direction.x +"y" + direction.y +"y" +direction.z );
		HashMap<OctNode,Float> _tempNodes = getRayNodes( start,  direction);
		//PApplet.println(_tempNodes);
		for (OctNode candidates : _tempNodes.keySet()) {
			if (this.nodeList.contains(candidates) && closest > _tempNodes.get(candidates))
			{
				closest = _tempNodes.get(candidates);
				_tempNode = candidates;
			}
		}
		//PApplet.println("selected"+_tempNode);
		return  _tempNode;
	}

	// *************************************************************************************
	// BOOLEAN AND LIST OPERATIONS
	// *************************************************************************************

	/**
	 * Boolean subtracts one node from the octree. Requires the min and max
	 * depth level to the octree to be set correctly.
	 */
	public void boolSub(OctNode subNode) {
		ArrayList<OctNode> listAdd = new ArrayList<OctNode>();
		ArrayList<OctNode> listDel = new ArrayList<OctNode>();
		if (nodeList.contains(subNode)) {
			nodeList.remove(subNode);
			return;
		}
		for (OctNode temp : subNode.getAllChildren(this.maxD - subNode.level)) {
			if (nodeList.contains(temp)) {
				listDel.add(temp);
			}
		}
		for (OctNode parent : subNode.getAllParents(subNode.level - this.minD)) {
			if (nodeList.contains(parent)) {
				listDel.add(parent);
				for (int i = subNode.level - parent.level; i >= 1; i--) {
					listAdd.addAll(subNode.getParent(i).getChildren(1));
					if (i > 1)
						listAdd.remove(subNode.getParent(i - 1));
				}
				listAdd.remove(subNode);
			}
		}
		nodeList.addAll(listAdd);
		nodeList.removeAll(listDel);
	}

	/**
	 * Boolean adds one node from the octree. Requires the min and max depth
	 * level to the octree to be set correctly.
	 */
	public void boolAdd(OctNode addNode) {
		ArrayList<OctNode> listAdd = new ArrayList<OctNode>();
		ArrayList<OctNode> listDel = new ArrayList<OctNode>();
		if (!nodeList.contains(addNode)) {
			listAdd.add(addNode);
		}
		for (OctNode temp : addNode.getAllChildren(this.maxD - addNode.level)) {
			if (nodeList.contains(temp)) {
				listDel.add(temp);
			}
		}
		for (OctNode temp : addNode.getAllParents(addNode.level - this.minD)) {
			if (nodeList.contains(temp)) {
				listDel.add(addNode);
			}
		}
		nodeList.addAll(listAdd);
		nodeList.removeAll(listDel);
	}

	// *************************************************************************************
	// NODE METHODS
	// *************************************************************************************

	/**
	 * Will return a node if there is a node in the octree with the same code
	 * and level as the _tempNode. Else returns null.
	 */
	public OctNode filterNode(OctNode _tempNode) {
		for (OctNode tempNode : nodeList) {
			if (nodeList.contains(_tempNode)) {
				return tempNode;
			}
		}
		return null;
	}

	/**
	 * Filters a list of node to keep only the nodes in the octree.
	 */
	public ArrayList<OctNode> filterNode(ArrayList<OctNode> _tempNodeList) {
		ArrayList<OctNode> t = new ArrayList<OctNode>();
		for (OctNode _tempNode : _tempNodeList) {
			if (nodeList.contains(_tempNode)) {
				t.add(_tempNode);
			}
		}
		return t;
	}

	/**
	 * Filters a set of node to keep only the nodes in the octree.
	 */
	public HashSet<OctNode> filterNode(HashSet<OctNode> _tempNodeList) {
		HashSet<OctNode> t = new HashSet<OctNode>();
		for (OctNode _tempNode : _tempNodeList) {
			if (nodeList.contains(_tempNode)) {
				t.add(_tempNode);
			}
		}
		return t;
	}

	// TODO filter to keep only this or that octant

	/**
	 * Subdivides a node by replacing it by its 8 children. Don't use it in a
	 * loop, use nodeSubdivide(ArrayList) or nodeSubdivide(HashSet).
	 */
	public void subdivideNode(OctNode _tempNode) {
		nodeList.remove(_tempNode);
		nodeList.addAll(_tempNode.getChildren(1));
	}

	/**
	 * Subdivides a list of nodes by replacing it by their 8 children.
	 */
	public void subdivideNode(ArrayList<OctNode> _tempList) {
		HashSet<OctNode> _tempNodeList = new HashSet<OctNode>(_tempList);
		HashSet<OctNode> listAdd = new HashSet<OctNode>();
		for (OctNode temp : _tempNodeList) {
			if (nodeList.contains(temp)) {
				nodeList.remove(temp);
				listAdd.addAll(temp.getChildren(1));
			}
		}
		nodeList.addAll(listAdd);
	}

	/**
	 * Subdivides a list of nodes by replacing it by their 8 children.
	 */
	public void subdivideNode(HashSet<OctNode> _tempNodeList) {
		HashSet<OctNode> listAdd = new HashSet<OctNode>();
		for (OctNode temp : _tempNodeList) {
			if (nodeList.contains(temp)) {
				nodeList.remove(temp);
				listAdd.addAll(temp.getChildren(1));
			}
		}
		nodeList.addAll(listAdd);
	}

	/**
	 * Replaces a node by its parent. Removes the siblings that are also in the
	 * parent, but not smaller nodes! Can create nested nodes!
	 */
	public void mergeNode(OctNode _tempNode) {
		nodeList.remove(_tempNode);
		nodeList.removeAll(_tempNode.getSiblings());
		nodeList.add(_tempNode.getParent(1));
	}

	/**
	 * Replaces list of nodes by their parent. Removes the siblings that are
	 * also in the parent, but not smaller nodes! Can create nested nodes!
	 */
	public void mergeNode(ArrayList<OctNode> _tempList) {
		HashSet<OctNode> _toAdd = new HashSet<OctNode>();
		HashSet<OctNode> _toRemove = new HashSet<OctNode>();
		for (OctNode temp : _tempList) {
			if (nodeList.contains(temp)) {
				_toRemove.add(temp);
				_toRemove.addAll(temp.getSiblings());
				_toAdd.add(temp.getParent(1));
			}
		}
		nodeList.removeAll(_toRemove);
		nodeList.addAll(_toAdd);
	}

	/**
	 * Replaces list of nodes by their parent. Removes the siblings that are
	 * also in the parent, but not smaller nodes! Can create nested nodes!
	 */
	public void mergeNode(HashSet<OctNode> _tempList) {
		HashSet<OctNode> _toAdd = new HashSet<OctNode>();
		HashSet<OctNode> _toRemove = new HashSet<OctNode>();
		for (OctNode temp : _tempList) {
			if (nodeList.contains(temp)) {
				_toRemove.add(temp);
				_toRemove.addAll(temp.getSiblings());
				_toAdd.add(temp.getParent(1));
			}
		}
		nodeList.removeAll(_toRemove);
		nodeList.addAll(_toAdd);
	}

	// *************************************************************************************
	// ALGORITHMS
	// *************************************************************************************

	/**
	 * Automatically set the minimum and maximum depth of the octree according
	 * to its content. Can be run before other algorithms to ensure maximum
	 * efficiency or a consistent result.
	 */
	public void algFixDepth() {
		int _minD = 127;
		int _maxD = 0;
		for (OctNode temp : nodeList) {
			if (_minD > temp.getLevel())
				_minD = temp.getLevel();
			if (_maxD < temp.getLevel())
				_maxD = temp.getLevel();
		}
		minD = _minD;
		maxD = _maxD;
		PApplet.println("RESULT of algDepth: min depth: " + minD
				+ ", max depth: " + maxD + ".");
	}

	/**
	 * Generates all possible nodes at the specified level and add them to the
	 * tree. Requires the min and max depth level to the octree to be set
	 * correctly.
	 */
	public void algGenerate(int _depth) {
		int depth = _depth;
		if (depth > maxD) {
			PApplet.println("ERROR in algGenerate: start depth " + depth
					+ " larger than max depth " + maxD
					+ ". Max depth used instead.");
			depth = maxD;
		}
		if (depth < minD) {
			PApplet.println("ERROR in algGenerate: start depth " + depth
					+ " smaller than min depth " + maxD
					+ ". Min depth used instead.");
			depth = minD;
		}
		for (int p = (int) (1 << depth); --p >= 0;) {
			for (int n = (int) (1 << depth); --n >= 0;) {
				for (int m = (int) (1 << depth); --m >= 0;) {
					OctNode tempNode = new OctNode(m, n, p, depth);
					this.addNode(tempNode);
				}
			}
		}
	}

	/**
	 * Generates all nodes where the center has a value smaller/equal/larger to
	 * the threshold. Requires the min and max depth level to the octree to be
	 * set correctly.
	 */
	public void algGenerateByCenter(OctFunction _f, float _threshold,
			int _depth, boolean _smaller, boolean _equal, boolean _larger) {
		float data;
		int depth = _depth;
		if (depth > maxD) {
			PApplet.println("ERROR in algGenerate: start depth " + depth
					+ " larger than max depth " + maxD
					+ ". Max depth used instead.");
			depth = maxD;
		}
		if (depth < minD) {
			PApplet.println("ERROR in algGenerate: start depth " + depth
					+ " smaller than min depth " + maxD
					+ ". Min depth used instead.");
			depth = minD;
		}
		for (int m = (int) (1 << depth); --m >= 0;) {
			for (int n = (int) (1 << depth); --n >= 0;) {
				for (int p = (int) (1 << depth); --p >= 0;) {
					OctNode tempNode = new OctNode(m, n, p, depth);
					data = _f.compute(tempNode.getCenter().toXYZ(this).x,
							tempNode.getCenter().toXYZ(this).y, tempNode
							.getCenter().toXYZ(this).z);
					if (_smaller && data < _threshold) {
						this.addNode((OctNode) tempNode);
					}
					if (_equal && data == _threshold) {
						this.addNode((OctNode) tempNode);
					}
					if (_larger && data > _threshold) {
						this.addNode((OctNode) tempNode);
					}
				}
			}
		}
	}

	// TODO generate nodes at the border of the octree to close the surface

	/**
	 * Generates all nodes that are below/crossing/above the threshold by
	 * looking at the 8 corners. Requires the min and max depth level to the
	 * octree to be set correctly. Use the same start and end to avoid
	 * recursion.
	 */
	public void algGenerateByCorners(OctFunction _f, float _threshold,
			int _start, int _end, boolean _below, boolean _equal, boolean _above) {
		int start = _start;
		int end = _end;
		if (start > end) {
			PApplet.println("ERROR in algGenerate: start depth " + start
					+ " larger than end depth " + end
					+ ". Values have been inverted.");
			int temp = end;
			end = start;
			start = temp;
		}
		if (start > maxD) {
			PApplet.println("ERROR in algGenerate: start depth " + start
					+ " larger than max depth " + maxD
					+ ". Max depth used instead.");
			start = maxD;
		}
		if (start < minD) {
			PApplet.println("ERROR in algGenerate: start depth " + start
					+ " smaller than min depth " + maxD
					+ ". Min depth used instead.");
			start = minD;
		}
		if (end > maxD) {
			PApplet.println("ERROR in algGenerate: end depth " + end
					+ " larger than max depth " + maxD
					+ ". Max depth used instead.");
			end = maxD;
		}
		if (end < minD) {
			PApplet.println("ERROR in algGenerate: end depth " + end
					+ " smaller than min depth " + maxD
					+ ". Min depth used instead.");
			end = minD;
		}

		HashSet<OctNode> toAdd = new HashSet<OctNode>();
		for (int m = (int) (1 << start); --m >= 0;) {
			for (int n = (int) (1 << start); --n >= 0;) {
				for (int p = (int) (1 << start); --p >= 0;) {
					OctNode tempNode = new OctNode(m, n, p, start);
					float data = 0;
					for (OCT_VERTEX v : OCT_VERTEX.getAll()) {
						if (_f.compute(tempNode.getVertex(v).toXYZ(this).x,
								tempNode.getVertex(v).toXYZ(this).y, tempNode
								.getVertex(v).toXYZ(this).z) > _threshold) {
							data++;
						}
					}
					if (_below && data == 0) {
						this.addNode(tempNode);
					}
					if ((data > 0 && data < 8)) {
						toAdd.add(tempNode);
					}
					if (_above && (data == 8)) {
						this.addNode(tempNode);
					}
				}
			}
		}

		while (start < end) {
			start += 1;
			HashSet<OctNode> toAdd2 = new HashSet<OctNode>();
			for (OctNode tempNode1 : toAdd) {
				toAdd2.addAll(tempNode1.getChildren());
			}
			toAdd = new HashSet<OctNode>();
			for (OctNode tempNode2 : toAdd2) {
				float data = 0;
				for (OCT_VERTEX v : OCT_VERTEX.getAll()) {
					if (_f.compute(tempNode2.getVertex(v).toXYZ(this).x,
							tempNode2.getVertex(v).toXYZ(this).y, tempNode2
							.getVertex(v).toXYZ(this).z) > _threshold)
						data++;
				}
				if (_below && data == 0) {
					this.addNode(tempNode2);
				}
				if ((data > 0 && data < 8)) {
					toAdd.add(tempNode2);
				}
				if (_above && (data == 8)) {
					this.addNode(tempNode2);
				}
			}
		}

		if (_equal) this.addNode(toAdd);
	}

	// TODO simplify except if crossing a certain function, so that the surface is always clean

	/**
	 * Simplifies the octree by maximizing the size of the nodes, grouping nodes
	 * in bigger nodes. Requires the min and max depth level to the octree to be
	 * set correctly.
	 * 
	 * @param _mindepth
	 *            minimum node size to keep, even if the 8 octants are full
	 */
	public void algSimplify(int _mindepth) {
		int depth = _mindepth;
		if (depth > maxD) {
			PApplet.println("ERROR in algSimplify: start depth " + depth
					+ " larger than max depth " + maxD
					+ ". Max depth used instead.");
			depth = maxD;
		}
		if (depth < minD) {
			PApplet.println("ERROR in algSimplify: start depth " + depth
					+ " smaller than min depth " + maxD
					+ ". Min depth used instead.");
			depth = minD;
		}
		boolean needToIterate = false;
		PApplet.println("ITERATING algSimplify:\t" + nodeList.size() + " nodes");
		HashSet<OctNode> tempToAdd = new HashSet<OctNode>();
		HashSet<OctNode> tempToDelete = new HashSet<OctNode>();
		for (OctNode tempNode : nodeList) {
			int count = 0;
			if (tempNode.getLevel() > depth) {
				for (OctNode tempSibling : tempNode.getSiblings()) {
					if (nodeList.contains(tempSibling)) {
						count++;
					}
				}
				if (count == 7) {
					OctNode tempParentNode = tempNode.getParent();
					tempToAdd.add(tempParentNode);
					tempToDelete.addAll(tempParentNode.getChildren(1));
					needToIterate = true;
				}
			}
		}
		this.deleteNode(tempToDelete);
		this.addNode(tempToAdd);
		if (needToIterate) {
			PApplet.println("ITERATING algSimplify:\t" + nodeList.size()
			+ " nodes");
			algSimplify(depth);
		}
	}

	/**
	 * Ensures that one node has no neighbor smaller that a certain size by
	 * subdividing nodes. Requires the min and max depth level to the octree to
	 * be set correctly.
	 * 
	 * @param _constraint
	 *            maxmimum level difference between two adjacent nodes
	 */

	public void algConstrain(int _constraint) {
		int constraint = _constraint;
		if (constraint < 0) {
			PApplet.println("ERROR in algConstraint: Constraint cannot be negative, 0 used instead.");
			constraint = 0;
		}
		if (constraint > maxD - minD) {
			PApplet.println("ERROR in algConstraint: Constraint too large, maxD - minD used instead.");
			constraint = maxD - minD;
		}
		boolean needToIterate = false;
		PApplet.println("ITERATING algConstrain:\t" + nodeList.size()
		+ " nodes");
		HashSet<OctNode> tempToAdd = new HashSet<OctNode>();
		HashSet<OctNode> tempToRemove = new HashSet<OctNode>();
		ArrayList<OCT_ENUM> enumlist = new ArrayList<OCT_ENUM>();
		enumlist.addAll(OCT_FACE.getAll());
		enumlist.addAll(OCT_EDGE.getAll());
		for (OctNode tempNode : nodeList) {
			for (OCT_ENUM e : enumlist) {
				ArrayList<OctNode> tempList = tempNode.getAllBiggerNbrs(
						tempNode.level - this.minD, e);
				tempList = this.filterNode(tempList);
				if (tempList.size() > 0) {
					for (OctNode tempNeighborNode : tempList) {
						if (tempNode.getLevel() - tempNeighborNode.getLevel() > constraint
								&& tempToRemove.contains(tempNeighborNode) == false) {
							tempToAdd.addAll(tempNeighborNode.getChildren());
							tempToRemove.add(tempNeighborNode);
							needToIterate = true;
						}
					}
				}
			}
		}
		nodeList.removeAll(tempToRemove);
		for (OctNode tempNode : tempToAdd) {
			nodeList.add(tempNode);
		}
		if (needToIterate) {
			PApplet.println("ITERATING algConstrain:\t" + nodeList.size()
			+ " nodes");
			algConstrain(constraint);
		}
	}

	/**
	 * Ensures that no node is nested in another node. Requires the min and max
	 * depth level to the octree to be set correctly.
	 */
	public void algCleanNested() {
		Iterator<OctNode> itr = nodeList.iterator();
		while (itr.hasNext()) {
			OctNode t = itr.next();
			boolean delete = false;
			for (OctNode p : t.getAllParents(t.level - this.minD)) {
				if (nodeList.contains(p))
					delete = true;
			}
			if (delete)
				itr.remove();
		}
	}

	/**
	 * Ensures that no node is outside of the bounds of the octree. Requires the
	 * min and max depth level to the octree to be set correctly.
	 */
	public void algCleanOutOfBounds() {
		Iterator<OctNode> itr = nodeList.iterator();
		while (itr.hasNext()) {
			OctNode t = itr.next();
			if (t.isCodeWB(this) != true)
				itr.remove();
		}
	}

	/**
	 * Ensures that no node is outside of the bounds of the octree. Requires the
	 * min and max depth level to the octree to be set correctly.
	 */
	public void algCleanWrongLevels() {
		Iterator<OctNode> itr = nodeList.iterator();
		while (itr.hasNext()) {
			OctNode t = itr.next();
			if (t.isLevelWB(this) != true)
				itr.remove();
		}
	}

	// *************************************************************************************
	// DRAW
	// *************************************************************************************

	public void drawAsCenters() {
		p5.beginShape(PConstants.POINTS);
		for (OctNode tempNode : nodeList) {
			PVector center = tempNode.getCenter().toXYZ(this).toPVector();
			p5.vertex(center.x, center.y, center.z);
		}
		p5.endShape();
	}

	public void drawAsVertices() {
		p5.beginShape(PConstants.POINTS);
		for (OctNode tempNode : nodeList) {
			for (OctRST t : tempNode.getVertices()) {
				p5.vertex(t.toXYZ(this).x, t.toXYZ(this).y, t.toXYZ(this).z);
			}
		}
		p5.endShape();
	}

	public void drawAsEdges() {
		for (OctNode tempNode : nodeList) {
			p5.beginShape(PConstants.LINES);
			for (OCT_EDGE e : OCT_EDGE.values()) {
				for (OctRST t : tempNode.getVertices(e.getVertices())) {
					p5.vertex(t.toXYZ(this).x, t.toXYZ(this).y, t.toXYZ(this).z);
				}
			}
			p5.endShape();
		}
	}

	public void drawAsFaces() {
		for (OctNode tempNode : nodeList) {
			p5.beginShape(PConstants.QUADS);
			for (OCT_FACE f : OCT_FACE.values()) {
				for (OctRST t : tempNode.getVertices(f.getVertices())) {
					p5.vertex(t.toXYZ(this).x, t.toXYZ(this).y, t.toXYZ(this).z);
				}
			}
			p5.endShape();
		}
		for (OctNode tempNode : selectedNodes) {
			p5.fill(255, 0, 0,100);
			p5.beginShape(PConstants.QUADS);
			for (OCT_FACE f : OCT_FACE.values()) {
				for (OctRST t : tempNode.getVertices(f.getVertices())) {
					p5.vertex(t.toXYZ(this).x, t.toXYZ(this).y, t.toXYZ(this).z);
				}
			}
			p5.endShape();
		}
	}

	/**
	 * Requires the min and max depth level to the octree to be set correctly.
	 * Does not work with nested nodes.
	 */
	public void drawAsFaceSkeleton() {
		ArrayList<OCT_ENUM> tempList = new ArrayList<OCT_ENUM>();
		tempList.add(OCT_FACE.L);
		tempList.add(OCT_FACE.B);
		tempList.add(OCT_FACE.D);
		for (OctNode tempNode : nodeList) {
			p5.beginShape(PConstants.LINES);
			for (OCT_ENUM e : tempList) {
				if (nodeList.contains(tempNode.getNbr(e))) {
					OctXYZ p1 = tempNode.getCenter().toXYZ(this);
					OctXYZ p2 = tempNode.getNbr(e).getCenter().toXYZ(this);
					p5.vertex(p1.x, p1.y, p1.z);
					p5.vertex(p2.x, p2.y, p2.z);
				}
			}
			for (OctNode tempN : tempNode.getAllBiggerNbrs(tempNode.level
					- this.minD, new ArrayList<OCT_ENUM>(OCT_FACE.getAll()))) {
				if (nodeList.contains(tempN)) {
					OctXYZ p1 = tempNode.getCenter().toXYZ(this);
					OctXYZ p2 = tempN.getCenter().toXYZ(this);
					p5.vertex(p1.x, p1.y, p1.z);
					p5.vertex(p2.x, p2.y, p2.z);
				}
			}
			p5.endShape();
		}
	}

	/**
	 * Requires the min and max depth level to the octree to be set correctly.
	 * Does not work with nested nodes.
	 */
	public void drawAsEdgeSkeleton() {
		ArrayList<OCT_ENUM> tempList = new ArrayList<OCT_ENUM>();
		tempList.add(OCT_EDGE.LB);
		tempList.add(OCT_EDGE.BD);
		tempList.add(OCT_EDGE.LD);
		for (OctNode tempNode : nodeList) {
			p5.beginShape(PConstants.LINES);
			for (OCT_ENUM e : tempList) {
				if (nodeList.contains(tempNode.getNbr(e))) {
					OctXYZ p1 = tempNode.getCenter().toXYZ(this);
					OctXYZ p2 = tempNode.getNbr(e).getCenter().toXYZ(this);
					p5.vertex(p1.x, p1.y, p1.z);
					p5.vertex(p2.x, p2.y, p2.z);
				}
			}
			for (OctNode tempN : tempNode.getAllBiggerNbrs(tempNode.level
					- this.minD, new ArrayList<OCT_ENUM>(OCT_EDGE.getAll()))) {
				if (nodeList.contains(tempN)) {
					OctXYZ p1 = tempNode.getCenter().toXYZ(this);
					OctXYZ p2 = tempN.getCenter().toXYZ(this);
					p5.vertex(p1.x, p1.y, p1.z);
					p5.vertex(p2.x, p2.y, p2.z);
				}
			}
			p5.endShape();
		}
	}

	public void drawBoundingBox() {
		p5.pushMatrix();
		OctXYZ tempV = origin.add(dimension.scale(0.5f));
		p5.translate(tempV.x, tempV.y, tempV.z);
		p5.box(dimension.x, dimension.y, dimension.z);
		p5.popMatrix();
	}

	public void drawOrigin() {
		p5.point(origin.x, origin.y, origin.z);
	}

	public void drawAxis() {
		p5.stroke(0, 255, 255);
		p5.line(origin.x, origin.y, origin.z, origin.x + (dimension.x / 10),
				origin.y, origin.z);
		p5.stroke(255, 0, 255);
		p5.line(origin.x, origin.y, origin.z, origin.x, origin.y
				+ (dimension.y / 10), origin.z);
		p5.stroke(255, 255, 0);
		p5.line(origin.x, origin.y, origin.z, origin.x, origin.y, origin.z
				+ (dimension.z / 10));
	}

} public class OctRST implements java.io.Serializable {

	public float r;
	public float s;
	public float t;
		
	public OctRST() {
		r = 0f;
		s = 0f;
		t = 0f;
	}
	
	public OctRST(float _r, float _s, float _t) {
		r = _r;
		s = _s;
		t = _t;
	}
	
	public OctRST(PVector _vertex) {		
		r = _vertex.x;
		s = _vertex.y;
		t = _vertex.z;
	}

	public OctRST(OctRST _vertex) {		
		r = _vertex.r;
		s = _vertex.s;
		t = _vertex.t;
	}
		
	public OctXYZ toXYZ(OctXYZ _origin, OctXYZ _scale) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = this.r * _scale.x + _origin.x;
		tempV.y = this.s * _scale.y + _origin.y;
		tempV.z = this.t * _scale.z + _origin.z;
		return tempV;
	}

	public OctXYZ toXYZ(OctOctree _octree) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = this.r * _octree.getDimension().x + _octree.getOrigin().x;
		tempV.y = this.s * _octree.getDimension().y + _octree.getOrigin().y;
		tempV.z = this.t * _octree.getDimension().z + _octree.getOrigin().z;
		return tempV;
	}
	
	public PVector toPVector() {		
		PVector tempV = new PVector();
		tempV.x = this.r;
		tempV.y = this.s;
		tempV.z = this.t;
		return tempV;
	}
	
	public void addSelf(float _r, float _s, float _t) {
		r = r+_r;
		s = s+_s;
		t = t+_t;		
	}
	
	public void subSelf(float _r, float _s, float _t) {
		r = r-_r;
		s = s-_s;
		t = t-_t;		
	}
	
	public void addSelf(OctRST _tempV) {
		r = r+_tempV.r;
		s = s+_tempV.s;
		t = t+_tempV.t;		
	}
	
	public void subSelf(OctRST _tempV) {
		r = r-_tempV.r;
		s = s-_tempV.s;
		t = t-_tempV.t;		
	}
	
	public void scaleSelf(float _f) {
		r = r*_f;
		s = s*_f;
		t = t*_f;		
	}
	
	public OctRST add(float _r, float _s, float _t) {
		OctRST tempV = new OctRST();
		tempV.r = r+_r;
		tempV.s = s+_s;
		tempV.t = t+_t;
		return tempV;
	}
	
	public OctRST sub(float _r, float _s, float _t) {
		OctRST tempV = new OctRST();
		tempV.r = r-_r;
		tempV.s = s-_s;
		tempV.t = t-_t;	
		return tempV;
	}
	
	public OctRST add(OctRST _tempV) {
		OctRST tempV = new OctRST();
		tempV.r = r+_tempV.r;
		tempV.s = s+_tempV.s;
		tempV.t = t+_tempV.t;
		return tempV;
	}
	
	public OctRST sub(OctRST _tempV) {
		OctRST tempV = new OctRST();
		tempV.r = r-_tempV.r;
		tempV.s = s-_tempV.s;
		tempV.t = t-_tempV.t;	
		return tempV;
	}
	
	public OctRST scale(float _f) {
		OctRST tempV = new OctRST();
		tempV.r = r*_f;
		tempV.s = s*_f;
		tempV.t = t*_f;
		return tempV;
	}
	
	/**
	 * Is the RST point in the node? Requires an octree.
	 */
	public boolean isInNode(OctNode _n) {
		boolean cont = true;
		if (!(this.r < _n.getVertex(OCT_VERTEX.get(1)).r && this.r > _n.getVertex(OCT_VERTEX.get(0)).r))
			cont = false;
		if (!(this.s < _n.getVertex(OCT_VERTEX.get(3)).s && this.s > _n.getVertex(OCT_VERTEX.get(0)).s))
			cont = false;
		if (!(this.t < _n.getVertex(OCT_VERTEX.get(4)).t && this.t > _n.getVertex(OCT_VERTEX.get(0)).t))
			cont = false;
		return cont;
	}
	 
    public float dot(OctRST vec){
        return (this.r * vec.r +
        		this.s * vec.s +
                this.t * vec.t);
    }

    public OctRST cross(OctRST vec){
        float newR = this.s*vec.t - this.t*vec.s;
        float newS = this.t*vec.r - this.r*vec.t;
        float newT = this.r*vec.s - this.s*vec.r;
        return new OctRST(newR, newS, newT);
    }
    
    public float getNorm() {
        return PApplet.sqrt (r * r + s * s + t * t);
    }
    
    public OctRST normalize(){
        return new OctRST(r/this.getNorm(), s/this.getNorm(), t/this.getNorm());
    }
	
	public void draw(OctOctree _octree) {
		_octree.p5.beginShape(PConstants.POINTS);
		OctRST t = new OctRST(this);
		_octree.p5.vertex(t.toXYZ(_octree).x, t.toXYZ(_octree).y,
				t.toXYZ(_octree).z);
		_octree.p5.endShape();
	}
	
	// TODO is in node
	
	@Override	
	public int hashCode() {		
		return new String("r" + r + "s" + s + "t" + t).hashCode();
	}

	@Override
	public boolean equals(Object _tempV) {
		OctRST tempV = (OctRST) _tempV;
		if (tempV.r == this.r && tempV.s == this.s && tempV.t == this.t)
			return true;
		else
			return false;
	}
	
	public void invertSelf() {
		r = -r;
		s = -s;
		t = -t;	
	}

} public class OctXYZ implements java.io.Serializable {

	public float x;
	public float y;
	public float z;
		
	public OctXYZ() {
		x = 0f;
		y = 0f;
		z = 0f;
	}
	
	public OctXYZ(float _x, float _y, float _z) {
		x = _x;
		y = _y;
		z = _z;
	}

	public OctXYZ(PVector _vertex) {		
		x = _vertex.x;
		y = _vertex.y;
		z = _vertex.z;
	}
	
	public OctXYZ(OctXYZ _vertex) {		
		x = _vertex.x;
		y = _vertex.y;
		z = _vertex.z;
	}
		
	public OctRST toRST(OctXYZ _origin, OctXYZ _scale) {
		OctRST tempV = new OctRST();
		tempV.r = (this.x- _origin.x) / _scale.x ;
		tempV.s = (this.y- _origin.y) / _scale.y ;
		tempV.t = (this.z- _origin.z) / _scale.z ;
		return tempV;
	}

	public OctRST toRST(OctOctree _octree) {
		OctRST tempV = new OctRST();
		tempV.r = (this.x - _octree.getOrigin().x) / _octree.getDimension().x ;
		tempV.s = (this.y - _octree.getOrigin().y) / _octree.getDimension().y ;
		tempV.t = (this.z - _octree.getOrigin().z) / _octree.getDimension().z ;
		return tempV;
	}
	
	//TODO
	public OctRST toRST(OctOctree _octree, int level) {
		OctRST tempV = new OctRST();
		tempV.r = (this.x - _octree.getOrigin().x) / _octree.getDimension().x * PApplet.pow(2,level);
		tempV.s = (this.y - _octree.getOrigin().y) / _octree.getDimension().y * PApplet.pow(2,level) ;
		tempV.t = (this.z - _octree.getOrigin().z) / _octree.getDimension().z * PApplet.pow(2,level) ;
		return tempV;
	}
	
	public PVector toPVector() {		
		PVector tempV = new PVector();
		tempV.x = this.x;
		tempV.y = this.y;
		tempV.z = this.z;
		return tempV;
	}
	
	public void addSelf(float _x, float _y, float _z) {
		x = x+_x;
		y = y+_y;
		z = z+_z;		
	}
	
	public void subSelf(float _x, float _y, float _z) {
		x = x-_x;
		y = y-_y;
		z = z-_z;		
	}
	
	public void addSelf(OctXYZ _tempV) {
		x = x+_tempV.x;
		y = y+_tempV.y;
		z = z+_tempV.z;		
	}
	
	public void subSelf(OctXYZ _tempV) {
		x = x-_tempV.x;
		y = y-_tempV.y;
		z = z-_tempV.z;		
	}
	
	public void scaleSelf(float _f) {
		x = x*_f;
		y = y*_f;
		z = z*_f;		
	}
	
	public OctXYZ add(float _x, float _y, float _z) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = x+_x;
		tempV.y = y+_y;
		tempV.z = z+_z;
		return tempV;
	}
	
	public OctXYZ sub(float _x, float _y, float _z) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = x-_x;
		tempV.y = y-_y;
		tempV.z = z-_z;	
		return tempV;
	}
	
	public OctXYZ add(OctXYZ _tempV) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = x+_tempV.x;
		tempV.y = y+_tempV.y;
		tempV.z = z+_tempV.z;
		return tempV;
	}
	
	public OctXYZ sub(OctXYZ _tempV) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = x-_tempV.x;
		tempV.y = y-_tempV.y;
		tempV.z = z-_tempV.z;	
		return tempV;
	}
	
	public OctXYZ scale(float _f) {
		OctXYZ tempV = new OctXYZ();
		tempV.x = x*_f;
		tempV.y = y*_f;
		tempV.z = z*_f;
		return tempV;
	}
	 
    public float dot(OctXYZ vec){
        return (this.x * vec.x +
        		this.y * vec.y +
                this.z * vec.z);
    }

    public OctXYZ cross(OctXYZ vec){
        float newX = this.y*vec.z - this.z*vec.y;
        float newY = this.z*vec.x - this.x*vec.z;
        float newZ = this.x*vec.y - this.y*vec.x;
        return new OctXYZ(newX, newY, newZ);
    }

    public float getNorm() {
        return PApplet.sqrt (x * x + y * y + z * z);
    }
    
    public OctXYZ normalize(){
        return new OctXYZ(x/this.getNorm(), y/this.getNorm(), z/this.getNorm());
    }
	
	/**
	 * Is the XYZ point in the node _n of the octree _o?
	 * Requires an octree to convert from XYZ to RST coordinates.
	 */
	public boolean isInNode(OctNode _n, OctOctree _o) {
		boolean cont = true;
		if (!(this.x <= _n.getVertex(OCT_VERTEX.get(1)).toXYZ(_o).x && this.x >= _n.getVertex(OCT_VERTEX.get(0)).toXYZ(_o).x))
			cont = false;
		if (!(this.y <= _n.getVertex(OCT_VERTEX.get(3)).toXYZ(_o).y && this.y >= _n.getVertex(OCT_VERTEX.get(0)).toXYZ(_o).y))
			cont = false;
		if (!(this.z <= _n.getVertex(OCT_VERTEX.get(4)).toXYZ(_o).z && this.z >= _n.getVertex(OCT_VERTEX.get(0)).toXYZ(_o).z))
			cont = false;
		return cont;
	}
	
	public void draw(OctOctree _octree) {
		_octree.p5.beginShape(PConstants.POINTS);
		OctXYZ t = new OctXYZ(this);
		_octree.p5.vertex(t.x, t.y, t.z);
		_octree.p5.endShape();
	}
	
	@Override	
	public int hashCode() {		
		return new String("x" + x + "y" + y + "z" + z).hashCode();
	}

	@Override
	public boolean equals(Object _tempV) {
		OctXYZ tempV = (OctXYZ) _tempV;
		if (tempV.x == this.x && tempV.y == this.y && tempV.z == this.z)
			return true;
		else
			return false;
	}

	public void invertSelf() {
		x = -x;
		y = -y;
		z = -z;	
	}

}
import java.util.ArrayList;
import java.util.HashMap;
import oct.enums.OCT_EDGE;
import oct.enums.OCT_ENUM;
import oct.enums.OCT_FACE;
import oct.enums.OCT_VERTEX;
import oct.math.OctFunction;
import oct.math.Jama.Matrix;
import oct.math.Jama.SingularValueDecomposition;
import oct.octree.OctNode;
import oct.octree.OctOctree;
import oct.octree.OctRST;
import oct.octree.OctXYZ;
import processing.core.PApplet;

				if (valueAt[e.getV0().getOrdinal()] * flip > threshold * flip) {
					bool1 = true;
				}
				if (valueAt[e.getV1().getOrdinal()] * flip > threshold * flip) {
					bool2 = true;
				}
				if (bool1 != bool2) {
					adapt = PApplet.map(threshold, valueAt[e.getV0()
							.getOrdinal()], valueAt[e.getV1().getOrdinal()],
							-1, 1);

					// 1 position inside the node
					OctRST vX = new OctRST((e.r) * nodeSizeX / 2, (e.s)
							* nodeSizeY / 2, (e.t) * nodeSizeZ / 2);

					// 2 smooth on the edge
					if (e.r == 0)
						vX.r += adapt * nodeSizeX / 2;
					if (e.s == 0)
						vX.s += adapt * nodeSizeY / 2;
					if (e.t == 0)
						vX.t += adapt * nodeSizeZ / 2;

					// 3 position inside the octree
					vX.addSelf((float) (posx + 0.5) * nodeSizeX,
							(float) (posy + 0.5) * nodeSizeY,
							(float) (posz + 0.5) * nodeSizeZ);

					OctXYZ g = getGradient(vX.toXYZ(myOctree));

					A.set(e.getOrdinal(), 0, g.x);
					A.set(e.getOrdinal(), 1, g.y);
					A.set(e.getOrdinal(), 2, g.z);
					// B.set(e.getOrdinal(), 0, g.dot(vX.toXYZ(myOctree).sub(masspoint)));
					vertexList[e.getOrdinal()] = vX.toXYZ(myOctree);
					gList[e.getOrdinal()] = g;
				}
			}

			masspoint = new OctXYZ();
			float count = 0;
			for (int a = 0; a < 12; a++) {
				if (vertexList[a] != null) {
					masspoint.addSelf(vertexList[a]);
					count++;
				}
			}
			masspoint.scaleSelf(1f / count);
			for (int a = 0; a < 12; a++) {
				if (vertexList[a] != null) {
					B.set(a, 0, gList[a].dot(vertexList[a].sub(masspoint)));
				}
			}

			// from
			// https://github.com/jarrettchisholm/glr/blob/master/src/terrain/dual_contouring/Qef.cpp

			SingularValueDecomposition mySvd;
			Matrix At = new Matrix(3, 12, 0);
			At = A.transpose();
			Matrix AtA = new Matrix(12, 12, 0);
			AtA = A.times(At);
			mySvd = new SingularValueDecomposition(AtA);
			OctXYZ vector = new OctXYZ((float) mySvd.getSingularValues()[0],
					(float) mySvd.getSingularValues()[1],
					(float) mySvd.getSingularValues()[2]);

			Matrix UU = new Matrix(12, 3);
			Matrix VV = new Matrix(3, 3);
			Matrix DD = new Matrix(3, 3);
			Matrix svdsolve = new Matrix(3, 1);

			UU = A.svd().getU();
			VV = A.svd().getV();
			DD = A.svd().getS();

			float param = qefClip;

			if (DD.get(2, 2) < param) {
				DD.set(2, 2, 0.0f);
			}
			if (DD.get(1, 1) < param) {
				DD.set(1, 1, 0.0f);
			}
			if (DD.get(0, 0) < param) {
				DD.set(0, 0, 0.0f);
			}

			int rows = 12;
			double w[] = new double[3];
			for (int ii = 0; ii < rows; ii++) {
				if (B.get(ii, 0) != 0.0)
					for (int jj = 0; jj < 3; jj++)
						w[jj] += B.get(ii, 0) * UU.get(ii, jj);

			}

			// // // introduce non-zero singular values in d into w
			for (int ii = 0; ii < 3; ii++) {
				if (DD.get(ii, ii) != 0.0)
					w[ii] /= DD.get(ii, ii);
			}

			// // // compute result vector x = V * w
			for (int ii = 0; ii < 3; ii++) {
				double tmp = 0.0;
				for (int jj = 0; jj < 3; jj++) {
					tmp += (w[jj] * VV.get(ii, jj));
				}
				svdsolve.set(ii, 0, tmp);
			}

			vector = new OctXYZ((float) svdsolve.get(0, 0), (float) svdsolve.get(1, 0), (float) svdsolve.get(2, 0));

			vector = (vector.scale(adaptMultiplier)).add(masspoint);
			
			qef.put(new OctNode(tempNode),vector);
			
		}


		ArrayList<OCT_ENUM> dirList = new ArrayList<OCT_ENUM>();
