(Mostly) working CSG plugin.

This commit is contained in:
Paolo Cignoni cignoni 2010-09-14 10:01:42 +00:00
parent 1be50d61ed
commit 5d19aaa9b1
5 changed files with 695 additions and 115 deletions

View File

@ -23,9 +23,14 @@
#include <Qt>
#include <QtGui>
#include "filter_csg.h"
#include "fixed.h"
#include <vcg/complex/trimesh/create/extended_marching_cubes.h>
#include <vcg/complex/trimesh/create/marching_cubes.h>
#include <fstream>
#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<vcg::math::fixed<8,vcg::math::base32>,scalar> intercept;
const float d = par.getFloat("Delta");
typedef float scalar;
//typedef Intercept<vcg::math::fixed<8,vcg::math::base32>,scalar> intercept;
typedef Intercept<fraction,scalar> intercept;
const scalar d = par.getFloat("Delta");
const Point3f delta(d, d, d);
InterceptVolume<intercept> v(InterceptSet3<intercept>(firstMesh->cm, delta)),
tmp(InterceptSet3<intercept>(secondMesh->cm, delta));
Log(0, "First volume");
InterceptVolume<intercept> v = InterceptSet3<intercept>(firstMesh->cm, delta);
ofstream out1("/Users/ranma42/Desktop/v1.intercepts.txt");
out1 << v;
out1.close();
Log(0, "Second volume");
InterceptVolume<intercept> tmp = InterceptSet3<intercept>(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<CMeshO, intercept> MyWalker;
typedef vcg::tri::ExtendedMarchingCubes<CMeshO, MyWalker> MyExtendedMarchingCubes;
typedef vcg::tri::MarchingCubes<CMeshO, MyWalker> MyMarchingCubes;
MyWalker walker;
MyMarchingCubes mc(mesh->cm, walker);
walker.BuildMesh<MyMarchingCubes>(mesh->cm, v, mc);
}
return true;

View File

@ -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

View File

@ -0,0 +1,49 @@
#ifndef BIGFRAC_H
#define BIGFRAC_H
namespace vcg {
namespace math {
template<typename T>
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

View File

@ -0,0 +1,60 @@
#ifndef BIGFRAC_H
#define BIGFRAC_H
#include <iostream>
#include <gmpxx.h>
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

View File

@ -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 <typename _dist_type, typename _scalar>
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<DistType, Scalar> &x) {
out << "Intercept[" << x._dist << "[" << x._sort_norm << "], (" << p3print(x._norm) << "), " << x._quality << "]";
return out;
}
};
template <typename InterceptType>
class InterceptRay
{
typedef typename InterceptType::DistType DistType;
typedef std::vector<InterceptType> 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<InterceptType> &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<std::vector<IRayType > > 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<InterceptType> &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<InterceptBeam<InterceptType> > 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 <const int coord>
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<InterceptType> & 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<InterceptType> &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<x.bbox.max[(coord+2)%3]+2; k+=1)
out << '+';
out << endl;
}
coord++;
}
out << "]";
return out;
}
const Point3x delta;
Box3i bbox;
private:
ContainerType beam;
};
@ -299,53 +539,98 @@ namespace vcg {
template <typename InterceptType>
class InterceptSet
{
typedef std::vector<InterceptType> ContainerType;
typedef InterceptRay<InterceptType> 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<InterceptType> &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<InterceptType> v;
ContainerType v;
};
template <typename InterceptType>
class InterceptSet1
{
typedef std::vector<InterceptSet<InterceptType> > ContainerType;
typedef std::vector<InterceptRay<InterceptType> > 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<InterceptType> &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<InterceptSet<InterceptType> > set;
ContainerType set;
};
template <typename InterceptType>
class InterceptSet2
{
typedef std::vector<InterceptSet1<InterceptType> > ContainerType;
typedef std::vector<std::vector<InterceptSet<InterceptType> > > NewContainerType;
typedef InterceptBeam<InterceptType> 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<InterceptType> &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<InterceptSet1<InterceptType> > set;
Box2i bbox;
ContainerType set;
};
template <typename InterceptType>
@ -357,6 +642,7 @@ namespace vcg {
typedef vcg::Point3<Scalar> Point3x;
typedef InterceptSet2<InterceptType> ISet2Type;
typedef InterceptVolume<InterceptType> SortedType;
typedef std::vector<ISet2Type> ContainerType;
template <const int CoordZ>
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 <class MeshType>
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<InterceptType> &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<ISet2Type> set;
ContainerType set;
};
template <typename MeshType, typename scalar>
template <typename MeshType, typename InterceptType>
class Walker
{
typedef typename MeshType::VertexPointer VertexPointer;
typedef Intercept<double,scalar> InterceptType;
typedef typename MeshType::CoordType CoordType;
public:
inline Walker(MeshType &mesh, const InterceptVolume<InterceptType> &volume) : m(mesh), v(volume) {}
template<typename EXTRACTOR_TYPE>
void BuildMesh(MeshType &mesh, InterceptVolume<InterceptType> &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 <const int coord>
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<MeshType>::AddVertices(m, 1);
p->P() = MeshType::CoordType(p1);
p->P().V(coord) = v.GetIntercept<coord>(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 <const int coord>
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<MeshType>::AddVertices(*_mesh, 1);
const InterceptType& i = _volume->GetIntercept<coord>(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<InterceptType> v;
InterceptVolume<InterceptType> *_volume;
MeshType *_mesh;
};
};
};