Skip to content
Snippets Groups Projects
Commit 9a5280a1 authored by Boris Sedji's avatar Boris Sedji
Browse files

No commit message 

No commit message
parent 58fbabbf
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
contrib/arc/Classes/QuadtreeLSImage.cpp
+ 469
− 136
View file @ 9a5280a1
...@@ -21,23 +21,42 @@ QuadtreeLSImage::QuadtreeLSImage(itk::Image< float, 2 >::Pointer image){ ...@@ -21,23 +21,42 @@ QuadtreeLSImage::QuadtreeLSImage(itk::Image< float, 2 >::Pointer image){
// Size will be a power of two larger than the image size // Size will be a power of two larger than the image size
_size[0] = 1; int size[3];
_size[1] = 1; size[0] = 1;
size[1] = 1;
while (_size[0]<region.GetSize(0)){ _ImageSize[0] = region.GetSize(0);
_size[0] = _size[0]*2; _ImageSize[1] = region.GetSize(1);
while (size[0]<region.GetSize(0)){
size[0] = size[0]*2;
} }
while (_size[1]<region.GetSize(1)){ while (size[1]<region.GetSize(1)){
_size[1] = _size[1]*2; size[1] = size[1]*2;
} }
BoxData *data = new BoxData; BoxData *data = new BoxData;
// The window of the image root // The window of the image root
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[0]=0;
data->boundingbox[1]=_size[0]; data->boundingbox[1]=size[1];
data->boundingbox[2]=0; data->boundingbox[2]=0;
data->boundingbox[3]=_size[1]; data->boundingbox[3]=size[1];
}
// initialisation of the quadtree // initialisation of the quadtree
_root = new Box(data,NULL,this); _root = new Box(data,NULL,this);
...@@ -168,10 +187,10 @@ void Box::PrintLevelSetValue(){ ...@@ -168,10 +187,10 @@ void Box::PrintLevelSetValue(){
} }
GModel *QuadtreeLSImage::CreateGModel(bool simplex, double facx, double facy,int & sizemax,int & sizemin){
GModel *QuadtreeLSImage::CreateGModel(){ double eps = 1e-6;
this->FillMeshInfo(sizemin);
this->FillMeshInfo();
std::map<int, MVertex*> vertexMap; std::map<int, MVertex*> vertexMap;
std::vector<MVertex*> vertexVector; std::vector<MVertex*> vertexVector;
...@@ -182,41 +201,43 @@ GModel *QuadtreeLSImage::CreateGModel(){ ...@@ -182,41 +201,43 @@ GModel *QuadtreeLSImage::CreateGModel(){
std::vector<int> physical; std::vector<int> physical;
std::vector<int> elementary; std::vector<int> elementary;
std::vector<int> partition; std::vector<int> partition;
std::vector<double> d;
std::map<int, std::vector<double> > data;
data.clear();
int numVertices; int numVertices;
numVertices = _ListNodes.size();
int minVertex = numVertices+1; std::vector<std::vector<int> >::iterator ite;
int maxVertex = -1;
std::vector<std::vector<int> >::iterator it;
it = _ListNodes.begin();
int i = 1; int i = 1;
while(it!=_ListNodes.end()){ int k = 1;
int num;
float xyz[3]; int a;
xyz[0] = ((*it)[0])*1;
xyz[1] = ((*it)[1])*1.2; int xlimit = _ImageSize[0] - _ImageSize[0]%sizemax; // lost of image pixels...
xyz[2] = 0; int ylimit = _ImageSize[1] - _ImageSize[1]%sizemax;
num = i;
MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num); std::cout<<"\nxlimit :"<< xlimit<<"\n";
vertexMap[num] = newVertex; std::cout<<"ylimit :"<< ylimit<<"\n";
it++;
i++; std::cout<<"Listnodes size :"<<_ListNodes.size()<<"\n";
}
minVertex = 1;
maxVertex = numVertices;
i = 1; i = 1;
it = _ListElements.begin(); ite = _ListElements.begin();
while(it!=_ListElements.end()) {
if (!simplex) // first try, to be improved
{
while(ite!=_ListElements.end()) {
int num, type, phys = 0, ele = 0, part = 0, numVertices; int num, type, phys = 0, ele = 0, part = 0, numVertices;
num = i; num = i;
type = 3; // MSH_QUA_4 type = 3; // MSH_QUA_4
ele = 26; ele = 1;
phys = 0; phys = 1;
part = 0; part = 0;
numVertices = MElement::getInfoMSH(type); numVertices = MElement::getInfoMSH(type);
std::vector <int> indices; std::vector <int> indices;
for(int j = 0; j < numVertices; j++){ for(int j = 0; j < numVertices; j++){
indices.push_back((*it)[j]); indices.push_back((*ite)[j]);
} }
numElement.push_back(i); numElement.push_back(i);
vertexIndices.push_back(indices); vertexIndices.push_back(indices);
...@@ -225,11 +246,324 @@ GModel *QuadtreeLSImage::CreateGModel(){ ...@@ -225,11 +246,324 @@ GModel *QuadtreeLSImage::CreateGModel(){
elementary.push_back(ele); elementary.push_back(ele);
partition.push_back(part); partition.push_back(part);
indices.clear(); indices.clear();
it++; 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++; 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); 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; return gmod;
...@@ -257,13 +591,13 @@ PView *QuadtreeLSImage::CreateLSPView(GModel* m){ ...@@ -257,13 +591,13 @@ PView *QuadtreeLSImage::CreateLSPView(GModel* m){
} }
void QuadtreeLSImage::FillMeshInfo(){ void QuadtreeLSImage::FillMeshInfo(int & pas){
if (!_ListNodes.empty()){ if (!_ListNodes.empty()){
return; return;
} }
_root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue); _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
} }
...@@ -272,14 +606,13 @@ void QuadtreeLSImage::FillMeshInfo(){ ...@@ -272,14 +606,13 @@ void QuadtreeLSImage::FillMeshInfo(){
// fonction par balayage de région // 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 size = _quadtree->GetRoot()->BoxSize();
int xpas = size[0]/size[1]*pas; int xpas =pas; // à modifier si sizeZ n'est pas le plus petit
int ypas = size[1]/size[1]*pas; // à modifier si sizeY n'est pas le plus petit int ypas = pas;
std::cout<<"\nxyzpas :"<< xpas << " " << ypas ; std::cout<<"\nxyzpas :"<< xpas << " " << ypas ;
...@@ -290,8 +623,8 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< ...@@ -290,8 +623,8 @@ void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector<
int iti=0; int iti=0;
for(int y = 0;y<=size[1];y+=ypas) for(int y = 0;y<=size;y+=ypas)
for(int x = 0;x<=size[0]; x+=xpas) for(int x = 0;x<=size; x+=xpas)
{ {
std::vector<Box*> Leafs; std::vector<Box*> Leafs;
Leafs.clear(); Leafs.clear();
...@@ -347,9 +680,9 @@ void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){ ...@@ -347,9 +680,9 @@ void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
_ElementNodes[2] = temp; _ElementNodes[2] = temp;
for (int i = 0 ; i < 4 ; i++) ElementNodes.push_back(_ElementNodes[i]); for (int i = 0 ; i < 4 ; i++) ElementNodes.push_back(_ElementNodes[i]);
ListElements.push_back(ElementNodes); ListElements.push_back(ElementNodes);
std::cout<<"\n Nodes :"; // std::cout<<"\n Nodes :";
for (int i = 0 ; i < 4 ; i++) // for (int i = 0 ; i < 4 ; i++)
std::cout<<" " << ElementNodes[i]; // std::cout<<" " << ElementNodes[i];
return; return;
} }
......
contrib/arc/Classes/QuadtreeLSImage.h
+ 11
− 11
View file @ 9a5280a1
...@@ -41,8 +41,8 @@ class QuadtreeLSImage ...@@ -41,8 +41,8 @@ class QuadtreeLSImage
Box* _root; Box* _root;
// Quadtree's image // Quadtree's image
itk::Image< float, 2 >::Pointer _image; itk::Image< float, 2 >::Pointer _image;
// Window size // Region size
int _size[2]; int _ImageSize[2];
// Mesh list of nodes // Mesh list of nodes
std::vector< std::vector <int> > _ListNodes; std::vector< std::vector <int> > _ListNodes;
// Mesh list of elements // Mesh list of elements
...@@ -53,7 +53,7 @@ class QuadtreeLSImage ...@@ -53,7 +53,7 @@ class QuadtreeLSImage
int _LeafNumber; int _LeafNumber;
// Count octree's leafs // Count octree's leafs
void SetLeafNumber(); void SetLeafNumber();
void FillMeshInfo(); void FillMeshInfo(int & pas);
public : public :
...@@ -69,10 +69,10 @@ class QuadtreeLSImage ...@@ -69,10 +69,10 @@ class QuadtreeLSImage
bool Smooth(); bool Smooth();
//std::vector< std::vector <int> > *GetListElements(); //std::vector< std::vector <int> > *GetListElements();
// Create GModel // Create GModel
GModel* CreateGModel(); GModel* CreateGModel(bool simplex, double facx, double facy,int & sizemax,int & sizemin);
// Create PView representation of the level set // Create PView representation of the level set
PView* CreateLSPView(GModel* m); PView* CreateLSPView(GModel* m);
int* GetSize(){return _size;} int* GetSize(){return _ImageSize;}
}; };
...@@ -113,7 +113,7 @@ class Box{ ...@@ -113,7 +113,7 @@ class Box{
// Give all the leafs containing this node // Give all the leafs containing this node
void GetLeafsWith(std::vector<Box*> &Leafs, int x, int y); void GetLeafsWith(std::vector<Box*> &Leafs, int x, int y);
// Recursive function to create the list of the mesh nodes, elements and levelset values // 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 // 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 SetElementNode(int pos,int iti){_ElementNodes[pos]=iti;}
void FillElementsNode(std::vector< std::vector<int> > &ListElements); void FillElementsNode(std::vector< std::vector<int> > &ListElements);
......
contrib/arc/ImportLS2dImage.cpp
+ 17
− 14
View file @ 9a5280a1
...@@ -66,17 +66,20 @@ int main( int argc, char *argv[] ) ...@@ -66,17 +66,20 @@ int main( int argc, char *argv[] )
k++; k++;
} }
bool simplex = 1;
double facx = 1;
double facy = 1;
// Create GModel with the octree mesh // Create GModel with the octree mesh
GModel *m = quadtree.CreateGModel(); GModel *m = quadtree.CreateGModel(simplex,facx,facy,sizemax, sizemin);
// Write a .msh file with the mesh // Write a .msh file with the mesh
std::string ModelName = "QuadtreeMesh.msh" ; 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 // Write a .msh file with the level set values as postview data
PView *pv = quadtree.CreateLSPView(m); // PView *pv = quadtree.CreateLSPView(m);
bool useadapt = true; // bool useadapt = true;
pv->getData(useadapt)->writeMSH("LSPView.msh", false); // pv->getData(useadapt)->writeMSH("LSPView.msh", false);
std::cout<<"\n"; std::cout<<"\n";
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment