Biomechanical Joint Model  Author: Anderson Maciel

---

comebiostructure.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2002 by Anderson Maciel                                 *
 *   andi.maciel@gmail.com                                                 *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 **************************************************************************/
//
//  PROJECT.....: CO-ME
//  RESPONSIBLE.: 
//
//  FILE........: comebiostructure.cpp
//  DESCRIPTION.: 
//                                 
//
//  AUTHOR......: Anderson Maciel
//  DATE........: May/27/2002
//  DESCRIPTION.: Class definition.
//

#include        <bio/comebiostructure.h>
#include        <general/comescenario.h>
#include        <algebra/comemesh.h>
#include        <RAPID/RAPID.H>

#include <math.h>
#include <time.h>


void
COME_BioStructure::setSystem( int systemN ){

        system = systemN;
}

void
COME_BioStructure::setSurface( COME_Mesh *surfaceN ){

        surface = surfaceN;
        surface->setParent( this );
        surface->setCollisionDetectableFaces();
}

void
COME_BioStructure::setFixed( bool flag ){

        fixed = flag;
}

void
COME_BioStructure::setLIM( const COME_Matrix &limN ){

        localInstanceMatrix = limN;
}

void
COME_BioStructure::setGIM( const COME_Matrix &gimN ){

        globalInstanceMatrix = gimN;
}

void
COME_BioStructure::setGlobalHooke( double hookeN ){

        globalHooke = hookeN;
}

void
COME_BioStructure::setCollisionRadius( double rad ){

        collisionRadius = rad;
}

void
COME_BioStructure::setType( int nType ){
        
        type = nType;
}

        
int
COME_BioStructure::getSystem() const {

        return system;
}

COME_Mesh*
COME_BioStructure::getSurface() {

        return surface;
}

COME_Mesh*
COME_BioStructure::getPointerToSurface() {

        return surface;
}

bool
COME_BioStructure::isFixed() const{

        return fixed;
}

const COME_Matrix&
COME_BioStructure::getLIM() const {

        return localInstanceMatrix;
}

const COME_Matrix&
COME_BioStructure::getGIM() const {

        return globalInstanceMatrix;
}


double
COME_BioStructure::getGlobalHooke() const {

        return globalHooke;
}

double
COME_BioStructure::getCollisionRadius() const {

        return collisionRadius;
}

double
COME_BioStructure::getTopStress() {

        return getTissue()->getTopStress();
}

void
COME_BioStructure::scale( double factorx, double factory, double factorz ){

        getTissue()->scale( factorx, factory, factorz );
        
        if( surface ){
                if( isFixed() ){
                        surface->scale( factorx, factory, factorz );
                } else {
                        getTissue()->makeAllLocalFrames();
                        surface->updateSkin();
                }
        }
        
}

void
COME_BioStructure::translate( double px, double py, double pz ){

        getTissue()->translate( px, py, pz );
        
        if( surface ){
                if( isFixed() ){
                        surface->translate( px, py, pz );
                } else {
                        getTissue()->makeAllLocalFrames();
                        surface->updateSkin();
                }
        }
}

void
COME_BioStructure::translate( const COME_Vector3D &disp ){
        
        translate( disp.getX(), disp.getY(), disp.getZ() );
}

void
COME_BioStructure::rotate( double rx, double ry, double rz ){

        getTissue()->rotate( rx, ry, rz );
        
        if( surface ){
                if( isFixed() ){
                        surface->rotate( rx, ry, rz );
                } else {
                        getTissue()->makeAllLocalFrames();
                        surface->updateSkin();
                }
        }
        
}

void
COME_BioStructure::scaleNominals( double factor ){

        if( !isFixed() ){
                getTissue()->scaleNominals( factor );
        }
}

void
COME_BioStructure::globalToLocal( COME_Matrix M, COME_Matrix MP ){

        if( getTissue() ){
                if( type == LIGAMENT ){
                        getTissue()->transform( M, MP );
                } else {
                        getTissue()->transform( M );
                }
        } else {
                surface->transform( M );
        }
}

void
COME_BioStructure::initializeHashStructure(){

        clock_t start,end;
        start = clock ();
        
        sUnit = COME_Vector3D( 0.0, 0.0, 1.0 ); 
        sSpherical = COME_Point2D( atan2( sUnit.y, sUnit.x ), acos( sUnit.z ) );
        hashSize = 20;
        
        if( hashTable ){
                delete [] hashTable;
        } else {
                hashTable = new COME_HashElement[hashSize*hashSize]; //malloc((unsigned) (nrow)*sizeof(double*))) == NULL)
        }
        
        for( int i = 0; i < hashSize; i++ ){
                for( int j = 0; j < hashSize; j++ ){
                        hashTable[elem(i,j)].face = NULL;
                        hashTable[elem(i,j)].distance = INFINITE;
                }
        }
        
        for( int iF = 0; iF < surface->getNumberFaces(); iF++ ){
                
                COME_Vector3D centerFace; centerFace.setVector3D( surface->getAFacePt(iF)->getCenterGlobalPosition() );
                centerFace.vpNormalize();
                COME_Point2D fSpherical = COME_Point2D( atan2( centerFace.y, centerFace.x ), acos( centerFace.z ) ) - sSpherical;
                int indexX = (int) ( ( fSpherical.getX() * ( ((double)(hashSize-1)) / ( 2.0 * M_PI ) ) ) + ( ((double)(hashSize-1)) / 2.0 ) );
                int indexY = (int) ( ( fSpherical.getY() * ( ((double)(hashSize-1)) / ( M_PI ) ) ) );//  + ( (double)(hashSize-1) / 2.0 ) );
                double distance = surface->getAFacePt(iF)->distancePointFaceGlobal( COME_Point3D() );
                //printf( "i: %d j: %d index: %d \n", indexX,indexY,  elem(indexX,indexY) );
                if( distance < hashTable[elem(indexX,indexY)].distance ){
                                if( ( surface->getAFacePt(iF)->getNormalGlobalPosition().vpDotProduct( centerFace ) < 0.0 ) ) { 
                                hashTable[elem(indexX,indexY)].face = surface->getAFacePt(iF);
                                hashTable[elem(indexX,indexY)].distance = distance;
                        }
                }
        }
        
        printf( "Hash completed with %d neighbor elements in 1st pass.\n", fillHashHoles() );
        printf( "Hash completed with %d neighbor elements in 2nd pass.\n", fillHashHoles() );
        
        for( int i = 0; i < hashSize; i++ ){
                for( int j = 0; j < hashSize; j++ ){
                        if( hashTable[elem(i,j)].distance > (INFINITE / 10.0) ) { 
                                
                                if( hashTable[elem(i,j)].face ){
                                        double cCollided[] = {0.1, 1.0, 0.1, 1.0};
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(0)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(1)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(2)->setColor( cCollided );
                                }
                                printf( "ZERO " );
                        } else{
                                if( hashTable[elem(i,j)].face ){
                                        double cCollided[] = {1.0, 1.0, 1.0, 1.0};
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(0)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(1)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(2)->setColor( cCollided );
                                }
                                printf( "%1.3f ", hashTable[elem(i,j)].distance );
                        }
                }
                printf( "\n" );
        }
        
        
        end = clock ();
        double dif = end - start;
        printf ("Init hash for collision time was %.8f seconds.\n", dif / (double) CLOCKS_PER_SEC );

}

void
COME_BioStructure::initializeHashStructureRayCasting(){

        clock_t start,end;
        start = clock ();
        
        sUnit = COME_Vector3D( 0.0, 0.0, 1.0 ); 
        sSpherical = COME_Point2D( atan2( sUnit.y, sUnit.x ), acos( sUnit.z ) );
        hashSize = 40;
        
        if( hashTable ){
                delete [] hashTable;
        } else {
                hashTable = new COME_HashElement[hashSize*hashSize]; //malloc((unsigned) (nrow)*sizeof(double*))) == NULL)
        }
        
        for( int i = 0; i < hashSize; i++ ){
                for( int j = 0; j < hashSize; j++ ){
                        
                        hashTable[elem(i,j)].face = NULL;
                        hashTable[elem(i,j)].distance = INFINITE;
                }
        }
        
        FILE * pFile;
        pFile = fopen ("saida.txt","w");
        
        for( int i = 0; i < hashSize; i++ ){
                for( int j = 0; j < hashSize; j++ ){
                        
                        // create ray
                        double theta = ( ( 2.0 * M_PI * i ) / ((double)( hashSize-1) ) ) - M_PI;
                        double phi = ( ( M_PI * j ) / ((double)( hashSize-1) ) );
                        
                        COME_Vector3D intersectRay( ( 1.0 * cos( theta ) * sin( phi ) ), ( 1.0 * sin( theta ) * sin( phi ) ), ( 1.0 * cos( phi) ) );
                        
                        //fprintf( pFile, "i: %d j: %d theta: %f phi: %f x: %f y: %f z: %f\n", i, j, theta, phi, intersectRay.x, intersectRay.y, intersectRay.z );
                        //printf( "i: %d j: %d theta: %f phi: %f x: %f y: %f z: %f\n", i, j, theta, phi, intersectRay.x, intersectRay.y, intersectRay.z );
                        
                        // for each face turned to center test intersection with ray
                        for( int iF = 0; iF < surface->getNumberFaces(); iF++ ){
                
                                COME_Vector3D centerFace( surface->getAFacePt(iF)->getVertexGlobalPositionPt( 0 )->x, surface->getAFacePt(iF)->getVertexGlobalPositionPt( 0 )->y, surface->getAFacePt(iF)->getVertexGlobalPositionPt( 0 )->z );
                                centerFace.vpNormalize();
                                if( ( surface->getAFacePt(iF)->getNormalGlobalPosition().vpDotProduct( centerFace ) < 0.0 ) ) { 
                                        COME_Point3D intersectPoint;
                                        if( surface->getAFacePt(iF)->intersectionRayTriangle( intersectRay, intersectPoint ) ) { 
                                                // put intersection onto the table with distance
                                                hashTable[elem(i,j)].face = surface->getAFacePt(iF);
                                                hashTable[elem(i,j)].distance = surface->getAFacePt(iF)->distancePointFaceGlobal( COME_Point3D() );
                                                break;
                                        }
                                }
                        }
                }
        }
        
        fclose (pFile);
        
        for( int i = 0; i < hashSize; i++ ){
                for( int j = 0; j < hashSize; j++ ){
                        if( hashTable[elem(i,j)].distance > (INFINITE / 10.0) ) { 
                                
                                if( hashTable[elem(i,j)].face ){
                                        double cCollided[] = {0.1, 1.0, 0.1, 1.0};
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(0)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(1)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(2)->setColor( cCollided );
                                }
                                printf( "ZERO " );
                        } else{
                                if( hashTable[elem(i,j)].face ){
                                        double cCollided[] = {1.0, 1.0, 1.0, 1.0};
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(0)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(1)->setColor( cCollided );
                                        hashTable[elem(i,j)].face->getVertexGlobalPositionPt(2)->setColor( cCollided );
                                }
                                printf( "%1.3f ", hashTable[elem(i,j)].distance );
                        }
                }
                printf( "\n" );
        }
        
        end = clock ();
        double dif = end - start;
        printf ("Init hash for collision time was %.8f seconds.\n", dif / (double) CLOCKS_PER_SEC );

}

COME_HashElement
COME_BioStructure::getHashElement( COME_Point2D p ){

        int indexX = (int) ( ( p.getX() * ( ((double)(hashSize-1)) / ( 2.0 * M_PI ) ) ) + ( ((double)(hashSize-1)) / 2.0 ) );
        int indexY = (int) ( ( p.getY() * ( ((double)(hashSize-1)) / (  M_PI ) ) ) );// + ( (double)(hashSize-1) / 2.0 ) );

        return hashTable[elem(indexX,indexY)];
}


int
COME_BioStructure::fillHashHoles(){

        int countAdditions = 0;
        COME_HashElement rowElement, columnElement;
        
        for( int i = 1; i < hashSize-1; i++ ){
                for( int j = 1; j < hashSize-1; j++ ){
                
                        rowElement.distance = columnElement.distance = INFINITE;
                        rowElement.face = columnElement.face = NULL;
        
                        if( !(hashTable[elem(i,j)].face) ){
                                for( int rowS = i-1; rowS >= 0; rowS-- ){
                                        if( hashTable[elem(rowS,j)].face ){
                                                for( int rowF = i+1; rowF < hashSize; rowF++ ){
                                                        if( hashTable[elem(rowF,j)].face ){
                                                                rowElement = ( hashTable[elem(rowS,j)].distance < hashTable[elem(rowF,j)].distance ) ? hashTable[elem(rowS,j)] : hashTable[elem(rowF,j)];
                                                                break;
                                                        }
                                                }
                                                break;
                                        }
                                }
                                
                                for( int colS = j-1; colS >= 0; colS-- ){
                                        if( hashTable[elem(i,colS)].face ){
                                                for( int colF = j+1; colF < hashSize; colF++ ){
                                                        if( hashTable[elem(i,colF)].face ){
                                                                columnElement = ( hashTable[elem(i,colS)].distance < hashTable[elem(i,colF)].distance ) ? hashTable[elem(i,colS)] : hashTable[elem(i,colF)];
                                                                break;
                                                        }
                                                }
                                                break;
                                        }
                                }
                                
                                if( rowElement.face || columnElement.face ) countAdditions++;
                                hashTable[elem(i,j)] = ( rowElement.distance < columnElement.distance ) ? rowElement : columnElement;
                        }
                }
        }
        return countAdditions;
}

---