#include "StdAfx.h"
#include "TH_MedialAxis.h"
void create_inner_sphere_shade(DS_cube* ds_cube, unsigned short z, unsigned short y, unsigned short x, unsigned short inn_diam, unsigned short out_diam, unsigned char shade)
{
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-1)/2;
unsigned short inn_rad = (inn_diam-1)/2;
unsigned int dz2, dz2_dy2, rad2 = rad*rad, rad_prev2 = inn_rad*inn_rad;
//odd diam
if (out_diam%2)
{
if (inn_diam == 0)
ds_cube->data[z][y][x] = shade;
for (dz = 0; dz <= rad; dz++)
{
dz2 = dz*dz;
z_p_dz = z+dz; z_m_dz = z-dz;
for (dy = 0; dy <= sqrt((float)rad2 - dz2); dy++)
{
//dy^2
dz2_dy2 = dz2 + dy*dy;
y_p_dy = y+dy; y_m_dy = y-dy;
dx_outer_lim = (short)sqrt((float)rad2 - dz2_dy2);
if (dz2_dy2 > rad_prev2)
dx_inner_lim = -1;
else
dx_inner_lim = (short)sqrt((float)rad_prev2 - 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->data[z_p_dz][y_p_dy][x_p_dx] = shade;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0))
ds_cube->data[z_p_dz][y_p_dy][x_m_dx] = shade;
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if (x_p_dx < ds_cube->wX)
ds_cube->data[z_p_dz][y_m_dy][x_p_dx] = shade;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0))
ds_cube->data[z_p_dz][y_m_dy][x_m_dx] = shade;
}
}
//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->data[z_m_dz][y_p_dy][x_p_dx] = shade;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0))
ds_cube->data[z_m_dz][y_p_dy][x_m_dx] = shade;
}
if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) )
{//neg y
if ( x_p_dx < ds_cube->wX )
ds_cube->data[z_m_dz][y_m_dy][x_p_dx] = shade;
if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0))
ds_cube->data[z_m_dz][y_m_dy][x_m_dx] = shade;
}
}
}
}
}
}
//even diam
else
{
if (inn_diam == 0)
{
ds_cube->data[z][y][x] = shade;
if ( (x+1) < ds_cube->wX)
ds_cube->data[z][y][x+1] = shade;
if ( (y+1) < ds_cube->wY)
ds_cube->data[z][y+1][x] = shade;
if ( ((x+1) < ds_cube->wX) && ((y+1) < ds_cube->wY) )
ds_cube->data[z][y+1][x+1] = shade;
if ( (z+1) < ds_cube->wZ)
{
ds_cube->data[z+1][y][x] = shade;
if ( (x+1) < ds_cube->wX)
ds_cube->data[z+1][y][x+1] = shade;
if ( (y+1) < ds_cube->wY)
ds_cube->data[z+1][y+1][x] = shade;
if ( ((x+1) < ds_cube->wX) && ((y+1) < ds_cube->wY) )
ds_cube->data[z+1][y+1][x+1] = shade;
}
}
for (dz = 0; dz <= rad; dz++)
{
dz2 = dz*dz;
z_p_dz = z+dz; z_m_dz = z-dz;
for (dy = 0; dy <= sqrt((float)rad2 - dz2); dy++)
{
//dy^2
dz2_dy2 = dz2 + dy*dy;
y_p_dy = y+dy; y_m_dy = y-dy;
dx_outer_lim = (short)sqrt((float)rad2 - dz2_dy2);
if (dz2_dy2 > rad_prev2)
dx_inner_lim = -1;
else
dx_inner_lim = (short)sqrt((float)rad_prev2 - dz2_dy2);
for (dx = dx_outer_lim; dx > dx_inner_lim; dx--)
{
x_p_dx = x+dx; x_m_dx = x-dx;
//pos z
if (z_p_dz < ds_cube->wZ)
{
if (y_p_dy < ds_cube->wY)
{
if (x_p_dx < ds_cube->wX)
ds_cube->data[z_p_dz+1][y_p_dy+1][x_p_dx+1] = shade;
if ( x_m_dx >= 0)
ds_cube->data[z_p_dz+1][y_p_dy+1][x_m_dx] = shade;
}
if ( y_m_dy >= 0 )
{//neg y
if (x_p_dx < ds_cube->wX)
ds_cube->data[z_p_dz+1][y_m_dy][x_p_dx+1] = shade;
if ( x_m_dx >= 0)
ds_cube->data[z_p_dz+1][y_m_dy][x_m_dx] = shade;
}
}
//neg z
if ( z_m_dz >= 0 )
{
if (y_p_dy < ds_cube->wY)
{
if ( x_p_dx < ds_cube->wX)
ds_cube->data[z_m_dz][y_p_dy+1][x_p_dx+1] = shade;
if ( x_m_dx >= 0)
ds_cube->data[z_m_dz][y_p_dy+1][x_m_dx] = shade;
}
if ( y_m_dy >= 0 )
{//neg y
if ( x_p_dx < ds_cube->wX )
ds_cube->data[z_m_dz][y_m_dy][x_p_dx+1] = shade;
if ( x_m_dx >= 0)
ds_cube->data[z_m_dz][y_m_dy][x_m_dx] = shade;
}
}
}
}
}
}
}
MedialAxis::MedialAxis(DS_1b_cube* ds_cube_ptr, void * medial_dialog_ptr)
{
ds_cube = ds_cube_ptr;
MedialDialog* MedialDlg = (MedialDialog*) medial_dialog_ptr;
progress_txt = &MedialDlg->txt_field_progress;
progress_bar = &MedialDlg->pbar_simulation_progress;
//sample creation
//DS_cube * test_cube = new DS_cube(101, 200, 200);
//for (int z = 0; z < test_cube->wZ; z++)
// for (int y = 0; y < test_cube->wY; y++)
// for (int x = 0; x < test_cube->wX; x++)
// test_cube->data[z][y][x] = 255;
//create_inner_sphere_shade(test_cube, 50, 50, 100, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 100, 50, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 100, 150, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 150, 100, 0, 79, 0);
//
//create_inner_sphere_shade(test_cube, 50, 65, 65, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 135, 65, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 135, 135, 0, 79, 0);
//create_inner_sphere_shade(test_cube, 50, 65, 135, 0, 79, 0);
//test_cube->write_to_RAW_STACK("test_sphere_circle_odd8(200_200_101).stack.raw");
}
MedialAxis::~MedialAxis(void)
{
}
DS_cube* MedialAxis::Thinning()
{
unsigned short x, y, z, iteration = 1;
unsigned char finished = 1;
// float mag;
// vect3f spot_dir;//, neigh_dir;
//vect3i delta;
//this will store the skeleton
skeleton = new DS_cube(ds_cube->wZ, ds_cube->wY, ds_cube->wX);
icube = new DS_cube(ds_cube->wZ, ds_cube->wY, ds_cube->wX);
//initialize work cubes
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))
{
icube->data[z][y][x] = 1;
skeleton->data[z][y][x] = 50;
}
else
{
icube->data[z][y][x] = 0;
skeleton->data[z][y][x] = 0;
}
}
}
}
//find medial surface
while(1)
{
{//log it
char txt[512];
sprintf(txt, "Now on iteration %d...", iteration);
progress_txt->SetWindowText(txt);
}
//for each spot
for (z = 0; z < ds_cube->wZ; z++)
{
for (y = 0; y < ds_cube->wY; y++)
{
for (x = 0; x < ds_cube->wX; x++)
{
//check if it is a pore
if (icube->data[z][y][x] == 0)
{
finished = 0;
//find neighbors
Mark_if_MS26(z, y, x, iteration);
}
}
}
}
++iteration;
if (finished)
break;
finished = 1;
}
////clear icube spots of MS for next step
//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 (skeleton->data[z][y][x] > 50)
// {
// icube->data[z][y][x] = 0;
// skeleton->data[z][y][x] = 0;
// }
// }
////find medial axis from surface
//while(1)
//{
// unsigned short num_found;
// {//log it
// char txt[512];
// sprintf(txt, "Now on iteration %d...", iteration);
// progress_txt->SetWindowText(txt);
// }
// //for each spot
// for (z = 0; z < ds_cube->wZ; z++)
// {
// for (y = 0; y < ds_cube->wY; y++)
// {
// for (x = 0; x < ds_cube->wX; x++)
// {
// //check if it is a pore
// if (icube->data[z][y][x] == 0)
// {
// finished = 0;
// //find neighbors
// num_found = Get_Pore_Vect(z, y, x, icube, iteration, spot_dir);
// if (num_found <= 6)
// {
// skeleton->data[z][y][x] = iteration+200;
// }
// //else
// //{
// icube->data[z][y][x] = iteration+1;
// //}
// //finished = 0;
// ////find neighbors
// //num_found = Get_Pore_Vect(z, y, x, icube, iteration, spot_dir);
// //if (num_found <= 2) //only 1 or two pores... must be MA already
// //{
// // icube->data[z][y][x] = iteration+1;
// // skeleton->data[z][y][x] = iteration+200;
// //}
// //else if (num_found <= 9)
// //{
// // icube->data[z][y][x] = iteration+1;
// // mag = sqrt(spot_dir[0]*spot_dir[0] + spot_dir[1]*spot_dir[1] + spot_dir[2]*spot_dir[2]);
// // Make_Unit_Vector(spot_dir, delta);
// // if (mag < 1.5)
// // {
// // skeleton->data[z][y][x] = iteration+200;
// // }
// // else
// // {
// // }
// //}
// }
// }
// }
// }
//
// ++iteration;
// if (finished)
// break;
// finished = 1;
//}
//clear background if user requests
if (1)
{
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) == 0) && (skeleton->data[z][y][x] < 150) )
skeleton->data[z][y][x] = 20;
else if (skeleton->data[z][y][x] == 50)
skeleton->data[z][y][x] = 0;
}
}
else
{
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))
skeleton->data[z][y][x] = 0;
}
}
icube->write_to_RAW_STACK("iteration_cube(200_200_101).stack.raw");
//clean up
delete icube;
return skeleton;
}
void MedialAxis::Mark_if_MS8(unsigned short z, unsigned short y, unsigned short x, unsigned short iteration)
{
char n=0, s=0, w=0, e=0, u=0, d=0;
char zpb=1, znb=1, ypb=1, ynb=1, xpb=1, xnb=1;
//bound checks
if ( (z+1) >= ds_cube->wZ )
zpb = 0;
else if ( (z-1) < 0 )
znb = 0;
if ( (y+1) >= ds_cube->wY )
ypb = 0;
else if ( (y-1) < 0 )
ynb = 0;
if ( (x+1) >= ds_cube->wX )
xpb = 0;
else if ( (x-1) < 0)
xnb = 0;
/*
odd skele check
*/
if ( zpb && (icube->data[z+1][y][x] == iteration) )
d = 1;
if ( ypb && (icube->data[z][y+1][x] == iteration) )
n = 1;
if ( xpb && (icube->data[z][y][x+1] == iteration) )
e = 1;
if ( znb && (icube->data[z-1][y][x] == iteration) )
u = 1;
if ( ynb && (icube->data[z][y-1][x] == iteration) )
s = 1;
if ( xnb && (icube->data[z][y][x-1] == iteration) )
w = 1;
//if we found neighbors, burn this spot
if (u || n || e || d || s || w)
icube->data[z][y][x] = (unsigned char)(iteration+1);
else
return;
//see if any opposite direction, mark skeleton if so
if (u && d)
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
else if (n && s)
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
else if (w && e)
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
//see if any triangles are formed
if (zpb && znb && ypb && ynb) //in zy plane
{
if ( (icube->data[z+1][y-1][x] == iteration) && (icube->data[z-1][y+1][x] == iteration) )
{//decreasing diag
if (n && u)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (s && d)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
if ( (icube->data[z+1][y+1][x] == iteration) && (icube->data[z-1][y-1][x] == iteration) )
{//increasing diag
if (u && s)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (d && n)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
}
if (xpb && xnb && ypb && ynb) //in xy plane
{
if ( (icube->data[z][y+1][x-1] == iteration) && (icube->data[z][y-1][x+1] == iteration) )
{//decreasing diag
if (n && e)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (w && s)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
if ( (icube->data[z][y+1][x+1] == iteration) && (icube->data[z][y-1][x-1] == iteration) )
{//increasing diag
if (n && w)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (s && e)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
}
if (xpb && xnb && zpb && znb) //in xz plane
{
if ( (icube->data[z+1][y][x-1] == iteration) && (icube->data[z-1][y][x+1] == iteration) )
{//decreasing diag
if (u && e)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (w && d)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
if ( (icube->data[z+1][y][x+1] == iteration) && (icube->data[z-1][y][x-1] == iteration) )
{//increasing diag
if (u && w)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
if (d && e)
{
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
goto end_of_diag;
}
}
}
//more diag stuff
//if (u && n && ((skeleton->data[z+1][y+1][x] > 50) || (icube->data[z+1][y+1][x] == 0)) )
// skeleton->data[z][y][x] = iteration+200;
//else if (u && s && ((skeleton->data[z+1][y-1][x] > 50) || (icube->data[z+1][y-1][x] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (d && s && ((skeleton->data[z-1][y-1][x] > 50) || (icube->data[z-1][y-1][x] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (d && n && ((skeleton->data[z-1][y+1][x] > 50) || (icube->data[z-1][y+1][x] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (n && e && ((skeleton->data[z][y+1][x+1] > 50) || (icube->data[z][y+1][x+1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (n && w && ((skeleton->data[z][y+1][x-1] > 50) || (icube->data[z][y+1][x-1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (s && w && ((skeleton->data[z][y-1][x-1] > 50) || (icube->data[z][y-1][x-1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (s && e && ((skeleton->data[z][y-1][x+1] > 50) || (icube->data[z][y-1][x+1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (u && e && ((skeleton->data[z+1][y][x+1] > 50) || (icube->data[z+1][y][x+1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (u && w && ((skeleton->data[z+1][y][x-1] > 50) || (icube->data[z+1][y][x-1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (d && w && ((skeleton->data[z-1][y][x-1] > 50) || (icube->data[z-1][y][x-1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
//else if (d && e && ((skeleton->data[z-1][y][x+1] > 50) || (icube->data[z-1][y][x+1] == 0)))
// skeleton->data[z][y][x] = iteration+200;
end_of_diag:
if ( (z < 2) && (z >= 2) )
return;
/*
Even skeleton check
*/
if ( d && (z >= 2) && (icube->data[z-1][y][x] == 0) && (icube->data[z-2][y][x] == iteration) ) //up check
{
icube->data[z-1][y][x] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z-1][y][x] = (unsigned char)(iteration+200);
}
if ( s && (y >= 2) && (icube->data[z][y-1][x] == 0) && (icube->data[z][y-2][x] == iteration) ) //north check
{
icube->data[z][y-1][x] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z][y-1][x] = (unsigned char)(iteration+200);
}
if ( e && (x >= 2) && (icube->data[z][y][x-1] == 0) && (icube->data[z][y][x-2] == iteration) ) //west check
{
icube->data[z][y][x-1] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z][y][x-1] = (unsigned char)(iteration+200);
}
if ( u && (z+2 < ds_cube->wZ) && (icube->data[z+1][y][x] == 0) && (icube->data[z+2][y][x] == iteration) ) //down check
{
icube->data[z+1][y][x] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z+1][y][x] = (unsigned char)(iteration+200);
}
if ( n && (y+2 < ds_cube->wY) && (icube->data[z][y+1][x] == 0) && (icube->data[z][y+2][x] == iteration) ) //south check
{
icube->data[z][y+1][x] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z][y+1][x] = (unsigned char)(iteration+200);
}
if ( w && (x+2 < ds_cube->wX) && (icube->data[z][y][x+1] == 0) && (icube->data[z][y][x+2] == iteration) ) //east check
{
icube->data[z][y][x+1] = (unsigned char)(iteration+1);
skeleton->data[z][y][x] = (unsigned char)(iteration+200);
skeleton->data[z][y][x+1] = (unsigned char)(iteration+200);
}
}
void MedialAxis::Mark_if_MS26(unsigned short z, unsigned short y, unsigned short x, unsigned short iteration)
{
int dx, dy, dz;
char mc[3][3][3];//cube that will represent neightboring voxels
char n=0, s=0, w=0, e=0, u=0, d=0, burn=0;
char zpb=1, znb=1, ypb=1, ynb=1, xpb=1, xnb=1;
//bound checks
if ( (z+1) >= ds_cube->wZ )
zpb = 0;
else if ( (z-1) < 0 )
znb = 0;
if ( (y+1) >= ds_cube->wY )
ypb = 0;
else if ( (y-1) < 0 )
ynb = 0;
if ( (x+1) >= ds_cube->wX )
xpb = 0;
else if ( (x-1) < 0)
xnb = 0;
/*
fill in all 26 neighbor voxels info
*/
for (dz = -1; dz <= 1; dz++)
{
for (dy = -1; dy <= 1; dy++)
{
for (dx = -1; dx <= 1; dx++)
{
if (
((dz+z) >= 0) && ((dz+z) < ds_cube->wZ) && /* check z, y, x bounds */
((dy+y) >= 0) && ((dy+y) < ds_cube->wY) &&
((dx+x) >= 0) && ((dx+x) < ds_cube->wX) &&
(icube->data[dz+z][dy+y][dx+x] == iteration)
)
{
mc[dz+1][dy+1][dx+1] = 1;
burn = 1;
}
else
mc[dz+1][dy+1][dx+1] = 0;
}
}
}
//if we found neighbors, burn this spot
if (burn)
icube->data[z][y][x] = iteration+1;
else //otherwise return cause we are in middle of a pore
return;
//mark variables to represent if there are any fibers on appropriate faces
for (dy = 0; dy <= 2; dy++) //up face
{
for (dx = 0; dx <= 2; dx++)
{
if (mc[0][dy][dx])
{
u = 1;
dx = 3; dy = 3;
}
}
}
for (dy = 0; dy <= 2; dy++) //down face
{
for (dx = 0; dx <= 2; dx++)
{
if (mc[2][dy][dx])
{
d = 1;
dx = 3; dy = 3;
}
}
}
for (dz = 0; dz <= 2; dz++) //north face
{
for (dx = 0; dx <= 2; dx++)
{
if (mc[dz][0][dx])
{
n = 1;
dz = 3; dx = 3;
}
}
}
for (dz = 0; dz <= 2; dz++) //south face
{
for (dx = 0; dx <= 2; dx++)
{
if (mc[dz][2][dx])
{
s = 1;
dz = 3; dx = 3;
}
}
}
for (dz = 0; dz <= 2; dz++) //west face
{
for (dy = 0; dy <= 2; dy++)
{
if (mc[dz][dy][0])
{
w = 1;
dz = 3; dy = 3;
}
}
}
for (dz = 0; dz <= 2; dz++) //east face
{
for (dy = 0; dy <= 2; dy++)
{
if (mc[dz][dy][2])
{
e = 1;
dz = 3; dy = 3;
}
}
}
//see if any opposite direction fibers, mark MS if so
if ( (mc[0][1][1] && d) || (mc[2][1][1] && u) ||
(mc[1][0][1] && s) || (mc[1][2][1] && n) ||
(mc[1][1][0] && e) || (mc[1][1][2] && w) )
skeleton->data[z][y][x] = iteration+200;
//see if any triangles are formed
//in zy plane
//
//
// if (zpb && znb && ypb && ynb) //in zy plane
// {
// if ( (icube->data[z+1][y-1][x] == iteration) && (icube->data[z-1][y+1][x] == iteration) )
// {//decreasing diag
// if (n && u)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (s && d)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// if ( (icube->data[z+1][y+1][x] == iteration) && (icube->data[z-1][y-1][x] == iteration) )
// {//increasing diag
// if (u && s)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (d && n)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// }
// if (xpb && xnb && ypb && ynb) //in xy plane
// {
// if ( (icube->data[z][y+1][x-1] == iteration) && (icube->data[z][y-1][x+1] == iteration) )
// {//decreasing diag
// if (n && e)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (w && s)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// if ( (icube->data[z][y+1][x+1] == iteration) && (icube->data[z][y-1][x-1] == iteration) )
// {//increasing diag
// if (n && w)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (s && e)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// }
// if (xpb && xnb && zpb && znb) //in xz plane
// {
// if ( (icube->data[z+1][y][x-1] == iteration) && (icube->data[z-1][y][x+1] == iteration) )
// {//decreasing diag
// if (u && e)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (w && d)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// if ( (icube->data[z+1][y][x+1] == iteration) && (icube->data[z-1][y][x-1] == iteration) )
// {//increasing diag
// if (u && w)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// if (d && e)
// {
// skeleton->data[z][y][x] = iteration+200;
// goto end_of_diag;
// }
// }
// }
//
//end_of_diag:
if ( (z < 2) && (z >= 2) )
return;
/*
Even skeleton check
*/
//if ( d && (z >= 2) && (icube->data[z-1][y][x] == 0) && (icube->data[z-2][y][x] == iteration) ) //up check
//{
// icube->data[z-1][y][x] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z-1][y][x] = iteration+200;
//}
//if ( s && (y >= 2) && (icube->data[z][y-1][x] == 0) && (icube->data[z][y-2][x] == iteration) ) //north check
//{
// icube->data[z][y-1][x] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z][y-1][x] = iteration+200;
//}
//if ( e && (x >= 2) && (icube->data[z][y][x-1] == 0) && (icube->data[z][y][x-2] == iteration) ) //west check
//{
// icube->data[z][y][x-1] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z][y][x-1] = iteration+200;
//}
//if ( u && (z+2 < ds_cube->wZ) && (icube->data[z+1][y][x] == 0) && (icube->data[z+2][y][x] == iteration) ) //down check
//{
// icube->data[z+1][y][x] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z+1][y][x] = iteration+200;
//}
//if ( n && (y+2 < ds_cube->wY) && (icube->data[z][y+1][x] == 0) && (icube->data[z][y+2][x] == iteration) ) //south check
//{
// icube->data[z][y+1][x] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z][y+1][x] = iteration+200;
//}
//if ( w && (x+2 < ds_cube->wX) && (icube->data[z][y][x+1] == 0) && (icube->data[z][y][x+2] == iteration) ) //east check
//{
// icube->data[z][y][x+1] = iteration+1;
// skeleton->data[z][y][x] = iteration+200;
// skeleton->data[z][y][x+1] = iteration+200;
//}
}
/***********************************
MEDIAL AXIS USING VECTORS
***********************************/
DS_cube* MedialAxis::Thinning_Vect()
{
unsigned short x, y, z, iteration = 1;
unsigned char finished = 1;
float mag;
vect3f spot_dir, neigh_dir;
vect3i delta;
//this will store the skeleton
skeleton = new DS_cube(ds_cube->wZ, ds_cube->wY, ds_cube->wX);
icube = new DS_cube(ds_cube->wZ, ds_cube->wY, ds_cube->wX);
//copy ds lightly over to skeleton so we can see original fibers, and initialize icube
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))
{
skeleton->data[z][y][x] = 50;
icube->data[z][y][x] = 1;
}
else
{
skeleton->data[z][y][x] = 0;
icube->data[z][y][x] = 0;
}
}
}
}
//iterate until we are done
while(1)
{
{
char txt[512];
sprintf(txt, "Now on iteration %d...", iteration);
progress_txt->SetWindowText(txt);
}
//for each spot
for (z = 1; z < ds_cube->wZ-1; z++)
{
for (y = 1; y < ds_cube->wY-1; y++)
{
for (x = 1; x < ds_cube->wX-1; x++)
{
//check if it is a pore
if ( (icube->data[z][y][x] == 0) || (icube->data[z][y][x] > iteration) )
{
//check if this spot is on the interface, and if so, get vector dir from fibers
if(Get_Fiber_Vect(z, y, x, iteration, spot_dir))
{
finished = 0;
//burn spot for next iteration
icube->data[z][y][x] = iteration+1;
//check if it is strong in one direction, and what direction
mag = sqrt(spot_dir[0]*spot_dir[0] + spot_dir[1]*spot_dir[1] + spot_dir[2]*spot_dir[2]);
Make_Unit_Vector(spot_dir, delta);
//if it is weak, or direction hits a pore/MA, or opp dir is fiber, mark this as MA
if ( (mag < 1.5) || (icube->data[z+delta[2]][y+delta[1]][x+delta[0]] == 0)
|| (skeleton->data[z+delta[2]][z+delta[1]][z+delta[0]] > 50) )
// || ( (icube->data[z-delta[2]][y-delta[1]][x-delta[0]] > 0)
// && (icube->data[z-delta[2]][y-delta[1]][x-delta[0]] <= iteration)) )
{
skeleton->data[z][y][x] = 200+iteration;
//icube->data[z][y][x] = 0;
}
else
{//probably not part of MA, but now we check for the 'even' case
//go in opposite direction of fibers and see if its a pore
if ( (icube->data[z-delta[2]][y-delta[1]][x-delta[0]] == 0) ||
(icube->data[z-delta[2]][y-delta[1]][x-delta[0]] > iteration) )
{
//it is a pore, so get dir vector of that spot
if (Get_Fiber_Vect(z-delta[2], y-delta[1], x-delta[0], iteration, neigh_dir))
{ //neighbor has a vector value, make sure they don't cancel each other
spot_dir[0] += neigh_dir[0];
spot_dir[1] += neigh_dir[1];
spot_dir[2] += neigh_dir[2];
mag = sqrt(spot_dir[0]*spot_dir[0] + spot_dir[1]*spot_dir[1] + spot_dir[2]*spot_dir[2]);
if (mag <= 1.5)
{
//if (!skeleton->data[z][y][x])
// icube->data[z][y][x] = 0;
skeleton->data[z][y][x] = 200+iteration;
}
//else
// skeleton->data[z][y][x] = 0;
}
}
}
}
}
}
}
}
if (finished)
break;
finished = 1;
++iteration;
}
icube->write_to_RAW_STACK("iteration_cube(200_200_101).stack.raw");
//clean up
delete icube;
return skeleton;
}
void MedialAxis::Make_Unit_Vector(vect3f vector, vect3i ret_vect)
{
float size = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
if (vector[0] < 0)
ret_vect[0] = (int)(vector[0]/size - .50f);
else
ret_vect[0] = (int)(vector[0]/size + .50f);
if (vector[1] < 0)
ret_vect[1] = (int)(vector[1]/size - .50f);
else
ret_vect[1] = (int)(vector[1]/size + .50f);
if (vector[2] < 0)
ret_vect[2] = (int)(vector[2]/size - .50f);
else
ret_vect[2] = (int)(vector[2]/size + .50f);
}
unsigned char MedialAxis::Get_Fiber_Vect(unsigned short z, unsigned short y, unsigned short x, unsigned short iteration, vect3f ret_vect)
{
static unsigned char mark_cube[3][3][3];
char dx, dy, dz, on_edge_flag = 0;
ret_vect[0] = 0.0;
ret_vect[1] = 0.0;
ret_vect[2] = 0.0;
//fill in mark_cube with surroundings
for (dz = -1; dz <= 1; dz++)
{
if ( (((int)z+dz) >= 0) && (((int)z+dz) < ds_cube->wZ) )
{// make sure we are in z cube bounds
for (dy = -1; dy <= 1; dy++)
{
if ( (((int)y+dy) >= 0) && (((int)y+dy) < ds_cube->wY) )
{//make sure we are in y cube bounds
for (dx = -1; dx <= 1; dx++)
{
//mark spot fiber if in bounds and is by interface
if ( (((int)x+dx) >= 0) && (((int)x+dx) < ds_cube->wX)
&& (icube->data[(int)z+dz][(int)y+dy][(int)x+dx] > 0)
&& (icube->data[(int)z+dz][(int)y+dy][(int)x+dx] <= iteration) )
{//all good
if (dz || dy || dx)
{// not center
on_edge_flag++;
if (dz && dy && dx)
{//corner
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
/*ret_vect[0] += 0.57735*dx;
ret_vect[1] += 0.57735*dy;
ret_vect[2] += 0.57735*dz;*/
}
else if ( (dy && dx) || (dz && dx) || (dz && dy))
{//edge
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
/*ret_vect[0] += 0.707107*dx;
ret_vect[1] += 0.707107*dy;
ret_vect[2] += 0.707107*dz;*/
}
else
{//direct neigh
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
}
}
}
}
}
}
}
}
//not on fiber interface, so don't do anything
if (on_edge_flag == 0)
return 0;
return 1;
}
unsigned char MedialAxis::Get_Pore_Vect(unsigned short z, unsigned short y, unsigned short x, DS_cube * cube, unsigned short iteration, vect3f ret_vect)
{
static unsigned char mark_cube[3][3][3];
char dx, dy, dz, pore_count = 0;
ret_vect[0] = 0.0;
ret_vect[1] = 0.0;
ret_vect[2] = 0.0;
//fill in mark_cube with surroundings
for (dz = -1; dz <= 1; dz++)
{
if ( (((int)z+dz) >= 0) && (((int)z+dz) < ds_cube->wZ) )
{// make sure we are in z cube bounds
for (dy = -1; dy <= 1; dy++)
{
if ( (((int)y+dy) >= 0) && (((int)y+dy) < ds_cube->wY) )
{//make sure we are in y cube bounds
for (dx = -1; dx <= 1; dx++)
{
//mark spot if it is a pore
if ( (((int)x+dx) >= 0) && (((int)x+dx) < ds_cube->wX) &&
((cube->data[(int)z+dz][(int)y+dy][(int)x+dx] == 0) ||
(cube->data[(int)z+dz][(int)y+dy][(int)x+dx] > iteration)))
{//all good
if (dz || dy || dx)
{// not center
pore_count++;
if (dz && dy && dx)
{//corner
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
/*ret_vect[0] += 0.57735*dx;
ret_vect[1] += 0.57735*dy;
ret_vect[2] += 0.57735*dz;*/
}
else if ( (dy && dx) || (dz && dx) || (dz && dy))
{//edge
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
/*ret_vect[0] += 0.707107*dx;
ret_vect[1] += 0.707107*dy;
ret_vect[2] += 0.707107*dz;*/
}
else
{//direct neigh
ret_vect[0] += dx;
ret_vect[1] += dy;
ret_vect[2] += dz;
}
}
}
}
}
}
}
}
//return number of pores found
return pore_count;
}