COME - Biomodel

00001 /***************************************************************************
00002  *   Copyright (C) 2002 by Anderson Maciel                                 *
00003  *   andi.maciel@gmail.com                                                 *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  **************************************************************************/
00020 
00022 //
00023 //  PROJECT.....: CO-ME
00024 //  RESPONSIBLE.: 
00025 //
00026 //  AUTHOR......: Anderson Maciel
00027 //  DATE........: June/05/2002
00028 //  DESCRIPTION.: Class definition.
00029 //
00031 
00032 #include        <general/comescenario.h>
00033 #include        <physics/comeforce.h>
00034 #include        <physics/comecollide.h>
00035 #include        <bio/comepatient.h>
00036 #include        <bio/comemolecule.h>
00037 #include        <bio/comemoleculescartilage.h>
00038 #include        <bio/comemoleculesbone.h>
00039 #include        <bio/comeclamp.h>
00040 #include        <general/comesimulator.h>
00041 #include        <general/comexml.h>
00042 #include        <algebra/comevector3d.h>
00043 #include        <kinematics/comepolyaxialjoint.h>
00044 #include        <kinematics/comemovement.h>
00045 #include        <math.h>
00046 
00047 #ifndef _WIN32
00048 #include        <sys/time.h>
00049 #else 
00050 #include        <time.h>
00051 #endif
00052 
00053 #include        <zthread/Mutex.h>
00054 #define _MECHATESTER_APP
00058 using namespace ZThread;
00059 
00060 static double accu;
00061 
00062 COME_Scenario::COME_Scenario(){
00063 
00064         availableMaterials = new vector<COME_Material*>;
00065         movement = NULL;
00066         gravity = COME_Vector3D( 0.0, 0.0, 0.0 );
00067         collisions = new COME_Collide( this );
00068         modified = false;
00069         patientList.push_back( new COME_Patient() );
00070         patientList.front()->setParent( this );
00071         mutex = new Mutex();
00072         accu = 0.0;
00073 }
00074 
00075 COME_Scenario::~COME_Scenario(){
00076 
00077 }
00078 
00079 
00083 void
00084 COME_Scenario::setPatientList( list<COME_Patient*> &listN ){
00085 
00086         patientList = listN;
00087 }
00088 
00089 void
00090 COME_Scenario::setGravity( const COME_Vector3D &gravityN ){
00091         
00092         gravity = gravityN;
00093 }
00094 
00098 COME_Patient*
00099 COME_Scenario::getPatient( int index ) const {
00100 
00101         list<COME_Patient*>::const_iterator iter;
00102         int i = 0;
00103         for (iter = patientList.begin(); !(iter == patientList.end()); iter++){
00104                 if( i == index ){
00105                         return (*iter);
00106                 }
00107                 i++;
00108         }
00109         return NULL;
00110 }
00111 
00112 list<COME_Patient*>*
00113 COME_Scenario::getPatientList() {
00114 
00115         return &patientList;
00116 }
00117 
00118 const COME_Vector3D
00119 COME_Scenario::getGravity() const {
00120 
00121         return gravity;
00122 }
00123 
00124 COME_Collide*
00125 COME_Scenario::getCollisionDetector() const {
00126 
00127         return collisions;
00128 }
00129 
00130 vector<COME_Material*>*
00131 COME_Scenario::getAvailableMaterialsList() const {
00132 
00133         return availableMaterials;
00134 }
00135 
00142 bool
00143 COME_Scenario::simulate( double timestep, double simClock, bool direction ){
00144 
00145         // Variables to save axes into files
00146         //double currFlex = 0.0, currAdduct = 0.0, currTwist = 0.0;
00147 
00148         //lock();
00149         #ifndef _WIN32
00150         struct timeval *Tps, *Tpf;
00151         Tps = (struct timeval*) malloc(sizeof(struct timeval));
00152         Tpf = (struct timeval*) malloc(sizeof(struct timeval));
00153         #endif
00154 
00155         // update kinematical matrices
00156 
00157         int indPat = patientList.size() - 1;
00158         
00159         if( movement ){
00160 
00161                 #ifndef _WIN32
00162                 gettimeofday (Tps, NULL );
00163                 #endif
00164                 
00165                 movement->setTime_Err( timestep / 2.0 );
00166                 
00167                 for( int i = 0; i < movement->getQtdMotion(); i++ ){
00168         
00169                         COME_Joint *jointAux = getPatient(indPat)->getJointByName( movement->getTimeline()[i].getJointName() );
00170         
00171                         if( ( movement->getTimeline()[i].getTimeIni() > simClock - movement->getTime_Err() ) && ( movement->getTimeline()[i].getTimeIni() <= simClock + movement->getTime_Err() ) ){
00172                                 //calculate deltaS
00173                                 double gds = 0;
00174                                 switch( movement->getTimeline()[i].getMotionType() ){
00175                                         case FLEX       :       gds = movement->getTimeline()[i].getParameter() - (((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX])->getCurrent();
00176                                                                 break;
00177                                         case ADDUCT     :       gds = movement->getTimeline()[i].getParameter() - (((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT])->getCurrent();
00178                                                                 break;
00179                                         case TWIST      :       gds = movement->getTimeline()[i].getParameter() - (((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST])->getCurrent();
00180                                                                 break;
00181                                 }
00182                                 double gdt = ( movement->getTimeline()[i].getTimeFin() - movement->getTimeline()[i].getTimeIni() ) / timestep;
00183                                 movement->getTimeline()[i].setDeltaS( gds / gdt );
00184                         }
00185 
00186                         // move to new position 
00187                         double newPos = 0.0;
00188                         switch( movement->getTimeline()[i].getMotionType() ){
00189                                 case FLEX   : newPos = (((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX])->getCurrent() + movement->getTimeline()[i].getDeltaS();
00190                                         ((COME_PolyaxialJoint*)jointAux )->setFlexionTo( newPos );
00191                                         //printf( "FLEX:%f\n", newPos );
00192                                         //currFlex = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMin() ) * newPos + ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMin();
00193                                         break;
00194                                 case ADDUCT : newPos = (((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT])->getCurrent() + movement->getTimeline()[i].getDeltaS();
00195                                         ((COME_PolyaxialJoint*)jointAux )->setAdductionTo( newPos );
00196                                         //printf( "ADDUCT:%f\n", newPos );
00197                                         //currAdduct = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMin() ) * newPos + ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMin();
00198                                         break;
00199                                 case TWIST  : newPos = (((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST])->getCurrent() + movement->getTimeline()[i].getDeltaS();
00200                                         ((COME_PolyaxialJoint*)jointAux )->setTwistTo( newPos );
00201                                         //printf( "TWIST:%f\n", newPos );
00202                                         //currTwist = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMin() ) * newPos + ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMin();
00203                                         break;
00204                         }
00205 
00206 
00207                         if( ( movement->getTimeline()[i].getTimeFin() > simClock - movement->getTime_Err() ) && ( movement->getTimeline()[i].getTimeFin() <= simClock + movement->getTime_Err() ) ){
00208                                 //reset deltaS
00209                                 movement->getTimeline()[i].setDeltaS( 0.0000000000 );
00210                         }
00211                 }
00212         }
00213 
00214         // Apply changes. Should be replaced by a Patient->Update to be executed only if there are changes
00215         getPatient(indPat)->getRootJoint()->vpMakeGims( COME_Matrix() );
00216         #ifndef _WIN32
00217         gettimeofday (Tpf, NULL); 
00218         ((COME_Simulator*)parent)->difKinematics += (Tpf->tv_sec-Tps->tv_sec)*1000000 + Tpf->tv_usec-Tps->tv_usec;
00219         gettimeofday (Tps, NULL );
00220         #endif
00221         
00223         if( COME::flagCollisionTreatment == NEIGHBORS ){
00224         
00225                 collisions->updateProximityStructure();
00226                 collisions->geometricResponse();
00227         } else if( COME::flagCollisionTreatment == SPHERICAL_SLIDING ){
00228         
00229                 collisions->sphericalSlidingResponse();
00230         }
00231         #ifndef _WIN32
00232         gettimeofday (Tpf, NULL); 
00233         ((COME_Simulator*)parent)->difCollisions += (Tpf->tv_sec-Tps->tv_sec)*1000000 + Tpf->tv_usec-Tps->tv_usec;
00234         #endif
00235 
00236         
00237         list<COME_Patient*>::iterator iter;
00238         for (iter = patientList.begin(); !(iter == patientList.end()); iter++){
00239                 #ifndef _WIN32
00240                 gettimeofday (Tps, NULL );
00241                 #endif
00242                 if( direction == FUTURE ){
00243                         if( !(*iter)->update( timestep, simClock ) ){
00244                                 return false;
00245                         }
00246                 } else {
00247                         if( !(*iter)->update( -timestep, simClock )){
00248                                 return false;
00249                         }
00250                 }
00251                 
00252                 #ifndef _WIN32
00253                 gettimeofday (Tpf, NULL); 
00254                 ((COME_Simulator*)parent)->difDeformation += (Tpf->tv_sec-Tps->tv_sec)*1000000 + Tpf->tv_usec-Tps->tv_usec;
00255                 #endif
00256 
00257                 #ifdef _MECHATESTER_APP
00258                 // Update mobile clamp position during force test
00259                 COME_Vector3D dispClamp;
00260                 double totalVectors = 0.0;
00261                 list<COME_BioStructure*>::const_iterator iterOrgans;
00262                 for( iterOrgans = ((list<COME_BioStructure *>*)getPatient(indPat)->getPtOrganList())->begin(), iterOrgans++, iterOrgans++; iterOrgans != ((list<COME_BioStructure *>*)getPatient(indPat)->getPtOrganList())->end(); iterOrgans++  ){
00263                 
00264                         list<COME_Molecule*>* molecules = ((COME_MoleculesCartilage*)(*iterOrgans))->getTissue()->getShape();
00265                         list<COME_Molecule*>::iterator iterMolecules;
00266                         for ( iterMolecules = molecules->begin(); iterMolecules != molecules->end(); iterMolecules++){
00267                         
00268                                 if( (*iterMolecules)->intersectClampMobile() ){
00269                                         dispClamp += (*iterMolecules)->getPosition() - (*iterMolecules)->getPositionMinusOne( );
00270                                         totalVectors+=1.0;
00271                                 } 
00272                         }
00273                 }
00274                 if ( totalVectors ){
00275                         getPatient(indPat)->getOrgan("clamp1")->translate( dispClamp / totalVectors);
00276                         getPatient(indPat)->getOrgan("clamp1")->updateSkin();
00277                 }
00278                 #endif // MECHATESTER
00279         }
00280         if( COME::flagCollisionTreatment == NEIGHBORS ){
00281                 #ifndef _WIN32
00282                 gettimeofday (Tps, NULL );
00283                 #endif
00284                 collisions->resetProximitiesDisps();
00285                 #ifndef _WIN32
00286                 gettimeofday (Tpf, NULL); 
00287                 ((COME_Simulator*)parent)->difCollisions += (Tpf->tv_sec-Tps->tv_sec)*1000000 + Tpf->tv_usec-Tps->tv_usec;
00288                 #endif
00289         }
00290         
00291         
00292         // Detect collisions between the objects of the scene
00293         if( COME::flagCollisionTreatment == SPHERES ){
00294         
00295                 collisions->treatCollisions();
00296 
00297         } else if( ( COME::flagCollisionTreatment == MESH ) || (COME::flagCollisionTreatment == HYBRID )){
00298                 
00299                 collisions->detectContacts();
00300 
00301                 vector<COME_Collision> c = collisions->getCollisions();
00302                 for (int k=0; k < c.size(); k++ ){
00303                         
00304                         COME_BioStructure *organ1 = c[k].getOrgan1();
00305                         COME_BioStructure *organ2 = c[k].getOrgan2();
00306                         
00307                         COME_Force f1 = c[k].getForceInFace( organ1 );
00308                         COME_Force f2 = c[k].getForceInFace( organ2 );
00309 
00310                         int face1 = c[k].getFace( organ1 );
00311                         int face2 = c[k].getFace( organ2 );
00312 
00313                         organ1->getSurface()->getAFace( face1 ).getNearestMolecule()->addCollisionForce( f1 );
00314                         organ2->getSurface()->getAFace( face2 ).getNearestMolecule()->addCollisionForce( f2 );
00315                         
00316                 }
00317         } else if( COME::flagCollisionTreatment == DISPLACEMENT ){
00318         
00319                 if( collisions->detectContactsDisplacement() ){
00320                 
00321                         list<COME_BioStructure*>::const_iterator iterOrgans;
00322                         COME_Patient *activePatient = getPatient(getPatientList()->size()-1);
00323                         for( iterOrgans = activePatient->getPtOrganList()->begin(); iterOrgans != activePatient->getPtOrganList()->end(); iterOrgans++ ){
00324                 
00325                                 if( (*iterOrgans)->isFixed() ){
00326                                         ((COME_MoleculesBone*)(*iterOrgans))->respondCollision();
00327                                 } else {
00328                                         ((COME_MoleculesCartilage*)(*iterOrgans))->respondCollision();
00329                                 }
00330                         }
00331                 }
00332         
00333         } else if( COME::flagCollisionTreatment == NEIGHBORS ){
00334         
00335                 // DO NOTHING: Code is integrated in updateEuler
00336         
00337         } else if( COME::flagCollisionTreatment == SPHERICAL_SLIDING ){
00338         
00339                 // DO NOTHING: Code is integrated in updateEuler
00340         
00341         } else {
00342                 // NO COLLISION DETECTION
00343         }
00344         
00345         accu += timestep;
00346         if( flagExportToPrecalculatedFile ){
00347                 if( accu >= COME::flagAnimationResolution ){
00348                         saveJointIntoFile( getPatient(0)->getRootJoint( ), UPDATE );
00349                         accu = 0.0;
00350                 }
00351         }
00352         
00353         //unlock();
00354 
00355         return true;
00356 }
00357 
00361 void
00362 COME_Scenario::createPrecalculatedFiles() {
00363 
00364         list<COME_Patient*>::iterator iter;
00365         for (iter = patientList.begin(); !(iter == patientList.end()); iter++){
00366                 list<COME_BioStructure*>::const_iterator iterOrgans;
00367                 for( iterOrgans = ((list<COME_BioStructure *>*)(*iter)->getPtOrganList())->begin(); iterOrgans != ((list<COME_BioStructure *>*)(*iter)->getPtOrganList())->end(); iterOrgans++  ){
00368                         //if( !((COME_MoleculesCartilage*)(*iterOrgans))->isFixed() )
00369                                 ((COME_MoleculesCartilage*)(*iterOrgans))->getTissue()->saveDeformationFile( NEW ); // to create empty file
00370                                 //((COME_MoleculesCartilage*)(*iterOrgans))->getTissue()->saveDeformationFile( UPDATE ); // to save zero position
00371                 }
00372         }
00373 
00374         if( getPatient(0) ){
00375                 if( getPatient(0)->getRootJoint( )->getDescription() == "root" ){
00376                         saveJointIntoFile( getPatient(0)->getRootJoint( ), NEW ); // to create empty file
00377                         //saveJointIntoFile( getPatient(0)->getRootJoint( ), UPDATE ); // to save zero position
00378                 }
00379         }
00380 
00381 }
00382 
00383 
00389 void
00390 COME_Scenario::saveJointIntoFile( COME_Joint *jointAux, int fileMode ) {
00391         
00392         // Save axes and angles into file
00393         FILE *file;
00394         string strMode;
00395         string fileName = COME::baseFolder + "/" + jointAux->getDescription() + "angles.dat";
00396         
00397         if( fileMode == NEW ){
00398                 strMode = "wb";// create new
00399                 fclose( fopen( fileName.c_str(), strMode.c_str() ) );
00400         } else {
00401                 strMode = "a+b"; // update at the end
00402 
00403                 if( file = fopen( fileName.c_str(), strMode.c_str() ) ){
00404                 
00405                         double *matrixAux[4];
00406                         matrixAux[0] = new double[4];
00407                         matrixAux[1] = new double[4];
00408                         matrixAux[2] = new double[4];
00409                         matrixAux[3] = new double[4];
00410                         
00411                         jointAux->getLim().getMatrixD( matrixAux ); 
00412                 
00413                         double array[16];
00414                         array[0] = ((COME_Simulator*)parent)->getClock();
00415                         
00416                         array[1] = matrixAux[0][0];
00417                         array[2] = matrixAux[1][0];
00418                         array[3] = matrixAux[2][0];
00419                 
00420                         array[4] = matrixAux[0][1];
00421                         array[5] = matrixAux[1][1];
00422                         array[6] = matrixAux[2][1];
00423                 
00424                         array[7] = matrixAux[0][2];
00425                         array[8] = matrixAux[1][2];
00426                         array[9] = matrixAux[2][2];
00427                 
00428                         array[10] = matrixAux[3][0];
00429                         array[11] = matrixAux[3][1];
00430                         array[12] = matrixAux[3][2];
00431                         
00432                         array[13] = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMin() ) * ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getCurrent() + ((COME_PolyaxialJoint*)jointAux )->getDofList()[FLEX]->getMin();
00433                 
00434                         array[14] = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMin() ) * ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getCurrent() + ((COME_PolyaxialJoint*)jointAux )->getDofList()[ADDUCT]->getMin();
00435                 
00436                         array[15] = ( ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMax() - ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMin() ) * ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getCurrent() + ((COME_PolyaxialJoint*)jointAux )->getDofList()[TWIST]->getMin();
00437                 
00438                         fwrite( array, sizeof(double), 16, file );
00439                         
00440                         fclose(file);
00441                 }
00442         }
00443         
00444         list<COME_Joint*>::iterator childAux;
00445         int childCount = 0;
00446         for( childAux = jointAux->getChildListPt()->begin(); childCount < jointAux->getChildListPt()->size();  childAux++, childCount++ ){
00447                 saveJointIntoFile( (*childAux), fileMode );
00448         }
00449 }
00450 
00455 bool
00456 COME_Scenario::loadFile( string file_name, COME_Simulator *simulator ){
00457 
00458         COME_Xml doc;
00459 
00460         // load scene 
00461         doc.loadSceneFile( file_name, patientList, simulator, this );
00462         
00463         vector<COME_Material*>* newMaterials = doc.getAllMaterials();
00464         for( int i = 0; i < newMaterials->size(); i++ ){
00465                 availableMaterials->push_back( (*newMaterials)[i] );
00466         }
00467                 
00468         // load movements
00469         movement = new COME_Movement();
00470         movement->vpLoadMotionFile( file_name.c_str() );
00471         movement->setTime_Err( simulator->getTimestep() / 2.0 );
00472         
00473         parent = simulator;
00474         simulator->setScene( this );
00475         
00477         for( int iS = 0; iS < getPatient(0)->getPtOrganList()->size(); iS++ ){
00478 
00479                 for( int iI = 0; iI < getPatient(0)->getOrgan(iS)->internalsDesc.size(); iI++ ){
00480                 
00481                         getPatient(0)->getOrgan(iS)->internals.push_back( getPatient(0)->getOrgan(getPatient(0)->getOrgan(iS)->internalsDesc[iI] ) );
00482                 }
00483         }
00484         
00485         modified = true;
00486 
00487         return true;
00488 }
00489 
00490 
00494 void
00495 COME_Scenario::initializeAllHashs(){
00496 
00497         for( int iS = 0; iS < getPatient(0)->getPtOrganList()->size(); iS++ ){
00498                 if( getPatient(0)->getOrgan(iS)->internalsDesc.size() > 0 ){
00499                         getPatient(0)->getOrgan(iS)->initializeHashStructureRayCasting();
00500                         printf( "Organ %s initialized for spherical sliding. \n", getPatient(0)->getOrgan(iS)->getDescription().c_str() );
00501                 }
00502         }
00503 }
00504 
00509 bool
00510 COME_Scenario::addOrganFromFile( string file_name, COME_Simulator *simulator ){
00511 
00512         COME_Xml doc( availableMaterials );
00513 
00514         // load scene 
00515         COME_MoleculesCartilage *newOrgan = (COME_MoleculesCartilage*)doc.loadOrganFile( file_name, patientList, simulator, this );
00516         
00517         //((COME_MoleculesCartilage*)newOrgan)->initializeSkinning();
00518         //printf( "skinning initialized \n" );
00519         //newOrgan->updateSkin();
00520         //printf( "skin updated \n" );
00521         
00522         // Following lines (adding new materials is not necessary because
00523         // availableMaterials passed above is already updated in loadOrganFile
00524         /*vector<COME_Material*>* newMaterials = doc.getAllMaterials();
00525         for( int i = 0; i < newMaterials->size(); i++ ){
00526                 availableMaterials->push_back( (*newMaterials)[i] );
00527         }*/
00528 
00529         newOrgan->getTissue()->setMoleculesMaterial( (*availableMaterials)[0] );
00530         getPatient(0)->addOrgan( newOrgan );
00531 
00532         parent = simulator;
00533         simulator->setScene( this );
00534         
00535         modified = true;
00536         return true;
00537 }
00538 
00543 bool
00544 COME_Scenario::addMaterialsFromFile( string file_name, COME_Simulator *simulator ){
00545 
00546         COME_Xml doc;
00547 
00548         // load scene 
00549         doc.loadMaterialsFile( file_name, patientList, simulator, this );
00550         
00551         vector<COME_Material*>* newMaterials = doc.getAllMaterials();
00552         for( int i = 0; i < newMaterials->size(); i++ ){
00553                 availableMaterials->push_back( (*newMaterials)[i] );
00554         }
00555 
00556         modified = true;
00557         return true;
00558 }
00559 
00564 bool
00565 COME_Scenario::saveFile( string file_name, COME_Simulator *simulator ){
00566 
00567         COME_Xml docS;
00568         docS.saveSceneFile( file_name, patientList, simulator );
00569 
00570         return true;
00571 }
00572 
00573 
00576 /*bool
00577 COME_Scenario::loadOrgans( vector<string> organFiles, int mode ){
00578 
00579         double xResolution = 0.01;
00580         double yResolution = 0.01;
00581         double zResolution = 0.01;
00582         double scaleFactorFile = 0.005; //-00
00583 
00584 
00586         double fracture_dist = (xResolution+yResolution+zResolution)/3 + (xResolution+yResolution+zResolution)/6;
00587         double friction_const = 0.3;
00588         double radius = fracture_dist/2;
00589 
00590         for( int iF = 0; iF < organFiles.size(); iF++ ){
00591 
00592                 COME_MoleculesCartilage *newOrgan = new COME_MoleculesCartilage();
00593                 newOrgan->setParent( patientList.front() );
00594                 
00595                 COME_Mesh *organMainMesh;
00596                 //organMainMesh.loadFile( organFiles[iF], scaleFactorFile );
00597                 organMainMesh->loadFile( organFiles[iF], scaleFactorFile, COME_Point3D( -97.9, 4.5, -25.7 ) );
00598                 
00599                 if( mode == MARCHING_CUBES ){
00600 
00601                         COME_Point3D mins, maxs;
00602                         organMainMesh->getEnvelop( mins, maxs );
00603                         // Marchig cubes to create balls
00604                         for( double xMcubes = mins.getX(); xMcubes <= maxs.getX(); xMcubes += xResolution ){
00605                                 for( double yMcubes = mins.getY(); yMcubes <= maxs.getY(); yMcubes += yResolution ){
00606                                         for( double zMcubes = mins.getZ(); zMcubes <= maxs.getZ(); zMcubes += zResolution ){
00607 
00608                                                 COME_Point3D currentPoint( xMcubes, yMcubes, zMcubes );
00609                                                 if( organMainMesh->isInside( currentPoint ) ){
00610                                                         
00611                                                         // Attention: Using allways the same material
00612                                                         COME_Molecule *molecule = new COME_Molecule( currentPoint, radius, friction_const, (*availableMaterials)[0], COME_Vector3D(), (COME_BioStructure*)newOrgan );
00613                                                         molecule->setFixed( true );
00614                                                         newOrgan->getTissue()->addMolecule( molecule );
00615 
00616                                                 }
00617                                                 printf( "." );
00618                                         }
00619                                         printf( "\n" );
00620                                 }
00621                                 printf( "\n" );
00622                         }
00623                 } else if( mode == VERTICES ){
00624 
00625                         double layersThickness = 0.045;
00626                         double minDist = 0.016;
00627                         radius = 0.02;
00628                         
00629                         int layersNumber = 0;
00630                         COME_Mesh *organMesh = organMainMesh;
00631                         
00632                         for( int numFaces = 1; numFaces >= 1;  ){
00633 
00634                                 //organMesh.scaleOnNormals( layersThickness );
00635                                 organMesh->scale( 0.85 ); //1.12 for pelvis cartilage
00636                                 //numFaces = organMesh.simplify( minDist );
00637                                 //printf( "Remeshed to %d faces.\n", numFaces );
00638                                 //getchar();
00639 
00640                                 if( numFaces > 0 ){
00641                                         // Create balls on each vertex of the mesh
00642                                         for( int iVert = 0; iVert < organMesh->getVertices().size(); iVert++ ){
00643 
00644                                                 COME_Point3D pointPosition = COME_Point3D( organMesh->getAVertex( iVert ));
00645                                                 COME_Molecule *molecule = new COME_Molecule( pointPosition, radius, friction_const, (*availableMaterials)[0], COME_Vector3D(), newOrgan );
00646                                                 if( layersNumber == 1 ){
00647 
00648                                                         molecule->setMaterial( (*availableMaterials)[1] );
00649 
00650                                                 }
00651                                                 newOrgan->getTissue()->addMolecule( molecule, minDist );
00652                                         }
00653                                 }
00654                                 layersNumber++;
00655                                 if( layersNumber == 2 )
00656                                         numFaces = 0;
00657                         }
00658                 }
00659 
00660                 newOrgan->setSurface( organMainMesh );
00661                 patientList.front()->addOrgan( newOrgan );
00662                 printf( "Finished loading organ.\n" );
00663         }
00664 
00665         return true;
00666 }
00667 */
00668 bool
00669 COME_Scenario::loadOrganMesh( string organFile, int type, double scaleFactor ){
00670 
00671         COME_BioStructure *newOrgan;
00672 
00673         switch( type ){
00674                 case CARTILAGE  : newOrgan = new COME_MoleculesCartilage(); newOrgan->setFixed( false ); break;
00675                 case LIGAMENT   : newOrgan = new COME_MoleculesCartilage(); newOrgan->setFixed( false ); break;
00676                 case BONE       : newOrgan = new COME_MoleculesBone(); newOrgan->setFixed( true ); break;
00677                 case CLAMP      : newOrgan = new COME_Clamp(); newOrgan->setFixed( true ); break;
00678                 //case LIGAMENT : newOrgan = new COME_MoleculesLigament(); break;
00679         }
00680         newOrgan->setType( type );
00681         newOrgan->setParent( patientList.front() );
00682         patientList.front()->addOrgan( newOrgan );
00683 
00684         COME_Mesh *organMainMesh = new COME_Mesh();
00685         // translation for pubofemoral to RNW04
00686         //organMainMesh->loadFile( organFile, scaleFactor, COME_Point3D( 21.0, -4.5, 8.76 ) );
00687         // translation for pubofemoral to DEMO
00688         //organMainMesh->loadFile( organFile, scaleFactor, COME_Point3D( 21.0, -4.5, 8.76 ) );
00689         
00690         // APPLY ALL TRANSLATIONS
00691         COME_Point3D oldhip( -97.9, 4.5, -25.7 );
00692         //COME_Point3D labrum( -119.8, 8.535, -34.84 );
00693         //COME_Point3D pubofemoral( -118.9, 9, -34.46 );
00694         COME_Point3D labrum( 21.9, -4.035, 9.14 );
00695         COME_Point3D pubofemoral( 21.0, -4.5, 8.76 );
00696         COME_Point3D comedemo( 0.0, 0.0, -2.5 );
00697         
00698         string descOrgan = "fuck"; //"comedemo"; //////////////////////////////////////////// CHANGE THE STRING HERE TO SWITCH TRANSLATIONS
00699         if( descOrgan == "labrum")
00700                 organMainMesh->loadFile( organFile, scaleFactor, labrum );
00701         else if( ( descOrgan == "pubofemoral") || ( descOrgan == "ischiofemoral") )
00702                 organMainMesh->loadFile( organFile, scaleFactor, pubofemoral );
00703         else if( descOrgan == "pelviscart")
00704                 organMainMesh->loadFile( organFile, scaleFactor, oldhip );
00705         else if( descOrgan == "pelvisbone")
00706                 organMainMesh->loadFile( organFile, scaleFactor, oldhip );
00707         else if( descOrgan == "femurcart")
00708                 organMainMesh->loadFile( organFile, scaleFactor, oldhip );
00709         else if( descOrgan == "femurbone")
00710                 organMainMesh->loadFile( organFile, scaleFactor, oldhip );
00711         else if( descOrgan == "comedemo")
00712                 organMainMesh->loadFile( organFile, scaleFactor, comedemo );
00713         else {
00714                 organMainMesh->loadFile( organFile, scaleFactor, COME_Point3D( 0, 0, 0 ) );
00715         }
00716         //organMainMesh->loadFile( organFile, scaleFactor, COME_Point3D( 0.0, 0.0, 0.0 ) );
00717         
00718         newOrgan->setSurface( organMainMesh );
00719         switch( type ){
00720                 case CARTILAGE  : ((COME_MoleculesCartilage*)newOrgan)->initializeSkinning(); break;
00721                 case LIGAMENT   : ((COME_MoleculesCartilage*)newOrgan)->initializeSkinning(); break;
00722                 case BONE       : ((COME_MoleculesBone*)newOrgan)->initializeSkinning();break;
00723                 case CLAMP      : ((COME_Clamp*)newOrgan)->initializeSkinning();break;
00724                 //case LIGAMENT : ((COME_MoleculesLigament*)newOrgan)->initializeSkinning() break;
00725         }
00726         printf( "skinning initialized \n" );
00727         newOrgan->updateSkin();
00728         printf( "skin updated \n" );
00729         
00730         return true;
00731 }
00732 
00738 void
00739 COME_Scenario::getEnvelop( COME_Point3D& mins, COME_Point3D& maxs ) const {
00740 
00741         patientList.front()->getEnvelop( mins, maxs );
00742 
00743         list<COME_Patient*>::const_iterator iter;
00744         for (iter = patientList.begin(); !(iter == patientList.end()); iter++){
00745                 
00746                 COME_Point3D minsP , maxsP;
00747                 (*iter)->getEnvelop( minsP, maxsP );
00748                 
00749                 if( minsP.getX() < mins.getX() ){
00750                         mins.setX( minsP.getX() );
00751                 }
00752                 if( minsP.getY() < mins.getY() ){
00753                         mins.setY( minsP.getY() );
00754                 }
00755                 if( minsP.getZ() < mins.getZ() ){
00756                         mins.setZ( minsP.getZ() );
00757                 }
00758 
00759                 if( maxsP.getX() > maxs.getX() ){
00760                         maxs.setX( maxsP.getX() );
00761                 }
00762                 if( maxsP.getY() > maxs.getY() ){
00763                         maxs.setY( maxsP.getY() );
00764                 }
00765                 if( maxsP.getZ() > maxs.getZ() ){
00766                         maxs.setZ( maxsP.getZ() );
00767                 }
00768         }
00769 
00770 }
comescenario.cpp
