#include "stdafx.h"
#include "Dispore_All.h"
using namespace std;
//constructor
Dispore_All::Dispore_All(DS_1b_cube* input_cube)
{
cDS_cube = input_cube;
wZ = input_cube->wZ;
wY = input_cube->wY;
wX = input_cube->wX;
slice = new DS_slice(wY, wX);
//these store individual pore results, most possible pores is wX*wY/2
porearea = new unsigned long int [wX*wY/2];
porebdry = new unsigned long int [wX*wY/2];
Rh_pore = new float [wX*wY/2];
//this is the value that is used for marking pores
mark_val = 0x80;
}
//destructor
Dispore_All::~Dispore_All()
{
delete slice;
delete porearea;
delete porebdry;
delete Rh_pore;
}
bool Dispore_All::range(unsigned short int i, unsigned short int j)
{
/*Return 1 if (n,m) are legal (row,column) indices for image X*/
if ( (i < 0) || (i >= wY))
return false;
if ( (j < 0) || (j >= wX))
return false;
return true;
}
//recursive, if large pore can use a lot of stack space
//if program crash, increase menu\project\poromedia properties\configuration properties\linker\system\stack reserve size (default 30 000 000, might need more for very large samples)
void Dispore_All::mark4_speed(unsigned short int yseed, unsigned short int xseed)
{
int temp;
//lower spot (y-1)
temp = yseed-1;
if ( ( temp >= min_Y) && (!slice->data[temp][xseed]) )
{
slice->data[temp][xseed] = mark_val;
mark4_speed((unsigned short)temp, xseed);
}
//left spot
temp = xseed-1;
if ( ( temp >= 0) && (!slice->data[yseed][temp]) )
{
slice->data[yseed][temp] = mark_val;
if (min_X > temp)
min_X = (unsigned short)temp;
mark4_speed(yseed, (unsigned short)temp);
}
//right spot
temp = xseed+1;
if ( ( temp < wX) && (!slice->data[yseed][temp]) )
{
slice->data[yseed][temp] = mark_val;
if (max_X < temp)
max_X = (unsigned short)temp;
mark4_speed(yseed, (unsigned short)temp);
}
//upper spot (y+1)
temp = yseed+1;
if ( ( temp < wY) && (!slice->data[temp][xseed]) )
{
slice->data[temp][xseed] = mark_val;
if (max_Y < temp)
max_Y = (unsigned short)temp;
mark4_speed((unsigned short)temp, xseed);
}
}
//slower (can go over the same spot many times), but uses constant memory
void Dispore_All::mark4_mem(unsigned short int yseed, unsigned short int xseed)
{
UNREFERENCED_PARAMETER(xseed); //to get rid of compiler warning message
bool y_nflag, x_nflag;
int temp;
for (unsigned short int y = yseed; y <= max_Y; y++)
{
for (unsigned short int x = min_X; x <= max_X; x++)
{
//make sure we are looking at a valid spot
if (slice->data[y][x] == mark_val)
{
y_nflag = x_nflag = false;
//y-1 spot
temp = y-1;
if ( (temp >= yseed) && (!slice->data[temp][x]) ) //check bounds, and make sure its a pore
{
slice->data[temp][x] = mark_val; //mark spot
y_nflag = true;
}
//left spot
temp = x-1;
if ( (temp >= 0) && (!slice->data[y][temp]) )
{
slice->data[y][temp] = mark_val; //mark spot
if (min_X > temp) min_X = (unsigned short)temp;
x_nflag = true;
}
//right spot
temp = x+1;
if ( (temp < wX) && (!slice->data[y][temp]) )
{
slice->data[y][temp] = mark_val; //mark spot
if (max_X < temp) max_X = (unsigned short)temp;
}
//upper spot
temp = y+1;
if ( (temp < wY) && (!slice->data[temp][x]) )
{
slice->data[temp][x] = mark_val; //mark spot
if (max_Y < temp) max_Y = (unsigned short)temp;
}
if (y_nflag || x_nflag)
{
//go to previous row if y_nflag
if (y_nflag) y--;
//if x_nflag set, the go back a column even if y_nflag not set
if (x_nflag) x-=2;
//if not, means ynflag HAS to be set, so just stay in same column
else x--;
}
}
//else go to next spot
}
}
}
//
////original mark4
//void Dispore_All::mark4_mem(unsigned short int iseed, unsigned short int jseed)
//{ //call: mark4(x,mark_val,ii,jj);
///*Mark a 4-/connected region with mark_val, starting at (iseed,jseed)*/
// //unsigned int k;
// unsigned short int imax, jmax, idec, jdec;
// register unsigned short int i,j;
// register short int n,m;
//
////determine starting "window"
// if (iseed > 0)
// iseed--;
// if (jseed > 0)
// jseed--;
// imax = iseed + 2;
// jmax = jseed + 2;
// if (imax >= wY)
// imax = wY - 1;
// if (jmax >= wX)
// jmax = wX - 1;
//
// for (i=iseed; i <= imax; i++)
// {
// for (j=jseed; j <= jmax; j++)
// {
// if (slice->data[i][j] == mark_val)
// {
// idec = jdec = 0;
// for (n=-1; n <= 1; n++)
// for (m=-1; m <= 1; m++)
// {
// if ( ((n*m) != 0) || (n==0 && m==0) )
// continue;
// if (range(i+n, j+m) == 0)
// continue;
// if ( !slice->data[i+n][j+m])
// {
// slice->data[i+n][j+m] = mark_val;
//
// //expand boundaries if needed
// if (min_X > j+m) min_X = j+m;
// if (max_X < j+m) max_X = j+m;
// if (max_Y < i+n) max_Y = i+n;
//
// if ((imax <= i+n) && (i+n+1 < wY))
// imax = i+n+1;
// if ((jmax <= j+m) && (j+m+1 < wX))
// jmax = j+m+1;
// if (n < 0)
// idec = 1;
// if (m < 0)
// jdec = 1;
// }
// }
// if ((i > iseed) * idec)
// i = i--;
// if ((j > 1) * jdec)
// j = j-2;
// if (j+1 < jseed)
// jseed = j+1;
// }
// }
// }
//}
//void mark8(struct image *x,int mark_val,int iseed,int jseed)
//{
// /*Mark an 8-connected region, beginning at (iseed, jseed), with mark_val*/
// int i,j,n,m,again;
// if(x->data[iseed][jseed]!=0)
// return;
// x->data[iseed][jseed]=(unsigned char)mark_val;
// do
// {
// again = 0;
// for(i=0;i<x->info->nr;i++)
// {
// for(j=0;j<x->info->nc;j++)
// {
// if(x->data[i][j]==mark_val)
// {
// for(n=i-1;n<=i+1;n++)
// {
// for(m=j-1;m<=j+1;m++)
// {
// if(range(x,n,m)==0)
// continue;
// if(x->data[n][m]==0)
// {
// x->data[n][m]=(unsigned char)mark_val;
// again=1;
// }
// }
// }
// }
// }
// }
// for(i=x->info->nr-1;i>=0;i--)
// {
// for(j=x->info->nc-1;j>=0;j--)
// {
// if(x->data[i][j]==mark_val)
// {
// for(n=i-1;n<=i+1;n++)
// {
// for(m=j-1;m<=j+1;m++)
// {
// if(range(x,n,m)==0)
// continue;
// if(x->data[n][m]==0)
// {
// x->data[n][m]=(unsigned char)mark_val;
// again = 1;
// }
// }
// }
// }
// }
// }
// }while(again);
//}
void Dispore_All::region_4(int *error_code)
{ //call: region_4(x,1,&err);
/*Locate a black (0) region and mark it with mark_val. 4-connected.*/
register unsigned short int i,j;
*error_code = 0;
for(i = min_Y; i < wY; i++)
{
for(j = ns_X; j < wX; j++)
if( !slice->data[i][j] )
goto pore_found;
ns_X = 0;
}
//this only happens when we are done with the current slice
*error_code = NO_REGION;
return;
pore_found:
ns_X = min_X = max_X = j;
min_Y = max_Y = i;
//mark this spot
slice->data[i][j] = mark_val;
mark4_speed(i, j);
}
//void region_8 (struct image *x, int mark_val, int *error_code)
//{
// /*Locate a black (0) region and mark it with mark_val mark_val. 8-conneceted*/
// int i,j,ii,jj;
// *error_code=0;
// ii=-1;
// jj=-1;
// for(i=0;i<x->info->nr;i++)
// {
// for(j=0;j<x->info->nc;j++)
// {
// if(x->data[i][j]==0)
// {
// ii=i;
// jj=j;
// break;
// }
// }
// if(ii>=0)
// break;
// }
// if(ii<0)
// {
// *error_code=NO_REGION;
// return;
// }
// mark8(x,mark_val,ii,jj);
//}
void Dispore_All::del_reg()
{
/*Change all pixels with mark_val mark_val to mark_val BACKGROUND*/
register unsigned int i,j;
for(i = min_Y; i <= max_Y; i++)
for(j = min_X; j <= max_X; j++)
if(slice->data[i][j] == mark_val)
slice->data[i][j] = 0xFF;
}
unsigned long int Dispore_All::area()
{
/*Count and return the number of pixels having mark_val VAL*/
register unsigned short int i,j;
unsigned long int k;
k=0;
for(i = min_Y; i <= max_Y; i++)
for(j = min_X; j <= max_X; j++)
if( slice->data[i][j] == mark_val)
k++;
return k;
}
unsigned char Dispore_All::nay4(unsigned short int i, unsigned short int j)
{
/*Return the number of 4-connected neighbors of (i,j) with mark_val VAL*/
// register short int n,m;
unsigned char k=0;
int temp;
//y-1
temp = i-1;
if ( (temp >= 0) && (slice->data[temp][j] == mark_val) )
k++;
//x-1
temp = j-1;
if ( (temp >= 0) && (slice->data[i][temp] == mark_val) )
k++;
//x+1
temp = j+1;
if ( (temp < wX) && (slice->data[i][temp] == mark_val) )
k++;
//y+1
temp = i+1;
if ( (temp < wY) && (slice->data[temp][j] == mark_val) )
k++;
return k;
// if( slice->data[i][j] != mark_val)
// return 0;
//
// k=0;
// for(n = -1; n <= 1; n++)
// for(m = -1; m <= 1; m++)
// {
// if ( ((n*m) != 0) || (n==0 && m==0) )
// continue;
// if(range(i+n, j+m))
// if( slice->data[i+n][j+m] == mark_val)
// k++;
// }
// return k; //k-1
}
//int nay8(struct image *x,int i,int j,int marked_val)
//{
// /*Return the number of 8-connected neighbors of (i,j) having mark_val marked_val*/
// int n,m,k;
// if(x->data[i][j]!=marked_val)
// return 0;
// k=0;
// for(n=-1;n<=1;n++)
// {
// for(m=-1;m<=1;m++)
// {
// if(range(x,i+n,j+m))
// if(x->data[i+n][j+m]==marked_val)
// k++;
// }
// }
// return k-1;
//}
unsigned long int Dispore_All::boundary()
{
/*Compute the boundary of the region(s) marked with marked_val*/
short int k;
register unsigned short int i,j;
unsigned long int p=0;
/*Remove all pixels except those having mark_val marked_val*/
for(i = min_Y; i <= max_Y; i++)
for(j = min_X; j <= max_X; j++)
if( slice->data[i][j] == mark_val)
{
k = nay4(i, j);/* How many neighbors are marked_val*/
if(k < 4)/*If not all, this is on boundary*/
p += (4-k); //p=p+(4-k)*(float)1.0;
}
return p;
}
//float perimeter(struct image *x,int marked_val,int *error_code)
//{
// /*Compute the perimeter of the region(s) marked with marked_val*/
// unsigned int k,t;
// register unsigned int i,j;
// float p;
// struct image *y;
// *error_code = 0;
// p = 0.0;
// y = 0;
// copy(x,&y,error_code);
// if(*error_code)
// return 0.0;
// /*Remove all pixels except those having mark_val marked_val*/
// for(i=0;i<y->info->nr;i++)
// {
// for(j=0;j<y->info->nc;j++)
// {
// if(x->data[i][j]!=marked_val)
// {
// y->data[i][j]=BACKGROUND;
// continue;
// }
// k=nay4(x,i,j,marked_val);/*How many neighbors are marked_val*/
// if(k<4)/*If not all, this is on perim*/
// y->data[i][j]=0;
// else
// y->data[i][j]=BACKGROUND;
// }
// }
// for(i=0;i<y->info->nr;i++)
// {
// for(j=0;j<y->info->nc;j++)
// {
// if(y->data[i][j]!=0)
// continue;
// /*Match one of the templates*/
// k=1;
// t=0;
// for(register int ii=-1;ii<=1;ii++)
// for(register int jj=-1;jj<=1;jj++)
// {
// if(ii==0&&jj==0)
// continue;
// if(y->data[i+ii][j+jj]==0)
// t=t+k;
// k=k<<1;
// }
///*
// Templates for 1.207:
// o o o o o # o # o o # o # o o o o # o o o # o o
// # # o # # o o # o o # o o # o o # o o # # o # #
// o o # o o o # o o o o # o # o o # o # o o o o o
// T= 210 014 042 202 101 104 060 021
//
// Templates for 1.414:
// # o o o o # # o o o o # o o o # o #
// o P o o P o o P o o P o o P o o P o
// o o # # o o # o o o o # # o # o o o
// T= 201 044 041 204 240 005
//
// Templates for 1.0:
//
// o o o o # o o o o o o o o # o o # o
// # # # o # o # # o o # # # # o o # #
// o o o o # o o # o o # o o o o o o o
// 030 102 72 80 10 18
//*/
// if(t==0210||t==014||t==042||t==0202||t==0101||t==0104||t==060||t==021)
// {
// p+=(float)1.207;
// continue;
// }
//
// if(t==0201||t==044||t==041||t==0204||t==0240||t==005)
// {
// p+=(float)1.414;
// continue;
// }
//
// if(t==030||t==0102||t==72||t==80||t==10||t==18)
// {
// p+=1.0;
// continue;
// }
//
// p+=(float)1.207;
// }
// }
// freeimage(y,error_code);
// return p;
//}
//returns:
// 0 if all good,
// 10 if stop received
int Dispore_All::run_dispore_all(char* out_prefix, char* plane, CProgressCtrl *pbar_ptr, bool *stop_req)
{
register unsigned int porenmb;
int err;
register unsigned short int index;
register unsigned long int ind;
unsigned long long sumarea, sumbdry; //__int64
char fnamew[512];
//progress control
pbar_ptr->SetPos(0);
pbar_ptr->SetStep(1);
pbar_ptr->SetRange(0, wZ);
{
char log_txt[512];
sprintf(log_txt, "Porosity: Dispore_all: plane %s, slices: %d", plane, wZ);
write_to_log(log_txt);
}
//creating the --_po_bdry.txt file
char po_bdry_out_file[512];
sprintf(po_bdry_out_file, "%s_%s_po_bdry.txt\0", out_prefix, plane);
ofstream fout(po_bdry_out_file);
if(!fout.is_open())
{
char log_txt[512];
sprintf(log_txt, "Dispore_All: cannot open file %s.", po_bdry_out_file);
return 1;
}
fout << "%index\t\tsumarea\t\tsumbdry\t\tratio" << endl;
//create the _po.dis file
sprintf(fnamew,"%s_%s.dis", out_prefix, plane);
ofstream ofilew(fnamew);
if(!ofilew.is_open())
{
char log_txt[512];
sprintf(log_txt, "Dispore All: Cannot open file %s.", fnamew);
write_to_log(log_txt);
return 1;
}
//allocating needed memory
//maxporenmb = 262000;
// porearea=(unsigned long int*)malloc(sizeof(unsigned long int)*maxporenmb);
// porebdry=(unsigned long int*)malloc(sizeof(unsigned long int)*maxporenmb);
// Rh_pore=(float*)malloc(sizeof(float)*maxporenmb);
// if( (!Rh_pore) | (!porebdry) | (!porearea) )
// {
// write_to_log("Dispore_All: out of memory -Rh_pore- or -porebdry- or -porearea-");
// if (porearea) free(porearea);
// if (porebdry) free(porebdry);
// if (Rh_pore) free(Rh_pore);
// fout.close();
// ofilew.close();
// return 2;
// }
for(index = 1; index <= wZ; index++)
{
if (index%10 == 0)
{
char log_txt[512];
sprintf(log_txt, "Porosity: Dispore_all: processing slice #%d", index);
write_to_log(log_txt);
}
//reset variables for new slice
sumbdry = sumarea = err = 0;
porenmb = 1;
min_Y = ns_X = 0;
cDS_cube->to_DS_slice(slice, index-1);
/* process slice */
do
{
region_4(&err); //locate pore
if(err)
break;
sumarea += porearea[porenmb] = area(); //area of that pore
sumbdry += porebdry[porenmb] = boundary(); //boundary (perimeter) of that pore
Rh_pore[porenmb] = (float)porearea[porenmb] / (float)porebdry[porenmb];
porenmb++;
del_reg(); //erase pore (-> fiber)
}
while (err == 0);
//write slice header
ofilew << "%%Slice number: " << index << endl;
ofilew << "%%Pore" << endl;
ofilew << "%%No.\tArea\tPerimeter\tRh";
//if no pores found, skip to next slice
if (porenmb == 1)
{
if (index != wZ) //to prevent blank lines at end of file
ofilew << endl << endl;
}
else
{ //write .dis file
ofilew << endl;
for(ind = 1; ind < (porenmb-1); ind++)
ofilew << ind << "\t" << porearea[ind] << "\t" << porebdry[ind] << "\t" << Rh_pore[ind] << endl;
ofilew << ind << "\t" << porearea[ind] << "\t" << porebdry[ind] << "\t" << Rh_pore[ind];
if (index != wZ) //to prevent blank lines at end of file
ofilew << endl << endl;
//write --_po_bdry.txt (sum area, sum boundary)
fout << index << "\t\t" << sumarea << "\t\t" << sumbdry << "\t\t" << (double)sumarea/(double)sumbdry << endl;
}
//update progress
pbar_ptr->StepIt();
if (*stop_req)
break;
}
fout.close();
ofilew.close();
if (*stop_req)
return 10;
else
return 0;
}