00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <CompuCell3D/Simulator.h>
00026 #include <CompuCell3D/Potts3D/Potts3D.h>
00027
00028
00029 #include <CompuCell3D/Potts3D/CellInventory.h>
00030 #include <CompuCell3D/Field3D/Field3D.h>
00031 #include <CompuCell3D/Field3D/WatchableField3D.h>
00032 #include <CompuCell3D/Boundary/BoundaryStrategy.h>
00033
00034 using namespace CompuCell3D;
00035
00036 #include <iostream>
00037 #include <cmath>
00038 using namespace std;
00039
00040 #define EXP_STL
00041 #include "NeighborTrackerPlugin.h"
00042
00043 NeighborTrackerPlugin::NeighborTrackerPlugin() :
00044 cellFieldG(0),
00045 checkSanity(false),
00046 checkFreq(1),
00047 periodicX(false),
00048 periodicY(false),
00049 periodicZ(false),
00050
00051 maxNeighborIndex(0),
00052 boundaryStrategy(0)
00053
00054 {
00055
00056 }
00057
00058 NeighborTrackerPlugin::~NeighborTrackerPlugin() {
00059 }
00060
00061 void NeighborTrackerPlugin::init(Simulator *_simulator, CC3DXMLElement *_xmlData) {
00062
00063
00064 if (_xmlData && _xmlData->findElement("CheckLatticeSanityFrequency")) {
00065 checkSanity=true;
00066 checkFreq=_xmlData->getFirstElement("CheckLatticeSanityFrequency")->getUInt();
00067 }
00068 cout<<"INITIALIZING CELL BOUNDARYTRACKER PLUGIN"<<endl;
00069 simulator=_simulator;
00070 Potts3D *potts = simulator->getPotts();
00071 cellFieldG = (WatchableField3D<CellG *> *)potts->getCellFieldG();
00072
00074 cellInventoryPtr=& potts->getCellInventory();
00075
00076
00077
00079 BasicClassAccessorBase * cellBTAPtr=&neighborTrackerAccessor;
00083 potts->getCellFactoryGroupPtr()->registerClass(cellBTAPtr);
00084
00085 potts->registerCellGChangeWatcher(this);
00086
00087 fieldDim=cellFieldG->getDim();
00088
00089 adjNeighbor = AdjacentNeighbor(fieldDim);
00090 if(potts->getBoundaryXName()=="Periodic"){
00091 adjNeighbor.setPeriodicX();
00092 periodicX=true;
00093 }
00094 if(potts->getBoundaryYName()=="Periodic"){
00095 adjNeighbor.setPeriodicY();
00096 periodicY=true;
00097 }
00098 if(potts->getBoundaryZName()=="Periodic"){
00099 adjNeighbor.setPeriodicZ();
00100 periodicZ=true;
00101 }
00102
00103
00104 boundaryStrategy=BoundaryStrategy::getInstance();
00105 maxNeighborIndex=boundaryStrategy->getMaxNeighborIndexFromNeighborOrder(1);
00106
00107
00108 }
00109
00111
00112 void NeighborTrackerPlugin::field3DChange(const Point3D &pt, CellG *newCell,
00113 CellG *oldCell) {
00114
00115
00116 if(newCell==oldCell){
00117 return;
00118 }
00119
00120 const Field3DIndex & field3DIndex=adjNeighbor.getField3DIndex();
00121
00122 long currentPtIndex=0;
00123 long adjNeighborIndex=0;
00124 long adjFace2FaceNeighborIndex=0;
00125
00126 CellG * currentCellPtr=0;
00127 CellG * adjCellPtr=0;
00128
00129 unsigned int token = 0;
00130 double distance;
00131 int oldDiff = 0;
00132 int newDiff = 0;
00133 Point3D n;
00134 Point3D ptAdj;
00135 CellG *nCell=0;
00136 Neighbor neighbor;
00137
00138 set<NeighborSurfaceData> * oldCellNeighborSurfaceDataSetPtr=0;
00139 set<NeighborSurfaceData> * newCellNeighborSurfaceDataSetPtr=0;
00140 pair<set<NeighborSurfaceData>::iterator,bool > set_NSD_itr_OK_Pair;
00141 set<NeighborSurfaceData>::iterator set_NSD_itr;
00142
00143
00144 if(newCell){
00145
00146 newCellNeighborSurfaceDataSetPtr = &neighborTrackerAccessor.get(newCell->extraAttribPtr)->cellNeighbors;
00147 }
00148
00149 if(oldCell){
00150
00151 oldCellNeighborSurfaceDataSetPtr = &neighborTrackerAccessor.get(oldCell->extraAttribPtr)->cellNeighbors;
00152 }
00153
00154
00155 currentPtIndex=field3DIndex.index(pt);
00156 currentCellPtr=cellFieldG->getByIndex(currentPtIndex);
00157
00158
00159 if(oldCell){
00160
00162 long temp_index;
00163 for(unsigned int nIdx=0 ; nIdx <= maxNeighborIndex ; ++nIdx ){
00164 neighbor=boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt),nIdx);
00165 if(!neighbor.distance){
00166
00167 continue;
00168 }
00169 adjCellPtr=cellFieldG->get(neighbor.pt);
00170
00171 if( adjCellPtr != oldCell ){
00172
00173
00174 set_NSD_itr = oldCellNeighborSurfaceDataSetPtr->find(NeighborSurfaceData(adjCellPtr));
00175 if( set_NSD_itr != oldCellNeighborSurfaceDataSetPtr->end() ){
00176 set_NSD_itr->decrementCommonSurfaceArea(*set_NSD_itr);
00177 if(set_NSD_itr->OKToRemove())
00178 oldCellNeighborSurfaceDataSetPtr->erase(set_NSD_itr);
00179
00180 }else{
00181 cerr<<"Could not find cell address in the boundary - set of cellNeighbors is corrupted. Exiting ..."<<endl;
00182 exit(0);
00183 }
00184
00185
00186 if(adjCellPtr){
00187 set<NeighborSurfaceData> &set_NSD_ref = neighborTrackerAccessor.get(adjCellPtr->extraAttribPtr)->cellNeighbors;
00188 set<NeighborSurfaceData>::iterator sitr;
00189 sitr=set_NSD_ref.find(oldCell);
00190 if(sitr!=set_NSD_ref.end()){
00191 sitr->decrementCommonSurfaceArea(*sitr);
00192 if(sitr->OKToRemove())
00193 set_NSD_ref.erase(sitr);
00194
00195 }
00196 }
00197 }
00198
00199
00200
00201 }
00202 }
00203
00204 if(newCell){
00205
00206
00208 for(unsigned int nIdx=0 ; nIdx <= maxNeighborIndex ; ++nIdx ){
00209 neighbor=boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt),nIdx);
00210 if(!neighbor.distance){
00211
00212 continue;
00213 }
00214 adjCellPtr=cellFieldG->get(neighbor.pt);
00215
00216
00217 if( adjCellPtr != newCell ){
00218
00219 set_NSD_itr_OK_Pair=newCellNeighborSurfaceDataSetPtr->insert(NeighborSurfaceData(adjCellPtr));
00220
00221
00222 set_NSD_itr=set_NSD_itr_OK_Pair.first;
00223 set_NSD_itr->incrementCommonSurfaceArea(*set_NSD_itr);
00224
00225 if(adjCellPtr){
00226 set<NeighborSurfaceData> &set_NSD_ref = neighborTrackerAccessor.get(adjCellPtr->extraAttribPtr)->cellNeighbors;
00227 pair<set<NeighborSurfaceData>::iterator,bool> sitr_OK_pair=set_NSD_ref.insert(newCell);
00228 set<NeighborSurfaceData>::iterator sitr=sitr_OK_pair.first;
00229 sitr->incrementCommonSurfaceArea(*sitr);
00230 }
00231
00232 }
00233
00234
00235
00236
00237 }
00238 }
00239
00240 token = 0;
00241 distance = 0;
00242
00243
00244
00245 if(!oldCell){
00246
00247
00249 for(unsigned int nIdx=0 ; nIdx <= maxNeighborIndex ; ++nIdx ){
00250 neighbor=boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt),nIdx);
00251 if(!neighbor.distance){
00252
00253 continue;
00254 }
00255 adjCellPtr=cellFieldG->get(neighbor.pt);
00256
00257
00258 if( adjCellPtr != oldCell && !(ptAdj == pt)){
00259
00260 if(adjCellPtr){
00261 set<NeighborSurfaceData> &set_NSD_ref = neighborTrackerAccessor.get(adjCellPtr->extraAttribPtr)->cellNeighbors;
00262 set<NeighborSurfaceData>::iterator sitr;
00263 sitr=set_NSD_ref.find(oldCell);
00264 if(sitr!=set_NSD_ref.end()){
00265 sitr->decrementCommonSurfaceArea(*sitr);
00266 if(sitr->OKToRemove()){
00267 set_NSD_ref.erase(sitr);
00268
00269
00270 }
00271
00272 }
00273 }
00274 }
00275
00276
00277
00278 }
00279
00280
00281 }
00282
00283 token = 0;
00284 distance = 0;
00285
00286
00287 if(!newCell){
00288
00289
00290 for(unsigned int nIdx=0 ; nIdx <= maxNeighborIndex ; ++nIdx ){
00291 neighbor=boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt),nIdx);
00292 if(!neighbor.distance){
00293
00294 continue;
00295 }
00296
00297 if(cellFieldG->isValid(neighbor.pt)){
00298
00299 adjCellPtr=cellFieldG->get(neighbor.pt);
00300
00301 if( adjCellPtr != newCell ){
00302
00303 if(adjCellPtr){
00304 set<NeighborSurfaceData> &set_NSD_ref = neighborTrackerAccessor.get(adjCellPtr->extraAttribPtr)->cellNeighbors;
00305 pair<set<NeighborSurfaceData>::iterator,bool> sitr_OK_pair=set_NSD_ref.insert(newCell);
00306 set<NeighborSurfaceData>::iterator sitr=sitr_OK_pair.first;
00307 sitr->incrementCommonSurfaceArea(*sitr);
00308 }
00309
00310 }
00311
00312 }
00313
00314
00315 }
00316
00317 }
00318
00319
00321
00322 if(checkSanity){
00323
00324 ++changeCounter;
00325
00326 if(!(changeCounter % checkFreq)){
00327
00328 cerr<<"ChangeCounter:"<<changeCounter<<endl;
00329 testLatticeSanityFull();
00330
00331 }
00332 }
00333
00334
00335 }
00336
00338
00339 double distance(double x1,double y1,double z1,double x2,double y2,double z2){
00340 return sqrt (
00341 (x1-x2)*(x1-x2)+
00342 (y1-y2)*(y1-y2)+
00343 (z1-z2)*(z1-z2)
00344 );
00345 }
00346
00347
00348 void NeighborTrackerPlugin::testLatticeSanityFull(){
00349
00350
00351
00352 Dim3D fieldDim=cellFieldG->getDim();
00353
00354 Point3D pt(0,0,0);
00355 Point3D ptAdj;
00356 Neighbor neighbor;
00357
00358
00359
00361
00362
00363 CellG * currentCellPtr;
00364 CellG * adjCellPtr;
00365
00366 map<CellG*,set<NeighborSurfaceData> > mapCellNeighborSurfaceData;
00367 map<CellG*,set<NeighborSurfaceData> >::iterator mitr;
00368
00370
00371 unsigned int token = 0;
00372 double distance;
00373
00374 set<NeighborSurfaceData> * set_NSD_ptr;
00375
00376 pair<set<NeighborSurfaceData>::iterator,bool > set_NSD_itr_OK_Pair;
00377 set<NeighborSurfaceData>::iterator set_NSD_itr;
00378
00379
00380 for(int z=0 ; z < fieldDim.z ; ++z)
00381 for(int y=0 ; y < fieldDim.y ; ++y)
00382 for(int x=0 ; x < fieldDim.x ; ++x){
00383 pt.x=x;
00384 pt.y=y;
00385 pt.z=z;
00386
00387 token=0;
00388 distance=0;
00389
00390
00391
00392 currentCellPtr=cellFieldG->get(pt);
00393 if(!currentCellPtr)
00394 continue;
00395 if(!isBoundaryPixel(pt))
00396 continue;
00397
00398 mitr=mapCellNeighborSurfaceData.find(currentCellPtr);
00399 if(mitr != mapCellNeighborSurfaceData.end() ){
00400 set_NSD_ptr = & (mitr->second) ;
00401 }else{
00402 mapCellNeighborSurfaceData.insert(make_pair(currentCellPtr,set<NeighborSurfaceData>()));
00403 mitr=mapCellNeighborSurfaceData.find(currentCellPtr);
00404 set_NSD_ptr = &(mitr->second);
00405 }
00406
00407 for(unsigned int nIdx=0 ; nIdx <= maxNeighborIndex ; ++nIdx ){
00408 neighbor=boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt),nIdx);
00409 if(!neighbor.distance){
00410
00411 continue;
00412 }
00413 adjCellPtr=cellFieldG->get(neighbor.pt);
00414
00415 if(adjCellPtr != currentCellPtr){
00416 set_NSD_itr_OK_Pair = set_NSD_ptr->insert(NeighborSurfaceData(adjCellPtr));
00417
00418
00419 set_NSD_itr=set_NSD_itr_OK_Pair.first;
00420 set_NSD_itr->incrementCommonSurfaceArea(*set_NSD_itr);
00421
00422 }
00423
00424 }
00425
00426 }
00427
00428
00429
00430 if(mapCellNeighborSurfaceData.size() != cellInventoryPtr->getCellInventorySize()){
00431 cerr<<"Number of cells in the mapCellNeighborSurfaceData = "<<mapCellNeighborSurfaceData.size()
00432 <<" is different than in cell inventory: "<< cellInventoryPtr->getCellInventorySize()<<endl;
00433 exit(0);
00434 }
00435 CellInventory::cellInventoryIterator cInvItr;
00436 CellG * cell;
00437 int counter=0;
00438
00439 for(cInvItr=cellInventoryPtr->cellInventoryBegin() ; cInvItr !=cellInventoryPtr->cellInventoryEnd() ;++cInvItr ){
00440
00441 cell=*cInvItr;
00442 mitr = mapCellNeighborSurfaceData.find(cell);
00443 if(mitr==mapCellNeighborSurfaceData.end()){
00444 cerr<<"Cell "<<cell<<" does not appear in the just initialized mapCellNeighborSurfaceData"<<endl;
00445 exit(0);
00446
00447 }
00448 set_NSD_ptr=&(mitr->second);
00449
00450 if(! (*set_NSD_ptr == neighborTrackerAccessor.get(cell->extraAttribPtr)->cellNeighbors)){
00451 cerr<<"Have checked "<<counter<<" cells"<<endl;
00452 cerr<<"set of NeighborSurfaceData do not match for cell: "<<cell<<endl;
00453 exit(0);
00454 }
00455
00456 ++counter;
00457 }
00458
00459 cerr<<"FULL TEST: LATTICE IS SANE!!!!!"<<endl;
00460
00461 }
00462
00464 bool NeighborTrackerPlugin::isBoundaryPixel(Point3D pt){
00465
00466
00467 CellG * currentCellPtr=cellFieldG->get(pt);;
00468
00469 CellG * nCell;
00470 unsigned int token = 0;
00471 double distance = 0;
00472 Point3D n;
00473
00474
00475
00476 while (true) {
00477 n = cellFieldG->getNeighbor(pt, token, distance, false);
00478 if (distance > 1) break;
00479 nCell = cellFieldG->get(n);
00480 if(nCell != currentCellPtr)
00481 return true;
00482 }
00483
00484 return false;
00485
00486 }
00487
00488
00489 std::string NeighborTrackerPlugin::toString(){
00490 return "NeighborTracker";
00491 }
