#include "StdAfx.h"
#include "SphereGrowth.h"
#include "DistanceTransform.h"
//------------------------------------------------------------------------------------
// SphereGrowth::RunAnalysis
//
// Description - Outputs analysis results from Sphere run, if it is done
//
// Parameters - if 'outfile' is NULL, file_prefix+"_Sphere_Analysis.csv" is used. Otherwise
// 'outfile' is used for output file.
// 'dist_offset' is used to indicate to start the distribution at a specific line number
// note: has to be after all header information, otherwise error occurs
//------------------------------------------------------------------------------------
sg_status SphereGrowth::RunAnalysis(float pixel_size, sg_options options, int dist_offset, char *outfile /*= NULL*/)
{
char update_txt[512];
sprintf(update_txt, "Performing analysis on Sphere Map...");
progress_txt->SetWindowText(update_txt);
// Make sure we have information needed to run analysis
if ( !(cur_state & SG_STATE_GROWTH_DONE) )
{
return SG_ERR_INVALID_STATE;
}
// Check if option to stop after Sphere Growth is set
if ( cur_setting & SG_STOP_AFTER_SPHERE )
{
return SG_NO_ERROR;
}
// Variables
unsigned short z, y, x, i;
sg_sphr_type diam_min = (sg_sphr_type)1; // We can probably assume smallest diameter is 1
sg_sphr_type diam_max = max_diameter;
sg_sphr_type dist_size = diam_max - diam_min + 1;
char excel_out_file[512];
ofstream fout;
// Statistics variables
unsigned long long sphere_voxels = 0;
unsigned long long pore_voxels = 0;
unsigned int error_vox_fiber = 0;
unsigned int error_vox_pore = 0;
// Allocate memory for distribution array
unsigned long long *diam_dist = new unsigned long long [dist_size];
if ( diam_dist == NULL )
{
write_to_log("SphereGrowth::RunAnalysis - Out of memory for diameter distribution.");
return SG_ERR_MEMORY_ALLOC;
}
// Clear distribution array
memset( diam_dist, 0, dist_size*sizeof(unsigned long long) );
// Fill out diameter distribution.
// At the same time find number of sphere voxels that overlapped with
// fibers, and number of pores that were untouched by spheres
FOR_wZwYwX_LOOP(ds_cube, z, y, x)
{
// Check if spot is fiber, and has a Sphere diameter greater than 0.
// If so, then a sphere has overlapped where it shouldn't have.
if ( ds_cube->get_spot(z, y, x) )
{
if ( sphere_map->datac(z, y, x) > 0 )
{
++ diam_dist[ (sphere_map->datac(z, y, x) - diam_min) ];
++ sphere_voxels;
++ error_vox_fiber;
}
}
else // Pore
{
// Check if diameter is zero, in which case a valid spot did not
// get filled by a grown sphere
if ( sphere_map->datac(z, y, x) == (sg_sphr_type) 0)
{
++ error_vox_pore;
}
else // Everything valid, fill in distribution.
{
++ diam_dist[ (sphere_map->datac(z, y, x) - diam_min) ];
++ sphere_voxels;
}
++ pore_voxels;
}
}
// Determine appropriate file name
if ( outfile )
{
sprintf(excel_out_file, "%s", outfile);
}
else
{
sprintf(excel_out_file, "%s_Sphere_Analysis.csv", file_prefix);
}
// Check if file should be opened in append mode, or normal mode
if ( options & SG_APPEND_RESULTS )
{ // Append mode
fout.open(excel_out_file, ios::out | ios::app);
}
else
{ // Normal mode
fout.open(excel_out_file, ios::out);
}
// Make sure file is open
if( !fout.is_open() )
{
write_to_log("SphereGrowth::RunAnalysis - Error opening file %s.", excel_out_file);
return SG_ERR_FILE_IO;
}
// Add header information if requested
if ( !(options & SG_DO_NOT_ADD_HEADER) )
{ // add in future
}
OutputSphereHeader(fout);
// Print out total voxel and error information.
fout << "Voxel Counters:" << endl;
fout << "Sphere Voxels,," << "Pore Voxels,," << "Fiber Overlap,," << "Missed Pores" << endl;
fout << sphere_voxels << ",," << pore_voxels << ",," << error_vox_fiber << ",," << error_vox_pore << endl << endl;
// Volume ratios
long long total_voxels = ds_cube->wX*ds_cube->wY*ds_cube->wZ;
fout << "Volume Ratios:" << endl;
fout << "Sphere To Total,," << "Pores to Total (porosity)" << endl;
fout << (double)sphere_voxels/total_voxels << ",," << (double)pore_voxels/total_voxels << endl << endl;
// Determine mean
double mean = 0.0;
for(i = 0; i < dist_size; i++)
{
mean = mean + i*diam_dist[i];
}
mean = mean/sphere_voxels;
fout << "mean," << (diam_min + mean)*pixel_size << ",(microns),with average #," << \
( diam_dist[(int)mean] + diam_dist[(int)mean+1] )/2 << endl;
// Median
unsigned long long median = 0;
for(i = 0; i < dist_size; i++)
{
median = median + diam_dist[i];
if ( (median*2) > sphere_voxels)
{
break;
}
}
fout << "av median," << (float)(diam_min + i - 0.5)*pixel_size << \
",(microns),with average #," << (diam_dist[i] + diam_dist[i-1])/2 << endl;
// Most common diameter
unsigned long long max= 0 ;
unsigned char maxi = 0;
for(i = 0; i < dist_size; i++)
{
if (max < diam_dist[i])
{
max = diam_dist[i];
maxi = (unsigned char) i;
}
}
fout << "max," << (float)(maxi + diam_min)*pixel_size << \
",(microns),with #," << diam_dist[maxi] << endl;
// If requested, write out empty lines before distribution. This is used to ensure that
// the distribution always starts on a particular line for easy copy/past
if ( dist_offset > 0 )
{
fout << endl << "***** See line #" << dist_offset << " for distribution. *****" << endl;
// Need to close file first
fout.close();
if ( OpenFileForAppendAtLine(dist_offset, excel_out_file, fout) != 0 )
{
write_to_log("SphereGrowth::RunAnalysis - Error reopening file '%s' after distribution offset.", excel_out_file);
return SG_ERR_FILE_IO;
}
}
// Distribution
fout << "Diameters (microns),Probability,CDF,Normalized Distribution(% per unit bucket)" << endl;
if ( (cur_setting & SG_ODD_DIAMETER_ONLY) && (cur_setting & SG_PAD_ODD_DISTRIBUTION) )
{ // If pad odd is enabled, cut distribution values in half and transfer other half to next
// even diameter. This helps comparison in with an odd+even graph.
for(i = 0; i < dist_size; i++)
{
if ( (i%2) == 0 )
diam_dist[i] /= 2;
else
diam_dist[i] = diam_dist[i-1];
}
}
// Write out to file % values
double cdf = 0.0;
for(i = 0; i < dist_size; i++)
{
cdf = cdf + (double)diam_dist[i]/(double)sphere_voxels;
// Skip every 'even' spot if ODD diameter only is valid
if ( (!(cur_setting & SG_ODD_DIAMETER_ONLY)) || ((i%2) == 0) || (cur_setting & SG_PAD_ODD_DISTRIBUTION) )
{
fout << (double)(diam_min + i)*pixel_size << ","; // Diameter in micrometers
fout << (double)diam_dist[i]/(double)sphere_voxels << ","; // Probability
fout << cdf << ","; // CDF
fout << (double)diam_dist[i]/(double)sphere_voxels/(double)pixel_size << endl; // Normalized distribution
}
}
// Clean up
fout.close();
delete[] diam_dist;
return SG_NO_ERROR;
}
//------------------------------------------------------------------------------------
// SphereGrowth::OutputSphereHeader
// Description - Outputs run options to a .csv formatted file stream
//------------------------------------------------------------------------------------
sg_status SphereGrowth::OutputSphereHeader(ofstream &file_ptr)
{
// Make sure file is open
if( !file_ptr.is_open() )
{
write_to_log("SphereGrowth::OutputSphereHeader - Called with a closed file.");
return SG_ERR_FILE_IO;
}
// Print out run options that were selected
file_ptr << "Sphere settings:" << endl;
// Odd diameter only, or even as well
if ( cur_setting&SG_ODD_DIAMETER_ONLY )
{
file_ptr << ",Odd diameter only." << endl;
}
else
{
file_ptr << ",Odd+Even diameters." << endl;
}
// Boundary or mirror for boarder
if ( cur_setting&SG_BOUNDARY_IS_SOLID )
{
file_ptr << ",Samples boundaries were solid." << endl;
}
else
{
file_ptr << ",Samples boundaries were mirrored." << endl;
}
// Conserve memory option
if ( cur_setting&SG_CONSERVE_MEMORY )
{
file_ptr << ",RAM was conserved." << endl;
}
else
{
file_ptr << ",Full RAM was used." << endl;
}
// Distance transform
if ( cur_setting&SG_QUICK_DIST_XFRM )
{
file_ptr << ",Distance Transform was used with F-E-C values (" <<
(int)cur_FEC.face << "-" << (int)cur_FEC.edge << "-" << (int)cur_FEC.corner << ")." << endl;
}
else
{
file_ptr << ",Sphere fitting method was used." << endl;
}
// Indicate what format cubes were saved in
file_ptr << ",The following stacks were output: ";
if ( (cur_setting&(SG_SAVE_DIST_GRAY|SG_SAVE_DIST_RGB)) )
{
file_ptr << "Diameter Map ( " << ((cur_setting&SG_SAVE_DIST_GRAY) ? ("RAW "):("")) <<
((cur_setting&SG_SAVE_DIST_RGB) ? ("RGB "):("")) << ") ";
}
if ( (cur_setting&(SG_SAVE_GROWN_GRAY|SG_SAVE_GROWN_RGB)) )
{
file_ptr << "Sphere Map ( " << ((cur_setting&SG_SAVE_GROWN_GRAY) ? ("RAW "):("")) <<
((cur_setting&SG_SAVE_GROWN_RGB) ? ("RGB "):("")) << ")";
}
if ( !(cur_setting&(SG_SAVE_DIST_GRAY|SG_SAVE_DIST_RGB|SG_SAVE_GROWN_GRAY|SG_SAVE_GROWN_RGB)) )
{
file_ptr << "(none)";
}
file_ptr << endl << "___________________________________________________" << endl;
return SG_NO_ERROR;
}