diff --git a/contrib/arc/Classes/QuadtreeLSImage.cpp b/contrib/arc/Classes/QuadtreeLSImage.cpp
index 70994c66755866e3297e596f1e3e12a9d1ffeff9..9a2c0e457e59634be3d6db788ea4c30ba9ffeb98 100644
--- a/contrib/arc/Classes/QuadtreeLSImage.cpp
+++ b/contrib/arc/Classes/QuadtreeLSImage.cpp
@@ -1,5 +1,5 @@
//
-// Description:
+// Description:
//
//
// Author: <Boris Sedji>, 12/2009
@@ -13,31 +13,50 @@
#include "MElement.h"
QuadtreeLSImage::QuadtreeLSImage(itk::Image< float, 2 >::Pointer image){
-
+
itk::Image< float, 2 >::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;
- while (_size[0]<region.GetSize(0)){
- _size[0] = _size[0]*2;
+ int size[3];
+ size[0] = 1;
+ size[1] = 1;
+
+ _ImageSize[0] = region.GetSize(0);
+ _ImageSize[1] = region.GetSize(1);
+
+ while (size[0]<region.GetSize(0)){
+ size[0] = size[0]*2;
}
- while (_size[1]<region.GetSize(1)){
- _size[1] = _size[1]*2;
+ while (size[1]<region.GetSize(1)){
+ size[1] = size[1]*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];
+
+ if (size[0]>= size[1])
+ {
+ // The window of the image root
+ data->boundingbox[0]=0;
+ data->boundingbox[1]=size[0];
+ data->boundingbox[2]=0;
+ data->boundingbox[3]=size[0];
+
+ }
+ if (size[1]>= size[0])
+ {
+ // The window of the image root
+ data->boundingbox[0]=0;
+ data->boundingbox[1]=size[1];
+ data->boundingbox[2]=0;
+ data->boundingbox[3]=size[1];
+
+ }
// initialisation of the quadtree
_root = new Box(data,NULL,this);
@@ -52,7 +71,7 @@ void QuadtreeLSImage::Mesh(int maxsize,int minsize){
// taillemax //nombre de pixels max dans une direction
-
+
if (this->_root->HaveChildren()){
return;
}
@@ -90,7 +109,7 @@ void Box::FillLevelSetValue(itk::Image< float, 2 >::Pointer image, BoxData *data
// fill with pixel value, and if out of region size, fill with value at largest region size
for (int i=0;i<2;i++){
for (int j=0;j<2;j++){
-
+
if (data->boundingbox[i] < region.GetSize(0)) pixelIndex[0] = data->boundingbox[i];
else pixelIndex[0] = region.GetSize(0)-1; // x position
@@ -106,7 +125,7 @@ void Box::FillLevelSetValue(itk::Image< float, 2 >::Pointer image, BoxData *data
}
int Box::BoxSize(){
-
+
int SizeX;
int SizeY;
int size;
@@ -142,7 +161,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
@@ -168,11 +187,11 @@ void Box::PrintLevelSetValue(){
}
+GModel *QuadtreeLSImage::CreateGModel(bool simplex, double facx, double facy,int & sizemax,int & sizemin){
+
+ double eps = 1e-6;
+ this->FillMeshInfo(sizemin);
-GModel *QuadtreeLSImage::CreateGModel(){
-
- this->FillMeshInfo();
-
std::map<int, MVertex*> vertexMap;
std::vector<MVertex*> vertexVector;
@@ -182,57 +201,372 @@ GModel *QuadtreeLSImage::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();
int numVertices;
- numVertices = _ListNodes.size();
- int minVertex = numVertices+1;
- int maxVertex = -1;
- std::vector<std::vector<int> >::iterator it;
- it = _ListNodes.begin();
- 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] = 0;
- num = i;
- MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
- vertexMap[num] = newVertex;
- it++;
- i++;
- }
- minVertex = 1;
- maxVertex = numVertices;
+
+ std::vector<std::vector<int> >::iterator ite;
+
+ int i = 1;
+ int k = 1;
+
+ int a;
+
+ int xlimit = _ImageSize[0] - _ImageSize[0]%sizemax; // lost of image pixels...
+ int ylimit = _ImageSize[1] - _ImageSize[1]%sizemax;
+
+ std::cout<<"\nxlimit :"<< xlimit<<"\n";
+ std::cout<<"ylimit :"<< ylimit<<"\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 = 3; // MSH_QUA_4
- 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 = 3; // MSH_QUA_4
+ ele = 1;
+ phys = 1;
+ 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<4;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 )
+ {
+ inImage = false;
+ }
+ }
+
+ //inImage = true;
+ if (inImage)
+ {
+
+ k++;
+ for (int j = 0;j<4;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] = 0;
+ 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;
+ numVertices = 4;
+ std::vector <int> indices;
+ std::vector<int> front_indices;
+
+ // --- domain elements - first triangle ---
+
+
+ indices.clear();
+ for(int j = 0; j < 3; j++){
+ indices.push_back((*ite)[j]);
+ }
+ type = 2; // MSH_TRI_3
+ ele = 5;
+ phys = 5;
+ 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 face 1 -> y = 0 ---
+
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->y() < eps*facy)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 1;
+ phys = 1;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ // --- frontier element face 2 -> x = 0 ---
+
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->x() < eps*facx)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 2;
+ phys = 2;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+
+ // - frontier element face 3 -> y = _size[1]
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->y() > (ylimit - eps)*facy)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 3;
+ phys = 3;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ // --- frontier element face 4 -> x = _size[0] ---
+
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 4;
+ phys = 4;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+
+ indices.clear();
+ // second triangle
+ for(int j = 2; j < 5; j++){
+ indices.push_back((*ite)[j%4]);
+ }
+ type = 2; // MSH_TRI_3
+ ele = 5;
+ phys = 5;
+ 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 face 1 -> y = 0 ---
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->y() < eps*facy)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 1;
+ phys = 1;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+ // --- frontier element face 2 -> x = 0 ---
+
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->x() < eps*facx)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 2;
+ phys = 2;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+
+ // --- frontier element face 3 -> y = _size[1] ---
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->y() > (ylimit - eps)*facy)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 3;
+ phys = 3;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+ front_indices.clear();
+
+
+ // --- frontier element face 4 -> x = _size[0] ---
+
+ front_indices.clear();
+ for (int j =0;j<3;j++)
+ {
+ if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
+ {
+ front_indices.push_back(indices[j]);
+ }
+ }
+ if (front_indices.size() == 2)
+ {
+ type = 1; // MSH_LINE
+ ele = 4;
+ phys = 4;
+ part = 0;
+ numVertices = 2;
+ numElement.push_back(i);
+ vertexIndices.push_back(front_indices);
+ elementType.push_back(type);
+ physical.push_back(phys);
+ elementary.push_back(ele);
+ partition.push_back(part);
+ i++;
+ }
+
+ }
+ ite++; // next quadtree leaf
+ }
}
-
+
+
+ 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";
+
return gmod;
-
+
}
PView *QuadtreeLSImage::CreateLSPView(GModel* m){
@@ -241,7 +575,7 @@ PView *QuadtreeLSImage::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)) ;
@@ -253,17 +587,17 @@ PView *QuadtreeLSImage::CreateLSPView(GModel* m){
std::string postTypeName = "LevelSet Value" ;
PView *pv = new PView (postTypeName, "NodeData", m, data, 0.0);
- return pv;
+ return pv;
}
-void QuadtreeLSImage::FillMeshInfo(){
+void QuadtreeLSImage::FillMeshInfo(int & pas){
if (!_ListNodes.empty()){
return;
}
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue);
+ _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
}
@@ -272,35 +606,34 @@ void QuadtreeLSImage::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 = _quadtree->GetSize();
- int xpas = size[0]/size[1]*pas;
- int ypas = size[1]/size[1]*pas; // Ã modifier si sizeY n'est pas le plus petit
+ int size = _quadtree->GetRoot()->BoxSize();
+
+ int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
+ int ypas = pas;
std::cout<<"\nxyzpas :"<< xpas << " " << ypas ;
-
+
std::vector<std::vector<int> >::iterator it;
it = ListNodes.begin();
- int iti=0;
+ int iti=0;
- for(int y = 0;y<=size[1];y+=ypas)
- for(int x = 0;x<=size[0]; x+=xpas)
- {
+ 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);
- 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);
@@ -308,18 +641,18 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector<
XYZ.push_back(0);
for (int i = 0 ; i<2;i++ )
for (int j = 0 ; j<2;j++ )
-
+
if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2] ){
ListLSValue.push_back(Leafs[l]->GetData()->LevelSetValue[i][j]);
ListNodes.push_back(XYZ);
added = true;
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 j = 0 ; j<2;j++ ){
int pos ;
@@ -327,13 +660,13 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector<
if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2]){
Leafs[l]->SetElementNode(pos,iti);
}
- }
+ }
}
}
}
std::cout<<"\nNumber Nodes : "<< ListNodes.size();
_quadtree->GetRoot()->FillElementsNode(ListElements);
-
+
}
void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
@@ -347,14 +680,14 @@ void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
_ElementNodes[2] = temp;
for (int i = 0 ; i < 4 ; i++) ElementNodes.push_back(_ElementNodes[i]);
ListElements.push_back(ElementNodes);
- std::cout<<"\n Nodes :";
- for (int i = 0 ; i < 4 ; i++)
- std::cout<<" " << ElementNodes[i];
+// std::cout<<"\n Nodes :";
+// for (int i = 0 ; i < 4 ; i++)
+// std::cout<<" " << ElementNodes[i];
return;
}
for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
+ this->_children[n]->FillElementsNode(ListElements);
}
}
@@ -412,7 +745,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);
}
}*/
@@ -437,7 +770,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->_quadtree->GetListNodes();
+// ListNodes = this->_quadtree->GetListNodes();
// it = ListNodes->begin();
// std::vector<int> NodesIn;
// int l = 1;
@@ -456,9 +789,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);
// }
}
@@ -469,7 +802,7 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
void QuadtreeLSImage::SetLeafNumber(){
_LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
+ _root->CountLeafNumber(_LeafNumber);
}
void Box::CountLeafNumber(int &LeafNumber){
@@ -480,7 +813,7 @@ void Box::CountLeafNumber(int &LeafNumber){
}
for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
+ this->_children[n]->CountLeafNumber(LeafNumber);
}
}
@@ -489,67 +822,67 @@ void Box::CountLeafNumber(int &LeafNumber){
void Box::Divide(){
-// box divide in 4 boxes with limits xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2
+// box divide in 4 boxes with limits xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2
int n;
-
- // Box 1 : < xl ; < yl
- n = 0;
- BoxData* data1 = new BoxData;
+ // Box 1 : < xl ; < yl
+
+ 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;
data1->boundingbox[2]=this->_data->boundingbox[2];
data1->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
-
+
this->_children[n] = new Box(data1,this,this->_quadtree);
this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 2 : > xl ; < yl
-
- n = 1;
- BoxData* data2 = new BoxData;
+ // Box 2 : > xl ; < yl
+
+ 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];
data2->boundingbox[2]=this->_data->boundingbox[2];
data2->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
-
+
this->_children[n] = new Box(data2,this,this->_quadtree);
this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 3 : < xl ; > yl
-
+ // Box 3 : < xl ; > yl
+
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;
data3->boundingbox[3]=this->_data->boundingbox[3];
-
+
this->_children[n] = new Box(data3,this,this->_quadtree);
this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 4 : > xl ; > yl
-
- n = 3;
- BoxData* data4 = new BoxData;
-
+ // Box 4 : > xl ; > yl
+
+ 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;
data4->boundingbox[3]=this->_data->boundingbox[3];
-
+
this->_children[n] = new Box(data4,this,this->_quadtree);
this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
-
+
}
bool QuadtreeLSImage::Smooth(){
@@ -561,7 +894,7 @@ bool QuadtreeLSImage::Smooth(){
bool Box::Smooth(){
-
+
if (!this->HaveChildren()) {
bool smoothed = true;
BoxData*data = this->GetData();
@@ -572,11 +905,11 @@ bool Box::Smooth(){
_quadtree->GetRoot()->GetLeafsWith(Leafs,data->boundingbox[i],data->boundingbox[2+j]);
if (Leafs.size()!=4){
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;
}
@@ -589,7 +922,7 @@ bool Box::Smooth(){
bool smoothed = true;
for (int n = 0 ; n < 4 ; n++){
- if(!this->_children[n]->Smooth()) smoothed = false;
+ if(!this->_children[n]->Smooth()) smoothed = false;
}
return smoothed;
}
@@ -605,15 +938,15 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y){
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) ){
_children[0]->GetLeafsWith(Leafs,x,y);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) ){
_children[1]->GetLeafsWith(Leafs,x,y);
}
@@ -621,10 +954,10 @@ void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y){
if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) ){
_children[2]->GetLeafsWith(Leafs,x,y);
}
-
+
if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2)){
_children[3]->GetLeafsWith(Leafs,x,y);
}
-
+
}
diff --git a/contrib/arc/Classes/QuadtreeLSImage.h b/contrib/arc/Classes/QuadtreeLSImage.h
index 1ab7495fd7e219334ee1ed5d0d57f24bcc326dd2..f956762c9fd2b0c86ce2cd21ad26bffd94e2a626 100644
--- a/contrib/arc/Classes/QuadtreeLSImage.h
+++ b/contrib/arc/Classes/QuadtreeLSImage.h
@@ -1,5 +1,5 @@
//
-// Description:
+// Description:
//
//
// Author: <Boris Sedji>, 12/2009
@@ -34,15 +34,15 @@ class Box;
class QuadtreeLSImage
{
-
+
private :
// Quadtree's root
Box* _root;
// Quadtree's image
itk::Image< float, 2 >::Pointer _image;
- // Window size
- int _size[2];
+ // Region size
+ int _ImageSize[2];
// Mesh list of nodes
std::vector< std::vector <int> > _ListNodes;
// Mesh list of elements
@@ -53,7 +53,7 @@ class QuadtreeLSImage
int _LeafNumber;
// Count octree's leafs
void SetLeafNumber();
- void FillMeshInfo();
+ void FillMeshInfo(int & pas);
public :
@@ -65,14 +65,14 @@ class QuadtreeLSImage
void Mesh(int maxsize, int minsize);
itk::Image< float, 2 >::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,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, 2 >::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);
// 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
diff --git a/contrib/arc/ImportLS2dImage.cpp b/contrib/arc/ImportLS2dImage.cpp
index 34cc9e82f6daaf25c108e7af5687672471ed545d..1ac56bda2caf203b5a378d384422d31ab4ccc6f9 100755
--- a/contrib/arc/ImportLS2dImage.cpp
+++ b/contrib/arc/ImportLS2dImage.cpp
@@ -35,19 +35,19 @@ int main( int argc, char *argv[] )
typedef itk::ImageFileReader< ImageTypeFloat > ImageReaderTypeFloat;
ImageReaderTypeFloat::Pointer reader = ImageReaderTypeFloat::New();
-
- reader->SetFileName(argv[1]);
+
+ reader->SetFileName(argv[1]);
ImageTypeFloat::Pointer image = reader->GetOutput();
image->Update();
ImageTypeFloat::RegionType region;
region = image->GetLargestPossibleRegion ();
-
+
std::cout<<"\nImage dimensions : " << region.GetSize(0) << " x " << region.GetSize(1);
QuadtreeLSImage quadtree(image);
-
+
int sizemax = atoi(argv[2]);
int sizemin = atoi(argv[3]);
@@ -65,19 +65,22 @@ int main( int argc, char *argv[] )
std::cout<<"\nLeaf Number : "<< (quadtree.GetLeafNumber())<<"\n";
k++;
}
-
+
+ bool simplex = 1;
+ double facx = 1;
+ double facy = 1;
// Create GModel with the octree mesh
- GModel *m = quadtree.CreateGModel();
-
- // Write a .msh file with the mesh
+ GModel *m = quadtree.CreateGModel(simplex,facx,facy,sizemax, sizemin);
+
+ // Write a .msh file with the mesh
std::string ModelName = "QuadtreeMesh.msh" ;
- m->writeMSH(ModelName,2.1,false,false);
-
+ m->writeMSH(ModelName,2.1,false,false);
+
// Write a .msh file with the level set values as postview data
- PView *pv = quadtree.CreateLSPView(m);
- bool useadapt = true;
- pv->getData(useadapt)->writeMSH("LSPView.msh", false);
-
+// PView *pv = quadtree.CreateLSPView(m);
+// bool useadapt = true;
+// pv->getData(useadapt)->writeMSH("LSPView.msh", false);
+
std::cout<<"\n";
return 0;