Commit 0012f61f authored by Jan Eichler's avatar Jan Eichler
Browse files

working on issue #1

parent 1fb78c90
......@@ -4,6 +4,7 @@
#include <fstream>
#include <string>
#include <sys/stat.h>
#include <ctime>
#include <vigra/stdimage.hxx>
#include <vigra/impex.hxx>
......@@ -22,6 +23,8 @@
#include "marray.hxx"
#include "seg.h"
std::string cAxes_version = "21.02.11 slim";
class Grains
{
public:
......@@ -82,7 +85,10 @@ vigra::BImage quality_image(1,1);
//grain statistics
int nr_grains, mean_size;
float sum_vector_norm, regelungsgrad, regelungsgrad_g, concentration_parameter, concentration_parameter_g, spherical_aperture, spherical_aperture_g, woodcock, woodcock_g, e1, e2, e3, e1_g, e2_g, e3_g;
float sum_vector_norm, regelungsgrad, regelungsgrad_g, concentration_parameter, concentration_parameter_g, spherical_aperture, spherical_aperture_g, woodcock, woodcock_g,
e1, e2, e3, e1_g, e2_g, e3_g;
float eigenv_a_az[3], eigenv_g_az[3], eigenv_a_co[3], eigenv_g_co[3]; // (az_e1, az_e2, az_e3) and so on
vigra::TinyVector <float, 3> eigenvector_a[3], eigenvector_g[3];
//boundary statistics
int nr_boundaries;
......@@ -2201,9 +2207,42 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
}
vigra::linalg::symmetricEigensystem(orientation_tensor, eigenvalues, eigenvectors);
e1 = eigenvalues(2,0)/nr_points;
e2 = eigenvalues(1,0)/nr_points;
e3 = eigenvalues(0,0)/nr_points;
e1 = eigenvalues(2)/nr_points;
e2 = eigenvalues(1)/nr_points;
e3 = eigenvalues(0)/nr_points;
// eigenvectors area weighted
for (int i=0; i<=2; i++) // read out matrix and convert eigenvectors in tinyvectors
{
eigenvector_a[0][i] = eigenvectors(i,2);
eigenvector_a[1][i] = eigenvectors(i,1);
eigenvector_a[2][i] = eigenvectors(i,0);
}
// Transform to geo coordinates in horizontal projection
for (int e=0; e<3; e++)
{
if (eigenvector_a[e][2]>0.0)
eigenvector_a[e]=-eigenvector_a[e]; //rotate all vectors in negative z-direction ("down") - projection onto the lower hemisphere
if (eigenvector_a[e][1]>0.0 && eigenvector_a[e][0]>=0.0)
eigenv_a_az[e]=atan(eigenvector_a[e][0]/eigenvector_a[e][1])*180/M_PI;
else if (eigenvector_a[e][1]<0.0)
eigenv_a_az[e]=(atan(eigenvector_a[e][0]/eigenvector_a[e][1])+M_PI)*180/M_PI;
else if (eigenvector_a[e][1]>0.0 && eigenvector_a[e][0]<0.0)
eigenv_a_az[e]=(atan(eigenvector_a[e][0]/eigenvector_a[e][1])+2*M_PI)*180/M_PI;
else if (eigenvector_a[e][0]>0.0) eigenv_a_az[e]=90;
else if (eigenvector_a[e][0]<0.0) eigenv_a_az[e]=270;
else eigenv_a_az[e]=0;
eigenv_a_co[e]=acos(-eigenvector_a[e][2]/sqrt(vigra::squaredNorm(eigenvector_a[e])))*180/M_PI; //acos(-z) because lower hemisphere projection
}
std::cout << std::endl;
for (int e=0; e<3; e++)
std::cout << eigenv_a_az[e] << "\t" << eigenv_a_co[e] << std::endl;
woodcock = log(e3/e2)/log(e2/e1);
......@@ -2231,17 +2270,53 @@ int calculate_statistics(vigra::IImage &grain_label_image, vigra::FVector3Image
e2_g = eigenvalues(1,0)/nr_grains;
e3_g = eigenvalues(0,0)/nr_grains;
// eigenvectors grain weighted
for (int i=0; i<=2; i++) // read out matrix and convert eigenvectors in tinyvectors
{
eigenvector_g[0][i] = eigenvectors(i,2);
eigenvector_g[1][i] = eigenvectors(i,1);
eigenvector_g[2][i] = eigenvectors(i,0);
}
// Transform to geo coordinates in horizontal projection
for (int e=0; e<3; e++)
{
if (eigenvector_g[e][2]>0.0)
eigenvector_g[e]=-eigenvector_g[e]; //rotate all vectors in negative z-direction ("down") - projection onto the lower hemisphere
if (eigenvector_g[e][1]>0.0 && eigenvector_g[e][0]>=0.0)
eigenv_g_az[e]=atan(eigenvector_g[e][0]/eigenvector_g[e][1])*180/M_PI;
else if (eigenvector_g[e][1]<0.0)
eigenv_g_az[e]=(atan(eigenvector_g[e][0]/eigenvector_g[e][1])+M_PI)*180/M_PI;
else if (eigenvector_g[e][1]>0.0 && eigenvector_g[e][0]<0.0)
eigenv_g_az[e]=(atan(eigenvector_g[e][0]/eigenvector_g[e][1])+2*M_PI)*180/M_PI;
else if (eigenvector_g[e][0]>0.0) eigenv_g_az[e]=90;
else if (eigenvector_g[e][0]<0.0) eigenv_g_az[e]=270;
else eigenv_g_az[e]=0;
eigenv_g_co[e]=acos(-eigenvector_g[e][2]/sqrt(vigra::squaredNorm(eigenvector_g[e])))*180/M_PI; //acos(-z) because lower hemisphere projection
}
std::cout << std::endl;
for (int e=0; e<3; e++)
std::cout << eigenv_g_az[e] << "\t" << eigenv_g_co[e] << std::endl;
woodcock_g = log(e3_g/e2_g)/log(e2_g/e1_g);
std::cout << " done" << std::endl;
return 0;
}
int write_grains_txt (std::vector<Grains> &Grain, std::string path)
int write_general_results (std::string path)
{
//write output file
std::ofstream output_file(path.c_str());
std::time_t now = std::time(0);
char* time_now = std::ctime(&now);
output_file << "# ++ DATA INFO ++\n"
<< "# instrument = " << instrument << "\n"
<< "# section folder = " << sourcepath << "\n"
......@@ -2257,25 +2332,54 @@ int write_grains_txt (std::vector<Grains> &Grain, std::string path)
<< "# scanned height = " << scanned_height << "mm\n";
output_file << "\n# ++ EVALUATION INFO ++\n"
<< "# data processed: " << time_now
<< "# cAxes version = " << cAxes_version << "\n"
<< "# minimal geometric quality [%] = " << gqmin << "\n"
<< "# minimal retardation quality [%] = " << rqmin << "\n"
<< "# misorientation threshold [degree] = " << maximal_misorientation << "\n"
<< "# minimal grain size [pixel] = " << minimal_grain_size << "\n"
<< "# tile rotation angle [degree] = " << rot_angle << "\n"
<< "# section geometry (0 horizontal, 1 vertical) = " << vertical_cut << "\n"
<< "# horizontal orientation angle [degree] " << rotate_north << "\n"
<< "\n# ++ GENERAL C-AXIS STATISTICS ++\n"
<< "\n# ++ GRAIN STATISTICS ++\n"
<< "# nr of grains = " << nr_grains << "\n"
<< "# mean grain size [pixel] = " << mean_size << "\n"
<< "\n# ++ C-AXIS ORIENTATION STATISTICS ++\n"
<< "# sum vector norm = " << sum_vector_norm << "\n"
<< "\n# weighted by grain area: \n"
<< "# regelungsgrad = " << regelungsgrad << "\n"
<< "# concentration parameter = " << concentration_parameter << "\n"
<< "# spherical aperture [degree]= " << spherical_aperture << "\n"
<< "# eigenvalue e1 = " << e1 << "\n"
<< "# eigenvalue e2 = " << e2 << "\n"
<< "# eigenvalue e3 = " << e3 << "\n"
<< "# spherical aperture [degree] = " << spherical_aperture << "\n"
<< "# eigenvalues (e1, e2, e3) = (" << e1 << ", " << e2 << ", " << e3 << ")\n"
<< "# eigenvector_1 (azimuth, latitude, colatitude) [degree] = (" << eigenv_a_az[0] << ", " << 90-eigenv_a_co[0] << ", " << eigenv_a_co[0] << ")\n"
<< "# eigenvector_2 (azimuth, latitude, colatitude) [degree] = (" << eigenv_a_az[1] << ", " << 90-eigenv_a_co[1] << ", " << eigenv_a_co[1] << ")\n"
<< "# eigenvector_3 (azimuth, latitude, colatitude) [degree] = (" << eigenv_a_az[2] << ", " << 90-eigenv_a_co[2] << ", " << eigenv_a_co[2] << ")\n"
<< "# woodcock parameter = " << woodcock << "\n"
<< "\n# ++ GRAIN PROPERTIES ++\n"
<< "\n# not weighted (per grain): \n"
<< "# regelungsgrad = " << regelungsgrad_g << "\n"
<< "# concentration parameter = " << concentration_parameter_g << "\n"
<< "# spherical aperture [degree] = " << spherical_aperture_g << "\n"
<< "# eigenvalues (e1, e2, e3) = (" << e1_g << ", " << e2_g << ", " << e3_g << ")\n"
<< "# eigenvector_1 (azimuth, latitude, colatitude) [degree] = (" << eigenv_g_az[0] << ", " << 90-eigenv_g_co[0] << ", " << eigenv_g_co[0] << ")\n"
<< "# eigenvector_2 (azimuth, latitude, colatitude) [degree] = (" << eigenv_g_az[1] << ", " << 90-eigenv_g_co[1] << ", " << eigenv_g_co[1] << ")\n"
<< "# eigenvector_3 (azimuth, latitude, colatitude) [degree] = (" << eigenv_g_az[2] << ", " << 90-eigenv_g_co[2] << ", " << eigenv_g_co[2] << ")\n"
<< "# woodcock parameter = " << woodcock_g << "\n";
output_file.close();
return 0;
}
int write_grains_txt (std::vector<Grains> &Grain, std::string path)
{
//write output file
std::ofstream output_file(path.c_str());
output_file << "# ++ GRAIN PROPERTIES ++\n"
<< "# c1 grain number\n"
<< "# c2 size [pixel]\n"
<< "# c3 center x-coordinate\n"
......
......@@ -2,7 +2,7 @@
#include "misorient.cpp"
#include <unistd.h>
std::string cAxes_version = "20.11.11 slim";
//functions
int analysis();
......@@ -12,7 +12,7 @@ int main(int argc, char ** argv)
if (argc==1) //command line modus
{
//hello
std::cout << "cAxes is a program for the processing of Fabric Analyser data. \n";
std::cout << "cAxes is a program for the processing of the Fabric Analyser data. \n";
std::cout << "Version: " << cAxes_version << "\n\n";
//input parameters
......@@ -26,9 +26,9 @@ int main(int argc, char ** argv)
std::cin >> minimal_grain_size;
std::cout << "Enter tiles rotation angle in degrees: ";
std::cin >> rot_angle;
std::cout << "Enter section type (0 = horizontal; 1 = vertical): ";
std::cout << "Enter section type (0 = horizontal; 1 = vertical): ";
std::cin >> vertical_cut;
std::cout << "Rotate North in degrees: ";
std::cout << "Rotate North in degrees: ";
std::cin >> rotate_north;
//std::cout << "Eliminate tile borders? (0 = no; 1 = yes): ";
//std::cin >> eliminate_borders;
......@@ -41,32 +41,63 @@ int main(int argc, char ** argv)
{
std::ofstream statistics("statistics.txt");
std::time_t now = std::time(0);
char* time_now = std::ctime(&now);
statistics << "# ++ EVALUATION INFO ++\n"
<< "# data processed: " << time_now
<< "# cAxes version = " << cAxes_version << "\n"
<< "# minimal geometric quality [%] = " << gqmin << "\n"
<< "# minimal retardation quality [%] = " << rqmin << "\n"
<< "# misorientation threshold [degree] = " << maximal_misorientation << "\n"
<< "# minimal grain size [pixel] = " << minimal_grain_size << "\n"
<< "# NOTE: the last three columns are eigenvalues grain weighted\n"
<< "# tile rotation angle [degree] = " << rot_angle << "\n"
<< "# section geometry (0 horizontal, 1 vertical) = " << vertical_cut << "\n"
<< "# horizontal orientation angle [degree] " << rotate_north << "\n"
<< "\n# ++ GENERAL C-AXIS STATISTICS ++\n"
<< "\n# ++ COLUMN HEADERS ++\n"
<< "# c1 data directory\n"
<< "# c2 nr of grains\n"
<< "# c3 mean grain size [pixel]\n"
<< "# c4 sum vector norm\n"
<< "### weighted by grain area:\n"
<< "# c5 regelungsgrad\n"
<< "# c6 regelungsgrad g-weighted\n"
<< "# c7 concentration_parameter\n"
<< "# c8 concentration_parameter g-weighted\n"
<< "# c9 spherical_aperture [degree]\n"
<< "# c10 spherical_aperture g-weighted [degree]\n"
<< "# c11 eigenvalue e1\n"
<< "# c12 eigenvalue e1 g-weighted\n"
<< "# c13 eigenvalue e2\n"
<< "# c14 eigenvalue e2 g-weighted\n"
<< "# c15 eigenvalue e3\n"
<< "# c16 eigenvalue e3 g-weighted\n"
<< "# c17 woodcock parameter\n"
<< "# c18 woodcock parameter g-weighted\n";
<< "# c6 concentration_parameter\n"
<< "# c7 spherical_aperture [degree]\n"
<< "# c8 eigenvalue e1\n"
<< "# c9 eigenvalue e2\n"
<< "# c10 eigenvalue e3\n"
<< "# c11 eigenvector_1 azimuth [degree]\n"
<< "# c12 eigenvector_1 latitude [degree]\n"
<< "# c13 eigenvector_1 colatitude [degree]\n"
<< "# c14 eigenvector_2 azimuth [degree]\n"
<< "# c15 eigenvector_2 latitude [degree]\n"
<< "# c16 eigenvector_2 colatitude [degree]\n"
<< "# c17 eigenvector_3 azimuth [degree]\n"
<< "# c18 eigenvector_3 latitude [degree]\n"
<< "# c19 eigenvector_3 colatitude [degree]\n"
<< "# c20 woodcock parameter\n"
<< "### not weighted (per grain):\n"
<< "# c21 regelungsgrad\n"
<< "# c22 concentration_parameter\n"
<< "# c23 spherical_aperture [degree]\n"
<< "# c24 eigenvalue e1\n"
<< "# c25 eigenvalue e2\n"
<< "# c26 eigenvalue e3\n"
<< "# c27 eigenvector_1 azimuth [degree]\n"
<< "# c28 eigenvector_1 latitude [degree]\n"
<< "# c29 eigenvector_1 colatitude [degree]\n"
<< "# c30 eigenvector_2 azimuth [degree]\n"
<< "# c31 eigenvector_2 latitude [degree]\n"
<< "# c32 eigenvector_2 colatitude [degree]\n"
<< "# c33 eigenvector_3 azimuth [degree]\n"
<< "# c34 eigenvector_3 latitude [degree]\n"
<< "# c35 eigenvector_3 colatitude [degree]\n"
<< "# c36 woodcock parameter\n"
<< "\n# ++ DATA ++\n";
statistics.close();
......@@ -82,18 +113,24 @@ int main(int argc, char ** argv)
<< mean_size << "\t"
<< sum_vector_norm << "\t"
<< regelungsgrad << "\t"
<< regelungsgrad_g << "\t"
<< concentration_parameter << "\t"
<< concentration_parameter_g << "\t"
<< spherical_aperture << "\t"
<< spherical_aperture_g << "\t"
<< e1 << "\t"
<< e1_g << "\t"
<< e2 << "\t"
<< e2_g << "\t"
<< e3 << "\t"
<< e3_g << "\t"
<< eigenv_a_az[0] << "\t" << 90-eigenv_a_co[0] << "\t" << eigenv_a_co[0] << "\t"
<< eigenv_a_az[1] << "\t" << 90-eigenv_a_co[1] << "\t" << eigenv_a_co[1] << "\t"
<< eigenv_a_az[2] << "\t" << 90-eigenv_a_co[2] << "\t" << eigenv_a_co[2] << "\t"
<< woodcock << "\t"
<< regelungsgrad_g << "\t"
<< concentration_parameter_g << "\t"
<< spherical_aperture_g << "\t"
<< e1_g << "\t"
<< e2_g << "\t"
<< e3_g << "\t"
<< eigenv_g_az[0] << "\t" << 90-eigenv_g_co[0] << "\t" << eigenv_g_co[0] << "\t"
<< eigenv_g_az[1] << "\t" << 90-eigenv_g_co[1] << "\t" << eigenv_g_co[1] << "\t"
<< eigenv_g_az[2] << "\t" << 90-eigenv_g_co[2] << "\t" << eigenv_g_co[2] << "\t"
<< woodcock_g << "\n";
}
else statistics << sourcepath << " error status " << status << "\n";
......@@ -186,6 +223,7 @@ int analysis()
std::string filepath_boundary_rgb = new_directory+"/boundary_rgb.png";
std::string filepath_boundary_trendw = new_directory+"/boundary_trendw.png";
std::string filepath_regions = new_directory+"/regions.bmp";
std::string filepath_general_txt = new_directory+"/general_results.txt";
std::string filepath_grains_txt = new_directory+"/grains.txt";
std::string filepath_grains2_txt = new_directory+"/grains2.txt";
std::string filepath_boundaries_txt = new_directory+"/boundaries.txt";
......@@ -198,6 +236,8 @@ int analysis()
/*********************************************************/
/* Import data and create vector_image and quality_image */
/*********************************************************/
std::cout << "\n";
int import_status;
import_status = import_data_g50AWI();
......@@ -545,6 +585,8 @@ int analysis()
/**************************************/
std::cout << "Write text files.";
write_general_results(filepath_general_txt); std::cout << ".";
write_grains_txt(Grain, filepath_grains_txt); 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