diff --git a/.gitignore b/.gitignore index a5e7200..307cf2a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *~ #*# __pycache__/ -.ipynb_checkpoints/* +.ipynb_checkpoints/ +.ccls-cache/ diff --git a/python_isosurfaces_2018_2019/build_osx.sh b/python_isosurfaces_2018_2019/build_osx.sh new file mode 100755 index 0000000..6346360 --- /dev/null +++ b/python_isosurfaces_2018_2019/build_osx.sh @@ -0,0 +1,10 @@ +#!/bin/sh + +g++ -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include \ + -L/usr/local/Cellar/gmp/6.1.2_2/lib \ + -L/usr/local/Cellar/mpfr/4.0.1/lib \ + -L/usr/local/Cellar/boost/1.67.0_1/lib \ + -DCGAL_CONCURRENT_MESH_3 \ + -lCGAL -lgmp -lmpfr -lboost_thread-mt \ + ./mesh_an_implicit_function.cpp \ + -o mesh_an_implicit_function.o diff --git a/python_isosurfaces_2018_2019/cgal_dabbling.org b/python_isosurfaces_2018_2019/cgal_dabbling.org new file mode 100755 index 0000000..a26ecb5 --- /dev/null +++ b/python_isosurfaces_2018_2019/cgal_dabbling.org @@ -0,0 +1,173 @@ +#+TITLE: CGAL dabbling +#+DATE: <2018-08-06 Mon> +#+AUTHOR: Hodapp +#+EMAIL: hodapp87@gmail.com +#+OPTIONS: ':nil *:t -:t ::t <:t H:3 \n:nil ^:t arch:headline +#+OPTIONS: author:t c:nil creator:comment d:(not "LOGBOOK") date:t +#+OPTIONS: e:t email:nil f:t inline:t num:t p:nil pri:nil stat:t +#+OPTIONS: tags:t tasks:t tex:t timestamp:t toc:t todo:t |:t +#+DESCRIPTION: +#+EXCLUDE_TAGS: noexport +#+KEYWORDS: +#+LANGUAGE: en +#+SELECT_TAGS: export + +# By default I do not want that source code blocks are evaluated on export. Usually +# I want to evaluate them interactively and retain the original results. +#+PROPERTY: header-args :eval never-export :export both + +- CGAL is one of the most insanely cryptic and impenetrable libraries + I have found. +- Where I am stuck now: + - I can use 1D features in + [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] but this is the + wrong sort of data (it's a 3D mesh, yes, but not a surface mesh) + and "medit" and "vtu" are all it can write. How do I extract a + surface mesh to get a Polyhedron_3 so that I can write an OBJ? + - Or: Is there any way to use 1D features with a surface mesh? + - Even if I manage to put this Mesh_3 into the right form I have no + guarantees that it is a good surface mesh. It uses a different + algorithm than for the surface mesh - I'm simply trying to take + the surface of the result. I should also expect this algorithm to + be something more like O(N^3) with mesh resolution because it must + fill the volume with tetrahedrons, and I then just throw away all + of these. +- This is based around the source for [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_implicit_sphere_8cpp-example.html][mesh_implicit_sphere.cpp]] from + [[https://doc.cgal.org/latest/Mesh_3/index.html][CGAL: 3D Mesh Generation]]. +- Why am I not using [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh__3Functions.html#ga68ca38989087644fb6b893eb566be9ea][facets_in_complex_3_to_triangle_mesh()]]? +- [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] shows the use of + [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Mesh__domain__with__polyline__features__3.html][Mesh_domain_with_polyline_features_3]] to explicitly give features to + preserve + - This uses the ~make_mesh_3~ call while I'm using + ~make_surface_mesh~... + - I can only call ~CGAL::print_polyhedron_wavefront~ if I have a + ~Polyhedron~ and I cannot figure out how to get one from a ~C3t3~. + So, perhaps I am stuck with this "medit" mesh format, or else VTU. + - I think the terminology is trying to tell me: this isn't a surface + mesh, it's a 3D mesh made of tetrahedra. It sounds like the type + of data from a Delaunay triangulation is inherently rather + different from a surface mesh. + - https://doc.cgal.org/latest/Mesh_3/index.html#title24 does the + opposite direction of what I need +- There are [[https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__detect__features__grp.html][Feature Detection Functions]], but ~detect_features~ is in + [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__complex__mesh__domain__3.html][Polyhedral_complex_mesh_domain_3]] and + [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__mesh__domain__with__features__3.html#a5a868ac7b8673436766d28b7a80d2826][Polyhedral_mesh_domain_with_features_3]] +- Domains: + - [[https://doc.cgal.org/latest/Mesh_3/classMeshDomain__3.html][MeshDomain_3]] concept + - [[https://doc.cgal.org/latest/Mesh_3/classMeshDomainWithFeatures__3.html][MeshDomainWithFeatures_3]] concept +- Why make_surface_mesh, Implicit_surface_3, + Complex_2_in_triangulation_3 vs. make_mesh_3, + Implicit_mesh_domain_3, Mesh_complex_3_in_triangulation_3? +- [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh3Parameters.html#ga0a990b28d55157c62d4bfd2624d408af][parameters::features()]] - can do 1-dimensional features + - now how do I use it? +- [[https://doc.cgal.org/latest/Surface_mesher/index.html][make_surface_mesh]] args: + - ~SurfaceMeshC2T3& c2t3~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshComplex__2InTriangulation__3.html][SurfaceMeshComplex_2InTriangulation_3]] + - ~SurfaceMeshTraits::Surface_3 surface~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurface__3.html][Surface_3]] + - ~SurfaceMeshTraits traits~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshTraits__3.html][SurfaceMeshTraits_3]] + - ~FacetsCriteria criteria~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshFacetsCriteria__3.html][SurfaceMeshFacetsCriteria_3]] + - ~Tag~ +- "This algorithm of CGAL::make_surface_mesh is designed for smooth + implicit surfaces. If your implicit surface is not smooth, then the + sharp features of the surface will not be meshed correctly." + +#+BEGIN_SRC elisp + (setq org-confirm-babel-evaluate nil) + (setq org-src-fontify-natively t) + (setq org-src-tab-acts-natively t) + (org-version) +#+END_SRC + +#+RESULTS: +: 9.1.9 + +#+BEGIN_SRC sh :results verbatim + gcc --version +#+END_SRC + +#+RESULTS: +: Apple LLVM version 9.1.0 (clang-902.0.39.2) +: Target: x86_64-apple-darwin17.7.0 +: Thread model: posix +: InstalledDir: /Library/Developer/CommandLineTools/usr/bin + +#+NAME: includes +#+BEGIN_SRC C++ + #include + + #include + #include + #include + + #include + #include +#+END_SRC + +#+NAME: typesDomain +#+BEGIN_SRC C++ + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + typedef K::FT FT; + typedef K::Point_3 Point; + typedef FT (Function)(const Point&); + typedef CGAL::Implicit_mesh_domain_3 Mesh_domain; + + #ifdef CGAL_CONCURRENT_MESH_3 + typedef CGAL::Parallel_tag Concurrency_tag; + #else + typedef CGAL::Sequential_tag Concurrency_tag; + #endif +#+END_SRC + +#+NAME: typesTriangulation +#+BEGIN_SRC C++ + typedef CGAL::Mesh_triangulation_3::type Tr; + typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; +#+END_SRC + +#+NAME: typesCriteria +#+BEGIN_SRC C++ + typedef CGAL::Mesh_criteria_3 Mesh_criteria; +#+END_SRC + +#+NAME: sphereFunction +#+BEGIN_SRC C++ + FT sphere_function (const Point& p) { + return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; + } +#+END_SRC + +#+BEGIN_SRC C++ :noweb yes :flags -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include -L/usr/local/Cellar/gmp/6.1.2_2/lib -L/usr/local/Cellar/mpfr/4.0.1/lib -L/usr/local/Cellar/boost/1.67.0_1/lib -DCGAL_CONCURRENT_MESH_3 -lCGAL -lgmp -lmpfr -lboost_thread-mt + <> + + using namespace CGAL::parameters; + + <> + <> + <> + <> + + int main() + { + // Domain (Warning: Sphere_3 constructor uses squared radius !) + Mesh_domain domain(sphere_function, + K::Sphere_3(CGAL::ORIGIN, 2.)); + // Mesh criteria + Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025, + cell_radius_edge_ratio=2, cell_size=0.1); + + std::cout << "Generating..." << std::endl; + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + + // Output + std::ofstream medit_file("out.mesh"); + c3t3.output_to_medit(medit_file); + + std::cout << "Done" << std::endl; + + return 0; + } +#+END_SRC + +#+RESULTS: +| Generating... | +| Done | diff --git a/python_isosurfaces_2018_2019/isosurfaces.pdf b/python_isosurfaces_2018_2019/isosurfaces.pdf new file mode 100755 index 0000000..9aae5b3 Binary files /dev/null and b/python_isosurfaces_2018_2019/isosurfaces.pdf differ diff --git a/python_isosurfaces_2018_2019/isosurfaces.wxmx b/python_isosurfaces_2018_2019/isosurfaces.wxmx new file mode 100755 index 0000000..0ef3073 Binary files /dev/null and b/python_isosurfaces_2018_2019/isosurfaces.wxmx differ diff --git a/python_isosurfaces_2018_2019/mesh_an_implicit_function.cpp b/python_isosurfaces_2018_2019/mesh_an_implicit_function.cpp new file mode 100755 index 0000000..a9e0cab --- /dev/null +++ b/python_isosurfaces_2018_2019/mesh_an_implicit_function.cpp @@ -0,0 +1,171 @@ +// Taken from: +// https://doc.cgal.org/latest/Surface_mesher/Surface_mesher_2mesh_an_implicit_function_8cpp-example.html +// https://doc.cgal.org/latest/Surface_mesher/index.html +// https://doc.cgal.org/latest/Mesh_3/index.html + +#include +#include +#include +#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H + #include +#else + // NixOS currently has CGAL 4.11, not 4.12. I guess 4.12 is needed + // for this. The #ifdef is unnecessary, but the header and call for + // below are deprecated. + #include +#endif +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +// Triangulation +typedef CGAL::Surface_mesh_default_triangulation_3 Tr; +typedef CGAL::Complex_2_in_triangulation_3 C2t3; + +// Domain? +typedef Tr::Geom_traits GT; +typedef GT::Sphere_3 Sphere_3; +typedef GT::Point_3 Point_3; +typedef GT::Vector_3 Vector_3; +typedef GT::FT FT; +typedef FT (*Function)(Point_3); +typedef CGAL::Implicit_surface_3 Surface_3; +// how does this differ from CGAL::Implicit_mesh_domain_3? + +typedef CGAL::Polyhedron_3 Polyhedron; + +FT sphere_function(Point_3 p) { + Point_3 p2(p.x() + 0.1 * cos(p.x() * 20), + p.y(), + p.z() + 0.1 * sin(p.z() * 4)); + const FT x2=p2.x()*p2.x(), y2=p2.y()*p2.y(), z2=p2.z()*p2.z(); + return x2+y2+z2-1; +} + +Vector_3 sphere_gradient(Point_3 p) { + float A = 0.1; + float B = 0.1; + float F1 = 20; + float F2 = 4; + return Vector_3(2*(A*cos(p.x()*F1) + p.x())*(1 - A*F1*sin(p.x()*F1)), + 2*p.y(), + 2*(B*sin(p.z()*F2) + p.z())*(1 + B*F2*cos(p.z()*F2))); +} + +FT spiral_function(Point_3 p) { + float outer = 2.0; + float inner = 0.4; // 0.9 + float freq = 20; // 10 + float phase = M_PI; // 3 * M_PI / 2; + float thresh = 0.3; + const FT d1 = p.y()*outer - inner * sin(p.x()*freq + phase); + const FT d2 = p.z()*outer - inner * cos(p.x()*freq + phase); + return std::max(sqrt(d1*d1 + d2*d2) - thresh, + p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - 1.9*1.9); +} + +Vector_3 spiral_gradient(Point_3 p) { + float outer = 2.0; + float inner = 0.4; + float freq = 20; + float phase = M_PI; + float thresh = 0.3; + // "block([%1,%2,%3,%4,%5,%6],%1:P+x*F,%2:cos(%1),%3:z*O-I*%2,%4:sin(%1),%5:y*O-I*%4,%6:1/sqrt(%5^2+%3^2),[((2*F*I*%3*%4-2*F*I*%2*%5)*%6)/2,O*%5*%6,O*%3*%6])" + float v1 = phase + p.x() * freq; + float v2 = cos(v1); + float v3 = p.z() * outer - inner * v2; + float v4 = sin(v1); + float v5 = p.y() * outer - inner * v4; + float v6 = 1.0 / sqrt(v5*v5 + v3*v3); + return Vector_3(((2*freq*inner*v3*v4-2*freq*inner*v2*v5)*v6)/2, + outer * v5 * v6, + outer * v3 * v6); +} + +int main() { + Tr tr; // 3D-Delaunay triangulation + C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation + + FT bounding_sphere_rad = 2.0; + + // defining the surface + Surface_3 surface(spiral_function, // pointer to function + Sphere_3(CGAL::ORIGIN, + bounding_sphere_rad*bounding_sphere_rad)); // bounding sphere + + std::string fname("spiral_thing4.obj"); + float angular_bound = 30; + float radius_bound = 0.02; + float distance_bound = 0.02; + + // Note that "2." above is the *squared* radius of the bounding sphere! + // defining meshing criteria + CGAL::Surface_mesh_default_criteria_3 criteria( + angular_bound, radius_bound, distance_bound); + + std::cout << "angular bound = " << angular_bound << ", " + << "radius bound = " << radius_bound << ", " + << "distance bound = " << distance_bound << std::endl; + + std::cout << "Making surface mesh..." << std::endl; + // meshing surface + CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag()); + std::cout << "Vertices: " << tr.number_of_vertices() << std::endl; + + // This didn't work on some calls instead of 'poly': + //CGAL::Surface_mesh sm; + Polyhedron poly; + std::cout << "Turning facets to triangle mesh..." << std::endl; + +#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H + CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, poly); +#else + CGAL::output_surface_facets_to_polyhedron(c2t3, poly); +#endif + + FT err = 0.0; + FT inf = std::numeric_limits::infinity(); + for (Polyhedron::Point_iterator it = poly.points_begin(); + it != poly.points_end(); + ++it) + { + + FT rate = 2e-6; + FT f0 = abs(spiral_function(*it)); + FT f; + for (int i = 0; i < 100; ++i) { + f = spiral_function(*it); + Vector_3 grad(spiral_gradient(*it)); + + *it -= grad * rate * (f > 0 ? 1 : -1); + /* + std::cout << "Iter " << i << ": " + << "F(" << it->x() << "," << it->y() << "," << it->z() + << ")=" << f << ", F'=" << grad << std::endl; + */ + } + //FT diff = (abs(f) - abs(f0)) / f0; + /* + std::cout << "F(" << it->x() << "," << it->y() << "," << it->z() + << "): " << f0 << " to " << f << std::endl; + */ + + err += f * f; + } + err = sqrt(err); + std::cout << "RMS isosurface distance: " << err << std::endl; + + std::cout << "Mesh vertices: " << poly.size_of_vertices() << ", " + << "facets: " << poly.size_of_facets() << std::endl; + + std::cout << "Writing to " << fname << "..." << std::endl; + std::ofstream ofs(fname); + CGAL::print_polyhedron_wavefront(ofs, poly); +} diff --git a/python_isosurfaces_2018_2019/mesh_implicit_sphere.cpp b/python_isosurfaces_2018_2019/mesh_implicit_sphere.cpp new file mode 100755 index 0000000..6b0d49b --- /dev/null +++ b/python_isosurfaces_2018_2019/mesh_implicit_sphere.cpp @@ -0,0 +1,55 @@ +#include + +#include +#include +#include + +#include +#include + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::FT FT; +typedef K::Point_3 Point; +typedef FT (Function)(const Point&); +typedef CGAL::Implicit_mesh_domain_3 Mesh_domain; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +// Function +FT sphere_function (const Point& p) +{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; } + +int main() +{ + // Domain (Warning: Sphere_3 constructor uses squared radius !) + Mesh_domain domain(sphere_function, + K::Sphere_3(CGAL::ORIGIN, 2.)); + // Mesh criteria + Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025, + cell_radius_edge_ratio=2, cell_size=0.1); + + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + + // Output + std::ofstream medit_file("out.mesh"); + c3t3.output_to_medit(medit_file); + + return 0; +} + diff --git a/python_isosurfaces_2018_2019/mesh_two_implicit_spheres_with_balls.cpp b/python_isosurfaces_2018_2019/mesh_two_implicit_spheres_with_balls.cpp new file mode 100755 index 0000000..c45e7b7 --- /dev/null +++ b/python_isosurfaces_2018_2019/mesh_two_implicit_spheres_with_balls.cpp @@ -0,0 +1,106 @@ +#include + +#include +#include +#include + +#include +#include +#include + +#include + +// Kernel +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +// Domain +typedef K::FT FT; +typedef K::Point_3 Point; +typedef FT (Function)(const Point&); +typedef CGAL::Mesh_domain_with_polyline_features_3< + CGAL::Labeled_mesh_domain_3 > Mesh_domain; + +// Polyline +typedef std::vector Polyline_3; +typedef std::list Polylines; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3< + Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +// Function +FT sphere_function1 (const Point& p) +{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-2; } + +FT sphere_function2 (const Point& p) +{ return CGAL::squared_distance(p, Point(2, 0, 0))-2; } + +FT sphere_function (const Point& p) +{ + if(sphere_function1(p) < 0 || sphere_function2(p) < 0) + return -1; + else + return 1; +} + +#include + +int main() +{ + // Domain (Warning: Sphere_3 constructor uses squared radius !) + Mesh_domain domain = + Mesh_domain::create_implicit_mesh_domain(sphere_function, + K::Sphere_3(Point(1, 0, 0), 6.)); + + // Mesh criteria + Mesh_criteria criteria(edge_size = 0.15, + facet_angle = 25, facet_size = 0.15, + cell_radius_edge_ratio = 2, cell_size = 0.15); + + // Create edge that we want to preserve + Polylines polylines (1); + Polyline_3& polyline = polylines.front(); + + for(int i = 0; i < 360; ++i) + { + Point p (1, std::cos(i*CGAL_PI/180), std::sin(i*CGAL_PI/180)); + polyline.push_back(p); + } + polyline.push_back(polyline.front()); // close the line + + // Insert edge in domain + domain.add_features(polylines.begin(), polylines.end()); + + // Mesh generation without feature preservation + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + CGAL::parameters::no_features()); + + std::ofstream medit_file("out-no-protection.vtu"); + c3t3.output_to_medit(medit_file); + medit_file.close(); + c3t3.clear(); + + // Mesh generation with feature preservation + c3t3 = CGAL::make_mesh_3(domain, criteria); + + // Output + medit_file.open("out-with-protection.vtu"); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return 0; +} diff --git a/python_isosurfaces_2018_2019/nix-shell.sh b/python_isosurfaces_2018_2019/nix-shell.sh new file mode 100755 index 0000000..656111b --- /dev/null +++ b/python_isosurfaces_2018_2019/nix-shell.sh @@ -0,0 +1,3 @@ +#!/bin/sh + +nix-shell -p python37Packages.numpy-stl diff --git a/python_isosurfaces_2018_2019/scratch_ply.py b/python_isosurfaces_2018_2019/scratch_ply.py new file mode 100755 index 0000000..fd21679 --- /dev/null +++ b/python_isosurfaces_2018_2019/scratch_ply.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python + +import math + +def f1(T, I, O, P, F): + return lambda x: (x, (T + I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O) + +def f2(T, I, O, P, F): + return lambda x: (x, -(T - I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O) + +f = f2(O=2.0, I=0.4, F=20, P=math.pi, T=0.3) +print("ply") +print("format ascii 1.0") +r = range(-400, 400) +print("element vertex %d" % (len(r))) +print("property float32 x") +print("property float32 y") +print("property float32 z") +print("end_header") +for xi in r: + print("%f %f %f" % f(float(xi) / 200)) diff --git a/python_isosurfaces_2018_2019/shell.nix b/python_isosurfaces_2018_2019/shell.nix new file mode 100755 index 0000000..88752cb --- /dev/null +++ b/python_isosurfaces_2018_2019/shell.nix @@ -0,0 +1,9 @@ +{ pkgs ? import {} }: + +let python_with_deps = pkgs.python3.withPackages + (ps: [ps.numpy ps.numpy-stl]); +in pkgs.stdenv.mkDerivation rec { + name = "gfx_scratch"; + + buildInputs = with pkgs; [ python_with_deps ]; +} diff --git a/python_isosurfaces_2018_2019/shell_cgal.nix b/python_isosurfaces_2018_2019/shell_cgal.nix new file mode 100755 index 0000000..88ffd2e --- /dev/null +++ b/python_isosurfaces_2018_2019/shell_cgal.nix @@ -0,0 +1,10 @@ +{ pkgs ? import {} }: + +let stdenv = pkgs.stdenv; +in stdenv.mkDerivation rec { + name = "cgal_scratch"; + + buildInputs = with pkgs; [ cgal boost gmp mpfr ]; +} + +# g++ -lCGAL -lmpfr -lgmp mesh_an_implicit_function.cpp -o mesh_an_implicit_function.o diff --git a/python_isosurfaces_2018_2019/spiral_parametric.py b/python_isosurfaces_2018_2019/spiral_parametric.py new file mode 100755 index 0000000..520472b --- /dev/null +++ b/python_isosurfaces_2018_2019/spiral_parametric.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 + +import sys +import numpy +import stl.mesh + +# TODO: +# - This is a very naive triangulation strategy. It needs fixing - the +# way it handles 'flatter' areas isn't optimal at all, even if the +# sharper areas are much better than from CGAL or libfive. +# - Generate just part of the mesh and then copy. It is rotationally +# symmetric, as well as translationally symmetric at its period. + +fname = "spiral_outer0.stl" +freq = 20 +phase = 0 +scale = 1/16 # from libfive +inner = 0.4 * scale +outer = 2.0 * scale +rad = 0.3 * scale + +angle = lambda z: freq*z + phase + +# z starting & ending point: +z0 = -20*scale +z1 = 20*scale +# Number of z divisions: +m = 1600 +# Number of circle points: +n = 1000 + +dz = (z1 - z0) / (m-1) + +data = numpy.zeros((m-1)*n*2 + 2*n, dtype=stl.mesh.Mesh.dtype) +# Vertex count: +# From z0 to z0+dz is n circle points joined with 2 triangles to next -> n*2 +# z0+dz to z0+dz*2 is likewise... up through (m-1) of these -> (m-1)*n*2 +# Two endcaps each have circle points & center point -> 2*n +# Thus: (m-1)*n*2 + 2*n +v = data["vectors"] +print("Vertex count: {}".format(m*n*2 + 2*n)) + +verts = numpy.zeros((n, 3), dtype=numpy.float32) + +# For every z cross-section... +for z_idx in range(m): + #sys.stdout.write(".") + # z value: + z = z0 + dz*z_idx + # Angle of center point of circle (radians): + # (we don't actually need to normalize this) + rad = angle(z) + c,s = numpy.cos(rad), numpy.sin(rad) + # Center point of circle: + cx, cy = (inner + outer)*numpy.cos(rad), (inner + outer)*numpy.sin(rad) + # For every division of the circular cross-section... + if z_idx == 0: + # Begin with z0 endcap as a special case: + verts_last = numpy.zeros((n, 3), dtype=numpy.float32) + verts_last[:, 0] = cx + verts_last[:, 1] = cy + verts_last[:, 2] = z + else: + verts_last = verts + verts = numpy.zeros((n, 3), dtype=numpy.float32) + for ang_idx in range(n): + # Step around starting angle (the 'far' intersection of the + # line at angle 'rad' and this circle): + rad2 = rad + 2*ang_idx*numpy.pi/n + # ...and generate points on the circle: + xi = cx + outer*numpy.cos(rad2) + yi = cy + outer*numpy.sin(rad2) + verts[ang_idx, :] = [xi, yi, z] + #print("i={}, z={}, rad={}, cx={}, cy={}, rad2={}, xi={}, yi={}".format(i,z,rad,cx,cy, rad2, xi, yi)) + if z_idx == 0: + for i in range(n): + v[i][0,:] = verts[(i + 1) % n,:] + v[i][1,:] = verts[i,:] + v[i][2,:] = verts_last[i,:] + #print("Write vertex {}".format(i)) + else: + for i in range(n): + # Vertex index: + vi = z_idx*n*2 + i*2 - n + v[vi][0,:] = verts[(i + 1) % n,:] + v[vi][1,:] = verts[i,:] + v[vi][2,:] = verts_last[i,:] + #print("Write vertex {}".format(vi)) + v[vi+1][0,:] = verts_last[(i + 1) % n,:] + v[vi+1][1,:] = verts[(i + 1) % n,:] + v[vi+1][2,:] = verts_last[i,:] + #print("Write vertex {} (2nd half)".format(vi+1)) + +# then handle z1 endcap: +for i in range(n): + # See vi definition above. z_idx ends at m-1, i ends at n-1, and + # so evaluate vi+1 (final index it wrote), add 1 for the next, and + # then use 'i' to step one at a time: + vi = (m-1)*n*2 + (n-1)*2 - n + 2 + i + v[vi][0,:] = verts[i,:] + v[vi][1,:] = verts[(i + 1) % n,:] + v[vi][2,:] = [cx, cy, z] + # Note winding order (1 & 2 flipped from other endcap) + #print("Write vertex {} (endcap)".format(vi)) + +print("Writing {}...".format(fname)) +mesh = stl.mesh.Mesh(data, remove_empty_areas=False) +mesh.save(fname) +print("Done.") diff --git a/python_isosurfaces_2018_2019/spiral_parametric2.py b/python_isosurfaces_2018_2019/spiral_parametric2.py new file mode 100755 index 0000000..0794a0d --- /dev/null +++ b/python_isosurfaces_2018_2019/spiral_parametric2.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 + +import sys +import numpy +import stl.mesh + +# TODO: +# - Fix correction around high curvatures. It has boundary issues +# between the functions. +# - Check every vertex point against the actual isosurface. +# - Check rotation direction +# - Fix phase, which only works if 0! + +fname = "spiral_inner0_one_period.stl" +freq = 20 +phase = 0 +scale = 1/16 # from libfive +inner = 0.4 * scale +outer = 2.0 * scale +rad = 0.3 * scale +ext_phase = 0 + +""" +fname = "spiral_outer90_one_period.stl" +freq = 10 +#phase = numpy.pi/2 +phase = 0 +scale = 1/16 # from libfive +inner = 0.9 * scale +outer = 2.0 * scale +rad = 0.3 * scale +ext_phase = numpy.pi/2 +""" + +def angle(z): + return freq*z + phase + +def max_z(): + # This value is the largest |z| for which 'radical' >= 0 + # (thus, for x_cross to have a valid solution) + return (numpy.arcsin(rad / inner) - phase) / freq + +def radical(z): + return rad*rad - inner*inner * (numpy.sin(angle(z)))**2 + +def x_cross(z, sign): + # Single cross-section point in XZ for y=0. Set sign for positive + # or negative solution. + n1 = numpy.sqrt(radical(z)) + n2 = inner * numpy.cos(angle(z)) + if sign > 0: + return (n2-n1) / outer + else: + return (n2+n1) / outer + +def curvature_cross(z, sign): + # Curvature at a given cross-section point. This is fugly because + # it was produced from Maxima's optimized expression. + a1 = 1/outer + a2 = freq**2 + a3 = phase + z*freq + a4 = numpy.cos(a3) + a5 = a4**2 + a6 = numpy.sin(a3) + a7 = a6**2 + a8 = inner**2 + a9 = numpy.sqrt(rad**2 - a8*a7) + a10 = -a2*(inner**4)*a5*a7 / (a9**3) + a11 = 1 / a9 + a12 = -a2*a8*a5*a11 + a13 = a2*a8*a7*a11 + a14 = 1/(outer**2) + a15 = -freq*a8*a4*a6*a11 + if sign > 0: + return -a1*(a13+a12+a10+a2*inner*a4) / ((a14*(a15+freq*inner*a6)**2 + 1)**(3/2)) + else: + return a1*(a13+a13+a10-a2*inner*a4) / ((a14*(a15-freq*inner*a6)**2 + 1)**(3/2)) + +def cross_section_xz(eps): + # Generate points for a cross-section in XZ. 'eps' is the maximum + # distance in either axis. + verts = [] + signs = [-1, 1] + z_start = [0, max_z()] + z_end = [max_z(), 0] + # Yes, this is clunky and numerical: + for sign, z0, z1 in zip(signs, z_start, z_end): + print("sign={} z0={} z1={}".format(sign, z0, z1)) + z = z0 + x = x_cross(z, sign) + while (sign*z) >= (sign*z1): + verts.append([x, 0, z]) + x_last = x + dz = -sign*min(eps, abs(z - z1)) + if abs(dz) < 1e-8: + break + x = x_cross(z + dz, sign) + #curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign))) + curvature = abs(curvature_cross((z + dz)/2, sign)) + dx = (x - x_last) * curvature + print("start x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature)) + while abs(dx) > eps: + dz *= 0.8 + x = x_cross(z + dz, sign) + curvature = abs(curvature_cross((z + dz)/2, sign)) + #curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign))) + dx = (x - x_last) * curvature + print("iter x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature)) + z = z + dz + print("finish x={} z={} curvature={}".format(x, z, curvature)) + n = len(verts) + data = numpy.zeros((n*2, 3)) + data[:n, :] = verts + data[n:, :] = verts[::-1] + data[n:, 2] = -data[n:, 2] + return data + +def turn(points, dz): + # Note one full revolution is dz = 2*pi/freq + # How far to turn in radians (determined by dz): + rad = angle(dz) + c, s = numpy.cos(rad), numpy.sin(rad) + mtx = numpy.array([ + [ c, s, 0], + [-s, c, 0], + [ 0, 0, 1], + ]) + return points.dot(mtx) + [0, 0, dz] + +def screw_360(eps, dz): + #z0 = -10 * 2*numpy.pi/freq / 2 + z0 = -5 * 2*numpy.pi/freq / 2 + z1 = z0 + 2*numpy.pi/freq + #z1 = 5 * 2*numpy.pi/freq / 2 + #z0 = 0 + #z1 = 2*numpy.pi/freq + init_xsec = cross_section_xz(eps) + num_xsec_steps = init_xsec.shape[0] + zs = numpy.arange(z0, z1, dz) + num_screw_steps = len(zs) + vecs = num_xsec_steps * num_screw_steps * 2 + 2*num_xsec_steps + print("Generating {} vertices...".format(vecs)) + data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype) + v = data["vectors"] + # First endcap: + center = init_xsec.mean(0) + for i in range(num_xsec_steps): + v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:] + v[i][1,:] = init_xsec[i,:] + v[i][2,:] = center + # Body: + verts = init_xsec + for i,z in enumerate(zs): + verts_last = verts + verts = turn(init_xsec, z-z0) + if i > 0: + for j in range(num_xsec_steps): + # Vertex index: + vi = num_xsec_steps + (i-1)*num_xsec_steps*2 + j*2 + v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:] + v[vi][1,:] = verts[j,:] + v[vi][2,:] = verts_last[j,:] + #print("Write vertex {}".format(vi)) + v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:] + v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:] + v[vi+1][2,:] = verts_last[j,:] + #print("Write vertex {} (2nd half)".format(vi+1)) + # Second endcap: + center = verts.mean(0) + for i in range(num_xsec_steps): + vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i + v[vi][0,:] = center + v[vi][1,:] = verts[i,:] + v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:] + v[:, :, 2] += z0 + ext_phase / freq + v[:, :, :] /= scale + mesh = stl.mesh.Mesh(data, remove_empty_areas=False) + print("Beginning z: {}".format(z0/scale)) + print("Ending z: {}".format(z1/scale)) + print("Period: {}".format((z1-z0)/scale)) + return mesh + +#print("Writing {}...".format(fname)) +#mesh = stl.mesh.Mesh(data, remove_empty_areas=False) +#mesh.save(fname) +#print("Done.") + +# What's up with this? Note the jump from z=0.0424 to z=0.037. +""" +finish x=0.13228756555322954 z=0.042403103949074046 curvature=2.451108140319032 +sign=1 z0=0.042403103949074046 z1=0 +__main__:75: RuntimeWarning: invalid value encountered in double_scalars +start x=0.0834189730812818 dx=nan z=0.042403103949074046 dz=-0.005 curvature=nan +finish x=0.0834189730812818 z=0.03740310394907405 curvature=nan +""" +# Is it because curvature is undefined there - thus the starting step +# size of 0.005 is fine? + +m = screw_360(0.01, 0.001) +print("Writing {}...".format(fname)) +m.save(fname) diff --git a/python_isosurfaces_2018_2019/spiral_parametric3.py b/python_isosurfaces_2018_2019/spiral_parametric3.py new file mode 100755 index 0000000..804e757 --- /dev/null +++ b/python_isosurfaces_2018_2019/spiral_parametric3.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 + +import sys +import numpy +import stl.mesh + +# TODO: +# - This still has some strange errors around high curvature. +# They are plainly visible in the cross-section. +# (errr... which errors were those? I can't see in the render) +# - Check rotation direction +# - Fix phase, which only works if 0 (due to how I work with y) +# Things don't seem to line up right. +# - Why is there still a gap when using Array modifier? +# Check beginning and ending vertices maybe +# - Organize this so that it generates both meshes when run +# - Use SymPy instead of all this hard-coded stuff? + +# This is all rather tightly-coupled. Almost everything is specific +# to the isosurface I was trying to generate. walk_curve may be able +# to generalize to some other shapes. +class ExplicitSurfaceThing(object): + def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase): + self.freq = freq + self.phase = phase + self.scale = scale + self.inner = inner + self.outer = outer + self.rad = rad + self.ext_phase = ext_phase + + def angle(self, z): + return self.freq*z + self.phase + + def max_z(self): + # This value is the largest |z| for which 'radical' >= 0 + # (thus, for x_cross to have a valid solution) + return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq + + def radical(self, z): + return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2 + + # Implicit curve function + def F(self, x, z): + return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2 + + # Partial 1st derivatives of F: + def F_x(self, x, z): + return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z)) + + def F_z(self, x, z): + return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z)) + + # Curvature: + def K(self, x, z): + a1 = self.outer**2 + a2 = x**2 + a3 = self.freq*z + self.phase + a4 = numpy.cos(a3) + a5 = self.inner**2 + a6 = a4**2 + a7 = self.freq**2 + a8 = numpy.sin(a3)**2 + a9 = self.outer**3 + a10 = self.inner**3 + return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4)) + + def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7): + x, z = x0, z0 + eps2 = eps*eps + verts = [] + iters = 0 + # Until we return to the same point at which we started... + while True: + iters += 1 + verts.append([x, 0, z]) + # ...walk around the curve by stepping perpendicular to the + # gradient by 'eps'. So, first find the gradient: + dx = self.F_x(x, z) + dz = self.F_z(x, z) + # Normalize it: + f = 1/numpy.sqrt(dx*dx + dz*dz) + nx, nz = dx*f, dz*f + # Find curvature at this point because it tells us a little + # about how far we can safely move: + K_val = abs(self.K(x, z)) + eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps) + # Scale by 'eps' and use (-dz, dx) as perpendicular: + px, pz = -nz*eps_corr, nx*eps_corr + # Walk in that direction: + x += px + z += pz + # Moving in that direction is only good locally, and we may + # have deviated off the curve slightly. The implicit function + # tells us (sort of) how far away we are, and the gradient + # tells us how to minimize that: + #print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format( + # x, z, dx, dz, px, pz, K_val, eps_corr)) + F_val = self.F(x, z) + count = 0 + while abs(F_val) > gd_thresh: + count += 1 + dx = self.F_x(x, z) + dz = self.F_z(x, z) + f = 1/numpy.sqrt(dx*dx + dz*dz) + nx, nz = dx*f, dz*f + # If F is negative, we want to increase it (thus, follow + # gradient). If F is positive, we want to decrease it + # (thus, opposite of gradient). + F_val = self.F(x, z) + x += -F_val*nx + z += -F_val*nz + # Yes, this is inefficient gradient-descent... + diff = numpy.sqrt((x-x0)**2 + (z-z0)**2) + #print("{} gradient-descent iters. diff = {}".format(count, diff)) + if iters > 100 and diff < thresh: + #print("diff < eps, quitting") + #verts.append([x, 0, z]) + break + data = numpy.array(verts) + return data + + def x_cross(self, z, sign): + # Single cross-section point in XZ for y=0. Set sign for positive + # or negative solution. + n1 = numpy.sqrt(self.radical(z)) + n2 = self.inner * numpy.cos(self.angle(z)) + if sign > 0: + return (n2-n1) / self.outer + else: + return (n2+n1) / self.outer + + def turn(self, points, dz): + # Note one full revolution is dz = 2*pi/freq + # How far to turn in radians (determined by dz): + rad = self.angle(dz) + c, s = numpy.cos(rad), numpy.sin(rad) + mtx = numpy.array([ + [ c, s, 0], + [-s, c, 0], + [ 0, 0, 1], + ]) + return points.dot(mtx) + [0, 0, dz] + + def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False): + #z0 = -10 * 2*numpy.pi/freq / 2 + z0 = z0_period_start * 2*numpy.pi/self.freq / 2 + z1 = z0 + 2*numpy.pi/self.freq + #z1 = 5 * 2*numpy.pi/freq / 2 + #z0 = 0 + #z1 = 2*numpy.pi/freq + init_xsec = self.walk_curve(x_init, z_init, eps, thresh) + num_xsec_steps = init_xsec.shape[0] + zs = numpy.append(numpy.arange(z0, z1, dz), z1) + num_screw_steps = len(zs) + vecs = num_xsec_steps * num_screw_steps * 2 + offset = 0 + if endcaps: + offset = num_xsec_steps + vecs += 2*num_xsec_steps + print("Generating {} vertices...".format(vecs)) + data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype) + v = data["vectors"] + # First endcap: + if endcaps: + center = init_xsec.mean(0) + for i in range(num_xsec_steps): + v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:] + v[i][1,:] = init_xsec[i,:] + v[i][2,:] = center + # Body: + verts = init_xsec + for i,z in enumerate(zs): + verts_last = verts + verts = self.turn(init_xsec, z-z0) + if i > 0: + for j in range(num_xsec_steps): + # Vertex index: + vi = offset + (i-1)*num_xsec_steps*2 + j*2 + v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:] + v[vi][1,:] = verts[j,:] + v[vi][2,:] = verts_last[j,:] + #print("Write vertex {}".format(vi)) + v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:] + v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:] + v[vi+1][2,:] = verts_last[j,:] + #print("Write vertex {} (2nd half)".format(vi+1)) + # Second endcap: + if endcaps: + center = verts.mean(0) + for i in range(num_xsec_steps): + vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i + v[vi][0,:] = center + v[vi][1,:] = verts[i,:] + v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:] + v[:, :, 2] += z0 + self.ext_phase / self.freq + v[:, :, :] /= self.scale + mesh = stl.mesh.Mesh(data, remove_empty_areas=False) + print("Beginning z: {}".format(z0/self.scale)) + print("Ending z: {}".format(z1/self.scale)) + print("Period: {}".format((z1-z0)/self.scale)) + return mesh + +surf1 = ExplicitSurfaceThing( + freq = 20, + phase = 0, + scale = 1/16, # from libfive + inner = 0.4 * 1/16, + outer = 2.0 * 1/16, + rad = 0.3 * 1/16, + ext_phase = 0) + +z_init = 0 +x_init = surf1.x_cross(z_init, 1) +mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4) +fname = "spiral_inner0_one_period.stl" +print("Writing {}...".format(fname)) +mesh1.save(fname) + +surf2 = ExplicitSurfaceThing( + freq = 10, + phase = 0, + scale = 1/16, # from libfive + inner = 0.9 * 1/16, + outer = 2.0 * 1/16, + rad = 0.3 * 1/16, + ext_phase = numpy.pi/2) + +z_init = 0 +x_init = surf2.x_cross(z_init, 1) +mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4) +fname = "spiral_outer90_one_period.stl" +print("Writing {}...".format(fname)) +mesh2.save(fname)