outmoduled spherical harmonics, including differential

src/main.cpp

@@ -2,7 +2,11 @@
 #include "utils.h"
 #define TINYPLY_IMPLEMENTATION
 #include "ply.h"
-#include "sh/spherical_harmonics.h"
+#include "shbasis.h"
+
+#define float double
+#include "nelder-mead-optimizer/optimizer.h"
+#undef float
 
 using namespace std;
 using namespace Eigen;
@@ -25,8 +29,12 @@ inline Vec diffOriginDistance(const Vec& n){
 }
 
 inline Vec contactPoint(const Vec& n, double r){
-	Vec dh = diffOriginDistance(n);
-	return (originDistance(n,r) - dh.dot(n))*n + dh;
+
+	return originDistance(n,r)*n + diffOriginDistance(n);
+
+	// if n and diffOriginDistance(n) are not guaranteed to be perpendicular, use:
+	// Vec dh = diffOriginDistance(n);
+	// return (originDistance(n,r) - dh.dot(n))*n + dh;
 }
 
 // curve radius of an edge acc. to https://computergraphics.stackexchange.com/a/1719
@@ -37,53 +45,51 @@ inline double edgeRadius(const Vec& p1, const Vec& n1, const Vec& p2, const Vec&
 int main()
 {
 
-	Idx shDeg = 8; // degree of spherical harmonics
-	Idx basisSize = (shDeg+1)*(shDeg+1); // number of basis functions
-
-	auto mesh = read_ply_file("ico8.ply");
+	// load sphere mesh
+	auto mesh = read_ply_file("ico6.ply");
 	const Idx nVerts = mesh.v.size();
 
 	double r = 4*sqrt(6)/9;
 
+	// normalize data
 	for(auto& v:mesh.v){v.normalize();} // float precision -> double precision
 	mesh.n = mesh.v; // valid for sphere
 
 	for(Idx i=0; i<nVerts; i++){
 		auto& v = mesh.v[i];
 		auto& n = mesh.n[i];
 
-		v = contactPoint(n,r);
+		//v = contactPoint(n,r);
 	}
 
+	// test curve radius accuracy
 	double coat = 1e300;
 	for(auto& e:mesh.e){
 		double er = edgeRadius(mesh.v[e[0]], mesh.n[e[0]], mesh.v[e[1]], mesh.n[e[1]]);
 		if(er < coat){coat = er;}
 	}
 
+	// random
+	default_random_engine rnd;
+	uniform_real_distribution<double> dist(-1,1);
 
-	MatrixXd basis(basisSize, nVerts);
-	struct lm{int l; int m;};
-	vector<lm> lmPairs(basisSize);
-	{	
-		Idx i=0;
-		for(int l=0; l<=shDeg; l++){
-			for(int m=-l; m<=l; m++){
-				lmPairs[i] = {l, m};
-				i++;
-			}
-		}
-	}
+	// precompute spherical harmonics basis
+	auto start = std::chrono::high_resolution_clock::now();
+
+	ShBasis shb(15, mesh.n);
+
+	auto finish = std::chrono::high_resolution_clock::now();
+	std::chrono::duration<double> elapsed = finish - start;
+	cout << elapsed.count() << "\n";
+
+	for(Idx ilm=0; ilm<shb.basisSize; ilm++){
+		double r = dist(rnd)*0.01;
 
-	Idx u=0;
-	for(Idx ilm=0; ilm<basisSize; ilm++){
-		for(Idx iv=0; iv<nVerts; iv++){
-			basis(ilm, iv) = EvalSH(lmPairs[ilm].l, lmPairs[ilm].m, mesh.n[iv]);
-			u++;
+		for(Idx i=0; i<nVerts; i++){
+			mesh.v[i] += mesh.n[i] * shb.basis(ilm, i) * r;
 		}
 	}
 
-	cout << u << "\n";
 
 	write_ply_file("output.ply", mesh);
 #ifdef _DEBUG

src/nelder-mead-optimizer/optimizer.h

@@ -0,0 +1,248 @@
+// https://github.com/blinry/nelder-mead-optimizer
+
+#include <ctime>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <algorithm>
+using namespace std;
+
+// Float vector with standard operations
+class Vector {
+public:
+    Vector() {
+    }
+    Vector(float c0, float c1) {
+        coords.push_back(c0);
+        coords.push_back(c1);
+    }
+    Vector(float c0, float c1, float c2) {
+        coords.push_back(c0);
+        coords.push_back(c1);
+        coords.push_back(c2);
+    }
+
+    // add more constructors when N gets > 3
+
+    float& operator[](int i) {
+        return coords[i];
+    }
+    float at(int i) const {
+        return coords[i];
+    }
+    int dimension() const {
+        return coords.size();
+    }
+    void prepare(int size) {
+        for (int i=0; i<size; i++) {
+            coords.push_back(0);
+        }
+    }
+    Vector operator+(Vector other) {
+        Vector result;
+        result.prepare(dimension());
+        for (int i=0; i<dimension(); i++) {
+            result[i] = coords[i] + other[i];
+        }
+        return result;
+    }
+    void operator+=(Vector other) {
+        for (int i=0; i<dimension(); i++) {
+            coords[i] += other[i];
+        }
+    }
+    Vector operator-(Vector other) {
+        Vector result;
+        result.prepare(dimension());
+        for (int i=0; i<dimension(); i++) {
+            result[i] = coords[i] - other[i];
+        }
+        return result;
+    }
+    bool operator==(Vector other) {
+        if (dimension() != other.dimension()) {
+            return false;
+        }
+        for (int i=0; i<dimension(); i++) {
+            if (other[i] != coords[i]) {
+                return false;
+            }
+        }
+        return true;
+    }
+    Vector operator*(float factor) {
+        Vector result;
+        result.prepare(dimension());
+        for (int i=0; i<dimension(); i++) {
+            result[i] = coords[i]*factor;
+        }
+        return result;
+    }
+    Vector operator/(float factor) {
+        Vector result;
+        result.prepare(dimension());
+        for (int i=0; i<dimension(); i++) {
+            result[i] = coords[i]/factor;
+        }
+        return result;
+    }
+    void operator/=(float factor) {
+        for (int i=0; i<dimension(); i++) {
+            coords[i] /= factor;
+        }
+    }
+    bool operator<(const Vector other) const {
+        for (int i=0; i<dimension(); i++) {
+            if (at(i) < other.at(i))
+                return false;
+            else if (at(i) > other.at(i))
+                return true;
+        }
+        return false;
+    }
+    float length() {
+        float sum = 0;
+        for (int i=0; i<dimension(); i++) {
+            sum += coords[i]*coords[i];
+        }
+        return pow(sum, 0.5f);
+    }
+private:
+    vector<float> coords;
+};
+
+// This class stores known values for vectors. It throws unknown vectors.
+class ValueDB {
+    public:
+        ValueDB() {
+        }
+        float lookup(Vector vec) {
+            if (!contains(vec)) {
+                throw vec;
+            } else {
+                return values[vec];
+            }
+        }
+        void insert(Vector vec, float value) {
+            values[vec] = value;
+        }
+    private:
+        bool contains(Vector vec) {
+            map<Vector, float>::iterator it = values.find(vec); // TODO add tolerance
+            return it != values.end();
+        }
+        map<Vector, float> values;
+};
+
+class NelderMeadOptimizer {
+    public:
+        NelderMeadOptimizer(int dimension, float termination_distance=0.001) {
+            this->dimension = dimension;
+            srand(time(NULL));
+            alpha = 1;
+            gamma = 2;
+            rho = -0.5;
+            sigma = 0.5;
+            this->termination_distance = termination_distance;
+        }
+        // used in `step` to sort the vectors
+        bool operator()(const Vector& a, const Vector& b) {
+            return db.lookup(a) < db.lookup(b);
+        }
+        // termination criteria: each pair of vectors in the simplex has to
+        // have a distance of at most `termination_distance`
+        bool done() {
+            if (vectors.size() < dimension) {
+                return false;
+            }
+            for (int i=0; i<dimension+1; i++) {
+                for (int j=0; j<dimension+1; j++) {
+                    if (i==j) continue;
+                    if ((vectors[i]-vectors[j]).length() > termination_distance) {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+        void insert(Vector vec) {
+            if (vectors.size() < dimension+1) {
+                vectors.push_back(vec);
+            }
+        }
+        Vector step(Vector vec, float score) {
+            db.insert(vec, score);
+            try {
+                if (vectors.size() < dimension+1) {
+                    vectors.push_back(vec);
+                }
+
+                // otherwise: optimize!
+                if (vectors.size() == dimension+1) {
+                    while(!done()) {
+                        sort(vectors.begin(), vectors.end(), *this);
+                        Vector cog; // center of gravity
+                        cog.prepare(dimension);
+                        for (int i = 1; i<=dimension; i++) {
+                            cog += vectors[i];
+                        }
+                        cog /= dimension;
+                        Vector best = vectors[dimension];
+                        Vector worst = vectors[0];
+                        Vector second_worst = vectors[1];
+                        // reflect
+                        Vector reflected = cog + (cog - worst)*alpha;
+                        if (f(reflected) > f(second_worst) && f(reflected) < f(best)) {
+                            vectors[0] = reflected;
+                        } else if (f(reflected) > f(best)) {
+                            // expand
+                            Vector expanded = cog + (cog - worst)*gamma;
+                            if (f(expanded) > f(reflected)) {
+                                vectors[0] = expanded;
+                            } else {
+                                vectors[0] = reflected;
+                            }
+                        } else {
+                            // contract
+                            Vector contracted = cog + (cog - worst)*rho;
+                            if (f(contracted) > f(worst)) {
+                                vectors[0] = contracted;
+                            } else {
+                                for (int i=0; i<dimension; i++) {
+                                    vectors[i] = best + (vectors[i] - best)*sigma;
+                                }
+                            }
+                        }
+                    }
+
+                    // algorithm is terminating, output: simplex' center of gravity
+                    Vector cog;
+                    for (int i = 0; i<=dimension; i++) {
+                        cog += vectors[i];
+                    }
+                    return cog/(dimension+1);
+                } else {
+                    // as long as we don't have enough vectors, request random ones,
+                    // with coordinates between 0 and 1. If you want other start vectors,
+                    // simply ignore these and use `step` on the vectors you want.
+                    Vector result;
+                    result.prepare(dimension);
+                    for (int i = 0; i<dimension; ++i) {
+                        result[i] = 0.001*(rand()%1000);
+                    }
+                    return result;
+                }
+            } catch (Vector v) {
+                return v;
+            }
+        }
+    private:
+        float f(Vector vec) {
+            return db.lookup(vec);
+        }
+        int dimension;
+        float alpha, gamma, rho, sigma;
+        float termination_distance;
+        vector<Vector> vectors;
+        ValueDB db;
+};