8
8
import org .joml .Vector3d ;
9
9
10
10
/**
11
- * Find the pixels visible from skeleton points for further use in fitting ellipsoids.
11
+ * Find the pixels visible from skeleton points for further use in fitting
12
+ * ellipsoids.
12
13
*
13
- * The only pixels relevant for fitting an ellipsoid at a given skeleton point are the
14
- * points visible from that skeleton point.
14
+ * The only pixels relevant for fitting an ellipsoid at a given skeleton point
15
+ * are the points visible from that skeleton point.
15
16
*
16
17
* @author Michael Doube
17
18
*
@@ -22,90 +23,87 @@ public class RayCaster {
22
23
public static final int FORE = -1 ;
23
24
/** Background value */
24
25
public static final int BACK = 0 ;
25
-
26
+
26
27
/** Nyquist sampling of the edge length */
27
- private static final double STEP_SIZE = 1 / 2.3 ;
28
-
28
+ private static final double STEP_SIZE = 1 / 2.3 ;
29
+
29
30
/**
30
- * Iterate over all skeleton points, checking whether each foreground point is visible.
31
- * will be very slow brute force.
31
+ * Iterate over all skeleton points, checking whether each foreground point is
32
+ * visible. will be very slow brute force.
32
33
*
33
34
* Multithreaded over z slices (usual pattern)
34
35
*
35
36
* @param pixels
36
37
* @return
37
38
*/
38
39
public static ArrayList <ArrayList <int []>> getVisibleClouds (final List <Vector3d > skeletonPoints ,
39
- final byte [][] pixels , final int w , final int h , final int d ){
40
+ final byte [][] pixels , final int w , final int h , final int d ) {
40
41
final Thread [] threads = Multithreader .newThreads ();
41
42
final int nThreads = threads .length ;
42
43
@ SuppressWarnings ("unchecked" )
43
44
final ArrayList <ArrayList <int []>>[] visibleCloudPointsArray = (ArrayList <ArrayList <int []>>[]) new Object [nThreads ];
44
45
45
- //probably same solution as usual, give each thread a copy of the below structure, then merge them at the end.
46
46
final int nSkeletonPoints = skeletonPoints .size ();
47
-
48
- //give each thread its own list of lists, will be merged later
47
+
48
+ // give each thread its own list of lists, will be merged later
49
49
for (int i = 0 ; i < nThreads ; i ++) {
50
- //List<ArrayList<int[]>> visibleCloudPoints = new ArrayList<ArrayList<int[]>>(nSkeletonPoints);
51
50
visibleCloudPointsArray [i ] = new ArrayList <ArrayList <int []>>(nSkeletonPoints );
52
51
}
53
-
54
-
55
- //multithread this as a stream or as a Multithreader pattern.
56
- //have to handle potential write collisions on visibleCloudPoints
57
- //if multiple threads wanted to write a new visible point at the same time.
52
+
58
53
AtomicInteger ai = new AtomicInteger (0 );
59
54
for (int thread = 0 ; thread < nThreads ; thread ++) {
60
55
final ArrayList <ArrayList <int []>> threadVisibleCloudPoints = visibleCloudPointsArray [thread ];
61
56
threads [thread ] = new Thread (() -> {
62
-
63
- // for (int z = 0; z < d; z++) {
64
- for (int z = ai .getAndIncrement (); z < d ; z = ai .getAndIncrement ()) {
65
- final byte [] slice = pixels [z ];
66
- for (int y = 0 ; y < h ; y ++) {
67
- final int offset = y * w ;
68
- for (int x = 0 ; x < w ; x ++) {
69
- final int value = slice [offset + x ];
70
- if (value == FORE ) {
71
- //continue if not a surface pixel because all 'middle' pixels are occluded.
72
- if (!isSurface (pixels , x , y , z , w , h , d ))
73
- continue ;
74
- //look for fastest List iterator
75
- for (int i = 0 ; i < nSkeletonPoints ; i ++) {
76
- Vector3d v = skeletonPoints .get (i );
77
- final int qx = (int ) v .x ;
78
- final int qy = (int ) v .y ;
79
- final int qz = (int ) v .z ;
80
- final boolean isAVisiblePoint = isVisible (x , y , z , qx , qy , qz , w , h , d , STEP_SIZE , pixels );
81
- if (isAVisiblePoint ) {
82
- //add this pixel to the list of points visible from this skeleton point
83
- threadVisibleCloudPoints .get (i ).add (new int [] {x , y , z });
57
+
58
+ // for (int z = 0; z < d; z++) {
59
+ for (int z = ai .getAndIncrement (); z < d ; z = ai .getAndIncrement ()) {
60
+ final byte [] slice = pixels [z ];
61
+ for (int y = 0 ; y < h ; y ++) {
62
+ final int offset = y * w ;
63
+ for (int x = 0 ; x < w ; x ++) {
64
+ final int value = slice [offset + x ];
65
+ if (value == FORE ) {
66
+ // continue if not a surface pixel because all 'middle' pixels are occluded.
67
+ if (!isSurface (pixels , x , y , z , w , h , d ))
68
+ continue ;
69
+ // look for fastest List iterator
70
+ for (int i = 0 ; i < nSkeletonPoints ; i ++) {
71
+ Vector3d v = skeletonPoints .get (i );
72
+ final int qx = (int ) v .x ;
73
+ final int qy = (int ) v .y ;
74
+ final int qz = (int ) v .z ;
75
+ final boolean isAVisiblePoint = isVisible (x , y , z , qx , qy , qz , w , h , d , STEP_SIZE ,
76
+ pixels );
77
+ if (isAVisiblePoint ) {
78
+ // add this pixel to the list of points visible from this skeleton point
79
+ threadVisibleCloudPoints .get (i ).add (new int [] { x , y , z });
80
+ }
81
+ }
84
82
}
85
83
}
86
84
}
87
85
}
88
- }
86
+ });
89
87
}
90
- });
91
- }
92
88
Multithreader .startAndJoin (threads );
93
-
94
- //merge all the lists into one
95
- //could be done in a parallel stream I guess
89
+
90
+ // merge all the lists into one
91
+ // could be done in a parallel stream I guess
96
92
ArrayList <ArrayList <int []>> visibleCloudPoints = new ArrayList <ArrayList <int []>>(nSkeletonPoints );
97
93
for (int i = 0 ; i < nThreads ; i ++) {
98
94
ArrayList <ArrayList <int []>> threadVisibleCloudPoints = visibleCloudPointsArray [i ];
99
95
for (int s = 0 ; s < nSkeletonPoints ; s ++) {
100
96
visibleCloudPoints .get (s ).addAll (threadVisibleCloudPoints .get (s ));
101
97
}
102
98
}
103
-
99
+
104
100
return visibleCloudPoints ;
105
101
}
106
-
102
+
107
103
/**
108
- * Check whether this pixel (x, y, z) has any background neighbours (26 neighbourhood)
104
+ * Check whether this pixel (x, y, z) has any background neighbours (26
105
+ * neighbourhood)
106
+ *
109
107
* @param x
110
108
* @param y
111
109
* @param z
@@ -114,9 +112,9 @@ public static ArrayList<ArrayList<int[]>> getVisibleClouds(final List<Vector3d>
114
112
* @param d
115
113
* @return true if a background neighbour pixel is found, false otherwise
116
114
*/
117
- private static boolean isSurface (final byte [][] image , final int x , final int y , final int z ,
118
- final int w , final int h , final int d ) {
119
-
115
+ private static boolean isSurface (final byte [][] image , final int x , final int y , final int z , final int w ,
116
+ final int h , final int d ) {
117
+
120
118
if (isBackgroundNeighbour (image , x - 1 , y - 1 , z - 1 , w , h , d ))
121
119
return true ;
122
120
if (isBackgroundNeighbour (image , x , y - 1 , z - 1 , w , h , d ))
@@ -169,11 +167,11 @@ private static boolean isSurface(final byte[][] image, final int x, final int y,
169
167
return true ;
170
168
if (isBackgroundNeighbour (image , x + 1 , y + 1 , z + 1 , w , h , d ))
171
169
return true ;
172
-
173
- //no background neighbours were found, x, y, z is not a surface pixel
170
+
171
+ // no background neighbours were found, x, y, z is not a surface pixel
174
172
return false ;
175
173
}
176
-
174
+
177
175
/**
178
176
*
179
177
* @param image
@@ -185,18 +183,19 @@ private static boolean isSurface(final byte[][] image, final int x, final int y,
185
183
* @param d
186
184
* @return
187
185
*/
188
- private static boolean isBackgroundNeighbour (final byte [][] image ,
189
- final int x , final int y , final int z , final int w , final int h , final int d ) {
190
- //if pixel is within bounds and has background value it is a background neighbour
186
+ private static boolean isBackgroundNeighbour (final byte [][] image , final int x , final int y , final int z ,
187
+ final int w , final int h , final int d ) {
188
+ // if pixel is within bounds and has background value it is a background
189
+ // neighbour
191
190
if (withinBounds (x , y , z , w , h , d ))
192
191
if (getPixel (image , x , y , z , w ) == BACK )
193
192
return true ;
194
-
193
+
195
194
return false ;
196
195
}
197
-
196
+
198
197
/**
199
- * Get pixel in 3D image - no bounds checking
198
+ * Get pixel in 3D image - no bounds checking
200
199
*
201
200
* @param image 3D image
202
201
* @param x x- coordinate
@@ -206,9 +205,9 @@ private static boolean isBackgroundNeighbour(final byte[][] image,
206
205
* @return corresponding pixel
207
206
*/
208
207
private static byte getPixel (final byte [][] image , final int x , final int y , final int z , final int w ) {
209
- return image [z ][x + y * w ];
208
+ return image [z ][x + y * w ];
210
209
}
211
-
210
+
212
211
/**
213
212
* checks whether a pixel at (m, n, o) is within the image boundaries
214
213
*
@@ -227,73 +226,73 @@ private static boolean withinBounds(final int m, final int n, final int o, final
227
226
}
228
227
229
228
/**
230
- * Check whether a straight line can be drawn between pixels p and q, which is not occluded by any foreground pixel
229
+ * Check whether a straight line can be drawn between pixels p and q, which is
230
+ * not occluded by any foreground pixel
231
231
*
232
232
* Note all units are in pixels, so need to decalibrate back to pixel integers.
233
233
*
234
- * @param px starting point x
234
+ * @param px starting point x
235
235
* @param py
236
236
* @param pz
237
237
* @param qx
238
238
* @param qy
239
239
* @param qz
240
- * @param w image width in pixels
240
+ * @param w image width in pixels
241
241
* @param stepSize increment distance along vector in pixels
242
242
* @param image
243
- * @return true if it is possible to pass in a straight line from p to q without hitting an occluding pixel.
243
+ * @return true if it is possible to pass in a straight line from p to q without
244
+ * hitting an occluding pixel.
244
245
*/
245
- private static boolean isVisible (final int px , final int py , final int pz ,
246
- final int qx , final int qy , final int qz ,
247
- final int w , final int h , final int d ,
248
- final double stepSize , final byte [][] pixels ) {
249
-
250
- //calculate unit vector between the two points
246
+ private static boolean isVisible (final int px , final int py , final int pz , final int qx , final int qy , final int qz ,
247
+ final int w , final int h , final int d , final double stepSize , final byte [][] pixels ) {
248
+
249
+ // calculate unit vector between the two points
251
250
final double dx = qx - px ;
252
251
final double dy = qy - py ;
253
252
final double dz = qz - pz ;
254
-
253
+
255
254
final double distance = Math .sqrt (dx * dx + dy * dy + dz * dz );
256
-
257
- //unit vector
255
+
256
+ // unit vector
258
257
final double ux = dx / distance ;
259
258
final double uy = dy / distance ;
260
259
final double uz = dz / distance ;
261
-
262
- //unit vector * step size = increment
260
+
261
+ // unit vector * step size = increment
263
262
final double stepX = ux * stepSize ;
264
263
final double stepY = uy * stepSize ;
265
264
final double stepZ = uz * stepSize ;
266
-
267
- //starting at p, step along the vector until we arrive at q
268
- //need to bump out of the starting pixel by a unit vector
269
- //otherwise starting pixel may trigger a (wrong) foreground hit
265
+
266
+ // starting at p, step along the vector until we arrive at q
267
+ // need to bump out of the starting pixel by a unit vector
268
+ // otherwise starting pixel may trigger a (wrong) foreground hit
270
269
double cursorX = px + ux ;
271
270
double cursorY = py + uy ;
272
271
double cursorZ = pz + uz ;
273
-
274
- //number of steps required
272
+
273
+ // number of steps required
275
274
final int nSteps = (int ) (distance / stepSize );
276
-
275
+
277
276
for (int n = 0 ; n < nSteps ; n ++) {
278
- //calculate precise position along ray
277
+ // calculate precise position along ray
279
278
cursorX += stepX ;
280
279
cursorY += stepY ;
281
280
cursorZ += stepZ ;
282
-
283
- //round it to integer discrete space
281
+
282
+ // round it to integer discrete space
284
283
final int x = (int ) Math .round (cursorX );
285
284
final int y = (int ) Math .round (cursorY );
286
285
final int z = (int ) Math .round (cursorZ );
287
-
288
- //check if the pixel is foreground
286
+
287
+ // check if the pixel is foreground
289
288
if (pixels [z ][y * w + x ] == FORE )
290
289
return false ;
291
- //TODO make sure logic excludes the possibility
292
- //of accidentally counting p or q as occluding pixels.
293
- //easy to do if out by 1.
290
+ // TODO make sure logic excludes the possibility
291
+ // of accidentally counting p or q as occluding pixels.
292
+ // easy to do if out by 1.
294
293
}
295
- //if we got to the end of the loop without hitting something else
296
- //p and q are visible to each other
294
+ // if we got to the end of the loop without hitting something else
295
+ // p and q are visible to each other
297
296
return true ;
298
297
}
299
298
}
0 commit comments