Commit 8a443162 authored by aderam's avatar aderam
Browse files

NEW: MML distance to triangle surface

NEW: brain32 and face examples 
FIXED: Sofa Simulator now works if indexes of atoms do not follow each other

git-svn-id: svn+ssh://scm.forge.imag.fr/var/lib/gforge/chroot/scmrepos/svn/camitk/trunk/camitk@199 ec899d31-69d1-42ba-9299-647d76f65fb3
parent 1ade3307
......@@ -234,6 +234,7 @@ set(monitoring_H
monitor/monitors/MonitorSurface.h
monitor/monitors/MonitorDisplacement.h
monitor/monitors/MonitorNormDisplacement.h
monitor/monitors/MonitorPointToTriangleMeshDistanceFinal.h
reference/Reference.h
tools/Chrono.h
tools/Tools.h
......@@ -287,6 +288,7 @@ set(monitoring_SRCS
monitor/monitors/MonitorSurface.cpp
monitor/monitors/MonitorDisplacement.cpp
monitor/monitors/MonitorNormDisplacement.cpp
monitor/monitors/MonitorPointToTriangleMeshDistanceFinal.cpp
reference/Reference.cpp
tools/Chrono.cpp
tools/Tools.cpp
......@@ -360,6 +362,7 @@ export_headers(
tools/Chrono.h
tools/Tools.h
tools/Macros.h
tools/AtomIterator.h
tools/SurfaceExtractor/Facet.h
tools/SurfaceExtractor/SurfaceExtractor.h
SUBDIRECTORY tools
......@@ -386,6 +389,7 @@ export_headers(
monitor/monitors/MonitorForce.h
monitor/monitors/MonitorPointSetDistance.h
monitor/monitors/MonitorPointFinalSetDistance.h
monitor/monitors/MonitorPointToTriangleMeshDistanceFinal.h
SUBDIRECTORY monitor/monitors
COMPONENT monitoring
)
......
......@@ -407,6 +407,7 @@ void MonitoringManager::simulate(){
init();
doMove();
while(!checkStop()){
std::cout << currentTime << std::endl;
doMove();
}
saveMonitors();
......
......@@ -47,6 +47,7 @@ Monitor* MonitorFactory::createMonitor(mml::Monitor* m,MonitoringManager* monito
case mml::MonitorType::PointFinalSetDistance: return (new MonitorPointFinalSetDistance(m,monitoringManager)); break;
case mml::MonitorType::Volume: return (new MonitorVolume(m,monitoringManager)); break;
case mml::MonitorType::Surface: return (new MonitorSurface(m,monitoringManager)); break;
case mml::MonitorType::DistanceToTriangularMeshFinal: return (new MonitorPointToTriangleMeshDistanceFinal(m,monitoringManager)); break;
default: std::cout << "Monitor type error" << std::endl;
}
return NULL;
......
......@@ -38,6 +38,7 @@
#include "monitor/monitors/MonitorSurface.h"
#include "monitor/monitors/MonitorDisplacement.h"
#include "monitor/monitors/MonitorNormDisplacement.h"
#include "monitor/monitors/MonitorPointToTriangleMeshDistanceFinal.h"
/**
* A factory to create monitors
......
/*****************************************************************************
$CAMITK_LICENCE_BEGIN$
CamiTK - Computer Assisted Medical Intervention ToolKit
Visit http://camitk.imag.fr for more information
Copyright (C) 2012 Celine Fouard, Emmanuel Promayon, Yannick Keraval
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301 USA
$CAMITK_LICENCE_END$
*****************************************************************************/
#include "MonitorPointToTriangleMeshDistanceFinal.h"
#include "tools/Tools.h"
#include "tools/AtomIterator.h"
// -------------------- constructor --------------------
MonitorPointToTriangleMeshDistanceFinal::MonitorPointToTriangleMeshDistanceFinal(mml::Monitor* m,MonitoringManager* monitoringManager): Monitor(m,monitoringManager,SCALARSET){}
// -------------------- destructor --------------------
MonitorPointToTriangleMeshDistanceFinal::~MonitorPointToTriangleMeshDistanceFinal(){}
// -------------------- destructor --------------------
void MonitorPointToTriangleMeshDistanceFinal::calculate()
{
values.clear();
double posSimul[3];
double posRef[3];
double realTime;
double dist=0;
//TODO indexes are in index vector: do not iterate and use indexes vector directly?
AtomIterator it=AtomIterator(monitoringManager->getPml(),target);
for (it.begin();!it.end();it.next()){
it.currentAtom()->getPosition(posSimul);
posSimul[0]=posSimul[0]+dx;
posSimul[1]=posSimul[1]+dy;
posSimul[2]=posSimul[2]+dz;
if (references[0]->getDistanceToTriangularMesh(posSimul,dist)){
values.push_back(dist);
}else{
cerr << "no 'Elements' component found of no triangles in 'Elements' " << endl;
}
}
write();
}
// -------------------- destructor --------------------
std::string MonitorPointToTriangleMeshDistanceFinal::getTypeName()
{
return "Point to triangle mesh distance";
}
/*****************************************************************************
$CAMITK_LICENCE_BEGIN$
CamiTK - Computer Assisted Medical Intervention ToolKit
Visit http://camitk.imag.fr for more information
Copyright (C) 2012 Celine Fouard, Emmanuel Promayon, Yannick Keraval
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301 USA
$CAMITK_LICENCE_END$
*****************************************************************************/
#ifndef MONITOR_MONITORS_MONITORPOINTTOTRIANGLEMESHDISTANCEFINAL_H
#define MONITOR_MONITORS_MONITORPOINTTOTRIANGLEMESHDISTANCEFINAL_H
#include <string>
#include "MonitorIn.hxx"
#include "monitor/Monitor.h"
/**
* A monitor that calculate the distance between a point a triangle mes
* This is an optimized version of the point to triangle mesh distance, that only calculate distance with the final reference state
*/
class MonitorPointToTriangleMeshDistanceFinal: public Monitor{
public:
/**
* constructor
* @param m the xsdcxx generated monitor
*/
MonitorPointToTriangleMeshDistanceFinal(mml::Monitor* m,MonitoringManager* monitoringManager);
/// destructor
~MonitorPointToTriangleMeshDistanceFinal();
/// calculate current followed data and store them in values vector
void calculate();
/// return a string relative to monitor type
std::string getTypeName();
};
#endif // MONITOR_MONITORS_MONITORPOINTTOTRIANGLEMESHDISTANCEFINAL_H
\ No newline at end of file
......@@ -28,7 +28,7 @@
#include <monitor/MonitorFactory.h>
#include "tools/AtomIterator.h"
#include "tools/Tools.h"
#include "pml/CellProperties.h"
#include <string.h>
......@@ -212,6 +212,70 @@ bool Reference::getMonitoredData(string type, double time, int index, double &re
return false;
}
// -------------------- getDistanceToTriangularMesh --------------------
bool Reference::getDistanceToTriangularMesh(double pos[3], double& dist)
{
//look for nearest atom
Atom* a;
double temp[3];
dist=-1;
double realTime;
int i=0;
AtomIterator it=AtomIterator(pml,target);
for (it.begin();!it.end();it.next()){
it.currentAtom()->getPosition(temp);
double d=distance(pos,temp);
if (i==0 || d<dist){
dist=d;
a=it.currentAtom();
i++;
}
}
//heuristic? look cell containing tthe nearest atom and compute distance to triangle
std::vector<StructuralComponent*> comps=a->getAllStructuralComponents();
for (int i=0;i<comps.size();i++){
StructuralComponent* str=comps[i];
Cell* c = dynamic_cast<Cell*>(str);
if (c) {
if (c->getProperties()->getType()==StructureProperties::TRIANGLE){
double v1[3];
double v2[3];
double v3[3];
((Atom*)(c->getStructure(0)))->getPosition(v1);
((Atom*)(c->getStructure(1)))->getPosition(v2);
((Atom*)(c->getStructure(2)))->getPosition(v3);
double d=distanceToTriangle(pos,v1,v2,v3);
if (d<dist)
dist=d;
}
}
}
return dist!=-1;
// StructuralComponent* str=(StructuralComponent*)(pml->getComponentByName(target));
// if (!str)
// return false;
// dist=0;
// bool found=false;
// for (unsigned int i=0; i< str->getNumberOfCells();i++){
// Cell* c=str->getCell(i);
// if (c->getProperties()->getType()==StructureProperties::TRIANGLE){
// double v1[3];
// double v2[3];
// double v3[3];
// ((Atom*)(c->getStructure(0)))->getPosition(v1);
// ((Atom*)(c->getStructure(1)))->getPosition(v2);
// ((Atom*)(c->getStructure(2)))->getPosition(v3);
// double d=distanceToTriangle(pos,v1,v2,v3);
// if (!found || d<dist){
// dist=d;
// found=true;
// }
// }
// }
// return found;
}
// -------------------- toString --------------------
std::string Reference::toString()
......
......@@ -63,6 +63,15 @@ class Reference{
*/
bool getNearest(double pos[3], double ref[3]);
/**
* get distance to a triagular mesh, the target of the reference msut contain triangles
* this is an optimized method to get distance to triangular mesh where there is only one final time step (real reference) using final pml
* @param pos atom's position (with eventual offset)
* @param dist the distance
* @return true if method succeded (a triangular mesh was found)
*/
bool getDistanceToTriangularMesh(double pos[3],double& dist);
/**
* get the the values of a given monitor which do not depend of time or an atom (e.g. geometrical data)
* @param type name of the monitor, must be a existing monitor in the mmlOut document
......
......@@ -42,10 +42,12 @@ std::string ParametersWriter::write(){
std::ostringstream os;
//iterate over the fields of the exclusive components and add them as parameters in batch file
Properties* prop=monitoringManager->getPml()->getExclusiveComponents()->getProperties();
for(unsigned int i=0;i<prop->numberOfFields();i++){
std::string field=prop->getField(i);
os << field << "=" << prop->getString(field) << std::endl;
}
if (prop){
for(unsigned int i=0;i<prop->numberOfFields();i++){
std::string field=prop->getField(i);
os << field << "=" << prop->getString(field) << std::endl;
}
return os.str();
}
}
\ No newline at end of file
......@@ -183,6 +183,8 @@ SofaSimulator::SofaSimulator(MonitoringManager* monitoringManager): InteractiveS
break;
case 4:
//tetra
scn << " <TetrahedronSetTopologyContainer name=\"topo\" />" << endl;
scn << " <TetrahedronSetGeometryAlgorithms template=\"Vec3d\" name=\"GeomAlgo\" />" << endl;
scn<< " <TetrahedralCorotationalFEMForceField template=\"Vec3d\" name=\"FEM\" method=\"large\" poissonRatio=\"0.45\" youngModulus=\"1440\" />" << endl;
break;
default:
......@@ -194,23 +196,25 @@ SofaSimulator::SofaSimulator(MonitoringManager* monitoringManager): InteractiveS
}
build();
// create the atom / DOF map
unsigned int atomId = 0; // TODO: Warning! the actual PML should be used to make
unsigned int id=0; // TODO: Warning! the actual PML should be used to make
// perfect correspondance with the DOF index (for e.g. using two matching positions)
StructuralComponent* sc=monitoringManager->getPml()->getAtoms();
for (unsigned int i = 0; i < mechanicalObjects.size(); i++) {
mechanicalObjectAtomDOFMap.push_back(new std::MechanicalObjectAtomDOFMap());
mechanicalObjectDOFAtomMap.push_back(new std::MechanicalObjectDOFAtomMap());
for (unsigned int j = 0; j < getMechanicalObjectDOFPosition(i).size(); j++) {
// insert in the global map
unsigned int atomId=((Atom*)(sc->getStructure(id)))->getIndex();
atomsToDOF.insert(std::pair<unsigned int, std::MechanicalObjectDOFIndex> (atomId, std::MechanicalObjectDOFIndex(i, j)));
// insert in the local map
mechanicalObjectAtomDOFMap.back()->insert(std::pair<unsigned int, unsigned int>(atomId, j));
mechanicalObjectDOFAtomMap.back()->insert(std::pair<unsigned int, unsigned int>(j, atomId));
// next id
atomId++;
id++;
}
}
......
......@@ -27,15 +27,44 @@
#include "tools/Tools.h"
#include <iostream>
#include "math.h"
// -------------------- distance --------------------
double distance(double pos[3], double pos2[3]){
return sqrt( (pos[0] - pos2[0])*(pos[0] - pos2[0])
+ (pos[1] - pos2[1])*(pos[1] - pos2[1])
+ (pos[2] - pos2[2])*(pos[2] - pos2[2]));
// -------------------- distanceToTriangle --------------------
double distanceToTriangle(double point[3], double tri1[3], double tri2[3], double tri3[3])
{
double norm[3];
double v1_2[3];
double v1_3[3];
double v1_0[3];
for (int i=0;i<3;i++){
v1_2[i]=tri2[i]-tri1[i];
v1_3[i]=tri3[i]-tri1[i];
v1_0[i]=point[i]-tri1[i];
}
crossProduct(v1_2,v1_3,norm);
//distance to projection on triangle, it is the seeked distance if the projection is inside the triangle
double dist=distance(point,tri1)*dotProduct(v1_0,norm)/(normOf(v1_0)*normOf(norm));
if (dist<0)
dist=-dist;
//approximation here, return min between dist calculated and distance with nodes of the triangle
//TODO find true distance (i.e. if the minimum diatance is between the point and one edge of the triangle)
/*double min=dist;
double dists[3];
dists[0]=distance(point,tri1);
dists[1]=distance(point,tri2);
dists[2]=distance(point,tri3);
for(int i=0;i<3;i++){
if (dists[i]>min)
min=dists[i];
}
return min;*/
return dist;
}
// -------------------- timeParameter2double --------------------
double timeParameter2double(mml::TimeParameter& t){
switch (t.unit()){
......
......@@ -27,11 +27,69 @@
#define TOOLS_TOOLS_H
#include "MonitoringModel.hxx"
#include "math.h"
/// compute euclidian distance
double distance(double pos[3], double pos2[3]);
/// compute the distance of a point to a triangle in 3D
/// @param point coords of the point
/// @param tri1 coords of the 1st point of the triangle
/// @param tri2 coords of the 2nd point of the triangle
/// @param tri3 coords of the 3rd point of the triangle
double distanceToTriangle(double point[3], double tri1[3], double tri2[3], double tri3[3]);
/// convert a TimeParameter (from xsd-cxx generetaed file) to double
double timeParameter2double(mml::TimeParameter& t);
/// compute dot product
double dotProduct(double vec1[3], double vec2[3]);
/// compute cross product
double crossProduct(double vec1[3], double vec2[3], double res[3]);
/// normalize vector
void normalize(double vec[3]);
/// norm of vector
double normOf(double vec[3]);
// -------------------- distance --------------------
inline double distance(double pos[3], double pos2[3]){
return sqrt( (pos[0] - pos2[0])*(pos[0] - pos2[0])
+ (pos[1] - pos2[1])*(pos[1] - pos2[1])
+ (pos[2] - pos2[2])*(pos[2] - pos2[2]));
}
// -------------------- dotProduct --------------------
inline double dotProduct(double vec1[3], double vec2[3])
{
return vec1[0]*vec2[0]+vec1[1]*vec2[1]+vec1[2]*vec2[2];
}
// -------------------- distance --------------------
inline double crossProduct(double vec1[3], double vec2[3], double res[3])
{
res[0]=vec1[1]*vec1[2]-vec1[2]*vec1[1];
res[1]=vec1[2]*vec1[0]-vec1[0]*vec1[2];
res[2]=vec1[0]*vec1[1]-vec1[1]*vec1[0];
}
// -------------------- normalize --------------------
inline void normalize(double vec[3])
{
double norm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
if (norm >0){
vec[0]=vec[0]/norm;
vec[1]=vec[1]/norm;
vec[2]=vec[2]/norm;
}
}
// -------------------- normOf --------------------
inline double normOf(double vec[3])
{
return sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
}
#endif // TOOLS_TOOLS_H
\ No newline at end of file
......@@ -238,6 +238,7 @@ $CAMITK_LICENCE_END$
<xsd:enumeration value="Surface"/>
<xsd:enumeration value="Displacement"/>
<xsd:enumeration value="NormDisplacement"/>
<xsd:enumeration value="DistanceToTriangularMeshFinal"/>
</xsd:restriction>
</xsd:simpleType>
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment