#include "stdafx.h"
#include "stats_pore_all.h"
//
//#include "mathtool.h"
using namespace std;
int posround(float number)
{
int a=(int)number;
int b=a+1;
if(fabs(a-number)<=fabs(b-number))
return a;
else
return b;
}
Val_instance::Val_instance()
{
value=0;
instance=0;
}
void bubblunsint(unsigned int* a, unsigned int n)
{ //bubble sort algorithm
register unsigned int i,j;
for(i=0;i<n;i++)
for(j=1;j<(n-i);j++)
if(a[j-1]>a[j])
{
unsigned int t;
t=a[j-1];
a[j-1]=a[j];
a[j]=t;
}
}
List<Val_instance>* tabulate(List<unsigned int>* inputlist)
{
List<Val_instance>* outputlist=new List<Val_instance>;
unsigned int *a=(unsigned int*)malloc(sizeof(unsigned int)*inputlist->getNumber());
if(!a)
{
cout << "out of memory -a-" << endl;
exit(1);
}
for(register int i=0;i<inputlist->getNumber();i++)
a[i]=*(inputlist->getWhich(i+1));
bubblunsint(a,inputlist->getNumber());
Val_instance* first=new Val_instance;
first->value=a[0];
first->instance=1;
outputlist->Add(*first);
for(register int i=1;i<inputlist->getNumber();i++)
if((int)a[i] == outputlist->getWhich(1)->value)
outputlist->getWhich(1)->instance++;
else
{
Val_instance* newsample=new Val_instance;
newsample->value=a[i];
newsample->instance=1;
outputlist->Add(*newsample);
}
return outputlist;
}
unsigned int mymax(List<unsigned int>* inputlist)
{
unsigned int result=0;
for(register int i=1; i<=inputlist->getNumber(); i++)
if(*(inputlist->getWhich(i)) > result)
result=*(inputlist->getWhich(i));
return result;
}
float fmax(List<float>* inputlist)
{
float result=0.0;
for(register int i=1;i<=inputlist->getNumber();i++)
if(*(inputlist->getWhich(i))>result)
result=*(inputlist->getWhich(i));
return result;
}
int stats_pore_all(DS_1b_cube* info, short bins, char* output_prefix, char* plane)
{
unsigned short int total_slices;
register unsigned short int i;//, line;
char* filehead;
float max_rh_array=0.0;
char stats[512];
char data_file_name[512];
unsigned long long sum_area = 0;
double mean_Rh=0.0;
float mid_Rh;
ifstream fp, fp2;
unsigned int* Rh_array=(unsigned int*)malloc(sizeof(unsigned int)*(bins+1));
//what plane to process
switch(*plane)
{
case'x':
total_slices=info->wZ;
break;
case'y':
total_slices=info->wX;
break;
case'z':
total_slices=info->wY;
break;
default:
total_slices = 0;
}
sprintf(data_file_name,"%s_%s.dis", output_prefix, plane);
fp.open(data_file_name, ios::in);
if(!fp)
{
char log_txt[512];
sprintf(log_txt, "Stats Pore All: Cannot open file %s.", data_file_name);
write_to_log(log_txt);
return 1;
}
{
char log_txt[512];
sprintf(log_txt, "Porosity: stats_pore_all: plane %s, slices: %d", plane, total_slices);
write_to_log(log_txt);
}
for(i=1; i<=total_slices; i++)
{
//read .dis file
//sprintf(data_file_name,"%s_%s_%d.dis", output_prefix/*info->bin_prefix*/,plane,i);
//ifstream fp(data_file_name,ios::in);
//if(!fp)
//{
// char log_txt[512];
// sprintf(log_txt, "Stats Pore All: Cannot open file %s.", data_file_name);
// write_to_log(log_txt);
// return 1;
//}
filehead=(char*)malloc(sizeof(char)*256);
if(!filehead)
{
write_to_log("Stats Pore All: out of memory -filehead-");
return 2;
}
//skip header (3 first lines)
for(unsigned int line=1; line<=3; line++)
fp.getline(filehead, 256, '\n');
if(filehead) free(filehead);
//read data: search MAX Rh
unsigned long int one,two,three;
float four;
while( (!fp.eof()) && (fp.get() != '%') )/*(.eof()*/
{
fp.unget();
if (fp.get() != '\n')
{
fp.unget();
fp >> one >> two >> three >> four;
if (four > max_rh_array) max_rh_array=four;
}
}
}
fp.close();
//create array for stats
if(!Rh_array)
{
write_to_log("Stats Pore All: out of memory -Rh_array-");
return 2;
}
//initialise stats variables
for(i=0; i<=bins; i++)
Rh_array[i]=0;
// sprintf(data_file_name,"%s_%s.dis",output_prefix, plane);
fp2.open(data_file_name, ios::in); //, ios_base::binary);
if(!fp2)
{
char log_txt[512];
sprintf(log_txt, "Stats Pore All: cannot open file %s", data_file_name);
write_to_log(log_txt);
return 1;
}
//fp.seekg(0, ios::beg);
//read data
for(i=1; i<=total_slices; i++)
{
//read single .dis file
//sprintf(data_file_name,"%s_%s_%d.dis",output_prefix, plane, i);
//ifstream fp(data_file_name); //, ios_base::binary);
//if(!fp)
//{
// char log_txt[512];
// sprintf(log_txt, "Stats Pore All: cannot open file %s", data_file_name);
// write_to_log(log_txt);
// return 1;
//}
filehead=(char*)malloc(sizeof(char)*256);
if(!filehead)
{
write_to_log("Stats Pore All: out of memory -filehead-");
return 2;
}
//skip header (3 first lines)
for(unsigned short int line=1; line<=3; line++)
fp2.getline(filehead, 256, '\n');
if (filehead)
free(filehead);
//unsigned int max_num=0;
//read data + generate histogram (Rh_array)
unsigned long int one,two,three;
float four;
while( (!fp2.eof()) && (fp2.get() != '%') )/*(.eof()*/
{
fp2.unget();
if (fp2.get() != '\n')
{
fp2.unget();
fp2 >> one >> two >> three >> four; //%%No. Area Perimeter Rh
Rh_array[(unsigned int)(four*bins/max_rh_array)] += two; //two = area
mean_Rh += two*four;
sum_area += two;
}
}
//////////////////////////////////
//while(!fp.eof())
//{
// unsigned int one, two, three; //two=area
// float four; //Rh
// fp>>one>>two>>three>>four;
// Rh_array[(unsigned int)(four*100/max_rh_array)] += two; //two = area
// mean_Rh += two*four;
// sum_area += two;
//}
}
fp2.close();
//stats results
mean_Rh /= sum_area;
//in array Rh_array subscript is Rh*100/max_rh_array and value is corresponding sum area
//merge Rh_array[bins-1] & Rh_array[bins]
Rh_array[bins-1] += Rh_array[bins];
//print results
sprintf(stats,"%s_%s_stats_rh-area.txt", output_prefix, plane);
ofstream fid_stats(stats, ios::out);
if(!fid_stats.is_open())
{
char log_txt[512];
sprintf(log_txt, "Stats Pore All : cannot open file %s", stats);
write_to_log(log_txt);
return 1;
}
fid_stats<<"%% Hydraulic Radii (per Area) distribution for the whole sample"<<endl;
fid_stats<<"%% Mean_Rh = "<<mean_Rh<<endl;
fid_stats<<"%% Sum_Area = "<<sum_area<<endl;
fid_stats<<"%% Rh %Sum_Area"<<endl;
for(i=0; i<bins; i++)
{
mid_Rh = ((float)i+0.5f)*max_rh_array/bins;
fid_stats << mid_Rh << " " << Rh_array[i]/(double)sum_area<<endl;
}
fid_stats.close();
//release memory
if (Rh_array)
free(Rh_array);
return 0;
}