#include "StdAfx.h"
#include "SphereGrowth.h"
#include "DistanceTransform.h"
using namespace std;
//////////////////////////////////
//
// GetDiameterMap()
// Creates 'distance_xfrm' cube, that for each Pore voxel, contains the max diameter (in voxels)
// sphere that can fit without overlapping a Fiber. Note, if quick distance xfrm is used there might be
// some overlap.
//
//////////////////////////////////
sg_status SphereGrowth::GetDiameterMap()
{
unsigned short z, y, x;
// Setup progress control and associated variables
int prog_inc_val = ds_cube->wZ/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, (( (ds_cube->wZ)>99 )?(99):(ds_cube->wZ)) );
// Indicate Xfrm and Growth is currently invalid
cur_state &= (~(SG_STATE_XFRM_DONE | SG_STATE_GROWTH_DONE));
// Calculate diameter of max fitting sphere at each Pore voxel
if ( cur_setting & SG_QUICK_DIST_XFRM )
{
// Use distance transform method
write_to_log("SphereGrowth::GetDiameterMap: Creating distance transform using quick method...");
progress_txt->SetWindowText("Obtaining radius estimates using quick distance transform method...");
// Create distance transform run object
DistanceXfrm<sg_xfrm_type> *dist_xfrm_obj = new DistanceXfrm<sg_xfrm_type>(ds_cube, progress_bar);
// Make sure correct face/edge/corner values are set, and run BS dist xfrm
dist_xfrm_obj->changeFECValues(cur_FEC.face, cur_FEC.edge, cur_FEC.corner);
dist_xfrm_obj->runBS(distance_xfrm, DT_INIT_INPUT | \
( (cur_setting & SG_BOUNDARY_IS_SOLID) ? (DT_BOUNDARY_IS_SOLID):(0) ) );
// Normalize and convert radius to diameter
progress_txt->SetWindowText("Normalizing distance transform and converting to diameters...");
// Need to go through entire xfrm cube and do the following at each pore location:
// 1) Divide by the xfrm FACE value to normalize and get true distance value in pixels from boarder
// 2) Multiply by 2 to get diameter
// 3) Subtract one, since we just doubled the center pixel as well when we shouldn't have
FOR_wZwYwX_LOOP(ds_cube, z, y, x)
{
if ( !ds_cube->get_spot(z, y, x) )
{
distance_xfrm->datac(z, y, x,
(sg_xfrm_type)(((unsigned int)2*distance_xfrm->datac(z, y, x))/cur_FEC.face) );
if (distance_xfrm->datac(z, y, x) > 0)
{
distance_xfrm->datac(z, y, x,
distance_xfrm->datac(z, y, x)-1 );
}
// If we are only doing Odd diameters, make sure we are still odd, if not adjust by one to make it so
if ( ( cur_setting & SG_ODD_DIAMETER_ONLY ) &&
((distance_xfrm->datac(z, y, x)%2) == 0) &&
( distance_xfrm->datac(z, y, x) > 0) )
{
distance_xfrm->datac(z, y, x,
distance_xfrm->datac(z, y, x)-1 );
}
}
}
// We are done, no longer need distance transform template object
delete dist_xfrm_obj;
}
else
{
// Actually try to fit a max sphere at each Pore spot
write_to_log("SphereGrowth::GetDiameterMap - Creating distance transform using sphere fitting method...");
progress_txt->SetWindowText("Obtaining diameter estimates using sphere fitting method...");
// Clear xfrm cube
distance_xfrm->clear();
for (z = 0; z < ds_cube->wZ; z++)
{
for(y = 0; y < ds_cube->wY; y++)
{
for(x = 0; x < ds_cube->wX; x++)
{
if ( !ds_cube->get_spot(z, y, x) )
{ // Get exact diameter at each pore, uses neighbor info to speed things up
distance_xfrm->datac(z, y, x, GetExactDiameter(z, y, x));
}
}
}
// Update progress bar when needed
if ( ++prog_cur == prog_next_update )
{
prog_next_update += prog_inc_val;
progress_bar->StepIt();
}
}
}
// Log completion of odd diameter stage
write_to_log("SphereGrowth::GetDiameterMap - Finished calculating all sphere diameters.");
progress_txt->SetWindowText("Finding maximum diameter...");
max_diameter = 0;
FOR_wZwYwX_LOOP(ds_cube, z, y, x)
{
if ( (sg_sphr_type)distance_xfrm->datac(z, y, x) > max_diameter)
{
max_diameter = (sg_sphr_type) distance_xfrm->datac(z, y, x);
}
}
// Output max diameter to log file
write_to_log("SphereGrowth::GetDiameterMap - Maximum sphere diameter found: %d", (int)max_diameter);
// Update state as Xfrm complete
cur_state |= SG_STATE_XFRM_DONE;
return SG_NO_ERROR;
}
////////////////
//function for finding a radius at a given spot... uses neighbors radius's to speed things up
sg_xfrm_type SphereGrowth::GetExactDiameter(unsigned short z, unsigned short y, unsigned short x)
{
sg_xfrm_type prev_max_diam = 0;
// Find the maximum neighboring diameter from previous locations
if ( (x > 0) && (prev_max_diam < distance_xfrm->datac(z, y, x-1)) )
{
prev_max_diam = distance_xfrm->datac(z, y, x-1);
}
if ( (y > 0) && (prev_max_diam < distance_xfrm->datac(z, y-1, x)))
{
prev_max_diam = distance_xfrm->datac(z, y-1, x);
}
if ( (z > 0) && (prev_max_diam < distance_xfrm->datac(z-1, y, x)))
{
prev_max_diam = distance_xfrm->datac(z-1, y, x);
}
// First make sure we don't create a sphere bigger than max allowed
if ( prev_max_diam >= inf_radius )
{
write_to_log("We might have a sphere bigger than max allowed radius at [%dx,%dy,%dz]", x, y, z);
return prev_max_diam;
}
// First try to varify the same size as greatest last space
if ( prev_max_diam > 1 )
{
if ( ValidateOuterSphereLayer(z, y, x, prev_max_diam-2, prev_max_diam) == SG_FIBER_HIT )
{
// If we fail the validation, diameter might be 1 smaller than last time
if ( prev_max_diam > 2 )
{
if ( ValidateOuterSphereLayer(z, y, x, prev_max_diam-3, prev_max_diam-1) == SG_FIBER_HIT )
{
// Must be two smaller
return ( prev_max_diam - 2 );
}
return ( prev_max_diam - 1 );
}
else
{
// We just failed verifying a diam of 2, so must be 1 (other wise it would be a fiber voxel)
return ( 1 );
}
}
}
// If previous max diameter was 0, find maximum we can fit (this is needed if borders are mirrored)
if ( prev_max_diam == 0 )
{
while ( ValidateOuterSphereLayer(z, y, x, prev_max_diam, prev_max_diam+1) != SG_FIBER_HIT )
{
++ prev_max_diam;
}
return ( prev_max_diam );
}
// If we made it this far, try one bigger than last diameter.
// -if we fail, diameter is same as last
if ( ValidateOuterSphereLayer(z, y, x, prev_max_diam-1, prev_max_diam+1) == SG_FIBER_HIT )
{
return ( prev_max_diam );
}
// -if we succeed, try two bigger than last
if ( ValidateOuterSphereLayer(z, y, x, prev_max_diam, prev_max_diam+2) == SG_FIBER_HIT )
{
return ( prev_max_diam + 1 );
}
return ( prev_max_diam + 2 );
}
////////////////
//function for checking if there are any intersections with fibers and outside sphere wall
sg_status SphereGrowth::ValidateOuterSphereLayer(unsigned short z, unsigned short y, unsigned short x, sg_xfrm_type inn_diam, sg_xfrm_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;
unsigned short rad = (out_diam)/2;
unsigned short inn_rad = (inn_diam)/2;
unsigned int dz2, dz2_dy2, rad2, inn_rad2;
if ( (out_diam%2) == 0 )
{
rad2 = rad*rad;
}
else
{
rad2 = (unsigned int)(((double)rad+0.5f)*((double)rad+0.5f));
}
inn_rad2 = inn_rad*inn_rad;
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++)
{
//dy^2
dz2_dy2 = dz2 + dy*dy;
y_p_dy = y+dy; y_m_dy = y-dy;
dx_outer_lim = GetSqrt(rad2 - dz2_dy2);
dx_inner_lim = ( dz2_dy2 > inn_rad2 ) ? (-1) : ( 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 < ds_cube->wZ)
{
if (y_p_dy < ds_cube->wY)
{
if ( (x_p_dx < ds_cube->wX) && (ds_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx)) )
return SG_FIBER_HIT;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && ds_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx))
return SG_FIBER_HIT;
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if ( (x_p_dx < ds_cube->wX) && (ds_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx)) )
return SG_FIBER_HIT;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && ds_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx))
return SG_FIBER_HIT;
}
}
//neg z
if ( (z_m_dz != z_p_dz) && (z_m_dz >= 0) )
{
if (y_p_dy < ds_cube->wY)
{
if ( (x_p_dx < ds_cube->wX) && (ds_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx)) )
return SG_FIBER_HIT;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && ds_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx))
return SG_FIBER_HIT;
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if ( (x_p_dx < ds_cube->wX) && (ds_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx)) )
return SG_FIBER_HIT;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && ds_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx))
return SG_FIBER_HIT;
}
}
}
}
}
return SG_NO_ERROR;
}