1 /*
2 * vp_octree.c
3 *
4 * Routines to create and destroy MinMaxOctree objects for fast classification.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.23 $
29 */
30
31 #include "vp_global.h"
32
33 /*
34 * OctantOrder: octant traversal order, depending on best_view_axis
35 */
36
37 static int OctantOrder[3][8] = {
38 { 0, 2, 4, 6, 1, 3, 5, 7 }, /* VP_X_AXIS */
39 { 0, 4, 1, 5, 2, 6, 3, 7 }, /* VP_Y_AXIS */
40 { 0, 1, 2, 3, 4, 5, 6, 7 } /* VP_Z_AXIS */
41 };
42
43 static void CreatePyramid ANSI_ARGS((vpContext *vpc,
44 void *mm_pyramid[VP_MAX_OCTREE_LEVELS]));
45 static void DescendPyramid ANSI_ARGS((vpContext *vpc,
46 void *mm_pyramid[VP_MAX_OCTREE_LEVELS], int level,
47 int x, int y, int z, int nodes_per_side, void *parent_node,
48 int *octree_offset));
49 static void Compute1DSummedAreaTable ANSI_ARGS((vpContext *vpc));
50 static void Compute2DSummedAreaTable ANSI_ARGS((vpContext *vpc));
51 static void ClassifyOctree1 ANSI_ARGS((vpContext *vpc));
52 static void ClassifyOctree2 ANSI_ARGS((vpContext *vpc));
53 static void ComputeOctreeMask ANSI_ARGS((vpContext *vpc, int level,
54 int xn, int yn, int zn, int node_size, void *parent_node,
55 unsigned char *array, int max_level));
56
57 /*
58 * vpCreateMinMaxOctree
59 *
60 * Create a MinMaxOctree representation of the volume for fast classification.
61 */
62
63 vpResult
vpCreateMinMaxOctree(vpc,root_node_size,base_node_size)64 vpCreateMinMaxOctree(vpc, root_node_size, base_node_size)
65 vpContext *vpc;
66 int root_node_size; /* ignored for now */
67 int base_node_size; /* controls level of detail of smallest nodes */
68 {
69 int max_dim, retcode, p, f;
70 int field_size;
71 int bytes_per_node;
72 int level_size, levels;
73 void *mm_pyramid[VP_MAX_OCTREE_LEVELS];
74 int octree_offset;
75
76 /* check for errors */
77 if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
78 return(retcode);
79 if (vpc->num_clsfy_params <= 0 ||
80 vpc->num_clsfy_params > vpc->num_voxel_fields)
81 return(VPSetError(vpc, VPERROR_BAD_VOXEL));
82 for (p = 0; p < vpc->num_clsfy_params; p++) {
83 f = vpc->param_field[p];
84 if (f < 0 || f >= vpc->num_voxel_fields)
85 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
86 if (p > 0 && f <= vpc->param_field[p-1])
87 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
88 }
89
90 max_dim = vpc->xlen;
91 if (vpc->ylen > max_dim)
92 max_dim = vpc->ylen;
93 if (vpc->zlen > max_dim)
94 max_dim = vpc->zlen;
95 switch (base_node_size) { /* must be a power of 2 */
96 case 1:
97 case 2:
98 case 4:
99 case 8:
100 case 16:
101 case 32:
102 case 64:
103 case 128:
104 case 256:
105 case 512:
106 break;
107 default:
108 return(VPSetError(vpc, VPERROR_BAD_VALUE));
109 }
110 for (p = 0; p < vpc->num_clsfy_params; p++) {
111 if (vpc->field_size[vpc->param_field[p]] == 4)
112 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
113 }
114
115 /* allocate mm_octree structure */
116 Alloc(vpc, vpc->mm_octree, MinMaxOctree *, sizeof(MinMaxOctree),
117 "MinMaxOctree");
118 bzero(vpc->mm_octree, sizeof(MinMaxOctree));
119 vpc->mm_octree->base_node_size = base_node_size;
120
121 /* compute field sizes */
122 bytes_per_node = 0;
123 for (p = 0; p < vpc->num_clsfy_params; p++) {
124 vpc->mm_octree->node_offset[p] = bytes_per_node;
125 bytes_per_node += 2 * vpc->field_size[vpc->param_field[p]];
126 }
127 vpc->mm_octree->range_bytes_per_node = bytes_per_node;
128 vpc->mm_octree->status_offset = bytes_per_node;
129 bytes_per_node++; /* add byte for status field */
130 bytes_per_node = (bytes_per_node + 1) & ~1; /* align to short boundary */
131 vpc->mm_octree->base_bytes_per_node = bytes_per_node;
132 bytes_per_node = (bytes_per_node + 3) & ~3; /* align to word boundary */
133 vpc->mm_octree->child_offset = bytes_per_node;
134 bytes_per_node += sizeof(unsigned); /* add word for child field */
135 vpc->mm_octree->nonbase_bytes_per_node = bytes_per_node;
136
137 /* compute number of levels */
138 levels = 1;
139 level_size = base_node_size;
140 while (level_size < max_dim) {
141 level_size *= 2;
142 levels++;
143 }
144 if (levels >= VP_MAX_OCTREE_LEVELS) {
145 vpDestroyMinMaxOctree(vpc);
146 return(VPSetError(vpc, VPERROR_LIMIT_EXCEEDED));
147 }
148 vpc->mm_octree->levels = levels;
149 vpc->mm_octree->root_node_size = level_size;
150
151 /* build a min-max pyramid representation of the volume */
152 CreatePyramid(vpc, mm_pyramid);
153
154 /* determine how much space is needed for the octree nodes */
155 octree_offset = vpc->mm_octree->nonbase_bytes_per_node; /* root node */
156 DescendPyramid(vpc, mm_pyramid, 0, 0, 0, 0, 1, NULL, &octree_offset);
157
158 /* create the min-max octree nodes */
159 Alloc(vpc, vpc->mm_octree->root, void *, octree_offset, "mm_octree");
160 vpc->mm_octree->octree_bytes = octree_offset;
161 octree_offset = vpc->mm_octree->nonbase_bytes_per_node;
162 Debug((vpc, VPDEBUG_OCTREE, "Octree:\n"));
163 DescendPyramid(vpc, mm_pyramid, 0, 0, 0, 0, 1, vpc->mm_octree->root,
164 &octree_offset);
165
166 /* clean up and return */
167 Dealloc(vpc, mm_pyramid[0]);
168 return(VP_OK);
169 }
170
171 /*
172 * vpDestroyMinMaxOctree
173 *
174 * Destroy the MinMaxOctree representation of the volume.
175 */
176
177 vpResult
vpDestroyMinMaxOctree(vpc)178 vpDestroyMinMaxOctree(vpc)
179 vpContext *vpc;
180 {
181 if (vpc->mm_octree != NULL) {
182 if (vpc->mm_octree->root != NULL) {
183 Dealloc(vpc, vpc->mm_octree->root);
184 vpc->mm_octree->root = NULL;
185 }
186 Dealloc(vpc, vpc->mm_octree);
187 vpc->mm_octree = NULL;
188 }
189 return(VP_OK);
190 }
191
192 /*
193 * CreatePyramid
194 *
195 * Create a min-max pyramid representation of the volume.
196 */
197
198 static void
CreatePyramid(vpc,mm_pyramid)199 CreatePyramid(vpc, mm_pyramid)
200 vpContext *vpc;
201 void *mm_pyramid[VP_MAX_OCTREE_LEVELS];
202 {
203 int pyr_size; /* size of pyramid in bytes */
204 int level, pyr_levels; /* current, total pyramid levels */
205 int level_offset; /* byte offset to beginning of level */
206 int nodes_per_side; /* nodes per side at current level */
207 int node_size; /* voxels per side in node */
208 char *pyr_node; /* current node of pyramid */
209 char *pyr_src; /* pyramid node being read */
210 char *pyr_dst; /* pyramid node being written */
211 char *voxel; /* voxel being read */
212 int x, y, z; /* coordinates of current pyramid node */
213 int nx, ny, nz; /* coordinates of voxel within node */
214 int xlen, ylen, zlen; /* size of volume */
215 int voxel_xstride; /* volume strides */
216 int voxel_ystride;
217 int voxel_zstride;
218 int param_size[VP_MAX_FIELDS]; /* size of each parameter */
219 int param_offset[VP_MAX_FIELDS];/* voxel offset of each parameter */
220 int node_offset[VP_MAX_FIELDS]; /* offset of parameter in octree node */
221 int max_value[VP_MAX_FIELDS]; /* max. value of each parameter in node */
222 int min_value[VP_MAX_FIELDS]; /* min. value of each parameter in node */
223 int num_params; /* number of params for classifier */
224 int p; /* parameter number */
225 int value; /* parameter value */
226 int pyr_bytes_per_node; /* size of node in bytes */
227 int pyr_offsets[8]; /* offsets from pyr_src to each of its
228 neighbors (0,+X,+Y,+XY,+Z,+XZ,+YZ,+XYZ) */
229 int elem; /* index into pyr_offsets */
230
231 /* allocate space for pyramid */
232 ASSERT(vpc->mm_octree != NULL);
233 ASSERT(vpc->mm_octree->levels > 0);
234 ASSERT(vpc->xlen > 0);
235 ASSERT(vpc->raw_voxels != NULL);
236 ASSERT(vpc->num_clsfy_params > 0);
237 pyr_levels = vpc->mm_octree->levels;
238 pyr_size = vpc->mm_octree->base_bytes_per_node;
239 pyr_bytes_per_node = vpc->mm_octree->range_bytes_per_node;
240 for (level = pyr_levels; level > 0; level--)
241 pyr_size = pyr_size*8 + pyr_bytes_per_node;
242 Alloc(vpc, mm_pyramid[0], void *, pyr_size, "mm_pyramid");
243 level_offset = pyr_bytes_per_node;
244 nodes_per_side = 1;
245 for (level = 1; level < vpc->mm_octree->levels; level++) {
246 mm_pyramid[level] = (char *)mm_pyramid[level-1] + level_offset;
247 level_offset *= 8;
248 nodes_per_side *= 2;
249 }
250
251 /* build the base level of the pyramid */
252 xlen = vpc->xlen;
253 ylen = vpc->ylen;
254 zlen = vpc->zlen;
255 voxel_xstride = vpc->xstride;
256 voxel_ystride = vpc->ystride;
257 voxel_zstride = vpc->zstride;
258 voxel = vpc->raw_voxels;
259 num_params = vpc->num_clsfy_params;
260 for (p = 0; p < num_params; p++) {
261 param_size[p] = vpc->field_size[vpc->param_field[p]];
262 param_offset[p] = vpc->field_offset[vpc->param_field[p]];
263 node_offset[p] = vpc->mm_octree->node_offset[p];
264 }
265 node_size = vpc->mm_octree->base_node_size;
266 pyr_dst = mm_pyramid[pyr_levels-1];
267 Debug((vpc, VPDEBUG_PYRAMID, "Pyramid Level %d:\n", pyr_levels-1));
268 for (z = 0; z < nodes_per_side; z++) {
269 ReportStatus(vpc, (double)z / (double)nodes_per_side);
270
271 for (y = 0; y < nodes_per_side; y++) {
272 for (x = 0; x < nodes_per_side; x++) {
273 /* clear the min/max values for the current node */
274 for (p = 0; p < num_params; p++) {
275 max_value[p] = -1;
276 min_value[p] = 65536;
277 }
278
279 /* loop over voxels in the node */
280 if (z * node_size >= zlen || y * node_size >= ylen ||
281 x * node_size >= xlen) {
282 for (p = 0; p < num_params; p++) {
283 max_value[p] = 0;
284 min_value[p] = 0;
285 }
286 voxel += node_size * voxel_zstride;
287 } else for (nz = 0; nz < node_size; nz++) {
288 if (z * node_size + nz >= zlen) {
289 voxel += (node_size - nz) * voxel_zstride;
290 break;
291 }
292 for (ny = 0; ny < node_size; ny++) {
293 if (y * node_size + ny >= ylen) {
294 voxel += (node_size - ny) * voxel_ystride;
295 break;
296 }
297 for (nx = 0; nx < node_size; nx++) {
298 if (x * node_size + nx >= xlen) {
299 voxel += (node_size - nx) * voxel_xstride;
300 break;
301 }
302
303 /* compare each field against current min/max */
304 for (p = 0; p < num_params; p++) {
305 ASSERT(voxel == (char *)vpc->raw_voxels +
306 (x*node_size+nx)*voxel_xstride +
307 (y*node_size+ny)*voxel_ystride +
308 (z*node_size+nz)*voxel_zstride);
309 value = VoxelField(voxel, param_offset[p],
310 param_size[p]);
311 if (value > max_value[p])
312 max_value[p] = value;
313 if (value < min_value[p])
314 min_value[p] = value;
315 }
316 voxel += voxel_xstride;
317 } /* for nx */
318 voxel += voxel_ystride - node_size*voxel_xstride;
319 } /* for ny */
320 voxel += voxel_zstride - node_size*voxel_ystride;
321 } /* for nz */
322
323 /* store the min/max values for this node */
324 Debug((vpc, VPDEBUG_PYRAMID, " Node %d,%d,%d:\n", x, y, z));
325 for (p = 0; p < num_params; p++) {
326 ASSERT(max_value[p] >= 0 && max_value[p] < 65536);
327 ASSERT(min_value[p] >= 0 && min_value[p] < 65536);
328 Debug((vpc, VPDEBUG_PYRAMID,
329 " Param %d: min %d, max %d\n",
330 p, min_value[p], max_value[p]));
331 if (param_size[p] == 1) {
332 ByteField(pyr_dst, node_offset[p]) = min_value[p];
333 ByteField(pyr_dst, node_offset[p]+1) = max_value[p];
334 } else {
335 ASSERT(param_size[p] == 2);
336 ShortField(pyr_dst, node_offset[p]) = min_value[p];
337 ShortField(pyr_dst, node_offset[p]+2) = max_value[p];
338 }
339 }
340 pyr_dst += pyr_bytes_per_node;
341 voxel += node_size * (voxel_xstride - voxel_zstride);
342 } /* for x */
343 voxel += node_size*(voxel_ystride - nodes_per_side*voxel_xstride);
344 } /* for y */
345 voxel += node_size*(voxel_zstride - nodes_per_side*voxel_ystride);
346 } /* for z */
347 ReportStatus(vpc, 1.0);
348
349 /* build the rest of the pyramid */
350 for (level = pyr_levels-2; level >= 0; level--) {
351 ReportStatus(vpc, 1. - (double)(level+1)/(double)(pyr_levels-1));
352 Debug((vpc, VPDEBUG_PYRAMID, "Pyramid Level %d:\n", level));
353 pyr_dst = mm_pyramid[level];
354 pyr_node = mm_pyramid[level+1];
355 pyr_offsets[0] = 0;
356 pyr_offsets[1] = pyr_bytes_per_node;
357 pyr_offsets[2] = nodes_per_side * pyr_bytes_per_node;
358 pyr_offsets[3] = pyr_offsets[2] + pyr_bytes_per_node;
359 pyr_offsets[4] = pyr_offsets[2] * nodes_per_side;
360 pyr_offsets[5] = pyr_offsets[4] + pyr_bytes_per_node;
361 pyr_offsets[6] = pyr_offsets[4] + pyr_offsets[2];
362 pyr_offsets[7] = pyr_offsets[6] + pyr_bytes_per_node;
363 node_size *= 2;
364 nodes_per_side /= 2;
365 for (z = 0; z < nodes_per_side; z++) {
366 for (y = 0; y < nodes_per_side; y++) {
367 for (x = 0; x < nodes_per_side; x++) {
368 /* clear the min/max values for the current node */
369 for (p = 0; p < num_params; p++) {
370 max_value[p] = -1;
371 min_value[p] = 65536;
372 }
373
374 /* loop over the eight children of this node */
375 for (elem = 0; elem < 8; elem++) {
376 pyr_src = pyr_node + pyr_offsets[elem];
377 /* compare min/max values of children with current
378 min/max values for the node */
379 for (p = 0; p < num_params; p++) {
380 value = VoxelField(pyr_src, node_offset[p],
381 param_size[p]);
382 if (value < min_value[p])
383 min_value[p] = value;
384 value = VoxelField(pyr_src, node_offset[p] +
385 param_size[p], param_size[p]);
386 if (value > max_value[p])
387 max_value[p] = value;
388 }
389 }
390
391 /* store the min/max values for this node */
392 Debug((vpc, VPDEBUG_PYRAMID, " Node %d,%d,%d:\n",x,y,z));
393 for (p = 0; p < num_params; p++) {
394 ASSERT(max_value[p] >= 0 && max_value[p] < 65536);
395 ASSERT(min_value[p] >= 0 && min_value[p] < 65536);
396 Debug((vpc, VPDEBUG_PYRAMID,
397 " Param %d: min %d, max %d\n",
398 p, min_value[p], max_value[p]));
399 if (param_size[p] == 1) {
400 ByteField(pyr_dst, node_offset[p]) = min_value[p];
401 ByteField(pyr_dst, node_offset[p]+1)=max_value[p];
402 } else {
403 ASSERT(param_size[p] == 2);
404 ShortField(pyr_dst, node_offset[p]) = min_value[p];
405 ShortField(pyr_dst, node_offset[p]+2)=max_value[p];
406 }
407 }
408
409 /* go on to the next node */
410 pyr_dst += pyr_bytes_per_node;
411 pyr_node += 2*pyr_bytes_per_node;
412 } /* for x */
413 pyr_node += (2*nodes_per_side)*pyr_bytes_per_node;
414 } /* for y */
415 pyr_node += (2*nodes_per_side)*(2*nodes_per_side)*
416 pyr_bytes_per_node;
417 } /* for z */
418 } /* for level */
419 ReportStatus(vpc, 1.0);
420 }
421
422 /*
423 * DescendPyramid
424 *
425 * Descend the pyramid recursively, either to count how many nodes will
426 * be copied to the octree (if parent_node == NULL) or to actually copy them.
427 */
428
429 static void
DescendPyramid(vpc,mm_pyramid,level,x,y,z,nodes_per_side,parent_node,octree_offset)430 DescendPyramid(vpc, mm_pyramid, level, x, y, z, nodes_per_side,
431 parent_node, octree_offset)
432 vpContext *vpc; /* context */
433 void *mm_pyramid[VP_MAX_OCTREE_LEVELS]; /* min-max pyramid */
434 int level; /* current level */
435 int x, y, z; /* current node coordinates (in coordinate system
436 of the current level) */
437 int nodes_per_side; /* # nodes at current level per side of volume */
438 void *parent_node; /* parent octree node (or NULL) */
439 int *octree_offset; /* bytes from root of octree to next free location */
440 {
441 char *pyr_ptr;
442 char *child_node;
443 int p;
444 MinMaxOctree *mm_octree;
445 int pyr_bytes_per_node;
446 int base_bytes_per_node;
447 int nonbase_bytes_per_node;
448 int child_bytes_per_node;
449 int field_size;
450 int field_offset;
451 int child_offset;
452 int range;
453 int subdivide;
454
455 ASSERT(vpc->mm_octree != NULL);
456 mm_octree = vpc->mm_octree;
457 pyr_bytes_per_node = mm_octree->range_bytes_per_node;
458 base_bytes_per_node = mm_octree->base_bytes_per_node;
459 nonbase_bytes_per_node = mm_octree->nonbase_bytes_per_node;
460 child_offset = mm_octree->child_offset;
461 pyr_ptr = (char *)mm_pyramid[level] + ((z*nodes_per_side + y) *
462 nodes_per_side + x) * pyr_bytes_per_node;
463
464 /* copy min/max data from pyramid node to octree node */
465 if (parent_node != NULL) {
466 Debug((vpc, VPDEBUG_OCTREE,
467 " Node at level %d, coords %d,%d,%d, addr 0x%08x\n",
468 level, x, y, z, parent_node));
469 for (p = 0; p < pyr_bytes_per_node; p++)
470 ByteField(parent_node, p) = ByteField(pyr_ptr, p);
471 }
472
473 /* descend to next level */
474 if (level < mm_octree->levels-1) {
475 /* check if we should subdivide node or not */
476 subdivide = 0;
477 for (p = 0; p < vpc->num_clsfy_params; p++) {
478 field_size = vpc->field_size[vpc->param_field[p]];
479 field_offset = mm_octree->node_offset[p];
480 if (field_size == 1) {
481 range = ByteField(pyr_ptr, field_offset+1) -
482 ByteField(pyr_ptr, field_offset);
483 } else {
484 ASSERT(field_size == 2);
485 range = ShortField(pyr_ptr, field_offset+2) -
486 ShortField(pyr_ptr, field_offset);
487 }
488 if (range > vpc->param_maxrange[p]) {
489 subdivide = 1;
490 break;
491 }
492 }
493
494 if (subdivide) {
495 /* store offset to child */
496 if (parent_node != NULL) {
497 child_node = (char *)mm_octree->root + *octree_offset;
498 IntField(parent_node, child_offset) = *octree_offset;
499 Debug((vpc, VPDEBUG_OCTREE,
500 " Storing children at offset = %d, addr = 0x%08x\n",
501 *octree_offset, child_node));
502 }
503 if (level == mm_octree->levels-2)
504 child_bytes_per_node = base_bytes_per_node;
505 else
506 child_bytes_per_node = nonbase_bytes_per_node;
507 *octree_offset += 8 * child_bytes_per_node;
508 if (parent_node == NULL) {
509 child_node = NULL;
510 child_bytes_per_node = 0;
511 }
512
513 /* visit children */
514 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2, z*2,
515 nodes_per_side*2, child_node, octree_offset);
516 child_node += child_bytes_per_node;
517 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2, z*2,
518 nodes_per_side*2, child_node, octree_offset);
519 child_node += child_bytes_per_node;
520 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2+1, z*2,
521 nodes_per_side*2, child_node, octree_offset);
522 child_node += child_bytes_per_node;
523 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2+1, z*2,
524 nodes_per_side*2, child_node, octree_offset);
525 child_node += child_bytes_per_node;
526 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2, z*2+1,
527 nodes_per_side*2, child_node, octree_offset);
528 child_node += child_bytes_per_node;
529 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2, z*2+1,
530 nodes_per_side*2, child_node, octree_offset);
531 child_node += child_bytes_per_node;
532 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2+1, z*2+1,
533 nodes_per_side*2, child_node, octree_offset);
534 child_node += child_bytes_per_node;
535 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1,y*2+1,z*2+1,
536 nodes_per_side*2, child_node, octree_offset);
537 } else {
538 /* node has no children; store NULL pointer */
539 Debug((vpc, VPDEBUG_OCTREE, " Not subdividing.\n"));
540 if (parent_node != NULL) {
541 IntField(parent_node, child_offset) = 0;
542 }
543 }
544 }
545 }
546
547 /*
548 * VPComputeSummedAreaTable
549 *
550 * Build the summed-area table for fast-classification.
551 */
552
553 void
VPComputeSummedAreaTable(vpc)554 VPComputeSummedAreaTable(vpc)
555 vpContext *vpc;
556 {
557 /* use a special-case version for lower dimensions
558 (faster since C optimizer does a better job) */
559 switch (vpc->num_clsfy_params) {
560 case 1:
561 Compute1DSummedAreaTable(vpc);
562 break;
563 case 2:
564 Compute2DSummedAreaTable(vpc);
565 break;
566 default:
567 /* XXX add code for ND classifiers */
568 VPBug("VPComputeSummedAreaTable can only handle 1D or 2D classifiers");
569 break;
570 }
571 }
572
573 /*
574 * Compute1DSummedAreaTable
575 *
576 * Build a 1D summed area table.
577 */
578
579 static void
Compute1DSummedAreaTable(vpc)580 Compute1DSummedAreaTable(vpc)
581 vpContext *vpc;
582 {
583 int p0max, p0value;
584 unsigned table_size;
585 float opacity, min_opacity, *p0table;
586 unsigned sum;
587 unsigned *entry;
588
589 p0max = vpc->field_max[vpc->param_field[0]];
590 table_size = (p0max+1) * sizeof(unsigned);
591 p0table = vpc->clsfy_table[0];
592 min_opacity = vpc->min_opacity;
593 if (vpc->sum_table == NULL || table_size != vpc->sum_table_size) {
594 if (vpc->sum_table != NULL)
595 Dealloc(vpc, vpc->sum_table);
596 Alloc(vpc, vpc->sum_table, unsigned *, table_size, "sum_table");
597 vpc->sum_table_size = table_size;
598 }
599 entry = vpc->sum_table;
600 for (p0value = 0; p0value <= p0max; p0value++) {
601 opacity = p0table[p0value];
602 if (opacity > min_opacity)
603 sum = 1;
604 else
605 sum = 0;
606 if (p0value > 0)
607 sum += entry[-1];
608 entry[0] = sum;
609 entry++;
610 }
611 }
612
613 /*
614 * Compute2DSummedAreaTable
615 *
616 * Build a 2D summed area table.
617 */
618
619 static void
Compute2DSummedAreaTable(vpc)620 Compute2DSummedAreaTable(vpc)
621 vpContext *vpc;
622 {
623 int p0max, p0value, p1max, p1value;
624 unsigned table_size;
625 float opacity, min_opacity, *p0table, *p1table;
626 unsigned sum;
627 unsigned *entry;
628
629 p0max = vpc->field_max[vpc->param_field[0]];
630 p1max = vpc->field_max[vpc->param_field[1]];
631 table_size = (p0max+1) * (p1max+1) * sizeof(unsigned);
632 p0table = vpc->clsfy_table[0];
633 p1table = vpc->clsfy_table[1];
634 min_opacity = vpc->min_opacity;
635 if (vpc->sum_table == NULL || table_size != vpc->sum_table_size) {
636 if (vpc->sum_table != NULL)
637 Dealloc(vpc, vpc->sum_table);
638 Alloc(vpc, vpc->sum_table, unsigned *, table_size, "sum_table");
639 vpc->sum_table_size = table_size;
640 }
641 entry = vpc->sum_table;
642 for (p0value = 0; p0value <= p0max; p0value++) {
643 for (p1value = 0; p1value <= p1max; p1value++) {
644 opacity = p0table[p0value] * p1table[p1value];
645 if (opacity > min_opacity)
646 sum = 1;
647 else
648 sum = 0;
649 if (p1value > 0) {
650 sum += entry[-1];
651 if (p0value > 0) {
652 sum += entry[-(p1max+1)];
653 sum -= entry[-(p1max+1)-1];
654 }
655 } else if (p0value > 0) {
656 sum += entry[-(p1max+1)];
657 }
658 entry[0] = sum;
659 entry++;
660 }
661 }
662 }
663
664 /*
665 * VPClassifyOctree
666 *
667 * Descend an octree and classify each node as full, empty or partfull.
668 */
669
670 void
VPClassifyOctree(vpc)671 VPClassifyOctree(vpc)
672 vpContext *vpc;
673 {
674 /* use a special-case version for lower dimensions
675 (faster since C optimizer does a better job) */
676 switch (vpc->num_clsfy_params) {
677 case 1:
678 ClassifyOctree1(vpc);
679 break;
680 case 2:
681 ClassifyOctree2(vpc);
682 break;
683 default:
684 /* XXX add code for ND classifiers */
685 VPBug("VPClassifyOctree can only handle 2D classifiers");
686 break;
687 }
688 }
689
690 /*
691 * ClassifyOctree1
692 *
693 * Descend an octree and classify each node as full, empty or partfull.
694 * Specialized for a 1 parameter classification function (1D summed
695 * area table).
696 */
697
698 static void
ClassifyOctree1(vpc)699 ClassifyOctree1(vpc)
700 vpContext *vpc;
701 {
702 char *node_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node addresses */
703 int count_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node child counts;
704 when count drops to zero, pop up a level */
705 int level; /* current octree level */
706 int max_level; /* highest octree level */
707 char *octree_root; /* root node of octree */
708 char *node; /* current octree node */
709 unsigned area; /* area computed from the summed-area table */
710 int status; /* classification status of current node */
711 unsigned *sum_table; /* summed area table */
712 int p0max, p0min; /* parameter 0 extrema */
713 int p0size; /* parameter size */
714 int child_offset; /* offset of child field in node */
715 int status_offset; /* offset of status field in node */
716 int base_bytes_per_node; /* size of base node in bytes */
717 int nonbase_bytes_per_node; /* size of nonbase node in bytes */
718 int child_count; /* children left at current level */
719
720 /* initialize */
721 ASSERT(vpc->sum_table != NULL);
722 ASSERT(vpc->mm_octree != NULL);
723 ASSERT(vpc->mm_octree->root != NULL);
724 ASSERT(vpc->sum_table_size == sizeof(unsigned) *
725 (vpc->field_max[vpc->param_field[0]]+1));
726 sum_table = vpc->sum_table;
727 max_level = vpc->mm_octree->levels - 1;
728 octree_root = vpc->mm_octree->root;
729 p0size = vpc->field_size[vpc->param_field[0]];
730 status_offset = vpc->mm_octree->status_offset;
731 child_offset = vpc->mm_octree->child_offset;
732 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
733 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
734 node = octree_root;
735 level = 0;
736
737 /* do a depth-first, preorder traversal of the octree */
738 Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
739 while (1) {
740 /* find min/max values for both parameters in this node */
741 if (p0size == 1) {
742 p0min = ByteField(node, 0)-1;
743 p0max = ByteField(node, 1);
744 } else {
745 p0min = ShortField(node, 0)-1;
746 p0max = ShortField(node, 2);
747 }
748
749 /* integrate the opacities in the node using the summed area table */
750 area = sum_table[p0max];
751 if (p0min >= 0)
752 area -= sum_table[p0min];
753
754 /* decide if node is full, empty or partfull */
755 if (area == 0) {
756 status = MM_EMPTY;
757 } else if (level != max_level && IntField(node, child_offset) != 0 &&
758 area != (p0max - p0min)) {
759 status = MM_PARTFULL;
760 } else {
761 status = MM_FULL;
762 }
763 ByteField(node, status_offset) = status;
764 Debug((vpc, VPDEBUG_CLSFYOCTREE,
765 " Level %d: node is %s (addr 0x%08x)\n", level,
766 status == MM_EMPTY ? "empty" : (status == MM_FULL ? "full" :
767 "partfull"), node));
768
769 /* move to next node in tree traversal */
770 if (status == MM_PARTFULL) {
771 /* move down to first child in next level */
772 node = octree_root + IntField(node, child_offset);
773 Debug((vpc, VPDEBUG_CLSFYOCTREE,
774 " Descending. Children at offset %d, addr 0x%08x\n",
775 IntField(node, child_offset), node));
776 node_stack[level] = node;
777 count_stack[level] = 7; /* number of remaining children */
778 level++;
779 ASSERT(level <= max_level);
780 } else {
781 do {
782 /* move up to a node with unvisited children */
783 Debug((vpc, VPDEBUG_CLSFYOCTREE, " Ascending.\n"));
784 level--;
785 if (level < 0)
786 break;
787 child_count = count_stack[level]--;
788 ASSERT(child_count >= 0 && child_count <= 7);
789 } while (child_count == 0);
790 if (level < 0)
791 break; /* traversal of octree is done! */
792
793 /* descend to the next child of this node */
794 if (level == max_level-1)
795 node = node_stack[level] + base_bytes_per_node;
796 else
797 node = node_stack[level] + nonbase_bytes_per_node;
798 Debug((vpc, VPDEBUG_CLSFYOCTREE,
799 " Descending to child at 0x%08x.\n", node));
800 node_stack[level] = node;
801 level++;
802 ASSERT(level <= max_level);
803 }
804 } /* while (1) */
805 }
806
807 /*
808 * ClassifyOctree2
809 *
810 * Descend an octree and classify each node as full, empty or partfull.
811 * Specialized for a 2 parameter classification function (2D summed
812 * area table).
813 */
814
815 static void
ClassifyOctree2(vpc)816 ClassifyOctree2(vpc)
817 vpContext *vpc;
818 {
819 char *node_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node addresses */
820 int count_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node child counts;
821 when count drops to zero, pop up a level */
822 int level; /* current octree level */
823 int max_level; /* highest octree level */
824 char *octree_root; /* root node of octree */
825 char *node; /* current octree node */
826 unsigned area; /* area computed from the summed-area table */
827 int status; /* classification status of current node */
828 unsigned *sum_table; /* summed area table */
829 int sum_table_dim1; /* size of last dimension of sum_table */
830 int p0max, p0min; /* parameter 0 extrema */
831 int p1max, p1min; /* parameter 1 extrema */
832 int p0size, p1size; /* parameter sizes */
833 int child_offset; /* offset of child field in node */
834 int status_offset; /* offset of status field in node */
835 int base_bytes_per_node; /* size of base node in bytes */
836 int nonbase_bytes_per_node; /* size of nonbase node in bytes */
837 int child_count; /* children left at current level */
838
839 /* initialize */
840 ASSERT(vpc->sum_table != NULL);
841 ASSERT(vpc->mm_octree != NULL);
842 ASSERT(vpc->mm_octree->root != NULL);
843 ASSERT(vpc->sum_table_size == sizeof(unsigned) *
844 (vpc->field_max[vpc->param_field[0]]+1) *
845 (vpc->field_max[vpc->param_field[1]]+1));
846 sum_table = vpc->sum_table;
847 max_level = vpc->mm_octree->levels - 1;
848 octree_root = vpc->mm_octree->root;
849 p0size = vpc->field_size[vpc->param_field[0]];
850 p1size = vpc->field_size[vpc->param_field[1]];
851 sum_table_dim1 = vpc->field_max[vpc->param_field[1]] + 1;
852 status_offset = vpc->mm_octree->status_offset;
853 child_offset = vpc->mm_octree->child_offset;
854 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
855 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
856 node = octree_root;
857 level = 0;
858
859 /* do a depth-first, preorder traversal of the octree */
860 Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
861 while (1) {
862 /* find min/max values for both parameters in this node */
863 if (p0size == 1) {
864 p0min = ByteField(node, 0)-1;
865 p0max = ByteField(node, 1);
866 } else {
867 p0min = ShortField(node, 0)-1;
868 p0max = ShortField(node, 2);
869 }
870 if (p1size == 1) {
871 p1min = ByteField(node, 2*p0size)-1;
872 p1max = ByteField(node, 2*p0size+1);
873 } else {
874 p1min = ShortField(node, 2*p0size)-1;
875 p1max = ShortField(node, 2*p0size+2);
876 }
877
878 /* integrate the opacities in the node using the summed area table */
879 area = sum_table[p0max * sum_table_dim1 + p1max];
880 if (p0min >= 0) {
881 if (p1min >= 0) {
882 area += sum_table[p0min * sum_table_dim1 + p1min];
883 area -= sum_table[p0max * sum_table_dim1 + p1min];
884 }
885 area -= sum_table[p0min * sum_table_dim1 + p1max];
886 } else {
887 if (p1min >= 0)
888 area -= sum_table[p0max * sum_table_dim1 + p1min];
889 }
890
891 /* decide if node is full, empty or partfull */
892 if (area == 0) {
893 status = MM_EMPTY;
894 } else if (level != max_level && IntField(node, child_offset) != 0 &&
895 area != (p1max - p1min)*(p0max - p0min)) {
896 status = MM_PARTFULL;
897 } else {
898 status = MM_FULL;
899 }
900 ByteField(node, status_offset) = status;
901 Debug((vpc, VPDEBUG_CLSFYOCTREE,
902 " Level %d: node is %s (addr 0x%08x)\n", level,
903 status == MM_EMPTY ? "empty" : (status == MM_FULL ? "full" :
904 "partfull"), node));
905
906 /* move to next node in tree traversal */
907 if (status == MM_PARTFULL) {
908 /* move down to first child in next level */
909 node = octree_root + IntField(node, child_offset);
910 Debug((vpc, VPDEBUG_CLSFYOCTREE,
911 " Descending. Children at offset %d, addr 0x%08x\n",
912 IntField(node, child_offset), node));
913 node_stack[level] = node;
914 count_stack[level] = 7; /* number of remaining children */
915 level++;
916 ASSERT(level <= max_level);
917 } else {
918 do {
919 /* move up to a node with unvisited children */
920 Debug((vpc, VPDEBUG_CLSFYOCTREE, " Ascending.\n"));
921 level--;
922 if (level < 0)
923 break;
924 child_count = count_stack[level]--;
925 ASSERT(child_count >= 0 && child_count <= 7);
926 } while (child_count == 0);
927 if (level < 0)
928 break; /* traversal of octree is done! */
929
930 /* descend to the next child of this node */
931 if (level == max_level-1)
932 node = node_stack[level] + base_bytes_per_node;
933 else
934 node = node_stack[level] + nonbase_bytes_per_node;
935 Debug((vpc, VPDEBUG_CLSFYOCTREE,
936 " Descending to child at 0x%08x.\n", node));
937 node_stack[level] = node;
938 level++;
939 ASSERT(level <= max_level);
940 }
941 } /* while (1) */
942 }
943
944 /*
945 * VPInitOctreeLevelStack
946 *
947 * Initialize an MMOctreeLevel stack.
948 */
949
950 void
VPInitOctreeLevelStack(vpc,level_stack,axis,k)951 VPInitOctreeLevelStack(vpc, level_stack, axis, k)
952 vpContext *vpc;
953 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS];
954 int axis; /* principle viewing axis */
955 int k; /* current slice number */
956 {
957 int max_level, level, last_node_size;
958 int child_octant, child_bytes_per_node;
959 int *octant_order;
960
961 ASSERT(vpc->mm_octree != NULL);
962 max_level = vpc->mm_octree->levels-1;
963 level_stack[max_level].level_size = vpc->mm_octree->base_node_size;
964 level_stack[max_level].child_octant = -1;
965 level_stack[max_level].child_offset1 = -1;
966 level_stack[max_level].child_offset2 = -1;
967 level_stack[max_level].child2 = NULL;
968 last_node_size = vpc->mm_octree->base_node_size;
969 octant_order = OctantOrder[axis];
970 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Octants for next scanline:\n"));
971 for (level = max_level-1; level >= 0; level--) {
972 level_stack[level].level_size = last_node_size * 2;
973 child_octant = ((k / last_node_size) & 1) * MM_K_BIT;
974 last_node_size *= 2;
975 level_stack[level].child_octant = child_octant;
976 if (level == max_level-1)
977 child_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
978 else
979 child_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
980 ASSERT(child_octant >= 0 && child_octant < 7);
981 ASSERT(octant_order[child_octant] >= 0 &&
982 octant_order[child_octant] < 8);
983 level_stack[level].child_offset1 = octant_order[child_octant] *
984 child_bytes_per_node;
985 level_stack[level].child_offset2 = octant_order[child_octant+1] *
986 child_bytes_per_node;
987 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Level %d: %d, then %d\n",
988 level,octant_order[child_octant],octant_order[child_octant+1]));
989 }
990 }
991
992 /*
993 * VPComputeScanRuns
994 *
995 * For a given voxel scanline, produce a sequence of run lengths
996 * which give a conservative estimate of the non-transparent portions
997 * of the scanline. The runs are computed by finding which nodes
998 * of the classified min-max octree are intersected by the scanline.
999 *
1000 * The return value is the number of scanlines for which this run data
1001 * is valid.
1002 */
1003
1004 int
VPComputeScanRuns(vpc,level_stack,run_lengths,axis,j,icount)1005 VPComputeScanRuns(vpc, level_stack, run_lengths, axis, j, icount)
1006 vpContext *vpc;
1007 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS]; /* saved state */
1008 unsigned char *run_lengths; /* storage for run lengths */
1009 int axis; /* principle viewing axis */
1010 int j; /* scanline number */
1011 int icount; /* scanline length */
1012 {
1013 int octree_maxlevel;
1014 int level;
1015 int max_level = vpc->mm_octree->levels-1;
1016 int child_octant, child_bytes_per_node;
1017 int base_bytes_per_node, nonbase_bytes_per_node;
1018 int i;
1019 char *octree_root, *node;
1020 int run_type;
1021 int run_length;
1022 int run_piece;
1023 int status_offset;
1024 int child_offset;
1025 int status;
1026 int *octant_order;
1027
1028 Debug((vpc, VPDEBUG_OCTREERUNS, "Runs for scanline %d:\n", j));
1029 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Traversal for scanline %d:\n", j));
1030 ASSERT(vpc->mm_octree != NULL);
1031 ASSERT(vpc->mm_octree->root != NULL);
1032 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
1033 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
1034 status_offset = vpc->mm_octree->status_offset;
1035 child_offset = vpc->mm_octree->child_offset;
1036 octree_maxlevel = -1;
1037 i = icount;
1038 octree_root = vpc->mm_octree->root;
1039 node = octree_root;
1040 level = 0;
1041 run_type = MM_EMPTY;
1042 run_length = 0;
1043 octant_order = OctantOrder[axis];
1044
1045 /* traverse the octree */
1046 while (1) {
1047 /* descend tree to next node which is not partfull */
1048 while (1) {
1049 status = ByteField(node, status_offset);
1050 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Node at level %d: %s\n",
1051 level, status == MM_PARTFULL ? "partfull" :
1052 (status == MM_FULL ? "full" : "empty")));
1053 ASSERT(status == MM_PARTFULL || status == MM_FULL ||
1054 status == MM_EMPTY);
1055 if (status != MM_PARTFULL)
1056 break;
1057 ASSERT(IntField(node, child_offset) != 0);
1058 Debug((vpc, VPDEBUG_OCTREETRAVERSE,
1059 " Children at base %d, offsets %d, %d; ",
1060 IntField(node, child_offset),
1061 level_stack[level].child_offset1,
1062 level_stack[level].child_offset2));
1063 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "status %d, %d\n",
1064 ByteField(octree_root + IntField(node, child_offset) +
1065 level_stack[level].child_offset1, status_offset),
1066 ByteField(octree_root + IntField(node, child_offset) +
1067 level_stack[level].child_offset2, status_offset)));
1068 node = octree_root + IntField(node, child_offset);
1069 level_stack[level].child2 = node+level_stack[level].child_offset2;
1070 node += level_stack[level].child_offset1;
1071 level++;
1072 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Descending.\n"));
1073 ASSERT(level < vpc->mm_octree->levels);
1074 }
1075 if (level > octree_maxlevel)
1076 octree_maxlevel = level;
1077
1078 /* add current node to the list of runs */
1079 run_piece = MIN(level_stack[level].level_size, i);
1080 i -= run_piece;
1081 if (status == run_type) {
1082 run_length += run_piece;
1083 } else {
1084 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " New run.\n"));
1085 while (run_length > 255) {
1086 Debug((vpc, VPDEBUG_OCTREERUNS, " 255 0"));
1087 *run_lengths++ = 255;
1088 *run_lengths++ = 0;
1089 run_length -= 255;
1090 }
1091 Debug((vpc, VPDEBUG_OCTREERUNS, " %d", run_length));
1092 *run_lengths++ = run_length;
1093 run_type ^= 1;
1094 run_length = run_piece;
1095 }
1096 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Added %d to run.\n",
1097 run_piece));
1098 if (i == 0)
1099 break; /* traversal is done */
1100
1101 /* move back up the tree to the next node with unvisited children */
1102 do {
1103 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Ascending.\n"));
1104 level--;
1105 ASSERT(level >= 0);
1106 } while (level_stack[level].child2 == NULL);
1107
1108 /* descend to next child */
1109 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Next child--descending.\n"));
1110 node = level_stack[level].child2;
1111 level_stack[level].child2 = NULL;
1112 level++;
1113 ASSERT(level < vpc->mm_octree->levels);
1114 } /* while (1) */
1115
1116 /* write out the last run */
1117 while (run_length > 255) {
1118 Debug((vpc, VPDEBUG_OCTREERUNS, " 255 0"));
1119 *run_lengths++ = 255;
1120 *run_lengths++ = 0;
1121 run_length -= 255;
1122 }
1123 Debug((vpc, VPDEBUG_OCTREERUNS, " %d", run_length));
1124 *run_lengths++ = run_length;
1125 if (run_type == MM_EMPTY) {
1126 Debug((vpc, VPDEBUG_OCTREERUNS, " 0"));
1127 *run_lengths++ = 0;
1128 }
1129 Debug((vpc, VPDEBUG_OCTREERUNS, "\n"));
1130
1131 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Covered %d scanlines.\n",
1132 level_stack[octree_maxlevel].level_size));
1133
1134 /* update state for next scanline: adjust child_octant and then
1135 use it to compute child_offset1 and child_offset2 */
1136 j += level_stack[octree_maxlevel].level_size;
1137 max_level = vpc->mm_octree->levels-1;
1138 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Octants for next scanline:\n"));
1139 for (level = max_level-1; level >= 0; level--) {
1140 child_octant = level_stack[level].child_octant;
1141 if (level >= octree_maxlevel)
1142 child_octant &= MM_K_BIT;
1143 else if ((j & (level_stack[level].level_size/2)) == 0)
1144 child_octant &= ~MM_J_BIT;
1145 else
1146 child_octant |= MM_J_BIT;
1147 level_stack[level].child_octant = child_octant;
1148
1149 if (level == max_level-1)
1150 child_bytes_per_node = base_bytes_per_node;
1151 else
1152 child_bytes_per_node = nonbase_bytes_per_node;
1153 level_stack[level].child_offset1 = octant_order[child_octant] *
1154 child_bytes_per_node;
1155 level_stack[level].child_offset2 = octant_order[child_octant+1] *
1156 child_bytes_per_node;
1157 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Level %d: %d, then %d\n",
1158 level,octant_order[child_octant],octant_order[child_octant+1]));
1159 }
1160
1161 /* return the number of scanlines for which the run lengths are valid
1162 (which is the size of the smallest octree node the scanline hit) */
1163 return(level_stack[octree_maxlevel].level_size);
1164 }
1165
1166 /*
1167 * vpOctreeMask
1168 *
1169 * Fill a 3D array with a mask computed from an octree.
1170 * Each array element is set to one of three values depending upon
1171 * the value of the corresponding voxel in the octree:
1172 * 0 voxel is definitely transparent
1173 * 255 voxel may be non-transparent
1174 * 128 voxel may be non-transparent, and more detailed information
1175 * is available at deeper levels of the octree which were not
1176 * visited
1177 */
1178
1179 vpResult
vpOctreeMask(vpc,array,array_size,max_level)1180 vpOctreeMask(vpc, array, array_size, max_level)
1181 vpContext *vpc; /* context */
1182 unsigned char *array; /* array for result */
1183 int array_size; /* size of array in bytes */
1184 int max_level;
1185 {
1186 int c;
1187 unsigned char *aptr;
1188 int retcode;
1189
1190 /* error checks */
1191 if (vpc->mm_octree == NULL)
1192 return(VPSetError(vpc, VPERROR_BAD_SIZE));
1193 if ((retcode = VPCheckClassifier(vpc)) == NULL)
1194 return(retcode);
1195 if (array_size != vpc->xlen*vpc->ylen*vpc->zlen)
1196 return(VPSetError(vpc, VPERROR_BAD_SIZE));
1197
1198 /* classify the octree */
1199 VPComputeSummedAreaTable(vpc);
1200 VPClassifyOctree(vpc);
1201 ComputeOctreeMask(vpc, 0, 0, 0, 0, vpc->mm_octree->root_node_size,
1202 vpc->mm_octree->root, array, max_level);
1203 return(VP_OK);
1204 }
1205
1206 /*
1207 * ComputeOctreeMask
1208 *
1209 * Recursive helper function for vpOctreeMask.
1210 */
1211
1212 static void
ComputeOctreeMask(vpc,level,xn,yn,zn,node_size,parent_node,array,max_level)1213 ComputeOctreeMask(vpc, level, xn, yn, zn, node_size, parent_node, array,
1214 max_level)
1215 vpContext *vpc; /* context */
1216 int level; /* current level */
1217 int xn, yn, zn; /* current node coordinates (in coordinate system
1218 of the current level) */
1219 int node_size; /* voxel per side of node at this level */
1220 void *parent_node; /* parent octree node */
1221 unsigned char *array; /* array for storing result */
1222 int max_level; /* deepest level of the tree to visit */
1223 {
1224 char *child_node, *octree_root;
1225 int child_bytes_per_node;
1226 int child_offset;
1227 int status_offset;
1228 int status, value;
1229 int x, y, z, x0, y0, z0, x1, y1, z1;
1230 int array_ystride, array_zstride;
1231
1232 /* initialize */
1233 status_offset = vpc->mm_octree->status_offset;
1234 child_offset = vpc->mm_octree->child_offset;
1235 if (level == vpc->mm_octree->levels-2)
1236 child_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
1237 else
1238 child_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
1239 octree_root = vpc->mm_octree->root;
1240
1241 /* base case */
1242 status = ByteField(parent_node, status_offset);
1243 if (level == max_level || status != MM_PARTFULL) {
1244 if (status == MM_EMPTY)
1245 value = 0;
1246 else if (status == MM_FULL)
1247 value = 255;
1248 else if (status == MM_PARTFULL)
1249 value = 128;
1250 else
1251 VPBug("bad status value in ComputeOctreeMask, nodeaddr = 0x%08x",
1252 parent_node);
1253 x0 = xn * node_size;
1254 y0 = yn * node_size;
1255 z0 = zn * node_size;
1256 x1 = MIN(x0 + node_size, vpc->xlen) - 1;
1257 y1 = MIN(y0 + node_size, vpc->ylen) - 1;
1258 z1 = MIN(z0 + node_size, vpc->zlen) - 1;
1259 array_ystride = vpc->xlen;
1260 array_zstride = vpc->xlen * vpc->ylen;
1261 for (z = z0; z <= z1; z++) {
1262 for (y = y0; y <= y1; y++) {
1263 for (x = x0; x <= x1; x++) {
1264 array[z*array_zstride + y*array_ystride + x] = value;
1265 }
1266 }
1267 }
1268 return;
1269 }
1270 ASSERT(IntField(parent_node, child_offset) != 0);
1271
1272 /* visit children */
1273 child_node = octree_root + IntField(parent_node, child_offset);
1274 ComputeOctreeMask(vpc, level+1, xn*2, yn*2, zn*2, node_size/2,
1275 child_node, array, max_level);
1276 child_node += child_bytes_per_node;
1277 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2, zn*2, node_size/2,
1278 child_node, array, max_level);
1279 child_node += child_bytes_per_node;
1280 ComputeOctreeMask(vpc, level+1, xn*2, yn*2+1, zn*2, node_size/2,
1281 child_node, array, max_level);
1282 child_node += child_bytes_per_node;
1283 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2+1, zn*2, node_size/2,
1284 child_node, array, max_level);
1285 child_node += child_bytes_per_node;
1286 ComputeOctreeMask(vpc, level+1, xn*2, yn*2, zn*2+1, node_size/2,
1287 child_node, array, max_level);
1288 child_node += child_bytes_per_node;
1289 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2, zn*2+1, node_size/2,
1290 child_node, array, max_level);
1291 child_node += child_bytes_per_node;
1292 ComputeOctreeMask(vpc, level+1, xn*2, yn*2+1, zn*2+1, node_size/2,
1293 child_node, array, max_level);
1294 child_node += child_bytes_per_node;
1295 ComputeOctreeMask(vpc, level+1, xn*2+1,yn*2+1,zn*2+1, node_size/2,
1296 child_node, array, max_level);
1297 }
1298
1299 #ifdef DEBUG
1300 /*
1301 * VPCheckRuns
1302 *
1303 * Check a scanline of run lengths for validity by comparing it to
1304 * the raw volume data. Return value is the number of voxels in
1305 * nonzero runs which are actually zero (due to conservative
1306 * approximations.) If an error is detected then VPBug is called.
1307 */
1308
1309 int
VPCheckRuns(vpc,run_lengths,axis,k,j)1310 VPCheckRuns(vpc, run_lengths, axis, k, j)
1311 vpContext *vpc;
1312 unsigned char *run_lengths;/* run lengths */
1313 int axis; /* principle viewing axis */
1314 int k; /* slice number */
1315 int j; /* scanline number */
1316 {
1317 char *voxel;
1318 int i;
1319 int icount;
1320 int count;
1321 int is_non_zero;
1322 int num_runs;
1323 float opacity;
1324 int badpredictions;
1325 int istride;
1326
1327 switch (axis) {
1328 case VP_X_AXIS:
1329 voxel = (char *)vpc->raw_voxels + k*vpc->xstride + j*vpc->zstride;
1330 istride = vpc->ystride;
1331 icount = vpc->ylen;
1332 break;
1333 case VP_Y_AXIS:
1334 voxel = (char *)vpc->raw_voxels + k*vpc->ystride + j*vpc->xstride;
1335 istride = vpc->zstride;
1336 icount = vpc->zlen;
1337 break;
1338 case VP_Z_AXIS:
1339 voxel = (char *)vpc->raw_voxels + k*vpc->zstride + j*vpc->ystride;
1340 istride = vpc->xstride;
1341 icount = vpc->xlen;
1342 break;
1343 default:
1344 VPBug("bad axis in VPCheckRuns");
1345 }
1346
1347 count = 0;
1348 is_non_zero = 1;
1349 num_runs = 0;
1350 badpredictions = 0;
1351 for (i = 0; i < icount; i++) {
1352 while (count == 0) {
1353 count = *run_lengths++;
1354 is_non_zero = !is_non_zero;
1355 if (++num_runs > icount)
1356 VPBug("runaway scanline detected by VPCheckRuns");
1357 }
1358 opacity = VPClassifyVoxel(vpc, voxel);
1359 if (opacity > vpc->min_opacity &&
1360 fabs(opacity - vpc->min_opacity) > 0.001) {
1361 if (!is_non_zero) {
1362 printf("\n");
1363 printf("VPCheckRuns: error on voxel (i,j,k)=(%d,%d,%d), ",
1364 i, j, k);
1365 printf("viewaxis %d\n", axis);
1366 printf("Actual opacity: %17.15f\n", opacity);
1367 printf("Threshold: %17.15f\n", vpc->min_opacity);
1368 VPDumpView(vpc);
1369 VPDumpClassifier(vpc);
1370 VPBug("nonzero voxel in zero run detected by VPCheckRuns");
1371 }
1372 } else {
1373 if (is_non_zero)
1374 badpredictions++;
1375 }
1376 voxel += istride;
1377 count--;
1378 }
1379 if (count != 0)
1380 VPBug("run that overshoots end of scanline detected by VPCheckRuns");
1381 if (!is_non_zero) {
1382 if (*run_lengths != 0)
1383 VPBug("missing 0 run at end of scanline detected by VPCheckRuns");
1384 }
1385 return(badpredictions);
1386 }
1387
1388 /*
1389 * VPTestMinMaxOctree
1390 *
1391 * Test out the MinMaxOctree routines.
1392 */
1393
1394 void
VPTestMinMaxOctree(vpc)1395 VPTestMinMaxOctree(vpc)
1396 vpContext *vpc;
1397 {
1398 int x, y, z;
1399 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS];
1400 unsigned char run_lengths[VP_MAX_VOLUME_DIM];
1401 int badpredictions;
1402 int scans_left;
1403 float accuracy;
1404
1405 ASSERT(vpc->mm_octree != NULL);
1406 VPComputeSummedAreaTable(vpc);
1407 VPClassifyOctree(vpc);
1408
1409 badpredictions = 0;
1410 printf("Checking +Z axis runs...\n");
1411 for (z = 0; z < vpc->zlen; z++) {
1412 ReportStatus(vpc, (double)z / (double)vpc->zlen);
1413 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", z));
1414 VPInitOctreeLevelStack(vpc, level_stack, VP_Z_AXIS, z);
1415 scans_left = 0;
1416 for (y = 0; y < vpc->ylen; y++) {
1417 if (scans_left == 0) {
1418 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1419 VP_Z_AXIS, y, vpc->xlen);
1420 }
1421 scans_left--;
1422 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Z_AXIS, z, y);
1423 }
1424 }
1425 ReportStatus(vpc, 1.0);
1426 accuracy = 1. - ((double)badpredictions /
1427 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1428 printf("VPTestMinMaxOctree: PASSED.\n");
1429 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1430 accuracy*100., badpredictions);
1431
1432 badpredictions = 0;
1433 printf("Checking +Y axis runs...\n");
1434 for (y = 0; y < vpc->ylen; y++) {
1435 ReportStatus(vpc, (double)y / (double)vpc->ylen);
1436 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", y));
1437 VPInitOctreeLevelStack(vpc, level_stack, VP_Y_AXIS, y);
1438 scans_left = 0;
1439 for (x = 0; x < vpc->xlen; x++) {
1440 if (scans_left == 0) {
1441 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1442 VP_Y_AXIS, x, vpc->zlen);
1443 }
1444 scans_left--;
1445 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Y_AXIS, y, x);
1446 }
1447 }
1448 ReportStatus(vpc, 1.0);
1449 accuracy = 1. - ((double)badpredictions /
1450 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1451 printf("VPTestMinMaxOctree: PASSED.\n");
1452 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1453 accuracy*100., badpredictions);
1454
1455 badpredictions = 0;
1456 printf("Checking +X axis runs...\n");
1457 for (x = 0; x < vpc->xlen; x++) {
1458 ReportStatus(vpc, (double)x / (double)vpc->xlen);
1459 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", x));
1460 VPInitOctreeLevelStack(vpc, level_stack, VP_X_AXIS, x);
1461 scans_left = 0;
1462 for (z = 0; z < vpc->zlen; z++) {
1463 if (scans_left == 0) {
1464 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1465 VP_X_AXIS, z, vpc->ylen);
1466 }
1467 scans_left--;
1468 badpredictions += VPCheckRuns(vpc, run_lengths, VP_X_AXIS, x, z);
1469 }
1470 }
1471 ReportStatus(vpc, 1.0);
1472 accuracy = 1. - ((double)badpredictions /
1473 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1474 printf("VPTestMinMaxOctree: PASSED.\n");
1475 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1476 accuracy*100., badpredictions);
1477
1478 badpredictions = 0;
1479 printf("Checking -Z axis runs...\n");
1480 for (z = vpc->zlen-1; z >= 0; z--) {
1481 ReportStatus(vpc, (double)(vpc->zlen-1-z) / (double)vpc->zlen);
1482 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", z));
1483 VPInitOctreeLevelStack(vpc, level_stack, VP_Z_AXIS, z);
1484 scans_left = 0;
1485 for (y = 0; y < vpc->ylen; y++) {
1486 if (scans_left == 0) {
1487 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1488 VP_Z_AXIS, y, vpc->xlen);
1489 }
1490 scans_left--;
1491 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Z_AXIS, z, y);
1492 }
1493 }
1494 ReportStatus(vpc, 1.0);
1495 accuracy = 1. - ((double)badpredictions /
1496 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1497 printf("VPTestMinMaxOctree: PASSED.\n");
1498 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1499 accuracy*100., badpredictions);
1500
1501 badpredictions = 0;
1502 printf("Checking -Y axis runs...\n");
1503 for (y = vpc->ylen-1; y >= 0; y--) {
1504 ReportStatus(vpc, (double)(vpc->ylen-1-y) / (double)vpc->ylen);
1505 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", y));
1506 VPInitOctreeLevelStack(vpc, level_stack, VP_Y_AXIS, y);
1507 scans_left = 0;
1508 for (x = 0; x < vpc->xlen; x++) {
1509 if (scans_left == 0) {
1510 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1511 VP_Y_AXIS, x, vpc->zlen);
1512 }
1513 scans_left--;
1514 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Y_AXIS, y, x);
1515 }
1516 }
1517 ReportStatus(vpc, 1.0);
1518 accuracy = 1. - ((double)badpredictions /
1519 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1520 printf("VPTestMinMaxOctree: PASSED.\n");
1521 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1522 accuracy*100., badpredictions);
1523
1524 badpredictions = 0;
1525 printf("Checking -X axis runs...\n");
1526 for (x = vpc->xlen-1; x >= 0; x--) {
1527 ReportStatus(vpc, (double)(vpc->xlen-1-x) / (double)vpc->xlen);
1528 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", x));
1529 VPInitOctreeLevelStack(vpc, level_stack, VP_X_AXIS, x);
1530 scans_left = 0;
1531 for (z = 0; z < vpc->zlen; z++) {
1532 if (scans_left == 0) {
1533 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
1534 VP_X_AXIS, z, vpc->ylen);
1535 }
1536 scans_left--;
1537 badpredictions += VPCheckRuns(vpc, run_lengths, VP_X_AXIS, x, z);
1538 }
1539 }
1540 ReportStatus(vpc, 1.0);
1541 accuracy = 1. - ((double)badpredictions /
1542 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
1543 printf("VPTestMinMaxOctree: PASSED.\n");
1544 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
1545 accuracy*100., badpredictions);
1546 }
1547 #else
1548
1549 int
VPCheckRuns(vpc,run_lengths,axis,k,j)1550 VPCheckRuns(vpc, run_lengths, axis, k, j)
1551 vpContext *vpc;
1552 unsigned char *run_lengths;/* run lengths */
1553 int axis; /* principle viewing axis */
1554 int k; /* slice number */
1555 int j; /* scanline number */
1556 {
1557 }
1558
1559 void
VPTestMinMaxOctree(vpc)1560 VPTestMinMaxOctree(vpc)
1561 vpContext *vpc;
1562 {
1563 }
1564
1565 #endif /* DEBUG */
1566