Commit 5b961c6c authored by Jan Eichler's avatar Jan Eichler
Browse files

working on issue 2

parent 33fde170
......@@ -23,15 +23,25 @@
#include "marray.hxx"
#include "seg.h"
std::string cAxes_version = "21.03.01 slim";
std::string cAxes_version = "21.04.15development slim";
class Grains
class Segment
{
public:
//Attributes *************************************************************
std::vector<point> pointer;
int phase; // background=0, ice=1
int size;
double xcenter;
double ycenter;
int child_index; // reference to the child, e.g. grain index
};
class Grains: public Segment //child object of the parent Segment
{
public:
//Attributes *************************************************************
int parent_index; // reference to the parent region
vigra::TinyVector<float, 3> axis;
float az;
float co;
......@@ -60,7 +70,7 @@ class Boundaries
float twist_to_boundary;
};
/*
/*vpn.awi.de
typedef struct
{
int size;
......@@ -1180,17 +1190,17 @@ int export_grains_bw(vigra::IImage grain_label_image, std::string path)
int export_grains_rgb(vigra::IImage grain_label_image, std::vector<Grains> &Grain, std::string path)
{
vigra::FRGBImage color_image(mWidth, mHeight);
vigra::FRGBImage color_image(mWidth, mHeight, black);
float azim,colat;
int x,y;
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
if (grain_label_image(x,y)==0)
color_image(x,y)=black;
else
for (int i=1; i<=nr_grains; i++)
for (int pixel_iterator=0; pixel_iterator<Grain[i].pointer.size(); pixel_iterator++)
{
azim=Grain[grain_label_image(x,y)].az;
colat=Grain[grain_label_image(x,y)].co;
azim=Grain[i].az;
colat=Grain[i].co;
x = Grain[i].pointer[pixel_iterator].x;
y = Grain[i].pointer[pixel_iterator].y;
if (azim>=0.0 && azim<120.0)
{
......@@ -1212,7 +1222,7 @@ int export_grains_rgb(vigra::IImage grain_label_image, std::vector<Grains> &Grai
}
color_image(x,y)+=pow(cos(colat*M_PI/180),5)*(white-color_image(x,y));
}
exportImage(srcImageRange(color_image), vigra::ImageExportInfo(path.c_str()));
return 0;
}
......@@ -1221,15 +1231,15 @@ int export_grains_trendw(vigra::IImage grain_label_image, std::vector<Grains> &G
{
vigra::FRGBImage color_image(mWidth, mHeight);
float azim,colat;
int x,y;
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
if (grain_label_image(x,y)==0)
color_image(x,y)=black;
else
for (int i=1; i<=nr_grains; i++)
for (int pixel_iterator=0; pixel_iterator<Grain[i].pointer.size(); pixel_iterator++)
{
azim=Grain[grain_label_image(x,y)].az;
colat=Grain[grain_label_image(x,y)].co;
azim=Grain[i].az;
colat=Grain[i].co;
x = Grain[i].pointer[pixel_iterator].x;
y = Grain[i].pointer[pixel_iterator].y;
if (azim>=0.0 && azim<30.0)
{
......@@ -1310,20 +1320,23 @@ int export_grains_trendw(vigra::IImage grain_label_image, std::vector<Grains> &G
return 0;
}
int export_boundary_rgb(vigra::IImage boundary_label_image, vigra::IImage regions, vigra::BImage bground_image, std::vector<Grains> &Grain, std::vector<Boundaries> &Boundary, std::string path)
int export_boundary_rgb(vigra::IImage boundary_label_image, vigra::IImage regions, vigra::BImage bground_image, std::vector<Segment> &Segments, std::vector<Grains> &Grain, std::vector<Boundaries> &Boundary, std::string path)
{
vigra::FRGBImage color_image(mWidth, mHeight);
vigra::FRGBImage color_image(mWidth, mHeight, black);
float azim,colat;
int grain_index, segment_index;
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
{
if (boundary_label_image(x,y)==0)
color_image(x,y)=black;
else
segment_index = boundary_label_image(x,y);
if (Segments[segment_index].phase == 1)
{
azim=Grain[boundary_label_image(x,y)].az;
colat=Grain[boundary_label_image(x,y)].co;
grain_index = Segments[segment_index].child_index;
azim=Grain[grain_index].az;
colat=Grain[grain_index].co;
if (azim>=0.0 && azim<120.0)
{
......@@ -1345,34 +1358,19 @@ int export_boundary_rgb(vigra::IImage boundary_label_image, vigra::IImage region
}
color_image(x,y)+=pow(cos(colat*M_PI/180),5)*(white-color_image(x,y));
}
if (bground_image(x,y)==0) color_image(x,y)=black; //border to background
}
int dd = 1; //grain boundary thickness in the image
int dd = 1; //grain boundary thickness in the image
for (int y=dd; y<mHeight-dd; y++)
for (int x=dd; x<mWidth-dd; x++)
{
if (regions(x,y)==0)
{
for (int dy=-dd; dy<=dd; dy++)
for (int dx=-dd; dx<=dd; dx++)
{
for (int dy=-dd; dy<=dd; dy++)
for (int dx=-dd; dx<=dd; dx++)
color_image(x+dx,y+dy)=black;
// color_image(x-dx,y+dy)=black;
// color_image(x-dx,y-dy)=black;
// color_image(x+dx,y-dy)=black;
// color_image(x+1,y)=black;
// color_image(x-1,y)=black;
// color_image(x,y+1)=black;
// color_image(x,y-1)=black;
// color_image(x+1,y+1)=black;
// color_image(x-1,y-1)=black;
// color_image(x+1,y-1)=black;
// color_image(x-1,y+1)=black;
}
}
}
......@@ -1380,20 +1378,22 @@ int export_boundary_rgb(vigra::IImage boundary_label_image, vigra::IImage region
return 0;
}
int export_boundary_trendw(vigra::IImage boundary_label_image, vigra::IImage regions, vigra::BImage bground_image, std::vector<Grains> &Grain, std::vector<Boundaries> &Boundary, std::string path)
int export_boundary_trendw(vigra::IImage boundary_label_image, vigra::IImage regions, vigra::BImage bground_image, std::vector<Segment> &Segments, std::vector<Grains> &Grain, std::vector<Boundaries> &Boundary, std::string path)
{
vigra::FRGBImage color_image(mWidth, mHeight);
float azim,colat;
int grain_index, segment_index;
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
{
if (boundary_label_image(x,y)==0)
color_image(x,y)=black;
else
{
segment_index = boundary_label_image(x,y);
if (Segments[segment_index].phase == 1)
{
azim=Grain[boundary_label_image(x,y)].az;
colat=Grain[boundary_label_image(x,y)].co;
grain_index = Segments[segment_index].child_index;
azim=Grain[grain_index].az;
colat=Grain[grain_index].co;
if (azim>=0.0 && azim<30.0)
{
......@@ -1472,40 +1472,25 @@ int export_boundary_trendw(vigra::IImage boundary_label_image, vigra::IImage reg
if (bground_image(x,y)==0) color_image(x,y)=black; //border to background
}
int dd = 1; //grain boundary thickness in the image
int dd = 1; //grain boundary thickness in the image
for (int y=dd; y<mHeight-dd; y++)
for (int x=dd; x<mWidth-dd; x++)
{
if (regions(x,y)==0)
{
for (int dy=-dd; dy<=dd; dy++)
for (int dx=-dd; dx<=dd; dx++)
{
for (int dy=-dd; dy<=dd; dy++)
for (int dx=-dd; dx<=dd; dx++)
color_image(x+dx,y+dy)=black;
// color_image(x-dx,y+dy)=black;
// color_image(x-dx,y-dy)=black;
// color_image(x+dx,y-dy)=black;
// color_image(x+1,y)=black;
// color_image(x-1,y)=black;
// color_image(x,y+1)=black;
// color_image(x,y-1)=black;
// color_image(x+1,y+1)=black;
// color_image(x-1,y-1)=black;
// color_image(x+1,y-1)=black;
// color_image(x-1,y+1)=black;
}
}
}
}
exportImage(srcImageRange(color_image), vigra::ImageExportInfo(path.c_str()));
return 0;
}
int export_regions(vigra::IImage regions, std::vector<Boundaries> &Boundary, std::string path)
int export_regions(vigra::IImage regions, std::vector<Boundaries> &Boundary, std::string path) //experiment
{
vigra::FRGBImage image(mWidth, mHeight);
int xx, yy;
......@@ -1965,27 +1950,110 @@ int create_topo(vigra::FVector3Image &vector_image, vigra::FImage &topo_image)
return 0;
}
int calculate_grains(vigra::IImage &grain_label_image, vigra::FVector3Image &vector_image, std::vector<Grains> &Grain)
int classify_segments(std::vector<Segment> &Segments, vigra::IImage &grain_label_image, vigra::BImage bground_image)
{
std::cout << "classify segments" << std::endl;
for(int i=1; i<=nr_regions; i++)
{
Segments[i].phase=0;
Segments[i].size=0;
Segments[i].xcenter=0;
Segments[i].ycenter=0;
Segments[i].pointer.resize(1);
Segments[i].pointer[0].x = Segments[i].pointer[0].y = 0;
}
//calculate center of mass positions, mean segment size
point point_xy;
for (int x=0; x<mWidth; x++)
for (int y=0; y<mHeight; y++)
{
//segment labeling starts with 1
if (grain_label_image(x,y)>0)
{
int i = grain_label_image(x,y);
if (bground_image(x,y)==1)
Segments[i].phase=1;
Segments[i].size++;
Segments[i].xcenter+=x/float(mWidth);
Segments[i].ycenter+=y/float(mHeight);
point_xy.x = x;
point_xy.y = y;
Segments[i].pointer.push_back(point_xy);
}
}
for(int i=1; i<=nr_regions; i++)
{
Segments[i].xcenter=Segments[i].xcenter/Segments[i].size*float(mWidth);
Segments[i].ycenter=Segments[i].ycenter/Segments[i].size*float(mHeight);
if (Segments[i].phase == 1)
nr_grains++;
}
std::cout<<"Number of grains: "<<nr_grains<<std::endl;
return 0;
}
int calculate_grains(std::vector<Grains> &Grain, std::vector<Segment> &Segments, vigra::FVector3Image &vector_image, vigra::IImage &grain_label_image)
{
std::cout << "calculate grain parameters.";
for(int i=1; i<=nr_grains; i++)
{
Grain[i].size=0;
Grain[i].xcenter=0;
Grain[i].ycenter=0;
Grain[i].axis[0]=Grain[i].axis[1]=Grain[i].axis[2]=0;
Grain[i].az=0;
Grain[i].co=0;
Grain[i].r=0;
}
int grain_count = 0;
int segment_count = 0;
for(int i=1; i<=nr_regions; i++)
if (Segments[i].phase == 1)
{
grain_count++;
Segments[i].child_index = grain_count;
Grain[grain_count].parent_index = i;
Grain[grain_count].phase = Segments[i].phase;
Grain[grain_count].size = Segments[i].size;
Grain[grain_count].xcenter = Segments[i].xcenter;
Grain[grain_count].ycenter = Segments[i].ycenter;
Grain[grain_count].pointer = Segments[i].pointer;
}
//calculate mean axis vector
for(int i=1; i<=nr_grains; i++)
{
int x = 0;
int y = 0;
for(int pointer_iterator = 0; pointer_iterator < Grain[i].pointer.size(); pointer_iterator++)
{
x = Grain[i].pointer[pointer_iterator].x;
y = Grain[i].pointer[pointer_iterator].y;
if (vigra::dot(vector_image(x,y), Grain[i].axis)<0.0) Grain[i].axis-=vector_image(x,y);
else Grain[i].axis+=vector_image(x,y);
}
Grain[i].r=(2*sqrt(vigra::squaredNorm(Grain[i].axis))/Grain[i].size-1.0)*100;
Grain[i].axis=Grain[i].axis/sqrt(vigra::squaredNorm(Grain[i].axis));
}
/*
//calculate center of mass positions, mean grain size and mean axis vector
std::cout << ".";
for (int x=0; x<mWidth; x++)
for (int y=0; y<mHeight; y++)
{
//grain labeling starts with 1
if (grain_label_image(x,y)>0)
segment_count = grain_label_image(x,y);
if (Segments[segment_count].phase == 1)
{
int i = grain_label_image(x,y);
Grain[i].size++;
......@@ -1993,6 +2061,8 @@ int calculate_grains(vigra::IImage &grain_label_image, vigra::FVector3Image &vec
Grain[i].ycenter+=y/float(mHeight);
if (vigra::dot(vector_image(x,y), Grain[i].axis)<0.0) Grain[i].axis-=vector_image(x,y);
else Grain[i].axis+=vector_image(x,y);
//if (bground_image(x,y)==0) Grain[i].phase=0;
}
}
......@@ -2004,7 +2074,7 @@ int calculate_grains(vigra::IImage &grain_label_image, vigra::FVector3Image &vec
Grain[i].r=(2*sqrt(vigra::squaredNorm(Grain[i].axis))/Grain[i].size-1.0)*100;
Grain[i].axis=Grain[i].axis/sqrt(vigra::squaredNorm(Grain[i].axis));
}
*/
// Transform to geo coordinates in horizontal projection
std::cout << ".";
for(int i=1; i<=nr_grains; i++)
......@@ -2067,7 +2137,7 @@ float get_length(std::vector<point> pixels)
return length;
}
int calculate_boundaries(std::vector<Boundaries> &Boundary, std::vector<Grains> &Grain, std::string path)
int calculate_boundaries(std::vector<Boundaries> &Boundary, std::vector<Grains> &Grain, std::vector<Segment> &Segments, std::string path)
{
std::cout << "calculate boundary parameters..." << std::endl;
......@@ -2093,7 +2163,7 @@ int calculate_boundaries(std::vector<Boundaries> &Boundary, std::vector<Grains>
vigra::TinyVector<float, 3> xzplane;
xzplane[0]=0.0; xzplane[1]=1.0; xzplane[2]=0.0;
int grain1, grain2, xdist, ydist;
int grain1, grain2, segment_index1, segment_index2, xdist, ydist;
double jjdistance;
for (int i=1; i<=nr_boundaries; i++)
{
......@@ -2101,8 +2171,16 @@ int calculate_boundaries(std::vector<Boundaries> &Boundary, std::vector<Grains>
Boundary[i].length = get_length(segment.arcs[i-1]);
Boundary[i].pointer = segment.arcs[i-1]; //all points of the arc
grain1 = segment.two_boundings(i-1,0);
grain2 = segment.two_boundings(i-1,1);
segment_index1 = segment.two_boundings(i-1,0);
segment_index2 = segment.two_boundings(i-1,1);
if (Segments[segment_index1].phase == 1)
grain1 = Segments[segment_index1].child_index;
else
grain1 = 0;
if (Segments[segment_index2].phase == 1)
grain2 = Segments[segment_index2].child_index;
else
grain2 = 0;
Boundary[i].grain.push_back(grain1);
Boundary[i].grain.push_back(grain2);
......@@ -2166,7 +2244,7 @@ int calculate_boundaries(std::vector<Boundaries> &Boundary, std::vector<Grains>
return 0;
}
int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image &vector_image, std::vector<Grains> &Grain)
int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image &vector_image, std::vector<Grains> &Grain, vigra::BImage bground_image)
{
//General Statistics
std::cout << "calculate c-axes statistics.";
......@@ -2189,7 +2267,7 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
for (int x=0; x<mWidth; x++)
for (int y=0; y<mHeight; y++)
{
if (grain_label_image(x,y)>0)
if (grain_label_image(x,y)>0 && bground_image(x,y)==1)
{
nr_points++;
if (vigra::dot(vector_image(x,y), sum_vector)<0.0) sum_vector-=vector_image(x,y);
......@@ -2207,7 +2285,7 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
for (int x=0; x<mWidth; x++)
for (int y=0; y<mHeight; y++)
{
if (grain_label_image(x,y)>0)
if (grain_label_image(x,y)>0 && bground_image(x,y)==1)
{
if (vector_image(x,y)[2]<0.0) point = -vector_image(x,y);
else point = vector_image(x,y);
......@@ -2269,6 +2347,8 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
for (int i=1; i<=nr_grains; i++)
{
if (Grain[i].phase==1)
{
if (Grain[i].axis[2]<0.0) point = -Grain[i].axis;
else point = Grain[i].axis;
orientation_tensor(0,0)+=point[0]*point[0];
......@@ -2277,6 +2357,7 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
orientation_tensor(0,1)=orientation_tensor(1,0)+=point[0]*point[1];
orientation_tensor(0,2)=orientation_tensor(2,0)+=point[0]*point[2];
orientation_tensor(2,1)=orientation_tensor(1,2)+=point[2]*point[1];
}
}
vigra::linalg::symmetricEigensystem(orientation_tensor, eigenvalues, eigenvectors);
......@@ -2401,7 +2482,8 @@ int write_grains_txt (std::vector<Grains> &Grain, std::string path)
<< "# c5 mean azimuth [degree]\n"
<< "# c6 mean latitude [degree]\n"
<< "# c7 mean colatitude [degree]\n"
<< "# c8 regelungsgrad [%]\n";
<< "# c8 regelungsgrad [%]\n"
<< "# c9 phase (0 background, 1 ice)\n";
for(int i=1; i<=nr_grains; i++)
{
......@@ -2412,7 +2494,8 @@ int write_grains_txt (std::vector<Grains> &Grain, std::string path)
<< Grain[i].az << "\t"
<< 90-Grain[i].co << "\t"
<< Grain[i].co << "\t"
<< Grain[i].r << "\n";
<< Grain[i].r << "\t"
<< Grain[i].phase << "\n";
}
output_file.close();
......
......@@ -278,7 +278,7 @@ int analysis()
vigra::FImage topo_image(mWidth, mHeight);
create_topo(vector_image, topo_image);
// erase bad quality and backtround pixels from topo_image - set their values to 100
// erase bad quality and background pixels from topo_image - set their values to 100
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
if(quality_image(x,y)==0 || bground_image(x,y)==0)
......@@ -412,34 +412,49 @@ int analysis()
// create grain label image
vigra::IImage grain_label_image(mWidth,mHeight);
nr_regions = vigra::labelImageWithBackground(vigra::srcImageRange(label_image), vigra::destImage(grain_label_image), false, 0);
nr_grains = nr_regions;
std::cout<<"Number of grains: "<<nr_grains<<std::endl;
// nr_grains = nr_regions;
std::cout<<"Number of regions: "<<nr_regions<<std::endl;
if (nr_grains == 0) return 2;
if (nr_regions == 0)
{
std::cout<<"segmentation terminated: the image contains no regions."<<std::endl;
return 2;
}
// perform region classifier to distinguish between different types of regions
std::vector<Segment> Segments(nr_regions+1); //region index starts with 1
classify_segments(Segments, grain_label_image, bground_image);
// calculate grains
std::vector<Grains> Grain(nr_grains+1);
calculate_grains(Grain, Segments, vector_image, grain_label_image);
/*********************************************/
/* Region growth and extraction of */
/* topological network */
/*********************************************/
std::cout<<"region growing...";
// create a statistics functor for region growing
vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<float> > gradstat(nr_grains);
vigra::ArrayOfRegionStatistics<vigra::SeedRgDirectValueFunctor<float> > gradstat(nr_regions);
// perform region growing
vigra::IImage boundary_label_image(mWidth,mHeight);
vigra::seededRegionGrowing(vigra::srcImageRange(grain_label_image), vigra::srcImage(grain_label_image), vigra::destImage(boundary_label_image),
gradstat, vigra::CompleteGrow);
gradstat, vigra::CompleteGrow);
std::cout<<"done"<<std::endl;
/*
// mask background in boundary label image
for (int y=0; y<mHeight; y++)
for (int x=0; x<mWidth; x++)
if (bground_image(x,y)==0)
boundary_label_image(x,y) = 0;
*/
//array to store regions
unsigned int** region_image=new unsigned int* [mHeight];
for (int i=0; i<mHeight; i++) region_image[i]=new unsigned int[mWidth];
for(int y=0;y<mHeight;y++)
for(int x=0;x<mWidth;x++)
region_image[y][x]=boundary_label_image(x,y);
*/
//visualisation of the regions
vigra::IImage regions(mWidth,mHeight);
......@@ -459,8 +474,19 @@ int analysis()
y_old[x]=boundary_label_image(x,y);
}
}
//array to store regions
unsigned int** region_image=new unsigned int* [mHeight];
for (int i=0; i<mHeight; i++) region_image[i]=new unsigned int[mWidth];
for(int y=0;y<mHeight;y++)
for(int x=0;x<mWidth;x++)
region_image[y][x]=boundary_label_image(x,y);
//creating the hdf5 file
std::cout<<"create hdf5 files..."<<std::endl;
hid_t file_save=H5Fcreate(filepath_h5file.c_str(), H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
//dataspace for region_image
......@@ -532,16 +558,14 @@ int analysis()
*/
/****************************************************************************/
/* Create grain and boundary data structure and calculate c-axis statistics */
/* Create boundary data structure and calculate c-axis statistics */
/****************************************************************************/
std::vector<Grains> Grain(nr_grains+1);
calculate_grains(grain_label_image, vector_image, Grain);
std::vector<Boundaries> Boundary;
calculate_boundaries(Boundary, Grain, filepath_objects);
calculate_boundaries(Boundary, Grain, Segments, filepath_objects);
calculate_statistics(grain_label_image, vector_image, Grain);
calculate_statistics(grain_label_image, vector_image, Grain, bground_image);
/**************************************/
......@@ -566,10 +590,10 @@ int analysis()
export_grains_trendw(grain_label_image, Grain, filepath_grains_trendw); std::cout << ".";
//export boundaries in rgb color-code
export_boundary_rgb(boundary_label_image, regions, bground_image, Grain, Boundary, filepath_boundary_rgb); std::cout << ".";
export_boundary_rgb(boundary_label_image, regions, bground_image, Segments, Grain, Boundary, filepath_boundary_rgb); std::cout << ".";
//export boundaries in trend-white color-code
export_boundary_trendw(boundary_label_image, regions, bground_image, Grain, Boundary, filepath_boundary_trendw); std::cout << ".";
export_boundary_trendw(boundary_label_image, regions, bground_image, Segments, Grain, Boundary, filepath_boundary_trendw); std::cout << ".";
//export subgrain misorientation topography (testing)
export_topo_image(topo_image, filepath_topo_image); std::cout << ".";
......
Markdown is supported
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