computing first constant-width-solid using spherical harmonics

mqnc · Sep 23, 2018 · 35b260a · 35b260a
1 parent 83bad5c
commit 35b260a
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 17 deletions.
diff --git a/src/main.cpp b/src/main.cpp
@@ -17,41 +17,51 @@ int main()
 {
 
 	// load sphere mesh
-	auto mesh = read_ply_file("ico6.ply");
+	auto mesh = read_ply_file("ico8.ply");
 	const Idx nVerts = mesh.v.size();
 
 	// normalize data
 	for(auto& v:mesh.v){v.normalize();} // float precision -> double precision
 	mesh.n = mesh.v; // valid for sphere
 
-	// 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);
 
 	// precompute spherical harmonics basis
 	auto start = std::chrono::high_resolution_clock::now();
 
-	ShBasis shb(7, mesh.n, ShBasis::BOTH);
+	ShBasis shb(30, mesh.n, ShBasis::ODD);
 
 	auto finish = std::chrono::high_resolution_clock::now();
 	std::chrono::duration<double> elapsed = finish - start;
 	cout << elapsed.count() << "\n";
 
-	for(Idx ih=0; ih<shb.nHarmonics; ih++){
-		double r = dist(rnd)*0.1;
+	vector<double> coeffs(shb.nHarmonics);
+	for(auto& c:coeffs){
+		c = dist(rnd) * 0.01;
+	}
 
-		for(Idx i=0; i<nVerts; i++){
-			mesh.v[i] += mesh.n[i] * shb.basis(ih, i) * r;
+	for(Idx i=0; i<nVerts; i++){
+		mesh.v[i] = mesh.n[i]*100;
+		for(Idx ih=0; ih<shb.nHarmonics; ih++){
+			mesh.v[i] += mesh.n[i] * shb.basis(ih, i) * coeffs[ih];
+			mesh.v[i] += shb.diffBasis(ih, i) * coeffs[ih];
 		}
 	}
 
+	// find smallest curvature radius
+	double rCurvMin = 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 < rCurvMin){rCurvMin = er;}
+	}
+
+	// shrink surface by that radius
+	for(Idx i=0; i<nVerts; i++){
+		mesh.v[i] -= mesh.n[i]*rCurvMin;
+	}
+
 	write_ply_file("output.ply", mesh);
 #ifdef _DEBUG
 	system("PAUSE");

diff --git a/src/shbasis.h b/src/shbasis.h
@@ -8,11 +8,17 @@
 // using MatrixXV = Eigen::Matrix<Vec, Eigen::Dynamic, Eigen::Dynamic>;
 
 Vec perp1(Vec v){
-	return (Vec(v.y(), v.z(), -v.x()).cross(v)).normalized();
+	if(fabs(v.x()) <= fabs(v.y()) && fabs(v.x()) <= fabs(v.z())){
+		return v.cross(Vec(1,0,0)).normalized();
+	}
+	if(fabs(v.y()) <= fabs(v.z())){
+		return v.cross(Vec(0,1,0)).normalized();
+	}
+	return v.cross(Vec(0,0,1)).normalized();
 }
 
 Vec perp2(Vec v){
-	return v.cross(perp1(v));
+	return v.cross(perp1(v)).normalized();
 }
 
 class ShBasis{
@@ -60,8 +66,8 @@ class ShBasis{
 				Vec t2 = perp2(dir);
 
 				double b0  = sh::EvalSH(lmPairs[ih].l, lmPairs[ih].m, dir.normalized());
-				double db1 = sh::EvalSH(lmPairs[ih].l, lmPairs[ih].m, (dir+eps*t1).normalized());
-				double db2 = sh::EvalSH(lmPairs[ih].l, lmPairs[ih].m, (dir+eps*t2).normalized());
+				double db1 = sh::EvalSH(lmPairs[ih].l, lmPairs[ih].m, (dir+eps*t1).normalized()) - b0;
+				double db2 = sh::EvalSH(lmPairs[ih].l, lmPairs[ih].m, (dir+eps*t2).normalized()) - b0;
 
 				basis(ih, iDir) = b0;
 				diffBasis(ih, iDir) = db1/eps*t1 + db2/eps*t2;