diff --git a/contrib/arc/Classes/OctreeLSImage.cpp b/contrib/arc/Classes/OctreeLSImage.cpp
index e4cd02f92ac1ac9dad6b7d221bfa8f4eff3ec046..f6839e02ae7e4a03a4f59f463f17d08e83302b5b 100644
--- a/contrib/arc/Classes/OctreeLSImage.cpp
+++ b/contrib/arc/Classes/OctreeLSImage.cpp
@@ -1,5 +1,5 @@
//
-// Description:
+// Description:
//
//
// Author: <Boris Sedji>, 12/2009
@@ -13,36 +13,70 @@
#include "MElement.h"
OctreeLSImage::OctreeLSImage(itk::Image< float, 3 >::Pointer image){
-
+
itk::Image< float, 3 >::RegionType region;
// Define a region to get size of the image
region = image->GetLargestPossibleRegion ();
-
+
// Size will be a power of two larger than the image size
-
- _size[0] = 1;
- _size[1] = 1;
- _size[2] = 1;
- while (_size[0]<region.GetSize(0)){
- _size[0] = _size[0]*2;
+
+
+ int size[3];
+ size[0] = 1;
+ size[1] = 1;
+ size[2] = 1;
+
+ _ImageSize[0] = region.GetSize(0);
+ _ImageSize[1] = region.GetSize(1);
+ _ImageSize[2] = region.GetSize(2);
+
+ while (size[0]<_ImageSize[0]){
+ size[0] = size[0]*2;
}
- while (_size[1]<region.GetSize(1)){
- _size[1] = _size[1]*2;
+ while (size[1]<_ImageSize[1]){
+ size[1] = size[1]*2;
}
- while (_size[2]<region.GetSize(2)){
- _size[2] = _size[2]*2;
+ while (size[2]<_ImageSize[2]){
+ size[2] = size[2]*2;
}
+
+
BoxData *data = new BoxData;
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=_size[0];
- data->boundingbox[2]=0;
- data->boundingbox[3]=_size[1];
- data->boundingbox[4]=0;
- data->boundingbox[5]=_size[2];
+ if (size[0]>= size[1] & size[0]>=size[2] )
+ {
+ // The window of the image root
+ data->boundingbox[0]=0;
+ data->boundingbox[1]=size[0];
+ data->boundingbox[2]=0;
+ data->boundingbox[3]=size[0];
+ data->boundingbox[4]=0;
+ data->boundingbox[5]=size[0];
+ }
+ if (size[1]>= size[0] & size[1]>=size[2] )
+ {
+ // The window of the image root
+ data->boundingbox[0]=0;
+ data->boundingbox[1]=size[1];
+ data->boundingbox[2]=0;
+ data->boundingbox[3]=size[1];
+ data->boundingbox[4]=0;
+ data->boundingbox[5]=size[1];
+ }
+
+ if (size[2]>= size[1] & size[2]>=size[0] )
+ {
+ // The window of the image root
+ data->boundingbox[0]=0;
+ data->boundingbox[1]=size[2];
+ data->boundingbox[2]=0;
+ data->boundingbox[3]=size[2];
+ data->boundingbox[4]=0;
+ data->boundingbox[5]=size[2];
+ }
+
// initialisation of the octree
_root = new Box(data,NULL,this);
@@ -50,6 +84,8 @@ OctreeLSImage::OctreeLSImage(itk::Image< float, 3 >::Pointer image){
_image = image;
_LeafNumber = 0;
+ std::cout<<"size root :"<<_root->BoxSize()<<"\n";
+
}
@@ -57,7 +93,7 @@ void OctreeLSImage::Mesh(int maxsize,int minsize){
// taillemax //nombre de pixels max dans une direction
-
+
if (this->_root->HaveChildren()){
return;
}
@@ -96,14 +132,14 @@ void Box::FillLevelSetValue(itk::Image< float, 3 >::Pointer image, BoxData *data
for (int i=0;i<2;i++){
for (int j=0;j<2;j++){
for (int k=0;k<2;k++){
-
+
if (data->boundingbox[i] < region.GetSize(0)) pixelIndex[0] = data->boundingbox[i];
else pixelIndex[0] = region.GetSize(0)-1; // x position
if (data->boundingbox[2+j] < region.GetSize(1)) pixelIndex[1] = data->boundingbox[2+j];
else pixelIndex[1] = region.GetSize(1)-1; // y position
- if (data->boundingbox[4+k] < region.GetSize(2)) pixelIndex[2] = data->boundingbox[4+k];
+ if (data->boundingbox[4+k] < region.GetSize(2)) pixelIndex[2] = data->boundingbox[4+k];
else pixelIndex[2] = region.GetSize(2)-1; // z position
pixelValue = image->GetPixel( pixelIndex );
@@ -116,7 +152,7 @@ void Box::FillLevelSetValue(itk::Image< float, 3 >::Pointer image, BoxData *data
}
int Box::BoxSize(){
-
+
int SizeX;
int SizeY;
int SizeZ;
@@ -158,7 +194,7 @@ bool Box::IsHomogeneous(){
// Recursive dividing function - eight nodes
void Box::Mesh(int & maxsize, int & minsize){
-
+
if(((this->BoxSize()<=maxsize) & (this->IsHomogeneous())) | (this->BoxSize()<=minsize)){ // (max size & homogenous) ou (min size) => dont divide
return; // dont divide
@@ -187,11 +223,17 @@ void Box::PrintLevelSetValue(){
-GModel *OctreeLSImage::CreateGModel(){
-
- this->FillMeshInfo();
-
+GModel *OctreeLSImage::CreateGModel(bool simplex, double facx, double facy, double facz,int & sizemax,int & sizemin){
+
+ this->FillMeshInfo(sizemin);
+
+ itk::Image< float, 3 >::RegionType region;
+ // Define a region to get size of the image
+ region = _image->GetLargestPossibleRegion ();
+
+
std::map<int, MVertex*> vertexMap;
+ std::map<int,int> vertexNum;
std::vector<MVertex*> vertexVector;
std::vector<int> numElement;
@@ -200,57 +242,522 @@ GModel *OctreeLSImage::CreateGModel(){
std::vector<int> physical;
std::vector<int> elementary;
std::vector<int> partition;
+ std::vector<double> d;
+ std::map<int, std::vector<double> > data;
+ data.clear();
+
+ double eps = 1e-6;
int numVertices;
- numVertices = _ListNodes.size();
- int minVertex = numVertices+1;
- int maxVertex = -1;
- std::vector<std::vector<int> >::iterator it;
- it = _ListNodes.begin();
+
+ std::vector<std::vector<int> >::iterator ite;
+
int i = 1;
- while(it!=_ListNodes.end()){
- int num;
- float xyz[3];
- xyz[0] = ((*it)[0])*1;
- xyz[1] = ((*it)[1])*1.2;
- xyz[2] = ((*it)[2])*1.5;
- num = i;
- MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
- vertexMap[num] = newVertex;
- it++;
- i++;
- }
- minVertex = 1;
- maxVertex = numVertices;
+ int k = 1;
+
+ int a;
+
+ int xlimit = _ImageSize[0] - _ImageSize[0]%sizemax; // lost of image pixels...
+ int ylimit = _ImageSize[1] - _ImageSize[1]%sizemax;
+ int zlimit = _ImageSize[2] - _ImageSize[2]%sizemax;
+
+ std::cout<<"\nxlimit :"<< xlimit<<"\n";
+ std::cout<<"ylimit :"<< ylimit<<"\n";
+ std::cout<<"zlimit :"<< zlimit<<"\n";
+ std::cout<<"Listnodes size :"<<_ListNodes.size()<<"\n";
+
i = 1;
- it = _ListElements.begin();
- while(it!=_ListElements.end()) {
- int num, type, phys = 0, ele = 0, part = 0, numVertices;
- num = i;
- type = 5;
- ele = 26;
- phys = 0;
- part = 0;
- numVertices = MElement::getInfoMSH(type);
- std::vector <int> indices;
- for(int j = 0; j < numVertices; j++){
- indices.push_back((*it)[j]);
- }
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- indices.clear();
- it++;
- i++;
- }
-
+ ite = _ListElements.begin();
+ if (!simplex) // first try, to be improved
+ {
+ while(ite!=_ListElements.end()) {
+ int num, type, phys = 0, ele = 0, part = 0, numVertices;
+ num = i;
+ type = 5;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numVertices = MElement::getInfoMSH(type);
+ std::vector <int> indices;
+ for(int j = 0; j < numVertices; j++){
+ indices.push_back((*ite)[j]);
+ }
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ indices.clear();
+ ite++;
+ i++;
+ }
+ }
+ else // if simplex
+ {
+ while(ite!=_ListElements.end()) {
+
+ bool inImage;
+ inImage = true;
+ for (int j = 0;j<8;j++){
+// std::cout<<"_[(*ite)[j]]"<< (*ite)[j]<<"\n";
+// std::cout<<"_ListNodes[(*ite)[j]][0]"<< _ListNodes[(*ite)[j]-1][0] <<"\n";
+
+ if ( _ListNodes[(*ite)[j]-1][0] > xlimit | _ListNodes[(*ite)[j]-1][1] > ylimit | _ListNodes[(*ite)[j]-1][2] > zlimit ){
+ inImage = false;
+ }
+ }
+
+ if (inImage)
+ {
+
+ k++;
+ for (int j = 0;j<8;j++)
+ {
+ if (vertexMap.find((*ite)[j]) == vertexMap.end()) // if not in
+ {
+ float xyz[3];
+ d.clear();
+ xyz[0] = (_ListNodes[(*ite)[j]-1][0])*facx;
+ xyz[1] = (_ListNodes[(*ite)[j]-1][1])*facy;
+ xyz[2] = (_ListNodes[(*ite)[j]-1][2])*facz;
+ MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, (*ite)[j]);
+ vertexMap[(*ite)[j]] = newVertex;
+ d.push_back(_ListLSValue[(*ite)[j]-1]);
+ data[(*ite)[j]]=d;
+ d.clear();
+ }
+ }
+
+ int num, type, phys = 0, ele = 0, part = 0, numVertices;
+ num = i;
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numVertices = 4; //
+ std::vector <int> indices;
+ std::vector <int> front_1_indices;
+ std::vector <int> front_2_indices;
+ std::vector <int> front_3_indices;
+ std::vector <int> front_4_indices;
+ std::vector <int> front_5_indices;
+ std::vector <int> front_6_indices;
+
+ // --- Domain element tetraedre 1 ---
+ indices.clear();
+ indices.push_back((*ite)[3]);
+ indices.push_back((*ite)[6]);
+ indices.push_back((*ite)[7]);
+ indices.push_back((*ite)[4]);
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+
+
+ // --- frontier element---
+ front_1_indices.clear();
+ front_2_indices.clear();
+ front_3_indices.clear();
+ front_4_indices.clear();
+ front_5_indices.clear();
+ front_6_indices.clear();
+ for (int j =0;j<4;j++)
+ {
+ if (vertexMap[indices[j]]->z() > (zlimit-eps)*facz)
+ {
+ front_6_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
+ {
+ front_4_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->y() < eps*facy)
+ {
+ front_2_indices.push_back(indices[j]);
+ }
+ }
+
+ if (front_2_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 2;
+ phys = 2;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_2_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+ if (front_4_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 4;
+ phys = 4;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_4_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+ if (front_6_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 6;
+ phys = 6;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_6_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+ // --- Domain element tetraedre 2 ---
+ indices.clear();
+ indices.push_back((*ite)[6]);
+ indices.push_back((*ite)[3]);
+ indices.push_back((*ite)[1]);
+ indices.push_back((*ite)[4]);
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+
+
+ // --- Domain element tetraedre 3 ---
+ indices.clear();
+ indices.push_back((*ite)[3]);
+ indices.push_back((*ite)[4]);
+ indices.push_back((*ite)[0]);
+ indices.push_back((*ite)[1]);
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+
+
+ // --- frontier element tetraedre 3---
+ front_1_indices.clear();
+ front_2_indices.clear();
+ front_3_indices.clear();
+ front_4_indices.clear();
+ front_5_indices.clear();
+ front_6_indices.clear();
+ for (int j =0;j<4;j++)
+ {
+ if (vertexMap[indices[j]]->x() < eps*facx)
+ {
+ front_1_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->y() < eps*facy)
+ {
+ front_2_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->z() < eps*facz)
+ {
+ front_3_indices.push_back(indices[j]);
+ }
+ }
+ if (front_1_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 1;
+ phys = 1;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_1_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_2_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 2;
+ phys = 2;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_2_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_3_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 3;
+ phys = 3;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_3_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ // --- Domain element tetraedre 4 ---
+ indices.clear();
+ indices.push_back((*ite)[6]);
+ indices.push_back((*ite)[3]);
+ indices.push_back((*ite)[2]);
+ indices.push_back((*ite)[1]);
+ type = 4;
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+
+ // --- frontier element tetraedre 4---
+ front_1_indices.clear();
+ front_2_indices.clear();
+ front_3_indices.clear();
+ front_4_indices.clear();
+ front_5_indices.clear();
+ front_6_indices.clear();
+ for (int j =0;j<4;j++)
+ {
+ if (vertexMap[indices[j]]->x() < eps*facx)
+ {
+ front_1_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->y() > (ylimit-eps)*facy)
+ {
+ front_5_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->z() > (zlimit-eps)*facz)
+ {
+ front_6_indices.push_back(indices[j]);
+ }
+ }
+ if (front_1_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 1;
+ phys = 1;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_1_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_5_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 5;
+ phys = 5;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_5_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_6_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 6;
+ phys = 6;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_6_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+ // --- Domain element tetraedre 5 ---
+ indices.clear();
+ indices.push_back((*ite)[4]);
+ indices.push_back((*ite)[1]);
+ indices.push_back((*ite)[6]);
+ indices.push_back((*ite)[5]);
+ type = 4;
+ ele = 7;
+ phys = 7;
+ part = 0;
+ numElement.push_back(i);
+ vertexIndices.push_back(indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+
+ // --- frontier element tetraedre 4---
+ front_1_indices.clear();
+ front_2_indices.clear();
+ front_3_indices.clear();
+ front_4_indices.clear();
+ front_5_indices.clear();
+ front_6_indices.clear();
+ for (int j =0;j<4;j++)
+ {
+ if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
+ {
+ front_4_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->y() > (ylimit-eps)*facy)
+ {
+ front_5_indices.push_back(indices[j]);
+ }
+ if (vertexMap[indices[j]]->z() < eps*facz)
+ {
+ front_3_indices.push_back(indices[j]);
+ }
+ }
+ if (front_4_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 4;
+ phys = 4;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_4_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_5_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 5;
+ phys = 5;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_5_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ if (front_3_indices.size() == 3)
+ {
+ type = 2; // MSH_TRI
+ ele = 3;
+ phys = 3;
+ part = 0;
+ numVertices = 3;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_3_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+ }
+ ite++; // next octree leaf
+ }
+ std::cout<<"k :"<<k<<"\n";
+ }
+
+ std::cout<<"numElement size :"<<numElement.size()<<"\n";
+ std::cout<<"vertexIndices size :"<<vertexIndices.size()<<"\n";
+
+
+ std::cout<<"data size :"<<data.size()<<"\n";
+ std::cout<<"Vertexmap size :"<<vertexMap.size()<<"\n";
GModel *gmod = GModel::createGModel(vertexMap,numElement,vertexIndices,elementType,physical,elementary,partition);
-
+ // Write a .msh file with the level set values as postview data
+ std::string postTypeName = "LevelSet Value" ;
+ PView *pv = new PView (postTypeName, "NodeData", gmod, data);
+ //PView *pv = octree.CreateLSPView(m);
+ bool useadapt = true;
+ pv->getData(useadapt)->writeMSH("LSPView.msh", false);
+
+ std::cout<<"getNumScalars :"<< pv->getData(useadapt)->getNumScalars()<<"\n";
+
+//
+// std::string LevelSetFileName = "LevelSetValue.msh" ;
+// FILE *f = fopen(LevelSetFileName.c_str(), "w");
+// char name[256] = "LSValues";
+// int num;
+// float val;
+// int numnodes;
+// fprintf(f, "%s\n", name);
+// fprintf(f, "%d\n", data.size());
+// i = 1;
+// while(data.find(i)!=data.end()){
+// fprintf(f,"%d %f\n",i,data[i][0]);
+// i++;
+// }
+// fclose(f);
+
return gmod;
-
+
}
PView *OctreeLSImage::CreateLSPView(GModel* m){
@@ -259,7 +766,7 @@ PView *OctreeLSImage::CreateLSPView(GModel* m){
std::vector<float>::iterator itf;
itf = _ListLSValue.begin();
int i = 1;
-
+
while (itf!=_ListLSValue.end()){
std::vector<double> d;
d.push_back((*itf)) ;
@@ -271,17 +778,17 @@ PView *OctreeLSImage::CreateLSPView(GModel* m){
std::string postTypeName = "LevelSet Value" ;
PView *pv = new PView (postTypeName, "NodeData", m, data, 0.0);
- return pv;
+ return pv;
}
-void OctreeLSImage::FillMeshInfo(){
+void OctreeLSImage::FillMeshInfo(int & pas){
if (!_ListNodes.empty()){
return;
}
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue);
+ _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
}
@@ -290,35 +797,34 @@ void OctreeLSImage::FillMeshInfo(){
// fonction par balayage de région
-void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue){
+void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue,int &pas){
- int pas = 1;
-
- int* size = _octree->GetSize();
- int xpas = size[0]/size[2]*pas; // Ã modifier si sizeZ n'est pas le plus petit
- int ypas = size[1]/size[2]*pas;
- int zpas = size[2]/size[2]*pas;
+ int size = _octree->GetRoot()->BoxSize();
+
+ int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
+ int ypas = pas;
+ int zpas = pas;
std::cout<<"\nxyzpas :"<< xpas << " " << ypas << " " << zpas;
-
+
std::vector<std::vector<int> >::iterator it;
it = ListNodes.begin();
- int iti=0;
+ int iti=0;
- for(int z = 0;z<=size[2];z+=zpas)
- for(int y = 0;y<=size[1];y+=ypas)
- for(int x = 0;x<=size[0]; x+=xpas)
- {
+ for(int z = 0;z<=size;z+=zpas)
+ for(int y = 0;y<=size;y+=ypas)
+ for(int x = 0;x<=size; x+=xpas)
+ {
std::vector<Box*> Leafs;
Leafs.clear();
GetLeafsWith(Leafs, x, y, z);
- std::map<int,std::vector< int > > ElementsNodes;
+ std::map<int,std::vector< int > > ElementsNodes;
std::vector<int> NodesIn;
bool added = false;
- for(unsigned int l=0;l<Leafs.size();l++){
+ for(unsigned int l=0;l<Leafs.size();l++){
if (!added){
std::vector<int> XYZ;
XYZ.push_back(x);
@@ -334,49 +840,67 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector<
iti++;
}
}
- }else break;
+ }else break;
}
if (added){
- for(unsigned int l=0;l<Leafs.size();l++){
+ for(unsigned int l=0;l<Leafs.size();l++){
for (int i = 0 ; i<2;i++ )
for (int k = 0 ; k<2;k++ )
for (int j = 0 ; j<2;j++ ){
int pos ;
- pos = i*4+k*2+j*1;
+ pos = i*4+k*2+j*1; // pour l'autre méthode
if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2] & z == Leafs[l]->GetData()->boundingbox[k+4]){
- Leafs[l]->SetElementNode(pos,iti);
+ Leafs[l]->SetElementNode(pos,iti); //différence avec l'autre méthode
}
- }
+ }
}
}
}
std::cout<<"\nNumber Nodes : "<< ListNodes.size();
_octree->GetRoot()->FillElementsNode(ListElements);
-
+
}
void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
-
- // if it s a leaf then fill list node if node dont exist in list
+ // if it s a leaf then fill list node if node dont exist in list
if (!this->HaveChildren()) {
std::vector<int> ElementNodes;
for (int i = 0 ; i < 8 ; i++) ElementNodes.push_back(_ElementNodes[i]);
+ int temp;
+ temp = ElementNodes[2];
+ ElementNodes[2] = ElementNodes[3];
+ ElementNodes[3] = temp;
+ temp = ElementNodes[6];
+ ElementNodes[6] = ElementNodes[7];
+ ElementNodes[7] = temp;
+
+// temp = ElementNodes[1];
+// ElementNodes[1] = ElementNodes[2];
+// ElementNodes[2] = temp;
+// temp = ElementNodes[2];
+// ElementNodes[2] = ElementNodes[3];
+// ElementNodes[3] = temp;
+//
+// temp = ElementNodes[5];
+// ElementNodes[5] = ElementNodes[6];
+// ElementNodes[6] = temp;
+// temp = ElementNodes[6];
+// ElementNodes[6] = ElementNodes[7];
+// ElementNodes[7] = temp;
+
ListElements.push_back(ElementNodes);
- std::cout<<"\n Nodes :";
- for (int i = 0 ; i < 8 ; i++)
- std::cout<<" " << ElementNodes[i];
+// std::cout<<"\n Nodes :";
+// for (int i = 0 ; i < 8 ; i++)
+// std::cout<<" " << ElementNodes[i];
return;
}
for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
+ this->_children[n]->FillElementsNode(ListElements);
}
-
}
-
-
/*
void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue){
@@ -427,7 +951,7 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector<
}
for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillMeshInfo(ListNodes,ListElements,ListLSValue);
+ this->_children[n]->FillMeshInfo(ListNodes,ListElements,ListLSValue);
}
}*/
@@ -452,7 +976,7 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
// if (!this->HaveChildren()) {
// std::vector<std::vector<int> >::iterator it;
// std::vector<std::vector<int> > *ListNodes;
-// ListNodes = this->_octree->GetListNodes();
+// ListNodes = this->_octree->GetListNodes();
// it = ListNodes->begin();
// std::vector<int> NodesIn;
// int l = 1;
@@ -471,9 +995,9 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
// ListElements.push_back(NodesIn);
// return;
// }
-//
+//
// for (int n = 0 ; n < 8 ; n++){
-// this->_children[n]->GetLeafElements(ListElements);
+// this->_children[n]->GetLeafElements(ListElements);
// }
}
@@ -484,7 +1008,7 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
void OctreeLSImage::SetLeafNumber(){
_LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
+ _root->CountLeafNumber(_LeafNumber);
}
void Box::CountLeafNumber(int &LeafNumber){
@@ -495,7 +1019,7 @@ void Box::CountLeafNumber(int &LeafNumber){
}
for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
+ this->_children[n]->CountLeafNumber(LeafNumber);
}
}
@@ -507,11 +1031,11 @@ void Box::Divide(){
// box divide in 8 boxes whit limite xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2 ; zl = (zmin + zmax)/2
int n;
-
- // Box 1 : < xl ; < yl ; < zl
- n = 0;
- BoxData* data1 = new BoxData;
+ // Box 1 : < xl ; < yl ; < zl
+
+ n = 0;
+ BoxData* data1 = new BoxData;
data1->boundingbox[0]=this->_data->boundingbox[0];
data1->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
@@ -524,9 +1048,9 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 2 : > xl ; < yl ; < zl
-
- n = 1;
- BoxData* data2 = new BoxData;
+
+ n = 1;
+ BoxData* data2 = new BoxData;
data2->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data2->boundingbox[1]=this->_data->boundingbox[1];
@@ -539,10 +1063,10 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 3 : < xl ; > yl ; < zl
-
+
n = 2;
- BoxData* data3 = new BoxData;
-
+ BoxData* data3 = new BoxData;
+
data3->boundingbox[0]=this->_data->boundingbox[0];
data3->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data3->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
@@ -554,10 +1078,10 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 4 : > xl ; > yl ; < zl
-
- n = 3;
- BoxData* data4 = new BoxData;
-
+
+ n = 3;
+ BoxData* data4 = new BoxData;
+
data4->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data4->boundingbox[1]=this->_data->boundingbox[1];
data4->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
@@ -569,10 +1093,10 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 5 : < xl ; < yl ; > zl
-
+
n = 4;
- BoxData* data5 = new BoxData;
-
+ BoxData* data5 = new BoxData;
+
data5->boundingbox[0]=this->_data->boundingbox[0];
data5->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data5->boundingbox[2]=this->_data->boundingbox[2];
@@ -584,10 +1108,10 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 6 : > xl ; < yl ; > zl
-
+
n = 5;
- BoxData* data6 = new BoxData;
-
+ BoxData* data6 = new BoxData;
+
data6->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data6->boundingbox[1]=this->_data->boundingbox[1];
data6->boundingbox[2]=this->_data->boundingbox[2];
@@ -597,12 +1121,12 @@ void Box::Divide(){
this->_children[n] = new Box(data6,this,this->_octree);
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
+
// Box 7 : < xl ; > yl ; > zl
-
+
n = 6;
- BoxData* data7 = new BoxData;
-
+ BoxData* data7 = new BoxData;
+
data7->boundingbox[0]=this->_data->boundingbox[0];
data7->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data7->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
@@ -614,9 +1138,9 @@ void Box::Divide(){
this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
// Box 8 : > xl ; > yl ; > zl
-
- n = 7;
- BoxData* data8 = new BoxData;
+
+ n = 7;
+ BoxData* data8 = new BoxData;
data8->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
data8->boundingbox[1]=this->_data->boundingbox[1];
@@ -638,7 +1162,7 @@ bool OctreeLSImage::Smooth(){
bool Box::Smooth(){
-
+
if (!this->HaveChildren()) {
bool smoothed = true;
BoxData*data = this->GetData();
@@ -650,11 +1174,11 @@ bool Box::Smooth(){
_octree->GetRoot()->GetLeafsWith(Leafs,data->boundingbox[i],data->boundingbox[2+j],data->boundingbox[4+k]);
if (Leafs.size()!=8){
int min = Leafs[0]->BoxSize();
- for(unsigned int l=0;l<Leafs.size();l++){
+ for(unsigned int l=0;l<Leafs.size();l++){
if(min>Leafs[l]->BoxSize()) min = Leafs[l]->BoxSize();
- }
- for(unsigned int l=0;l<Leafs.size();l++){
- if(min<=(Leafs[l]->BoxSize())/4){
+ }
+ for(unsigned int l=0;l<Leafs.size();l++){
+ if(min<=(Leafs[l]->BoxSize())/4){
Leafs[l]->Divide();
smoothed = false;
}
@@ -668,7 +1192,7 @@ bool Box::Smooth(){
bool smoothed = true;
for (int n = 0 ; n < 8 ; n++){
- if(!this->_children[n]->Smooth()) smoothed = false;
+ if(!this->_children[n]->Smooth()) smoothed = false;
}
return smoothed;
}
@@ -684,15 +1208,15 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z){
data = parent->GetParent()->GetData();
parent=parent->GetParent();
}
- Leafs.push_back(this);
- return;
+ Leafs.push_back(this);
+ return;
}
-
+
if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[0]->GetLeafsWith(Leafs,x,y,z);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[1]->GetLeafsWith(Leafs,x,y,z);
}
@@ -700,7 +1224,7 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z){
if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[2]->GetLeafsWith(Leafs,x,y,z);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[3]->GetLeafsWith(Leafs,x,y,z);
}
@@ -708,7 +1232,7 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z){
if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[4]->GetLeafsWith(Leafs,x,y,z);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[5]->GetLeafsWith(Leafs,x,y,z);
}
@@ -716,9 +1240,9 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z){
if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[6]->GetLeafsWith(Leafs,x,y,z);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
_children[7]->GetLeafsWith(Leafs,x,y,z);
}
-
+
}
diff --git a/contrib/arc/Classes/OctreeLSImage.h b/contrib/arc/Classes/OctreeLSImage.h
index c4acd5e08e69f89f5c60a48809b5d7a905cbec0b..3b732115f69c59912e58e82eb14608b65091774e 100644
--- a/contrib/arc/Classes/OctreeLSImage.h
+++ b/contrib/arc/Classes/OctreeLSImage.h
@@ -1,5 +1,5 @@
//
-// Description:
+// Description:
//
//
// Author: <Boris Sedji>, 12/2009
@@ -34,15 +34,15 @@ class Box;
class OctreeLSImage
{
-
+
private :
// Octree's root
Box* _root;
// Octree's image
itk::Image< float, 3 >::Pointer _image;
- // Window size
- int _size[3];
+ // Region size
+ int _ImageSize[3];
// Mesh list of nodes
std::vector< std::vector <int> > _ListNodes;
// Mesh list of elements
@@ -53,7 +53,7 @@ class OctreeLSImage
int _LeafNumber;
// Count octree's leafs
void SetLeafNumber();
- void FillMeshInfo();
+ void FillMeshInfo(int & pas);
public :
@@ -65,14 +65,14 @@ class OctreeLSImage
void Mesh(int maxsize, int minsize);
itk::Image< float, 3 >::Pointer GetImage(){return _image;}
int GetLeafNumber(){this->SetLeafNumber();return _LeafNumber;}
- // Refine if too much generations between adjacent leafs
+ // Refine if too much generations between adjacent leafs
bool Smooth();
//std::vector< std::vector <int> > *GetListElements();
// Create GModel
- GModel* CreateGModel();
+ GModel* CreateGModel(bool simplex, double facx, double facy, double facz,int & sizemax,int & sizemin);
// Create PView representation of the level set
PView* CreateLSPView(GModel* m);
- int* GetSize(){return _size;}
+ int* GetSize(){return _ImageSize;}
};
@@ -93,7 +93,7 @@ class Box{
BoxData* GetData(){return _data;}
// Does it have children
bool HaveChildren();
- // Fill level set value in data with image values
+ // Fill level set value in data with image values
void FillLevelSetValue(itk::Image< float, 3 >::Pointer image, BoxData *data);
// The smallest length of the box
int BoxSize();
@@ -113,10 +113,10 @@ class Box{
// Give all the leafs containing this node
void GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z);
// Recursive function to create the list of the mesh nodes, elements and levelset values
- void FillMeshInfo(std::vector<std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector< float > &ListLSValue);
+ void FillMeshInfo(std::vector<std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector< float > &ListLSValue,int &pas);
// Recursive function to create the list of mesh elements in relation with the list of nodes
void SetElementNode(int pos,int iti){_ElementNodes[pos]=iti;}
void FillElementsNode(std::vector< std::vector<int> > &ListElements);
};
-#endif
\ No newline at end of file
+#endif