dist and coputeRT removed (see distance.h and updateEdges)

This commit is contained in:
ganovelli 2004-05-12 18:49:05 +00:00
parent 61c7d41569
commit 2fe139f9ac
1 changed files with 3 additions and 200 deletions

View File

@ -24,6 +24,9 @@
History
$Log: not supported by cvs2svn $
Revision 1.12 2004/05/12 14:43:36 cignoni
removed warning of unused variables
Revision 1.11 2004/05/12 12:50:20 turini
include color4
@ -1158,206 +1161,6 @@ void Swap ( const int z )
Plane3<ScalarType> plane;
#endif
void ComputeRT()
{
#ifdef __VCGLIB_FACE_RT
// Primo calcolo degli edges
edge[0] = V(1)->P(); edge[0] -= V(0)->P();
edge[1] = V(2)->P(); edge[1] -= V(1)->P();
edge[2] = V(0)->P(); edge[2] -= V(2)->P();
// Calcolo di plane
plane.SetDirection(edge[0]^edge[1]);
plane.SetOffset(plane.Direction() * V(0)->P());
plane.Normalize();
// Calcolo migliore proiezione
ScalarType nx = math::Abs(plane.Direction()[0]);
ScalarType ny = math::Abs(plane.Direction()[1]);
ScalarType nz = math::Abs(plane.Direction()[2]);
ScalarType d;
if(nx>ny && nx>nz) { _flags |= NORMX; d = 1/plane.Direction()[0]; }
else if(ny>nz) { _flags |= NORMY; d = 1/plane.Direction()[1]; }
else { _flags |= NORMZ; d = 1/plane.Direction()[2]; }
// Scalatura spigoli
edge[0] *= d;
edge[1] *= d;
edge[2] *= d;
#endif
}
/*
Point face distance
trova il punto <p> sulla faccia piu' vicino a <q>, con possibilità di
rejection veloce su se la distanza trovata è maggiore di <rejdist>
Commenti del 12/11/02
Funziona solo se la faccia e di quelle di tipo E (con edge e piano per faccia gia' calcolati)
algoritmo:
1) si calcola la proiezione <p> di q sul piano della faccia
2) se la distanza punto piano e' > rejdist ritorna
3) si lavora sul piano migliore e si cerca di capire se il punto sta dentro il triangolo:
a) prodotto vettore tra edge triangolo (v[i+1]-v[i]) e (p-v[i])
b) se il risultato e' negativo (gira in senso orario) allora il punto
sta fuori da quella parte e si fa la distanza punto segmento.
c) se il risultato sempre positivo allora sta dentro il triangolo
4) e si restituisce la distanza punto /piano gia` calcolata
Note sulla robustezza:
il calcolo del prodotto vettore e` la cosa piu` delicata:
possibili fallimenti quando a^b ~= 0
1) doveva essere <= 0 e viene positivo (q era fuori o sulla linea dell'edge)
allora capita che si faccia la distanza punto piano anziche` la distanza punto seg
2) doveva essere > 0 e viene <=0 (q era dentro il triangolo)
*/
bool Dist( const Point3<ScalarType> & q, ScalarType & dist, Point3<ScalarType> & p )
{
#ifdef __VCGLIB_FACE_RT
//const ScalarType EPSILON = ScalarType( 0.000001);
const ScalarType EPSILON = 0.00000001;
ScalarType b,b0,b1,b2;
// Calcolo distanza punto piano
ScalarType d = Distance( plane, q );
if( d>dist || d<-dist ) // Risultato peggiore: niente di fatto
return false;
// Calcolo del punto sul piano
// NOTA: aggiunto un '-d' in fondo Paolo C.
Point3<ScalarType> t = plane.Direction();
t[0] *= -d;
t[1] *= -d;
t[2] *= -d;
p = q; p += t;
#define PP(i) (v[i]->cP())
#define E(i) (edge[i])
switch( _flags & (NORMX|NORMY|NORMZ) )
{
case NORMX:
b0 = E(1)[1]*(p[2] - PP(1)[2]) - E(1)[2]*(p[1] - PP(1)[1]);
if(b0<=0)
{
b0 = PSDist(q,V(1)->cP(),V(2)->cP(),p);
if(dist>b0) { dist = b0; return true; }
else return false;
}
b1 = E(2)[1]*(p[2] - PP(2)[2]) - E(2)[2]*(p[1] - PP(2)[1]);
if(b1<=0)
{
b1 = PSDist(q,V(2)->cP(),V(0)->cP(),p);
if(dist>b1) { dist = b1; return true; }
else return false;
}
b2 = E(0)[1]*(p[2] - PP(0)[2]) - E(0)[2]*(p[1] - PP(0)[1]);
if(b2<=0)
{
b2 = PSDist(q,V(0)->cP(),V(1)->cP(),p);
if(dist>b2) { dist = b2; return true; }
else return false;
}
// sono tutti e tre > 0 quindi dovrebbe essere dentro;
// per sicurezza se il piu' piccolo dei tre e' < epsilon (scalato rispetto all'area della faccia
// per renderlo dimension independent.) allora si usa ancora la distanza punto
// segmento che e' piu robusta della punto piano, e si fa dalla parte a cui siamo piu'
// vicini (come prodotto vettore)
// Nota: si potrebbe rendere un pochino piu' veloce sostituendo Area()
// con il prodotto vettore dei due edge in 2d lungo il piano migliore.
if( (b=min(b0,min(b1,b2))) < EPSILON*Area())
{
ScalarType bt;
if(b==b0) bt = PSDist(q,V(1)->cP(),V(2)->cP(),p);
else if(b==b1) bt = PSDist(q,V(2)->cP(),V(0)->cP(),p);
else if(b==b2) bt = PSDist(q,V(0)->cP(),V(1)->cP(),p);
//printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
if(dist>bt) { dist = bt; return true; }
else return false;
}
break;
case NORMY:
b0 = E(1)[2]*(p[0] - PP(1)[0]) - E(1)[0]*(p[2] - PP(1)[2]);
if(b0<=0)
{
b0 = PSDist(q,V(1)->cP(),V(2)->cP(),p);
if(dist>b0) { dist = b0; return true; }
else return false;
}
b1 = E(2)[2]*(p[0] - PP(2)[0]) - E(2)[0]*(p[2] - PP(2)[2]);
if(b1<=0)
{
b1 = PSDist(q,V(2)->cP(),V(0)->cP(),p);
if(dist>b1) { dist = b1; return true; }
else return false;
}
b2 = E(0)[2]*(p[0] - PP(0)[0]) - E(0)[0]*(p[2] - PP(0)[2]);
if(b2<=0)
{
b2 = PSDist(q,V(0)->cP(),V(1)->cP(),p);
if(dist>b2) { dist = b2; return true; }
else return false;
}
if( (b=min(b0,min(b1,b2))) < EPSILON*Area())
{
ScalarType bt;
if(b==b0) bt = PSDist(q,V(1)->cP(),V(2)->cP(),p);
else if(b==b1) bt = PSDist(q,V(2)->cP(),V(0)->cP(),p);
else if(b==b2) bt = PSDist(q,V(0)->cP(),V(1)->cP(),p);
//printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
if(dist>bt) { dist = bt; return true; }
else return false;
}
break;
case NORMZ:
b0 = E(1)[0]*(p[1] - PP(1)[1]) - E(1)[1]*(p[0] - PP(1)[0]);
if(b0<=0)
{
b0 = PSDist(q,V(1)->cP(),V(2)->cP(),p);
if(dist>b0) { dist = b0; return true; }
else return false;
}
b1 = E(2)[0]*(p[1] - PP(2)[1]) - E(2)[1]*(p[0] - PP(2)[0]);
if(b1<=0)
{
b1 = PSDist(q,V(2)->cP(),V(0)->cP(),p);
if(dist>b1) { dist = b1; return true; }
else return false;
}
b2 = E(0)[0]*(p[1] - PP(0)[1]) - E(0)[1]*(p[0] - PP(0)[0]);
if(b2<=0)
{
b2 = PSDist(q,V(0)->cP(),V(1)->cP(),p);
if(dist>b2) { dist = b2; return true; }
else return false;
}
if( (b=min(b0,min(b1,b2))) < EPSILON*Area())
{
ScalarType bt;
if(b==b0) bt = PSDist(q,V(1)->cP(),V(2)->cP(),p);
else if(b==b1) bt = PSDist(q,V(2)->cP(),V(0)->cP(),p);
else if(b==b2) bt = PSDist(q,V(0)->cP(),V(1)->cP(),p);
//printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
if(dist>bt) { dist = bt; return true; }
else return false;
}
break;
}
#undef E
#undef PP
dist = ScalarType(fabs(d));
//dist = Distance(p,q);
#endif
return true;
}
/// return the index [0..2] of a vertex in a face
inline int VertexIndex( const FVTYPE * w ) const
{