rtsSkeleton.h
3.32 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include <vector>
using namespace std;
#include "rtsPoint3d.h"
#include "rts_itkVolume.h"
typedef vector<point3D<double>> filament;
#define BACKGROUND_NODE 0
#define SKELETON_NODE 255
#define ENDPOINT_NODE 128
#define VISITED_NODE 64
class rtsSkeleton
{
rts_itkVolume<unsigned char> implicit_skeleton;
vector<filament> explicit_filaments;
vector<point3D<unsigned int>> explicit_endpoints;
void ComputeEndPoints();
public:
void ConstructExplicitSkeleton();
void SmoothExplicitSkeleton(int iterations);
void SaveExplicitSkeleton(const char* filename);
}
void rtsSkeleton::ComputeEndPoints()
{
/*This function finds all of the filament end points in the implicit skeleton.
The end points are stored in the explicit_endpoints vector.
An end point is defined as a node with (n == 1) or (n > 2) connected components.
*/
unsigned int max_x = implicit_skeleton.DimX();
unsigned int max_y = implicit_skeleton.DimY();
unsigned int max_z = implicit_skeleton.DimZ();
unsigned int conn;
indextype x, y, z;
for(x=0; x<max_x; x++) //go through each voxel in the data set
for(y=0; y<max_y; y++)
for(z=0; z<max_z; z++)
{
if(implicit_skeleton.GetPixel(x, y, z) == SKELETON_NODE) //if the voxel is part of the skeleton
{
conn = Neighbors26(&implicit_skeleton, x, y, z, 1); //test the number of connected components
if(conn == 1 || conn > 2) //if the number of connected components meets the
//criteria, store the node in the explicit_endpoints vector.
{
explicit_endpoints.push_back(point3D<unsigned int>(x, y, z));
}
}
}
//tag the endpoints in the implicit function
vector<point3D<unsigned int>>::iterator i;
for(i=explicit_endpoints.begin(); i!=explicit_endpoints.end(); i++)
implicit_skeleton.xyz(i->x, i->y, i->z) = ENDPOINT_NODE;
//cout<<"Number of end points: "<<explicit_endpoints.size()<<endl;
}
void rtsSkeleton::ConstructExplicitSkeleton()
{
/*The implicit skeleton is constructed as follows:
1) Label all branch nodes.
2) For each branch, find the number of fibers connected to that branch.
a) Create each fiber as a seperate structure and store them in a queue.
b) Iteratively track each fiber.
*/
ComputeEndPoints(); //compute the endpoints for the explicit skeleton and
//label them in the implicit skeleton
//iterate through each end point
vector<point3D<unsigned int>>::iterator i;
rtsFunction3D<unsigned char> local;
unsigned int x, y, z;
for(i = explicit_endpoints.begin(); i != explicit_endpoints.end(); i++)
{
//cout<<"Tracing Filament..."<<endl;
//i->print();
local = getLocalRegion(i->x, i->y, i->z);
local(1, 1, 1) = 0;
//local.toConsole();
for(x=0; x<3; x++)
for(y=0; y<3; y++)
for(z=0; z<3; z++)
if(local.xyz(x, y, z) > VISITED_NODE)
{
//cout<<(int)local.xyz(x, y, z)<<endl;
//cout<<(int)implicit_skeleton.xyz(i->x - 1 + x,
// i->y - 1 + y,
// i->z - 1 + z)<<endl;
TraceFilament((*i), point3D<unsigned int>(i->x - 1 + x,
i->y - 1 + y,
i->z - 1 + z));
}
//since we have tracked all filaments attached to an endpoint
//set the endpoint as visited
implicit_skeleton.xyz(i->x, i->y, i->z) = VISITED_NODE;
}
//cout<<"Tracing complete"<<endl;
}