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 |