MA_vc_Main.cpp #2

#include "StdAfx.h"
#include "MA_Voxel_Coding.h"
#include <math.h>



// ----------------------------------------------------------------------------------------
// MA_VoxelCoding
//
//	Description -
// ----------------------------------------------------------------------------------------
MA_VoxelCoding::MA_VoxelCoding(DS_1b_cube *ds_cube_ptr, CEdit *prog_text, CProgressCtrl *prog_bar)
{
	queue_mark_coord.from_zyx(0xFFFF,0xFFFF,0xFFFF);
	run_settings = 0;
	// Set default pixel size and cube sizes
	pixel_size_in_um = 0.7f;
	node_distribution_cube_size_in_um = 35.0f;

	// Other parameter defaults
	delete_sp_below_size = 20;
	merge_nodes_bs_multiplier = 2.0f;

	memory_limit_total_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_bs_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_ss_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_cl_bytes = INFINITE_MEMORY_64BIT;

	BS_cube = NULL;
	SS_cube = NULL;
	CL_cube = NULL;
	
	for ( unsigned short int i = 0; i < TOTAL_MARK_CUBES; i++ )
	{
		MARK_cubes[i] = NULL;
	}

	// Set external pointers
	boundary_cube = ds_cube_ptr;
	progress_win_ptr = prog_text;
	progress_bar_ptr = prog_bar;
	current_recursion = 0;

	// Allocate items that do not use much memory, and remain the same over multiple runs.
	BS_Xfrm = new DistanceXfrm<bs_t>( boundary_cube, progress_bar_ptr );
	SS_Xfrm = new DistanceXfrm<ss_t>( boundary_cube );
	SS_Xfrm->changeFECValues( 1, 2, 3 );

	ss_infinity = SS_Xfrm->getInfinityValue();
	ss_medial_path = ss_infinity - 1;

#if MA_DO_FULL_TRACE
	// This can be used to enable/disable full trace on the fly during debugging
	full_trace_enabled = false;
	trace_text_buffer[0] = 0;
 #if MA_DISABLE_LOG_FLUSH
	disable_log_flush();
 #endif
#endif
}

MA_VoxelCoding::~MA_VoxelCoding(void)
{
	delete	BS_Xfrm;
	delete	SS_Xfrm;

#if MA_DISABLE_LOG_FLUSH
	enable_log_flush();
#endif
}


// --------------------------------------------------------------------------------------------
// 
//	Description -
// --------------------------------------------------------------------------------------------
ma_status	MA_VoxelCoding::runMedialAxis( ma_options run_opts, char *output_prefix )
{
	result_file_prefix = output_prefix;
	run_options = run_opts;

	// Allocate memory and initialize other variables
	prepareForRun();

	// ----------------------------------------------------------------------------------------
	// BS Distance transform
	if ( !(run_settings & BS_IS_PRELOADED) )
	{
		progress_win_ptr->SetWindowText("Now creating distance transform (BS Field)...");

		// Initialize FEC values
		BS_Xfrm->changeFECValues( 3, 4, 5 );

		// In case caching is used, transfer all memory to BS cube for the transform
		CL_cube->resizeCacheBuffer( 0 );
		SS_cube->resizeCacheBuffer( 0 );
		BS_cube->resizeCacheBuffer( memory_limit_total_bytes );
		
		// Do requested type of BS distance transform (solid or mirrored boundary)
		if ( run_options & MA_BOUNDARY_IS_SOLID )
		{
			write_to_log("MA_VoxelCoding::RunMedialAxis - Creating BS cube using walls as boundaries.");
			BS_Xfrm->runBS( BS_cube, DT_INIT_INPUT | DT_BOUNDARY_IS_SOLID );
		}
		else
		{
			write_to_log("MA_VoxelCoding::RunMedialAxis - Creating BS cube using mirrors as boundaries.");
			BS_Xfrm->runBS( BS_cube, DT_INIT_INPUT );
		}

		// Save BS cube if requested.	
		writeToFileBS( output_prefix, run_options );
	}

	// Skip SS and CL cube generation if MA cube is preloaded
	if ( !(run_settings & MA_IS_PRELOADED) )
	{
		// ----------------------------------------------------------------------------------------
		// SS Transform over entire sample
		if ( !(run_settings & SS_IS_PRELOADED) )
		{
			// Attempt to generate all SS fields, *progress kept within function*
			generateFullSS();

			// Save SS cube if requested.
			writeToFileSS( output_prefix, run_options );
		}


		// ----------------------------------------------------------------------------------------
		// Cluster cube and graph generation
		if ( !(run_settings & CL_IS_PRELOADED) )
		{
			// Create Cluster cube
			generateClusterCube();
		
			// Save CL cube (or variant) if requested.
			writeToFileCL( output_prefix, run_options );
		}

		// Now generate a graph structure using the CL cube, *progress kept within function*
		generateClusterGraph();


		// ----------------------------------------------------------------------------------------
		// Shortest paths and medial axis
		traceShortestPaths();

		// Save to file if requested
		writeToFileShortestPaths( output_prefix, run_options );
		traceMedialAxis();

		// Mark the medial axis in FINAL_MA cube, this is what the network graph is generated from
		MARK_cubes[FINAL_MA]->setall();
		for ( unsigned int j = 0; j < number_of_volumes; j++ )
		{
			for ( unsigned int i = 0; i < medial_path_volumes[j].number_of_medial_paths; i++ )
			{
				markPathInCube( medial_path_volumes[j].medial_paths_array[i], MARK_cubes[FINAL_MA], false );
			}
		}

		writeToFileMedialPaths( output_prefix, run_options );
	}

	// ----------------------------------------------------------------------------------------
	// Run analysis (no longer need CL cube or graph)
	delete	CL_cube;
	CL_cube = NULL;
	delete[]	cluster_nodes_array;
	cluster_nodes_array = NULL;

	MARK_cubes[TRAVELED]->clear();
	performAnalysis( output_prefix );

	// ----------------------------------------------------------------------------------------
	// Clean up and return
	cleanUpAfterRun();

	return	MA_NO_ERRORS;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
ma_status	MA_VoxelCoding::preloadBSCube( char *file_name )
{
	progress_win_ptr->SetWindowText("Loading BS cube...");
	write_to_log("MA_VoxelCoding::preloadBSCube - Preloading BS cube from %s.", file_name);

	BS_cube = new DSt_cube<bs_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_bs_bytes);

	if ( BS_cube->readFromFile(file_name) == DST_NO_ERROR )
	{
		run_settings |= BS_IS_PRELOADED;
		return	MA_NO_ERRORS;
	}

	delete	BS_cube;
	BS_cube = NULL;
	write_to_log("MA_VoxelCoding::preloadBSCube - Error occured during loading of file.", file_name);
	return	MA_FILE_OPEN_ERROR;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
ma_status	MA_VoxelCoding::preloadSSCube( char *file_name )
{
	progress_win_ptr->SetWindowText("Loading SS cube...");
	write_to_log("MA_VoxelCoding::preloadSSCube - Preloading SS cube from %s.", file_name);

	SS_cube = new DSt_cube<ss_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_ss_bytes);

	if ( SS_cube->readFromFile(file_name) == DST_NO_ERROR )
	{
		run_settings |= SS_IS_PRELOADED;
		return	MA_NO_ERRORS;
	}

	delete	SS_cube;
	SS_cube = NULL;
	write_to_log("MA_VoxelCoding::preloadSSCube - Error occured during loading of file.", file_name);
	return	MA_FILE_OPEN_ERROR;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
ma_status	MA_VoxelCoding::preloadCLCube( char *file_name )
{
	progress_win_ptr->SetWindowText("Loading CL cube...");
	write_to_log("MA_VoxelCoding::preloadCLCube - Preloading CL cube from %s.", file_name);

	CL_cube = new DSt_cube<cl_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_total_bytes);

	if ( CL_cube->readFromFile(file_name) == DST_NO_ERROR )
	{
		Vect3usi	cur_coord;

		number_of_clusters = 0;
	    
		progress_bar_ptr->SetRange(0, boundary_cube->wZ-1);
		progress_bar_ptr->SetPos(0);
		progress_bar_ptr->SetStep(1);

		progress_win_ptr->SetWindowText("Determining number of clusters...");

		// Find largest cluster (same as number of clusters)
		for ( cur_coord.z = 0;  cur_coord.z < boundary_cube->wZ;  cur_coord.z++)
		{
			for ( cur_coord.y = 0;  cur_coord.y < boundary_cube->wY;  cur_coord.y++)
				for ( cur_coord.x = 0;  cur_coord.x < boundary_cube->wX;  cur_coord.x++)
				{
					// Check if this is a higher numbered cluster (note: cluster IDs
					// start from '1', so no need to add one.
					if ( CL_cube->datac(cur_coord) > number_of_clusters )
					{
						number_of_clusters = CL_cube->datac(cur_coord);
					}
				}

			progress_bar_ptr->StepIt();
		}

		write_to_log("MA_VoxelCoding::preloadCLCube - Found %d clusters.", number_of_clusters);
		run_settings |= CL_IS_PRELOADED;
		return	MA_NO_ERRORS;
	}

	delete	CL_cube;
	CL_cube = NULL;
	write_to_log("MA_VoxelCoding::preloadCLCube - Error occured during loading of file.", file_name);
	return	MA_FILE_OPEN_ERROR;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
ma_status	MA_VoxelCoding::preloadMACube( char *file_name )
{	
	progress_win_ptr->SetWindowText("Loading MA cube...");
	write_to_log("MA_VoxelCoding::preloadMACube - Preloading MA cube from %s.", file_name);

	for ( unsigned short int i = 0; i < TOTAL_MARK_CUBES; i++ )
	{
		MARK_cubes[i] = new DS_1b_cube(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX);
	}

	if ( MARK_cubes[FINAL_MA]->readFromFile(file_name) )
	{
		// When MA_cube is written out, Medial axis is '1', but needs to be other way around
		// for analysis so invert.
		MARK_cubes[FINAL_MA]->invert();
		run_settings |= MA_IS_PRELOADED;
		return	MA_NO_ERRORS;
	}

	for ( unsigned short int i = 0; i < TOTAL_MARK_CUBES; i++ )
	{
		delete	MARK_cubes[i];
	}
	write_to_log("MA_VoxelCoding::preloadMACube - Error occured during loading of file.", file_name);
	return	MA_FILE_OPEN_ERROR;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
void	MA_VoxelCoding::enableCachingFromDisc( unsigned int memory_limit_in_mb )
{
	write_to_log("MA_VoxelCoding::enableCachingFromDisc - Caching from disc enabled, limit %d MB.", memory_limit_in_mb);

	// Determine memory limit in bytes (only used if MA_ENABLE_DISC_CACHE is enabled)
	memory_limit_total_bytes = (UINT64)memory_limit_in_mb * 1048576;

	unsigned int	total_bytes_per_spot = sizeof(bs_t)+sizeof(ss_t)+sizeof(cl_t);
	memory_limit_bs_bytes = memory_limit_total_bytes * sizeof(bs_t)/total_bytes_per_spot;
	memory_limit_ss_bytes = memory_limit_total_bytes * sizeof(ss_t)/total_bytes_per_spot;
	memory_limit_cl_bytes = memory_limit_total_bytes * sizeof(cl_t)/total_bytes_per_spot;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
void	MA_VoxelCoding::changeNodeDistributionCubeSize( float size_in_um )
{
	write_to_log("MA_VoxelCoding::changeNodeDistributionCubeSize - Node distribution cube size changed to %f micrometers.", size_in_um);
	node_distribution_cube_size_in_um = size_in_um;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
void	MA_VoxelCoding::changePixelSize( float size_in_um )
{
	write_to_log("MA_VoxelCoding::changePixelSize - Pixel size changed to %f micrometers.", size_in_um);
	pixel_size_in_um = size_in_um;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
void	MA_VoxelCoding::prepareForRun( )
{
	// Allocate required cubes (note by default cache is disabled)
	if ( !(run_settings & BS_IS_PRELOADED) )
	{
		BS_cube = new DSt_cube<bs_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_bs_bytes);
	}
	if ( !(run_settings & SS_IS_PRELOADED) )
	{
		SS_cube = new DSt_cube<ss_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_ss_bytes);
	}
	// Do not need CL cube at all if medial axis is preloaded
	if ( !(run_settings & CL_IS_PRELOADED) && !(run_settings & MA_IS_PRELOADED) )
	{
		number_of_clusters = 0;
		CL_cube = new DSt_cube<cl_t>(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX, CACHE_SUB_ROW, memory_limit_cl_bytes);
	}

	// Allocate binary cubes
	if ( !(run_settings & MA_IS_PRELOADED) )
	{
		for ( unsigned short int i = 0; i < TOTAL_MARK_CUBES; i++ )
		{
			MARK_cubes[i] = new DS_1b_cube(boundary_cube->wZ, boundary_cube->wY, boundary_cube->wX);
		}
	}

	// Smoothing will use '1' for pores
	SS_Xfrm_Smooth = new DistanceXfrm<ss_t>( MARK_cubes[SMOOTHING] );
	SS_Xfrm_Smooth->changeFECValues( 1, 2, 3 );
	SS_Xfrm_Smooth->changePoreValue( true );
	SS_Xfrm_Smooth->changeEndValueForSS( ss_medial_path - 2 );

	SS_Xfrm_Axis = new DistanceXfrm<ss_t>( MARK_cubes[FINAL_MA] );
	SS_Xfrm_Axis->changeFECValues( 1, 2, 3 );

	// Multi use group of queues of type Vect3usi, preallocated.  If a function uses this queue,
	// it should be emptied at finish.  Each function should find the first group with an empty list
	// to avoid hard to debug conflicts.  First group is id '0', each additional one being incremented.
	recursion_queue = new GroupedLists<int,Vect3usi>( 0, MA_PREALLOC_GENERAL_QUEUE );

	// Used in trace medpaths for backing up path coordinates when walling off clusters, and during analysis
	coordinate_queue = new LinkedList<Vect3usi>( 40 );

	// Used in analysis graph generation
	graph_node_queue = new LinkedList<NetworkGraphNode*>;

	// Used in Link functions for temporary max cluster point storage
	max_points_list = new LinkedList<Vect3usi>( 40 );
	clusters_max_points_lists = new LinkedList<LinkedList<Vect3usi>>(50);

	// These do not get allocated until we know more information.
	medial_path_volumes = NULL;
	cluster_nodes_array = NULL;

	// Initialize statistical counters
	number_of_volumes = 0;
	number_of_medial_paths_total = 0;
	number_of_bs_maximums = 0;
	number_of_bs_maximums_removed = 0;
	number_of_bs_maximums_removed_s2 = 0;

	number_of_local_max_clusters = 0;
	number_of_merging_clusters = 0;
	number_of_merging_clusters_removed = 0;
	number_of_dividing_clusters = 0;
	number_of_dividing_clusters_removed = 0;
	number_of_rp_clusters = 0;

	number_of_failed_ma_links = 0;
	number_of_failed_ma_path_links = 0;
}


// ----------------------------------------------------------------------------------------
// 
//	Description -
// ----------------------------------------------------------------------------------------
void	MA_VoxelCoding::cleanUpAfterRun( )
{
	run_settings = 0;
	memory_limit_total_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_bs_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_ss_bytes = INFINITE_MEMORY_64BIT;
	memory_limit_cl_bytes = INFINITE_MEMORY_64BIT;

	if ( BS_cube )
	{
		delete	BS_cube;
		BS_cube = NULL;
	}
	if ( CL_cube )
	{
		delete	CL_cube;
		CL_cube = NULL;
	}
	if ( SS_cube )
	{
		delete	SS_cube;
		SS_cube = NULL;
	}

	for ( unsigned short int i = 0; i < TOTAL_MARK_CUBES; i++ )
	{
		delete	MARK_cubes[i];
	}

	delete	SS_Xfrm_Smooth;
	delete	SS_Xfrm_Axis;

	delete	recursion_queue;
	delete	coordinate_queue;

	delete	graph_node_queue;

	delete	max_points_list;
	delete	clusters_max_points_lists;

	if ( medial_path_volumes )
	{
		delete[]	medial_path_volumes;
	}

	if ( cluster_nodes_array )
	{
		delete[]	cluster_nodes_array;
	}
}