porositySSA.cpp #2

#include "stdafx.h"
#include "porositySSA.h"

using namespace std;

/*
 if return is:
 2, stop received
 1, could not open output file
 0, everything good
*/

short int porositySSA(DS_1b_cube * DS_cube, char* output_prefix, void* Dlg)
{
	PorositySSADialog * PoroDlg = (PorositySSADialog *) Dlg;

	int threshold =	PoroDlg->GetIntValFromField(&PoroDlg->field_manual_threshold);
	float pixsize = PoroDlg->field_pixel_size_float;
//init progress bar
	PoroDlg->pbar_simulation_progress.SetPos(0);
	PoroDlg->pbar_simulation_progress.SetStep(1);
	PoroDlg->pbar_simulation_progress.SetRange(0, DS_cube->wZ);

//create SSA~sample_name~.txt
    char filename[512];
	sprintf(filename,"%s_SSA.txt", output_prefix);
    ofstream fp(filename);
    if(!fp.is_open())
    {
		char log_txt[512];
		sprintf(log_txt, "porositySSA: cannot open file %s", filename);
        write_to_log(log_txt);
        return 1;
    }
	fp<<"Sample "<< output_prefix <<endl;
   
	PoroDlg->txt_field_progress.SetWindowText("Running Porosity...");
//counting pores (black,0) 
	unsigned long long int pore=0;
	register unsigned short int i, j, k;
	unsigned long poro_table [4096];	//poro per slice (max 4096 slices)
	for (i = 0; i < DS_cube->wZ; i++)
	{
		unsigned long long int pore_temp = pore;
		for (j = 0; j < DS_cube->wY; j++)
			for (k = 0; k < DS_cube->wX; k++)
				if (!DS_cube->get_spot(i, j, k))//ds[aa]==0)	//if pore (0, black)
					pore++;
		poro_table[i] = (unsigned long)(pore - pore_temp);
        PoroDlg->pbar_simulation_progress.StepIt();
	}

//find interface pixels
    unsigned long long int interfaces=0;
    
	
	PoroDlg->txt_field_progress.SetWindowText("Running SSA...");
	//changed 2 to 1 as starting
	for (i = 1; i <= DS_cube->wZ - 2; i++)
	{
        for (j=1; j <= DS_cube->wY - 2; j++)
		{
            for (k=1; k <= DS_cube->wX - 2; k++)
			{
				if (DS_cube->get_spot(i, j, k)) //if fiber(white)
				{
					if(!DS_cube->get_spot(i-1, j, k) )//*(ds+(i-2)*info->sizeofx*info->sizeofy+(j-1)*info->sizeofx+(k-1)))==0)	//z-
                        interfaces++;
                    if(!DS_cube->get_spot(i, j-1, k) )//(int)(*(ds+(i-1)*info->sizeofx*info->sizeofy+(j-2)*info->sizeofx+(k-1)))==0)	//y-
                        interfaces++;
                    if(!DS_cube->get_spot(i, j, k-1) )//(int)(*(ds+(i-1)*info->sizeofx*info->sizeofy+(j-1)*info->sizeofx+(k-2)))==0)	//x-
                        interfaces++;
                    if(!DS_cube->get_spot(i+1, j, k) )//(int)(*(ds+  (i)*info->sizeofx*info->sizeofy+(j-1)*info->sizeofx+(k-1)))==0)	//z+
                        interfaces++;
                    if(!DS_cube->get_spot(i, j+1, k) )//(int)(*(ds+(i-1)*info->sizeofx*info->sizeofy+  (j)*info->sizeofx+(k-1)))==0)	//y+
                        interfaces++;
                    if(!DS_cube->get_spot(i, j, k+1) )//(int)(*(ds+(i-1)*info->sizeofx*info->sizeofy+(j-1)*info->sizeofx+  (k)))==0)	//x+
                        interfaces++;
				} 
			}
		}
		PoroDlg->pbar_simulation_progress.StepIt();

		if (PoroDlg->stop_request)
			return 2;
	}

//write results in .txt file
	double poros;
	unsigned long int totalVolume;
	totalVolume = (DS_cube->wX)*(DS_cube->wY)*(DS_cube->wZ);
	poros = (double)(pore)/(double)totalVolume;
	fp<<"Porosity "<<poros<<endl;

	double SSA;
	unsigned long int parsed_volume=((DS_cube->wX - 2)*(DS_cube->wY - 2)*(DS_cube->wZ - 2));
	SSA = (double)(interfaces)/((double)(parsed_volume)*pixsize*0.000001);
    fp<<"SSA "<<SSA<<" m2/m3"<<endl;
    
	SSA = (double)interfaces / ((1-poros) * (double)DENSITY_SSA * (double)(parsed_volume));
    fp<<"SSA "<<SSA<<" m2/g"<<endl;	

    fp<<"Total_Volume "<<totalVolume<<" pixelcube"<<endl;
    fp<<"Total_Pixel_Volume_Sampled "<<parsed_volume<<" pixelcube"<<endl;
	fp<<"Total_Volume_pore "<<pore<<" pixelcube"<<endl;
	fp<<"Total_Interface_Area "<<interfaces<<" pixelsquare"<<endl;
	fp<<"Pixel_Size "<<pixsize<<" microns"<<endl;
	fp<<"Threshold "<<threshold<<endl;

	//per slice results
	fp<<endl<<endl<<"slice#,porosity"<<endl;
	for (i = 0; i < DS_cube->wZ; i++)
		fp<<i<<","<< (float)poro_table[i] / (float)DS_cube->wX /(float)DS_cube->wY <<endl;

    fp.close();

	return 0;
}