#include "StdAfx.h"
#include "SphereGrowth.h"
#include "DistanceTransform.h"
using namespace std;
// -----------------------------------------------------------------
// GrowAllFromMinToMaxDiameter:
// Description - Grows spheres from a diameter of 1 up to 'max_diameter', one sphere size per pass.
// This results in larger spheres overwriting smaller spheres.
// Parameters - 'slice_to_snap' index of slice to take a snapshot of after each diameter, only useful if
// SG_SAVE_STAGE_CUBE is specified 'in cur_setting'. Saves this cube to file at end of this
// function.
// -----------------------------------------------------------------
sg_status SphereGrowth::GrowAllFromMinToMaxDiameter(unsigned short slice_to_snap/*=0*/)
{
// Make sure state is valid... Xfrm cube should be filled out
if ( !(cur_state & SG_STATE_XFRM_DONE) )
{
return SG_ERR_INVALID_STATE;
}
// Update state
cur_state &= (~SG_STATE_GROWTH_DONE);
char temp_file_name[512];
char text[512];
bool neigh[3][3][3];
sg_sphr_type cur_diam = 1, prev_max_diam;
unsigned short z, y, x;
unsigned short z_cube; // Used for flexibility in using a buffer or entire xfrm cube
char lz, uz, ly, uy, lx, ux, dx, dy, dz;
dst_status return_val;
// Reset progress control
int prog_inc_val = max_diameter/100 + 1;
int prog_next_update = prog_inc_val;
int prog_cur = 0;
progress_bar->SetPos(0);
progress_bar->SetStep(1);
progress_bar->SetRange(0, ((max_diameter>99)?(99):(max_diameter)) );
// If conserve memory option is set, we want to write out the xfrm cube to file and only
// use a 3 slices at a time
if ( cur_setting & SG_CONSERVE_MEMORY )
{
// Invalidate xfrm cube
cur_state &= (~SG_STATE_XFRM_DONE);
// Write out diameter map to temp file
sprintf(temp_file_name, "%s__TEMP__xfrm_buff", output_prefix);
if ( distance_xfrm->writeToFile(temp_file_name, DST_UPDATE_FILE_NAME) != DST_NO_ERROR )
{
return SG_ERR_FILE_IO;
}
// Can now delete Xfrm cube, and create a buffer of only 3 slices
delete distance_xfrm;
distance_xfrm = new DSt_cube<sg_xfrm_type>(3, ds_cube->wY, ds_cube->wX);
// Also need to create Sphere map now that we have free room
sphere_map = new DSt_cube<sg_sphr_type>(ds_cube->wZ, ds_cube->wY, ds_cube->wX);
if (sphere_map == NULL)
{ // Error occured during memory allocation
write_to_log("SphereGrowth::GrowAllFromMinToMaxDiameter - Error allocating 'sphere_map' cube of size (%d, %d, %d).", ds_cube->wX, ds_cube->wY, ds_cube->wZ);
return SG_ERR_MEMORY_ALLOC;
}
}
// Zero out sphere map, if even cubes are used the even_flag cube was cleared during memory setup
sphere_map->clear();
// Create cube for saving a slice after each growth stage
DSt_cube<sg_sphr_type> *snapshot = NULL;
unsigned short cur_snap_indx = 1;
if ( cur_setting & (SG_SAVE_STAGE_CUBE_GREY | SG_SAVE_STAGE_CUBE_RGB) )
{
// If odd diameter only, we need half the snapshot cube
if (cur_setting & SG_ODD_DIAMETER_ONLY)
snapshot = new DSt_cube<sg_xfrm_type>( max_diameter/2 + 2, ds_cube->wY, ds_cube->wX );
else
snapshot = new DSt_cube<sg_xfrm_type>( max_diameter + 1, ds_cube->wY, ds_cube->wX );
snapshot->clear();
}
// Now start growing Spheres, from smallest diameter until max diameter
while(cur_diam <= max_diameter)
{
// If we are only doing Odd diameters, we can skip all even diameters
if ( !((cur_setting & SG_ODD_DIAMETER_ONLY) && ((cur_diam%2) == 0)) )
{
sprintf(text, "Now spreading spheres with a diameter of %d (out of %d)...", cur_diam, max_diameter);
progress_txt->SetWindowText(text);
// Find every spot matching cur_diam, and grow a sphere there
for (z = 0; z < ds_cube->wZ; z++)
{
z_cube = z;
if ( cur_setting & SG_CONSERVE_MEMORY )
{
z_cube = 1;
// Read surrounding slices into memory buffer
if (z == 0)
{// For first slice, we need to offset by one, and make z_cube point to '0'
z_cube = 0;
return_val = distance_xfrm->readFromFile(temp_file_name, 0);
}
else if (z == (ds_cube->wZ-1))
{// For last slice, we need to offset by negative two, and make z_cube point to '2'
z_cube = 2;
return_val = distance_xfrm->readFromFile(temp_file_name, z-2);
}
else
{// Everything in between, we can just read from a negative one offset
return_val = distance_xfrm->readFromFile(temp_file_name, z-1);
}
if ( return_val != DST_NO_ERROR )
{ // Make sure there was no file reading error
return SG_ERR_FILE_IO;
}
}
// Set up Z limits (if on boundary)
lz = -1;
uz = 1;
if ( z == 0 )
{
lz = 0;
for (dy = 0; dy < 3; dy++)
for (dx = 0; dx < 3; dx++)
neigh[0][dy][dx] = false;
}
if ( z == (ds_cube->wZ-1) )
{
uz = 0;
for (dy = 0; dy < 3; dy++)
for (dx = 0; dx < 3; dx++)
neigh[2][dy][dx] = false;
}
for (y = 0; y < ds_cube->wY; y++)
{
// Set up Y limits (if on boundary)
ly = -1;
uy = 1;
if ( y == 0 )
{
ly = 0;
for (dz = 0; dz < 3; dz++)
for (dx = 0; dx < 3; dx++)
neigh[dz][0][dx] = false;
}
if ( y == (ds_cube->wY-1) )
{
uy = 0;
for (dz = 0; dz < 3; dz++)
for (dx = 0; dx < 3; dx++)
neigh[dz][2][dx] = false;
}
for (x = 0; x < ds_cube->wX; x++)
{
// Make sure current spot has the diameter we are currently growing
if (distance_xfrm->datac(z_cube, y, x) == cur_diam)
{
// Set up X limits (if on boundary)
lx = -1;
ux = 1;
if ( x == 0 )
{
lx = 0;
for (dz = 0; dz < 3; dz++)
for (dy = 0; dy < 3; dy++)
neigh[dz][dy][0] = false;
}
if ( x == (ds_cube->wX-1) )
{
ux = 0;
for (dz = 0; dz < 3; dz++)
for (dy = 0; dy < 3; dy++)
neigh[dz][dy][2] = false;
}
// Find the maximum neighboring radius in 6 directions
prev_max_diam = 0;
if (x > 0 && (prev_max_diam < distance_xfrm->datac(z_cube, y, x-1)))
prev_max_diam = distance_xfrm->datac(z_cube, y, x-1);
if (y > 0 && (prev_max_diam < distance_xfrm->datac(z_cube, y-1, x)))
prev_max_diam = distance_xfrm->datac(z_cube, y-1, x);
if (z > 0 && (prev_max_diam < distance_xfrm->datac(z_cube-1, y, x)))
prev_max_diam = distance_xfrm->datac(z_cube-1, y, x);
// For forward going
if ((x+1 < ds_cube->wX) && (prev_max_diam < distance_xfrm->datac(z_cube, y, x+1)) && (cur_diam < distance_xfrm->datac(z_cube, y, x+1)))
prev_max_diam = distance_xfrm->datac(z_cube, y, x+1);
if ((y+1 < ds_cube->wY) && (prev_max_diam < distance_xfrm->datac(z_cube, y+1, x)) && (cur_diam < distance_xfrm->datac(z_cube, y+1, x)))
prev_max_diam = distance_xfrm->datac(z_cube, y+1, x);
if ((z+1 < ds_cube->wZ) && (prev_max_diam < distance_xfrm->datac(z_cube+1, y, x)) && (cur_diam < distance_xfrm->datac(z_cube+1, y, x)))
prev_max_diam = distance_xfrm->datac(z_cube+1, y, x);
// We have three cases:
// 1) diam+1 < neighbor diam... then do nothing, cause neighbor will cover this sphere in full
// 2) diam+1 == neighbor diam, or
// diam = neighbor diam... then we spread radius from diam-2 to diam
// 3) diam > neighbor diam... then we should spread radius across entire sphere
// Case 1, we don't do anything!
// Case 2
if ( (cur_diam == prev_max_diam) || ((cur_diam+1) == prev_max_diam) )
{
if (cur_diam > 1)
GrowSingleSphere(z, y, x, cur_diam-2, cur_diam);
else
GrowSingleSphere(z, y, x, 0, cur_diam);
}
// Case 3
else if (cur_diam > prev_max_diam)
{
GrowSingleSphere(z, y, x, 0, cur_diam);
}
}//End IF CURRENT DIAMETER
}//End X loop
}//End Y loop
}//End Z loop
// Save snapshot slice to cube
if ( cur_setting & (SG_SAVE_STAGE_CUBE_GREY | SG_SAVE_STAGE_CUBE_RGB) )
{
// Copy entire sphere cube slice to snapshot cube
// WARNING: this will not work with a file cached data cube!!!
memcpy( &snapshot->data[cur_snap_indx][0][0], \
&sphere_map->data[slice_to_snap][0][0], \
sizeof(sg_sphr_type) * ds_cube->wX * ds_cube->wY );
// Increment index
++cur_snap_indx;
}
} // END if odd only skip even
// Go on to next diameter size
++cur_diam;
// Update progress bar when needed
if ( ++prog_cur == prog_next_update )
{
prog_next_update += prog_inc_val;
progress_bar->StepIt();
}
}
if ( cur_setting & SG_CONSERVE_MEMORY )
{
remove( temp_file_name );
}
// Save snapshot cube to file
if ( cur_setting & (SG_SAVE_STAGE_CUBE_GREY | SG_SAVE_STAGE_CUBE_RGB) )
{
char file_str[512];
// Before outputting to file, we want to mark all of the slices to contain all white
// (i.e above max value) for spots that are pores and have a snapshot value of 0.
// This will make it so the snapshot cube is white where no sphere has been grown.
FOR_wYwX_LOOP(ds_cube, y, x)
{
if ( !ds_cube->get_spot(slice_to_snap, y, x) )
{ // Spot is a pore, go through entire snapshot cubes Z index and mark any
// spots with a value of 0 as 'max_diamter' + 1
for (z = 0; z < snapshot->wZ; z++)
{
if ( snapshot->data[z][y][x] == 0 )
snapshot->data[z][y][x] = max_diameter+1;
}
}
}
if ( cur_setting & SG_SAVE_STAGE_CUBE_GREY )
{ // Save as grey
sprintf(file_str, "%s_Sphere_Map_Snapshots", output_prefix);
if ( snapshot->writeToFile(file_str, DST_UPDATE_FILE_NAME) != DST_NO_ERROR )
{ // This is not a fatal error, simply log it
write_to_log("SphereGrowth::GrowAllFromMinToMaxDiameter - Error outputting sphere map RAW snapshots to %s.", file_str);
}
}
if ( cur_setting & SG_SAVE_STAGE_CUBE_RGB )
{ // Save as RGB
sprintf(file_str, "%s_Sphere_Map_Snapshots", output_prefix);
if ( snapshot->writeToFileRGB(file_str, max_diameter, DST_UPDATE_FILE_NAME) != DST_NO_ERROR )
{ // This is not a fatal error, simply log it
write_to_log("SphereGrowth::GrowAllFromMinToMaxDiameter - Error outputting sphere map RGB snapshots to %s.", file_str);
}
}
delete snapshot;
} // END save snapshot cube
// Update state
cur_state |= SG_STATE_GROWTH_DONE;
return SG_NO_ERROR;
}
// -----------------------------------------------------------------
// GrowSingleSphere:
// Description - Will grow a sphere shell with an inner diameter of 'inn_diam' and an outer diameter
// of 'out_diam' at with center at location (z, y, x). If 'inn_diam' is 0, full sphere
// will be grown, otherwise just a shell.
// -----------------------------------------------------------------
void SphereGrowth::GrowSingleSphere(unsigned short z, unsigned short y, unsigned short x,
sg_sphr_type inn_diam, sg_sphr_type out_diam)
{
short int dz, dy;
int z_m_dz, z_p_dz, y_m_dy, y_p_dy, x_m_dx, x_p_dx;
short int dx, dx_inner_lim, dx_outer_lim;
sg_sphr_type rad = (out_diam)/2;
sg_sphr_type inn_rad = (inn_diam)/2;
unsigned int rad2;
unsigned int inn_rad2;
unsigned int dz2, dz2_dy2;
if ( (out_diam%2) == 0 )
{
rad2 = rad*rad;
}
else
{
rad2 = (unsigned int)(((double)rad+0.5f)*((double)rad+0.5f));
// inn_rad2 = (unsigned int)(((double)inn_rad+0.5f)*((double)inn_rad+0.5f));
}
inn_rad2 = inn_rad*inn_rad;
// Set initial spot to final diameter if spreading entire sphere
if ( inn_diam == 0 )
{
sphere_map->datac(z, y, x, out_diam);
}
for (dz = 0; dz <= rad; dz++)
{
dz2 = dz*dz;
z_p_dz = z+dz;
z_m_dz = z-dz;
for (dy = 0; dy <= GetSqrt(rad2 - dz2); dy++)
{
dz2_dy2 = dz2 + dy*dy;
y_p_dy = y+dy;
y_m_dy = y-dy;
dx_outer_lim = GetSqrt(rad2 - dz2_dy2);
if (dz2_dy2 > inn_rad2)
dx_inner_lim = -1;
else
dx_inner_lim = GetSqrt(inn_rad2 - dz2_dy2);
for (dx = dx_outer_lim; dx > dx_inner_lim; dx--)
{
x_p_dx = x+dx;
x_m_dx = x-dx;
//pos z, pos y
if (z_p_dz < sphere_map->wZ)
{
if (y_p_dy < sphere_map->wY)
{
if (x_p_dx < sphere_map->wX)
sphere_map->datac((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx, out_diam);
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) )
sphere_map->datac((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx, out_diam);
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if (x_p_dx < sphere_map->wX)
sphere_map->datac((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx, out_diam);
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) )
sphere_map->datac((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx, out_diam);
}
}
//neg z
if ( (z_m_dz != z_p_dz) && (z_m_dz >= 0) )
{//pos y
if (y_p_dy < sphere_map->wY)
{
if (x_p_dx < sphere_map->wX)
sphere_map->datac((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx, out_diam);
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) )
sphere_map->datac((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx, out_diam);
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if (x_p_dx < sphere_map->wX)
sphere_map->datac((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx, out_diam);
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) )
sphere_map->datac((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx, out_diam);
}
}
}//end FOR DX
}//end FOR DY
}//end FOR DZ
}