BOOKS i'm reading
 |
/*
* Module ID: chains.cpp
* Titre : Definition des classes de chaines.
* Utilite : Structure de donnee pour la location de point
* comme decrit dans 'Computational Geometry ' de
* Franco P. Preparata & Michael Ian Shamos (P. 48)
*
* Note : Les chaines monotones utilisees sur un diagramme de Voronoi
* permettent de trouver le plus proche voisin d'un point x en
* O(nlogn) comparativement a O(n^2) avec une methode brutale.
*
* Auteur : Olivier Langlois <olivier@olivierlanglois.net>
* Date : 24 Mars 1998
*
* Revision :
*
* 001 27-Mar-1998 - Olivier Langlois
* Les fstreams ne fonctionnent pas avec le compilateur GNU 2.7.2.1
* lorsque compile avec l'option -frtti.
*/
#ifdef DEBUG
#ifndef __GNUC__
#include "diagwin.h"
#include <typeinfo.h>
#else
#include <std/typeinfo.h>
#include "diagdos.h"
#endif
#endif
#include "collect/colvar.h"
#include "collect/dict.h"
#include "../tools/collect/dynastck.cpp"
#include "../tools/collect/cirqueue.cpp"
#include "../tools/collect/bintree.cpp"
#include "../tools/collect/rbtree.cpp"
#include "../tools/collect/pq.cpp"
#include "geo/point.h"
#include "geo/dcel.h"
#include "geo/chains.h"
#ifdef __GNUG__
#include "umath.h"
#else
#include <math.h>
#endif
#include "mathc.h"
DEFINE_COLLECTABLE( ChainMesh, __CHAINMESH )
DEFINE_COLLECTABLE( Chain , __CHAIN )
DEFINE_COLLECTABLE( ChainsSet, __CHAINSSET )
DCEL *Chain::edgesPtr = NULL;
#ifdef __GNUG__
void *restorePtr;
void ChainMesh::saveGuts( FILE *os )
#else
void ChainMesh::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
fwrite( &myAtom, sizeof(ClassID), 1, os );
edgeIdx.saveGuts(os);
fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
os.write((char *)&myAtom, sizeof(ClassID));
edgeIdx.saveGuts(os);
os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void ChainMesh::restoreGuts( FILE *is )
#else
void ChainMesh::restoreGuts( ifstream &is )
#endif
{
ClassID id;
// Verifie la validite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw int(0);
edgeIdx.restoreGuts(is);
// Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw int(0);
}
int ChainMesh::operator==(const ChainMesh &b) const
{
vertex *v1 = edgesPtr->getEdge(&edgeIdx)->getV1();
vertex *v2 = edgesPtr->getEdge(&edgeIdx)->getV2();
vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();
vertex *bv2 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();
return ((*v1 == *bv1 && *v2 == *bv2) ||
(*v1 == *bv2 && *v2 == *bv1) );
}
int ChainMesh::operator<(const ChainMesh &b) const
{
vertex *v1 = edgesPtr->getEdge(&edgeIdx)->getV1();
vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();
if( v1->isEqual(bv1) == BOOL_TRUE )
bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();
return( point::cmp_yx(*v1,*bv1)<0);
}
int ChainMesh::operator>(const ChainMesh &b) const
{
vertex *v1 = edgesPtr->getEdge(&edgeIdx)->getV1();
vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();
if( v1->isEqual(bv1) == BOOL_TRUE )
bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();
return( point::cmp_yx(*v1,*bv1)>0);
}
/******************************************************************************
*
* Nom : isInMesh
*
* Utilite : Verifie si le point est a la meme hauteur (ycoord) que le
* maillon.
*
* Reference : chains.h
*
****************************************************************************/
int ChainMesh::isInMesh( const point &p )
{
vertex *v1 = edgesPtr->getEdge(&edgeIdx)->getV1();
vertex *v2 = edgesPtr->getEdge(&edgeIdx)->getV2();
double resultV1 = point::cmp_yx(p, *v1);
double resultV2 = point::cmp_yx(p, *v2);
if( resultV1 != resultV2 ) return 0;
else if( resultV1 < 0 && resultV2 < 0 ) return -1;
else return 1;
}
/******************************************************************************
*
* Nom : findMesh
*
* Utilite : Trouve le maillon qui est a la meme hauteur que p.
*
* Reference : chains.h
*
****************************************************************************/
CollectableULong *Chain::findMesh( const point &p )
{
TreeNode< RBData<ChainMesh> > *n = Root;
while(n != Sentinel)
{
int result = n->Data->RealData->isInMesh(p);
if (result < 0)
n = n->Less;
else if ( result > 0 )
n = n->More;
else
break;
}
if( n != Sentinel ) return n->Data->RealData->getEdgeIdx();
else
{
return NULL;
}
}
#ifdef __GNUG__
void Chain::saveGuts( FILE *os )
#else
void Chain::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
fwrite( &myAtom, sizeof(ClassID), 1, os );
fwrite( &chainNum, sizeof(unsigned long), 1, os );
RedBlackTree<ChainMesh>::saveGuts(os);
fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
os.write((char *)&myAtom, sizeof(ClassID));
os.write((char *)&chainNum, sizeof(unsigned long));
RedBlackTree<ChainMesh>::saveGuts(os);
os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void Chain::restoreGuts( FILE *is )
#else
void Chain::restoreGuts( ifstream &is )
#endif
{
ClassID id;
// Verifie la validite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw int(0);
#ifdef __GNUG__
fread( &chainNum, sizeof(unsigned long), 1, is );
#else
is.read((char *)&chainNum, sizeof(unsigned long));
#endif
RedBlackTree<ChainMesh>::restoreGuts(is);
// Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw int(0);
}
void Chain::saveChainMesh( const ChainMesh &item, void *param )
{
#ifdef __GNUG__
FILE *f = (FILE *)param;
#else
ofstream &f = *(ofstream *)(ostream *)param;
#endif
const_cast< ChainMesh & >(item).saveGuts( f );
}
void Chain::restoreChainMesh( const ChainMesh &item, void *param )
{
#ifdef __GNUG__
FILE *f = (FILE *)param;
#else
ifstream &f = *(ifstream *)(istream *)param;
#endif
const_cast<ChainMesh &>(item).restoreGuts( f );
const_cast<ChainMesh &>(item).edgesPtr = edgesPtr;
}
ChainsSet::ChainsSet()
: BinaryTree<Chain>(&ChainsSet::saveChain,&ChainsSet::restoreChain),
edgesPtr(NULL) {}
ChainsSet::ChainsSet( DCEL *list )
: BinaryTree<Chain>(&ChainsSet::saveChain,&ChainsSet::restoreChain),
edgesPtr(list) {init();}
// Classe utilisee par ChainsSet::init
class weight : public Collectable
{
public:
weight() : w(1), min(0), max(0) {}
size_t w, min, max;
};
void ChainsSet::init()
{
Dictionary weightTable;
IDList sortedVertices;
CollectableULong ymin = edgesPtr->Vertices().getyminIdx();
CollectableULong ymax = edgesPtr->Vertices().getymaxIdx();
// Initialisation de weightTable
DictIterator it(edgesPtr->Edges());
while( ++it == BOOL_TRUE )
weightTable.Insert( it.getKey()->copy(), new weight );
DictIterator VerticesIt(edgesPtr->Vertices());
while( ++VerticesIt == BOOL_TRUE )
{
// Insere toutes les vertices dans la liste sauf ymin & ymax
if( VerticesIt.getKey()->isEqual(&ymin) == BOOL_FALSE &&
VerticesIt.getKey()->isEqual(&ymax) == BOOL_FALSE )
// Ici, on suppose que les vertices ne se trouve pas dans une autre liste
sortedVertices.insert((point *)VerticesIt.getValue());
}
sortedVertices.sort(sortSites);
unsigned long Win, Wout;
/*
* d1 : Arete sortante la plus a gauche
*/
edge *d1;
segment sd1;
// Premiere passe
segment curSeg;
IDListIterator svit(sortedVertices);
vertex *curVertex;
edge *curEdge;
CollectableULong *idx;
D( while( (curVertex = (vertex *)++svit) ) cout << *curVertex << endl; )
D( svit.reset(); )
while( (curVertex = (vertex *)++svit) )
{
IDListIterator eit(*curVertex->getEdgeIdxList());
Wout = 0;
Win = 0;
d1 = NULL;
while ( (idx = (CollectableULong *)++eit) )
{
curEdge = edgesPtr->getEdge(idx);
vertex *secondVertex = curEdge->getV1();
if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
curSeg = segment(*curVertex,*curEdge->getV2());
else
curSeg = segment(*curVertex,*secondVertex);
D( cout << "Angle :" << curSeg.angle() << endl; )
if( curSeg.angle() >= -0.00001 && curSeg.angle() < M_PI-0.00001 )
{
// Arete sortante
if( !d1 || sd1.angle() < curSeg.angle() )
{
d1 = curEdge;
sd1 = curSeg;
}
Wout += ((weight *)
weightTable.getValue(idx))->w;
}
else
{
// Arete rentrante
Win += ((weight *)
weightTable.getValue(idx))->w;
}
} // WHILE(edge)
// A la place d'une exception, une arete de correction pourrait etre
// cree.
if( !Win || ! Wout ) throw GeoEx(GX_REGPSLG);
if( Win > Wout )
{
((weight *)
weightTable.getValue(d1->getKey()))->w = Win-Wout+1;
}
} // WHILE(vertices)
// A la fin de la premiere passe, on peut connaitre le nombre de chaines.
PQ<edgeInfo> edgePQ;
edgeInfo Sentinel;
// Une clef valide ne peut pas etre negative.
Sentinel.Key = -100;
#ifdef __BORLANDC__
PQ<edgeInfo>::setSentinel(&Sentinel);
#else
edgePQ.setSentinel(&Sentinel);
#endif
edgeInfo *curEdgeInfo;
size_t numChains = 0;
int curChain = 0;
curVertex = edgesPtr->getVertex(&ymax);
IDListIterator ymaxEdgeit
(*curVertex->getEdgeIdxList());
while( (idx = (CollectableULong *)++ymaxEdgeit) )
{
numChains += ((weight *)weightTable.getValue(idx))->w;
curEdge = edgesPtr->getEdge(idx);
vertex *secondVertex = curEdge->getV1();
if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
curSeg = segment(*curVertex,*curEdge->getV2());
else
curSeg = segment(*curVertex,*secondVertex);
double Angle = curSeg.angle();
curEdgeInfo = new edgeInfo;
// Si l'angle est PI, il est change pour -PI
// Cela est necessaire pour l'ordre des elements de edgePQ
curEdgeInfo->Key = ((Angle > 0)?-Angle:Angle);
curEdgeInfo->e = curEdge;
edgePQ.insert(curEdgeInfo);
}
D( cout << "numChains :" << numChains << endl; )
for( unsigned long i = 0; i < numChains; i++ )
Insert(new Chain(i));
balance();
while( (curEdgeInfo = edgePQ.removeFirst()) != &Sentinel )
{
weight *w = (weight *)weightTable.getValue(curEdgeInfo->e->getKey());
w->min = curChain;
w->max = curChain + w->w - 1;
curChain += w->w;
delete curEdgeInfo;
}
// Deuxieme passe
while( (curVertex = (vertex *)--svit) )
{
IDListIterator eit(*curVertex->getEdgeIdxList());
Wout = 0;
Win = 0;
curChain = numChains;
while ( (idx = (CollectableULong *)++eit) )
{
curEdge = edgesPtr->getEdge(idx);
vertex *secondVertex = curEdge->getV1();
if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
curSeg = segment(*curVertex,*curEdge->getV2());
else
curSeg = segment(*curVertex,*secondVertex);
D( cout << "Angle :" << curSeg.angle() << endl; )
// La constante M_PI est plus precise que le nombre retourne par
// segment::angle() lorsque la reponse est PI c'est pourquoi le test
// == n'est pas utilise
if( curSeg.angle() >= -0.00001 && curSeg.angle() < M_PI-0.00001 )
{
// Arete sortante
weight *w = (weight *)
weightTable.getValue(idx);
Wout += w->w;
insert( new ChainMesh(curEdge->getKey(), edgesPtr), w->min, w->max );
if( w->min < curChain )
curChain = w->min;
}
else
{
// Arete rentrante
double Angle = curSeg.angle();
curEdgeInfo = new edgeInfo;
// Si l'angle est PI, il est change pour -PI
// Cela est necessaire pour l'ordre des elements de edgePQ
curEdgeInfo->Key = ((Angle > 0)?-Angle:Angle);
curEdgeInfo->e = curEdge;
edgePQ.insert(curEdgeInfo);
Win += ((weight *)
weightTable.getValue(idx))->w;
}
} // WHILE(edge)
if( Win < Wout )
{
weight *w = (weight *)weightTable.getValue(edgePQ.first()->e->getKey());
D( cout << "Win " << Win << " Wout " << Wout << " Wd2 " << w->w << endl; )
w->w = Wout-Win+w->w;
}
// Assigne les min/max aretes entrante
while( (curEdgeInfo = edgePQ.removeFirst()) != &Sentinel )
{
#ifdef __WIN32__
D( if( curChain < 0 || curChain > numChains) DiagOutWin("Erreur").DisplayMsg("curChain Corrompu", DIAG_WARNING); )
#endif
D( cout << "Segment :" << curEdgeInfo->e->getSegment() << endl; )
weight *w = (weight *)weightTable.getValue(curEdgeInfo->e->getKey());
w->min = curChain;
w->max = curChain + w->w - 1;
curChain += w->w;
delete curEdgeInfo;
}
} // WHILE(vertices)
// Insere les aretes connectees a ymin.
curVertex = edgesPtr->getVertex(&ymin);
IDListIterator eit(*curVertex->getEdgeIdxList());
while ( (idx = (CollectableULong *)++eit) )
{
weight *w = (weight *)
weightTable.getValue(idx);
insert( new ChainMesh(idx, edgesPtr), w->min, w->max );
}
// Vide la liste de ses vertices pour eviter qu'ils ne soient effaces
sortedVertices.clearWithoutDestroy();
}
/******************************************************************************
*
* Nom : insert
*
* Utilite : Insere une maille.
*
* Reference : chains.h
*
****************************************************************************/
void ChainsSet::insert( ChainMesh *m, unsigned long min, unsigned long max)
{
TreeNode<Chain> *n = Root;
D( cout << "min :" << min << " max :" << max << endl; )
while(n != Sentinel &&
(n->Data->getChainNum() < min || n->Data->getChainNum() > max) )
{
D( cout << "chainNum :" << n->Data->getChainNum() << endl; )
if (max < n->Data->getChainNum())
n = n->Less;
else
n = n->More;
}
if( n != Sentinel ) n->Data->Insert(m);
else
{
D( cout << "Chaine non trouvee" << endl; )
delete m;
}
}
/******************************************************************************
*
* Nom : location
*
* Utilite : Trouve la location de p.
*
* Reference : chains.h
*
****************************************************************************/
face *ChainsSet::location( const point &p )
{
TreeNode<Chain> *n = Root;
edge *e = NULL;
vertex *v1;
vertex *v2;
Boolean atRight;
face *loc;
while(1)
{
CollectableULong *result = n->Data->findMesh(p);
if( result )
{
e = edgesPtr->getEdge(result);
v1 = e->getV1();
v2 = e->getV2();
if( point::cmp_yx( *v1, *v2 ) > 0 )
{
vertex *tmp = v1;
v1 = v2;
v2 = tmp;
}
if( right_turn( *v1, *v2, p ) == BOOL_TRUE )
atRight = BOOL_TRUE;
// Le point est contenu dans le segment, Le point n'est dans aucune
// face mais les 2 faces les plus proches sont e.F1 et e.F2 et puisque
// dans la methode DCEL::setBox la face outside est toujour F1 et dans
// la majorite des applications elle moins utile que les autres faces,
// F2 est retournee dans cette situation.
else if( segment( *v1, *v2 ).contains(p) == BOOL_TRUE )
return e->getF2();
else
atRight = BOOL_FALSE;
if ( atRight == BOOL_FALSE && n->Less != Sentinel )
n = n->Less;
else if ( atRight == BOOL_TRUE && n->More != Sentinel )
n = n->More;
else
{
face *f1 = e->getF1();
face *f2 = e->getF2();
if( f1->getLocation().distance(p) < f2->getLocation().distance(p) )
loc = f1;
else
loc = f2;
break;
}
} // IF(result)
// Sinon continue a chercher l'arbre primaire
else if( e )
{
if ( atRight == BOOL_FALSE && n->Less != Sentinel )
n = n->Less;
else if ( atRight == BOOL_TRUE && n->More != Sentinel )
n = n->More;
else
{
face *f1 = e->getF1();
face *f2 = e->getF2();
if( f1->getLocation().distance(p) < f2->getLocation().distance(p) )
loc = f1;
else
loc = f2;
break;
}
}
// Pour eviter un loop infini mais cela ne devrait pas arriver.
else return NULL;
}
return loc;
}
#ifdef __GNUG__
void ChainsSet::saveGuts( FILE *os )
#else
void ChainsSet::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
fwrite( &myAtom, sizeof(ClassID), 1, os );
edgesPtr->saveGuts(os);
BinaryTree<Chain>::saveGuts(os);
fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
os.write((char *)&myAtom, sizeof(ClassID));
edgesPtr->saveGuts(os);
BinaryTree<Chain>::saveGuts(os);
os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void ChainsSet::restoreGuts( FILE *is )
#else
void ChainsSet::restoreGuts( ifstream &is )
#endif
{
ClassID id;
// Verifie la validite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw TreeEx(BTX_FILE);
delete edgesPtr;
edgesPtr = new DCEL;
edgesPtr->restoreGuts(is);
Chain::edgesPtr = edgesPtr;
BinaryTree<Chain>::restoreGuts(is);
// Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
fread( &id, sizeof(ClassID), 1, is );
#else
is.read( (char *)&id, sizeof(ClassID) );
#endif
if( id != myAtom ) throw TreeEx(BTX_FILE);
}
void ChainsSet::saveChain( const Chain &item, void *param )
{
#ifdef __GNUG__
FILE *f = (FILE *)param;
#else
ofstream &f = *(ofstream *)(ostream *)param;
#endif
const_cast< Chain & >(item).saveGuts( f );
}
void ChainsSet::restoreChain( const Chain &item, void *param )
{
#ifdef __GNUG__
FILE *f = (FILE *)param;
#else
ifstream &f = *(ifstream *)(istream *)param;
#endif
const_cast< Chain & >(item).restoreGuts( f );
}
#ifdef __GNUG__
template class PQ<edgeInfo>;
template class BinaryTree<Chain>;
template class BinaryTree<RBData<ChainMesh> >;
template class RBData<ChainMesh>;
template class RedBlackTree<ChainMesh>;
template class Queue<TreeNode<Chain> *>;
template class cirQueue<TreeNode<Chain> *>;
template class Queue<TreeNode<RBData<ChainMesh> > *>;
template class cirQueue<TreeNode<RBData<ChainMesh> > *>;
template class TreeNode< Chain >;
template class TreeNode< RBData<ChainMesh> >;
typedef TreeNode<Chain>* nodeChainPtr;
typedef TreeNode<RBData<ChainMesh> >* nodeChainMeshPtr;
template void saveNodeToFile( const nodeChainPtr &,void * );
template void saveNodeToFile( const nodeChainMeshPtr &,void * );
template void saveBTNodeToFile(
const Chain &item,
void *param
);
template void saveBTNodeToFile(
const RBData<ChainMesh> &item,
void *param
);
template void saveRBTDataToFile( const ChainMesh &, void * );
template void saveRBTNodeToFile( const RBData<ChainMesh> &item, void *param );
template void restoreRBTNodeFromFile( const RBData<ChainMesh> &item,
void *param );
template void ProcessData( const RBData<ChainMesh> &item, void *param );
#endif
|