diff --git a/src/libslic3r/Triangulation.cpp b/src/libslic3r/Triangulation.cpp index fc811277e5..9119df7d9b 100644 --- a/src/libslic3r/Triangulation.cpp +++ b/src/libslic3r/Triangulation.cpp @@ -26,7 +26,9 @@ inline void insert_edge(Triangulation::HalfEdges &edges, uint32_t &offset, const } // namespace Private -Triangulation::Indices Triangulation::triangulate(const Points &points, const HalfEdges &half_edges) +Triangulation::Indices Triangulation::triangulate(const Points & points, + const HalfEdges &half_edges, + bool allow_opposit_edge) { // IMPROVE use int point insted of float !!! @@ -67,12 +69,11 @@ Triangulation::Indices Triangulation::triangulate(const Points &points, const Ha for (size_t i = 0; i < 3; ++i) pi[i] = map[face->vertex(i)]; // Do not use triangles with opposit edges - if (half_edges.find(std::make_pair(pi[1], pi[0])) != half_edges.end()) - continue; - if (half_edges.find(std::make_pair(pi[2], pi[1])) != half_edges.end()) - continue; - if (half_edges.find(std::make_pair(pi[0], pi[2])) != half_edges.end()) - continue; + if (!allow_opposit_edge) { + if (half_edges.find(std::make_pair(pi[1], pi[0])) != half_edges.end()) continue; + if (half_edges.find(std::make_pair(pi[2], pi[1])) != half_edges.end()) continue; + if (half_edges.find(std::make_pair(pi[0], pi[2])) != half_edges.end()) continue; + } indices.emplace_back(pi[0], pi[1], pi[2]); } diff --git a/src/libslic3r/Triangulation.hpp b/src/libslic3r/Triangulation.hpp index a143587232..ccc08b52a6 100644 --- a/src/libslic3r/Triangulation.hpp +++ b/src/libslic3r/Triangulation.hpp @@ -26,9 +26,12 @@ public: /// /// Points to connect /// Constraint for edges, pair is from point(first) to - /// point(second) Triangles + /// point(second) + /// Flag for filtration result indices by opposit half edge + /// Triangles static Indices triangulate(const Points & points, - const HalfEdges &half_edges); + const HalfEdges &half_edges, + bool allow_opposit_edge = false); static Indices triangulate(const Polygon &polygon); static Indices triangulate(const Polygons &polygons); static Indices triangulate(const ExPolygons &expolygons); diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt index 05898db28a..3f732d1a39 100644 --- a/tests/libslic3r/CMakeLists.txt +++ b/tests/libslic3r/CMakeLists.txt @@ -22,6 +22,9 @@ add_executable(${_TEST_NAME}_tests test_optimizers.cpp test_png_io.cpp test_timeutils.cpp + test_quadric_edge_collapse.cpp + test_triangulation.cpp + test_emboss.cpp test_indexed_triangle_set.cpp ../libnest2d/printer_parts.cpp ) diff --git a/tests/libslic3r/Simplify.h b/tests/libslic3r/Simplify.h new file mode 100644 index 0000000000..24c126e34e --- /dev/null +++ b/tests/libslic3r/Simplify.h @@ -0,0 +1,1036 @@ +///////////////////////////////////////////// +// +// Mesh Simplification Tutorial +// +// (C) by Sven Forstmann in 2014 +// +// License : MIT +// http://opensource.org/licenses/MIT +// +//https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification +// +// 5/2016: Chris Rorden created minimal version for OSX/Linux/Windows compile + +//#include +//#include +//#include +//#include +//#include +#include +//#include +//#include +#include +#include +#include +#include +#include +#include +#include //FLT_EPSILON, DBL_EPSILON + +#define loopi(start_l,end_l) for ( int i=start_l;i1) input=1; + return (double) acos ( input ); + } + + inline double angle2( const vec3f& v , const vec3f& w ) + { + vec3f a = v , b= *this; + double dot = a.x*b.x + a.y*b.y + a.z*b.z; + double len = a.length() * b.length(); + if(len==0)len=1; + + vec3f plane; plane.cross( b,w ); + + if ( plane.x * a.x + plane.y * a.y + plane.z * a.z > 0 ) + return (double) -acos ( dot / len ); + + return (double) acos ( dot / len ); + } + + inline vec3f rot_x( double a ) + { + double yy = cos ( a ) * y + sin ( a ) * z; + double zz = cos ( a ) * z - sin ( a ) * y; + y = yy; z = zz; + return *this; + } + inline vec3f rot_y( double a ) + { + double xx = cos ( -a ) * x + sin ( -a ) * z; + double zz = cos ( -a ) * z - sin ( -a ) * x; + x = xx; z = zz; + return *this; + } + inline void clamp( double min, double max ) + { + if (xmax) x=max; + if (y>max) y=max; + if (z>max) z=max; + } + inline vec3f rot_z( double a ) + { + double yy = cos ( a ) * y + sin ( a ) * x; + double xx = cos ( a ) * x - sin ( a ) * y; + y = yy; x = xx; + return *this; + } + inline vec3f invert() + { + x=-x;y=-y;z=-z;return *this; + } + inline vec3f frac() + { + return vec3f( + x-double(int(x)), + y-double(int(y)), + z-double(int(z)) + ); + } + + inline vec3f integer() + { + return vec3f( + double(int(x)), + double(int(y)), + double(int(z)) + ); + } + + inline double length() const + { + return (double)sqrt(x*x + y*y + z*z); + } + + inline vec3f normalize( double desired_length = 1 ) + { + double square = sqrt(x*x + y*y + z*z); + /* + if (square <= 0.00001f ) + { + x=1;y=0;z=0; + return *this; + }*/ + //double len = desired_length / square; + x/=square;y/=square;z/=square; + + return *this; + } + static vec3f normalize( vec3f a ); + + static void random_init(); + static double random_double(); + static vec3f random(); + + static int random_number; + + double random_double_01(double a){ + double rnf=a*14.434252+a*364.2343+a*4213.45352+a*2341.43255+a*254341.43535+a*223454341.3523534245+23453.423412; + int rni=((int)rnf)%100000; + return double(rni)/(100000.0f-1.0f); + } + + vec3f random01_fxyz(){ + x=(double)random_double_01(x); + y=(double)random_double_01(y); + z=(double)random_double_01(z); + return *this; + } + +}; + +vec3f barycentric(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c){ + vec3f v0 = b-a; + vec3f v1 = c-a; + vec3f v2 = p-a; + double d00 = v0.dot(v0); + double d01 = v0.dot(v1); + double d11 = v1.dot(v1); + double d20 = v2.dot(v0); + double d21 = v2.dot(v1); + double denom = d00*d11-d01*d01; + double v = (d11 * d20 - d01 * d21) / denom; + double w = (d00 * d21 - d01 * d20) / denom; + double u = 1.0 - v - w; + return vec3f(u,v,w); +} + +vec3f interpolate(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c, const vec3f attrs[3]) +{ + vec3f bary = barycentric(p,a,b,c); + vec3f out = vec3f(0,0,0); + out = out + attrs[0] * bary.x; + out = out + attrs[1] * bary.y; + out = out + attrs[2] * bary.z; + return out; +} + +double min(double v1, double v2) { + return fmin(v1,v2); +} + + +class SymetricMatrix { + + public: + + // Constructor + + SymetricMatrix(double c=0) { loopi(0,10) m[i] = c; } + + SymetricMatrix( double m11, double m12, double m13, double m14, + double m22, double m23, double m24, + double m33, double m34, + double m44) { + m[0] = m11; m[1] = m12; m[2] = m13; m[3] = m14; + m[4] = m22; m[5] = m23; m[6] = m24; + m[7] = m33; m[8] = m34; + m[9] = m44; + } + + // Make plane + + SymetricMatrix(double a,double b,double c,double d) + { + m[0] = a*a; m[1] = a*b; m[2] = a*c; m[3] = a*d; + m[4] = b*b; m[5] = b*c; m[6] = b*d; + m[7 ] =c*c; m[8 ] = c*d; + m[9 ] = d*d; + } + + double operator[](int c) const { return m[c]; } + + // Determinant + + double det( int a11, int a12, int a13, + int a21, int a22, int a23, + int a31, int a32, int a33) + { + double det = m[a11]*m[a22]*m[a33] + m[a13]*m[a21]*m[a32] + m[a12]*m[a23]*m[a31] + - m[a13]*m[a22]*m[a31] - m[a11]*m[a23]*m[a32]- m[a12]*m[a21]*m[a33]; + return det; + } + + const SymetricMatrix operator+(const SymetricMatrix& n) const + { + return SymetricMatrix( m[0]+n[0], m[1]+n[1], m[2]+n[2], m[3]+n[3], + m[4]+n[4], m[5]+n[5], m[6]+n[6], + m[ 7]+n[ 7], m[ 8]+n[8 ], + m[ 9]+n[9 ]); + } + + SymetricMatrix& operator+=(const SymetricMatrix& n) + { + m[0]+=n[0]; m[1]+=n[1]; m[2]+=n[2]; m[3]+=n[3]; + m[4]+=n[4]; m[5]+=n[5]; m[6]+=n[6]; m[7]+=n[7]; + m[8]+=n[8]; m[9]+=n[9]; + return *this; + } + + double m[10]; +}; +/////////////////////////////////////////// + +namespace Simplify +{ + // Global Variables & Strctures + enum Attributes { + NONE, + NORMAL = 2, + TEXCOORD = 4, + COLOR = 8 + }; + struct Triangle { int v[3];double err[4];int deleted,dirty,attr;vec3f n;vec3f uvs[3];int material; }; + struct Vertex { vec3f p;int tstart,tcount;SymetricMatrix q;int border;}; + struct Ref { int tid,tvertex; }; + std::vector triangles; + std::vector vertices; + std::vector refs; + std::string mtllib; + std::vector materials; + + // Helper functions + + double vertex_error(SymetricMatrix q, double x, double y, double z); + double calculate_error(int id_v1, int id_v2, vec3f &p_result); + bool flipped(vec3f p,int i0,int i1,Vertex &v0,Vertex &v1,std::vector &deleted); + void update_uvs(int i0,const Vertex &v,const vec3f &p,std::vector &deleted); + void update_triangles(int i0,Vertex &v,std::vector &deleted,int &deleted_triangles); + void update_mesh(int iteration); + void compact_mesh(); + // + // Main simplification function + // + // target_count : target nr. of triangles + // agressiveness : sharpness to increase the threshold. + // 5..8 are good numbers + // more iterations yield higher quality + // + + void simplify_mesh(int target_count, double agressiveness=7, bool verbose=false) + { + // init + loopi(0,triangles.size()) + { + triangles[i].deleted=0; + } + + // main iteration loop + int deleted_triangles=0; + std::vector deleted0,deleted1; + int triangle_count=triangles.size(); + //int iteration = 0; + //loop(iteration,0,100) + for (int iteration = 0; iteration < 100; iteration ++) + { + if(triangle_count-deleted_triangles<=target_count)break; + + // update mesh once in a while + if(iteration%5==0) + { + update_mesh(iteration); + } + + // clear dirty flag + loopi(0,triangles.size()) triangles[i].dirty=0; + + // + // All triangles with edges below the threshold will be removed + // + // The following numbers works well for most models. + // If it does not, try to adjust the 3 parameters + // + double threshold = 0.000000001*pow(double(iteration+3),agressiveness); + + // target number of triangles reached ? Then break + if ((verbose) && (iteration%5==0)) { + printf("iteration %d - triangles %d threshold %g\n",iteration,triangle_count-deleted_triangles, threshold); + } + + // remove vertices & mark deleted triangles + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + if(t.err[3]>threshold) continue; + if(t.deleted) continue; + if(t.dirty) continue; + + loopj(0,3)if(t.err[j] deleted0,deleted1; + int triangle_count=triangles.size(); + //int iteration = 0; + //loop(iteration,0,100) + for (int iteration = 0; iteration < 9999; iteration ++) + { + // update mesh constantly + update_mesh(iteration); + // clear dirty flag + loopi(0,triangles.size()) triangles[i].dirty=0; + // + // All triangles with edges below the threshold will be removed + // + // The following numbers works well for most models. + // If it does not, try to adjust the 3 parameters + // + double threshold = DBL_EPSILON; //1.0E-3 EPS; + if (verbose) { + printf("lossless iteration %d\n", iteration); + } + + // remove vertices & mark deleted triangles + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + if(t.err[3]>threshold) continue; + if(t.deleted) continue; + if(t.dirty) continue; + + loopj(0,3)if(t.err[j] &deleted) + { + + loopk(0,v0.tcount) + { + Triangle &t=triangles[refs[v0.tstart+k].tid]; + if(t.deleted)continue; + + int s=refs[v0.tstart+k].tvertex; + int id1=t.v[(s+1)%3]; + int id2=t.v[(s+2)%3]; + + if(id1==i1 || id2==i1) // delete ? + { + + deleted[k]=1; + continue; + } + vec3f d1 = vertices[id1].p-p; d1.normalize(); + vec3f d2 = vertices[id2].p-p; d2.normalize(); + if(fabs(d1.dot(d2))>0.999) return true; + vec3f n; + n.cross(d1,d2); + n.normalize(); + deleted[k]=0; + if(n.dot(t.n)<0.2) return true; + } + return false; + } + + // update_uvs + + void update_uvs(int i0,const Vertex &v,const vec3f &p,std::vector &deleted) + { + loopk(0,v.tcount) + { + Ref &r=refs[v.tstart+k]; + Triangle &t=triangles[r.tid]; + if(t.deleted)continue; + if(deleted[k])continue; + vec3f p1=vertices[t.v[0]].p; + vec3f p2=vertices[t.v[1]].p; + vec3f p3=vertices[t.v[2]].p; + t.uvs[r.tvertex] = interpolate(p,p1,p2,p3,t.uvs); + } + } + + // Update triangle connections and edge error after a edge is collapsed + + void update_triangles(int i0,Vertex &v,std::vector &deleted,int &deleted_triangles) + { + vec3f p; + loopk(0,v.tcount) + { + Ref &r=refs[v.tstart+k]; + Triangle &t=triangles[r.tid]; + if(t.deleted)continue; + if(deleted[k]) + { + t.deleted=1; + deleted_triangles++; + continue; + } + t.v[r.tvertex]=i0; + t.dirty=1; + t.err[0]=calculate_error(t.v[0],t.v[1],p); + t.err[1]=calculate_error(t.v[1],t.v[2],p); + t.err[2]=calculate_error(t.v[2],t.v[0],p); + t.err[3]=min(t.err[0],min(t.err[1],t.err[2])); + refs.push_back(r); + } + } + + // compact triangles, compute edge error and build reference list + + void update_mesh(int iteration) + { + if(iteration>0) // compact triangles + { + int dst=0; + loopi(0,triangles.size()) + if(!triangles[i].deleted) + { + triangles[dst++]=triangles[i]; + } + triangles.resize(dst); + } + // + // Init Quadrics by Plane & Edge Errors + // + // required at the beginning ( iteration == 0 ) + // recomputing during the simplification is not required, + // but mostly improves the result for closed meshes + // + if( iteration == 0 ) + { + loopi(0,vertices.size()) + vertices[i].q=SymetricMatrix(0.0); + + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + vec3f n,p[3]; + loopj(0,3) p[j]=vertices[t.v[j]].p; + n.cross(p[1]-p[0],p[2]-p[0]); + n.normalize(); + t.n=n; + loopj(0,3) vertices[t.v[j]].q = + vertices[t.v[j]].q+SymetricMatrix(n.x,n.y,n.z,-n.dot(p[0])); + } + loopi(0,triangles.size()) + { + // Calc Edge Error + Triangle &t=triangles[i];vec3f p; + loopj(0,3) t.err[j]=calculate_error(t.v[j],t.v[(j+1)%3],p); + t.err[3]=min(t.err[0],min(t.err[1],t.err[2])); + } + } + + // Init Reference ID list + loopi(0,vertices.size()) + { + vertices[i].tstart=0; + vertices[i].tcount=0; + } + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + loopj(0,3) vertices[t.v[j]].tcount++; + } + int tstart=0; + loopi(0,vertices.size()) + { + Vertex &v=vertices[i]; + v.tstart=tstart; + tstart+=v.tcount; + v.tcount=0; + } + + // Write References + refs.resize(triangles.size()*3); + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + loopj(0,3) + { + Vertex &v=vertices[t.v[j]]; + refs[v.tstart+v.tcount].tid=i; + refs[v.tstart+v.tcount].tvertex=j; + v.tcount++; + } + } + + // Identify boundary : vertices[].border=0,1 + if( iteration == 0 ) + { + std::vector vcount,vids; + + loopi(0,vertices.size()) + vertices[i].border=0; + + loopi(0,vertices.size()) + { + Vertex &v=vertices[i]; + vcount.clear(); + vids.clear(); + loopj(0,v.tcount) + { + int k=refs[v.tstart+j].tid; + Triangle &t=triangles[k]; + loopk(0,3) + { + int ofs=0,id=t.v[k]; + while(ofs try to find best result + vec3f p1=vertices[id_v1].p; + vec3f p2=vertices[id_v2].p; + vec3f p3=(p1+p2)/2; + double error1 = vertex_error(q, p1.x,p1.y,p1.z); + double error2 = vertex_error(q, p2.x,p2.y,p2.z); + double error3 = vertex_error(q, p3.x,p3.y,p3.z); + error = min(error1, min(error2, error3)); + if (error1 == error) p_result=p1; + if (error2 == error) p_result=p2; + if (error3 == error) p_result=p3; + } + return error; + } + + char *trimwhitespace(char *str) + { + char *end; + + // Trim leading space + while(isspace((unsigned char)*str)) str++; + + if(*str == 0) // All spaces? + return str; + + // Trim trailing space + end = str + strlen(str) - 1; + while(end > str && isspace((unsigned char)*end)) end--; + + // Write new null terminator + *(end+1) = 0; + + return str; + } + + //Option : Load OBJ + void load_obj(const char* filename, bool process_uv=false){ + vertices.clear(); + triangles.clear(); + //printf ( "Loading Objects %s ... \n",filename); + FILE* fn; + if(filename==NULL) return ; + if((char)filename[0]==0) return ; + if ((fn = fopen(filename, "rb")) == NULL) + { + printf ( "File %s not found!\n" ,filename ); + return; + } + char line[1000]; + memset ( line,0,1000 ); + int vertex_cnt = 0; + int material = -1; + std::map material_map; + std::vector uvs; + std::vector > uvMap; + + while(fgets( line, 1000, fn ) != NULL) + { + Vertex v; + vec3f uv; + + if (strncmp(line, "mtllib", 6) == 0) + { + mtllib = trimwhitespace(&line[7]); + } + if (strncmp(line, "usemtl", 6) == 0) + { + std::string usemtl = trimwhitespace(&line[7]); + if (material_map.find(usemtl) == material_map.end()) + { + material_map[usemtl] = materials.size(); + materials.push_back(usemtl); + } + material = material_map[usemtl]; + } + + if ( line[0] == 'v' && line[1] == 't' ) + { + if ( line[2] == ' ' ) + if(sscanf(line,"vt %lf %lf", + &uv.x,&uv.y)==2) + { + uv.z = 0; + uvs.push_back(uv); + } else + if(sscanf(line,"vt %lf %lf %lf", + &uv.x,&uv.y,&uv.z)==3) + { + uvs.push_back(uv); + } + } + else if ( line[0] == 'v' ) + { + if ( line[1] == ' ' ) + if(sscanf(line,"v %lf %lf %lf", + &v.p.x, &v.p.y, &v.p.z)==3) + { + vertices.push_back(v); + } + } + int integers[9]; + if ( line[0] == 'f' ) + { + Triangle t; + bool tri_ok = false; + bool has_uv = false; + + if(sscanf(line,"f %d %d %d", + &integers[0],&integers[1],&integers[2])==3) + { + tri_ok = true; + }else + if(sscanf(line,"f %d// %d// %d//", + &integers[0],&integers[1],&integers[2])==3) + { + tri_ok = true; + }else + if(sscanf(line,"f %d//%d %d//%d %d//%d", + &integers[0],&integers[3], + &integers[1],&integers[4], + &integers[2],&integers[5])==6) + { + tri_ok = true; + }else + if(sscanf(line,"f %d/%d/%d %d/%d/%d %d/%d/%d", + &integers[0],&integers[6],&integers[3], + &integers[1],&integers[7],&integers[4], + &integers[2],&integers[8],&integers[5])==9) + { + tri_ok = true; + has_uv = true; + }else // Add Support for v/vt only meshes + if (sscanf(line, "f %d/%d %d/%d %d/%d", + &integers[0], &integers[6], + &integers[1], &integers[7], + &integers[2], &integers[8]) == 6) + { + tri_ok = true; + has_uv = true; + } + else + { + printf("unrecognized sequence\n"); + printf("%s\n",line); + while(1); + } + if ( tri_ok ) + { + t.v[0] = integers[0]-1-vertex_cnt; + t.v[1] = integers[1]-1-vertex_cnt; + t.v[2] = integers[2]-1-vertex_cnt; + t.attr = 0; + + if ( process_uv && has_uv ) + { + std::vector indices; + indices.push_back(integers[6]-1-vertex_cnt); + indices.push_back(integers[7]-1-vertex_cnt); + indices.push_back(integers[8]-1-vertex_cnt); + uvMap.push_back(indices); + t.attr |= TEXCOORD; + } + + t.material = material; + //geo.triangles.push_back ( tri ); + triangles.push_back(t); + //state_before = state; + //state ='f'; + } + } + } + + if ( process_uv && uvs.size() ) + { + loopi(0,triangles.size()) + { + loopj(0,3) + triangles[i].uvs[j] = uvs[uvMap[i][j]]; + } + } + + fclose(fn); + + //printf("load_obj: vertices = %lu, triangles = %lu, uvs = %lu\n", vertices.size(), triangles.size(), uvs.size() ); + } // load_obj() + + // Optional : Store as OBJ + + void write_obj(const char* filename) + { + FILE *file=fopen(filename, "w"); + int cur_material = -1; + bool has_uv = (triangles.size() && (triangles[0].attr & TEXCOORD) == TEXCOORD); + + if (!file) + { + printf("write_obj: can't write data file \"%s\".\n", filename); + exit(0); + } + if (!mtllib.empty()) + { + fprintf(file, "mtllib %s\n", mtllib.c_str()); + } + loopi(0,vertices.size()) + { + //fprintf(file, "v %lf %lf %lf\n", vertices[i].p.x,vertices[i].p.y,vertices[i].p.z); + fprintf(file, "v %g %g %g\n", vertices[i].p.x,vertices[i].p.y,vertices[i].p.z); //more compact: remove trailing zeros + } + if (has_uv) + { + loopi(0,triangles.size()) if(!triangles[i].deleted) + { + fprintf(file, "vt %g %g\n", triangles[i].uvs[0].x, triangles[i].uvs[0].y); + fprintf(file, "vt %g %g\n", triangles[i].uvs[1].x, triangles[i].uvs[1].y); + fprintf(file, "vt %g %g\n", triangles[i].uvs[2].x, triangles[i].uvs[2].y); + } + } + int uv = 1; + loopi(0,triangles.size()) if(!triangles[i].deleted) + { + if (triangles[i].material != cur_material) + { + cur_material = triangles[i].material; + fprintf(file, "usemtl %s\n", materials[triangles[i].material].c_str()); + } + if (has_uv) + { + fprintf(file, "f %d/%d %d/%d %d/%d\n", triangles[i].v[0]+1, uv, triangles[i].v[1]+1, uv+1, triangles[i].v[2]+1, uv+2); + uv += 3; + } + else + { + fprintf(file, "f %d %d %d\n", triangles[i].v[0]+1, triangles[i].v[1]+1, triangles[i].v[2]+1); + } + //fprintf(file, "f %d// %d// %d//\n", triangles[i].v[0]+1, triangles[i].v[1]+1, triangles[i].v[2]+1); //more compact: remove trailing zeros + } + fclose(file); + } +}; +/////////////////////////////////////////// diff --git a/tests/libslic3r/test_emboss.cpp b/tests/libslic3r/test_emboss.cpp new file mode 100644 index 0000000000..ff427cce56 --- /dev/null +++ b/tests/libslic3r/test_emboss.cpp @@ -0,0 +1,181 @@ +#include + +#include +#include // only debug visualization + +#include +#include + +using namespace Slic3r; + +namespace Private{ + +// calculate multiplication of ray dir to intersect - inspired by +// segment_segment_intersection when ray dir is normalized retur distance from +// ray point to intersection No value mean no intersection +std::optional ray_segment_intersection(const Vec2d &r_point, + const Vec2d &r_dir, + const Vec2d &s0, + const Vec2d &s1) +{ + auto denominate = [](const Vec2d &v0, const Vec2d &v1) -> double { + return v0.x() * v1.y() - v1.x() * v0.y(); + }; + + Vec2d segment_dir = s1 - s0; + double d = denominate(segment_dir, r_dir); + if (std::abs(d) < std::numeric_limits::epsilon()) + // Line and ray are collinear. + return {}; + + Vec2d s12 = s0 - r_point; + double s_number = denominate(r_dir, s12); + bool change_sign = false; + if (d < 0.) { + change_sign = true; + d = -d; + s_number = -s_number; + } + + if (s_number < 0. || s_number > d) + // Intersection outside of segment. + return {}; + + double r_number = denominate(segment_dir, s12); + if (change_sign) r_number = -r_number; + + if (r_number < 0.) + // Intersection before ray start. + return {}; + + return r_number / d; +} + +Vec2d get_intersection(const Vec2d & point, + const Vec2d & dir, + const std::array &triangle) +{ + std::optional t; + for (size_t i = 0; i < 3; ++i) { + size_t i2 = i + 1; + if (i2 == 3) i2 = 0; + if (!t.has_value()) { + t = ray_segment_intersection(point, dir, triangle[i], + triangle[i2]); + continue; + } + + // small distance could be preccission inconsistance + std::optional t2 = ray_segment_intersection(point, dir, + triangle[i], + triangle[i2]); + if (t2.has_value() && *t2 > *t) t = t2; + } + assert(t.has_value()); // Not found intersection. + return point + dir * (*t); +} + +Vec3d calc_hit_point(const igl::Hit & h, + const Vec3i & triangle, + const std::vector &vertices) +{ + double c1 = h.u; + double c2 = h.v; + double c0 = 1.0 - c1 - c2; + Vec3d v0 = vertices[triangle[0]].cast(); + Vec3d v1 = vertices[triangle[1]].cast(); + Vec3d v2 = vertices[triangle[2]].cast(); + return v0 * c0 + v1 * c1 + v2 * c2; +} + +Vec3d calc_hit_point(const igl::Hit &h, indexed_triangle_set &its) +{ + return calc_hit_point(h, its.indices[h.id], its.vertices); +} +} // namespace Private + +TEST_CASE("Emboss text", "[Emboss]") +{ + const char *font_name = "C:/windows/fonts/arialbd.ttf"; + char letter = '%'; + float flatness = 2.; + + std::optional font = Emboss::load_font(font_name); + REQUIRE(font.has_value()); + + std::optional glyph = Emboss::letter2glyph(*font, letter, flatness); + REQUIRE(glyph.has_value()); + + ExPolygons shape = glyph->shape; + REQUIRE(!shape.empty()); + + float z_depth = 1.f; + Emboss::ProjectZ projection(z_depth); + indexed_triangle_set its = Emboss::polygons2model(shape, projection); + + CHECK(!its.indices.empty()); +} + +TEST_CASE("Test hit point", "[AABBTreeIndirect]") +{ + indexed_triangle_set its; + its.vertices = { + Vec3f(1, 1, 1), + Vec3f(2, 10, 2), + Vec3f(10, 0, 2), + }; + its.indices = {Vec3i(0, 2, 1)}; + auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set( + its.vertices, its.indices); + + Vec3d ray_point(8, 1, 0); + Vec3d ray_dir(0, 0, 1); + igl::Hit hit; + AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices, tree, + ray_point, ray_dir, hit); + Vec3d hp = Private::calc_hit_point(hit, its); + CHECK(abs(hp.x() - ray_point.x()) < .1); + CHECK(abs(hp.y() - ray_point.y()) < .1); +} + +TEST_CASE("ray segment intersection", "[MeshBoolean]") +{ + Vec2d r_point(1, 1); + Vec2d r_dir(1, 0); + + // colinear + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 0), Vec2d(2, 0)).has_value()); + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(0, 0)).has_value()); + + // before ray + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 0), Vec2d(0, 2)).has_value()); + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 2), Vec2d(0, 0)).has_value()); + + // above ray + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 2), Vec2d(2, 3)).has_value()); + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 3), Vec2d(2, 2)).has_value()); + + // belove ray + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(2, -1)).has_value()); + CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, -1), Vec2d(2, 0)).has_value()); + + // intersection at [2,1] distance 1 + auto t1 = Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(2, 2)); + REQUIRE(t1.has_value()); + auto t2 = Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 2), Vec2d(2, 0)); + REQUIRE(t2.has_value()); + + CHECK(abs(*t1 - *t2) < std::numeric_limits::epsilon()); +} + + + +TEST_CASE("triangle intersection", "[]") +{ + Vec2d point(1, 1); + Vec2d dir(-1, 0); + std::array triangle = {Vec2d(0, 0), Vec2d(5, 0), Vec2d(0, 5)}; + Vec2d i = Private::get_intersection(point, dir, triangle); + CHECK(abs(i.x()) < std::numeric_limits::epsilon()); + CHECK(abs(i.y() - 1.) < std::numeric_limits::epsilon()); +} \ No newline at end of file diff --git a/tests/libslic3r/test_indexed_triangle_set.cpp b/tests/libslic3r/test_indexed_triangle_set.cpp index 6a8cdfc05b..91b43f8e83 100644 --- a/tests/libslic3r/test_indexed_triangle_set.cpp +++ b/tests/libslic3r/test_indexed_triangle_set.cpp @@ -101,220 +101,3 @@ TEST_CASE("Split two watertight meshes", "[its_split][its]") { debug_write_obj(res, "parts_watertight"); } - -#include -static float triangle_area(const Vec3f &v0, const Vec3f &v1, const Vec3f &v2) -{ - Vec3f ab = v1 - v0; - Vec3f ac = v2 - v0; - return ab.cross(ac).norm() / 2.f; -} - -static float triangle_area(const Vec3crd &triangle_inices, const std::vector &vertices) -{ - return triangle_area(vertices[triangle_inices[0]], - vertices[triangle_inices[1]], - vertices[triangle_inices[2]]); -} - -static std::mt19937 create_random_generator() { - std::random_device rd; - std::mt19937 gen(rd()); - return gen; -} - -std::vector its_sample_surface(const indexed_triangle_set &its, - double sample_per_mm2, - std::mt19937 random_generator = create_random_generator()) -{ - std::vector samples; - std::uniform_real_distribution rand01(0.f, 1.f); - for (const auto &triangle_indices : its.indices) { - float area = triangle_area(triangle_indices, its.vertices); - float countf; - float fractional = std::modf(area * sample_per_mm2, &countf); - int count = static_cast(countf); - - float generate = rand01(random_generator); - if (generate < fractional) ++count; - if (count == 0) continue; - - const Vec3f &v0 = its.vertices[triangle_indices[0]]; - const Vec3f &v1 = its.vertices[triangle_indices[1]]; - const Vec3f &v2 = its.vertices[triangle_indices[2]]; - for (int c = 0; c < count; c++) { - // barycentric coordinate - Vec3f b; - b[0] = rand01(random_generator); - b[1] = rand01(random_generator); - if ((b[0] + b[1]) > 1.f) { - b[0] = 1.f - b[0]; - b[1] = 1.f - b[1]; - } - b[2] = 1.f - b[0] - b[1]; - Vec3f pos; - for (int i = 0; i < 3; i++) { - pos[i] = b[0] * v0[i] + b[1] * v1[i] + b[2] * v2[i]; - } - samples.push_back(pos); - } - } - return samples; -} - - -#include "libslic3r/AABBTreeIndirect.hpp" - -struct CompareConfig -{ - float max_distance = 3.f; - float max_average_distance = 2.f; -}; - -bool is_similar(const indexed_triangle_set &from, - const indexed_triangle_set &to, - const CompareConfig &cfg) -{ - // create ABBTree - auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set( - from.vertices, from.indices); - float sum_distance = 0.f; - float max_distance = 0.f; - - auto collect_distances = [&](const Vec3f &surface_point) { - size_t hit_idx; - Vec3f hit_point; - float distance2 = - AABBTreeIndirect::squared_distance_to_indexed_triangle_set( - from.vertices, from.indices, tree, surface_point, hit_idx, hit_point); - float distance = sqrt(distance2); - if (max_distance < distance) max_distance = distance; - sum_distance += distance; - }; - - for (const Vec3f &vertex : to.vertices) { - collect_distances(vertex); - } - - for (const Vec3i &t : to.indices) { - Vec3f center(0,0,0); - for (size_t i = 0; i < 3; ++i) { - center += to.vertices[t[i]] / 3; - } - collect_distances(center); - } - - size_t count = to.vertices.size() + to.indices.size(); - float avg_distance = sum_distance / count; - if (avg_distance > cfg.max_average_distance || - max_distance > cfg.max_distance) - return false; - return true; -} - -TEST_CASE("Reduce one edge by Quadric Edge Collapse", "[its]") -{ - indexed_triangle_set its; - its.vertices = {Vec3f(-1.f, 0.f, 0.f), Vec3f(0.f, 1.f, 0.f), - Vec3f(1.f, 0.f, 0.f), Vec3f(0.f, 0.f, 1.f), - // vertex to be removed - Vec3f(0.9f, .1f, -.1f)}; - its.indices = {Vec3i(1, 0, 3), Vec3i(2, 1, 3), Vec3i(0, 2, 3), - Vec3i(0, 1, 4), Vec3i(1, 2, 4), Vec3i(2, 0, 4)}; - // edge to remove is between vertices 2 and 4 on trinagles 4 and 5 - - indexed_triangle_set its_ = its; // copy - // its_write_obj(its, "tetrhedron_in.obj"); - uint32_t wanted_count = its.indices.size() - 1; - its_quadric_edge_collapse(its, wanted_count); - // its_write_obj(its, "tetrhedron_out.obj"); - CHECK(its.indices.size() == 4); - CHECK(its.vertices.size() == 4); - - for (size_t i = 0; i < 3; i++) { - CHECK(its.indices[i] == its_.indices[i]); - } - - for (size_t i = 0; i < 4; i++) { - if (i == 2) continue; - CHECK(its.vertices[i] == its_.vertices[i]); - } - - const Vec3f &v = its.vertices[2]; // new vertex - const Vec3f &v2 = its_.vertices[2]; // moved vertex - const Vec3f &v4 = its_.vertices[4]; // removed vertex - for (size_t i = 0; i < 3; i++) { - bool is_between = (v[i] < v4[i] && v[i] > v2[i]) || - (v[i] > v4[i] && v[i] < v2[i]); - CHECK(is_between); - } - CompareConfig cfg; - cfg.max_average_distance = 0.014f; - cfg.max_distance = 0.75f; - - CHECK(is_similar(its, its_, cfg)); - CHECK(is_similar(its_, its, cfg)); -} - -#include "test_utils.hpp" -TEST_CASE("Simplify mesh by Quadric edge collapse to 5%", "[its]") -{ - TriangleMesh mesh = load_model("frog_legs.obj"); - double original_volume = its_volume(mesh.its); - uint32_t wanted_count = mesh.its.indices.size() * 0.05; - REQUIRE_FALSE(mesh.empty()); - indexed_triangle_set its = mesh.its; // copy - float max_error = std::numeric_limits::max(); - its_quadric_edge_collapse(its, wanted_count, &max_error); - //its_write_obj(its, "frog_legs_qec.obj"); - CHECK(its.indices.size() <= wanted_count); - double volume = its_volume(its); - CHECK(fabs(original_volume - volume) < 33.); - - CompareConfig cfg; - cfg.max_average_distance = 0.043f; - cfg.max_distance = 0.32f; - - CHECK(is_similar(mesh.its, its, cfg)); - CHECK(is_similar(its, mesh.its, cfg)); -} - -bool exist_triangle_with_twice_vertices(const std::vector& indices) -{ - for (const auto &face : indices) - if (face[0] == face[1] || - face[0] == face[2] || - face[1] == face[2]) return true; - return false; -} - -TEST_CASE("Simplify trouble case", "[its]") -{ - TriangleMesh tm = load_model("simplification.obj"); - REQUIRE_FALSE(tm.empty()); - float max_error = std::numeric_limits::max(); - uint32_t wanted_count = 0; - its_quadric_edge_collapse(tm.its, wanted_count, &max_error); - CHECK(!exist_triangle_with_twice_vertices(tm.its.indices)); -} - -TEST_CASE("Simplified cube should not be empty.", "[its]") -{ - auto its = its_make_cube(1, 2, 3); - float max_error = std::numeric_limits::max(); - uint32_t wanted_count = 0; - its_quadric_edge_collapse(its, wanted_count, &max_error); - CHECK(!its.indices.empty()); -} - -TEST_CASE("Neighbors in cube", "[its]") -{ - auto its = its_make_cube(1, 2, 3); - auto neighbors = its_face_neighbors(its); - int no_value = -1; - for (auto &neighbor : neighbors) { - for (size_t i = 0; i < 3; i++) { - CHECK(neighbor[i] != no_value); - } - } -} diff --git a/tests/libslic3r/test_meshboolean.cpp b/tests/libslic3r/test_meshboolean.cpp index 3f9d648278..4dce6d2862 100644 --- a/tests/libslic3r/test_meshboolean.cpp +++ b/tests/libslic3r/test_meshboolean.cpp @@ -39,7 +39,7 @@ Vec3d calc_normal(const Vec3i &triangle, const std::vector &vertices) TEST_CASE("Add TriangleMeshes", "[MeshBoolean]") { TriangleMesh tm1 = make_sphere(1.6, 1.6); - Vec3f move(5, -3, 7); + Vec3f move(5, -3, 7); move.normalize(); tm1.translate(0.3 * move); its_write_obj(tm1.its, "tm1.obj"); @@ -48,493 +48,3 @@ TEST_CASE("Add TriangleMeshes", "[MeshBoolean]") MeshBoolean::cgal::plus(tm1, tm2); its_write_obj(tm1.its, "test_add.obj"); } - -#include "libslic3r/Emboss.hpp" - -ExPolygons ttf2polygons(const char * font_name, char letter, float flatness = 1.f) { - auto font = Emboss::load_font(font_name); - if (!font.has_value()) return ExPolygons(); - return Emboss::letter2glyph(*font, letter, flatness)->shape; -} - -#include "libslic3r/SVG.hpp" -void store_to_svg(Polygons polygons,std::string file_name = "letter.svg") -{ - double scale = 1e6; - BoundingBox bb; - for (auto& p : polygons) { - p.scale(scale); - bb.merge(p.points); - } - SVG svg(file_name, bb); - svg.draw(polygons); -} - -struct Plane -{ - Vec3d point = Vec3d(0., 0., 0.); // lay on plane - define zero position - Vec3d normal = Vec3d(0., 0., 1.);// [unit vector] - define orientation -}; - -struct OrientedPlane: public Plane -{ - // Must be perpendiculat to normal - Vec3d up = Vec3d(0., 1., 0.); // [unit vector] -}; - -struct EmbossConfig -{ - // emboss plane must be above surface - // point define zero for 2d polygon and must be out of model - // normal is direction to model - // up define orientation of polygon - OrientedPlane projection; - - - // Move surface distance - // Positive value out of model (Raised) - // Negative value into model (Engraved) - float height = 1.; // [in milimeters] -}; - -#include -#include -Vec3d calc_hit_point(const igl::Hit & h, - const Vec3i & triangle, - const std::vector &vertices) -{ - double c1 = h.u; - double c2 = h.v; - double c0 = 1.0 - c1 - c2; - Vec3d v0 = vertices[triangle[0]].cast(); - Vec3d v1 = vertices[triangle[1]].cast(); - Vec3d v2 = vertices[triangle[2]].cast(); - return v0 * c0 + v1 * c1 + v2 * c2; -} - -Vec3d calc_hit_point(const igl::Hit &h, indexed_triangle_set &its) -{ - return calc_hit_point(h, its.indices[h.id], its.vertices); -} - -TEST_CASE("Test hit point", "[AABBTreeIndirect]") -{ - indexed_triangle_set its; - its.vertices = { - Vec3f(1,1,1), - Vec3f(2, 10, 2), - Vec3f(10, 0, 2), - }; - its.indices = {Vec3i(0, 2, 1)}; - auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set( - its.vertices, its.indices); - - Vec3d ray_point(8, 1, 0); - Vec3d ray_dir(0,0,1); - igl::Hit hit; - AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices, tree, - ray_point, ray_dir, hit); - Vec3d hp = calc_hit_point(hit, its); - CHECK(abs(hp.x() - ray_point.x()) < .1); - CHECK(abs(hp.y() - ray_point.y()) < .1); -} - -// represents triangle extend by seam -struct TrianglePath -{ - // input edge, output edge, when cross triangle border - std::optional> edges; - - // when edges has value than first and last point lay on triangle edge - std::vector points; - - // first point has index offset_id in result vertices - uint32_t offset_id; -}; -using TrianglePaths = std::vector; - -// create transformation matrix to convert direction vectors -// do not care about up vector -// Directions are normalized -Eigen::Matrix3d create_transformation(const Vec3d &from_dir, const Vec3d &to_dir) -{ - Vec3d axis = from_dir.cross(to_dir); - axis.normalize(); - double angle = acos(from_dir.dot(to_dir)); - auto rotation = Eigen::AngleAxisd(angle, axis); - return rotation.matrix(); -} - -TEST_CASE("Transformation matrix", "[]") { - Vec3d d1(3, -7, 13); - Vec3d d2(-9, 5, 1); - d1.normalize(); - d2.normalize(); - auto tr_mat = create_transformation(d1, d2); - - Vec3d d1_tr = tr_mat * d1; - Vec3d diff = d1_tr - d2; - double eps = 1e-12; - for (double d : diff) - CHECK(abs(d) < std::numeric_limits::epsilon()); -} - -// calculate multiplication of ray dir to intersect - inspired by segment_segment_intersection -// when ray dir is normalized retur distance from ray point to intersection -// No value mean no intersection -std::optional ray_segment_intersection( - const Vec2d& r_point, const Vec2d& r_dir, const Vec2d& s0, const Vec2d& s1){ - - auto denominate = [](const Vec2d& v0,const Vec2d& v1)->double{ - return v0.x() * v1.y() - v1.x() * v0.y(); - }; - - Vec2d segment_dir = s1 - s0; - double d = denominate(segment_dir, r_dir); - if (std::abs(d) < std::numeric_limits::epsilon()) - // Line and ray are collinear. - return {}; - - Vec2d s12 = s0 - r_point; - double s_number = denominate(r_dir, s12); - bool change_sign = false; - if (d < 0.) { - change_sign = true; - d = -d; - s_number = -s_number; - } - - if (s_number < 0.|| s_number > d) - // Intersection outside of segment. - return {}; - - double r_number = denominate(segment_dir, s12); - if (change_sign) - r_number = - r_number; - - if(r_number < 0.) - // Intersection before ray start. - return {}; - - return r_number / d; -} - -TEST_CASE("ray segment intersection", "[MeshBoolean]") -{ - Vec2d r_point(1,1); - Vec2d r_dir(1,0); - - // colinear - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,0), Vec2d(2,0)).has_value()); - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(0,0)).has_value()); - - // before ray - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,0), Vec2d(0,2)).has_value()); - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,2), Vec2d(0,0)).has_value()); - - // above ray - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,2), Vec2d(2,3)).has_value()); - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,3), Vec2d(2,2)).has_value()); - - // belove ray - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(2,-1)).has_value()); - CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,-1), Vec2d(2,0)).has_value()); - - // intersection at [2,1] distance 1 - auto t1 = ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(2,2)); - REQUIRE(t1.has_value()); - auto t2 = ray_segment_intersection(r_point, r_dir, Vec2d(2,2), Vec2d(2,0)); - REQUIRE(t2.has_value()); - - CHECK(abs(*t1 - *t2) < std::numeric_limits::epsilon()); -} - -Vec2d get_intersection(const Vec2d& point, const Vec2d& dir, const std::array& triangle) -{ - std::optional t; - for (size_t i = 0; i < 3; ++i) { - size_t i2 = i+1; - if(i2 == 3) i2 = 0; - if (!t.has_value()) { - t = ray_segment_intersection(point, dir, triangle[i], triangle[i2]); - continue; - } - - // small distance could be preccission inconsistance - std::optional t2 = ray_segment_intersection( - point, dir, triangle[i], triangle[i2]); - if (t2.has_value() && *t2 > *t) t = t2; - - } - assert(t.has_value()); //Not found intersection. - return point + dir * (*t); -} - -TEST_CASE("triangle intersection", "[]") -{ - Vec2d point(1,1); - Vec2d dir(-1, 0); - std::array triangle = { - Vec2d(0, 0), - Vec2d(5, 0), - Vec2d(0, 5) - }; - Vec2d i = get_intersection(point, dir, triangle); - CHECK(abs(i.x()) < std::numeric_limits::epsilon()); - CHECK(abs(i.y() - 1.) < std::numeric_limits::epsilon()); -} - -#include -#include - -indexed_triangle_set emboss3d(const Polygon &shape, float height) { - // CW order of triangle indices - std::vector shape_triangles = Triangulation::triangulate(shape); - - indexed_triangle_set result; - const Points &pts = shape.points; - size_t count_point = pts.size(); - result.vertices.reserve(2 * count_point); - // top points - std::transform(pts.begin(), pts.end(), std::back_inserter(result.vertices), - [](const Point &p) { return Vec3f(p.x(), p.y(), 0); }); - // bottom points - std::transform(pts.begin(), pts.end(), std::back_inserter(result.vertices), - [height](const Point &p) { return Vec3f(p.x(), p.y(), height); }); - - result.indices.reserve(shape_triangles.size() * 2 + count_point*2); - // top triangles - change to CCW - for (const Vec3i &t : shape_triangles) - result.indices.emplace_back(t.x(), t.z(), t.y()); - // bottom triangles - use CW - for (const Vec3i &t : shape_triangles) - result.indices.emplace_back( - t.x() + count_point, - t.y() + count_point, - t.z() + count_point - ); - - // quads around - zig zag by triangles - for (uint32_t i = 0; i < count_point; ++i) { - // previous index - uint32_t ip = (i == 0) ? (count_point - 1) : (i - 1); - // bottom indices - uint32_t i2 = i + count_point; - uint32_t ip2 = ip + count_point; - - result.indices.emplace_back(i, i2, ip); - result.indices.emplace_back(ip2, ip, i2); - } - return result; -} - -#include -#include -#include -#include - -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Exact_predicates_exact_constructions_kernel EK; -typedef CGAL::Surface_mesh Mesh; -namespace PMP = CGAL::Polygon_mesh_processing; -namespace params = PMP::parameters; - -// is it realy neccessary to copy all? -Mesh its_to_mesh(const indexed_triangle_set &its) -{ - Mesh mesh; - size_t vertices_count = its.vertices.size(); - size_t edges_count = (its.indices.size() * 3) / 2; - size_t faces_count = its.indices.size(); - mesh.reserve(vertices_count, edges_count, faces_count); - - for (auto &v : its.vertices) - mesh.add_vertex(typename Mesh::Point{v.x(), v.y(), v.z()}); - - using VI = typename Mesh::Vertex_index; - for (auto &f : its.indices) - mesh.add_face(VI(f(0)), VI(f(1)), VI(f(2))); - - return mesh; -} - -indexed_triangle_set mesh_to_its(const Mesh &mesh) -{ - indexed_triangle_set res; - res.vertices.reserve(mesh.num_vertices()); - res.indices.reserve(mesh.num_faces()); - - for (auto &vi : mesh.vertices()) { - auto &v = mesh.point(vi); // index addresing, to compare same float point - res.vertices.emplace_back(v.x(),v.y(), v.z()); - } - - for (auto &face : mesh.faces()) { - auto vtc = mesh.vertices_around_face(mesh.halfedge(face)); - int i = 0; - Vec3i facet; - for (auto v : vtc) { - if (i > 2) { - i = 0; - break; - } - facet(i++) = v; - } - if (i == 3) res.indices.emplace_back(facet); - } - return res; -} - -void emboss3d_(const Polygon& shape, const EmbossConfig &cfg, indexed_triangle_set &its) -{ - indexed_triangle_set shape3d = emboss3d(shape, cfg.height); - - //its_write_obj(shape3d,"shape3d.obj"); - - Mesh m1 = its_to_mesh(its); - Mesh m2 = its_to_mesh(shape3d); - - size_t in_vertices_count = its.vertices.size(); - size_t in_indices_count = its.indices.size(); - - auto tt = m1.property_map("e:removed"); - Mesh::Property_map ecm1 = tt.first; // hope in copy - - PMP::corefine(m1, m2 - , PMP::parameters::edge_is_constrained_map(ecm1) - //, PMP::parameters::do_not_modify(true) - ); - //PMP::Corefinement::Intersection_of_triangle_meshes() - - size_t count_true = 0; - for (const Mesh::Edge_index &e : edges(m1)) - if (ecm1[e]) { - ++count_true; - const Mesh::Halfedge_index &h = e.halfedge(); - const Mesh::Halfedge_index &h2 = m1.opposite(h); - const Mesh::Vertex_index &v = m1.target(h); - std::cout << - "edge " << e << - " halfedge " << h << - " target " << v << - " index0 " << v.idx() << - " index1 " << m1.target(h2).idx() << - std::endl; - } - - - its = mesh_to_its(m1); - - // twice and move additional vertices - Vec3f move = (cfg.projection.normal * cfg.height).cast(); - size_t new_add_vertice = its.vertices.size() - in_vertices_count; - its.vertices.reserve(its.vertices.size() + new_add_vertice); - for (size_t i = in_vertices_count; i < its.vertices.size(); ++i) - its.vertices.emplace_back(its.vertices[i] + move); - - // zig zag edge collected edge - for (size_t i = in_indices_count; i < its.indices.size(); ++i) { - // new added triangle - Vec3i &t = its.indices[i]; - - std::array is_new; - bool is_inner = true; - for (size_t i = 0; i < 3; ++i) { - bool is_n = t[0] >= in_vertices_count; - is_inner |= is_n; - is_new[i] = is_n; - } - - } - - its_write_obj(mesh_to_its(m1), "corefine1.obj"); - its_write_obj(mesh_to_its(m2), "corefine2.obj"); - - /*MeshBoolean::cgal::CGALMesh cm1 = m1->m; - - My_visitor sm_v; - - CGAL::Polygon_mesh_processing::corefine(m1->m, m2->m, - CGAL::Polygon_mesh_processing::parameters::visitor(sm_v) - , CGAL::parameters::do_not_modify(true) - );*/ - - //TriangleMesh tm1 = cgal_to_triangle_mesh(*m1); - //store_obj("tm1.obj", &tm1); - //TriangleMesh tm2 = cgal_to_triangle_mesh(*m2); - //store_obj("tm2.obj", &tm1); - - its = mesh_to_its(m1); - return; -} - -void emboss3d(const Polygon& shape, const EmbossConfig &cfg, indexed_triangle_set &its) -{ - indexed_triangle_set shape3d = emboss3d(shape, abs(cfg.height)); - //its_write_obj(shape3d,"shape3d.obj"); - auto m1 = MeshBoolean::cgal::triangle_mesh_to_cgal(its); - auto m2 = MeshBoolean::cgal::triangle_mesh_to_cgal(shape3d); - - bool is_plus = shape.is_counter_clockwise() == (cfg.height > 0.f); - if (is_plus) { - MeshBoolean::cgal::plus(*m1, *m2); - } else { - MeshBoolean::cgal::minus(*m1, *m2); - } - - TriangleMesh tm1 = cgal_to_triangle_mesh(*m1); - store_obj("tm1.obj", &tm1); - //TriangleMesh tm2 = cgal_to_triangle_mesh(*m2); - //store_obj("tm2.obj", &tm1); - - its = tm1.its; // copy -} - -TEST_CASE("Emboss polygon example", "[MeshBoolean]") -{ - const char *font_name = "C:/windows/fonts/arialbd.ttf"; - char letter = '%'; - float flatness = 2.; - ExPolygons espolygons = ttf2polygons(font_name, letter, flatness); - - //TriangleMesh tm = make_sphere(1., 1.); tm.scale(10.f); - - TriangleMesh tm = make_cube(10., 5., 2.); - tm.translate(Vec3f(0, 0, 1.7)); - Polygons polygons; - polygons = {Polygon({{1, 1}, {1, 2}, {2, 2}, {2, 1}})}; // rectangle CW - polygons = {Polygon({{1, 1}, {2, 1}, {2, 2}, {1, 2}})}; // rectangle CCW - - EmbossConfig ec; - ec.height = 3; - indexed_triangle_set its = tm.its; // copy - for (const auto& polygon : polygons) { - // TODO: differ CW and CCW - emboss3d_(polygon, ec, its); - } -} - -TEST_CASE("Triangulate by cgal", "[Triangulation]") -{ - Points points = {Point(1, 1), Point(2, 1), Point(2, 2), Point(1, 2)}; - Triangulation::HalfEdges edges1 = {{1, 3}}; - std::vector indices1 = Triangulation::triangulate(points, edges1); - - auto check = [](int i1, int i2, Vec3i t)->bool { return true; - return (t[0] == i1 || t[1] == i1 || t[2] == i1) && - (t[0] == i2 || t[1] == i2 || t[2] == i2); - }; - REQUIRE(indices1.size() == 2); - int i1 = edges1.begin()->first, - i2 = edges1.begin()->second; - CHECK(check(i1, i2, indices1[0])); - CHECK(check(i1, i2, indices1[1])); - - Triangulation::HalfEdges edges2 = {{0, 2}}; - std::vector indices2 = Triangulation::triangulate(points, edges2); - REQUIRE(indices2.size() == 2); - i1 = edges2.begin()->first; - i2 = edges2.begin()->second; - CHECK(check(i1, i2, indices2[0])); - CHECK(check(i1, i2, indices2[1])); -} diff --git a/tests/libslic3r/test_quadric_edge_collapse.cpp b/tests/libslic3r/test_quadric_edge_collapse.cpp new file mode 100644 index 0000000000..ef71c98b0b --- /dev/null +++ b/tests/libslic3r/test_quadric_edge_collapse.cpp @@ -0,0 +1,242 @@ +#include +#include + +#include +#include // its - indexed_triangle_set +#include // no priority queue + +#include "libslic3r/AABBTreeIndirect.hpp" // is similar + +using namespace Slic3r; + +namespace Private { + +struct Similarity +{ + float max_distance = 0.f; + float average_distance = 0.f; + + Similarity() = default; + Similarity(float max_distance, float average_distance) + : max_distance(max_distance), average_distance(average_distance) + {} +}; + +// border for our algorithm with frog_leg model and decimation to 5% +Similarity frog_leg_5(0.32f, 0.043f); + +Similarity get_similarity(const indexed_triangle_set &from, + const indexed_triangle_set &to) +{ + // create ABBTree + auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set( + from.vertices, from.indices); + float sum_distance = 0.f; + + float max_distance = 0.f; + auto collect_distances = [&](const Vec3f &surface_point) { + size_t hit_idx; + Vec3f hit_point; + float distance2 = + AABBTreeIndirect::squared_distance_to_indexed_triangle_set( + from.vertices, from.indices, tree, surface_point, hit_idx, + hit_point); + float distance = sqrt(distance2); + if (max_distance < distance) max_distance = distance; + sum_distance += distance; + }; + + for (const Vec3f &vertex : to.vertices) { collect_distances(vertex); } + for (const Vec3i &t : to.indices) { + Vec3f center(0, 0, 0); + for (size_t i = 0; i < 3; ++i) { center += to.vertices[t[i]] / 3; } + collect_distances(center); + } + + size_t count = to.vertices.size() + to.indices.size(); + float average_distance = sum_distance / count; + + std::cout << "max_distance = " << max_distance << ", average_distance = " << average_distance << std::endl; + return Similarity(max_distance, average_distance); +} + +void is_better_similarity(const indexed_triangle_set &its_first, + const indexed_triangle_set &its_second, + const Similarity & compare) +{ + Similarity s1 = get_similarity(its_first, its_second); + Similarity s2 = get_similarity(its_second, its_first); + + CHECK(s1.average_distance < compare.average_distance); + CHECK(s1.max_distance < compare.max_distance); + CHECK(s2.average_distance < compare.average_distance); + CHECK(s2.max_distance < compare.max_distance); +} + +void is_worse_similarity(const indexed_triangle_set &its_first, + const indexed_triangle_set &its_second, + const Similarity & compare) +{ + Similarity s1 = get_similarity(its_first, its_second); + Similarity s2 = get_similarity(its_second, its_first); + + if (s1.max_distance < compare.max_distance && + s2.max_distance < compare.max_distance) + CHECK(false); +} + +bool exist_triangle_with_twice_vertices(const std::vector &indices) +{ + for (const auto &face : indices) + if (face[0] == face[1] || face[0] == face[2] || face[1] == face[2]) + return true; + return false; +} + +} // namespace Private + +TEST_CASE("Reduce one edge by Quadric Edge Collapse", "[its]") +{ + indexed_triangle_set its; + its.vertices = {Vec3f(-1.f, 0.f, 0.f), Vec3f(0.f, 1.f, 0.f), + Vec3f(1.f, 0.f, 0.f), Vec3f(0.f, 0.f, 1.f), + // vertex to be removed + Vec3f(0.9f, .1f, -.1f)}; + its.indices = {Vec3i(1, 0, 3), Vec3i(2, 1, 3), Vec3i(0, 2, 3), + Vec3i(0, 1, 4), Vec3i(1, 2, 4), Vec3i(2, 0, 4)}; + // edge to remove is between vertices 2 and 4 on trinagles 4 and 5 + + indexed_triangle_set its_ = its; // copy + // its_write_obj(its, "tetrhedron_in.obj"); + uint32_t wanted_count = its.indices.size() - 1; + its_quadric_edge_collapse(its, wanted_count); + // its_write_obj(its, "tetrhedron_out.obj"); + CHECK(its.indices.size() == 4); + CHECK(its.vertices.size() == 4); + + for (size_t i = 0; i < 3; i++) { + CHECK(its.indices[i] == its_.indices[i]); + } + + for (size_t i = 0; i < 4; i++) { + if (i == 2) continue; + CHECK(its.vertices[i] == its_.vertices[i]); + } + + const Vec3f &v = its.vertices[2]; // new vertex + const Vec3f &v2 = its_.vertices[2]; // moved vertex + const Vec3f &v4 = its_.vertices[4]; // removed vertex + for (size_t i = 0; i < 3; i++) { + bool is_between = (v[i] < v4[i] && v[i] > v2[i]) || + (v[i] > v4[i] && v[i] < v2[i]); + CHECK(is_between); + } + Private::Similarity max_similarity(0.75f, 0.014f); + Private::is_better_similarity(its, its_, max_similarity); +} + +TEST_CASE("Simplify frog_legs.obj to 5% by Quadric edge collapse", "[its][quadric_edge_collapse]") +{ + TriangleMesh mesh = load_model("frog_legs.obj"); + double original_volume = its_volume(mesh.its); + uint32_t wanted_count = mesh.its.indices.size() * 0.05; + REQUIRE_FALSE(mesh.empty()); + indexed_triangle_set its = mesh.its; // copy + float max_error = std::numeric_limits::max(); + its_quadric_edge_collapse(its, wanted_count, &max_error); + // its_write_obj(its, "frog_legs_qec.obj"); + CHECK(its.indices.size() <= wanted_count); + double volume = its_volume(its); + CHECK(fabs(original_volume - volume) < 33.); + + Private::is_better_similarity(mesh.its, its, Private::frog_leg_5); +} + +#include +#include "Simplify.h" +TEST_CASE("Simplify frog_legs.obj to 5% by IGL/qslim", "[]") +{ + std::string obj_filename = "frog_legs.obj"; + TriangleMesh mesh = load_model(obj_filename); + REQUIRE_FALSE(mesh.empty()); + indexed_triangle_set &its = mesh.its; + double original_volume = its_volume(its); + uint32_t wanted_count = its.indices.size() * 0.05; + + Eigen::MatrixXd V(its.vertices.size(), 3); + Eigen::MatrixXi F(its.indices.size(), 3); + for (size_t j = 0; j < its.vertices.size(); ++j) { + Vec3d vd = its.vertices[j].cast(); + for (int i = 0; i < 3; ++i) V(j, i) = vd(i); + } + + for (size_t j = 0; j < its.indices.size(); ++j) { + const auto &f = its.indices[j]; + for (int i = 0; i < 3; ++i) F(j, i) = f(i); + } + + size_t max_m = wanted_count; + Eigen::MatrixXd U; + Eigen::MatrixXi G; + Eigen::VectorXi J, I; + CHECK(igl::qslim(V, F, max_m, U, G, J, I)); + + // convert to its + indexed_triangle_set its_out; + its_out.vertices.reserve(U.size()/3); + its_out.indices.reserve(G.size()/3); + for (size_t i = 0; i < U.size()/3; i++) + its_out.vertices.emplace_back(U(i, 0), U(i, 1), U(i, 2)); + for (size_t i = 0; i < G.size()/3; i++) + its_out.indices.emplace_back(G(i, 0), G(i, 1), G(i, 2)); + + // check if algorithm is still worse than our + Private::is_worse_similarity(its_out, its, Private::frog_leg_5); + // its_out, its --> avg_distance: 0.0351217, max_distance 0.364316 + // its, its_out --> avg_distance: 0.0412358, max_distance 0.238913 +} + +TEST_CASE("Simplify frog_legs.obj to 5% by simplify", "[]") { + std::string obj_filename = "frog_legs.obj"; + TriangleMesh mesh = load_model(obj_filename); + uint32_t wanted_count = mesh.its.indices.size() * 0.05; + Simplify::load_obj((TEST_DATA_DIR PATH_SEPARATOR + obj_filename).c_str()); + Simplify::simplify_mesh(wanted_count, 5, true); + + // convert to its + indexed_triangle_set its_out; + its_out.vertices.reserve(Simplify::vertices.size()); + its_out.indices.reserve(Simplify::triangles.size()); + for (size_t i = 0; i < Simplify::vertices.size(); i++) { + const Simplify::Vertex &v = Simplify::vertices[i]; + its_out.vertices.emplace_back(v.p.x, v.p.y, v.p.z); + } + for (size_t i = 0; i < Simplify::triangles.size(); i++) { + const Simplify::Triangle &t = Simplify::triangles[i]; + its_out.indices.emplace_back(t.v[0], t.v[1], t.v[2]); + } + + // check if algorithm is still worse than our + Private::is_worse_similarity(its_out, mesh.its, Private::frog_leg_5); + // its_out, mesh.its --> max_distance = 0.700494, average_distance = 0.0902524 + // mesh.its, its_out --> max_distance = 0.393184, average_distance = 0.0537392 +} + +TEST_CASE("Simplify trouble case", "[its]") +{ + TriangleMesh tm = load_model("simplification.obj"); + REQUIRE_FALSE(tm.empty()); + float max_error = std::numeric_limits::max(); + uint32_t wanted_count = 0; + its_quadric_edge_collapse(tm.its, wanted_count, &max_error); + CHECK(!Private::exist_triangle_with_twice_vertices(tm.its.indices)); +} + +TEST_CASE("Simplified cube should not be empty.", "[its]") +{ + auto its = its_make_cube(1, 2, 3); + float max_error = std::numeric_limits::max(); + uint32_t wanted_count = 0; + its_quadric_edge_collapse(its, wanted_count, &max_error); + CHECK(!its.indices.empty()); +} diff --git a/tests/libslic3r/test_triangulation.cpp b/tests/libslic3r/test_triangulation.cpp new file mode 100644 index 0000000000..bc0f371d74 --- /dev/null +++ b/tests/libslic3r/test_triangulation.cpp @@ -0,0 +1,136 @@ +#include + +#include +#include // only debug visualization + +using namespace Slic3r; + +namespace Private{ +void store_trinagulation(const ExPolygons & shape, + const std::vector &triangles, + const char* file_name = "Triangulation_visualization.svg", + double scale = 1e5) +{ + BoundingBox bb; + for (const auto &expoly : shape) bb.merge(expoly.contour.points); + bb.scale(scale); + SVG svg_vis(file_name, bb); + svg_vis.draw(shape, "gray", .7f); + + size_t count = count_points(shape); + Points points; + points.reserve(count); + auto insert_point = [](const Slic3r::Polygon &polygon, Points &points) { + for (const Point &p : polygon.points) points.emplace_back(p); + }; + for (const ExPolygon &expolygon : shape) { + insert_point(expolygon.contour, points); + for (const Slic3r::Polygon &hole : expolygon.holes) + insert_point(hole, points); + } + + for (const auto &t : triangles) { + Polygon triangle({points[t[0]], points[t[1]], points[t[2]]}); + triangle.scale(scale); + svg_vis.draw(triangle, "green"); + } + + // prevent visualization in test + CHECK(false); +} +} // namespace + + +TEST_CASE("Triangulate rectangle with restriction on edge", "[Triangulation]") +{ + // 0 1 2 3 + Points points = {Point(1, 1), Point(2, 1), Point(2, 2), Point(1, 2)}; + Triangulation::HalfEdges edges1 = {{1, 3}}; + std::vector indices1 = Triangulation::triangulate(points, edges1, true); + + auto check = [](int i1, int i2, Vec3i t) -> bool { + return true; + return (t[0] == i1 || t[1] == i1 || t[2] == i1) && + (t[0] == i2 || t[1] == i2 || t[2] == i2); + }; + REQUIRE(indices1.size() == 2); + int i1 = edges1.begin()->first, i2 = edges1.begin()->second; + CHECK(check(i1, i2, indices1[0])); + CHECK(check(i1, i2, indices1[1])); + + Triangulation::HalfEdges edges2 = {{0, 2}}; + std::vector indices2 = Triangulation::triangulate(points, edges2, true); + REQUIRE(indices2.size() == 2); + i1 = edges2.begin()->first; + i2 = edges2.begin()->second; + CHECK(check(i1, i2, indices2[0])); + CHECK(check(i1, i2, indices2[1])); +} + +TEST_CASE("Triangulation polygon", "[triangulation]") +{ + Points points = {Point(416, 346), Point(445, 362), Point(463, 389), + Point(469, 427), Point(445, 491)}; + + Polygon polygon(points); + Polygons polygons({polygon}); + ExPolygon expolygon(points); + ExPolygons expolygons({expolygon}); + + std::vector tp = Triangulation::triangulate(polygon); + std::vector tps = Triangulation::triangulate(polygons); + std::vector tep = Triangulation::triangulate(expolygon); + std::vector teps = Triangulation::triangulate(expolygons); + + //Private::store_trinagulation(expolygons, teps); + + CHECK(tp.size() == tps.size()); + CHECK(tep.size() == teps.size()); + CHECK(tp.size() == tep.size()); + CHECK(tp.size() == 3); +} + +TEST_CASE("Triangulation M shape polygon", "[triangulation]") +{ + // 0 1 2 3 4 + Polygon shape_M = {Point(0, 0), Point(2, 0), Point(2, 2), Point(1, 1), Point(0, 2)}; + + std::vector triangles = Triangulation::triangulate(shape_M); + + // Check outer triangle is not contain + std::set outer_triangle = {2, 3, 4}; + bool is_in = false; + for (const Vec3i &t : triangles) { + for (size_t i = 0; i < 3; i++) { + int index = t[i]; + if (outer_triangle.find(index) == outer_triangle.end()) { + is_in = false; + break; + } else { + is_in = true; + } + } + if (is_in) break; + } + + //Private::store_trinagulation({ExPolygon(shape_M)}, triangles); + CHECK(triangles.size() == 3); + CHECK(!is_in); +} + +// same point in triangulation are not Supported +//TEST_CASE("Triangulation 2 polygons with same point", "[triangulation][not supported]") +//{ +// Slic3r::Polygon polygon1 = { +// Point(416, 346), Point(445, 362), +// Point(463, 389), Point(469, 427) /* This point */, +// Point(445, 491) +// }; +// Slic3r::Polygon polygon2 = { +// Point(495, 488), Point(469, 427) /* This point */, +// Point(495, 364) +// }; +// ExPolygons shape2d = {ExPolygon(polygon1), ExPolygon(polygon2)}; +// std::vector shape_triangles = Triangulation::triangulate(shape2d); +// store_trinagulation(shape2d, shape_triangles); +//}