[FIX] Now it should compile somewhere other than here too

src/modules/lic/fibernavigator/Vector.cpp

+// Vector.cpp: implementation of the Vector class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "Vector.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Vector::Vector()
+{
+	x = 0;
+	y = 0;
+	z = 0;
+}
+
+Vector::Vector(double newX, double newY, double newZ)
+{
+	x = newX;
+	y = newY;
+	z = newZ;
+}
+
+Vector::Vector(double array[])
+{
+	x = array[0];
+	y = array[1];
+	z = array[2];
+}
+
+Vector::Vector(float array[])
+{
+	x = array[0];
+	y = array[1];
+	z = array[2];
+}
+
+Vector::Vector(int array[])
+{
+	x = array[0];
+	y = array[1];
+	z = array[2];
+}
+
+Vector::~Vector()
+{
+
+}
+
+Vector Vector::operator+(const Vector& rhs) const
+{
+    return Vector(x + rhs.x, y + rhs.y, z + rhs.z);
+}
+
+Vector Vector::operator-(const Vector& rhs) const
+{
+    return Vector(x - rhs.x, y - rhs.y, z - rhs.z);
+}
+
+Vector Vector::operator*(const double factor) const
+{
+    return Vector(x * factor, y * factor, z * factor);
+}
+
+Vector Vector::operator/(const double factor) const
+{
+    return Vector(x / factor, y / factor, z / factor);
+}
+
+double Vector::operator[](const int index) const
+{
+	switch (index) {
+		case VECTOR_X:
+			return x;
+			break;
+		case VECTOR_Y:
+			return y;
+			break;
+		default:
+			return z;
+			break;
+	}
+}
+
+double& Vector::operator[](const int index)
+{
+	switch (index) {
+		case VECTOR_X:
+			return x;
+			break;
+		case VECTOR_Y:
+			return y;
+			break;
+		default:
+			return z;
+			break;
+	}
+}
+
+void Vector::translateBy(const Vector& rhs)
+{
+	x += rhs.x;
+	y += rhs.y;
+	z += rhs.z;
+}
+
+void Vector::scaleBy(const double factor)
+{
+	x *= factor;
+	y *= factor;
+	z *= factor;
+}
+
+void Vector::normalize()
+{
+	double length = sqrt(x * x + y * y + z * z);
+	if (length > 0.000001) {
+		x /= length;
+		y /= length;
+		z /= length;
+	}
+}
+
+double Vector::normalizeAndReturn()
+{
+	double length = sqrt(x * x + y * y + z * z);
+	if (length > 0.0001) {
+		x /= length;
+		y /= length;
+		z /= length;
+	}
+	return length;
+}
+
+void Vector::zero()
+{
+	x = y = z = 0.0;
+}
+
+double Vector::getLength() const
+{
+	return sqrt(x * x + y * y + z * z);
+}
+
+double Vector::getSquaredLength() const
+{
+	return (x * x + y * y + z * z);
+}
+
+double Vector::Dot(const Vector& rhs) const
+{
+	return (x * rhs.x + y * rhs.y + z * rhs.z);
+}
+
+Vector Vector::Cross(const Vector& rhs) const
+{
+	return Vector((y * rhs.z) - (z * rhs.y),
+		          (z * rhs.x) - (x * rhs.z),
+				  (x * rhs.y) - (y * rhs.x));
+}
+
+void Vector::toArray(double array[]) const
+{
+	array[0] = x;
+	array[1] = y;
+	array[2] = z;
+}
+
+void Vector::toArray(float array[]) const
+{
+	array[0] = (float)x;
+	array[1] = (float)y;
+	array[2] = (float)z;
+}
+
+void Vector::fromArray(double array[])
+{
+	x = array[0];
+	y = array[1];
+	z = array[2];
+}
+
+void Vector::fromArray(float array[])
+{
+	x = array[0];
+	y = array[1];
+	z = array[2];
+}
+
+void Vector::rotateX(const double degrees) {
+	double radians = degrees * 3.14159265358979 / 180.0;
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origY = y;
+	y =	y * cosAngle - z * sinAngle;
+	z = origY * sinAngle + z * cosAngle;
+}
+
+void Vector::rotateY(const double degrees) {
+	double radians = degrees * 3.14159265358979 / 180.0;
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origX = x;
+	x =	x * cosAngle + z * sinAngle;
+	z = z * cosAngle - origX * sinAngle;
+}
+
+void Vector::rotateZ(const double degrees) {
+	double radians = degrees * 3.14159265358979 / 180.0;
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origX = x;
+	x =	x * cosAngle - y * sinAngle;
+	y = origX * sinAngle + y * cosAngle;
+}
+
+void Vector::radianRotateX(const double radians) {
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origY = y;
+	y =	y * cosAngle - z * sinAngle;
+	z = origY * sinAngle + z * cosAngle;
+}
+
+void Vector::radianRotateY(const double radians) {
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origX = x;
+	x =	x * cosAngle + z * sinAngle;
+	z = z * cosAngle - origX * sinAngle;
+}
+
+void Vector::radianRotateZ(const double radians) {
+	double cosAngle = cos(radians);
+	double sinAngle = sin(radians);
+	double origX = x;
+	x =	x * cosAngle - y * sinAngle;
+	y = origX * sinAngle + y * cosAngle;
+}
+
+void Vector::rotateAxis(const Vector& axis, double degrees)
+{
+	Vector rotator;
+	Vector rx;
+	Vector ry;
+	Vector rz;
+
+	// Normalize the axis and store in rz for later
+	// potential use.
+	rz = axis / axis.getLength();
+
+	// If this vector is parallel to the axis (Now in rz),
+	// don't bother rotating it. Check with Cauchy-Schwartz
+	// u dot v == length of u * length of v iff u is linearly
+	// dependant on v. Length of rz is one, because it is
+	// normalized.
+	if (!(fabs(this->Dot(rz)) == this->getLength())) {
+	    // If we're not already rotating around the Z, transform to Z
+	    if (axis.x == 0 && axis.y == 0) {
+            // In this case, the axis is along Z already, so we
+			// wont bother rotating our axis of rotation to Z.
+			// Since we're checking directly for 0, this is going
+			// to happen extremely rarely. However, if it does
+			// happen, some of the below math falls apart. So,
+			// we check.
+			rotator = *this;
+		} else {
+	        // Build the rotation matrix
+	        // rz was assigned already while normalizing the axis.
+	        rx = this->Cross(axis);
+	        rx.normalize();
+	        ry = rz.Cross(rx);
+
+	        // Move this vector such that the axis would be in Z
+	        rotator = Vector(rx.Dot(*this), ry.Dot(*this), rz.Dot(*this));
+		}
+
+	    // Rotate this vector around Z
+		rotator.rotateZ(degrees);
+
+	    if (axis.x == 0 && axis.y == 0) {
+		    *this = rotator;
+		} else {
+    	    // Move back so axis is in original location.
+	        this->x = rotator.x * rx.x +
+		              rotator.y * ry.x +
+			          rotator.z * rz.x;
+	        this->y = rotator.x * rx.y +
+		              rotator.y * ry.y +
+			          rotator.z * rz.y;
+	        this->z = rotator.x * rx.z +
+		              rotator.y * ry.z +
+		    	      rotator.z * rz.z;
+		}
+	}
+}
+
+Vector Vector::interpolate1(const Vector& endPoint, const double t) const
+{
+	return Vector((1.0 - t) * x + t * endPoint.x,
+		          (1.0 - t) * y + t * endPoint.y,
+				  (1.0 - t) * z + t * endPoint.z);
+}
+
+Vector Vector::interpolate2(const Vector& midControl, const Vector& endControl, const double t) const
+{
+    Vector left = this->interpolate1(midControl, t);
+	Vector right = midControl.interpolate1(endControl, t);
+	return left.interpolate1(right, t);
+}
+
+Vector Vector::interpolate3(const Vector& leftControl, const Vector& rightControl, const Vector& endControl, const double t) const
+{
+    Vector begin = this->interpolate1(leftControl, t);
+	Vector mid = leftControl.interpolate1(rightControl, t);
+	Vector end = rightControl.interpolate1(endControl, t);
+	return begin.interpolate2(mid, end, t);
+}

src/modules/lic/fibernavigator/Vector.h

+// Vector.h: interface for the Vector class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#if !defined(AFX_VECTOR_H__0AFDDCFD_6A7A_4943_8068_CBA5CEC834C5__INCLUDED_)
+#define AFX_VECTOR_H__0AFDDCFD_6A7A_4943_8068_CBA5CEC834C5__INCLUDED_
+
+#if _MSC_VER > 1000
+#pragma once
+#endif // _MSC_VER > 1000
+
+#include <math.h>
+
+#define VECTOR_X 0
+#define VECTOR_Y 1
+#define VECTOR_Z 2
+
+class Vector  
+{
+public:
+	Vector();
+	Vector(double x, double y, double z);
+	Vector(double[]);
+	Vector(float[]);
+	Vector(int[]);
+	virtual ~Vector();
+
+	double x;
+	double y;
+	double z;
+
+	/* Operators for vector addition. Note that these
+	   create a new Vector.
+	*/
+	Vector operator+(const Vector& rhs) const;
+	Vector operator-(const Vector& rhs) const;
+
+	/* Operators for scalar multiplication. Note that
+	   these create a new Vector
+	*/
+	Vector operator*(const double factor) const;
+	Vector operator/(const double factor) const;
+
+	/* Array operator allows access to x, y and z
+	   members through indices VECTOR_X, VECTOR_Y,
+	   and VECTOR_Z, respectively (0, 1, and 2).
+	   Only these three values should be used. If
+	   another index is used, the Z-value is returned.
+	   (RHS version)
+	*/
+	double operator[](const int index) const;
+
+	/* Array operator allows access to x, y and z
+	   members through indices VECTOR_X, VECTOR_Y,
+	   and VECTOR_Z, respectively (0, 1, and 2).
+	   Only these three values should be used. If
+	   another index is used, the Z-value is returned.
+	   (LHS version)
+	*/
+	double& operator[](const int index);
+
+	/* Operation for vector addition. This operation
+	   changes the original Vector. If the original is
+	   not needed after the operation, then this is faster
+	   than creating a new Vector object.
+	*/
+    void translateBy(const Vector& rhs);
+
+	/* Operation for scalar multiplication. This operation
+	   changes the original Vector. If the original is
+	   not needed after the operation, then this is faster
+	   than creating a new Vector object.
+	*/
+	void scaleBy(const double factor);
+
+	/* Normalizes the vector, if the vector is not the
+	   zero vector. If it is the zero vector, then it is
+	   unchanged.
+	*/
+	void normalize();
+
+	/* Normalizes the vector, if the vector is not the
+	   zero vector. If it is the zero vector, then it is
+	   unchanged. Returns the length the vector has before
+	   normalization.
+	*/
+	double normalizeAndReturn();
+
+	/* Zeros out the vector.
+	*/
+	void zero();
+
+	/* Return the length of the vector
+	*/
+	double getLength() const;
+	/* Return the squared length of the vector. Slightly
+	   faster, and we may only need to compare relative lengths.
+	*/
+	double getSquaredLength() const;
+
+	/* Return the dot product of this vector with the
+	   specified right-hand side vector.
+	*/
+	double Dot(const Vector& rhs) const;
+
+	/* Return the cross product of this vector with the
+	   specified right-hand side vector.
+    */
+	Vector Cross(const Vector& rhs) const;
+
+	/* Put this vector in the specified float array. The
+	   array must be at least 3 elements long.
+	*/
+	void toArray(double array[]) const;
+	/* Put this vector in the specified float array. The
+	   array must be at least 3 elements long. Each component
+	   is cast to a float, so loss of precision is likely.
+	*/
+	void toArray(float array[]) const;
+
+	/* Set this vector to the contents of the first three
+	   elements of the specified double array. The
+	   array must be at least 3 elements long.
+	*/
+	void fromArray(double array[]);
+	/* Set this vector to the contents of the first three
+	   elements of the specified float array. The
+	   array must be at least 3 elements long.
+	*/
+	void fromArray(float array[]);
+
+	/* Rotate this vector about an axis. The rotation is specified
+	   in degrees.
+	*/
+	void rotateX(const double degrees);
+	void rotateY(const double degrees);
+	void rotateZ(const double degrees);
+
+	/* Rotate this vector about an axis. The rotation is specified
+	   in radians. This is faster than the rotations specified in degrees, if
+	   you already know the radians.
+	*/
+	void radianRotateX(const double radians);
+	void radianRotateY(const double radians);
+	void radianRotateZ(const double radians);
+
+	/* Rotate this vector about an arbitrary axis. The rotation
+	   is specified in degrees.
+	*/
+	void rotateAxis(const Vector& axis, const double degrees);
+
+	/* Return a linearly interpolated vector between this
+	   vector and an endpoint. t should vary between 0 and 1.
+	*/
+	Vector interpolate1(const Vector& endPoint, const double t) const;
+
+	/* Return a qudratic Bezier interpolated vector with the
+	   three controls points this vector, midControl, and endControl.
+	   t should vary between 0 and 1.
+	*/
+	Vector interpolate2(const Vector& midControl, const Vector& endControl, const double t) const;
+
+	/* Return a cubic Bezier interpolated vector with the four
+	   controls points this vector, leftControl, rightControl,
+	   and endControl. t should vary between 0 and 1.
+	*/
+	Vector interpolate3(const Vector& leftControl, const Vector& rightControl, const Vector& endControl, const double t) const;
+
+};
+
+#endif // !defined(AFX_VECTOR_H__0AFDDCFD_6A7A_4943_8068_CBA5CEC834C5__INCLUDED_)

src/modules/lic/fibernavigator/loopSubD.cpp

+// loopSubD.cpp: implementation of the loopSubD class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "loopSubD.h"
+#include "triangleMesh.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+loopSubD::loopSubD()
+{
+
+}
+
+loopSubD::~loopSubD()
+{
+
+}
+
+loopSubD::loopSubD(TriangleMesh* nTriMesh){
+
+	triMesh = nTriMesh;
+	Vector newEdge;
+
+	numTriVerts = triMesh->getNumVertices();
+	numTriFaces = triMesh->getNumTriangles();
+
+	triMesh->reserveVerts(numTriVerts*4);
+	triMesh->reserveTriangles(numTriFaces*4);
+
+	//wxDateTime dt = wxDateTime::Now();
+	//printf("[%02d:%02d:%02d] ", dt.GetHour(), dt.GetMinute(), dt.GetSecond());
+
+	Vector* newVertexPositions = new Vector[numTriVerts];
+	//printf("start loop subdivision on %d vertices , %d triangles\n", numTriVerts, numTriFaces);
+	for(int i=0; i<numTriVerts; i++){
+		newVertexPositions[i] = calcNewPosition(i);
+	}
+
+	//printf("loop subdivision pass 2\n");
+	for(int i=0; i<numTriFaces; i++){
+		insertCenterTriangle(i);
+	}
+
+	//printf("loop subdivision pass 3\n");
+	for(int i=0; i<numTriFaces; i++){
+		insertCornerTriangles(i);
+	}
+
+	//printf("loop subdivision pass 4\n");
+	for(int i=0; i<numTriVerts; i++){
+		triMesh->setVertex(i, newVertexPositions[i]);
+	}
+
+	delete[] newVertexPositions;
+
+	//wxDateTime dt1 = wxDateTime::Now();
+	//printf("[%02d:%02d:%02d] ", dt1.GetHour(), dt1.GetMinute(), dt1.GetSecond());
+	//printf("loop subdivision done, mesh now has %d vertices , %d triangles\n",triMesh->getNumVertices(),triMesh->getNumTriangles());
+	//triMesh->printInfo();
+
+}
+
+Vector loopSubD::calcNewPosition(unsigned int vertNum)
+{
+	std::vector<unsigned int> starP = triMesh->getStar(vertNum);
+	int starSize = starP.size();
+
+
+	Vector oldPos = triMesh->getVertex(vertNum);
+	double alpha = getAlpha(starSize);
+	oldPos.scaleBy(1.0 - ((double) starSize * alpha));
+
+	Vector newPos;
+	int edgeV = 0;
+	for(int i=0; i<starSize; i++){
+		edgeV = triMesh->getNextVertex(starP[i], vertNum);
+		newPos.translateBy(triMesh->getVertex(edgeV));
+	}
+	newPos.scaleBy(alpha);
+
+	return oldPos + newPos;
+}
+
+void loopSubD::insertCenterTriangle(unsigned int triNum){
+
+	Triangle intP = triMesh->getTriangle(triNum);
+	int edgeVerts[3];
+
+	for(int i=0; i<3; i++){
+		edgeVerts[i] = calcEdgeVert(triNum, intP.pointID[i], intP.pointID[(i + 1) % 3], intP.pointID[(i + 2) % 3]);
+	}
+	triMesh->addTriangle(edgeVerts[0], edgeVerts[1], edgeVerts[2], triMesh->getTriangleTensor(triNum));
+}
+
+void loopSubD::insertCornerTriangles(unsigned int triNum){
+
+	// comment:		center are twisted from the orignal vertices.
+	// original:	0, 1, 2
+	// center:		a, b, c
+	// reAsgnOrig:	0, a, c
+	// addTris:		1, b, a
+	// addTris:		2, c, b
+	//
+	Triangle originalTri = triMesh->getTriangle(triNum);
+	Triangle centerTri   = triMesh->getTriangle(triNum + numTriFaces);
+
+	triMesh->addTriangle(originalTri.pointID[1], centerTri.pointID[1], centerTri.pointID[0], triMesh->getTriangleTensor(triNum));
+	triMesh->addTriangle(originalTri.pointID[2], centerTri.pointID[2], centerTri.pointID[1], triMesh->getTriangleTensor(triNum));
+	triMesh->setTriangle(triNum, originalTri.pointID[0], centerTri.pointID[0], centerTri.pointID[2]);
+}
+
+int loopSubD::calcEdgeVert(int triNum, unsigned int edgeV1, unsigned int edgeV2, unsigned int V3){
+
+