00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "vp_global.h"
00032
00033
00034
00035
00036
00037 static int OctantOrder[3][8] = {
00038 { 0, 2, 4, 6, 1, 3, 5, 7 },
00039 { 0, 4, 1, 5, 2, 6, 3, 7 },
00040 { 0, 1, 2, 3, 4, 5, 6, 7 }
00041 };
00042
00043 static void CreatePyramid ANSI_ARGS((vpContext *vpc,
00044 void *mm_pyramid[VP_MAX_OCTREE_LEVELS]));
00045 static void DescendPyramid ANSI_ARGS((vpContext *vpc,
00046 void *mm_pyramid[VP_MAX_OCTREE_LEVELS], int level,
00047 int x, int y, int z, int nodes_per_side, void *parent_node,
00048 int *octree_offset));
00049 static void Compute1DSummedAreaTable ANSI_ARGS((vpContext *vpc));
00050 static void Compute2DSummedAreaTable ANSI_ARGS((vpContext *vpc));
00051 static void ClassifyOctree1 ANSI_ARGS((vpContext *vpc));
00052 static void ClassifyOctree2 ANSI_ARGS((vpContext *vpc));
00053 static void ComputeOctreeMask ANSI_ARGS((vpContext *vpc, int level,
00054 int xn, int yn, int zn, int node_size, void *parent_node,
00055 unsigned char *array, int max_level));
00056
00057
00058
00059
00060
00061
00062
00063 vpResult
00064 vpCreateMinMaxOctree(vpc, root_node_size, base_node_size)
00065 vpContext *vpc;
00066 int root_node_size;
00067 int base_node_size;
00068 {
00069 int max_dim, retcode, p, f;
00070 int field_size;
00071 int bytes_per_node;
00072 int level_size, levels;
00073 void *mm_pyramid[VP_MAX_OCTREE_LEVELS];
00074 int octree_offset;
00075
00076
00077 if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
00078 return(retcode);
00079 if (vpc->num_clsfy_params <= 0 ||
00080 vpc->num_clsfy_params > vpc->num_voxel_fields)
00081 return(VPSetError(vpc, VPERROR_BAD_VOXEL));
00082 for (p = 0; p < vpc->num_clsfy_params; p++) {
00083 f = vpc->param_field[p];
00084 if (f < 0 || f >= vpc->num_voxel_fields)
00085 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
00086 if (p > 0 && f <= vpc->param_field[p-1])
00087 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
00088 }
00089
00090 max_dim = vpc->xlen;
00091 if (vpc->ylen > max_dim)
00092 max_dim = vpc->ylen;
00093 if (vpc->zlen > max_dim)
00094 max_dim = vpc->zlen;
00095 switch (base_node_size) {
00096 case 1:
00097 case 2:
00098 case 4:
00099 case 8:
00100 case 16:
00101 case 32:
00102 case 64:
00103 case 128:
00104 case 256:
00105 case 512:
00106 break;
00107 default:
00108 return(VPSetError(vpc, VPERROR_BAD_VALUE));
00109 }
00110 for (p = 0; p < vpc->num_clsfy_params; p++) {
00111 if (vpc->field_size[vpc->param_field[p]] == 4)
00112 return(VPSetError(vpc, VPERROR_BAD_CLASSIFIER));
00113 }
00114
00115
00116 Alloc(vpc, vpc->mm_octree, MinMaxOctree *, sizeof(MinMaxOctree),
00117 "MinMaxOctree");
00118 bzero(vpc->mm_octree, sizeof(MinMaxOctree));
00119 vpc->mm_octree->base_node_size = base_node_size;
00120
00121
00122 bytes_per_node = 0;
00123 for (p = 0; p < vpc->num_clsfy_params; p++) {
00124 vpc->mm_octree->node_offset[p] = bytes_per_node;
00125 bytes_per_node += 2 * vpc->field_size[vpc->param_field[p]];
00126 }
00127 vpc->mm_octree->range_bytes_per_node = bytes_per_node;
00128 vpc->mm_octree->status_offset = bytes_per_node;
00129 bytes_per_node++;
00130 bytes_per_node = (bytes_per_node + 1) & ~1;
00131 vpc->mm_octree->base_bytes_per_node = bytes_per_node;
00132 bytes_per_node = (bytes_per_node + 3) & ~3;
00133 vpc->mm_octree->child_offset = bytes_per_node;
00134 bytes_per_node += sizeof(unsigned);
00135 vpc->mm_octree->nonbase_bytes_per_node = bytes_per_node;
00136
00137
00138 levels = 1;
00139 level_size = base_node_size;
00140 while (level_size < max_dim) {
00141 level_size *= 2;
00142 levels++;
00143 }
00144 if (levels >= VP_MAX_OCTREE_LEVELS) {
00145 vpDestroyMinMaxOctree(vpc);
00146 return(VPSetError(vpc, VPERROR_LIMIT_EXCEEDED));
00147 }
00148 vpc->mm_octree->levels = levels;
00149 vpc->mm_octree->root_node_size = level_size;
00150
00151
00152 CreatePyramid(vpc, mm_pyramid);
00153
00154
00155 octree_offset = vpc->mm_octree->nonbase_bytes_per_node;
00156 DescendPyramid(vpc, mm_pyramid, 0, 0, 0, 0, 1, NULL, &octree_offset);
00157
00158
00159 Alloc(vpc, vpc->mm_octree->root, void *, octree_offset, "mm_octree");
00160 vpc->mm_octree->octree_bytes = octree_offset;
00161 octree_offset = vpc->mm_octree->nonbase_bytes_per_node;
00162 Debug((vpc, VPDEBUG_OCTREE, "Octree:\n"));
00163 DescendPyramid(vpc, mm_pyramid, 0, 0, 0, 0, 1, vpc->mm_octree->root,
00164 &octree_offset);
00165
00166
00167 Dealloc(vpc, mm_pyramid[0]);
00168 return(VP_OK);
00169 }
00170
00171
00172
00173
00174
00175
00176
00177 vpResult
00178 vpDestroyMinMaxOctree(vpc)
00179 vpContext *vpc;
00180 {
00181 if (vpc->mm_octree != NULL) {
00182 if (vpc->mm_octree->root != NULL) {
00183 Dealloc(vpc, vpc->mm_octree->root);
00184 vpc->mm_octree->root = NULL;
00185 }
00186 Dealloc(vpc, vpc->mm_octree);
00187 vpc->mm_octree = NULL;
00188 }
00189 return(VP_OK);
00190 }
00191
00192
00193
00194
00195
00196
00197
00198 static void
00199 CreatePyramid(vpc, mm_pyramid)
00200 vpContext *vpc;
00201 void *mm_pyramid[VP_MAX_OCTREE_LEVELS];
00202 {
00203 int pyr_size;
00204 int level, pyr_levels;
00205 int level_offset;
00206 int nodes_per_side;
00207 int node_size;
00208 char *pyr_node;
00209 char *pyr_src;
00210 char *pyr_dst;
00211 char *voxel;
00212 int x, y, z;
00213 int nx, ny, nz;
00214 int xlen, ylen, zlen;
00215 int voxel_xstride;
00216 int voxel_ystride;
00217 int voxel_zstride;
00218 int param_size[VP_MAX_FIELDS];
00219 int param_offset[VP_MAX_FIELDS];
00220 int node_offset[VP_MAX_FIELDS];
00221 int max_value[VP_MAX_FIELDS];
00222 int min_value[VP_MAX_FIELDS];
00223 int num_params;
00224 int p;
00225 int value;
00226 int pyr_bytes_per_node;
00227 int pyr_offsets[8];
00228
00229 int elem;
00230
00231
00232 ASSERT(vpc->mm_octree != NULL);
00233 ASSERT(vpc->mm_octree->levels > 0);
00234 ASSERT(vpc->xlen > 0);
00235 ASSERT(vpc->raw_voxels != NULL);
00236 ASSERT(vpc->num_clsfy_params > 0);
00237 pyr_levels = vpc->mm_octree->levels;
00238 pyr_size = vpc->mm_octree->base_bytes_per_node;
00239 pyr_bytes_per_node = vpc->mm_octree->range_bytes_per_node;
00240 for (level = pyr_levels; level > 0; level--)
00241 pyr_size = pyr_size*8 + pyr_bytes_per_node;
00242 Alloc(vpc, mm_pyramid[0], void *, pyr_size, "mm_pyramid");
00243 level_offset = pyr_bytes_per_node;
00244 nodes_per_side = 1;
00245 for (level = 1; level < vpc->mm_octree->levels; level++) {
00246 mm_pyramid[level] = (char *)mm_pyramid[level-1] + level_offset;
00247 level_offset *= 8;
00248 nodes_per_side *= 2;
00249 }
00250
00251
00252 xlen = vpc->xlen;
00253 ylen = vpc->ylen;
00254 zlen = vpc->zlen;
00255 voxel_xstride = vpc->xstride;
00256 voxel_ystride = vpc->ystride;
00257 voxel_zstride = vpc->zstride;
00258 voxel = vpc->raw_voxels;
00259 num_params = vpc->num_clsfy_params;
00260 for (p = 0; p < num_params; p++) {
00261 param_size[p] = vpc->field_size[vpc->param_field[p]];
00262 param_offset[p] = vpc->field_offset[vpc->param_field[p]];
00263 node_offset[p] = vpc->mm_octree->node_offset[p];
00264 }
00265 node_size = vpc->mm_octree->base_node_size;
00266 pyr_dst = mm_pyramid[pyr_levels-1];
00267 Debug((vpc, VPDEBUG_PYRAMID, "Pyramid Level %d:\n", pyr_levels-1));
00268 for (z = 0; z < nodes_per_side; z++) {
00269 ReportStatus(vpc, (double)z / (double)nodes_per_side);
00270
00271 for (y = 0; y < nodes_per_side; y++) {
00272 for (x = 0; x < nodes_per_side; x++) {
00273
00274 for (p = 0; p < num_params; p++) {
00275 max_value[p] = -1;
00276 min_value[p] = 65536;
00277 }
00278
00279
00280 if (z * node_size >= zlen || y * node_size >= ylen ||
00281 x * node_size >= xlen) {
00282 for (p = 0; p < num_params; p++) {
00283 max_value[p] = 0;
00284 min_value[p] = 0;
00285 }
00286 voxel += node_size * voxel_zstride;
00287 } else for (nz = 0; nz < node_size; nz++) {
00288 if (z * node_size + nz >= zlen) {
00289 voxel += (node_size - nz) * voxel_zstride;
00290 break;
00291 }
00292 for (ny = 0; ny < node_size; ny++) {
00293 if (y * node_size + ny >= ylen) {
00294 voxel += (node_size - ny) * voxel_ystride;
00295 break;
00296 }
00297 for (nx = 0; nx < node_size; nx++) {
00298 if (x * node_size + nx >= xlen) {
00299 voxel += (node_size - nx) * voxel_xstride;
00300 break;
00301 }
00302
00303
00304 for (p = 0; p < num_params; p++) {
00305 ASSERT(voxel == (char *)vpc->raw_voxels +
00306 (x*node_size+nx)*voxel_xstride +
00307 (y*node_size+ny)*voxel_ystride +
00308 (z*node_size+nz)*voxel_zstride);
00309 value = VoxelField(voxel, param_offset[p],
00310 param_size[p]);
00311 if (value > max_value[p])
00312 max_value[p] = value;
00313 if (value < min_value[p])
00314 min_value[p] = value;
00315 }
00316 voxel += voxel_xstride;
00317 }
00318 voxel += voxel_ystride - node_size*voxel_xstride;
00319 }
00320 voxel += voxel_zstride - node_size*voxel_ystride;
00321 }
00322
00323
00324 Debug((vpc, VPDEBUG_PYRAMID, " Node %d,%d,%d:\n", x, y, z));
00325 for (p = 0; p < num_params; p++) {
00326 ASSERT(max_value[p] >= 0 && max_value[p] < 65536);
00327 ASSERT(min_value[p] >= 0 && min_value[p] < 65536);
00328 Debug((vpc, VPDEBUG_PYRAMID,
00329 " Param %d: min %d, max %d\n",
00330 p, min_value[p], max_value[p]));
00331 if (param_size[p] == 1) {
00332 ByteField(pyr_dst, node_offset[p]) = min_value[p];
00333 ByteField(pyr_dst, node_offset[p]+1) = max_value[p];
00334 } else {
00335 ASSERT(param_size[p] == 2);
00336 ShortField(pyr_dst, node_offset[p]) = min_value[p];
00337 ShortField(pyr_dst, node_offset[p]+2) = max_value[p];
00338 }
00339 }
00340 pyr_dst += pyr_bytes_per_node;
00341 voxel += node_size * (voxel_xstride - voxel_zstride);
00342 }
00343 voxel += node_size*(voxel_ystride - nodes_per_side*voxel_xstride);
00344 }
00345 voxel += node_size*(voxel_zstride - nodes_per_side*voxel_ystride);
00346 }
00347 ReportStatus(vpc, 1.0);
00348
00349
00350 for (level = pyr_levels-2; level >= 0; level--) {
00351 ReportStatus(vpc, 1. - (double)(level+1)/(double)(pyr_levels-1));
00352 Debug((vpc, VPDEBUG_PYRAMID, "Pyramid Level %d:\n", level));
00353 pyr_dst = mm_pyramid[level];
00354 pyr_node = mm_pyramid[level+1];
00355 pyr_offsets[0] = 0;
00356 pyr_offsets[1] = pyr_bytes_per_node;
00357 pyr_offsets[2] = nodes_per_side * pyr_bytes_per_node;
00358 pyr_offsets[3] = pyr_offsets[2] + pyr_bytes_per_node;
00359 pyr_offsets[4] = pyr_offsets[2] * nodes_per_side;
00360 pyr_offsets[5] = pyr_offsets[4] + pyr_bytes_per_node;
00361 pyr_offsets[6] = pyr_offsets[4] + pyr_offsets[2];
00362 pyr_offsets[7] = pyr_offsets[6] + pyr_bytes_per_node;
00363 node_size *= 2;
00364 nodes_per_side /= 2;
00365 for (z = 0; z < nodes_per_side; z++) {
00366 for (y = 0; y < nodes_per_side; y++) {
00367 for (x = 0; x < nodes_per_side; x++) {
00368
00369 for (p = 0; p < num_params; p++) {
00370 max_value[p] = -1;
00371 min_value[p] = 65536;
00372 }
00373
00374
00375 for (elem = 0; elem < 8; elem++) {
00376 pyr_src = pyr_node + pyr_offsets[elem];
00377
00378
00379 for (p = 0; p < num_params; p++) {
00380 value = VoxelField(pyr_src, node_offset[p],
00381 param_size[p]);
00382 if (value < min_value[p])
00383 min_value[p] = value;
00384 value = VoxelField(pyr_src, node_offset[p] +
00385 param_size[p], param_size[p]);
00386 if (value > max_value[p])
00387 max_value[p] = value;
00388 }
00389 }
00390
00391
00392 Debug((vpc, VPDEBUG_PYRAMID, " Node %d,%d,%d:\n",x,y,z));
00393 for (p = 0; p < num_params; p++) {
00394 ASSERT(max_value[p] >= 0 && max_value[p] < 65536);
00395 ASSERT(min_value[p] >= 0 && min_value[p] < 65536);
00396 Debug((vpc, VPDEBUG_PYRAMID,
00397 " Param %d: min %d, max %d\n",
00398 p, min_value[p], max_value[p]));
00399 if (param_size[p] == 1) {
00400 ByteField(pyr_dst, node_offset[p]) = min_value[p];
00401 ByteField(pyr_dst, node_offset[p]+1)=max_value[p];
00402 } else {
00403 ASSERT(param_size[p] == 2);
00404 ShortField(pyr_dst, node_offset[p]) = min_value[p];
00405 ShortField(pyr_dst, node_offset[p]+2)=max_value[p];
00406 }
00407 }
00408
00409
00410 pyr_dst += pyr_bytes_per_node;
00411 pyr_node += 2*pyr_bytes_per_node;
00412 }
00413 pyr_node += (2*nodes_per_side)*pyr_bytes_per_node;
00414 }
00415 pyr_node += (2*nodes_per_side)*(2*nodes_per_side)*
00416 pyr_bytes_per_node;
00417 }
00418 }
00419 ReportStatus(vpc, 1.0);
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429 static void
00430 DescendPyramid(vpc, mm_pyramid, level, x, y, z, nodes_per_side,
00431 parent_node, octree_offset)
00432 vpContext *vpc;
00433 void *mm_pyramid[VP_MAX_OCTREE_LEVELS];
00434 int level;
00435 int x, y, z;
00436
00437 int nodes_per_side;
00438 void *parent_node;
00439 int *octree_offset;
00440 {
00441 char *pyr_ptr;
00442 char *child_node;
00443 int p;
00444 MinMaxOctree *mm_octree;
00445 int pyr_bytes_per_node;
00446 int base_bytes_per_node;
00447 int nonbase_bytes_per_node;
00448 int child_bytes_per_node;
00449 int field_size;
00450 int field_offset;
00451 int child_offset;
00452 int range;
00453 int subdivide;
00454
00455 ASSERT(vpc->mm_octree != NULL);
00456 mm_octree = vpc->mm_octree;
00457 pyr_bytes_per_node = mm_octree->range_bytes_per_node;
00458 base_bytes_per_node = mm_octree->base_bytes_per_node;
00459 nonbase_bytes_per_node = mm_octree->nonbase_bytes_per_node;
00460 child_offset = mm_octree->child_offset;
00461 pyr_ptr = (char *)mm_pyramid[level] + ((z*nodes_per_side + y) *
00462 nodes_per_side + x) * pyr_bytes_per_node;
00463
00464
00465 if (parent_node != NULL) {
00466 Debug((vpc, VPDEBUG_OCTREE,
00467 " Node at level %d, coords %d,%d,%d, addr 0x%08x\n",
00468 level, x, y, z, parent_node));
00469 for (p = 0; p < pyr_bytes_per_node; p++)
00470 ByteField(parent_node, p) = ByteField(pyr_ptr, p);
00471 }
00472
00473
00474 if (level < mm_octree->levels-1) {
00475
00476 subdivide = 0;
00477 for (p = 0; p < vpc->num_clsfy_params; p++) {
00478 field_size = vpc->field_size[vpc->param_field[p]];
00479 field_offset = mm_octree->node_offset[p];
00480 if (field_size == 1) {
00481 range = ByteField(pyr_ptr, field_offset+1) -
00482 ByteField(pyr_ptr, field_offset);
00483 } else {
00484 ASSERT(field_size == 2);
00485 range = ShortField(pyr_ptr, field_offset+2) -
00486 ShortField(pyr_ptr, field_offset);
00487 }
00488 if (range > vpc->param_maxrange[p]) {
00489 subdivide = 1;
00490 break;
00491 }
00492 }
00493
00494 if (subdivide) {
00495
00496 if (parent_node != NULL) {
00497 child_node = (char *)mm_octree->root + *octree_offset;
00498 IntField(parent_node, child_offset) = *octree_offset;
00499 Debug((vpc, VPDEBUG_OCTREE,
00500 " Storing children at offset = %d, addr = 0x%08x\n",
00501 *octree_offset, child_node));
00502 }
00503 if (level == mm_octree->levels-2)
00504 child_bytes_per_node = base_bytes_per_node;
00505 else
00506 child_bytes_per_node = nonbase_bytes_per_node;
00507 *octree_offset += 8 * child_bytes_per_node;
00508 if (parent_node == NULL) {
00509 child_node = NULL;
00510 child_bytes_per_node = 0;
00511 }
00512
00513
00514 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2, z*2,
00515 nodes_per_side*2, child_node, octree_offset);
00516 child_node += child_bytes_per_node;
00517 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2, z*2,
00518 nodes_per_side*2, child_node, octree_offset);
00519 child_node += child_bytes_per_node;
00520 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2+1, z*2,
00521 nodes_per_side*2, child_node, octree_offset);
00522 child_node += child_bytes_per_node;
00523 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2+1, z*2,
00524 nodes_per_side*2, child_node, octree_offset);
00525 child_node += child_bytes_per_node;
00526 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2, z*2+1,
00527 nodes_per_side*2, child_node, octree_offset);
00528 child_node += child_bytes_per_node;
00529 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1, y*2, z*2+1,
00530 nodes_per_side*2, child_node, octree_offset);
00531 child_node += child_bytes_per_node;
00532 DescendPyramid(vpc, mm_pyramid, level+1, x*2, y*2+1, z*2+1,
00533 nodes_per_side*2, child_node, octree_offset);
00534 child_node += child_bytes_per_node;
00535 DescendPyramid(vpc, mm_pyramid, level+1, x*2+1,y*2+1,z*2+1,
00536 nodes_per_side*2, child_node, octree_offset);
00537 } else {
00538
00539 Debug((vpc, VPDEBUG_OCTREE, " Not subdividing.\n"));
00540 if (parent_node != NULL) {
00541 IntField(parent_node, child_offset) = 0;
00542 }
00543 }
00544 }
00545 }
00546
00547
00548
00549
00550
00551
00552
00553 void
00554 VPComputeSummedAreaTable(vpc)
00555 vpContext *vpc;
00556 {
00557
00558
00559 switch (vpc->num_clsfy_params) {
00560 case 1:
00561 Compute1DSummedAreaTable(vpc);
00562 break;
00563 case 2:
00564 Compute2DSummedAreaTable(vpc);
00565 break;
00566 default:
00567
00568 VPBug("VPComputeSummedAreaTable can only handle 1D or 2D classifiers");
00569 break;
00570 }
00571 }
00572
00573
00574
00575
00576
00577
00578
00579 static void
00580 Compute1DSummedAreaTable(vpc)
00581 vpContext *vpc;
00582 {
00583 int p0max, p0value;
00584 unsigned table_size;
00585 float opacity, min_opacity, *p0table;
00586 unsigned sum;
00587 unsigned *entry;
00588
00589 p0max = vpc->field_max[vpc->param_field[0]];
00590 table_size = (p0max+1) * sizeof(unsigned);
00591 p0table = vpc->clsfy_table[0];
00592 min_opacity = vpc->min_opacity;
00593 if (vpc->sum_table == NULL || table_size != vpc->sum_table_size) {
00594 if (vpc->sum_table != NULL)
00595 Dealloc(vpc, vpc->sum_table);
00596 Alloc(vpc, vpc->sum_table, unsigned *, table_size, "sum_table");
00597 vpc->sum_table_size = table_size;
00598 }
00599 entry = vpc->sum_table;
00600 for (p0value = 0; p0value <= p0max; p0value++) {
00601 opacity = p0table[p0value];
00602 if (opacity > min_opacity)
00603 sum = 1;
00604 else
00605 sum = 0;
00606 if (p0value > 0)
00607 sum += entry[-1];
00608 entry[0] = sum;
00609 entry++;
00610 }
00611 }
00612
00613
00614
00615
00616
00617
00618
00619 static void
00620 Compute2DSummedAreaTable(vpc)
00621 vpContext *vpc;
00622 {
00623 int p0max, p0value, p1max, p1value;
00624 unsigned table_size;
00625 float opacity, min_opacity, *p0table, *p1table;
00626 unsigned sum;
00627 unsigned *entry;
00628
00629 p0max = vpc->field_max[vpc->param_field[0]];
00630 p1max = vpc->field_max[vpc->param_field[1]];
00631 table_size = (p0max+1) * (p1max+1) * sizeof(unsigned);
00632 p0table = vpc->clsfy_table[0];
00633 p1table = vpc->clsfy_table[1];
00634 min_opacity = vpc->min_opacity;
00635 if (vpc->sum_table == NULL || table_size != vpc->sum_table_size) {
00636 if (vpc->sum_table != NULL)
00637 Dealloc(vpc, vpc->sum_table);
00638 Alloc(vpc, vpc->sum_table, unsigned *, table_size, "sum_table");
00639 vpc->sum_table_size = table_size;
00640 }
00641 entry = vpc->sum_table;
00642 for (p0value = 0; p0value <= p0max; p0value++) {
00643 for (p1value = 0; p1value <= p1max; p1value++) {
00644 opacity = p0table[p0value] * p1table[p1value];
00645 if (opacity > min_opacity)
00646 sum = 1;
00647 else
00648 sum = 0;
00649 if (p1value > 0) {
00650 sum += entry[-1];
00651 if (p0value > 0) {
00652 sum += entry[-(p1max+1)];
00653 sum -= entry[-(p1max+1)-1];
00654 }
00655 } else if (p0value > 0) {
00656 sum += entry[-(p1max+1)];
00657 }
00658 entry[0] = sum;
00659 entry++;
00660 }
00661 }
00662 }
00663
00664
00665
00666
00667
00668
00669
00670 void
00671 VPClassifyOctree(vpc)
00672 vpContext *vpc;
00673 {
00674
00675
00676 switch (vpc->num_clsfy_params) {
00677 case 1:
00678 ClassifyOctree1(vpc);
00679 break;
00680 case 2:
00681 ClassifyOctree2(vpc);
00682 break;
00683 default:
00684
00685 VPBug("VPClassifyOctree can only handle 2D classifiers");
00686 break;
00687 }
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698 static void
00699 ClassifyOctree1(vpc)
00700 vpContext *vpc;
00701 {
00702 char *node_stack[VP_MAX_OCTREE_LEVELS];
00703 int count_stack[VP_MAX_OCTREE_LEVELS];
00704
00705 int level;
00706 int max_level;
00707 char *octree_root;
00708 char *node;
00709 unsigned area;
00710 int status;
00711 unsigned *sum_table;
00712 int p0max, p0min;
00713 int p0size;
00714 int child_offset;
00715 int status_offset;
00716 int base_bytes_per_node;
00717 int nonbase_bytes_per_node;
00718 int child_count;
00719
00720
00721 ASSERT(vpc->sum_table != NULL);
00722 ASSERT(vpc->mm_octree != NULL);
00723 ASSERT(vpc->mm_octree->root != NULL);
00724 ASSERT(vpc->sum_table_size == sizeof(unsigned) *
00725 (vpc->field_max[vpc->param_field[0]]+1));
00726 sum_table = vpc->sum_table;
00727 max_level = vpc->mm_octree->levels - 1;
00728 octree_root = vpc->mm_octree->root;
00729 p0size = vpc->field_size[vpc->param_field[0]];
00730 status_offset = vpc->mm_octree->status_offset;
00731 child_offset = vpc->mm_octree->child_offset;
00732 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
00733 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
00734 node = octree_root;
00735 level = 0;
00736
00737
00738 Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
00739 while (1) {
00740
00741 if (p0size == 1) {
00742 p0min = ByteField(node, 0)-1;
00743 p0max = ByteField(node, 1);
00744 } else {
00745 p0min = ShortField(node, 0)-1;
00746 p0max = ShortField(node, 2);
00747 }
00748
00749
00750 area = sum_table[p0max];
00751 if (p0min >= 0)
00752 area -= sum_table[p0min];
00753
00754
00755 if (area == 0) {
00756 status = MM_EMPTY;
00757 } else if (level != max_level && IntField(node, child_offset) != 0 &&
00758 area != (p0max - p0min)) {
00759 status = MM_PARTFULL;
00760 } else {
00761 status = MM_FULL;
00762 }
00763 ByteField(node, status_offset) = status;
00764 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00765 " Level %d: node is %s (addr 0x%08x)\n", level,
00766 status == MM_EMPTY ? "empty" : (status == MM_FULL ? "full" :
00767 "partfull"), node));
00768
00769
00770 if (status == MM_PARTFULL) {
00771
00772 node = octree_root + IntField(node, child_offset);
00773 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00774 " Descending. Children at offset %d, addr 0x%08x\n",
00775 IntField(node, child_offset), node));
00776 node_stack[level] = node;
00777 count_stack[level] = 7;
00778 level++;
00779 ASSERT(level <= max_level);
00780 } else {
00781 do {
00782
00783 Debug((vpc, VPDEBUG_CLSFYOCTREE, " Ascending.\n"));
00784 level--;
00785 if (level < 0)
00786 break;
00787 child_count = count_stack[level]--;
00788 ASSERT(child_count >= 0 && child_count <= 7);
00789 } while (child_count == 0);
00790 if (level < 0)
00791 break;
00792
00793
00794 if (level == max_level-1)
00795 node = node_stack[level] + base_bytes_per_node;
00796 else
00797 node = node_stack[level] + nonbase_bytes_per_node;
00798 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00799 " Descending to child at 0x%08x.\n", node));
00800 node_stack[level] = node;
00801 level++;
00802 ASSERT(level <= max_level);
00803 }
00804 }
00805 }
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815 static void
00816 ClassifyOctree2(vpc)
00817 vpContext *vpc;
00818 {
00819 char *node_stack[VP_MAX_OCTREE_LEVELS];
00820 int count_stack[VP_MAX_OCTREE_LEVELS];
00821
00822 int level;
00823 int max_level;
00824 char *octree_root;
00825 char *node;
00826 unsigned area;
00827 int status;
00828 unsigned *sum_table;
00829 int sum_table_dim1;
00830 int p0max, p0min;
00831 int p1max, p1min;
00832 int p0size, p1size;
00833 int child_offset;
00834 int status_offset;
00835 int base_bytes_per_node;
00836 int nonbase_bytes_per_node;
00837 int child_count;
00838
00839
00840 ASSERT(vpc->sum_table != NULL);
00841 ASSERT(vpc->mm_octree != NULL);
00842 ASSERT(vpc->mm_octree->root != NULL);
00843 ASSERT(vpc->sum_table_size == sizeof(unsigned) *
00844 (vpc->field_max[vpc->param_field[0]]+1) *
00845 (vpc->field_max[vpc->param_field[1]]+1));
00846 sum_table = vpc->sum_table;
00847 max_level = vpc->mm_octree->levels - 1;
00848 octree_root = vpc->mm_octree->root;
00849 p0size = vpc->field_size[vpc->param_field[0]];
00850 p1size = vpc->field_size[vpc->param_field[1]];
00851 sum_table_dim1 = vpc->field_max[vpc->param_field[1]] + 1;
00852 status_offset = vpc->mm_octree->status_offset;
00853 child_offset = vpc->mm_octree->child_offset;
00854 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
00855 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
00856 node = octree_root;
00857 level = 0;
00858
00859
00860 Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
00861 while (1) {
00862
00863 if (p0size == 1) {
00864 p0min = ByteField(node, 0)-1;
00865 p0max = ByteField(node, 1);
00866 } else {
00867 p0min = ShortField(node, 0)-1;
00868 p0max = ShortField(node, 2);
00869 }
00870 if (p1size == 1) {
00871 p1min = ByteField(node, 2*p0size)-1;
00872 p1max = ByteField(node, 2*p0size+1);
00873 } else {
00874 p1min = ShortField(node, 2*p0size)-1;
00875 p1max = ShortField(node, 2*p0size+2);
00876 }
00877
00878
00879 area = sum_table[p0max * sum_table_dim1 + p1max];
00880 if (p0min >= 0) {
00881 if (p1min >= 0) {
00882 area += sum_table[p0min * sum_table_dim1 + p1min];
00883 area -= sum_table[p0max * sum_table_dim1 + p1min];
00884 }
00885 area -= sum_table[p0min * sum_table_dim1 + p1max];
00886 } else {
00887 if (p1min >= 0)
00888 area -= sum_table[p0max * sum_table_dim1 + p1min];
00889 }
00890
00891
00892 if (area == 0) {
00893 status = MM_EMPTY;
00894 } else if (level != max_level && IntField(node, child_offset) != 0 &&
00895 area != (p1max - p1min)*(p0max - p0min)) {
00896 status = MM_PARTFULL;
00897 } else {
00898 status = MM_FULL;
00899 }
00900 ByteField(node, status_offset) = status;
00901 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00902 " Level %d: node is %s (addr 0x%08x)\n", level,
00903 status == MM_EMPTY ? "empty" : (status == MM_FULL ? "full" :
00904 "partfull"), node));
00905
00906
00907 if (status == MM_PARTFULL) {
00908
00909 node = octree_root + IntField(node, child_offset);
00910 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00911 " Descending. Children at offset %d, addr 0x%08x\n",
00912 IntField(node, child_offset), node));
00913 node_stack[level] = node;
00914 count_stack[level] = 7;
00915 level++;
00916 ASSERT(level <= max_level);
00917 } else {
00918 do {
00919
00920 Debug((vpc, VPDEBUG_CLSFYOCTREE, " Ascending.\n"));
00921 level--;
00922 if (level < 0)
00923 break;
00924 child_count = count_stack[level]--;
00925 ASSERT(child_count >= 0 && child_count <= 7);
00926 } while (child_count == 0);
00927 if (level < 0)
00928 break;
00929
00930
00931 if (level == max_level-1)
00932 node = node_stack[level] + base_bytes_per_node;
00933 else
00934 node = node_stack[level] + nonbase_bytes_per_node;
00935 Debug((vpc, VPDEBUG_CLSFYOCTREE,
00936 " Descending to child at 0x%08x.\n", node));
00937 node_stack[level] = node;
00938 level++;
00939 ASSERT(level <= max_level);
00940 }
00941 }
00942 }
00943
00944
00945
00946
00947
00948
00949
00950 void
00951 VPInitOctreeLevelStack(vpc, level_stack, axis, k)
00952 vpContext *vpc;
00953 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS];
00954 int axis;
00955 int k;
00956 {
00957 int max_level, level, last_node_size;
00958 int child_octant, child_bytes_per_node;
00959 int *octant_order;
00960
00961 ASSERT(vpc->mm_octree != NULL);
00962 max_level = vpc->mm_octree->levels-1;
00963 level_stack[max_level].level_size = vpc->mm_octree->base_node_size;
00964 level_stack[max_level].child_octant = -1;
00965 level_stack[max_level].child_offset1 = -1;
00966 level_stack[max_level].child_offset2 = -1;
00967 level_stack[max_level].child2 = NULL;
00968 last_node_size = vpc->mm_octree->base_node_size;
00969 octant_order = OctantOrder[axis];
00970 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Octants for next scanline:\n"));
00971 for (level = max_level-1; level >= 0; level--) {
00972 level_stack[level].level_size = last_node_size * 2;
00973 child_octant = ((k / last_node_size) & 1) * MM_K_BIT;
00974 last_node_size *= 2;
00975 level_stack[level].child_octant = child_octant;
00976 if (level == max_level-1)
00977 child_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
00978 else
00979 child_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
00980 ASSERT(child_octant >= 0 && child_octant < 7);
00981 ASSERT(octant_order[child_octant] >= 0 &&
00982 octant_order[child_octant] < 8);
00983 level_stack[level].child_offset1 = octant_order[child_octant] *
00984 child_bytes_per_node;
00985 level_stack[level].child_offset2 = octant_order[child_octant+1] *
00986 child_bytes_per_node;
00987 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Level %d: %d, then %d\n",
00988 level,octant_order[child_octant],octant_order[child_octant+1]));
00989 }
00990 }
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004 int
01005 VPComputeScanRuns(vpc, level_stack, run_lengths, axis, j, icount)
01006 vpContext *vpc;
01007 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS];
01008 unsigned char *run_lengths;
01009 int axis;
01010 int j;
01011 int icount;
01012 {
01013 int octree_maxlevel;
01014 int level;
01015 int max_level = vpc->mm_octree->levels-1;
01016 int child_octant, child_bytes_per_node;
01017 int base_bytes_per_node, nonbase_bytes_per_node;
01018 int i;
01019 char *octree_root, *node;
01020 int run_type;
01021 int run_length;
01022 int run_piece;
01023 int status_offset;
01024 int child_offset;
01025 int status;
01026 int *octant_order;
01027
01028 Debug((vpc, VPDEBUG_OCTREERUNS, "Runs for scanline %d:\n", j));
01029 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Traversal for scanline %d:\n", j));
01030 ASSERT(vpc->mm_octree != NULL);
01031 ASSERT(vpc->mm_octree->root != NULL);
01032 base_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
01033 nonbase_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
01034 status_offset = vpc->mm_octree->status_offset;
01035 child_offset = vpc->mm_octree->child_offset;
01036 octree_maxlevel = -1;
01037 i = icount;
01038 octree_root = vpc->mm_octree->root;
01039 node = octree_root;
01040 level = 0;
01041 run_type = MM_EMPTY;
01042 run_length = 0;
01043 octant_order = OctantOrder[axis];
01044
01045
01046 while (1) {
01047
01048 while (1) {
01049 status = ByteField(node, status_offset);
01050 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Node at level %d: %s\n",
01051 level, status == MM_PARTFULL ? "partfull" :
01052 (status == MM_FULL ? "full" : "empty")));
01053 ASSERT(status == MM_PARTFULL || status == MM_FULL ||
01054 status == MM_EMPTY);
01055 if (status != MM_PARTFULL)
01056 break;
01057 ASSERT(IntField(node, child_offset) != 0);
01058 Debug((vpc, VPDEBUG_OCTREETRAVERSE,
01059 " Children at base %d, offsets %d, %d; ",
01060 IntField(node, child_offset),
01061 level_stack[level].child_offset1,
01062 level_stack[level].child_offset2));
01063 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "status %d, %d\n",
01064 ByteField(octree_root + IntField(node, child_offset) +
01065 level_stack[level].child_offset1, status_offset),
01066 ByteField(octree_root + IntField(node, child_offset) +
01067 level_stack[level].child_offset2, status_offset)));
01068 node = octree_root + IntField(node, child_offset);
01069 level_stack[level].child2 = node+level_stack[level].child_offset2;
01070 node += level_stack[level].child_offset1;
01071 level++;
01072 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Descending.\n"));
01073 ASSERT(level < vpc->mm_octree->levels);
01074 }
01075 if (level > octree_maxlevel)
01076 octree_maxlevel = level;
01077
01078
01079 run_piece = MIN(level_stack[level].level_size, i);
01080 i -= run_piece;
01081 if (status == run_type) {
01082 run_length += run_piece;
01083 } else {
01084 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " New run.\n"));
01085 while (run_length > 255) {
01086 Debug((vpc, VPDEBUG_OCTREERUNS, " 255 0"));
01087 *run_lengths++ = 255;
01088 *run_lengths++ = 0;
01089 run_length -= 255;
01090 }
01091 Debug((vpc, VPDEBUG_OCTREERUNS, " %d", run_length));
01092 *run_lengths++ = run_length;
01093 run_type ^= 1;
01094 run_length = run_piece;
01095 }
01096 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Added %d to run.\n",
01097 run_piece));
01098 if (i == 0)
01099 break;
01100
01101
01102 do {
01103 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Ascending.\n"));
01104 level--;
01105 ASSERT(level >= 0);
01106 } while (level_stack[level].child2 == NULL);
01107
01108
01109 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Next child--descending.\n"));
01110 node = level_stack[level].child2;
01111 level_stack[level].child2 = NULL;
01112 level++;
01113 ASSERT(level < vpc->mm_octree->levels);
01114 }
01115
01116
01117 while (run_length > 255) {
01118 Debug((vpc, VPDEBUG_OCTREERUNS, " 255 0"));
01119 *run_lengths++ = 255;
01120 *run_lengths++ = 0;
01121 run_length -= 255;
01122 }
01123 Debug((vpc, VPDEBUG_OCTREERUNS, " %d", run_length));
01124 *run_lengths++ = run_length;
01125 if (run_type == MM_EMPTY) {
01126 Debug((vpc, VPDEBUG_OCTREERUNS, " 0"));
01127 *run_lengths++ = 0;
01128 }
01129 Debug((vpc, VPDEBUG_OCTREERUNS, "\n"));
01130
01131 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Covered %d scanlines.\n",
01132 level_stack[octree_maxlevel].level_size));
01133
01134
01135
01136 j += level_stack[octree_maxlevel].level_size;
01137 max_level = vpc->mm_octree->levels-1;
01138 Debug((vpc, VPDEBUG_OCTREETRAVERSE, "Octants for next scanline:\n"));
01139 for (level = max_level-1; level >= 0; level--) {
01140 child_octant = level_stack[level].child_octant;
01141 if (level >= octree_maxlevel)
01142 child_octant &= MM_K_BIT;
01143 else if ((j & (level_stack[level].level_size/2)) == 0)
01144 child_octant &= ~MM_J_BIT;
01145 else
01146 child_octant |= MM_J_BIT;
01147 level_stack[level].child_octant = child_octant;
01148
01149 if (level == max_level-1)
01150 child_bytes_per_node = base_bytes_per_node;
01151 else
01152 child_bytes_per_node = nonbase_bytes_per_node;
01153 level_stack[level].child_offset1 = octant_order[child_octant] *
01154 child_bytes_per_node;
01155 level_stack[level].child_offset2 = octant_order[child_octant+1] *
01156 child_bytes_per_node;
01157 Debug((vpc, VPDEBUG_OCTREETRAVERSE, " Level %d: %d, then %d\n",
01158 level,octant_order[child_octant],octant_order[child_octant+1]));
01159 }
01160
01161
01162
01163 return(level_stack[octree_maxlevel].level_size);
01164 }
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179 vpResult
01180 vpOctreeMask(vpc, array, array_size, max_level)
01181 vpContext *vpc;
01182 unsigned char *array;
01183 int array_size;
01184 int max_level;
01185 {
01186 int c;
01187 unsigned char *aptr;
01188 int retcode;
01189
01190
01191 if (vpc->mm_octree == NULL)
01192 return(VPSetError(vpc, VPERROR_BAD_SIZE));
01193 if ((retcode = VPCheckClassifier(vpc)) == NULL)
01194 return(retcode);
01195 if (array_size != vpc->xlen*vpc->ylen*vpc->zlen)
01196 return(VPSetError(vpc, VPERROR_BAD_SIZE));
01197
01198
01199 VPComputeSummedAreaTable(vpc);
01200 VPClassifyOctree(vpc);
01201 ComputeOctreeMask(vpc, 0, 0, 0, 0, vpc->mm_octree->root_node_size,
01202 vpc->mm_octree->root, array, max_level);
01203 return(VP_OK);
01204 }
01205
01206
01207
01208
01209
01210
01211
01212 static void
01213 ComputeOctreeMask(vpc, level, xn, yn, zn, node_size, parent_node, array,
01214 max_level)
01215 vpContext *vpc;
01216 int level;
01217 int xn, yn, zn;
01218
01219 int node_size;
01220 void *parent_node;
01221 unsigned char *array;
01222 int max_level;
01223 {
01224 char *child_node, *octree_root;
01225 int child_bytes_per_node;
01226 int child_offset;
01227 int status_offset;
01228 int status, value;
01229 int x, y, z, x0, y0, z0, x1, y1, z1;
01230 int array_ystride, array_zstride;
01231
01232
01233 status_offset = vpc->mm_octree->status_offset;
01234 child_offset = vpc->mm_octree->child_offset;
01235 if (level == vpc->mm_octree->levels-2)
01236 child_bytes_per_node = vpc->mm_octree->base_bytes_per_node;
01237 else
01238 child_bytes_per_node = vpc->mm_octree->nonbase_bytes_per_node;
01239 octree_root = vpc->mm_octree->root;
01240
01241
01242 status = ByteField(parent_node, status_offset);
01243 if (level == max_level || status != MM_PARTFULL) {
01244 if (status == MM_EMPTY)
01245 value = 0;
01246 else if (status == MM_FULL)
01247 value = 255;
01248 else if (status == MM_PARTFULL)
01249 value = 128;
01250 else
01251 VPBug("bad status value in ComputeOctreeMask, nodeaddr = 0x%08x",
01252 parent_node);
01253 x0 = xn * node_size;
01254 y0 = yn * node_size;
01255 z0 = zn * node_size;
01256 x1 = MIN(x0 + node_size, vpc->xlen) - 1;
01257 y1 = MIN(y0 + node_size, vpc->ylen) - 1;
01258 z1 = MIN(z0 + node_size, vpc->zlen) - 1;
01259 array_ystride = vpc->xlen;
01260 array_zstride = vpc->xlen * vpc->ylen;
01261 for (z = z0; z <= z1; z++) {
01262 for (y = y0; y <= y1; y++) {
01263 for (x = x0; x <= x1; x++) {
01264 array[z*array_zstride + y*array_ystride + x] = value;
01265 }
01266 }
01267 }
01268 return;
01269 }
01270 ASSERT(IntField(parent_node, child_offset) != 0);
01271
01272
01273 child_node = octree_root + IntField(parent_node, child_offset);
01274 ComputeOctreeMask(vpc, level+1, xn*2, yn*2, zn*2, node_size/2,
01275 child_node, array, max_level);
01276 child_node += child_bytes_per_node;
01277 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2, zn*2, node_size/2,
01278 child_node, array, max_level);
01279 child_node += child_bytes_per_node;
01280 ComputeOctreeMask(vpc, level+1, xn*2, yn*2+1, zn*2, node_size/2,
01281 child_node, array, max_level);
01282 child_node += child_bytes_per_node;
01283 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2+1, zn*2, node_size/2,
01284 child_node, array, max_level);
01285 child_node += child_bytes_per_node;
01286 ComputeOctreeMask(vpc, level+1, xn*2, yn*2, zn*2+1, node_size/2,
01287 child_node, array, max_level);
01288 child_node += child_bytes_per_node;
01289 ComputeOctreeMask(vpc, level+1, xn*2+1, yn*2, zn*2+1, node_size/2,
01290 child_node, array, max_level);
01291 child_node += child_bytes_per_node;
01292 ComputeOctreeMask(vpc, level+1, xn*2, yn*2+1, zn*2+1, node_size/2,
01293 child_node, array, max_level);
01294 child_node += child_bytes_per_node;
01295 ComputeOctreeMask(vpc, level+1, xn*2+1,yn*2+1,zn*2+1, node_size/2,
01296 child_node, array, max_level);
01297 }
01298
01299 #ifdef DEBUG
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309 int
01310 VPCheckRuns(vpc, run_lengths, axis, k, j)
01311 vpContext *vpc;
01312 unsigned char *run_lengths;
01313 int axis;
01314 int k;
01315 int j;
01316 {
01317 char *voxel;
01318 int i;
01319 int icount;
01320 int count;
01321 int is_non_zero;
01322 int num_runs;
01323 float opacity;
01324 int badpredictions;
01325 int istride;
01326
01327 switch (axis) {
01328 case VP_X_AXIS:
01329 voxel = (char *)vpc->raw_voxels + k*vpc->xstride + j*vpc->zstride;
01330 istride = vpc->ystride;
01331 icount = vpc->ylen;
01332 break;
01333 case VP_Y_AXIS:
01334 voxel = (char *)vpc->raw_voxels + k*vpc->ystride + j*vpc->xstride;
01335 istride = vpc->zstride;
01336 icount = vpc->zlen;
01337 break;
01338 case VP_Z_AXIS:
01339 voxel = (char *)vpc->raw_voxels + k*vpc->zstride + j*vpc->ystride;
01340 istride = vpc->xstride;
01341 icount = vpc->xlen;
01342 break;
01343 default:
01344 VPBug("bad axis in VPCheckRuns");
01345 }
01346
01347 count = 0;
01348 is_non_zero = 1;
01349 num_runs = 0;
01350 badpredictions = 0;
01351 for (i = 0; i < icount; i++) {
01352 while (count == 0) {
01353 count = *run_lengths++;
01354 is_non_zero = !is_non_zero;
01355 if (++num_runs > icount)
01356 VPBug("runaway scanline detected by VPCheckRuns");
01357 }
01358 opacity = VPClassifyVoxel(vpc, voxel);
01359 if (opacity > vpc->min_opacity &&
01360 fabs(opacity - vpc->min_opacity) > 0.001) {
01361 if (!is_non_zero) {
01362 printf("\n");
01363 printf("VPCheckRuns: error on voxel (i,j,k)=(%d,%d,%d), ",
01364 i, j, k);
01365 printf("viewaxis %d\n", axis);
01366 printf("Actual opacity: %17.15f\n", opacity);
01367 printf("Threshold: %17.15f\n", vpc->min_opacity);
01368 VPDumpView(vpc);
01369 VPDumpClassifier(vpc);
01370 VPBug("nonzero voxel in zero run detected by VPCheckRuns");
01371 }
01372 } else {
01373 if (is_non_zero)
01374 badpredictions++;
01375 }
01376 voxel += istride;
01377 count--;
01378 }
01379 if (count != 0)
01380 VPBug("run that overshoots end of scanline detected by VPCheckRuns");
01381 if (!is_non_zero) {
01382 if (*run_lengths != 0)
01383 VPBug("missing 0 run at end of scanline detected by VPCheckRuns");
01384 }
01385 return(badpredictions);
01386 }
01387
01388
01389
01390
01391
01392
01393
01394 void
01395 VPTestMinMaxOctree(vpc)
01396 vpContext *vpc;
01397 {
01398 int x, y, z;
01399 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS];
01400 unsigned char run_lengths[VP_MAX_VOLUME_DIM];
01401 int badpredictions;
01402 int scans_left;
01403 float accuracy;
01404
01405 ASSERT(vpc->mm_octree != NULL);
01406 VPComputeSummedAreaTable(vpc);
01407 VPClassifyOctree(vpc);
01408
01409 badpredictions = 0;
01410 printf("Checking +Z axis runs...\n");
01411 for (z = 0; z < vpc->zlen; z++) {
01412 ReportStatus(vpc, (double)z / (double)vpc->zlen);
01413 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", z));
01414 VPInitOctreeLevelStack(vpc, level_stack, VP_Z_AXIS, z);
01415 scans_left = 0;
01416 for (y = 0; y < vpc->ylen; y++) {
01417 if (scans_left == 0) {
01418 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01419 VP_Z_AXIS, y, vpc->xlen);
01420 }
01421 scans_left--;
01422 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Z_AXIS, z, y);
01423 }
01424 }
01425 ReportStatus(vpc, 1.0);
01426 accuracy = 1. - ((double)badpredictions /
01427 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01428 printf("VPTestMinMaxOctree: PASSED.\n");
01429 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01430 accuracy*100., badpredictions);
01431
01432 badpredictions = 0;
01433 printf("Checking +Y axis runs...\n");
01434 for (y = 0; y < vpc->ylen; y++) {
01435 ReportStatus(vpc, (double)y / (double)vpc->ylen);
01436 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", y));
01437 VPInitOctreeLevelStack(vpc, level_stack, VP_Y_AXIS, y);
01438 scans_left = 0;
01439 for (x = 0; x < vpc->xlen; x++) {
01440 if (scans_left == 0) {
01441 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01442 VP_Y_AXIS, x, vpc->zlen);
01443 }
01444 scans_left--;
01445 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Y_AXIS, y, x);
01446 }
01447 }
01448 ReportStatus(vpc, 1.0);
01449 accuracy = 1. - ((double)badpredictions /
01450 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01451 printf("VPTestMinMaxOctree: PASSED.\n");
01452 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01453 accuracy*100., badpredictions);
01454
01455 badpredictions = 0;
01456 printf("Checking +X axis runs...\n");
01457 for (x = 0; x < vpc->xlen; x++) {
01458 ReportStatus(vpc, (double)x / (double)vpc->xlen);
01459 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", x));
01460 VPInitOctreeLevelStack(vpc, level_stack, VP_X_AXIS, x);
01461 scans_left = 0;
01462 for (z = 0; z < vpc->zlen; z++) {
01463 if (scans_left == 0) {
01464 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01465 VP_X_AXIS, z, vpc->ylen);
01466 }
01467 scans_left--;
01468 badpredictions += VPCheckRuns(vpc, run_lengths, VP_X_AXIS, x, z);
01469 }
01470 }
01471 ReportStatus(vpc, 1.0);
01472 accuracy = 1. - ((double)badpredictions /
01473 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01474 printf("VPTestMinMaxOctree: PASSED.\n");
01475 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01476 accuracy*100., badpredictions);
01477
01478 badpredictions = 0;
01479 printf("Checking -Z axis runs...\n");
01480 for (z = vpc->zlen-1; z >= 0; z--) {
01481 ReportStatus(vpc, (double)(vpc->zlen-1-z) / (double)vpc->zlen);
01482 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", z));
01483 VPInitOctreeLevelStack(vpc, level_stack, VP_Z_AXIS, z);
01484 scans_left = 0;
01485 for (y = 0; y < vpc->ylen; y++) {
01486 if (scans_left == 0) {
01487 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01488 VP_Z_AXIS, y, vpc->xlen);
01489 }
01490 scans_left--;
01491 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Z_AXIS, z, y);
01492 }
01493 }
01494 ReportStatus(vpc, 1.0);
01495 accuracy = 1. - ((double)badpredictions /
01496 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01497 printf("VPTestMinMaxOctree: PASSED.\n");
01498 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01499 accuracy*100., badpredictions);
01500
01501 badpredictions = 0;
01502 printf("Checking -Y axis runs...\n");
01503 for (y = vpc->ylen-1; y >= 0; y--) {
01504 ReportStatus(vpc, (double)(vpc->ylen-1-y) / (double)vpc->ylen);
01505 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", y));
01506 VPInitOctreeLevelStack(vpc, level_stack, VP_Y_AXIS, y);
01507 scans_left = 0;
01508 for (x = 0; x < vpc->xlen; x++) {
01509 if (scans_left == 0) {
01510 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01511 VP_Y_AXIS, x, vpc->zlen);
01512 }
01513 scans_left--;
01514 badpredictions += VPCheckRuns(vpc, run_lengths, VP_Y_AXIS, y, x);
01515 }
01516 }
01517 ReportStatus(vpc, 1.0);
01518 accuracy = 1. - ((double)badpredictions /
01519 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01520 printf("VPTestMinMaxOctree: PASSED.\n");
01521 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01522 accuracy*100., badpredictions);
01523
01524 badpredictions = 0;
01525 printf("Checking -X axis runs...\n");
01526 for (x = vpc->xlen-1; x >= 0; x--) {
01527 ReportStatus(vpc, (double)(vpc->xlen-1-x) / (double)vpc->xlen);
01528 Debug((vpc, VPDEBUG_OCTREERUNS, "*** Slice %d ***\n", x));
01529 VPInitOctreeLevelStack(vpc, level_stack, VP_X_AXIS, x);
01530 scans_left = 0;
01531 for (z = 0; z < vpc->zlen; z++) {
01532 if (scans_left == 0) {
01533 scans_left = VPComputeScanRuns(vpc, level_stack, run_lengths,
01534 VP_X_AXIS, z, vpc->ylen);
01535 }
01536 scans_left--;
01537 badpredictions += VPCheckRuns(vpc, run_lengths, VP_X_AXIS, x, z);
01538 }
01539 }
01540 ReportStatus(vpc, 1.0);
01541 accuracy = 1. - ((double)badpredictions /
01542 (double)(vpc->xlen*vpc->ylen*vpc->zlen));
01543 printf("VPTestMinMaxOctree: PASSED.\n");
01544 printf("Prediction accuracy: %.1f%% (%d bad predictions)\n",
01545 accuracy*100., badpredictions);
01546 }
01547 #else
01548
01549 int
01550 VPCheckRuns(vpc, run_lengths, axis, k, j)
01551 vpContext *vpc;
01552 unsigned char *run_lengths;
01553 int axis;
01554 int k;
01555 int j;
01556 {
01557 return 0 ;
01558 }
01559
01560 void
01561 VPTestMinMaxOctree(vpc)
01562 vpContext *vpc;
01563 {
01564 }
01565
01566 #endif