From 8ae1fbf6127650e5dc2b5e9a4ed3c857366a461b Mon Sep 17 00:00:00 2001 From: Trammell hudson Date: Sat, 3 Mar 2018 14:14:06 -0500 Subject: [PATCH] hidden wireframe works, change v3 api a little bit --- camera.c | 78 +++--- hiddenwire.c | 35 +-- seg.h | 12 +- stl_3d.c | 8 +- svg.h | 41 +++ test1.scad | 8 +- test1.stl | Bin 684 -> 5795 bytes tri.c | 685 ++++++++++++++++++++++++++++++--------------------- tri.h | 40 ++- unfold.c | 8 +- v3.h | 96 +++++++- wireframe.c | 6 +- 12 files changed, 637 insertions(+), 380 deletions(-) create mode 100644 svg.h diff --git a/camera.c b/camera.c index d0e0469..0e17873 100644 --- a/camera.c +++ b/camera.c @@ -13,7 +13,7 @@ struct _camera_t { float near; float far; - float r[4][4]; + m44_t r; }; @@ -53,7 +53,7 @@ camera_setup( // and the "up" normal v3_t v = v3_norm(v3_cross(w, u)); - float cam[4][4] = { + m44_t cam = {{ #if 0 { u.p[0], v.p[0], w.p[0], 0 }, { u.p[1], v.p[1], w.p[1], 0 }, @@ -65,59 +65,62 @@ camera_setup( { w.p[0], w.p[1], w.p[2], -v3_dot(w,eye) }, { 0, 0, 0, 1 }, #endif - }; + }}; fprintf(stderr, "Camera:\n"); for(int i = 0 ; i < 4 ; i++) { for(int j = 0 ; j < 4 ; j++) - fprintf(stderr, " %+5.3f", cam[i][j]); + fprintf(stderr, " %+5.3f", cam.m[i][j]); fprintf(stderr, "\n"); } // now compute the perspective projection matrix +if(0) { float s = 1.0 / tan(fov * M_PI / 180 / 2); - c->near = 4.0; - c->far = 200; + c->near = 1; + c->far = 2; float f1 = - c->far / (c->far - c->near); float f2 = - c->far * c->near / (c->far - c->near); - float pers[4][4] = { + m44_t pers = {{ { s, 0, 0, 0 }, { 0, s, 0, 0 }, - { 0, 0, f1, -1 }, - { 0, 0, f2, 0 }, - }; + { 0, 0, f2, -1 }, + { 0, 0, f1, 0 }, + }}; fprintf(stderr, "Perspective:\n"); for(int i = 0 ; i < 4 ; i++) { for(int j = 0 ; j < 4 ; j++) - fprintf(stderr, " %+5.3f", pers[i][j]); + fprintf(stderr, " %+5.3f", pers.m[i][j]); fprintf(stderr, "\n"); } - // and apply it to the camera matrix to generate transform - for(int i = 0 ; i < 4 ; i++) - { - for(int j = 0 ; j < 4 ; j++) - { - float d = 0; - for(int k = 0 ; k < 4 ; k++) - d += pers[i][k] * cam[k][j]; - c->r[i][j] = d; - } - } + m44_mult(&c->r, &cam, &pers); +} else { + // no perspective + m44_t pers = {{ + { 1, 0, 0, 0 }, + { 0, 1, 0, 0 }, + { 0, 0, 1, 0 }, + { 0, 0, 0, 1 }, + }}; + // and apply it to the camera matrix to generate transform + m44_mult(&c->r, &cam, &pers); +} fprintf(stderr, "Cam*Pers\n"); for(int i = 0 ; i < 4 ; i++) { for(int j = 0 ; j < 4 ; j++) - fprintf(stderr, " %+5.3f", c->r[i][j]); + fprintf(stderr, " %+5.3f", c->r.m[i][j]); fprintf(stderr, "\n"); } + } @@ -133,32 +136,29 @@ camera_project( v3_t * const v_out ) { - float v[4] = { v_in->p[0], v_in->p[1], v_in->p[2], 1 }; - float p[4] = { 0, 0, 0, 0}; + v4_t v = {{ v_in->p[0], v_in->p[1], v_in->p[2], 1 }}; + v4_t p = m44_multv(&c->r, &v); - for (int i = 0 ; i < 4 ; i++) - for (int j = 0 ; j < 4 ; j++) - p[i] += c->r[i][j] * v[j]; + p.p[2] *= -1; // what if p->p[4] == 0? // pz < 0 == The point is behind us; do not display? //if (p[2] < c->near || p[2] > c->far) - if (p[2] < 0) + if (p.p[2] <= 0) return 0; - for (int i = 0 ; i < 3 ; i++) - p[i] /= p[3]; - - if(0) fprintf(stderr, "%.2f %.2f %.2f -> %.2f %.2f %.2f %.2f\n", - v[0], v[1], v[2], - p[0], p[1], p[2], p[3] - ); + // shrink by the distance + p.p[0] *= 1000 / p.p[2]; + p.p[1] *= 1000 / p.p[2]; + //p[2] /= 1000; // Transform to screen coordinate frame, // and return it to the caller - v_out->p[0] = p[0]; - v_out->p[1] = p[1]; - v_out->p[2] = p[2]; + v_out->p[0] = p.p[0]; + v_out->p[1] = p.p[1]; + v_out->p[2] = p.p[2]; + + v3_print(*v_out); return 1; } diff --git a/hiddenwire.c b/hiddenwire.c index 68b3921..35c36b0 100644 --- a/hiddenwire.c +++ b/hiddenwire.c @@ -18,6 +18,7 @@ static int debug = 0; #include "v3.h" #include "tri.h" #include "camera.h" +#include "svg.h" static const char usage[] = @@ -77,25 +78,6 @@ stl_face_t; -void -svg_line( - const char * color, - const float * p1, - const float * p2, - float thick -) -{ - // invert the sense of y - printf("\n", - p1[0], - -p1[1], - p2[0], - -p2[1], - color, - thick - ); -} - static inline int v2_eq( @@ -302,11 +284,6 @@ int main( behind++; goto reject_early; } - - // scale to the image size - s[j].p[0] *= width; - s[j].p[1] *= width; - s[j].p[2] *= width; } if(debug >= 2) @@ -326,7 +303,7 @@ int main( tri_t * const tri = tri_new(s, stl->p); // reject this face if any of the vertices are behind us - if (tri->min[2] < 0) + if (tri->min.p[2] < 0) { behind++; goto reject; @@ -342,6 +319,7 @@ int main( } // if it has any off-screen coords, reject it +/* if (!onscreen(&tri->p[0], width, height) || !onscreen(&tri->p[1], width, height) || !onscreen(&tri->p[2], width, height)) @@ -350,6 +328,7 @@ tri_print(tri); offscreen++; goto reject; } +*/ // prune the small triangles in the screen space if (tri_area_2d(tri) < prune) @@ -378,6 +357,8 @@ reject: reject_early: continue; } + for(tri_t * t = zlist ; t ; t = t->next) + tri_print(t); if (debug) fprintf(stderr, "Retained %d triangles, rejected %d behind, %d offscreen, %d backface, %d small\n", retained, behind, offscreen, backface, small_area); @@ -385,6 +366,7 @@ reject_early: // drop any triangles that are totally occluded by another // triangle. this reduces the amount of work for later rejected = 0; +#if 0 for(tri_t * t = zlist ; t ; t = t->next) { tri_t * t2_next; @@ -405,6 +387,7 @@ reject_early: } if (debug) fprintf(stderr, "Rejected %d fully occluded triangles\n", rejected); +#endif // generate a list of segments, dropping any coplanar ones @@ -469,7 +452,7 @@ reject_early: seg_t * s = slist; slist = s->next; - tri_seg_intersect(zlist, s, &slist_visible); + tri_seg_hidden(zlist, s, &slist_visible); } } else { // don't do any intersection tests diff --git a/seg.h b/seg.h index 042eb4f..d013b9a 100644 --- a/seg.h +++ b/seg.h @@ -34,15 +34,17 @@ seg_print( const seg_t * const s ) { - fprintf(stderr, "%.0f,%.0f -> %.0f,%.0f (was %.0f,%.0f -> %.0f,%.0f)\n", + fprintf(stderr, "%+5.1f,%+5.1f,%5.1f -> %+5.1f,%+5.1f,%5.1f\n", s->p[0].p[0], s->p[0].p[1], + s->p[0].p[2], s->p[1].p[0], s->p[1].p[1], - s->src[0].p[0], - s->src[0].p[1], - s->src[1].p[0], - s->src[1].p[1] + s->p[1].p[2] + //s->src[0].p[0], + //s->src[0].p[1], + //s->src[1].p[0], + //s->src[1].p[1] ); } diff --git a/stl_3d.c b/stl_3d.c index f9416e3..b9aa550 100644 --- a/stl_3d.c +++ b/stl_3d.c @@ -38,7 +38,7 @@ stl_vertex_find( { stl_vertex_t * const v = &vertices[x]; - if (v3_eq(&v->p, p)) + if (v3_eq(v->p, *p)) return v; } @@ -95,11 +95,11 @@ stl_angle( for (int i = 0 ; i < 3 ; i++) { x4 = f2->vertex[i]->p; - if (v3_eq(&x1, &x4)) + if (v3_eq(x1, x4)) continue; - if (v3_eq(&x2, &x4)) + if (v3_eq(x2, x4)) continue; - if (v3_eq(&x3, &x4)) + if (v3_eq(x3, x4)) continue; break; } diff --git a/svg.h b/svg.h new file mode 100644 index 0000000..2750f38 --- /dev/null +++ b/svg.h @@ -0,0 +1,41 @@ +#ifndef _svg_h_ +#define _svg_h_ + +static inline void +svg_line( + const char * color, + const float * p1, + const float * p2, + float thick +) +{ + // invert the sense of y + printf("\n", + p1[0], + -p1[1], + p2[0], + -p2[1], + color, + thick + ); +} + +static inline void +svg_circle( + const char * color, + const float p1, + const float p2, + float radius +) +{ + // invert the sense of y + printf("\n", + p1, + -p2, + radius, + color, + 1.0 + ); +} + +#endif diff --git a/test1.scad b/test1.scad index 904e475..52ab31f 100644 --- a/test1.scad +++ b/test1.scad @@ -1 +1,7 @@ -cube([10,20,30], center=true); +render() difference() +{ + cube([10,20,25], center=true); + cube([20,15,20], center=true); +} + +translate([3,2,4]) cube([25,8,8], center=true); diff --git a/test1.stl b/test1.stl index a017da2edb6a17c95861912088004677a7a593ea..536ed94b483c704255326ba8c8dfeb791a19bf11 100644 GIT binary patch literal 5795 zcmcIo%Wm5+5WMFr_6HPL0hBh!pt%M;^aF~ZR)GMH6*x}O&o7@VlIC)UhI&g@lxB8j zxT5v39A5YI{e3?E`1bhp^Lg3NhZ~~bFTdt99hcMF%YlZGcC<59qUG~EpXjhG@9tls zKlAB4|D|D~ac36Z-A)%5zT8gciod_Of!EK+eR!z(%Yko9Vvty*b%@0PXd$8*9R!u! zEkPxS2STHbCYZ9ixD^kpu@ZnZ7~Pg2Qe44dYc~W6HWR!NQXSIzR#u2$Y5r0Yd13(? z#ta&PVR6~A8xs{N+-%|HIt879gyhE^ttz9#Zep>)nY;*PZ zW4)R#-YOS0B#D=)C(A(w>ISioG{(NUO+e?fs{z3 zrZc%bI2Z+eYg!7Z5!c9o%!jHBS?*Yb-Z)~Dp(cWq1rNAooW;M-WN-SMwo&kU>=uy26n zJR{dV7ceDQ83@rvxWn-)D&AV()c7eE!Br!OQt2if2%jue^f{&CsYxPdw$)Do^19e>giZN5o0 KC0Lo#yuSh15^V_p diff --git a/tri.c b/tri.c index c43a7bc..0d8cd6e 100644 --- a/tri.c +++ b/tri.c @@ -34,11 +34,8 @@ tri_new( // compute the bounding box for the triangle in camera space - for(int j = 0 ; j < 3 ; j++) - { - t->min[j] = min(min(t->p[0].p[j], t->p[1].p[j]), t->p[2].p[j]); - t->max[j] = max(max(t->p[0].p[j], t->p[1].p[j]), t->p[2].p[j]); - } + t->min = v3_min(v3_min(t->p[0], t->p[1]), t->p[2]); + t->max = v3_max(v3_max(t->p[0], t->p[1]), t->p[2]); return t; } @@ -59,7 +56,7 @@ tri_insert( // check to see if our new triangle is closer than // the current triangle - if(iter->min[2] > t->min[2]) + if(iter->min.p[2] > t->min.p[2]) break; zlist = &(iter->next); @@ -110,7 +107,7 @@ tri_print( const tri_t * const t ) { - fprintf(stderr, "%.0f,%.0f,%.0f %.0f,%.0f,%.0f %.0f,%.0f,%.0f norm %.3f,%.3f,%.3f\n", + fprintf(stderr, "{{{%+5.1f,%+5.1f,%5.1f }},{{%+5.1f,%+5.1f,%5.1f }},{{%+5.1f,%+5.1f,%5.1f }}}\n", // norm %.3f %.3f %.3f\n", t->p[0].p[0], t->p[0].p[1], t->p[0].p[2], @@ -119,10 +116,10 @@ tri_print( t->p[1].p[2], t->p[2].p[0], t->p[2].p[1], - t->p[2].p[2], - t->normal.p[0], - t->normal.p[1], - t->normal.p[2] + t->p[2].p[2] + //t->normal.p[0], + //t->normal.p[1], + //t->normal.p[2] ); } @@ -149,7 +146,7 @@ tri_coplanar( { for(int j = 0 ; j < 3 ; j++) { - if (!v3_eq(&t0->p[i], &t1->p[j])) + if (!v3_eq(t0->p[i], t1->p[j])) continue; matches |= 1 << i; break; @@ -162,10 +159,14 @@ tri_coplanar( case 0x6: return 1; case 0x5: return 2; case 0x7: - fprintf(stderr, "uh, three points match?\n"); - tri_print(t0); - tri_print(t1); - return -1; + // these are likely small triangles that can be ignored + if (tri_debug > 3) + { + fprintf(stderr, "uh, three points match?\n"); + tri_print(t0); + tri_print(t1); + } + return 0; default: // no shared edge return -1; @@ -212,6 +213,60 @@ tri_find_z( } +/* + * Find the barycentric coordinates for point of an XY coordinate in a triangle. + * + * p can be written as a combination of t01 and t02, + * p - t0 = a * (t1 - t0) + b * (t2 - t0) + * setting t0 to 0, this becomes: + * p = a * t1 + b * t2 + * which is two equations with two unknowns + * + * The x and y coordinates are based on the two sides of the triangle + * and the z coordinate is the screen coordinate z in the triangle. + */ +v3_t +tri_bary_coord( + const tri_t * const t, + const v3_t * const p +) +{ + const float t1x = t->p[1].p[0] - t->p[0].p[0]; + const float t1y = t->p[1].p[1] - t->p[0].p[1]; + const float t1z = t->p[1].p[2] - t->p[0].p[2]; + const float t2x = t->p[2].p[0] - t->p[0].p[0]; + const float t2y = t->p[2].p[1] - t->p[0].p[1]; + const float t2z = t->p[2].p[2] - t->p[0].p[2]; + const float px = p->p[0] - t->p[0].p[0]; + const float py = p->p[1] - t->p[0].p[1]; + + float a = (px * t2y - py * t2x) / (t1x * t2y - t2x * t1y); + float b = (px * t1y - py * t1x) / (t2x * t1y - t1x * t2y); + v3_t v = {{ + a, + b, + t->p[0].p[2] + a * t1z + b * t2z, + }}; + + //v.p[2] = 1.0 - v.p[0] - v.p[1]; + + return v; +} + +// Returns true if the bary centry point is inside the triangle +// which means that a and b are non-negative and sum to less than 1 +int +tri_bary_inside( + const v3_t p +) +{ + const float a = p.p[0]; + const float b = p.p[1]; + + return 0 <= a && 0 <= b && a + b <= 1; +} + + /** Compute the points of intersection for two segments in 2d, and z points. * * This is a specialized ray intersection algorithm for the @@ -235,18 +290,22 @@ hidden_intersect( const float p0_x = p0->p[0]; const float p0_y = p0->p[1]; const float p0_z = p0->p[2]; + const float p1_x = p1->p[0]; const float p1_y = p1->p[1]; const float p1_z = p1->p[2]; + const float p2_x = p2->p[0]; const float p2_y = p2->p[1]; const float p2_z = p2->p[2]; + const float p3_x = p3->p[0]; const float p3_y = p3->p[1]; const float p3_z = p3->p[2]; const float s1_x = p1_x - p0_x; const float s1_y = p1_y - p0_y; + const float s2_x = p3_x - p2_x; const float s2_y = p3_y - p2_y; @@ -279,15 +338,12 @@ if(0) fprintf(stderr, "collision: %.0f,%.0f,%.0f->%.0f,%.0f,%.0f %.0f,%.0f,%.0f- r1 ); - const float ix = p0_x + (r0 * s1_x); - const float iy = p0_y + (r0 * s1_y); - // compute the z intercept for each on the two different coordinates if(l0_int) { *l0_int = (v3_t){{ - ix, - iy, + p0_x + r0 * s1_x, + p0_y + r0 * s1_y, p0_z + r0 * (p1_z - p0_z) }}; } @@ -295,8 +351,8 @@ if(0) fprintf(stderr, "collision: %.0f,%.0f,%.0f->%.0f,%.0f,%.0f %.0f,%.0f,%.0f- if(l1_int) { *l1_int = (v3_t){{ - ix, - iy, + p2_x + r1 * s2_x, + p2_y + r1 * s2_y, p2_z + r1 * (p3_z - p2_z) }}; } @@ -305,270 +361,7 @@ if(0) fprintf(stderr, "collision: %.0f,%.0f,%.0f->%.0f,%.0f,%.0f %.0f,%.0f,%.0f- } -/* - * Recursive algorithm: - * Given a line segment and a list of triangles, - * find if the line segment crosses any triangle. - * If it crosses a triangle the segment will be shortened - * and an additional one might be created. - * Recusively try intersecting the new segment (starting at the same triangle) - * and then continue trying the shortened segment. - */ - -void -tri_seg_intersect( - const tri_t * zlist, - seg_t * s, - seg_t ** slist_visible -) -{ - const float p0z = s->p[0].p[2]; - const float p1z = s->p[1].p[2]; - const float seg_max_z = max(p0z, p1z); - - // avoid processing empty segments - const float seg_len = v3_len(&s->p[0], &s->p[1]); - if (seg_len < EPS) - return; - -static int recursive; -recursive++; - -//fprintf(stderr, "%d: processing segment ", recursive); seg_print(s); - fprintf(stderr, "--- recursive %d\n", recursive); - seg_print(s); - - for( const tri_t * t = zlist ; t ; t = t->next ) - { - // if the segment is closer than the triangle, - // then we no longer have to check any further into - // the zlist (it is sorted by depth). - if (seg_max_z <= t->min[2]) - break; - #if 0 - // make sure that we're not comparing to our own triangle - // or one that shares an edge with us (which might be in - // a different order) - if (v2_eq(s->src[0].p, t->p[0].p, 0.0005) - && v2_eq(s->src[1].p, t->p[1].p, 0.0005)) - continue; - if (v2_eq(s->src[0].p, t->p[1].p, 0.0005) - && v2_eq(s->src[1].p, t->p[2].p, 0.0005)) - continue; - if (v2_eq(s->src[0].p, t->p[2].p, 0.0005) - && v2_eq(s->src[1].p, t->p[0].p, 0.0005)) - continue; - if (v2_eq(s->src[0].p, t->p[1].p, 0.0005) - && v2_eq(s->src[1].p, t->p[0].p, 0.0005)) - continue; - if (v2_eq(s->src[0].p, t->p[2].p, 0.0005) - && v2_eq(s->src[1].p, t->p[1].p, 0.0005)) - continue; - if (v2_eq(s->src[0].p, t->p[0].p, 0.0005) - && v2_eq(s->src[1].p, t->p[2].p, 0.0005)) - continue; -#endif - - if (tri_debug >= 2) - tri_print(t); -/* - // if the segment is co-linear to any of the - // triangle edges, include it - for(int i = 0 ; i < 3 ; i++) - { - if (parallel( - &s->p[0], &s->p[1], - &t->p[i], &t->p[(i+1)%3] - )) - goto next_segment; - } -*/ - - float z0, z1; - int inside0 = tri_find_z(t, &s->p[0], &z0); - int inside1 = tri_find_z(t, &s->p[1], &z1); - - if (tri_debug >= 2 && (inside0 || inside1)) - { - fprintf(stderr, "inside %d %d\n", inside0, inside1); - } - - // if both are inside but the segment is infront of the - // triangle, then we retain the segment. - // otherwies we discard the segment - if (inside0 && inside1) - { - if (s->p[0].p[2] <= z0 - && s->p[1].p[2] <= z1) - continue; - if (tri_debug >= 2) - fprintf(stderr, "BOTH INSIDE\n"); - recursive--; - return; - } - - // split the segment for each intersection with the - // triangle segments and add it to the work queue. - int intersections = 0; - v3_t is[3] = {}; // 3d point of segment intercept - v3_t it[3] = {}; // 3d point of triangle intercept - - for(int j = 0 ; j < 3 ; j++) - { - float ratio = hidden_intersect( - &s->p[0], &s->p[1], - &t->p[j], &t->p[(j+1)%3], - &is[intersections], - &it[intersections] - ); - - if (ratio < 0) - continue; - - if (tri_debug >= 2) - fprintf(stderr, "%d ratio=%.2f\n", j, ratio); - intersections++; - } - - - // if none of them intersect, we keep looking - if (intersections == 0) - continue; - - if (tri_debug >= 2) - fprintf(stderr, "%d intersections\n", intersections); - - if (intersections == 3) - { - // this likely means that the triangle is very, very - // small, so let's just ignore this triangle - if (tri_debug >= 2) - fprintf(stderr, "Three intersections\n"); - continue; - } - - - if (intersections == 2) - { - // figure out how far it is to each of the intersections - const float d00 = v3_len(&s->p[0], &is[0]); - const float d01 = v3_len(&s->p[0], &is[1]); - const float d10 = v3_len(&s->p[1], &is[0]); - const float d11 = v3_len(&s->p[1], &is[1]); - - if (tri_debug >= 2) - fprintf(stderr, "Two intersections\n"); - - // discard segments that have two interesections that match - // the segment exactly (distance from segment ends to - // intersection point close enough to zero). - if (d00 < EPS && d11 < EPS) - { - recursive--; - return; - } - if (d01 < EPS && d10 < EPS) - { - recursive--; - return; - } - - // if the segment intersection is closer than the triangle, - // then we do nothing. degenerate cases are not handled - if (d00 <= d01 - && is[0].p[2] <= it[0].p[2] - && is[1].p[2] <= it[1].p[2]) - continue; - if (d00 > d01 - && is[1].p[2] <= it[0].p[2] - && is[0].p[2] <= it[1].p[2]) - continue; - - // segment is behind the triangle, - // we have to create a new segment - // and shorten the existing segment - // find the two intersections that we have - // update the src field - - // we need to create a new segment - seg_t * news; - if (d00 < d01) - { - // split from p0 to ix0 - news = seg_new(s->p[0], is[0]); - news->src[0] = s->src[0]; - news->src[1] = s->src[1]; - s->p[0] = is[1]; - } else { - // split from p0 to ix1 - news = seg_new(s->p[0], is[1]); - news->src[0] = s->src[0]; - news->src[1] = s->src[1]; - s->p[0] = is[0]; - } - - // recursively start splitting the new segment - // starting at the next triangle down the z-depth - tri_seg_intersect(zlist->next, news, slist_visible); - - // continue splitting our current segment - continue; - } - - if (intersections == 1) - { - // if there is an intersection, but the segment intercept - // is closer than the triangle intercept, then no problem. - // we do not bother with degenerate cases of intersecting - // triangles - if (is[0].p[2] <= it[0].p[2] - && is[1].p[2] <= it[0].p[2]) - { - //svg_line("#00FF00", s->p[0].p, s->p[1].p, 10); - continue; - } - - if (inside0) - { - // shorten it on the 0 side - s->p[0] = is[0]; -// huh? shouldn't we process this one? -return; - continue; - } else - if (inside1) - { - // shorten it on the 1 side - s->p[1] = is[0]; -// huh? shouldn't we process this one? -return; - continue; - } else { - // both outside, but an intersection? - // split at that point and hope for the best - seg_t * const news = seg_new(s->p[0], is[0]); - news->src[0] = s->src[0]; - news->src[1] = s->src[1]; - s->p[0] = is[0]; - - tri_seg_intersect(zlist->next, news, slist_visible); - // continue splitting our current segment - continue; - } - } - -next_segment: - continue; - } - - // if we've reached here the segment is visible - // and should be added to the visible list - s->next = *slist_visible; - *slist_visible = s; - recursive--; -} - /* * Fast check to see if t2 is entire occluded by t. @@ -612,3 +405,323 @@ tri_behind( // they are all on the same side return 0; } +#endif + + + +/* + tri_no_intersection, // nothing changed + tri_infront, // segment is in front of the triangle + tri_hidden, // segment is completely occluded + tri_clipped, // segment is partially occluded on one end + tri_split, // segment is partially occluded in the middle +*/ +tri_intersect_t +tri_seg_intersect( + const tri_t * t, + seg_t * s, + seg_t ** new_seg // only if tri_split +) +{ + // avoid processing nearly empty segments + const float seg_len = v3_len(&s->p[0], &s->p[1]); + if (seg_len < EPS) + return tri_hidden; + + const v3_t p_max = v3_max(s->p[0], s->p[1]); + const v3_t p_min = v3_min(s->p[0], s->p[1]); + + // if the segment is closer than the triangle, + // then we no longer have to check any further into + // the zlist (it is sorted by depth). + if (p_max.p[2] <= t->min.p[2]) + return tri_infront; + + // check for four quadrant outside the bounding box + // of the triangle min/max, which would have no chance + // of intersecting with the triangle + if (p_min.p[0] < t->min.p[0] + && p_max.p[0] < t->min.p[0]) + return tri_no_intersection; + if (p_min.p[1] < t->min.p[1] + && p_max.p[1] < t->min.p[1]) + return tri_no_intersection; + if (p_min.p[0] > t->max.p[0] + && p_max.p[0] > t->max.p[0]) + return tri_no_intersection; + if (p_min.p[1] > t->max.p[1] + && p_max.p[1] > t->max.p[1]) + return tri_no_intersection; + + // there is a possibility that this line crosses the triangle + // compute the coordinates in triangle space + const v3_t tp0 = tri_bary_coord(t, &s->p[0]); + const v3_t tp1 = tri_bary_coord(t, &s->p[1]); + + // if both are inside and not both on the same edge of + // the triangle, then the segment is totally hidden. + if (tri_bary_inside(tp0) && tri_bary_inside(tp1)) + { + // if the segment z is closer than the triangle z + // then the segment is in front of the triangle + if (s->p[0].p[2] < tp0.p[2] && s->p[1].p[2] < tp1.p[2]) + return tri_no_intersection; + + // if the barycentric coord is 0 for the same edge + // for both points, then it is on the original line + if (tp0.p[0] < EPS && tp1.p[0] < EPS) + return tri_no_intersection; + if (tp0.p[1] < EPS && tp1.p[1] < EPS) + return tri_no_intersection; + + // compute the third barycentric coordinate and check + float c0 = 1.0 - tp0.p[0] - tp0.p[1]; + float c1 = 1.0 - tp1.p[0] - tp1.p[1]; + if (c0 < EPS && c1 < EPS) + return tri_no_intersection; + + // it is not on an edge and not infront of the triangle + // so the segment is totally occluded + return tri_hidden; + } + + // find the intersection point for each of the three + // sides of the triangle + v3_t is[3] = {}; // 3d point of segment intercept + v3_t it[3] = {}; // 3d point of triangle intercept + float ratio[3] = {}; // length along the line + int intersections = 0; + + for(int j = 0 ; j < 3 ; j++) + { + ratio[intersections] = hidden_intersect( + &s->p[0], &s->p[1], + &t->p[j], &t->p[(j+1)%3], + &is[intersections], + &it[intersections] + ); + + if (ratio[intersections] < 0) + continue; + + // if the segment intersection is closer than the + // triangle intersection, this does not count as + // an intersection and we can ignore it. + if (is[intersections].p[2] < it[intersections].p[2]) + continue; + + if (tri_debug >= 2) + { + fprintf(stderr, "%d ratio=%.2f %+6.1f", j, ratio[intersections], it[intersections].p[2]); + v3_print(is[intersections]); + } + intersections++; + } + + // check for duplicate intersections, which happens if + // the lines go through at precisely the corners + // this might mean that we hit exactly at one + // point and two of the points are the same + if (intersections == 3) + { + if (v3_eq(is[0], is[2])) + intersections--; + else + if (v3_eq(is[1], is[2])) + intersections--; + else + if (v3_eq(is[0], is[1])) + { + intersections--; + is[1] = is[2]; + it[1] = it[2]; + } + } + + if (intersections == 2 && v3_eq(is[0], is[1])) + intersections--; + + // no intersections? there is nothing to do + if (intersections == 0) + return tri_no_intersection; + + // three intersections? maybe a very small triangle + if (intersections == 3) + { + fprintf(stderr, "THREE INTERSECTIONS?\n"); + fprintf(stderr, "is0="); v3_print(is[0]); + fprintf(stderr, "is1="); v3_print(is[1]); + fprintf(stderr, "is2="); v3_print(is[2]); + svg_line("#0000FF", s->p[0].p, s->p[1].p, 8); + return tri_no_intersection; + } + + if (intersections == 1) + { + if (tri_bary_inside(tp0)) + { + // if the intercept point on the segment is + // closer than the intercept point on the triangle edge, + // then there is no occlusion + if (is[0].p[2] <= it[0].p[2]) + return tri_no_intersection; + + // clipped from intersection to p1 + s->p[0] = is[0]; + return tri_clipped; + } + + if (tri_bary_inside(tp1)) + { + // if the intercept point on the segment is + // closer than the intercept point on the triangle edge, + // then there is no occlusion + if (is[0].p[2] < it[0].p[2]) + return tri_no_intersection; + + // clipped from p0 to intersection + s->p[1] = is[0]; + return tri_clipped; + } + + // something isn't right. maybe we have a small triangle? + fprintf(stderr, "ONE INTERSECTION?"); +/* + svg_line("#FFFF00", s->p[0].p, s->p[1].p, 20); + svg_line("#000080", t->p[0].p, t->p[1].p, 8); + svg_line("#000080", t->p[1].p, t->p[2].p, 8); + svg_line("#000080", t->p[2].p, t->p[0].p, 8); +*/ + seg_print(s); + v3_print(tp0); + v3_print(tp1); + tri_print(t); + *new_seg = seg_new(is[0], s->p[1]); + s->p[1] = is[0]; + return tri_split; + } + + // two intersections: find the one that is closer to p0 + // modify the existing segment and create a new segment + const float d00 = v3_mag2(v3_sub(is[0], s->p[0])); + const float d01 = v3_mag2(v3_sub(is[1], s->p[0])); + const float d10 = v3_mag2(v3_sub(is[0], s->p[1])); + const float d11 = v3_mag2(v3_sub(is[1], s->p[1])); + + // if any of the intersections points are zero from an + // end point on the segment, then skip that part + if (tri_debug > 4) + { + seg_print(s); + tri_print(t); + fprintf(stderr, "d: %f %f %f %f\n", d00, d01, d10, d11); + } + + if (d00 < EPS && d11 < EPS) + return tri_hidden; + if (d01 < EPS && d10 < EPS) + return tri_hidden; + + if (d00 < EPS) + { + s->p[0] = is[1]; + return tri_clipped; + } else + if (d01 < EPS) + { + s->p[0] = is[0]; + return tri_clipped; + } else + if (d10 < EPS) + { + s->p[1] = is[1]; + return tri_clipped; + } else + if (d11 < EPS) + { + s->p[1] = is[0]; + return tri_clipped; + } + + // neither end points match, so create a new segment + // that excludes the space covered by the triangle. + // determine which is closer to point is[0] + if (d00 < d01) + { + // p0 is closer to is0, so new segment is is1 to p1 + *new_seg = seg_new(is[1], s->p[1]); + s->p[1] = is[0]; + } else { + // p0 is closer to is1, so new segment is is0 to p1 + *new_seg = seg_new(is[0], s->p[1]); + s->p[1] = is[1]; + } + +fprintf(stderr, "SPLIT: "); +seg_print(*new_seg); + return tri_split; +} + + +int +tri_seg_hidden( + const tri_t * zlist, + seg_t * s, + seg_t ** slist_visible +) +{ + int count = 0; + fprintf(stderr, "TEST: "); + seg_print(s); + + for( const tri_t * t = zlist ; t ; t = t->next ) + { + seg_t * new_seg = NULL; + //tri_print(t); + tri_intersect_t type = tri_seg_intersect(t, s, &new_seg); + //fprintf(stderr, "rc=%d\n", type); + + // if there is no intersection or if the segment has + // been clipped on one side, keep looking + if (type == tri_no_intersection) + continue; + if (type == tri_clipped) + { + //fprintf(stderr, "CLIP: "); + seg_print(s); + continue; + } + + // if this segment is infront of this triangle then we can + // stop searching + if (type == tri_infront) + break; + + // if this segment is totally occluded, we're done + if (type == tri_hidden) + return count; + + // if this line has been split into two, process the + // new segment starting at the next triangle since it + // has already intersected this one + if (type == tri_split) + { + static int recursive; + if (tri_debug > 4) fprintf(stderr, "RECURSIVE %d\n", recursive++); + int new_count = tri_seg_hidden(t->next, new_seg, slist_visible); + if (tri_debug > 4) fprintf(stderr, "END %d: %d segments\n", --recursive, new_count); + if (tri_debug > 4) fprintf(stderr, "CLIP: "); + if (tri_debug > 4) seg_print(s); + count += new_count; + continue; + } + + fprintf(stderr, "unknown type %d\n", type); + return -1; + } + + // we've reached the end and it is still visible + s->next = *slist_visible; + *slist_visible = s; + return ++count; +} diff --git a/tri.h b/tri.h index e51d5b1..96e0392 100644 --- a/tri.h +++ b/tri.h @@ -6,6 +6,7 @@ #include "v3.h" #include "seg.h" +#include "svg.h" extern int tri_debug; @@ -15,8 +16,8 @@ struct _tri_t v3_t p[3]; // camera space v3_t normal; // camera space v3_t normal_xyz; // original xyz space - float min[3]; // camera space - float max[3]; // camera space + v3_t min; // camera space + v3_t max; // camera space tri_t * next; tri_t ** prev; }; @@ -108,17 +109,18 @@ hidden_intersect( /* - * Recursive algorithm: * Given a line segment and a list of triangles, * find if the line segment crosses any triangle. * If it crosses a triangle the segment will be shortened * and an additional one might be created. * Recusively try intersecting the new segment (starting at the same triangle) * and then continue trying the shortened segment. + * + * Line segments will be added to the visible list. + * Returns the number of new elements created */ - -void -tri_seg_intersect( +int +tri_seg_hidden( const tri_t * zlist, seg_t * s, seg_t ** slist_visible @@ -134,4 +136,30 @@ tri_behind( const tri_t * const t2 ); + +/* + * There are four possible line/triangle intersections. + */ +typedef enum { + tri_no_intersection, // nothing changed + tri_infront, // segment is in front of the triangle + tri_hidden, // segment is completely occluded + tri_clipped, // segment is partially occluded on one end + tri_split, // segment is partially occluded in the middle +} tri_intersect_t; + + +tri_intersect_t +tri_seg_intersect( + const tri_t * tri, + seg_t * s, + seg_t ** new_seg // only if tri_split +); + +v3_t +tri_bary_coord( + const tri_t * const t, + const v3_t * const p +); + #endif diff --git a/unfold.c b/unfold.c index 5b675a1..81cf73b 100644 --- a/unfold.c +++ b/unfold.c @@ -92,7 +92,7 @@ edge_eq2( const v3_t * const v10 = &t1->p[e1]; const v3_t * const v11 = &t1->p[(e1+1) % 3]; - if (v3_eq(v00, v11) && v3_eq(v01, v10)) + if (v3_eq(*v00, *v11) && v3_eq(*v01, *v10)) return 1; return 0; @@ -613,11 +613,11 @@ coplanar_check( for (int i = 0 ; i < 3 ; i++) { x4 = f2->p[i]; - if (v3_eq(&x1, &x4)) + if (v3_eq(x1, x4)) continue; - if (v3_eq(&x2, &x4)) + if (v3_eq(x2, x4)) continue; - if (v3_eq(&x3, &x4)) + if (v3_eq(x3, x4)) continue; break; } diff --git a/v3.h b/v3.h index a5a15ce..da97862 100644 --- a/v3.h +++ b/v3.h @@ -5,9 +5,10 @@ #define _papercraft_v3_h_ #include +#include #include -#define EPS 0.001 +#define EPS 0.00001 #ifndef M_PI #define M_PI 3.1415926535897932384 @@ -52,13 +53,13 @@ typedef struct static inline int v3_eq( - const v3_t * v1, - const v3_t * v2 + const v3_t v1, + const v3_t v2 ) { - float dx = v1->p[0] - v2->p[0]; - float dy = v1->p[1] - v2->p[1]; - float dz = v1->p[2] - v2->p[2]; + float dx = v1.p[0] - v2.p[0]; + float dy = v1.p[1] - v2.p[1]; + float dz = v1.p[2] - v2.p[2]; if (-EPS < dx && dx < EPS && -EPS < dy && dy < EPS @@ -209,6 +210,36 @@ v3_cross( } +static inline v3_t +v3_min( + const v3_t a, + const v3_t b +) +{ + v3_t c = { { + min(a.p[0], b.p[0]), + min(a.p[1], b.p[1]), + min(a.p[2], b.p[2]), + }}; + return c; +} + + +static inline v3_t +v3_max( + const v3_t a, + const v3_t b +) +{ + v3_t c = { { + max(a.p[0], b.p[0]), + max(a.p[1], b.p[1]), + max(a.p[2], b.p[2]), + }}; + return c; +} + + // Compute the length of a line in screen space, ignoring Z static inline float v3_dist_2d( @@ -223,4 +254,57 @@ v3_dist_2d( } +static inline void +v3_print(const v3_t p) +{ + fprintf(stderr, "%+6.1f %+6.1f %+6.1f\n", + p.p[0], + p.p[1], + p.p[2] + ); +} + +typedef struct { + float p[4]; +} v4_t; + +typedef struct { + float m[4][4]; +} m44_t; + + +static inline void +m44_mult( + m44_t * r, + const m44_t * a, + const m44_t * b +) +{ + for(int i = 0 ; i < 4 ; i++) + { + for(int j = 0 ; j < 4 ; j++) + { + float d = 0; + for(int k = 0 ; k < 4 ; k++) + d += a->m[i][k] * b->m[k][j]; + r->m[i][j] = d; + } + } +} + + +static inline v4_t +m44_multv( + const m44_t * const m, + const v4_t * const v +) +{ + v4_t p = {}; + for (int i = 0 ; i < 4 ; i++) + for (int j = 0 ; j < 4 ; j++) + p.p[i] += m->m[i][j] * v->p[j]; + + return p; +} + #endif diff --git a/wireframe.c b/wireframe.c index 13559ea..3f69dec 100644 --- a/wireframe.c +++ b/wireframe.c @@ -62,7 +62,7 @@ coplanar_check( for (int i = 0 ; i < 3 ; i++) { for (int j = 0 ; j < 3 ; j++) - if (v3_eq(&f1->p[i], &f2->p[j])) + if (v3_eq(f1->p[i], f2->p[j])) match[i] = 1; } @@ -132,7 +132,7 @@ coplanar_check( return mask; #else // if the normals are close enough, then it is coplanner - if (v3_eq(&f1->normal, &f2->normal)) + if (v3_eq(f1->normal, f2->normal)) return mask; else return 0; @@ -184,7 +184,7 @@ stl_vertex_find( { stl_vertex_t * const v = vertices[x]; - if (v3_eq(&v->p, p)) + if (v3_eq(v->p, *p)) return v; }