temp_rtsSignedDistance.h
11.8 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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
#include "rtsImplicit3D.h"
#include <iostream>
using namespace std;
#define DIST_MAX 255
#define DIST_UNSIGNED 0
#define DIST_SIGNED 1
float ComputeSurfaceDistance(point3D<float> p0, point3D<float> p1, float val0, float val1, float s)
{
/*This function computes the distance from p0 to the surface, given two points p0 and p1
on either side of the surface (with values v0 and v1 respectively). surface specifies
the value at the surface.
*/
//compute the normalized position of the surface between p0 and p1
float s_norm_pos = (s - val0) / (val1 - val0);
//compute the actual position of the surface
point3D<float> s_pos = p0 + s_norm_pos * (p1 - p0);
//compute the distance from p0 to the surface
float result = (s_pos - p0).Length();
//cout<<"distance: "<<result<<endl;
return result;
}
void CreateBoundaryConditions(rtsImplicit3D<unsigned char>* function, float threshold,
rtsImplicit3D<float>* &result, rtsImplicit3D<bool>* &mask)
{
/*This function creates an initial signed distance function from a threshold image.
All voxels adjacent to the surface specified by the threshold are initialized with a
distance value. Low values are inside, high values are outside.
*/
//current and neighboring voxel flags (false = inside, true = outside)
bool c, x_p, x_n, y_p, y_n, z_p, z_n;
float d_xp, d_xn, d_yp, d_yn, d_zp, d_zn;
float in_out = 1;
//boundary condition function and the mask
result = new rtsImplicit3D<float>(function->DimX(), function->DimY(), function->DimZ());
//get the parameterization
point3D<float> min_domain= function->getMinDomain();
point3D<float> max_domain = function->getMaxDomain();
result->Parameterize(min_domain.x, max_domain.x, min_domain.y, max_domain.y, min_domain.z, max_domain.z);
(*result) = DIST_MAX;
//create a mask
mask = new rtsImplicit3D<bool>(function->DimX(), function->DimY(), function->DimZ());
(*mask) = false;
cout<<"done making boundary condition function"<<endl;
//for each voxel
vector3D<int> size(function->DimX(), function->DimY(), function->DimZ()); //get the function size
int x, y, z;
for(x=0; x<size.x; x++)
for(y=0; y<size.y; y++)
for(z=0; z<size.z; z++)
{
//reset flags
c=x_p=x_n=y_p=y_n=z_p=z_n=true;
in_out = 1.0;
//look at the current voxel
if((*function)(x, y, z) < threshold)
c=false;
else c=true;
//look at each neighboring voxel
if(x-1 < 0) x_n = c; //X
else if((*function)(x-1, y, z) < threshold) x_n = false;
if(x+1 >= size.x) x_p = c;
else if((*function)(x+1, y, z) < threshold) x_p = false;
if(y-1 < 0) y_n = c; //Y
else if((*function)(x, y-1, z) < threshold) y_n = false;
if(y+1 >= size.y) y_p = c;
else if((*function)(x, y+1, z) < threshold) y_p = false;
if(z-1 < 0) z_n = c; //Z
else if((*function)(x, y, z-1) < threshold) z_n = false;
if(z+1 >= size.z) z_p = c;
else if((*function)(x, y, z+1) < threshold) z_p = false;
//set the distance from the isosurface
if(c == false)
in_out = -1.0;
if(x_n != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z),
(*function).getParameter(x-1, y, z),
(*function)(x, y, z),
(*function)(x-1, y, z),
threshold) * in_out);
if(x_p != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z),
(*function).getParameter(x+1, y, z),
(*function)(x, y, z),
(*function)(x+1, y, z),
threshold) * in_out);
if(y_n != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z),
(*function).getParameter(x, y-1, z),
(*function)(x, y, z),
(*function)(x, y-1, z),
threshold) * in_out);
if(y_p != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z),
(*function).getParameter(x, y+1, z),
(*function)(x, y, z),
(*function)(x, y+1, z),
threshold) * in_out);
if(z_n != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z-1),
(*function).getParameter(x, y, z),
(*function)(x, y, z),
(*function)(x, y, z-1),
threshold) * in_out);
if(z_p != c)
(*result)(x, y, z) = min((*result)(x,y,z),
ComputeSurfaceDistance((*function).getParameter(x, y, z),
(*function).getParameter(x, y, z+1),
(*function)(x, y, z),
(*function)(x, y, z+1),
threshold) * in_out);
//set the mask to 1 if the voxel is on an edge node
if(x_n != c || x_p != c || y_n != c || y_p != c || z_n != c || z_p != c)
(*mask)(x, y, z) = true;
}
//if a line between the two voxels crosses the surface
//find the distance between the voxel center and the surface
cout<<"done computing boundary conditions"<<endl;
}
double ManhattanDistance(rtsImplicit3D<float>* function,
point3D<indextype> point,
vector3D<float> voxelsize,
int type = DIST_UNSIGNED)
{
/*This function updates the manhattan distance from a surface using the manhattan
distance of its neighboring points.
*/
indextype x, y, z;
x=point.x; y=point.y, z=point.z;
int sign = 1;
float result = DIST_MAX;
float near_value; //the value of the neighbor being considered
float possible_value;
if(x!=0)
{
near_value = (*function)(x-1, y, z);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.x);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.x);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
if(x!=function->DimX()-1)
{
near_value = (*function)(x+1, y, z);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.x);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.x);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
if(y!=0)
{
near_value = (*function)(x, y-1, z);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.y);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.y);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
if(y!=function->DimY()-1)
{
near_value = (*function)(x, y+1, z);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.y);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.y);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
if(z!=0)
{
near_value = (*function)(x, y, z-1);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.z);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.z);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
if(z!=function->DimZ()-1)
{
near_value = (*function)(x, y, z+1);
if(type == DIST_UNSIGNED)
result = min(result, near_value + voxelsize.z);
else if(type == DIST_SIGNED)
{
if(near_value<0) sign = -1; else sign = 1; //determine if the value is inside or outside
possible_value = sign*(fabs(near_value) + voxelsize.z);
if(fabs(possible_value) < fabs(result))
result = possible_value;
}
}
return result;
}
void Eikonal_Manhattan(rtsImplicit3D<float>* &function,
rtsImplicit3D<bool>* mask,
int type = DIST_UNSIGNED)
{
/*This function estimates the Eikonal equation based on the manhattan distance. The
function constantly increases values outwards from the boundary conditions.
*/
//compute the distance between two voxels
vector3D<float> voxel_size;
voxel_size.x = fabs(function->getParameter(0, 0, 0).x - function->getParameter(1, 0, 0).x);
voxel_size.y = fabs(function->getParameter(0, 0, 0).y - function->getParameter(0, 1, 0).y);
voxel_size.z = fabs(function->getParameter(0, 0, 0).z - function->getParameter(0, 0, 1).z);
//use fast sweeping to compute the manhattan distance
//0:X 0:Y 0:Z
cout<<"first iteration..."<<endl;
int x,y,z;
for(x=0; x<function->DimX(); x++)
for(y=0; y<function->DimY(); y++)
for(z=0; z<function->DimZ(); z++)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"second iteration..."<<endl;
//0:X 0:Y Z:0
for(x=0; x<function->DimX(); x++)
for(y=0; y<function->DimY(); y++)
for(z=function->DimZ()-1; z>=0; z--)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"third iteration..."<<endl;
//0:X Y:0 0:Z
for(x=0; x<function->DimX(); x++)
for(y=function->DimY()-1; y>=0; y--)
for(z=0; z<function->DimZ(); z++)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"fourth iteration..."<<endl;
//0:X Y:0 Z:0
for(x=0; x<function->DimX(); x++)
for(y=function->DimY()-1; y>=0; y--)
for(z=function->DimZ()-1; z>=0; z--)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"fifth iteration..."<<endl;
//X:0 0:Y 0:Z
for(x=function->DimX()-1; x>=0; x--)
for(y=0; y<function->DimY(); y++)
for(z=0; z<function->DimZ(); z++)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"sixth iteration..."<<endl;
//X:0 0:Y Z:0
for(x=function->DimX()-1; x>=0; x--)
for(y=0; y<function->DimY(); y++)
for(z=function->DimZ()-1; z>=0; z--)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"seventh iteration..."<<endl;
//X:0 Y:0 0:Z
for(x=function->DimX()-1; x>=0; x--)
for(y=function->DimY()-1; y>=0; y--)
for(z=0; z<function->DimZ(); z++)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
cout<<"eighth iteration..."<<endl;
//X:0 Y:0 Z:0
for(x=function->DimX()-1; x>=0; x--)
for(y=function->DimY()-1; y>=0; y--)
for(z=function->DimZ()-1; z>=0; z--)
//if the current point is not a boundary value
if(!(*mask)(x, y, z))
(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
cout<<"done."<<endl;
}