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........: August/05/2002
00028 //  DESCRIPTION.: Class declaration.
00029 //
00031 
00032 
00033 #include <bio/comemoleculesbone.h>
00034 #include <bio/comematerial.h>
00035 #include <math.h>
00036 #include <algebra/comemesh.h>
00037 
00038 //#include <LinAlg/linalg.h>
00039 //using namespace LinAlg;
00040 //#include <MarchingCubes/marching_cubes.h>
00041 //#include <MarchingCubes/blob.h>
00042 
00043 
00044 
00048 
00049 COME_MoleculesBone::COME_MoleculesBone(){
00050 
00051         tissue = new COME_MoleculesTissue();
00052         tissue->setParent(this);
00053         surface = NULL;
00054         hashTable = NULL;
00055         updateMassCenter();
00056         globalHooke = 0.0;
00057         collisionRadius = 1.0;
00058 }
00059 
00060 COME_MoleculesBone::~COME_MoleculesBone(){
00061 
00062         
00063 }
00064 
00068 void
00069 COME_MoleculesBone::setTissue( COME_MoleculesTissue *tissueN ){
00070 
00071         tissue = tissueN;
00072 }
00073 
00074 void
00075 COME_MoleculesBone::setMassCenter( const COME_Point3D& massC ){
00076 
00077         massCenter = massC;
00078 }
00079 
00080 COME_MoleculesTissue*
00081 COME_MoleculesBone::getTissue(){
00082 
00083         return tissue;
00084 }
00085 
00086 COME_Point3D
00087 COME_MoleculesBone::getMassCenter() const{
00088 
00089         return massCenter;
00090 }
00091 
00092 //Recalculates the Mesh based on the new positions of anchoring the molecules
00093 void 
00094 COME_MoleculesBone:: updateSkin(){
00095         if( surface )
00096                 surface->updateSkin();
00097 }
00098 
00102 void
00103 COME_MoleculesBone::updateSurface( bool regenerate ){
00104 
00105 //      if( regenerate ){
00106 //              updateSurface( 1.0, 150.0, 0.005 ); //1.0  120.0 0.005  
00107 //              surface.estabilishLinkVerticesFaces();
00108 
00109                 // para os pequenos 1.0, 80.0, 0.005
00110                 // para os grandes 1.0, 160.0, 0.01
00111                 //collision:   //2.3  (15)
00112 
00113 //      } else {
00114                 //generateEnvelopSurface();
00115 //      }
00116   
00117 
00118         if( regenerate ){
00119                 COME_Point3D mins, maxs;
00120                 getEnvelop( mins, maxs );
00121                 double volume = ( maxs.getX() - mins.getX() ) * ( maxs.getY() - mins.getY() ) * ( maxs.getZ() - mins.getZ() );
00122                 double radiusAvg = tissue->getRadiusAverage();
00123 
00124                 updateSurface( 1.0, 1/(radiusAvg/2), sqrt( volume )/20.0 + 0.005 );
00125                 surface->estabilishLinkVerticesFaces();
00126 
00127         } else {
00128                 surface->update();
00129         }
00130 }
00131 
00132 
00137 void
00138 COME_MoleculesBone::updateSurface( double skinSphere, double stiff, double resolution ){
00139 
00140 /*
00141         // Metaballs represented as the spheres with Raquel's density function
00142         MetaballsSphere metaSphere( skinSphere );
00143         
00144         vector<Primitive*>                      primitives;
00145         VEC                                                     C(3);           // sphere center
00146 
00147         list<COME_Molecule*>::iterator iter;    
00148         for( iter = tissue->getShape()->begin(); iter != tissue->getShape()->end(); iter++ ){
00149                 //read spheres
00150                 COME_Point3D point      = (*iter)->getPosition();
00151 
00152                 C[0] = point.getX();
00153                 C[1] = point.getY();
00154                 C[2] = point.getZ();
00155                 
00156                 //add sphere  
00157                 metaSphere.addSphere( Sphere( C, (*iter)->getRadius(), stiff ) );
00158         }
00159   
00160         //add to primitives
00161         primitives.push_back( &metaSphere );
00162 
00163         //create surface
00164         marchingCubes surfaceTriangles( primitives, resolution );
00165         
00166 
00167         // copy vertices and triangle indices into the Class Mesh       
00168         
00169         int iPos;
00170 
00171         vector<VEC>     vertexList      = surfaceTriangles.getVertices();
00172         vector<VEC> normalList  = surfaceTriangles.getNormals();
00173 
00174         // copy all vertices and the normals of the triangulated surface 
00175         surface->resetMesh();
00176         for ( iPos = 0; iPos < vertexList.size(); iPos++ ){
00177                 
00178                 COME_Vector3D normal = COME_Vector3D( normalList[ iPos ][ 0 ], normalList[ iPos ][ 1 ], normalList[ iPos ][ 2 ] );
00179                 surface->addVertex( COME_Vertex3D( vertexList[ iPos ][ 0 ], vertexList[ iPos ][ 1 ], vertexList[ iPos ][ 2 ], normal ) );
00180         }
00181 
00182         //copy all indices of vertices 
00183         vector<Triangle>        trianglesList   = surfaceTriangles.getTriangles();
00184         for (vector<Triangle>::const_iterator tri=trianglesList.begin(); tri!=trianglesList.end(); tri++){     
00185                 
00186                 vector<int>                     indicesTriangles;
00187 
00188                 indicesTriangles.push_back( tri->indices[0] );
00189                 indicesTriangles.push_back( tri->indices[1] );
00190                 indicesTriangles.push_back( tri->indices[2] );
00191 
00192                 surface->addFaceAsIndices( indicesTriangles );
00193         }
00194 */
00195 }
00196 
00201 void
00202 COME_MoleculesBone::generateEnvelopSurface(){
00203         
00204         surface->resetMesh();
00205 
00206         COME_Point3D mins, maxs;
00207         getEnvelop( mins, maxs );
00208 
00209         COME_Vertex3D v0 = COME_Vertex3D( mins.getX(), mins.getY(), mins.getZ() );
00210         COME_Vertex3D v1 = COME_Vertex3D( mins.getX(), maxs.getY(), mins.getZ() );
00211         COME_Vertex3D v2 = COME_Vertex3D( mins.getX(), mins.getY(), maxs.getZ() );
00212         COME_Vertex3D v3 = COME_Vertex3D( mins.getX(), maxs.getY(), maxs.getZ() );
00213         COME_Vertex3D v4 = COME_Vertex3D( maxs.getX(), mins.getY(), maxs.getZ() );
00214         COME_Vertex3D v5 = COME_Vertex3D( maxs.getX(), maxs.getY(), maxs.getZ() );
00215         COME_Vertex3D v6 = COME_Vertex3D( maxs.getX(), mins.getY(), mins.getZ() );
00216         COME_Vertex3D v7 = COME_Vertex3D( maxs.getX(), maxs.getY(), mins.getZ() );
00217 
00218         v0.setNormal( COME_Vector3D( v0 - v5 ).vpNormalize() );
00219         v1.setNormal( COME_Vector3D( v1 - v4 ).vpNormalize() );
00220         v2.setNormal( COME_Vector3D( v2 - v7 ).vpNormalize() );
00221         v3.setNormal( COME_Vector3D( v3 - v6 ).vpNormalize() );
00222         v4.setNormal( COME_Vector3D( v4 - v1 ).vpNormalize() );
00223         v5.setNormal( COME_Vector3D( v5 - v0 ).vpNormalize() );
00224         v6.setNormal( COME_Vector3D( v6 - v3 ).vpNormalize() );
00225         v7.setNormal( COME_Vector3D( v7 - v2 ).vpNormalize() );
00226         
00227         surface->addVertex( v0 );
00228         surface->addVertex( v1 );
00229         surface->addVertex( v2 );
00230         surface->addVertex( v3 );
00231         surface->addVertex( v4 );
00232         surface->addVertex( v5 );
00233         surface->addVertex( v6 );
00234         surface->addVertex( v7 );
00235 
00236         vector<int>     indTri;
00237 
00238         indTri.push_back( 0 ); indTri.push_back( 3 ); indTri.push_back( 1 );
00239         surface->addFaceAsIndices( indTri );
00240         indTri.clear();
00241         indTri.push_back( 0 ); indTri.push_back( 2 ); indTri.push_back( 3 );
00242         surface->addFaceAsIndices( indTri );
00243         indTri.clear();
00244         indTri.push_back( 2 ); indTri.push_back( 5 ); indTri.push_back( 3 );
00245         surface->addFaceAsIndices( indTri );
00246         indTri.clear();
00247         indTri.push_back( 2 ); indTri.push_back( 4 ); indTri.push_back( 5 );
00248         surface->addFaceAsIndices( indTri );
00249         indTri.clear();
00250         indTri.push_back( 4 ); indTri.push_back( 7 ); indTri.push_back( 5 );
00251         surface->addFaceAsIndices( indTri );
00252         indTri.clear();
00253         indTri.push_back( 4 ); indTri.push_back( 6 ); indTri.push_back( 7 );
00254         surface->addFaceAsIndices( indTri );
00255         indTri.clear();
00256         indTri.push_back( 6 ); indTri.push_back( 1 ); indTri.push_back( 7 );
00257         surface->addFaceAsIndices( indTri );
00258         indTri.clear();
00259         indTri.push_back( 6 ); indTri.push_back( 0 ); indTri.push_back( 1 );
00260         surface->addFaceAsIndices( indTri );
00261         indTri.clear();
00262         indTri.push_back( 1 ); indTri.push_back( 5 ); indTri.push_back( 7 );
00263         surface->addFaceAsIndices( indTri );
00264         indTri.clear();
00265         indTri.push_back( 1 ); indTri.push_back( 3 ); indTri.push_back( 5 );
00266         surface->addFaceAsIndices( indTri );
00267         indTri.clear();
00268         indTri.push_back( 6 ); indTri.push_back( 2 ); indTri.push_back( 0 );
00269         surface->addFaceAsIndices( indTri );
00270         indTri.clear();
00271         indTri.push_back( 6 ); indTri.push_back( 4 ); indTri.push_back( 2 );
00272         surface->addFaceAsIndices( indTri );
00273         indTri.clear();
00274 
00275 }
00276 
00280 void
00281 COME_MoleculesBone::getEnvelop( COME_Point3D& mins, COME_Point3D& maxs ){
00282 
00283         tissue->getEnvelop( mins, maxs );
00284 }
00285 
00295 
00296 bool
00297 COME_MoleculesBone::update( double timestep, double simClock ){
00298 
00299 /*
00300         double xh , hh, h6;
00301         COME_Vector3D *dymV, *dymF, *dymT, *dymW;
00302         
00303         COME_Vector3D *dytV, *dytF, *dytT, *dytW;
00304         
00305         COME_Vector3D *ytX, *ytL, *ytP;
00306         COME_Matrix *ytR;
00307 
00308         dymV = new COME_Vector3D();
00309         dymF = new COME_Vector3D();
00310         dymT = new COME_Vector3D();
00311         dymW = new COME_Vector3D();
00312         
00313         dytV = new COME_Vector3D();
00314         dytF = new COME_Vector3D();
00315         dytT = new COME_Vector3D();
00316         dytW = new COME_Vector3D();
00317         
00318         ytX = new COME_Vector3D();
00319         ytL = new COME_Vector3D();
00320         ytP = new COME_Vector3D();
00321         ytR = new COME_Matrix();
00322 
00323         COME_Vector3D momentumTemp;
00324         COME_Vector3D w;
00325 
00326         //COME_Matrix I = COME_Matrix( rotationMatrix * inertiaTensor * rotationMatrix.getTransposed() );
00327 
00328         hh = timestep * 0.5;
00329         h6 = timestep / 6.0;
00330         xh = simClock + hh;
00331         
00332         // First step //////////////////////////////////////////////////////
00333         
00334         ytX->setVector3D( ( linearVelocity * hh ) + massCenter );
00335         ytL->setVector3D( rotationMatrix * inertiaTensor * rotationMatrix.getTransposed() * angularVelocity );
00336         momentumTemp.setVector3D( *ytL ); // save initial momentum
00337         ytP->setVector3D( linearVelocity * tissue->getMass() );
00338         ytR->setMatrix( rotationMatrix );
00339         
00340         // Second step //////////////////////////////////////////////////////
00341         derivs( xh, hh, ytX, ytL, ytP, ytR, dytV, dytF, dytT, dytW ); 
00342 
00343         ytX->setVector3D( ( *dytV * hh ) + massCenter );
00344         ytL->setVector3D( *dytW * *ytR * inertiaTensor * (*dytW * ytR->getTransposed() ) * *dytW );
00345         ytP->setVector3D( *dytV * tissue->getMass() );
00346         ytR->setMatrix( rotationMatrix );//*ytR += (*dytW); // probably non-correct
00347 
00348         // Third step //////////////////////////////////////////////////////
00349         derivs( xh, hh, ytX, ytL, ytP, ytR, dymV, dymF, dymT, dymW ); 
00350 
00351         ytX->setVector3D( ( *dymV * timestep ) + massCenter );
00352         ytL->setVector3D( *dymW * *ytR * inertiaTensor * (*dymW * ytR->getTransposed() ) * *dymW );
00353         ytP->setVector3D( *dymV * tissue->getMass() );
00354         ytR->setMatrix( rotationMatrix );//ytR->setMatrix( *( dymWR->add( *ytR ) ) );
00355 
00356         *dymV += *dytV;
00357         *dymF += *dytF;
00358         *dymT += *dytT;
00359         *dymW += *dytW;
00360 
00361 
00362         // Fourth step //////////////////////////////////////////////////////
00363         derivs( simClock + timestep, timestep, ytX, ytL, ytP, ytR, dytV, dytF, dytT, dytW ); 
00364 
00365         massCenter = ( ( linearVelocity + *dytV + ( *dymV * 2.0 ) ) * h6 ) + massCenter;
00366         COME_Vector3D rotatedVec = COME_Vector3D( angularVelocity * rotationMatrix );
00367         rotationMatrix = ( ( angularVelocity + (*dytW * *ytR) + ( *dymW * *ytR * 2.0 ) ) * h6 ).transMult( rotatedVec );
00368         angularVelocity =  rotationMatrix * inertiaTensor.getInverse() * rotationMatrix.getTransposed() * momentumTemp;
00369         linearVelocity = ( ( linearVelocity + *dytV + ( *dymV * 2.0 ) ) * h6 ) * tissue->getMass();
00370                 
00371         delete dymV;
00372         delete dymF;
00373         delete dymT;
00374         delete dymW;
00375         
00376         delete dytV; 
00377         delete dytF;
00378         delete dytT;
00379         delete dytW;
00380 
00381         delete ytX;
00382         delete ytL;
00383         delete ytP;
00384         delete ytR;
00385 */
00386         updateSkin();
00387         return true;
00388 }
00389 
00397 void
00398 COME_MoleculesBone::derivs( double currTime, double timestep,
00399                                                     COME_Point3D *X, COME_Vector3D *L, COME_Vector3D *P, COME_Matrix *R,
00400                                                         COME_Vector3D *V, COME_Vector3D *F, COME_Vector3D *T, COME_Vector3D *W ){
00401         
00402         if( !isFixed() ){
00403                 
00404                 V->setVector3D( *P / tissue->getMass() );
00405                 //F = addAllForces();
00406                 //T = ( forcePoint - massCenter ) * forceVector;
00407                 W->setVector3D(  (*R) * inertiaTensor.getInverse() * R->getTransposed() * (*L) );
00408         }
00409 }
00410 
00411 
00417 void
00418 COME_MoleculesBone::updateInertiaTensor( double discretization ){
00419 
00420 
00421         inertiaTensor.reset();
00422         int hits = 0;
00423         
00424         COME_Point3D mins, maxs;
00425         getEnvelop( mins, maxs );
00426 
00427         for( double x = mins.getX(); x < maxs.getX(); x += discretization *  ( maxs.getX() - mins.getX() ) ){
00428                 for( double y = mins.getY(); y < maxs.getY(); y += discretization *  ( maxs.getY() - mins.getY() ) ){
00429                         for( double z = mins.getZ(); z < maxs.getZ(); z += discretization *  ( maxs.getZ() - mins.getZ() ) ){
00430 
00431                                 if( surface->isInside( COME_Point3D(x,y,z) ) ){
00432 
00433                                         inertiaTensor.setValueAt( 0, 0, inertiaTensor.getValueAt( 0, 0 ) + ( y * y + z * z ) );
00434                                         inertiaTensor.setValueAt( 0, 1, inertiaTensor.getValueAt( 0, 1 ) + ( -x * y ) );
00435                                         inertiaTensor.setValueAt( 0, 2, inertiaTensor.getValueAt( 0, 2 ) + ( -x * z ) );
00436 
00437                                         inertiaTensor.setValueAt( 1, 0, inertiaTensor.getValueAt( 1, 0 ) + ( -x * y ) );
00438                                         inertiaTensor.setValueAt( 1, 1, inertiaTensor.getValueAt( 1, 1 ) + ( x * x + z * z ) );
00439                                         inertiaTensor.setValueAt( 1, 2, inertiaTensor.getValueAt( 1, 2 ) + ( -y * z ) );
00440 
00441                                         inertiaTensor.setValueAt( 2, 0, inertiaTensor.getValueAt( 2, 0 ) + ( -x * z ) );
00442                                         inertiaTensor.setValueAt( 2, 1, inertiaTensor.getValueAt( 2, 1 ) + ( -y * z ) );
00443                                         inertiaTensor.setValueAt( 2, 2, inertiaTensor.getValueAt( 2, 2 ) + ( x * x + y * y ) );
00444 
00445                         hits++;
00446                                 }
00447                         }
00448                 }
00449         }
00450 
00451         inertiaTensor.multiplyScalar( 1/hits );
00452         inertiaTensor.multiplyScalar( tissue->getMass() );
00453 
00454 }
00455 
00460 void
00461 COME_MoleculesBone::updateMassCenter(){
00462 
00463         massCenter = tissue->getMassCenter();
00464 }
00465 
00469 void
00470 COME_MoleculesBone::initializeSkinning(){
00471 
00472         if( this->getSurface() == NULL ){
00473                 surface = new COME_Mesh();
00474                 surface->setParent( this );
00475                 updateSurface( IMPLICIT );
00476                 printf( "update IMPLICIT BONE done \n" );
00477         }
00478         this->getSurface()->anchorToMolecules( this->getTissue(), 4 );
00479         cout << "initialising Faces Neighbours...";
00480         this->getPointerToSurface()->initFacesNeighbours();
00481         cout << "done \ninitialising Vertices Neighbours...";
00482         this->getPointerToSurface()->initVerticesNeighbours();
00483         cout << "done\n";
00484 }
00485 
00486 void
00487 COME_MoleculesBone::respondCollision(){
00488 
00489         // do nothing for now
00490         //printf( "Entrou bone::respondcollision \n" );
00491 }
comemoleculesbone.cpp
