commit 6857de388e82da7b0563665e2749d2f8353a3f0d Author: andreika-git Date: Wed Oct 4 23:03:24 2023 +0300 Initial release circa 2015 diff --git a/src/mh.cpp b/src/mh.cpp new file mode 100644 index 0000000..b968808 --- /dev/null +++ b/src/mh.cpp @@ -0,0 +1,200 @@ +/* +Marching Hedrons LUT generator +author: andrei (c) 2015 + +*/ + +#include +#include +#include + +#include "mh.h" + +#define OUTPUT_BIN + +//#define OUTPUT_TEXT + +const int deviceIdx = 0; +//const Distribution dist = UniformDistribution; + +////////////////////////////////////////////////////////////////////////// + + +void addFace(FaceHVec *faces, int v0, int v1, int v2, int skipVert) +{ + if (v0 == skipVert || v1 == skipVert || v2 == skipVert) + return; + faces->push_back(Face::getFace(v0, v1, v2)); +} + + +//void setPoints(Point3HVec *pointNum, int *idx, int mask) +void setPoints(Vec3HVec *pointNum, int *idx, int mask, int inv, int index) +{ + for (int i = 0, j = 0; i < 26; i++) + { + int e0 = mask & (1 << MH::cell_edges[i][0]); + int e1 = mask & (1 << MH::cell_edges[i][1]); + if (MH::cell_edges[i][0] == 14) // special case for central edges + { + if (!e0) + continue; + } + if ((e0 && !e1) || (!e0 && e1)) + { + Vec3 p; + if (index >= 0 && i < 8) + { + float coeff; + coeff = MH::edge_indices[index][i & 3] == 0 ? 0.2f : 0.8f; + p = MH::cell_vertices[MH::cell_edges[i][0]].interp(MH::cell_vertices[MH::cell_edges[i][1]], coeff); + } + else + { + int k; + if (inv == 0) + k = 0; + else if (i >= 8) + k = 1; + else + k = 0; + if (i < 12) + { + if ((mask & (1 << MH::edgeFacesCentralVertices[i][0])) && (mask & (1 << MH::edgeFacesCentralVertices[i][1]))) + { + if (inv == 0) + k = 1; + else if (i >= 8) + k = 0; + else + k = 1; + } + } + p.x = MH::cell_hedron_vertices[k][i][0]; + p.y = MH::cell_hedron_vertices[k][i][1]; + p.z = MH::cell_hedron_vertices[k][i][2]; + } + + + idx[j++] = i; + pointNum->push_back(p); + } + } +} + +////////////////////////////////////////////////////// + +int main(int argc, const char* argv[]) +{ +#ifdef OUTPUT_BIN + FILE *fpi = fopen("mh_idx.bin", "wb"); + FILE *fpi0 = fopen("mh_offs.bin", "wb"); + FILE *fpih = fopen("mh_h_idx.bin", "wb"); + FILE *fp = fopen("mh.bin", "wb"); + int offs = 0; +#endif + + MH::initConvexHull(); + + std::vector faceGroupStore[7]; + int totalIdx = 0; + + for (int inv = 0; inv < 2; inv++) + { + for (int idx = 0; idx < 16; idx++) + { + if (inv == 0) + std::cout << "Starting Hedrons #" << idx << "...\n"; + else + std::cout << "Starting Holes #" << idx << "...\n"; + + for (int i = 0; i < (1 << 14); i++) + { + Vec3HVec pointVec; + int indices[26]; + + if (i > 0 && (i % 1000) == 0) + { + std::cout << "// " << i << std::endl; + } + + int mask = i | (1 << 14); + + setPoints(&pointVec, indices, mask, inv, idx); + + FaceHVec faces; + FaceHVec faceGroups[7]; + + if (pointVec.size() == 3) + { + addFace(&faces, 0, 1, 2, -1); + faceGroups[0].push_back(faces[0]); + } + else + { + FaceHVec fac; + MH::computeConvexHull(pointVec, fac); + MH::removeSpecialFaces(indices, fac, faces); + MH::FixNormals(pointVec, faces); + MH::SortOutFaces(indices, faces, faceGroups); + } + +#ifdef OUTPUT_TEXT + std::cout << " /*" << i << " (N=" << faces.size() << ") */ { "; + for (int j = 0; j < faces.size(); j++) + { + std::cout << "{ " << indices[faces[j].v[0]] << " /*" << faces[j].v[0] << "*/, " << indices[faces[j].v[1]] << "/*" << faces[j].v[1] << "*/, " << indices[faces[j].v[2]] << " /*" << faces[j].v[2] << "*/ }, "; + } + + std::cout << "{ -1, -1, -1 } }, " << std::endl; +#endif + +#ifdef OUTPUT_BIN + for (int f = 0; f < 7; f++) + { + IntHVec offsList; + bool addedNew = false; + int fidx = MH::AddUniqueFaceGroups(faceGroupStore[f], indices, faceGroups[f], &addedNew); + + FaceRecord &fr = faceGroupStore[f][fidx]; + if (addedNew) // new record added + { + fr.totalIdx = totalIdx; + fwrite(&offs, sizeof(offs), 1, fpi0); + fputc(fr.faces.size(), fp); + offs ++; + for (size_t j = 0; j < fr.faces.size(); j++) + { + fputc(fr.faces[j].v[0], fp); + fputc(fr.faces[j].v[1], fp); + fputc(fr.faces[j].v[2], fp); + offs += 3; + } + totalIdx++; + } + + fputc(fr.totalIdx & 0xff, inv == 0 ? fpi : fpih); + fputc((fr.totalIdx >> 8) & 0xff, inv == 0 ? fpi : fpih); + //fputc((fr.totalIdx >> 16) & 0xff, inv == 0 ? fpi : fpih); + } +#endif + } + } + } + +#ifdef OUTPUT_BIN + fclose(fp); + fclose(fpi); + fclose(fpi0); + fclose(fpih); +#endif + + std::cout << "NUM: " << faceGroupStore[0].size() << " " << faceGroupStore[1].size() << " " << faceGroupStore[2].size() + << " " << faceGroupStore[3].size() << " " << faceGroupStore[4].size() << " " << faceGroupStore[5].size() + << " " << faceGroupStore[6].size() << std::endl; + std::cout << "TOTALIDX = " << totalIdx << std::endl; + + return 0; +} + + diff --git a/src/mh.h b/src/mh.h new file mode 100644 index 0000000..ca1379a --- /dev/null +++ b/src/mh.h @@ -0,0 +1,118 @@ +/* +Marching Hedrons LUT generator +author: andrei (c) 2015 + +*/ + +#ifndef _MH_H_ +#define _MH_H_ + +#include +#include "vec.h" + +struct Face +{ + int v[3]; + + static Face getFace(int v0, int v1, int v2) + { + Face f; + // 012 + f.v[0] = v0; f.v[1] = v1; f.v[2] = v2; + return f; + } + + static Face getSortedFace(int v0, int v1, int v2) + { + Face f; + // 012 + if (v0 < v1 && v1 < v2) + { + f.v[0] = v0; f.v[1] = v1; f.v[2] = v2; + } + // 021 + if (v0 < v2 && v2 < v1) + { + f.v[0] = v0; f.v[1] = v2; f.v[2] = v1; + } + // 102 + if (v1 < v0 && v0 < v2) + { + f.v[0] = v1; f.v[1] = v0; f.v[2] = v2; + } + // 120 + if (v1 < v2 && v2 < v0) + { + f.v[0] = v1; f.v[1] = v2; f.v[2] = v0; + } + // 201 + if (v2 < v0 && v0 < v1) + { + f.v[0] = v2; f.v[1] = v0; f.v[2] = v1; + } + // 210 + if (v2 < v1 && v1 < v0) + { + f.v[0] = v2; f.v[1] = v1; f.v[2] = v0; + } + return f; + } + + friend bool operator== (const Face &n1, const Face &n2) + { + return Face::equal(n1, n2); + } + + friend bool operator!= (const Face &n1, const Face &n2) + { + return !Face::equal(n1, n2); + }; + +private: + static bool equal(const Face & f0, const Face & f1) + { + return (f0.v[0] == f1.v[0] && f0.v[1] == f1.v[1] && f0.v[2] == f1.v[2]) || + (f0.v[0] == f1.v[0] && f0.v[1] == f1.v[2] && f0.v[2] == f1.v[1]) || + (f0.v[0] == f1.v[1] && f0.v[1] == f1.v[0] && f0.v[2] == f1.v[2]) || + (f0.v[0] == f1.v[1] && f0.v[1] == f1.v[2] && f0.v[2] == f1.v[0]) || + (f0.v[0] == f1.v[2] && f0.v[1] == f1.v[1] && f0.v[2] == f1.v[0]) || + (f0.v[0] == f1.v[2] && f0.v[1] == f1.v[0] && f0.v[2] == f1.v[1]); + } +}; + +typedef std::vector FaceHVec; + +typedef std::vector Vec3HVec; + +typedef std::vector IntHVec; + +struct FaceRecord +{ + FaceHVec faces; + int idx; + int totalIdx; +}; + +class MH +{ +public: + static const int cell_edges[26][2]; + static const Vec3 cell_vertices[8]; + static const float cell_hedron_vertices[2][26][3]; + static const int faceVertices[6][9]; + static const int edgeFacesCentralVertices[12][2]; + static int edge_indices[16][4]; + +public: + + static bool initConvexHull(); + static bool computeConvexHull(const Vec3HVec & pointVec, FaceHVec & faces); + static bool removeSpecialFaces(int *idx, FaceHVec & faces, FaceHVec & outfaces); + static void FixNormals(const Vec3HVec & pointVec, FaceHVec & outfaces); + static void SortOutFaces(int *idx, FaceHVec & faces, FaceHVec faceGroups[7]); + static int AddUniqueFaceGroups(std::vector &faceStore, int *idx, FaceHVec &faceGroups, bool *addedNew); + + static int tr_tri_intersect3D(const Vec3 & C1, const Vec3 & P1, const Vec3 & P2, const Vec3 & D1, const Vec3 & Q1, const Vec3 & Q2); +}; + +#endif // of _MH_H_ diff --git a/src/mh_convex.cpp b/src/mh_convex.cpp new file mode 100644 index 0000000..6207329 --- /dev/null +++ b/src/mh_convex.cpp @@ -0,0 +1,255 @@ +/* +Marching Hedrons LUT generator +author: andrei (c) 2015 + +Marching hedron convex test + +*/ + +#include + +#include "mh.h" + +#define EPS 1e-6 + + +bool MH::initConvexHull() +{ + for (int i = 0; i < 16; i++) + { + for (int j = 0; j < 4; j++) + edge_indices[i][j] = (i >> j) & 1; + } + return true; +} + +bool MH::computeConvexHull(const Vec3HVec & pointVec, FaceHVec & faces) +{ + int s = pointVec.size(); + FaceHVec allfaces; + for (int i0 = 0; i0 < s; i0++) + { + for (int i1 = i0 + 1; i1 < s; i1++) + { + for (int i2 = i1 + 1; i2 < s; i2++) + { + allfaces.push_back(Face::getFace(i0, i1, i2)); + } + } + } + + for (size_t i = 0; i < allfaces.size(); i++) + { + Face &f = allfaces[i]; + Vec3 n = (pointVec[f.v[2]] - pointVec[f.v[0]]).cross(pointVec[f.v[1]] - pointVec[f.v[0]]); + bool pos = false; + bool neg = false; + for (int j = 0; j < s; j++) + { + if (j == f.v[0] || j == f.v[1] || j == f.v[2]) + continue; + Vec3 p = pointVec[j] - pointVec[f.v[0]]; + float d = p.dot(n); + if (d < -EPS) + neg = true; + else if (d > EPS) + pos = true; + } + // if there are vertices on both sides of a face plane + if (pos && neg) + continue; + + // test if this face intersects with already added faces + bool intersect = false; + for (size_t k = 0; k < faces.size(); k++) + { + Face &f2 = faces[k]; + + int ret = tr_tri_intersect3D(pointVec[f.v[0]], pointVec[f.v[1]] - pointVec[f.v[0]], pointVec[f.v[2]] - pointVec[f.v[0]], + pointVec[f2.v[0]], pointVec[f2.v[1]] - pointVec[f2.v[0]], pointVec[f2.v[2]] - pointVec[f2.v[0]]); + if (ret == 1) + { + intersect = true; + break; + } + } + if (!intersect) + faces.push_back(f); + } + + return true; +} + +bool MH::removeSpecialFaces(int *idx, FaceHVec & faces, FaceHVec & outfaces) +{ + // if there are special faces (faces that have all 3 verts of 12 'mc' edges), then remove them all - IF there are any other (non-special faces) + size_t i; + for (i = 0; i < faces.size(); i++) + { + if (idx[faces[i].v[0]] >= 12 || idx[faces[i].v[1]] >= 12 || idx[faces[i].v[2]] >= 12) + { + // we found non-special face + for (i = 0; i < faces.size(); i++) + { + if (idx[faces[i].v[0]] >= 12 || idx[faces[i].v[1]] >= 12 || idx[faces[i].v[2]] >= 12) + { + outfaces.push_back(faces[i]); + } + } + return true; + } + } + + return false; +} + +void MH::FixNormals(const Vec3HVec & pointVec, FaceHVec & faces) +{ + static const Vec3 c(0.5f, 0.5f, 0.5f), one_third(0.333333333333f, 0.333333333333f, 0.333333333333f); + + for (size_t i = 0; i < faces.size(); i++) + { + Face &f = faces[i]; + Vec3 cf = (pointVec[f.v[0]] + pointVec[f.v[1]] + pointVec[f.v[2]]) * one_third; + Vec3 n = (pointVec[f.v[2]] - pointVec[f.v[0]]).cross(pointVec[f.v[1]] - pointVec[f.v[0]]); + if (n.dot(cf - c) > 0) + { + // swap + int v = f.v[1]; + f.v[1] = f.v[0]; + f.v[0] = v; + } + } +} + +void MH::SortOutFaces(int *idx, FaceHVec & faces, FaceHVec faceGroups[7]) +{ + for (size_t i = 0; i < faces.size(); i++) + { + Face &f = faces[i]; + bool inGroup = false; + for (int j = 0; j < 6; j++) + { + bool fv[3] = { false, false, false }; + for (int k = 0; k < 9; k++) + { + for (int vi = 0; vi < 3; vi++) + { + if (idx[f.v[vi]] == faceVertices[j][k]) + fv[vi] = true; + } + } + // if all 3 face verts are from the same cube face, then put them into face group + if (fv[0] && fv[1] && fv[2]) + { + faceGroups[j].push_back(f); + inGroup = true; + break; + } + } + if (!inGroup) + faceGroups[6].push_back(f); + } +} + +int MH::AddUniqueFaceGroups(std::vector &faceStore, int *idx, FaceHVec &faces, bool *addedNew) +{ + for (size_t j = 0; j < faceStore.size(); j++) + { + FaceRecord & fsvec = faceStore[j]; + if (fsvec.faces.size() != faces.size()) + continue; + bool notFound = false; + for (size_t i = 0; i < faces.size(); i++) + { + int a0 = idx[faces[i].v[0]]; + int a1 = idx[faces[i].v[1]]; + int a2 = idx[faces[i].v[2]]; + bool found = false; + for (size_t k = 0; k < fsvec.faces.size(); k++) + { + int b0 = fsvec.faces[k].v[0]; + int b1 = fsvec.faces[k].v[1]; + int b2 = fsvec.faces[k].v[2]; + if (a0 == b0) + { + if (a1 == b1) + { + if (a2 == b2) + { + found = true; + break; + } + } + else if (a1 == b2) + { + if (a2 == b1) + { + found = true; + break; + } + } + } + else if (a0 == b1) + { + if (a1 == b0) + { + if (a2 == b2) + { + found = true; + break; + } + } + else if (a1 == b2) + { + if (a2 == b0) + { + found = true; + break; + } + } + } + else if (a0 == b2) + { + if (a1 == b0) + { + if (a2 == b1) + { + found = true; + break; + } + } + else if (a1 == b1) + { + if (a2 == b0) + { + found = true; + break; + } + } + } + } + if (!found) + { + notFound = true; + break; + } + } + if (notFound) + continue; + + // ok, we found it + return j; + } + + // not found, add it to the list + FaceRecord newrec; + for (size_t k = 0; k < faces.size(); k++) + { + newrec.faces.push_back(Face::getFace(idx[faces[k].v[0]], idx[faces[k].v[1]], idx[faces[k].v[2]])); + } + newrec.idx = faceStore.size(); + faceStore.push_back(newrec); + *addedNew = true; + return newrec.idx; +} diff --git a/src/mh_tables.cpp b/src/mh_tables.cpp new file mode 100644 index 0000000..9ea20e0 --- /dev/null +++ b/src/mh_tables.cpp @@ -0,0 +1,92 @@ +/* +Marching Hedrons LUT generator +author: andrei (c) 2015 + +*/ + +#include "mh.h" + +const int MH::cell_edges[26][2] = +{ + // the first 12 edges are the same as in MC + { 0, 1 }, { 1, 5 }, { 5, 4 }, { 4, 0 }, + { 3, 2 }, { 2, 6 }, { 6, 7 }, { 7, 3 }, + { 0, 3 }, { 1, 2 }, { 5, 6 }, { 4, 7 }, + + // there are 14 more edges coming from the central point: + // - 8 to corner + { 14, 0 }, { 14, 1 }, { 14, 2 }, { 14, 3 }, + { 14, 4 }, { 14, 5 }, { 14, 6 }, { 14, 7 }, + // - and 6 to facenormal (center of each cube face) + { 14, 8 }, { 14, 9 }, { 14, 10 }, { 14, 11 }, { 14, 12 }, { 14, 13 }, +}; + + +#define E01 0.05 +#define E09 0.95 +#define V00 0.01 +#define V10 0.99 + +const Vec3 MH::cell_vertices[8] = +{ + Vec3(V00, V00, V00), Vec3(V10, V00, V00), Vec3(V10, V10, V00), Vec3(V00, V10, V00), + Vec3(V00, V00, V10), Vec3(V10, V00, V10), Vec3(V10, V10, V10), Vec3(V00, V10, V10), +}; + +#define A05 0.2 + +const float MH::cell_hedron_vertices[2][26][3] = +{ + // low + { + // 12 cube edges + // ++-- + // ++-- + { A05, V00, V00 }, { V10, V00, A05 }, { A05, V00, V10 }, { V00, V00, A05 }, + { A05, V10, V00 }, { V10, V10, A05 }, { A05, V10, V10 }, { V00, V10, A05 }, + { V00, A05, V00 }, { V10, A05, V00 }, { V10, A05, V10 }, { V00, A05, V10 }, + + // 8 center-to-corner edges + { E01, E01, E01 }, { E09, E01, E01 }, { E09, E09, E01 }, { E01, E09, E01 }, + { E01, E01, E09 }, { E09, E01, E09 }, { E09, E09, E09 }, { E01, E09, E09 }, + + // 6 center-to-facenormal edges + { 0.5, 1.0, 0.5 }, { 0.5, 0.5, 0.0 }, { 1.0, 0.5, 0.5 }, { 0.5, 0.5, 1.0 }, { 0.0, 0.5, 0.5 }, { 0.5, 0.0, 0.5 }, + }, + +#undef A05 +#define A05 0.8 + // high + { + // 12 cube edges + { A05, V00, V00 }, { V10, V00, A05 }, { A05, V00, V10 }, { V00, V00, A05 }, + { A05, V10, V00 }, { V10, V10, A05 }, { A05, V10, V10 }, { V00, V10, A05 }, + { V00, A05, V00 }, { V10, A05, V00 }, { V10, A05, V10 }, { V00, A05, V10 }, + + // 8 center-to-corner edges + { E01, E01, E01 }, { E09, E01, E01 }, { E09, E09, E01 }, { E01, E09, E01 }, + { E01, E01, E09 }, { E09, E01, E09 }, { E09, E09, E09 }, { E01, E09, E09 }, + + // 6 center-to-facenormal edges + { 0.5, 1.0, 0.5 }, { 0.5, 0.5, 0.0 }, { 1.0, 0.5, 0.5 }, { 0.5, 0.5, 1.0 }, { 0.0, 0.5, 0.5 }, { 0.5, 0.0, 0.5 }, + } +}; + +const int MH::faceVertices[6][9] = +{ + { 4, 5, 6, 7, 14, 15, 18, 19, 20 }, + { 0, 8, 4, 9, 12, 13, 14, 15, 21 }, + { 1, 9, 5, 10, 13, 14, 18, 17, 22 }, + { 2, 10, 6, 11, 16, 17, 18, 19, 23 }, + { 3, 8, 7, 11, 12, 15, 19, 16, 24 }, + { 0, 1, 2, 3, 12, 13, 17, 16, 25 }, +}; + +const int MH::edgeFacesCentralVertices[12][2] = +{ + { 13, 9 }, { 13, 10 }, { 13, 11 }, { 13, 12 }, + { 8, 9 }, { 8, 10 }, { 8, 11 }, { 8, 12 }, + { 9, 12 }, { 9, 10 }, { 10, 11 }, { 11, 12 }, +}; + +int MH::edge_indices[16][4]; diff --git a/src/mh_tritri.cpp b/src/mh_tritri.cpp new file mode 100644 index 0000000..a03af75 --- /dev/null +++ b/src/mh_tritri.cpp @@ -0,0 +1,356 @@ +/* + +Marching Hedrons LUT generator +author: andrei (c) 2015 + +Come code taken from: + +"A fast triangle to triangle intersection test for collision detection" +Oren Tropp, Ayellet Tal, Ilan Shimshoni +Computer Animation and Virtual Worlds 17(5) 2006, pp 527-535. + +You are free to use the code but cite the paper. + + +The following code tests for 3D triangle triangle intersection. +Main procedures: + +int tr_tri_intersect3D (float *C1, float *P1, float *P2, +float *D1, float *Q1, float *Q2); + +int coplanar_tri_tri(float N[3],float V0[3],float V1[3],float V2[3], +float U0[3],float U1[3],float U2[3]); + + +tr_tri_intersect3D - C1 is a vertex of triangle A. P1,P2 are the two edges originating from this vertex. +D1 is a vertex of triangle B. P1,P2 are the two edges originating from this vertex. +Returns zero for disjoint triangles and non-zero for intersection. + +coplanar_tri_tri - This procedure for testing coplanar triangles for intersection is +taken from Tomas Moller's algorithm. +See article "A Fast Triangle-Triangle Intersection Test", +Journal of Graphics Tools, 2(2), 1997 +V1,V2,V3 are vertices of one triangle with normal N. U1,U2,U3 are vertices of the other +triangle. + +*/ +#include +#include +#include +#include +#include +#include +#ifndef UNIX +#define drand48() (rand()*1.0/RAND_MAX) +#endif + + +#include "vec.h" +#include "mh.h" + +int coplanar_tri_tri(float N[3], float V0[3], float V1[3], float V2[3], + float U0[3], float U1[3], float U2[3]); + +// some vector macros +#define CROSS(dest,v1,v2) \ + dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \ + dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \ + dest[2]=v1[0]*v2[1]-v1[1]*v2[0]; + + + +#define sVpsV_2( Vr, s1, V1,s2, V2);\ +{\ + Vr[0] = s1*V1[0] + s2*V2[0];\ + Vr[1] = s1*V1[1] + s2*V2[1];\ +}\ + +#define myVpV(g,v2,v1);\ +{\ + g[0] = v2[0]+v1[0];\ + g[1] = v2[1]+v1[1];\ + g[2] = v2[2]+v1[2];\ +}\ + +#define myVmV(g,v2,v1);\ +{\ + g[0] = v2[0]-v1[0];\ + g[1] = v2[1]-v1[1];\ + g[2] = v2[2]-v1[2];\ +}\ + +// 2D intersection of segment and triangle. +#define seg_collide3( q, r)\ +{\ + p1[0]=SF*P1[0];\ + p1[1]=SF*P1[1];\ + p2[0]=SF*P2[0];\ + p2[1]=SF*P2[1];\ + det1 = p1[0]*q[1]-q[0]*p1[1];\ + gama1 = (p1[0]*r[1]-r[0]*p1[1])*det1;\ + alpha1 = (r[0]*q[1] - q[0]*r[1])*det1;\ + alpha1_legal = (alpha1>=0) && (alpha1<=(det1*det1) && (det1!=0));\ + det2 = p2[0]*q[1] - q[0]*p2[1];\ + alpha2 = (r[0]*q[1] - q[0]*r[1]) *det2;\ + gama2 = (p2[0]*r[1] - r[0]*p2[1]) * det2;\ + alpha2_legal = (alpha2>=0) && (alpha2<=(det2*det2) && (det2 !=0));\ + det3=det2-det1;\ + gama3=((p2[0]-p1[0])*(r[1]-p1[1]) - (r[0]-p1[0])*(p2[1]-p1[1]))*det3;\ + if (alpha1_legal)\ +{\ + if (alpha2_legal)\ +{\ + if ( ((gama1<=0) && (gama1>=-(det1*det1))) || ((gama2<=0) && (gama2>=-(det2*det2))) || (gama1*gama2<0)) return 12;\ +}\ + else\ +{\ + if ( ((gama1<=0) && (gama1>=-(det1*det1))) || ((gama3<=0) && (gama3>=-(det3*det3))) || (gama1*gama3<0)) return 13;\ +}\ +}\ + else\ + if (alpha2_legal)\ +{\ + if ( ((gama2<=0) && (gama2>=-(det2*det2))) || ((gama3<=0) && (gama3>=-(det3*det3))) || (gama2*gama3<0)) return 23;\ +}\ + return 0;\ +} + + +#define FABS(x) ((x)>=0?(x):-(x)) /* implement as is fastest on your machine */ + +#define EQU(V0, V1) (FABS((V0) - (V1)) < 1e-6) +#define ZERO(V) (FABS(V) < 1e-6) + +//main procedure + +int tr_tri_intersect3D(float *C1, float *P1, float *P2, float *D1, float *Q1, float *Q2) +{ + float t[3], p1[3], p2[3], r[3], r4[3]; + float beta1, beta2, beta3; + float gama1, gama2, gama3; + float det1, det2, det3; + float dp0, dp1, dp2; + float dq1, dq2, dq3, dr, dr3; + float alpha1, alpha2; + bool alpha1_legal, alpha2_legal; + float SF; + bool beta1_legal, beta2_legal; + + myVmV(r, D1, C1); + // determinant computation + dp0 = P1[1] * P2[2] - P2[1] * P1[2]; + dp1 = P1[0] * P2[2] - P2[0] * P1[2]; + dp2 = P1[0] * P2[1] - P2[0] * P1[1]; + dq1 = Q1[0] * dp0 - Q1[1] * dp1 + Q1[2] * dp2; + dq2 = Q2[0] * dp0 - Q2[1] * dp1 + Q2[2] * dp2; + dr = -r[0] * dp0 + r[1] * dp1 - r[2] * dp2; + + + + beta1 = dr*dq2; // beta1, beta2 are scaled so that beta_i=beta_i*dq1*dq2 + beta2 = dr*dq1; + beta1_legal = (beta2 >= 0) && (beta2 <= dq1*dq1) && (dq1 != 0); + beta2_legal = (beta1 >= 0) && (beta1 <= dq2*dq2) && (dq2 != 0); + + dq3 = dq2 - dq1; + dr3 = +dr - dq1; // actually this is -dr3 + + + if (ZERO(dq1) && ZERO(dq2)) + { + if (dr != 0) return 0; // triangles are on parallel planes + else + { // triangles are on the same plane + float C2[3], C3[3], D2[3], D3[3], N1[3]; + // We use the coplanar test of Moller which takes the 6 vertices and 2 normals + //as input. + myVpV(C2, C1, P1); + myVpV(C3, C1, P2); + myVpV(D2, D1, Q1); + myVpV(D3, D1, Q2); + CROSS(N1, P1, P2); + return coplanar_tri_tri(N1, C1, C2, C3, D1, D2, D3); + } + } + + else if (!beta2_legal && !beta1_legal) return 0;// fast reject-all vertices are on + // the same side of the triangle plane + + else if (beta2_legal && beta1_legal) //beta1, beta2 + { + SF = dq1*dq2; + sVpsV_2(t, beta2, Q2, (-beta1), Q1); + } + + else if (beta1_legal && !beta2_legal) //beta1, beta3 + { + SF = dq1*dq3; + beta1 = beta1 - beta2; // all betas are multiplied by a positive SF + beta3 = dr3*dq1; + sVpsV_2(t, (SF - beta3 - beta1), Q1, beta3, Q2); + } + + else if (beta2_legal && !beta1_legal) //beta2, beta3 + { + SF = dq2*dq3; + beta2 = beta1 - beta2; // all betas are multiplied by a positive SF + beta3 = dr3*dq2; + sVpsV_2(t, (SF - beta3), Q1, (beta3 - beta2), Q2); + Q1 = Q2; + beta1 = beta2; + } + sVpsV_2(r4, SF, r, beta1, Q1); + seg_collide3(t, r4); // calculates the 2D intersection + return 0; +} + + + + +/* this edge to edge test is based on Franlin Antonio's gem: +"Faster Line Segment Intersection", in Graphics Gems III, +pp. 199-202 */ + +int EDGE_EDGE_TEST(float *V0, float *V1, float *U0, float *U1, int i0, int i1, float Ax, float Ay) +{ + float Bx, By, Cx, Cy, e, d, f; + Bx = U0[i0] - U1[i0]; + By = U0[i1] - U1[i1]; + Cx = V0[i0] - U0[i0]; + Cy = V0[i1] - U0[i1]; + + // if edges are adjacent + if (EQU(V0[i0], U0[i0]) && EQU(V0[i1], U0[i1])) + return 0; + if (EQU(V0[i0], U1[i0]) && EQU(V0[i1], U1[i1])) + return 0; + if (EQU(V1[i0], U0[i0]) && EQU(V1[i1], U0[i1])) + return 0; + if (EQU(V1[i0], U1[i0]) && EQU(V1[i1], U1[i1])) + return 0; + + f = Ay*Bx - Ax*By; + d = By*Cx - Bx*Cy; + if ((f > 0 && d >= 0 && d <= f) || (f<0 && d <= 0 && d >= f)) + { + e = Ax*Cy - Ay*Cx; + if (f>0) + { + if (e >= 0 && e <= f) return 1; + } + else + { + if (e <= 0 && e >= f) return 1; + } + } + return 0; +} + +int EDGE_AGAINST_TRI_EDGES(float *V0, float *V1, float *U0, float *U1, float *U2, int i0, int i1) +{ + float Ax,Ay; + Ax=V1[i0]-V0[i0]; + Ay=V1[i1]-V0[i1]; + /* test edge U0,U1 against V0,V1 */ + if (EDGE_EDGE_TEST(V0, V1, U0, U1, i0, i1, Ax, Ay) != 0) + return 1; + /* test edge U1,U2 against V0,V1 */ + if (EDGE_EDGE_TEST(V0, V1, U1, U2, i0, i1, Ax, Ay) != 0) + return 1; + /* test edge U2,U1 against V0,V1 */ + if (EDGE_EDGE_TEST(V0, V1, U2, U0, i0, i1, Ax, Ay) != 0) + return 1; + return 0; +} + +#define POINT_IN_TRI(V0,U0,U1,U2) \ +{ \ + float a,b,c,d0,d1,d2; \ + /* is T1 completly inside T2? */ \ + /* check if V0 is inside tri(U0,U1,U2) */ \ + a=U1[i1]-U0[i1]; \ + b=-(U1[i0]-U0[i0]); \ + c=-a*U0[i0]-b*U0[i1]; \ + d0=a*V0[i0]+b*V0[i1]+c; \ + \ + a=U2[i1]-U1[i1]; \ + b=-(U2[i0]-U1[i0]); \ + c=-a*U1[i0]-b*U1[i1]; \ + d1=a*V0[i0]+b*V0[i1]+c; \ + \ + a=U0[i1]-U2[i1]; \ + b=-(U0[i0]-U2[i0]); \ + c=-a*U2[i0]-b*U2[i1]; \ + d2=a*V0[i0]+b*V0[i1]+c; \ + if(d0*d1>0.0) \ +{ \ + if(d0*d2>0.0) return 1; \ +} \ +} + +//This procedure testing for intersection between coplanar triangles is taken +// from Tomas Moller's +//"A Fast Triangle-Triangle Intersection Test",Journal of Graphics Tools, 2(2), 1997 +int coplanar_tri_tri(float N[3], float V0[3], float V1[3], float V2[3], + float U0[3], float U1[3], float U2[3]) +{ + float A[3]; + short i0, i1; + /* first project onto an axis-aligned plane, that maximizes the area */ + /* of the triangles, compute indices: i0,i1. */ + A[0] = FABS(N[0]); + A[1] = FABS(N[1]); + A[2] = FABS(N[2]); + if (A[0]>A[1]) + { + if (A[0]>A[2]) + { + i0 = 1; /* A[0] is greatest */ + i1 = 2; + } + else + { + i0 = 0; /* A[2] is greatest */ + i1 = 1; + } + } + else /* A[0]<=A[1] */ + { + if (A[2]>A[1]) + { + i0 = 0; /* A[2] is greatest */ + i1 = 1; + } + else + { + i0 = 0; /* A[1] is greatest */ + i1 = 2; + } + } + + /* test all edges of triangle 1 against the edges of triangle 2 */ + if (EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2, i0, i1) != 0) + return 1; + if (EDGE_AGAINST_TRI_EDGES(V1, V2, U0, U1, U2, i0, i1) != 0) + return 1; + if (EDGE_AGAINST_TRI_EDGES(V2, V0, U0, U1, U2, i0, i1) != 0) + return 1; + + /* finally, test if tri1 is totally contained in tri2 or vice versa */ + POINT_IN_TRI(V0, U0, U1, U2); + POINT_IN_TRI(U0, V0, V1, V2); + + return 0; +} + +int MH::tr_tri_intersect3D(const Vec3 & C1, const Vec3 & P1, const Vec3 & P2, const Vec3 & D1, const Vec3 & Q1, const Vec3 & Q2) +{ + float *fC1 = (float *)(Vec3 *)&C1; + float *fP1 = (float *)(Vec3 *)&P1; + float *fP2 = (float *)(Vec3 *)&P2; + float *fD1 = (float *)(Vec3 *)&D1; + float *fQ1 = (float *)(Vec3 *)&Q1; + float *fQ2 = (float *)(Vec3 *)&Q2; + return ::tr_tri_intersect3D(fC1, fP1, fP2, fD1, fQ1, fQ2); +} + diff --git a/src/vec.h b/src/vec.h new file mode 100644 index 0000000..09bfc0f --- /dev/null +++ b/src/vec.h @@ -0,0 +1,139 @@ +#ifndef VEC_H +#define VEC_H + +class Vec2 +{ +public: + float x, y; + +public: + Vec2() : x(0.0f), y(0.0f) {} + Vec2(float _x, float _y) : x(_x), y(_y) {} + Vec2(const Vec2 & v) : x(v.x), y(v.y) {} +}; + +class Vec3 +{ +public: + float x, y, z; + +public: + Vec3() : x(0.0f), y(0.0f), z(0.0f) {} + Vec3(float _x, float _y, float _z) : x(_x), y(_y), z(_z) {} + Vec3(const Vec3 & v) : x(v.x), y(v.y), z(v.z) {} + +public: + Vec3 normalize() + { + float l = sqrtf(x * x + y * y + z * z); + if (l != 0.0f) + { + x /= l; + y /= l; + z /= l; + } + return *this; + } + + Vec3 cross(const Vec3 & vec) const + { + return Vec3(y * vec.z - z * vec.y, z * vec.x - x * vec.z, x * vec.y - y * vec.x); + } + + float dot(const Vec3 & vec) const + { + return x * vec.x + y * vec.y + z * vec.z; + } + + Vec3 interp(const Vec3 & vec, float coeff) const + { + float coeff_1 = 1.0f - coeff; + return Vec3(x * coeff_1 + vec.x * coeff, y * coeff_1 + vec.y * coeff, z * coeff_1 + vec.z * coeff); + } + + float squaredDist() const + { + return x * x + y * y + z * z; + } + + friend Vec3 operator + (const Vec3 & u, const Vec3 & v) + { + return Vec3(u.x + v.x, u.y + v.y, u.z + v.z); + } + + friend Vec3 operator - (const Vec3 & u, const Vec3 & v) + { + return Vec3(u.x - v.x, u.y - v.y, u.z - v.z); + } + + friend Vec3 operator * (const Vec3 & u, const Vec3 & v) + { + return Vec3(u.x * v.x, u.y * v.y, u.z * v.z); + } + + friend Vec3 operator / (const Vec3 & u, const Vec3 & v) + { + return Vec3(u.x / v.x, u.y / v.y, u.z / v.z); + } + + friend bool operator == (const Vec3 & u, const Vec3 & v) + { + return u.x == v.x && u.y == v.y && u.z == v.z; + } + + friend bool operator != (const Vec3 & u, const Vec3 & v) + { + return u.x != v.x || u.y != v.y || u.z != v.z; + } + +public: + static const Vec3 Zero() + { + return Vec3(); + } +}; + +class Vec4 +{ +public: + float x, y, z, w; + +public: + Vec4() : x(0.0f), y(0.0f), z(0.0f), w(0.0f) {} + Vec4(float _x, float _y, float _z, float _w) : x(_x), y(_y), z(_z), w(_w) {} + Vec4(const Vec4 & v) : x(v.x), y(v.y), z(v.z), w(v.w) {} + +public: + float dot(const Vec4 & vec) const + { + return x * vec.x + y * vec.y + z * vec.z + w * vec.w; + } + + Vec4 interp(const Vec4 & vec, float coeff) const + { + float coeff_1 = 1.0f - coeff; + return Vec4(x * coeff_1 + vec.x * coeff, y * coeff_1 + vec.y * coeff, z * coeff_1 + vec.z * coeff, w * coeff_1 + vec.w * coeff); + } + + friend Vec4 operator + (const Vec4 & u, const Vec4 & v) + { + return Vec4(u.x + v.x, u.y + v.y, u.z + v.z, u.w + v.w); + } + + friend Vec4 operator - (const Vec4 & u, const Vec4 & v) + { + return Vec4(u.x - v.x, u.y - v.y, u.z - v.z, u.w - v.w); + } + + friend Vec4 operator * (const Vec4 & u, const Vec4 & v) + { + return Vec4(u.x * v.x, u.y * v.y, u.z * v.z, u.w * v.w); + } + + friend Vec4 operator / (const Vec4 & u, const Vec4 & v) + { + return Vec4(u.x / v.x, u.y / v.y, u.z / v.z, u.w / v.w); + } +}; + +#endif // of VEC_H