diff --git a/src/fgt/filter_csg/filter_csg.cpp b/src/fgt/filter_csg/filter_csg.cpp index 6863e4e66..40bb11505 100644 --- a/src/fgt/filter_csg/filter_csg.cpp +++ b/src/fgt/filter_csg/filter_csg.cpp @@ -23,9 +23,14 @@ #include #include + #include "filter_csg.h" -#include "fixed.h" +#include +#include + +#include +#include "gmpfrac.h" #include "intercept.h" using namespace std; @@ -110,12 +115,26 @@ bool FilterCSG::applyFilter(QAction *filter, MeshDocument &md, RichParameterSet firstMesh->updateDataMask(MeshModel::MM_FACENORMAL); secondMesh->updateDataMask(MeshModel::MM_FACENORMAL); - typedef CMeshO::ScalarType scalar; - typedef Intercept,scalar> intercept; - const float d = par.getFloat("Delta"); + typedef float scalar; + //typedef Intercept,scalar> intercept; + typedef Intercept intercept; + const scalar d = par.getFloat("Delta"); const Point3f delta(d, d, d); - InterceptVolume v(InterceptSet3(firstMesh->cm, delta)), - tmp(InterceptSet3(secondMesh->cm, delta)); + Log(0, "First volume"); + InterceptVolume v = InterceptSet3(firstMesh->cm, delta); + + ofstream out1("/Users/ranma42/Desktop/v1.intercepts.txt"); + out1 << v; + out1.close(); + + Log(0, "Second volume"); + InterceptVolume tmp = InterceptSet3(secondMesh->cm, delta); + + ofstream out2("/Users/ranma42/Desktop/v2.intercepts.txt"); + out2 << tmp; + out2.close(); + + Log(0, "Volumes rasterized"); MeshModel *mesh; switch(par.getEnum("Operator")){ @@ -139,6 +158,17 @@ bool FilterCSG::applyFilter(QAction *filter, MeshDocument &md, RichParameterSet return true; } + ofstream out3("/Users/ranma42/Desktop/out.intercepts.txt"); + out3 << v; + out3.close(); + + Log("[MARCHING CUBES] Building mesh..."); + typedef vcg::intercept::Walker MyWalker; + typedef vcg::tri::ExtendedMarchingCubes MyExtendedMarchingCubes; + typedef vcg::tri::MarchingCubes MyMarchingCubes; + MyWalker walker; + MyMarchingCubes mc(mesh->cm, walker); + walker.BuildMesh(mesh->cm, v, mc); } return true; diff --git a/src/fgt/filter_csg/filter_csg.pro b/src/fgt/filter_csg/filter_csg.pro index 34a9cbc0e..36b26bcd1 100644 --- a/src/fgt/filter_csg/filter_csg.pro +++ b/src/fgt/filter_csg/filter_csg.pro @@ -1,13 +1,20 @@ include (../../shared.pri) +HEADERS += filter_csg.h \ + bigint.h \ + fixed.h \ + intercept.h \ + fraction.h \ + gmpfrac.h -HEADERS += filter_csg.h \ - bigint.h \ - fixed.h \ - intercept.h +SOURCES += filter_csg.cpp +TARGET = filter_csg +TEMPLATE = lib +QT += opengl +CONFIG += plugin -SOURCES += filter_csg.cpp +QMAKE_INCDIR += gmp-5.0.1/ -TARGET = filter_csg -TEMPLATE = lib -QT += opengl -CONFIG += plugin +#QMAKE_LIBS += -lgmpxx -lgmp +#QMAKE_LIBDIR += gmp-5.0.1/.libs/ + +OBJECTS += gmp-5.0.1/.libs/libgmpxx.a gmp-5.0.1/.libs/libgmp.a diff --git a/src/fgt/filter_csg/fraction.h b/src/fgt/filter_csg/fraction.h new file mode 100644 index 000000000..0c6d2c965 --- /dev/null +++ b/src/fgt/filter_csg/fraction.h @@ -0,0 +1,49 @@ +#ifndef BIGFRAC_H +#define BIGFRAC_H + +namespace vcg { + namespace math { + template + class fraction { + public: + fraction() : num(0), den(1) { } + fraction(T n) : num(n), den(1) {} + fraction(T n, T d) : num(n), den(d) {} + fraction(const fraction &from) : num(from.num), den(from.den) {} + ~fraction() {} + + fraction& operator=(const fraction &from) { num = from.num; den = from.den; return *this; } + fraction& operator+=(const fraction &b) { + num = num * y.den + y.num * den; + den *= y.den; + return *this; + } + + inline fraction& operator-=(const fraction &y) { + num = num * y.den - y.num * den; + den *= y.den; + return *this; + } + + inline fraction& operator*=(const fraction &y) { num *= y.num; den *= y.den; return *this; } + inline fraction& operator/=(const fraction &y) { num *= y.den; den *= y.num; return *this; } + + inline fraction operator+(const fraction &y) const { return fraction(*this) += y; } + inline fraction operator-(const fraction &y) const { return fraction(*this) -= y; } + inline fraction operator*(const fraction &y) const { return fraction(*this) *= y; } + inline fraction operator/(const fraction &y) const { return fraction(*this) /= y; } + + inline bool operator==(const fraction &y) const { return num * y.den == den * y.num; } + inline bool operator!=(const fraction &y) const { return num * y.den != den * y.num; } + inline bool operator<=(const fraction &y) const { return num * y.den <= den * y.num; } + inline bool operator>=(const fraction &y) const { return num * y.den >= den * y.num; } + inline bool operator<(const fraction &y) const { return num * y.den < den * y.num; } + inline bool operator>(const fraction &y) const { return num * y.den > den * y.num; } + + private: + T num, den; + }; + } +} + +#endif // BIGFRAC_H diff --git a/src/fgt/filter_csg/gmpfrac.h b/src/fgt/filter_csg/gmpfrac.h new file mode 100644 index 000000000..1676effb5 --- /dev/null +++ b/src/fgt/filter_csg/gmpfrac.h @@ -0,0 +1,60 @@ +#ifndef BIGFRAC_H +#define BIGFRAC_H + +#include +#include + +namespace vcg { + namespace math { + class fraction { + public: + fraction() : f(0) { } + fraction(int n) : f(n, 1) {} + fraction(int n, int d) : f(n, d) { f.canonicalize(); } + fraction(const fraction &from) : f(from.f) {} + ~fraction() {} + + inline fraction& operator=(const fraction &y) { f = y.f; return *this; } + inline fraction& operator+=(const fraction &y) { f += y.f; return *this; } + inline fraction& operator-=(const fraction &y) { f -= y.f; return *this; } + inline fraction& operator*=(const fraction &y) { f *= y.f; return *this; } + inline fraction& operator/=(const fraction &y) { f /= y.f; return *this; } + + inline fraction operator+(const fraction &y) const { return fraction(*this) += y; } + inline fraction operator-(const fraction &y) const { return fraction(*this) -= y; } + inline fraction operator*(const fraction &y) const { return fraction(*this) *= y; } + inline fraction operator/(const fraction &y) const { return fraction(*this) /= y; } + + inline bool operator==(const fraction &y) const { return f == y.f; } + inline bool operator!=(const fraction &y) const { return f != y.f; } + inline bool operator<=(const fraction &y) const { return f <= y.f; } + inline bool operator>=(const fraction &y) const { return f >= y.f; } + inline bool operator<(const fraction &y) const { return f < y.f; } + inline bool operator>(const fraction &y) const { return f > y.f; } + + inline double toFloat() const { return f.get_d(); } + + friend std::ostream& operator<<(std::ostream &out, const fraction &x) { + //return out << x.f.get_d() << "[" << x.f << "]"; + return out << x.f; + } + + friend inline long floor(const fraction &x) { + mpz_class q; + mpz_fdiv_q (q.get_mpz_t(), x.f.get_num_mpz_t(), x.f.get_den_mpz_t()); + return q.get_si(); + } + + friend inline long ceil(const fraction &x) { + mpz_class q; + mpz_cdiv_q (q.get_mpz_t(), x.f.get_num_mpz_t(), x.f.get_den_mpz_t()); + return q.get_si(); + } + + private: + mpq_class f; + }; + } +} + +#endif // BIGFRAC_H diff --git a/src/fgt/filter_csg/intercept.h b/src/fgt/filter_csg/intercept.h index 268e43529..253f8a9f6 100644 --- a/src/fgt/filter_csg/intercept.h +++ b/src/fgt/filter_csg/intercept.h @@ -8,67 +8,184 @@ #include "fixed.h" +#undef OLD + +#define p2print(point) ((point).X()) << ", " << ((point).Y()) +#define p3print(point) p2print(point) << ", " << ((point).Z()) + namespace vcg { namespace intercept { using namespace vcg::math; + + const bool debugRaster = true; + ofstream debugRasterOut("/Users/ranma42/Desktop/out.raster.txt"); + ofstream debugInOut("/Users/ranma42/Desktop/out.in.txt"); template class Intercept { - typedef vcg::Point3<_scalar> Point3x; - - private: - Point3x norm; - _scalar quality; - public: + typedef vcg::Point3<_scalar> Point3x; typedef _dist_type DistType; typedef _scalar Scalar; - DistType dist; + private: + Point3x _norm; + _dist_type _dist; + _scalar _sort_norm; + _scalar _quality; + + public: + inline Intercept () { } + + inline Intercept (const DistType &dist, const Scalar &sort_norm) : _dist(dist), _sort_norm(sort_norm) { } - inline Intercept () {}; + inline Intercept (const DistType &dist, const Scalar &sort_norm, const Point3x &norm, _scalar quality) : + _norm(norm), _dist(dist), _sort_norm(sort_norm), _quality(quality) { } + + inline Intercept operator -() const { return Intercept(_dist, -_sort_norm, -_norm, _quality); } + +#if OLD + inline bool operator <(const Intercept &other) const { return _dist < other._dist || (_dist == other.dist() && _sort_norm < other._sort_norm); } +#else + inline bool operator <(const Intercept &other) const { return _dist < other._dist; } +#endif + //inline bool operator <(const DistType &other) const { return _dist < other; } + + //inline friend bool operator<(const fraction &_dist, const Intercept &other) { return _dist < other._dist; } - inline Intercept (DistType dist, const Point3x &norm, Scalar quality) : norm(norm), quality(quality), dist(dist) {}; + inline const DistType& dist() const { return _dist; } + + inline const Scalar& sort_norm() const { return _sort_norm; } - inline Intercept operator -() const { return Intercept(dist, -norm, quality); } - - inline bool operator <(const Intercept &other) const { return dist < other.dist; } - - inline bool operator <(DistType other) const { return dist < other; } + inline const Scalar& quality() const { return _quality; } + + inline const Point3x& norm() const { return _norm; } + + friend ostream& operator<<(ostream &out, const Intercept &x) { + out << "Intercept[" << x._dist << "[" << x._sort_norm << "], (" << p3print(x._norm) << "), " << x._quality << "]"; + return out; + } }; - + template class InterceptRay { typedef typename InterceptType::DistType DistType; - + typedef std::vector ContainerType; - + public: - inline InterceptRay() {}; - + inline InterceptRay() { } + inline InterceptRay(const ContainerType &set) : v(set) { std::sort(v.begin(), v.end()); - assert(v.size() % 2 == 0); + cleanup(); + assert (isValid()); + } + +#if OLD + inline void cleanup() { + v.resize(v.size()); + } +#else + inline void cleanup() { + ContainerType newv; + newv.reserve(v.size()); + + typename ContainerType::const_iterator from = v.begin(); + while (from != v.end()) { + typename ContainerType::const_iterator to = from + 1; + while (to != v.end() && from->dist() == to->dist()) + to++; + + size_t pcond = false, ncond = false; + typename ContainerType::const_iterator pos, neg; + while (from != to) { + if (from->sort_norm() > 0) { + pcond ^= true; + pos = from; + } else { + ncond ^= true; + neg = from; + } + from++; + } + + if (pcond != ncond) + newv.push_back(pcond ? *pos : *neg); + } + newv.resize(newv.size()); + v = newv; + } +#endif + + inline bool isValid() const { + if (v.empty()) + return true; + + if (v.size() % 2 != 0) { + cerr << "Not a solid! (size: " << v.size() << ")" << endl; + return false; + } + + typename ContainerType::const_iterator curr = v.begin(); + typename ContainerType::const_iterator next = curr+1; + while (next != v.end()) { + if (!(*curr < *next)) { + cerr << "Not sorted! (" << *curr << " >= " << *next << ")" << endl; + return false; + } + curr = next; + next++; + } + + return true; } - inline const InterceptType& GetIntercept(DistType s) const { - // TODO: s==p.dist - typename ContainerType::const_iterator p = std::lower_bound(v.begin(), v.end(), s); + inline const InterceptType& GetIntercept(const DistType& s, const DistType& u) const { + typename ContainerType::const_iterator p = std::lower_bound(v.begin(), v.end(), InterceptType(s,-2)); + assert (IsInExt(s) != IsInExt(u)); assert (p != v.end()); + // TODO + //if (p->dist() == s && p->sort_norm() < 0) + //if (p->dist() == s) + // p++; + assert (InterceptType(s,-2) < *p && *p < InterceptType(u,2)); return *p; } - inline bool IsIn(DistType s) const { - // TODO: s==p.dist - typename ContainerType::const_iterator p = std::lower_bound(v.begin(), v.end(), s); + inline int IsInExtDebug(const DistType& s) const { + typename ContainerType::const_iterator p = std::lower_bound(v.begin(), v.end(), InterceptType(s,-2)); + if (p == v.end()) + return -1; + else if (p->dist() == s) + return 0; + else + return 2*((p - v.begin()) & 1) - 1; + } + + inline int IsInExt(const DistType& s) const { + int r = IsInExtDebug(s); + debugInOut << *this << " contains " << s << "? " << r << endl; + return r; + } + + inline bool IsInOld(const DistType& s) const { + typename ContainerType::const_iterator p = std::lower_bound(v.begin(), v.end(), InterceptType(s,-2)); + //if (p != v.end() && p->dist() == s && p->sort_norm() < 0) + if (p != v.end() && p->dist() == s) + //p++; + p = v.end(); + debugInOut << *this << " contains " << s << "? " << (p != v.end() && (p - v.begin()) & 1) << endl; return p != v.end() && (p - v.begin()) & 1; } inline InterceptRay operator &(const InterceptRay & other) const { - ContainerType newv(v.size() + other.v.size()); - typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = v.end(); + typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = other.v.end(); + ContainerType newv; + + newv.reserve(v.size() + other.v.size()); while (i != endi && j != endj) { if (*j < *i) { std::swap(i, j); @@ -89,13 +206,14 @@ namespace vcg { } i += 2; } - newv.resize(newv.size()); return InterceptRay(newv); } inline InterceptRay operator |(const InterceptRay & other) const { - ContainerType newv(v.size() + other.v.size()); - typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = v.end(); + typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = other.v.end(); + ContainerType newv; + + newv.reserve(v.size() + other.v.size()); while (i != endi && j != endj) { if (*j < *i) { std::swap(i, j); @@ -118,13 +236,14 @@ namespace vcg { } newv.insert(newv.end(), i, endi); newv.insert(newv.end(), j, endj); - newv.resize(newv.size()); return InterceptRay(newv); } inline InterceptRay operator -(const InterceptRay & other) const { - ContainerType newv(v.size() + other.v.size()); - typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = v.end(); + typename ContainerType::const_iterator i = v.begin(), j = other.v.begin(), endi = v.end(), endj = other.v.end(); + ContainerType newv; + + newv.reserve(v.size() + other.v.size()); while (i != endi && j != endj) { while (j != endj && !(*i < *(j+1))) // j < J <= i j += 2; @@ -155,10 +274,20 @@ namespace vcg { i += 2; } newv.insert(newv.end(), i, endi); - newv.resize(newv.size()); return InterceptRay(newv); } + friend ostream& operator<<(ostream &out, const InterceptRay &x) { + typename ContainerType::const_iterator i; + typename ContainerType::const_iterator end = x.v.end(); + out << "InterceptRay["; + for (i = x.v.begin(); i != end; ++i) + out << *i; + out << "]"; + assert (x.isValid()); + return out; + } + private: ContainerType v; }; @@ -173,30 +302,41 @@ namespace vcg { public: typedef std::vector > ContainerType; - inline InterceptBeam(const vcg::Box2i &box, const ContainerType &rays) : bbox(box), ray(rays) {}; + inline InterceptBeam(const vcg::Box2i &box, const ContainerType &rays) : bbox(box), ray(rays) { } inline const IRayType& GetInterceptRay (const vcg::Point2i &p) const { assert(bbox.IsIn(p)); vcg::Point2i c = p - bbox.min; + assert(c.X() < ray.size() && c.Y() < ray[c.X()].size()); return ray[c.X()][c.Y()]; } - inline bool IsIn (const vcg::Point2i &p, DistType s) const { - return bbox.IsIn(p) && GetInterceptRay(p).IsIn(s); + inline int IsInExt (const vcg::Point2i &p, const DistType& s) const { + int r; + if (bbox.IsIn(p)) + r = GetInterceptRay(p).IsInExt(s); + else + r = -1; + assert (-1 <= r && r <= 1); + return r; + } + + inline bool IsInOld (const vcg::Point2i &p, DistType s) const { + return bbox.IsIn(p) && GetInterceptRay(p).IsInOld(s); } inline InterceptBeam &operator &=(const InterceptBeam & other) { vcg::Box2i newbbox(bbox); newbbox.Intersect(other.bbox); - for(int i = 0; i < newbbox.DimX(); ++i) { - for(int j = 0; j < newbbox.DimY(); ++j) { + for(int i = 0; i <= newbbox.DimX(); ++i) { + for(int j = 0; j <= newbbox.DimY(); ++j) { vcg::Point2i p = newbbox.min + vcg::Point2i(i,j); ray[i][j] = GetInterceptRay(p) & other.GetInterceptRay(p); } - ray[i].resize(newbbox.DimY()); + ray[i].resize(newbbox.DimY() + 1); } - ray.resize(newbbox.DimX()); + ray.resize(newbbox.DimX() + 1); bbox = newbbox; return *this; } @@ -232,6 +372,19 @@ namespace vcg { return *this; } + friend ostream& operator<<(ostream &out, const InterceptBeam &x) { + out << "InterceptBeam[" << p2print(x.bbox.min) << " - " << p2print(x.bbox.max) << "][" << endl; + + for(int i = x.bbox.min.X(); i <= x.bbox.max.X(); ++i) { + for(int j = x.bbox.min.Y(); j <= x.bbox.max.Y(); ++j) { + vcg::Point2i p(i,j); + out << p2print(p) << ": " << x.GetInterceptRay(p) << endl; + } + } + out << "]"; + return out; + } + private: vcg::Box2i bbox; ContainerType ray; @@ -251,12 +404,13 @@ namespace vcg { public: typedef typename std::vector > ContainerType; - inline InterceptVolume(const Point3x &d, const ContainerType &beams) : delta(d), beam(beams) { assert(beams.size() == 3); }; + inline InterceptVolume(const Box3i &b, const Point3x &d, const ContainerType &beams) : delta(d), bbox(b), beam(beams) { assert(beams.size() == 3); }; inline InterceptVolume & operator &=(const InterceptVolume & other) { assert(checkConsistency(other)); for (int i = 0; i < 3; ++i) beam[i] &= other.beam[i]; + bbox.Intersect(other.bbox); return *this; } @@ -264,6 +418,7 @@ namespace vcg { assert(checkConsistency(other)); for (int i = 0; i < 3; ++i) beam[i] |= other.beam[i]; + bbox.Add(other.bbox); return *this; } @@ -275,22 +430,107 @@ namespace vcg { } template - inline const InterceptType& GetIntercept (const vcg::Point3i &p) const { + inline const InterceptType& GetIntercept (const Point3i &p1, const Point3i &p2) const { + assert(0 <= coord && coord < 3); + assert(p2 == p1 + Point3i(coord == 0, coord == 1, coord == 2)); + assert(IsInExt(p1) != IsInExt(p2)); + + const int c1 = (coord + 1) % 3; + const int c2 = (coord + 2) % 3; + return beam[coord].GetInterceptRay(Point2i(p1.V(c1), p1.V(c2))).GetIntercept(p1.V(coord), p2.V(coord)); + } + + inline const InterceptRay & GetInterceptRay (int coord, const vcg::Point3i &p) const { assert(0 <= coord && coord < 3); const int c1 = (coord + 1) % 3; const int c2 = (coord + 2) % 3; - return beam[coord].GetInterceptRay(vcg::Point2i(p.V(c1), p.V(c2))).GetIntercept(p.V(coord)); + return beam[coord].GetInterceptRay(vcg::Point2i(p.V(c1), p.V(c2))); } - inline bool IsIn (const vcg::Point3i &p) const { + inline int IsInExt (const vcg::Point3i &p) const { + int r[3]; // TODO: simplify me! + for (int i = 0; i < 3; ++i) + r[i] = beam[i].IsInExt(Point2i(p.V((i+1)%3), p.V((i+2)%3)), p.V(i)); + + if (r[0] == 0) + r[0] += r[1] + r[2]; + if (r[1] == 0) + r[1] += r[0] + r[1]; + if (r[2] == 0) + r[2] += r[2] + r[0]; + + if (r[0]>0 && r[1]>0 && r[2]>0) + return 1; + else if ((r[0]<0 && r[1]<0 && r[2]<0) || + r[0]==0 && r[1]==0 && r[2] == 0) + return -1; + + cerr << "Inconsistency: " << p3print(p) << p3print(delta) << endl; + for (int i = 0; i < 3; ++i) { + cerr << beam[i].IsInExt(Point2i(p.V((i+1)%3), p.V((i+2)%3)), p.V(i)); + cerr << ": " << beam[i].GetInterceptRay(Point2i(p.V((i+1)%3), p.V((i+2)%3))) << endl; + } + + assert (false); + return 0; + } + + inline bool IsInOld (const vcg::Point3i &p) const { bool r[3]; // TODO: simplify me! for (int i = 0; i < 3; ++i) - r[i] = beam[i].IsIn(p); - assert(r[0] == r[1] && r[1] == r[2]); - return r[0]; + r[i] = beam[i].IsInOld(Point2i(p.V((i+1)%3), p.V((i+2)%3)), p.V(i)); + if (!(r[0] == r[1] && r[1] == r[2])) { + cerr << "Inconsistency: " << p3print(p) << p3print(delta) << endl; + for (int i = 0; i < 3; ++i) + cerr << r[i] << ": " << beam[i].GetInterceptRay(Point2i(p.V((i+1)%3), p.V((i+2)%3))) << endl; + } + //assert(r[0] == r[1] && r[1] == r[2]); + return r[0] & r[1] & r[2]; + } + + friend ostream& operator<<(ostream &out, const InterceptVolume &x) { + out << "InterceptVolume[" << p3print(x.delta) << "][" << endl; + int coord = 0; + for(typename ContainerType::const_iterator iter = x.beam.begin(); iter != x.beam.end(); ++iter) { + out << *iter << endl; + out << "Beam " << coord << endl; + + for (int i=x.bbox.min[coord]; i<=x.bbox.max[coord]; i+=1) { + out << i << endl; + + for (int k=x.bbox.min[(coord+2)%3]; k<=x.bbox.max[(coord+2)%3]+2; k+=1) + out << '+'; + out << endl; + + for (int j=x.bbox.min[(coord+1)%3]; j<=x.bbox.max[(coord+1)%3]; j+=1) { + out << '+'; + for (int k=x.bbox.min[(coord+2)%3]; k<=x.bbox.max[(coord+2)%3]; k+=1) { + Point3i p(i,j,k); + int in = iter->IsInExt(Point2i(j, k), i); + char c = '?'; + if (in < 0) + c = ' '; + else if (in > 0) + c = '#'; + out << p3print(Point3i(i,j,k)) << " -> " << in << "[" << c << "]" << endl; + //out << c; + } + out << '+' << endl; + } + + for (int k=x.bbox.min[(coord+2)%3]; k class InterceptSet { + typedef std::vector ContainerType; typedef InterceptRay SortedType; public: + inline InterceptSet() { } + inline void AddIntercept (const InterceptType &x) { v.push_back(x); } inline operator SortedType() const { return SortedType(v); } + friend ostream& operator<<(ostream &out, const InterceptSet &x) { + typename ContainerType::const_iterator i; + typename ContainerType::const_iterator end = x.v.end(); + out << "InterceptSet["; + for (i = x.v.begin(); i != end; ++i) + out << *i << endl; + out << "]"; + return out; + } + private: - std::vector v; + ContainerType v; }; template class InterceptSet1 { + typedef std::vector > ContainerType; typedef std::vector > SortedType; public: - inline InterceptSet1() {}; + inline InterceptSet1() { } inline operator SortedType() const { return SortedType(set.begin(), set.end()); } + inline void resize(size_t size) { set.resize(size); } + inline void AddIntercept (const int i, const InterceptType &x) { + assert(i >= 0 && i < set.size()); set[i].AddIntercept(x); } + friend ostream& operator<<(ostream &out, const InterceptSet1 &x) { + typename ContainerType::const_iterator i; + typename ContainerType::const_iterator end = x.set.end(); + out << "InterceptSet1["; + for (i = x.set.begin(); i != end; ++i) + out << *i << endl; + out << "]InterceptSet1"; + return out; + } + private: - std::vector > set; + ContainerType set; }; template class InterceptSet2 { + typedef std::vector > ContainerType; + typedef std::vector > > NewContainerType; typedef InterceptBeam SortedType; public: - inline InterceptSet2(const vcg::Box2i &box) : bbox(box) { set.reserve(bbox.DimX()); }; + inline InterceptSet2(const Box2i &box) : bbox(box), set(box.DimX() + 1) { + typename ContainerType::iterator i; + typename ContainerType::iterator end = set.end(); + for (i = set.begin(); i != end; ++i) + i->resize(box.DimY() + 1); + } inline operator SortedType() const { return SortedType(bbox, typename SortedType::ContainerType(set.begin(), set.end())); } inline void AddIntercept (const vcg::Point2i &p, const InterceptType &x) { assert(bbox.IsIn(p)); vcg::Point2i c = p - bbox.min; + assert(c.X() >= 0 && c.X() < set.size()); set[c.X()].AddIntercept(c.Y(), x); } + friend ostream& operator<<(ostream &out, const InterceptSet2 &x) { + typename ContainerType::const_iterator i; + typename ContainerType::const_iterator end = x.set.end(); + out << "InterceptSet2["; + for (i = x.set.begin(); i != end; ++i) + out << *i << endl; + out << "]InterceptSet2"; + return out; + } + private: - vcg::Box2i bbox; - std::vector > set; + Box2i bbox; + ContainerType set; }; template @@ -357,6 +642,7 @@ namespace vcg { typedef vcg::Point3 Point3x; typedef InterceptSet2 ISet2Type; typedef InterceptVolume SortedType; + typedef std::vector ContainerType; template void RasterFace(const Point3dt &v0, const Point3dt &v1, const Point3dt &v2, @@ -369,18 +655,97 @@ namespace vcg { const Point3dt d21 = v2 - v1; const Point3dt d02 = v0 - v2; + const DistType det0 = d21[crd2] * d02[crd1] - d21[crd1] * d02[crd2]; + const DistType det1 = d21[crd0] * d02[crd2] - d21[crd2] * d02[crd0]; + const DistType det2 = d21[crd1] * d02[crd0] - d21[crd0] * d02[crd1]; + + if (debugRaster) { + debugRasterOut << "Face [" << crd0 << "] "; + debugRasterOut << "(" << p3print(v0) << "), "; + debugRasterOut << "(" << p3print(v1) << "), "; + debugRasterOut << "(" << p3print(v2) << ")" << endl; + + debugRasterOut << p3print(d10) << endl; + } for(int x = ibox.min[crd1]; x <= ibox.max[crd1]; ++x) { for(int y = ibox.min[crd2]; y <= ibox.max[crd2]; ++y) { - /* - scalar n0 = (x-v0[crd1])*d10[crd2] - (y-v0[crd2])*d10[crd1]; - scalar n1 = (x-v1[crd1])*d21[crd2] - (y-v1[crd2])*d21[crd1]; - scalar n2 = (x-v2[crd1])*d02[crd2] - (y-v2[crd2])*d02[crd1]; + DistType n0 = (v1[crd1]-x)*d21[crd2] - (v1[crd2]-y)*d21[crd1]; + DistType n1 = (v2[crd1]-x)*d02[crd2] - (v2[crd2]-y)*d02[crd1]; + DistType n2 = (v0[crd1]-x)*d10[crd2] - (v0[crd2]-y)*d10[crd1]; - if( (n0>-EPS && n1>-EPS && n2>-EPS) || - (n0< EPS && n1< EPS && n2< EPS)) - set[crd0].AddIntercept(vcg::Point2i(x, y), - InterceptType((dist - x*norm[crd1] - y*norm[crd2]) / norm[crd0], norm, quality)); - */ + if (debugRaster) + debugRasterOut << x << "," << y << ": " << p3print(Point3dt(n0,n1,n2)) << " -> "; + + /* + if (n0>0 || n1>0 || n2>0) { + if (n0 == 0 && + (d21[crd1] > 0 || d21[crd2] < 0)) + n0 = 1; + + if (n1 == 0 && + (d02[crd1] > 0 || d02[crd2] < 0)) + n1 = 1; + + if (n2 == 0 && + (d10[crd1] > 0 || d10[crd2] < 0)) + n2 = 1; + } else { + if (n0 == 0 && + (d21[crd1] < 0 || d21[crd2] > 0)) + n0 = -1; + + if (n1 == 0 && + (d02[crd1] < 0 || d02[crd2] > 0)) + n1 = -1; + + if (n2 == 0 && + (d10[crd1] < 0 || d10[crd2] > 0)) + n2 = -1; + } + */ + + if (crd1 > crd2) { + if (n2 == 0) + n2 = d10[crd1]; + if (n2 == 0) + n2 -= d10[crd2]; + + if (n0 == 0) + n0 = d21[crd1]; + if (n0 == 0) + n0 -= d21[crd2]; + + if (n1 == 0) + n1 = d02[crd1]; + if (n1 == 0) + n1 -= d02[crd2]; + } else { + if (n2 == 0) + n2 -= d10[crd2]; + if (n2 == 0) + n2 = d10[crd1]; + + if (n0 == 0) + n0 -= d21[crd2]; + if (n0 == 0) + n0 = d21[crd1]; + + if (n1 == 0) + n1 -= d02[crd2]; + if (n1 == 0) + n1 = d02[crd1]; + } + + if (debugRaster) + debugRasterOut << p3print(Point3dt(n0,n1,n2)) << endl; + + if((n0>0 && n1>0 && n2>0) || (n0<0 && n1<0 && n2<0)) { + DistType d = (v0[crd2] - y) * det2 + (v0[crd1] - x) * det1; + d /= det0; + d += v0[crd0]; + assert(d >= ibox.min[crd0] && d <= ibox.max[crd0]); + set[crd0].AddIntercept(vcg::Point2i(x, y), InterceptType(d, norm[crd0], norm, quality)); + } } } } @@ -400,71 +765,140 @@ namespace vcg { RasterFace<2>(v0, v1, v2, ibox, norm, quality); } - void ScanFace(const Point3x &v0, const Point3x &v1, const Point3x &v2, - Scalar quality, const Point3x &norm) { - ScanFace(Point3dt(v0.X(), v0.Y(), v0.Z()), - Point3dt(v1.X(), v1.Y(), v1.Z()), - Point3dt(v2.X(), v2.Y(), v2.Z()), - quality, norm); - } - public: template - inline InterceptSet3(const MeshType &m, const Point3x &d) : delta(d) { + inline InterceptSet3(const MeshType &m, const Point3x &d) : delta(d), + bbox(Point3i(floor(m.bbox.min.X() / d.X()) -1, + floor(m.bbox.min.Y() / d.Y()) -1, + floor(m.bbox.min.Z() / d.Z()) -1), + Point3i(ceil(m.bbox.max.X() / d.X()) + 0, + ceil(m.bbox.max.Y() / d.Y()) + 0, + ceil(m.bbox.max.Z() / d.Z()) + 0)) + { const Point3x invDelta(Scalar(1) / delta.X(), Scalar(1) / delta.Y(), Scalar(1) / delta.Z()); - set.push_back(ISet2Type(vcg::Box2i(/*TODO*/))); - set.push_back(ISet2Type(vcg::Box2i(/*TODO*/))); - set.push_back(ISet2Type(vcg::Box2i(/*TODO*/))); + vcg::Box2i xy, yz, zx; + yz.Set(bbox.min.Y(), bbox.min.Z(), bbox.max.Y(), bbox.max.Z()); + zx.Set(bbox.min.Z(), bbox.min.X(), bbox.max.Z(), bbox.max.X()); + xy.Set(bbox.min.X(), bbox.min.Y(), bbox.max.X(), bbox.max.Y()); + + set.push_back(ISet2Type(yz)); + set.push_back(ISet2Type(zx)); + set.push_back(ISet2Type(xy)); typename MeshType::ConstFaceIterator i, end = m.face.end(); - for (i = m.face.begin(); i != end; ++i) - ScanFace (Point3x(i->V(0)->P()).Scale(invDelta), - Point3x(i->V(1)->P()).Scale(invDelta), - Point3x(i->V(2)->P()).Scale(invDelta), + for (i = m.face.begin(); i != end; ++i) { + Point3x v0(i->V(0)->P()), v1(i->V(1)->P()), v2(i->V(2)->P()); + v0.Scale(invDelta); + v1.Scale(invDelta); + v2.Scale(invDelta); + ScanFace (Point3dt(v0.X(), v0.Y(), v0.Z()), + Point3dt(v1.X(), v1.Y(), v1.Z()), + Point3dt(v2.X(), v2.Y(), v2.Z()), 0, Point3x(i->cN()).Scale(delta).Normalize()); + } } - inline operator SortedType() const { return SortedType(delta, typename SortedType::ContainerType(set.begin(), set.end())); } + inline operator SortedType() const { return SortedType(bbox, delta, typename SortedType::ContainerType(set.begin(), set.end())); } + + friend ostream& operator<<(ostream &out, const InterceptSet3 &x) { + typename ContainerType::const_iterator i; + typename ContainerType::const_iterator end = x.set.end(); + out << "InterceptSet3["; + for (i = x.set.begin(); i != end; ++i) + out << *i << endl; + out << "]InterceptSet3"; + return out; + } const Point3x delta; + const Box3i bbox; private: - std::vector set; + ContainerType set; }; - template + template class Walker { typedef typename MeshType::VertexPointer VertexPointer; - typedef Intercept InterceptType; + typedef typename MeshType::CoordType CoordType; public: - inline Walker(MeshType &mesh, const InterceptVolume &volume) : m(mesh), v(volume) {} + template + void BuildMesh(MeshType &mesh, InterceptVolume &volume, EXTRACTOR_TYPE &extractor) + { + //Init(volume); + _volume = &volume; + _mesh = &mesh; + _mesh->Clear(); + vcg::Point3i p1, p2; - const scalar V(int i, int j, int k) const { return scalar(v.IsIn(vcg::Point3i(i, j, k)) ? 1.0 : -1.0); } - - template - void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer p) { - assert(p1 == p2 + vcg::Point3i(coord == 0, coord == 1, coord == 2)); - p = m.vert.end(); - // MeshType::Allocator::AddVertices(m, 1); - p->P() = MeshType::CoordType(p1); - p->P().V(coord) = v.GetIntercept(p1).dist; - assert(p1.V(coord) <= p->P().V(coord) && p->P().V(coord) <= p2.V(coord)); + //Begin(); + extractor.Initialize(); + for (int j=_volume->bbox.min.Y(); j<=_volume->bbox.max.Y(); ++j) + { + for (int i=_volume->bbox.min.X(); i<=_volume->bbox.max.X(); ++i) + for (int k=_volume->bbox.min.Z(); k<=_volume->bbox.max.Z(); ++k) + extractor.ProcessCell(vcg::Point3i(i,j,k), + vcg::Point3i(i+1,j+1,k+1)); + //NextSlice(); + } + extractor.Finalize(); + _volume = NULL; + _mesh = NULL; } - void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer p) { GetIntercept<0>(p1, p2, v); } + const float V(int i, int j, int k) const { return _volume->IsInExt(vcg::Point3i(i, j, k)); } - void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer p) { GetIntercept<1>(p1, p2, v); } + template + void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer& p) { + assert(p2 == p1 + vcg::Point3i(coord == 0, coord == 1, coord == 2)); + assert(_volume->IsInExt(p1) != _volume->IsInExt(p2)); + p = &*vcg::tri::Allocator::AddVertices(*_mesh, 1); + const InterceptType& i = _volume->GetIntercept(p1, p2); + p->P().V(coord) = i.dist().toFloat(); + p->P().V((coord+1)%3) = p1[(coord+1)%3]; + p->P().V((coord+2)%3) = p1[(coord+2)%3]; + if (!(p1.V(coord) <= p->P().V(coord) && p->P().V(coord) <= p2.V(coord))) { + cerr << "Inconsistency: " << p3print(p1) << " <= " << p3print(p->P()) << " <= " << p3print(p2) << endl; + cerr << (_volume->GetInterceptRay(coord, p2)) << endl; + } + assert(p1.V(coord) <= p->P().V(coord) && p->P().V(coord) <= p2.V(coord)); - void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer p) { GetIntercept<2>(p1, p2, v); } + p->P().Scale(_volume->delta); + p->N() = i.norm(); + p->Q() = i.quality(); + } + + void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer& p) { GetIntercept<0>(p1, p2, p); } + + void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer& p) { GetIntercept<1>(p1, p2, p); } + + void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer& p) { GetIntercept<2>(p1, p2, p); } + + bool Exist(const Point3i &p1, const Point3i &p2, VertexPointer& p) { + bool in1 = V(p1.X(), p1.Y(), p1.Z()); + bool in2 = V(p2.X(), p2.Y(), p2.Z()); + if (in1 == in2) + return false; + + Point3i d = p2 - p1; + if (d.X()) + GetXIntercept(p1, p2, p); + else if (d.Y()) + GetYIntercept(p1, p2, p); + else if (d.Z()) + GetZIntercept(p1, p2, p); + + return true; + } private: - MeshType &m; - InterceptVolume v; + InterceptVolume *_volume; + MeshType *_mesh; }; }; };