00001 /*
00002  * vp_octree.c
00003  *
00004  * Routines to create and destroy MinMaxOctree objects for fast classification.
00005  *
00006  * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
00007  * Junior University.  All rights reserved.
00008  *
00009  * Permission to use, copy, modify and distribute this software and its
00010  * documentation for any purpose is hereby granted without fee, provided
00011  * that the above copyright notice and this permission notice appear in
00012  * all copies of this software and that you do not sell the software.
00013  * Commercial licensing is available by contacting the author.
00014  * 
00015  * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
00016  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
00017  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
00018  *
00019  * Author:
00020  *    Phil Lacroute
00021  *    Computer Systems Laboratory
00022  *    Electrical Engineering Dept.
00023  *    Stanford University
00024  */
00025 
00026 /*
00027  * $Date: 2001/12/17 16:16:23 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 /*
00034  * OctantOrder: octant traversal order, depending on best_view_axis
00035  */
00036 
00037 static int OctantOrder[3][8] = {
00038     { 0, 2, 4, 6, 1, 3, 5, 7 }, /* VP_X_AXIS */
00039     { 0, 4, 1, 5, 2, 6, 3, 7 }, /* VP_Y_AXIS */
00040     { 0, 1, 2, 3, 4, 5, 6, 7 }  /* VP_Z_AXIS */
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  * vpCreateMinMaxOctree
00059  *
00060  * Create a MinMaxOctree representation of the volume for fast classification.
00061  */
00062 
00063 vpResult
00064 vpCreateMinMaxOctree(vpc, root_node_size, base_node_size)
00065 vpContext *vpc;
00066 int root_node_size;     /* ignored for now */
00067 int base_node_size;     /* controls level of detail of smallest nodes */
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     /* check for errors */
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) {   /* must be a power of 2 */
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     /* allocate mm_octree structure */
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     /* compute field sizes */
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++;                           /* add byte for status field */
00130     bytes_per_node = (bytes_per_node + 1) & ~1; /* align to short boundary */
00131     vpc->mm_octree->base_bytes_per_node = bytes_per_node;
00132     bytes_per_node = (bytes_per_node + 3) & ~3; /* align to word boundary */
00133     vpc->mm_octree->child_offset = bytes_per_node;
00134     bytes_per_node += sizeof(unsigned);         /* add word for child field */
00135     vpc->mm_octree->nonbase_bytes_per_node = bytes_per_node;
00136 
00137     /* compute number of levels */
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     /* build a min-max pyramid representation of the volume */
00152     CreatePyramid(vpc, mm_pyramid);
00153 
00154     /* determine how much space is needed for the octree nodes */
00155     octree_offset = vpc->mm_octree->nonbase_bytes_per_node; /* root node */
00156     DescendPyramid(vpc, mm_pyramid, 0, 0, 0, 0, 1, NULL, &octree_offset);
00157 
00158     /* create the min-max octree nodes */
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     /* clean up and return */
00167     Dealloc(vpc, mm_pyramid[0]);
00168     return(VP_OK);
00169 }
00170 
00171 /*
00172  * vpDestroyMinMaxOctree
00173  *
00174  * Destroy the MinMaxOctree representation of the volume.
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  * CreatePyramid
00194  *
00195  * Create a min-max pyramid representation of the volume.
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;               /* size of pyramid in bytes */
00204     int level, pyr_levels;      /* current, total pyramid levels */
00205     int level_offset;           /* byte offset to beginning of level */
00206     int nodes_per_side;         /* nodes per side at current level */
00207     int node_size;              /* voxels per side in node */
00208     char *pyr_node;             /* current node of pyramid */
00209     char *pyr_src;              /* pyramid node being read */
00210     char *pyr_dst;              /* pyramid node being written */
00211     char *voxel;                /* voxel being read */
00212     int x, y, z;                /* coordinates of current pyramid node */
00213     int nx, ny, nz;             /* coordinates of voxel within node */
00214     int xlen, ylen, zlen;       /* size of volume */
00215     int voxel_xstride;          /* volume strides */
00216     int voxel_ystride;
00217     int voxel_zstride;
00218     int param_size[VP_MAX_FIELDS]; /* size of each parameter */
00219     int param_offset[VP_MAX_FIELDS];/* voxel offset of each parameter */
00220     int node_offset[VP_MAX_FIELDS]; /* offset of parameter in octree node */
00221     int max_value[VP_MAX_FIELDS]; /* max. value of each parameter in node */
00222     int min_value[VP_MAX_FIELDS]; /* min. value of each parameter in node */
00223     int num_params;             /* number of params for classifier */
00224     int p;                      /* parameter number */
00225     int value;                  /* parameter value */
00226     int pyr_bytes_per_node;     /* size of node in bytes */
00227     int pyr_offsets[8];         /* offsets from pyr_src to each of its
00228                                    neighbors (0,+X,+Y,+XY,+Z,+XZ,+YZ,+XYZ) */
00229     int elem;                   /* index into pyr_offsets */
00230 
00231     /* allocate space for pyramid */
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     /* build the base level of the pyramid */
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                 /* clear the min/max values for the current node */
00274                 for (p = 0; p < num_params; p++) {
00275                     max_value[p] = -1;
00276                     min_value[p] = 65536;
00277                 }
00278 
00279                 /* loop over voxels in the node */
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                             /* compare each field against current min/max */
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                         } /* for nx */
00318                         voxel += voxel_ystride - node_size*voxel_xstride;
00319                     } /* for ny */
00320                     voxel += voxel_zstride - node_size*voxel_ystride;
00321                 } /* for nz */
00322 
00323                 /* store the min/max values for this node */
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             } /* for x */
00343             voxel += node_size*(voxel_ystride - nodes_per_side*voxel_xstride);
00344         } /* for y */
00345         voxel += node_size*(voxel_zstride - nodes_per_side*voxel_ystride);
00346     } /* for z */
00347     ReportStatus(vpc, 1.0);
00348 
00349     /* build the rest of the pyramid */
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                     /* clear the min/max values for the current node */
00369                     for (p = 0; p < num_params; p++) {
00370                         max_value[p] = -1;
00371                         min_value[p] = 65536;
00372                     }
00373 
00374                     /* loop over the eight children of this node */
00375                     for (elem = 0; elem < 8; elem++) {
00376                         pyr_src = pyr_node + pyr_offsets[elem];
00377                         /* compare min/max values of children with current
00378                            min/max values for the node */
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                     /* store the min/max values for this node */
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                     /* go on to the next node */
00410                     pyr_dst += pyr_bytes_per_node;
00411                     pyr_node += 2*pyr_bytes_per_node;
00412                 } /* for x */
00413                 pyr_node += (2*nodes_per_side)*pyr_bytes_per_node;
00414             } /* for y */
00415             pyr_node += (2*nodes_per_side)*(2*nodes_per_side)*
00416                         pyr_bytes_per_node;
00417         } /* for z */
00418     } /* for level */
00419     ReportStatus(vpc, 1.0);
00420 }
00421 
00422 /*
00423  * DescendPyramid
00424  *
00425  * Descend the pyramid recursively, either to count how many nodes will
00426  * be copied to the octree (if parent_node == NULL) or to actually copy them.
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;         /* context */
00433 void *mm_pyramid[VP_MAX_OCTREE_LEVELS]; /* min-max pyramid */
00434 int level;              /* current level */
00435 int x, y, z;            /* current node coordinates (in coordinate system
00436                            of the current level) */
00437 int nodes_per_side;     /* # nodes at current level per side of volume */
00438 void *parent_node;      /* parent octree node (or NULL) */
00439 int *octree_offset;     /* bytes from root of octree to next free location */
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     /* copy min/max data from pyramid node to octree node */
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     /* descend to next level */
00474     if (level < mm_octree->levels-1) {
00475         /* check if we should subdivide node or not */
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             /* store offset to child */
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             /* visit children */
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             /* node has no children; store NULL pointer */
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  * VPComputeSummedAreaTable
00549  *
00550  * Build the summed-area table for fast-classification.
00551  */
00552 
00553 void
00554 VPComputeSummedAreaTable(vpc)
00555 vpContext *vpc;
00556 {
00557     /* use a special-case version for lower dimensions
00558        (faster since C optimizer does a better job) */
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         /* XXX add code for ND classifiers */
00568         VPBug("VPComputeSummedAreaTable can only handle 1D or 2D classifiers");
00569         break;
00570     }
00571 }
00572 
00573 /*
00574  * Compute1DSummedAreaTable
00575  *
00576  * Build a 1D summed area table.
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  * Compute2DSummedAreaTable
00615  *
00616  * Build a 2D summed area table.
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  * VPClassifyOctree
00666  *
00667  * Descend an octree and classify each node as full, empty or partfull.
00668  */
00669 
00670 void
00671 VPClassifyOctree(vpc)
00672 vpContext *vpc;
00673 {
00674     /* use a special-case version for lower dimensions
00675        (faster since C optimizer does a better job) */
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         /* XXX add code for ND classifiers */
00685         VPBug("VPClassifyOctree can only handle 2D classifiers");
00686         break;
00687     }
00688 }
00689 
00690 /*
00691  * ClassifyOctree1
00692  *
00693  * Descend an octree and classify each node as full, empty or partfull.
00694  * Specialized for a 1 parameter classification function (1D summed
00695  * area table).
00696  */
00697 
00698 static void
00699 ClassifyOctree1(vpc)
00700 vpContext *vpc;
00701 {
00702     char *node_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node addresses */
00703     int count_stack[VP_MAX_OCTREE_LEVELS];  /* stack of node child counts;
00704                                    when count drops to zero, pop up a level */
00705     int level;                  /* current octree level */
00706     int max_level;              /* highest octree level */
00707     char *octree_root;          /* root node of octree */
00708     char *node;                 /* current octree node */
00709     unsigned area;              /* area computed from the summed-area table */
00710     int status;                 /* classification status of current node */
00711     unsigned *sum_table;        /* summed area table */
00712     int p0max, p0min;           /* parameter 0 extrema */
00713     int p0size;                 /* parameter size */
00714     int child_offset;           /* offset of child field in node */
00715     int status_offset;          /* offset of status field in node */
00716     int base_bytes_per_node;    /* size of base node in bytes */
00717     int nonbase_bytes_per_node; /* size of nonbase node in bytes */
00718     int child_count;            /* children left at current level */
00719 
00720     /* initialize */
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     /* do a depth-first, preorder traversal of the octree */
00738     Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
00739     while (1) {
00740         /* find min/max values for both parameters in this node */
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         /* integrate the opacities in the node using the summed area table */
00750         area = sum_table[p0max];
00751         if (p0min >= 0)
00752             area -= sum_table[p0min];
00753 
00754         /* decide if node is full, empty or partfull */
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         /* move to next node in tree traversal */
00770         if (status == MM_PARTFULL) {
00771             /* move down to first child in next level */
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;     /* number of remaining children */
00778             level++;
00779             ASSERT(level <= max_level);
00780         } else {
00781             do {
00782                 /* move up to a node with unvisited children */
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;  /* traversal of octree is done! */
00792 
00793             /* descend to the next child of this node */
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     } /* while (1) */
00805 }
00806 
00807 /*
00808  * ClassifyOctree2
00809  *
00810  * Descend an octree and classify each node as full, empty or partfull.
00811  * Specialized for a 2 parameter classification function (2D summed
00812  * area table).
00813  */
00814 
00815 static void
00816 ClassifyOctree2(vpc)
00817 vpContext *vpc;
00818 {
00819     char *node_stack[VP_MAX_OCTREE_LEVELS]; /* stack of node addresses */
00820     int count_stack[VP_MAX_OCTREE_LEVELS];  /* stack of node child counts;
00821                                    when count drops to zero, pop up a level */
00822     int level;                  /* current octree level */
00823     int max_level;              /* highest octree level */
00824     char *octree_root;          /* root node of octree */
00825     char *node;                 /* current octree node */
00826     unsigned area;              /* area computed from the summed-area table */
00827     int status;                 /* classification status of current node */
00828     unsigned *sum_table;        /* summed area table */
00829     int sum_table_dim1;         /* size of last dimension of sum_table */
00830     int p0max, p0min;           /* parameter 0 extrema */
00831     int p1max, p1min;           /* parameter 1 extrema */
00832     int p0size, p1size;         /* parameter sizes */
00833     int child_offset;           /* offset of child field in node */
00834     int status_offset;          /* offset of status field in node */
00835     int base_bytes_per_node;    /* size of base node in bytes */
00836     int nonbase_bytes_per_node; /* size of nonbase node in bytes */
00837     int child_count;            /* children left at current level */
00838 
00839     /* initialize */
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     /* do a depth-first, preorder traversal of the octree */
00860     Debug((vpc, VPDEBUG_CLSFYOCTREE, "Classifying octree:\n"));
00861     while (1) {
00862         /* find min/max values for both parameters in this node */
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         /* integrate the opacities in the node using the summed area table */
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         /* decide if node is full, empty or partfull */
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         /* move to next node in tree traversal */
00907         if (status == MM_PARTFULL) {
00908             /* move down to first child in next level */
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;     /* number of remaining children */
00915             level++;
00916             ASSERT(level <= max_level);
00917         } else {
00918             do {
00919                 /* move up to a node with unvisited children */
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;  /* traversal of octree is done! */
00929 
00930             /* descend to the next child of this node */
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     } /* while (1) */
00942 }
00943 
00944 /*
00945  * VPInitOctreeLevelStack
00946  *
00947  * Initialize an MMOctreeLevel stack.
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;       /* principle viewing axis */
00955 int k;          /* current slice number */
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  * VPComputeScanRuns
00994  *
00995  * For a given voxel scanline, produce a sequence of run lengths
00996  * which give a conservative estimate of the non-transparent portions
00997  * of the scanline.  The runs are computed by finding which nodes
00998  * of the classified min-max octree are intersected by the scanline.
00999  *
01000  * The return value is the number of scanlines for which this run data
01001  * is valid.
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]; /* saved state */
01008 unsigned char *run_lengths; /* storage for run lengths */
01009 int axis;               /* principle viewing axis */
01010 int j;                  /* scanline number */
01011 int icount;             /* scanline length */
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     /* traverse the octree */
01046     while (1) {
01047         /* descend tree to next node which is not partfull */
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         /* add current node to the list of runs */
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;      /* traversal is done */
01100 
01101         /* move back up the tree to the next node with unvisited children */
01102         do {
01103             Debug((vpc, VPDEBUG_OCTREETRAVERSE, "    Ascending.\n"));
01104             level--;
01105             ASSERT(level >= 0);
01106         } while (level_stack[level].child2 == NULL);
01107 
01108         /* descend to next child */
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     } /* while (1) */
01115 
01116     /* write out the last run */
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     /* update state for next scanline: adjust child_octant and then
01135        use it to compute child_offset1 and child_offset2 */
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     /* return the number of scanlines for which the run lengths are valid
01162        (which is the size of the smallest octree node the scanline hit) */
01163     return(level_stack[octree_maxlevel].level_size);
01164 }
01165 
01166 /*
01167  * vpOctreeMask
01168  *
01169  * Fill a 3D array with a mask computed from an octree. 
01170  * Each array element is set to one of three values depending upon
01171  * the value of the corresponding voxel in the octree:
01172  *    0         voxel is definitely transparent
01173  *    255       voxel may be non-transparent
01174  *    128       voxel may be non-transparent, and more detailed information
01175  *              is available at deeper levels of the octree which were not
01176  *              visited
01177  */
01178 
01179 vpResult
01180 vpOctreeMask(vpc, array, array_size, max_level)
01181 vpContext *vpc;         /* context */
01182 unsigned char *array;   /* array for result */
01183 int array_size;         /* size of array in bytes */
01184 int max_level;
01185 {
01186     int c;
01187     unsigned char *aptr;
01188     int retcode;
01189 
01190     /* error checks */
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     /* classify the octree */
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  * ComputeOctreeMask
01208  *
01209  * Recursive helper function for vpOctreeMask.
01210  */
01211 
01212 static void
01213 ComputeOctreeMask(vpc, level, xn, yn, zn, node_size, parent_node, array,
01214                   max_level)
01215 vpContext *vpc;         /* context */
01216 int level;              /* current level */
01217 int xn, yn, zn;         /* current node coordinates (in coordinate system
01218                            of the current level) */
01219 int node_size;          /* voxel per side of node at this level */
01220 void *parent_node;      /* parent octree node */
01221 unsigned char *array;   /* array for storing result */
01222 int max_level;          /* deepest level of the tree to visit */
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     /* initialize */
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     /* base case */
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     /* visit children */
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  * VPCheckRuns
01302  *
01303  * Check a scanline of run lengths for validity by comparing it to
01304  * the raw volume data.  Return value is the number of voxels in
01305  * nonzero runs which are actually zero (due to conservative
01306  * approximations.)  If an error is detected then VPBug is called.
01307  */
01308 
01309 int
01310 VPCheckRuns(vpc, run_lengths, axis, k, j)
01311 vpContext *vpc;
01312 unsigned char *run_lengths;/* run lengths */
01313 int axis;               /* principle viewing axis */
01314 int k;                  /* slice number */
01315 int j;                  /* scanline number */
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  * VPTestMinMaxOctree
01390  *
01391  * Test out the MinMaxOctree routines.
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;/* run lengths */
01553 int axis;               /* principle viewing axis */
01554 int k;                  /* slice number */
01555 int j;                  /* scanline number */
01556 {
01557   return 0 ;  /* RWCox */
01558 }
01559 
01560 void
01561 VPTestMinMaxOctree(vpc)
01562 vpContext *vpc;
01563 {
01564 }
01565 
01566 #endif /* DEBUG */