00001 /*
00002  * vp_rle.c
00003  *
00004  * Routines for run-length encoding classified volume data.
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 static void EncodeScanline ANSI_ARGS((vpContext *vpc, void *voxels,
00034     MinMaxOctree *octree, MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS]));
00035 static void InitRLE ANSI_ARGS((vpContext *vpc));
00036 static void CacheVoxel ANSI_ARGS((vpContext *vpc, double opacity,
00037     void *rawvoxel));
00038 static void CacheLength ANSI_ARGS((vpContext *vpc, int length));
00039 static void CountNonZeroVoxel ANSI_ARGS((RunData *rundata, int index,
00040     int end_of_scan, vpContext *vpc));
00041 static void RepackClassifiedVolume ANSI_ARGS((vpContext *vpc));
00042 static RLEVoxels *CreateEmptyRLEVoxels ANSI_ARGS((vpContext *vpc,
00043     int ilen, int jlen, int klen));
00044 static RLEVoxels *CreateRLEVoxelsFromRunData ANSI_ARGS((vpContext *vpc,
00045     int ilen, int jlen, int klen, int non_zero_count, RunData *run_data,
00046     int rle_bytes_per_voxel));
00047 static void StoreNonZeroVoxel ANSI_ARGS((void *voxel, int rle_bytes_per_voxel,
00048     void *data, unsigned char *lengths, RunData *rundata, int index));
00049 static void PadScanlines ANSI_ARGS((int ilen, int jlen, int klen,
00050     RunData *run_data, unsigned char *lengths));
00051 static ConstructionBuffer *CreateConstructionBuffer ANSI_ARGS((
00052     vpContext *vpc));
00053 static void DestroyConstructionBuffer ANSI_ARGS((vpContext *vpc,
00054     ConstructionBuffer *cbuf));
00055 static GBuffer *CreateGBuffer ANSI_ARGS((vpContext *vpc));
00056 static void DestroyGBuffer ANSI_ARGS((vpContext *vpc, GBuffer *gbuf));
00057 #ifdef INDEX_VOLUME
00058 static vpResult ComputeIndex ANSI_ARGS((vpContext *vpc,
00059     RLEVoxels *rle_voxels));
00060 #endif
00061 static vpResult ComputeScanOffsets ANSI_ARGS((vpContext *vpc,
00062     RLEVoxels *rle_voxels));
00063 #ifdef DEBUG
00064 static void ValidateRLEVoxels ANSI_ARGS((vpContext *vpc, RLEVoxels *rle,
00065     int istride, int jstride, int kstride, int axis));
00066 #endif
00067 
00068 /*
00069  * vpClassifyScalars
00070  *
00071  * Classify an array of scalars and store the result in the classified volume.
00072  */
00073 
00074 vpResult
00075 vpClassifyScalars(vpc, scalar_data, length, scalar_field, grad_field,
00076                   norm_field)
00077 vpContext *vpc;         /* context */
00078 unsigned char *scalar_data; /* 3D array of scalar data */
00079 int length;             /* number of scalars in scalar_data */
00080 int scalar_field;       /* voxel field for scalar, or VP_SKIP_FIELD */
00081 int grad_field;         /* voxel field for gradient, or VP_SKIP_FIELD */
00082 int norm_field;         /* voxel field for normal, or VP_SKIP_FIELD */
00083 {
00084     int xlen, ylen, zlen;       /* volume dimensions */
00085     int y, z;                   /* loop indices */
00086     unsigned char *scalar;      /* pointer to current scalar */
00087     int scalar_ystride;         /* stride to next scalar scanline */
00088     int scalar_zstride;         /* stride to next scalar slice */
00089     char *voxel;                /* pointer to current voxel */
00090     unsigned char *s_py, *s_my, *s_pz, *s_mz; /* ptrs to adjacent scans */
00091     int retcode;                /* return code from vpScanlineNormals */
00092     void *voxel_scan;           /* buffer for storing one scan of raw voxels */
00093 
00094     /* check for errors */
00095     xlen = vpc->xlen;
00096     ylen = vpc->ylen;
00097     zlen = vpc->zlen;
00098     if (xlen == 0 || ylen == 0 || zlen == 0 || vpc->raw_bytes_per_voxel == 0)
00099         return(VPSetError(vpc, VPERROR_BAD_VOLUME));
00100     if (xlen * ylen * zlen != length)
00101         return(VPSetError(vpc, VPERROR_BAD_SIZE));
00102 
00103     /* initialize */
00104     scalar = scalar_data;
00105     scalar_ystride = xlen;
00106     scalar_zstride = xlen*ylen;
00107     Alloc(vpc, voxel_scan, void *, xlen*vpc->raw_bytes_per_voxel,"voxel_scan");
00108 
00109     /* compute volume data */
00110     for (z = 0; z < zlen; z++) {
00111         ReportStatus(vpc, (double)z / (double)zlen);
00112         for (y = 0; y < ylen; y++) {
00113             s_my = (y == 0)      ? NULL : scalar - scalar_ystride;
00114             s_py = (y == ylen-1) ? NULL : scalar + scalar_ystride;
00115             s_mz = (z == 0)      ? NULL : scalar - scalar_zstride;
00116             s_pz = (z == zlen-1) ? NULL : scalar + scalar_zstride;
00117             voxel = voxel_scan;
00118             retcode = vpScanlineNormals(vpc, xlen, scalar, s_my, s_py,
00119                                         s_mz, s_pz, voxel, scalar_field,
00120                                         grad_field, norm_field);
00121             if (retcode != VP_OK) {
00122                 Dealloc(vpc, voxel_scan);
00123                 return(retcode);
00124             }
00125             retcode = vpClassifyScanline(vpc, voxel_scan);
00126             if (retcode != VP_OK) {
00127                 Dealloc(vpc, voxel_scan);
00128                 return(retcode);
00129             }
00130             scalar += scalar_ystride;
00131         }
00132         scalar += scalar_zstride - ylen*scalar_ystride;
00133     }
00134     ReportStatus(vpc, 1.0);
00135     Dealloc(vpc, voxel_scan);
00136     return(VP_OK);
00137 }
00138 
00139 /*
00140  * vpClassifyVolume
00141  *
00142  * Classify the current raw volume and store the result in the
00143  * classified volume.
00144  */
00145 
00146 vpResult
00147 vpClassifyVolume(vpc)
00148 vpContext *vpc;         /* context */
00149 {
00150     int xlen, ylen, zlen;       /* volume dimensions */
00151     int y, z;                   /* loop indices */
00152     char *voxel;                /* pointer to current voxel */
00153     int voxel_ystride;          /* stride to next voxel scanline */
00154     int voxel_zstride;          /* stride to next voxel slice */
00155     int retcode;                /* return code */
00156     MinMaxOctree *octree;       /* octree for fast classification */
00157     MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS]; /* stack for octree */
00158 
00159     /* check for errors */
00160     if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
00161         return(retcode);
00162     if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
00163         return(retcode);
00164 
00165     /* initialize */
00166     vpDestroyClassifiedVolume(vpc);
00167     InitRLE(vpc);
00168     xlen = vpc->xlen;
00169     ylen = vpc->ylen;
00170     zlen = vpc->zlen;
00171     voxel = vpc->raw_voxels;
00172     voxel_ystride = vpc->ystride;
00173     voxel_zstride = vpc->zstride;
00174     octree = vpc->mm_octree;
00175     if (octree != NULL) {
00176         VPComputeSummedAreaTable(vpc);
00177         VPClassifyOctree(vpc);
00178     }
00179 
00180     /* compute volume data */
00181     for (z = 0; z < zlen; z++) {
00182         ReportStatus(vpc, (double)z / (double)zlen);
00183         if (octree != NULL)
00184             VPInitOctreeLevelStack(vpc, level_stack, VP_Z_AXIS, z);
00185         for (y = 0; y < ylen; y++) {
00186             EncodeScanline(vpc, voxel, octree, level_stack);
00187             voxel += voxel_ystride;
00188         }
00189         voxel += voxel_zstride - ylen*voxel_ystride;
00190     }
00191     ReportStatus(vpc, 1.0);
00192 
00193     return(VP_OK);
00194 }
00195 
00196 /*
00197  * vpClassifyScanline
00198  *
00199  * Apply the classification function to a scanline of raw voxels and append
00200  * it to the classified volume.
00201  */
00202 
00203 vpResult
00204 vpClassifyScanline(vpc, voxels)
00205 vpContext *vpc;         /* context */
00206 void *voxels;           /* voxel scanline */
00207 {
00208     int retcode;
00209 
00210     /* initialize if this is the first scanline */
00211     if (vpc->cbuf == NULL) {
00212         if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
00213             return(retcode);
00214         vpDestroyClassifiedVolume(vpc);
00215         InitRLE(vpc);
00216     }
00217 
00218     /* encode scanline */
00219     EncodeScanline(vpc, voxels, NULL, NULL);
00220 
00221     return(VP_OK);
00222 }
00223 
00224 /*
00225  * EncodeScanline
00226  *
00227  * Classify and run-length encode one scanline of voxels.
00228  */
00229 
00230 static void
00231 EncodeScanline(vpc, voxels, octree, level_stack)
00232 vpContext *vpc;         /* context */
00233 void *voxels;           /* voxel scanline */
00234 MinMaxOctree *octree;   /* octree for fast classification */
00235 MMOctreeLevel level_stack[VP_MAX_OCTREE_LEVELS]; /* stack for octree */
00236 {
00237     ConstructionBuffer *cbuf;   /* state preserved between calls */
00238     RunData rundata_x;          /* statistics for current x-axis run */
00239     RunData *rundata_y;         /* statistics for all y-axis runs */
00240     RunData *rundata_z;         /* statistics for all z-axis runs */
00241     int skip_rle_x;             /* if true, do not compute rle_x */
00242     int skip_rle_y;             /* if true, do not compute rle_y */
00243     int skip_rle_z;             /* if true, do not compute rle_z */
00244     int y, z;                   /* y and z coordinates of scanline */
00245     int x;                      /* index of current voxel in scanline */
00246     int xlen, ylen, zlen;       /* volume dimensions */
00247     float opacity;              /* current value of the opacity (0.0-1.0) */
00248     float min_opacity;          /* low opacity threshold */
00249     int raw_bytes_per_voxel;    /* bytes in unclassified voxel */
00250     int run_length;             /* length of last run */
00251     char *rawvoxel;             /* current unclassified voxel */
00252     unsigned char *octree_run_ptr;      /* pointer to current run */
00253     int voxel_count;            /* voxels remaining in current run */
00254     int retcode;
00255 
00256     /* initialize */
00257     cbuf = vpc->cbuf;
00258     xlen = vpc->xlen;
00259     ylen = vpc->ylen;
00260     zlen = vpc->zlen;
00261     bzero(&rundata_x, sizeof(RunData));
00262     rundata_y = &cbuf->rundata_y[cbuf->next_z];
00263     rundata_z = &cbuf->rundata_z[cbuf->next_y * xlen];
00264     skip_rle_x = vpc->skip_rle_x;
00265     skip_rle_y = vpc->skip_rle_y;
00266     skip_rle_z = vpc->skip_rle_z;
00267     min_opacity = vpc->min_opacity;
00268     raw_bytes_per_voxel = vpc->raw_bytes_per_voxel;
00269     rawvoxel = voxels;
00270     y = cbuf->next_y;
00271     z = cbuf->next_z;
00272     if (octree != NULL) {
00273         if (cbuf->octree_scans_left == 0) {
00274             cbuf->octree_scans_left = VPComputeScanRuns(vpc, level_stack,
00275                 cbuf->octree_runs, VP_Z_AXIS, y, xlen);
00276         }
00277         cbuf->octree_scans_left--;
00278         octree_run_ptr = cbuf->octree_runs;
00279     }
00280 
00281     /* loop over voxels in the scanline */
00282     x = 0;
00283     while (x < xlen) {
00284         if (octree == NULL) {
00285             /* no octree available, so process all of the voxels in the scan */
00286             voxel_count = xlen;
00287         } else do {
00288             /* skip over a run of zero voxels */
00289             voxel_count = *octree_run_ptr++;
00290             rundata_y += zlen * voxel_count;
00291             rundata_z += voxel_count;
00292             rawvoxel += raw_bytes_per_voxel * voxel_count;
00293             x += voxel_count;
00294 
00295             /* get length of nonzero voxel run */
00296             voxel_count = *octree_run_ptr++;
00297         } while (voxel_count == 0 && x < xlen);
00298 
00299         /* process the voxels in the nonzero run */
00300         while (voxel_count-- > 0) {
00301             /* compute opacity */
00302             opacity = VPClassifyVoxel(vpc, rawvoxel);
00303 
00304             /* compare opacity to threshold */
00305             if (opacity > min_opacity) {
00306                 /* voxel is non-transparent, so save it */
00307                 CacheVoxel(vpc, opacity, rawvoxel);
00308                 if (!skip_rle_x) {
00309                     rundata_y->p.p1.non_zero_count++;
00310                     CountNonZeroVoxel(rundata_y, y, 0, NULL);
00311                 }
00312                 if (!skip_rle_y) {
00313                     rundata_z->p.p1.non_zero_count++;
00314                     CountNonZeroVoxel(rundata_z, z, 0, NULL);
00315                 }
00316                 rundata_x.p.p1.non_zero_count++;
00317                 CountNonZeroVoxel(&rundata_x, x, 0, vpc);
00318             }
00319 
00320             rundata_y += zlen;
00321             rundata_z++;
00322             rawvoxel += raw_bytes_per_voxel;
00323             x++;
00324         } /* while (voxel_count) */
00325     } /* for x */
00326 
00327     /* finish off the statistics for the scanline */
00328     CountNonZeroVoxel(&rundata_x, x, 1, vpc);
00329 
00330     /* update saved state */
00331     cbuf->non_zero_count += rundata_x.p.p1.non_zero_count;
00332     cbuf->x_run_count += rundata_x.p.p1.run_count;
00333 
00334     /* check if this is the last scanline in the volume */
00335     if (++cbuf->next_y == ylen) {
00336         cbuf->next_y = 0;
00337         cbuf->octree_scans_left = 0;
00338         if (++cbuf->next_z == zlen) {
00339             RepackClassifiedVolume(vpc);
00340             DestroyConstructionBuffer(vpc, vpc->cbuf);
00341             vpc->cbuf = NULL;
00342 #ifdef DEBUG
00343             printf("\r");
00344             if (!skip_rle_x) {
00345                 printf("Checking X scanline offsets....\n");
00346                 VPCheckScanOffsets(vpc->rle_x, vpc->rle_bytes_per_voxel);
00347             }
00348             if (!skip_rle_y) {
00349                 printf("Checking Y scanline offsets....\n");
00350                 VPCheckScanOffsets(vpc->rle_y, vpc->rle_bytes_per_voxel);
00351             }
00352             if (!skip_rle_z) {
00353                 printf("Checking Z scanline offsets....\n");
00354                 VPCheckScanOffsets(vpc->rle_z, vpc->rle_bytes_per_voxel);
00355             }
00356             VPValidateClassifiedVolume(vpc);
00357 #endif
00358         }
00359     }
00360 }
00361 
00362 /*
00363  * InitRLE
00364  *
00365  * Initialize in preparation for creating a new run-length encoded volume.
00366  */
00367 
00368 static void
00369 InitRLE(vpc)
00370 vpContext *vpc;
00371 {
00372     int f;
00373     int rle_bytes_per_voxel, size, offset;
00374     int maxsize = 0;
00375 
00376     /* find out how many bytes of the raw voxel are used for shading */
00377     rle_bytes_per_voxel = 0;
00378     for (f = 0; f < vpc->num_shade_fields; f++) {
00379         size = vpc->field_size[f];
00380         offset = vpc->field_offset[f] + size;
00381         if (offset > rle_bytes_per_voxel)
00382             rle_bytes_per_voxel = offset;
00383         if (size > maxsize)
00384             maxsize = size;
00385     }
00386 
00387     /* add one byte for opacity and then pad to the byte boundary of
00388        the largest field in the voxel; this ensures alignment; the
00389        opacity is always stored in the last byte (so the padding
00390        is in between the shading fields and the opacity field) */
00391     rle_bytes_per_voxel++;
00392     rle_bytes_per_voxel = (rle_bytes_per_voxel + maxsize-1) & ~(maxsize-1);
00393     vpc->rle_bytes_per_voxel = rle_bytes_per_voxel;
00394 
00395     /* initialize construction buffer */
00396     vpc->cbuf = CreateConstructionBuffer(vpc);
00397 }
00398 
00399 /*
00400  * CacheVoxel
00401  *
00402  * Cache one voxel's data in the ConstructionBuffer.
00403  */
00404 
00405 static void
00406 CacheVoxel(vpc, opacity, rawvoxel)
00407 vpContext *vpc;                 /* context */
00408 double opacity;                 /* voxel's opacity */
00409 void *rawvoxel;                 /* raw voxel data */
00410 {
00411     ConstructionBuffer *cbuf;   /* state during construction of volume */
00412     GBuffer *data_buf;          /* storage for cached voxels */
00413     void *data_ptr;             /* pointer to current voxel's storage */
00414     int rle_bytes_per_voxel;    /* bytes per voxel after classification */
00415     int opc_int;                /* quantized opacity */
00416 
00417     /* initialize */
00418     cbuf = vpc->cbuf;
00419     data_buf = cbuf->data_buf_tail;
00420     rle_bytes_per_voxel = vpc->rle_bytes_per_voxel;
00421 
00422     /* allocate more memory if necessary */
00423     if (data_buf->bytes_left < rle_bytes_per_voxel) {
00424         /* allocate more memory */
00425         data_buf->next = CreateGBuffer(vpc);
00426         data_buf = data_buf->next;
00427         cbuf->data_buf_tail = data_buf;
00428     }
00429     data_ptr = data_buf->data_ptr;
00430 
00431     /* copy the voxel fields required for shading */
00432     bcopy(rawvoxel, data_ptr, rle_bytes_per_voxel-1);
00433 
00434     /* quantize and store the opacity */
00435     opc_int = opacity*255.;
00436     if (opc_int > 255)
00437         opc_int = 255;
00438     else if (opc_int < 0)
00439         opc_int = 0;
00440     ByteField(data_ptr, rle_bytes_per_voxel-1) = opc_int;
00441 
00442     data_buf->data_ptr += rle_bytes_per_voxel;
00443     data_buf->bytes_left -= rle_bytes_per_voxel;
00444 }
00445 
00446 /*
00447  * CacheLength
00448  *
00449  * Cache one run length in the ConstructionBuffer.
00450  */
00451 
00452 static void
00453 CacheLength(vpc, length)
00454 vpContext *vpc;
00455 int length;
00456 {
00457     GBuffer *lengths_buf;
00458 
00459     lengths_buf = vpc->cbuf->lengths_buf_tail;
00460     if (lengths_buf->bytes_left == 0) {
00461         /* allocate more memory */
00462         lengths_buf->next = CreateGBuffer(vpc);
00463         lengths_buf = lengths_buf->next;
00464         vpc->cbuf->lengths_buf_tail = lengths_buf;
00465     }
00466     *(lengths_buf->data_ptr)++ = length;
00467     lengths_buf->bytes_left--;
00468 }
00469 
00470 /*
00471  * CountNonZeroVoxel
00472  *
00473  * Update the run count and nonzero voxel count for a voxel scanline.
00474  * This routine adds one non-zero voxel to the scanline.  Index
00475  * indicates the position of the voxel in the scanline.  If that
00476  * position is not immediately adjacent to the last non-zero voxel then
00477  * a run of zero voxels is added as well.
00478  *
00479  * If the vpc argument is non-NULL then the lengths of any completed
00480  * runs are written out to the run length buffer.
00481  */
00482 
00483 static void
00484 CountNonZeroVoxel(rundata, index, end_of_scan, vpc)
00485 RunData *rundata;       /* statistics for the scanline */
00486 int index;              /* index of voxel in scanline */
00487 int end_of_scan;        /* if true then finish the scanline instead of
00488                            adding a voxel */
00489 vpContext *vpc;         /* context in which run lengths should be stored */
00490 {
00491     int run_length;
00492 
00493     if (rundata->next_index != index) {
00494         /* a run of zero voxels intervenes between the current index
00495            and the last nonzero voxel that was processed */
00496 
00497         if (rundata->next_index != 0) {
00498             /* add the last nonzero run to the statistics */
00499             run_length = rundata->run_length;
00500             while (run_length > 255) {
00501                 /* run is too long, so split it */
00502                 run_length -= 255;
00503                 rundata->p.p1.run_count += 2;
00504                 if (vpc != NULL) {
00505                     CacheLength(vpc, 255);
00506                     CacheLength(vpc, 0);
00507                 }
00508             }
00509             rundata->p.p1.run_count++;
00510             if (vpc != NULL)
00511                 CacheLength(vpc, run_length);
00512         }
00513 
00514         /* add the last zero run to the statistics */
00515         run_length = index - rundata->next_index;
00516         while (run_length > 255) {
00517             /* run is too long, so split it */
00518             run_length -= 255;
00519             rundata->p.p1.run_count += 2;
00520             if (vpc != NULL) {
00521                 CacheLength(vpc, 255);
00522                 CacheLength(vpc, 0);
00523             }
00524         }
00525         rundata->p.p1.run_count++;
00526         if (vpc != NULL)
00527             CacheLength(vpc, run_length);
00528 
00529         if (end_of_scan) {
00530             /* add a zero-length nonzero run to finish the scanline */
00531             rundata->p.p1.run_count++;
00532             if (vpc != NULL)
00533                 CacheLength(vpc, 0);
00534         } else {
00535             /* start the new run */
00536             rundata->run_length = 1;
00537             rundata->next_index = index + 1;
00538         }
00539     } else if (!end_of_scan) {
00540         /* add a nonzero voxel to the current run */
00541         if (rundata->next_index == 0) {
00542             rundata->p.p1.run_count++;  /* count initial zero run */
00543             if (vpc != NULL)
00544                 CacheLength(vpc, 0);
00545         }
00546         rundata->run_length++;
00547         rundata->next_index = index + 1;
00548     } else {
00549         /* scanline ends with a nonzero voxel run */
00550         run_length = rundata->run_length;
00551         while (run_length > 255) {
00552             /* run is too long, so split it */
00553             run_length -= 255;
00554             rundata->p.p1.run_count += 2;
00555             if (vpc != NULL) {
00556                 CacheLength(vpc, 255);
00557                 CacheLength(vpc, 0);
00558             }
00559         }
00560         rundata->p.p1.run_count++;
00561         if (vpc != NULL)
00562             CacheLength(vpc, run_length);
00563     }
00564 }
00565 
00566 /*
00567  * RepackClassifiedVolume
00568  *
00569  * Repack the data in the ConstructionBuffer after the last call to
00570  * vpClassifyScanline.  This procedure creates the three run-length
00571  * encoded copies of the classified voxels.
00572  */
00573 
00574 static void
00575 RepackClassifiedVolume(vpc)
00576 vpContext *vpc;
00577 {
00578     int xlen, ylen, zlen;       /* volume dimensions */
00579     int x, y, z;                /* voxel coordinates */
00580     int non_zero_count;         /* number of nonzero voxels in volume */
00581     int rle_bytes_per_voxel;    /* bytes per classified voxel */
00582     int skip_rle_x;             /* if true, compute rle_x */
00583     int skip_rle_y;             /* if true, compute rle_y */
00584     int skip_rle_z;             /* if true, compute rle_z */
00585     char *x_data;               /* voxel data for x viewpoint */
00586     char *y_data;               /* voxel data for y viewpoint */
00587     char *z_data;               /* voxel data for z viewpoint */
00588     unsigned char *x_lengths;   /* run length for x viewpoint */
00589     unsigned char *y_lengths;   /* run length for y viewpoint */
00590     unsigned char *z_lengths;   /* run length for z viewpoint */
00591     int z_data_offset;          /* offset to current data value in z volume */
00592     int z_length_offset;        /* offset to current length value in z volume*/
00593     GBuffer *data_buf;          /* next GBuffer containing voxel data */
00594     char *data;                 /* pointer to next voxel */
00595     int data_bytes_left;        /* bytes of data left in data buffer */
00596     GBuffer *lengths_buf;       /* next GBuffer containing length data */
00597     unsigned char *lengths;     /* pointer to next length */
00598     int lengths_bytes_left;     /* bytes of data left in lengths buffer */
00599     int x_run_length;           /* length of current x-scanline run */
00600     int is_non_zero;            /* true if current x-scanline run is nonzero */
00601     RunData *rundata_y;         /* statistics for y-axis runs */
00602     RunData *rundata_z;         /* statistics for z-axis runs */
00603 
00604     /* initialize */
00605     xlen = vpc->xlen;
00606     ylen = vpc->ylen;
00607     zlen = vpc->zlen;
00608     non_zero_count = vpc->cbuf->non_zero_count;
00609     rle_bytes_per_voxel = vpc->rle_bytes_per_voxel;
00610     skip_rle_x = vpc->skip_rle_x;
00611     skip_rle_y = vpc->skip_rle_y;
00612     skip_rle_z = vpc->skip_rle_z;
00613 
00614     /* check for empty volume */
00615     if (non_zero_count == 0) {
00616         if (!skip_rle_x)
00617             vpc->rle_x = CreateEmptyRLEVoxels(vpc, ylen, zlen, xlen);
00618         if (!skip_rle_y)
00619             vpc->rle_y = CreateEmptyRLEVoxels(vpc, zlen, xlen, ylen);
00620         if (!skip_rle_z)
00621             vpc->rle_z = CreateEmptyRLEVoxels(vpc, xlen, ylen, zlen);
00622         return;
00623     }
00624 
00625     /* allocate space for y-axis runs (used for the x viewing axis) */
00626     if (!skip_rle_x) {
00627         vpc->rle_x = CreateRLEVoxelsFromRunData(vpc, ylen, zlen, xlen,
00628             non_zero_count, vpc->cbuf->rundata_y, rle_bytes_per_voxel);
00629         x_data = vpc->rle_x->data;
00630         x_lengths = vpc->rle_x->run_lengths;
00631     }
00632 
00633     /* allocate space for z-axis runs (used for the y viewing axis) */
00634     if (!skip_rle_y) {
00635         vpc->rle_y = CreateRLEVoxelsFromRunData(vpc, zlen, xlen, ylen, 
00636             non_zero_count, vpc->cbuf->rundata_z, rle_bytes_per_voxel);
00637         y_data = vpc->rle_y->data;
00638         y_lengths = vpc->rle_y->run_lengths;
00639     }
00640 
00641     /* allocate space for x-axis runs (used for the z viewing axis) */
00642     if (!skip_rle_z) {
00643         vpc->rle_z = VPCreateRLEVoxels(vpc, xlen, ylen, zlen, non_zero_count,
00644             vpc->cbuf->x_run_count, rle_bytes_per_voxel);
00645         Alloc(vpc, vpc->rle_z->scan_offsets, ScanOffset *,
00646               zlen*sizeof(ScanOffset), "scan_offsets");
00647         vpc->rle_z->scan_offsets_per_slice = 1;
00648         z_data = vpc->rle_z->data;
00649         z_lengths = vpc->rle_z->run_lengths;
00650         z_data_offset = 0;
00651         z_length_offset = 0;
00652     }
00653 
00654     /* copy data into the three RLEVoxels structures */
00655     data_buf = vpc->cbuf->data_buf_head;
00656     data = NULL;
00657     data_bytes_left = 0;
00658     lengths_buf = vpc->cbuf->lengths_buf_head;
00659     lengths = NULL;
00660     lengths_bytes_left = 0;
00661     x_run_length = 0;
00662     is_non_zero = 1;
00663     for (z = 0; z < zlen; z++) {
00664         ReportStatus(vpc, (double)z / (double)zlen);
00665         if (!skip_rle_z) {
00666             vpc->rle_z->scan_offsets[z].first_data = z_data_offset;
00667             vpc->rle_z->scan_offsets[z].first_len =
00668                 (z_length_offset & 0x1) ? z_length_offset + 1 :
00669                 z_length_offset;
00670         }
00671         rundata_z = vpc->cbuf->rundata_z;
00672         for (y = 0; y < ylen; y++) {
00673             rundata_y = &vpc->cbuf->rundata_y[z];
00674             for (x = 0; x < xlen; x++) {
00675                 while (x_run_length == 0) {
00676                     /* find length of next run */
00677                     if (lengths_bytes_left <= 0) {
00678                         /* go to next lengths buffer */
00679                         lengths = (unsigned char *)lengths_buf->data;
00680                         lengths_bytes_left = GBUFFER_SIZE -
00681                                              lengths_buf->bytes_left;
00682                         lengths_buf = lengths_buf->next;
00683                         if (!skip_rle_z) {
00684                             bcopy(lengths, z_lengths, lengths_bytes_left);
00685                             z_lengths += lengths_bytes_left;
00686                         }
00687                     }
00688                     x_run_length = *lengths++;
00689                     lengths_bytes_left--;
00690                     is_non_zero = !is_non_zero;
00691                     z_length_offset++;
00692                 }
00693                 x_run_length--; /* consume one voxel */
00694                 if (is_non_zero) {
00695                     /* find the data for this voxel */
00696                     if (data_bytes_left <= 0) {
00697                         data = data_buf->data;
00698                         data_bytes_left = GBUFFER_SIZE - data_buf->bytes_left;
00699                         data_buf = data_buf->next;
00700                         if (!skip_rle_z) {
00701                             bcopy(data, z_data, data_bytes_left);
00702                             z_data = (char *)z_data + data_bytes_left;
00703                         }
00704                     }
00705 
00706                     /* store voxel */
00707                     if (!skip_rle_x) {
00708                         StoreNonZeroVoxel(data, rle_bytes_per_voxel, x_data,
00709                                           x_lengths, rundata_y, y);
00710                     }
00711                     if (!skip_rle_y) {
00712                         StoreNonZeroVoxel(data, rle_bytes_per_voxel, y_data,
00713                                           y_lengths, rundata_z, z);
00714                     }
00715                     data += rle_bytes_per_voxel;
00716                     data_bytes_left -= rle_bytes_per_voxel;
00717                     z_data_offset += rle_bytes_per_voxel;
00718                 }
00719                 rundata_y += zlen;
00720                 rundata_z++;
00721             } /* for x */
00722         } /* for y */
00723     } /* for z */
00724     ReportStatus(vpc, 1.0);
00725 
00726     if (!skip_rle_x)
00727         PadScanlines(ylen, zlen, xlen, vpc->cbuf->rundata_y, x_lengths);
00728     if (!skip_rle_y)
00729         PadScanlines(zlen, xlen, ylen, vpc->cbuf->rundata_z, y_lengths);
00730 }
00731 
00732 /*
00733  * CreateEmptyRLEVoxels
00734  *
00735  * Create an empty RLEVoxels object (all voxels transparent).
00736  */
00737 
00738 static RLEVoxels *
00739 CreateEmptyRLEVoxels(vpc, ilen, jlen, klen)
00740 vpContext *vpc;
00741 int ilen, jlen, klen;
00742 {
00743     RLEVoxels *rle_voxels;
00744     int j, k;
00745     unsigned char *run_lengths;
00746     ScanOffset *scan_offsets;
00747 
00748     rle_voxels = VPCreateRLEVoxels(vpc, ilen, jlen, klen, 1, 2*jlen*klen, 1);
00749     Alloc(vpc, rle_voxels->scan_offsets, ScanOffset *, klen*sizeof(ScanOffset),
00750           "scan_offsets");
00751     rle_voxels->scan_offsets_per_slice = 1;
00752     run_lengths = rle_voxels->run_lengths;
00753     scan_offsets = rle_voxels->scan_offsets;
00754     for (k = 0; k < klen; k++) {
00755         scan_offsets->first_len = k*jlen*2;
00756         scan_offsets->first_data = 0;
00757         scan_offsets++;
00758         for (j = 0; j < jlen; j++) {
00759             *run_lengths++ = ilen;
00760             *run_lengths++ = 0;
00761         }
00762     }
00763     return(rle_voxels);
00764 }
00765 
00766 /*
00767  * CreateRLEVoxelsFromRunData
00768  *
00769  * Allocate an RLEVoxels structure using the data in a RunData array
00770  * in order to determine the required size.  Also reinitialize the RunData
00771  * array with pointers to the RLEVoxels data for the beginning of
00772  * each scanline.
00773  */
00774 
00775 static RLEVoxels *
00776 CreateRLEVoxelsFromRunData(vpc, ilen, jlen, klen, non_zero_count, run_data,
00777                            rle_bytes_per_voxel)
00778 vpContext *vpc;         /* context */
00779 int ilen, jlen, klen;   /* size of volume in rotated object space */
00780 int non_zero_count;     /* number of nonzero voxels in volume */
00781 RunData *run_data;      /* array of run statistics (jlen*klen entries) */
00782 int rle_bytes_per_voxel;/* number of bytes to allocate for each voxel */
00783 {
00784     int j, k;                   /* scanline and slice number */
00785     int scan_run_count;         /* runs in current scanline */
00786     int run_count;              /* runs in entire volume */
00787     int scan_non_zero_count;    /* nonzero voxels in scanline */
00788     int data_offset;            /* scanline's offset in RLEVoxels->data */
00789     int length_offset;          /* scanline's offset in
00790                                    RLEVoxels->run_lengths */
00791     ScanOffset *slice_offset;   /* offsets for each slice */
00792     RLEVoxels *rle_voxels;      /* return value */
00793 
00794     Alloc(vpc, slice_offset, ScanOffset *, klen*sizeof(ScanOffset),
00795           "scan_offsets");
00796 
00797     /* accumulate the statistics for the last run in each scanline,
00798        count the total number of runs, and store the data and length
00799        offsets for the beginning of the scanline */
00800     data_offset = 0;
00801     length_offset = 0;
00802     run_count = 0;
00803     for (k = 0; k < klen; k++) {
00804         slice_offset[k].first_data = data_offset;
00805         slice_offset[k].first_len = length_offset;
00806         for (j = 0; j < jlen; j++) {
00807             CountNonZeroVoxel(run_data, ilen, 1, NULL);
00808             scan_non_zero_count = run_data->p.p1.non_zero_count;
00809             scan_run_count = run_data->p.p1.run_count;
00810             run_data->run_length = 0;
00811             run_data->next_index = 0;
00812             run_data->p.p2.data_offset = data_offset;
00813             run_data->p.p2.length_offset = length_offset;
00814             data_offset += scan_non_zero_count * rle_bytes_per_voxel;
00815             length_offset += scan_run_count * sizeof(unsigned char);
00816             run_count += scan_run_count;
00817             run_data++;
00818         }
00819     }
00820 
00821     /* allocate space */
00822     rle_voxels = VPCreateRLEVoxels(vpc, ilen, jlen, klen, non_zero_count,
00823                                    run_count, rle_bytes_per_voxel);
00824     rle_voxels->scan_offsets_per_slice = 1;
00825     rle_voxels->scan_offsets = slice_offset;
00826     return(rle_voxels);
00827 }
00828 
00829 /*
00830  * StoreNonZeroVoxel
00831  *
00832  * Store a nonzero voxel in an RLEVoxels object.  This function is
00833  * just like CountNonZeroVoxel except that it actually stores voxel data.
00834  */
00835 
00836 static void
00837 StoreNonZeroVoxel(voxel, rle_bytes_per_voxel, data, lengths, rundata, index)
00838 void *voxel;            /* input voxel data */
00839 int rle_bytes_per_voxel;/* size of voxel */
00840 void *data;             /* location to store voxel */
00841 unsigned char *lengths; /* location to store run lengths */
00842 RunData *rundata;       /* run length statistics for current voxel scanline */
00843 int index;              /* index of voxel in scanline */
00844 {
00845     int run_length;
00846 
00847     /* store the voxel */
00848     if (voxel != NULL) {
00849         bcopy(voxel, (char *)data + rundata->p.p2.data_offset,
00850               rle_bytes_per_voxel);
00851         rundata->p.p2.data_offset += rle_bytes_per_voxel;
00852     }
00853 
00854     /* update run lengths */
00855     if (rundata->next_index != index) {
00856         /* a run of zero voxels intervenes between the current index
00857            and the last nonzero voxel that was processed */
00858 
00859         if (rundata->next_index != 0) {
00860             /* add the last nonzero run to the statistics */
00861             run_length = rundata->run_length;
00862             while (run_length > 255) {
00863                 /* run is too long, so split it */
00864                 run_length -= 255;
00865                 lengths[rundata->p.p2.length_offset++] = 255;
00866                 lengths[rundata->p.p2.length_offset++] = 0;
00867             }
00868             lengths[rundata->p.p2.length_offset++] = run_length;
00869         }
00870 
00871         /* add the last zero run to the statistics */
00872         run_length = index - rundata->next_index;
00873         while (run_length > 255) {
00874             /* run is too long, so split it */
00875             run_length -= 255;
00876             lengths[rundata->p.p2.length_offset++] = 255;
00877             lengths[rundata->p.p2.length_offset++] = 0;
00878         }
00879         lengths[rundata->p.p2.length_offset++] = run_length;
00880 
00881         if (voxel == NULL) {
00882             /* add a zero-length nonzero run to finish the scanline */
00883             lengths[rundata->p.p2.length_offset++] = 0;
00884         } else {
00885             /* start the new run */
00886             rundata->run_length = 1;
00887             rundata->next_index = index + 1;
00888         }
00889     } else if (voxel != NULL) {
00890         /* add a nonzero voxel to the current run */
00891         if (rundata->next_index == 0) {
00892             lengths[rundata->p.p2.length_offset++] = 0;
00893         }
00894         rundata->run_length++;
00895         rundata->next_index = index + 1;
00896     } else {
00897         /* scanline ends with a nonzero voxel run */
00898         run_length = rundata->run_length;
00899         while (run_length > 255) {
00900             /* run is too long, so split it */
00901             run_length -= 255;
00902             lengths[rundata->p.p2.length_offset++] = 255;
00903             lengths[rundata->p.p2.length_offset++] = 0;
00904         }
00905         lengths[rundata->p.p2.length_offset++] = run_length;
00906     }
00907 }
00908 
00909 /*
00910  * PadScanlines
00911  *
00912  * Make sure each scanline has an even number of runs.
00913  */
00914 
00915 static void
00916 PadScanlines(ilen, jlen, klen, run_data, lengths)
00917 int ilen, jlen, klen;   /* size of volume in rotated object space */
00918 RunData *run_data;      /* array of run statistics (jlen*klen entries) */
00919 unsigned char *lengths; /* beginning of run lengths array */
00920 {
00921     int scan_count;             /* number of scanlines */
00922     int scan_run_count;         /* number of runs in scanline */
00923     int scan;                   /* current scanline number */
00924 
00925     scan_count = jlen * klen;
00926     for (scan = 0; scan < scan_count; scan++) {
00927         StoreNonZeroVoxel(NULL, 0, NULL, lengths, run_data, ilen);
00928         run_data++;
00929     }
00930 }
00931 
00932 /*
00933  * VPCreateRLEVoxels
00934  *
00935  *
00936  * Allocate a new RLEVoxels object.
00937  */
00938 
00939 RLEVoxels *
00940 VPCreateRLEVoxels(vpc, ilen, jlen, klen, data_count, run_count,
00941                   rle_bytes_per_voxel)
00942 vpContext *vpc;         /* context */
00943 int ilen, jlen, klen;   /* dimensions in rotated object space */
00944 int data_count;         /* number of nonzero voxels */
00945 int run_count;          /* number of runs */
00946 int rle_bytes_per_voxel;/* bytes of storage for one voxel */
00947 {
00948     RLEVoxels *rle_voxels;
00949 
00950     Alloc(vpc, rle_voxels, RLEVoxels *, sizeof(RLEVoxels), "RLEVoxels");
00951     rle_voxels->ilen = ilen;
00952     rle_voxels->jlen = jlen;
00953     rle_voxels->klen = klen;
00954     rle_voxels->run_count = run_count;
00955     if (run_count > 0) {
00956         Alloc(vpc, rle_voxels->run_lengths, unsigned char *, run_count,
00957               "run_lengths");
00958     } else {
00959         rle_voxels->run_lengths = NULL;
00960     }
00961     rle_voxels->data_count = data_count;
00962     if (data_count > 0) {
00963         Alloc(vpc, rle_voxels->data, void *, data_count * rle_bytes_per_voxel,
00964               "voxeldata");
00965     } else {
00966         rle_voxels->data = NULL;
00967     }
00968     rle_voxels->scan_offsets_per_slice = 0;
00969     rle_voxels->scan_offsets = NULL;
00970     rle_voxels->mmapped = 0;
00971 #ifdef INDEX_VOLUME
00972     rle_voxels->voxel_index = NULL;
00973 #endif
00974     return(rle_voxels);
00975 }
00976 
00977 /*
00978  * VPDestroyRLEVoxels
00979  *
00980  * Destroy an RLEVoxels object.
00981  */
00982 
00983 void
00984 VPDestroyRLEVoxels(vpc, rle_voxels)
00985 vpContext *vpc;
00986 RLEVoxels *rle_voxels;
00987 {
00988     if (!rle_voxels->mmapped) {
00989         if (rle_voxels->run_lengths != NULL)
00990             Dealloc(vpc, rle_voxels->run_lengths);
00991         if (rle_voxels->data != NULL)
00992             Dealloc(vpc, rle_voxels->data);
00993         if (rle_voxels->scan_offsets != NULL)
00994             Dealloc(vpc, rle_voxels->scan_offsets);
00995     }
00996 #ifdef INDEX_VOLUME
00997     if (rle_voxels->voxel_index != NULL)
00998         Dealloc(vpc, rle_voxels->voxel_index);
00999 #endif
01000     Dealloc(vpc, rle_voxels);
01001 }
01002 
01003 /*
01004  * CreateConstructionBuffer
01005  *
01006  * Create a ConstructionBuffer object.
01007  */
01008 
01009 static ConstructionBuffer *
01010 CreateConstructionBuffer(vpc)
01011 vpContext *vpc;
01012 {
01013     ConstructionBuffer *cbuf;
01014     int xlen, ylen, zlen;
01015 
01016     xlen = vpc->xlen;
01017     ylen = vpc->ylen;
01018     zlen = vpc->zlen;
01019     Alloc(vpc, cbuf, ConstructionBuffer *, sizeof(ConstructionBuffer),
01020           "ConstructionBuffer");
01021     Alloc(vpc, cbuf->rundata_y, RunData *, xlen*zlen*sizeof(RunData),
01022           "rundata_y");
01023     Alloc(vpc, cbuf->rundata_z, RunData *, ylen*xlen*sizeof(RunData),
01024           "rundata_z");
01025     bzero(cbuf->rundata_y, xlen*zlen*sizeof(RunData));
01026     bzero(cbuf->rundata_z, ylen*xlen*sizeof(RunData));
01027     cbuf->data_buf_head = CreateGBuffer(vpc);
01028     cbuf->data_buf_tail = cbuf->data_buf_head;
01029     cbuf->lengths_buf_head = CreateGBuffer(vpc);
01030     cbuf->lengths_buf_tail = cbuf->lengths_buf_head;
01031     cbuf->non_zero_count = 0;
01032     cbuf->x_run_count = 0;
01033     cbuf->octree_scans_left = 0;
01034     cbuf->next_z = 0;
01035     cbuf->next_y = 0;
01036     return(cbuf);
01037 }
01038 
01039 /*
01040  * DestroyConstructionBuffer
01041  *
01042  * Destroy a ConstructionBuffer object.
01043  */
01044 
01045 static void
01046 DestroyConstructionBuffer(vpc, cbuf)
01047 vpContext *vpc;
01048 ConstructionBuffer *cbuf;
01049 {
01050     GBuffer *gbuf, *next_gbuf;
01051 
01052     Dealloc(vpc, cbuf->rundata_y);
01053     Dealloc(vpc, cbuf->rundata_z);
01054     for (gbuf = cbuf->data_buf_head; gbuf != NULL; gbuf = next_gbuf) {
01055         next_gbuf = gbuf->next;
01056         DestroyGBuffer(vpc, gbuf);
01057     }
01058     for (gbuf = cbuf->lengths_buf_head; gbuf != NULL; gbuf = next_gbuf) {
01059         next_gbuf = gbuf->next;
01060         DestroyGBuffer(vpc, gbuf);
01061     }
01062     Dealloc(vpc, cbuf);
01063 }
01064 
01065 /*
01066  * CreateGBuffer
01067  *
01068  * Create a GBuffer object.
01069  */
01070 
01071 static GBuffer *
01072 CreateGBuffer(vpc)
01073 vpContext *vpc;
01074 {
01075     GBuffer *gbuf;
01076 
01077     Alloc(vpc, gbuf, GBuffer *, sizeof(GBuffer), "GBuffer");
01078     gbuf->bytes_left = GBUFFER_SIZE;
01079     gbuf->data_ptr = gbuf->data;
01080     gbuf->next = NULL;
01081     return(gbuf);
01082 }
01083 
01084 /*
01085  * DestroyGBuffer
01086  *
01087  * Destroy a GBuffer.
01088  */
01089 
01090 static void
01091 DestroyGBuffer(vpc, gbuf)
01092 vpContext *vpc;
01093 GBuffer *gbuf;
01094 {
01095     Dealloc(vpc, gbuf);
01096 }
01097 
01098 /*
01099  * vpDestroyClassifiedVolume
01100  *
01101  * Free all memory associated with a classified volume.
01102  */
01103 
01104 vpResult
01105 vpDestroyClassifiedVolume(vpc)
01106 vpContext *vpc;
01107 {
01108     if (vpc->cbuf != NULL) {
01109         DestroyConstructionBuffer(vpc, vpc->cbuf);
01110         vpc->cbuf = NULL;
01111     }
01112     if (vpc->rle_x != NULL) {
01113         VPDestroyRLEVoxels(vpc, vpc->rle_x);
01114         vpc->rle_x = NULL;
01115     }
01116     if (vpc->rle_y != NULL) {
01117         VPDestroyRLEVoxels(vpc, vpc->rle_y);
01118         vpc->rle_y = NULL;
01119     }
01120     if (vpc->rle_z != NULL) {
01121         VPDestroyRLEVoxels(vpc, vpc->rle_z);
01122         vpc->rle_z = NULL;
01123     }
01124     return(VP_OK);
01125 }
01126 
01127 #ifdef INDEX_VOLUME
01128 /*
01129  * vpComputeRLEIndex
01130  *
01131  * Compute indexes for the classified volume data in a context.
01132  */
01133 
01134 vpResult
01135 vpComputeRLEIndex(vpc)
01136 vpContext *vpc;
01137 {
01138     vpResult result;
01139 
01140     if ((result = VPComputeRLEScanOffsets(vpc)) != VP_OK)
01141         return(result);
01142     if (vpc->rle_x != NULL) {
01143         if ((result = ComputeIndex(vpc, vpc->rle_x)) != VP_OK)
01144             return(result);
01145     }
01146     if (vpc->rle_y != NULL) {
01147         if ((result = ComputeIndex(vpc, vpc->rle_y)) != VP_OK)
01148             return(result);
01149     }
01150     if (vpc->rle_z != NULL) {
01151         if ((result = ComputeIndex(vpc, vpc->rle_z)) != VP_OK)
01152             return(result);
01153     }
01154     return(VP_OK);
01155 }
01156 
01157 /*
01158  * ComputeIndex
01159  *
01160  * Compute an index that maps 3D voxel coordinates to the RLE run data
01161  * for the corresponding voxel.  The values stored in the index are
01162  * byte offsets to the beginning of the run containing the voxel,
01163  * plus a count indicating the position of the voxel in the run.
01164  * Return value is a result code.
01165  */
01166 
01167 static vpResult
01168 ComputeIndex(vpc, rle_voxels)
01169 vpContext *vpc;
01170 RLEVoxels *rle_voxels;
01171 {
01172     int ilen, jlen, klen;       /* size of volume */
01173     unsigned char *RLElen;      /* pointer to current run length */
01174     VoxelLocation *index;       /* pointer to current index entry */
01175     int i, j, k;                /* current voxel coordinates */
01176     unsigned len_offset;        /* offset in bytes from beginning of
01177                                    scanline to current run length */
01178     unsigned data_offset;       /* offset in bytes from beginning of
01179                                    scanline to current voxel data */
01180     int run_is_zero;            /* true if current run is a zero run */
01181     int run_count;              /* voxels left in current run */
01182     int voxel_size;             /* size of a voxel in bytes */
01183 
01184     ilen = rle_voxels->ilen;
01185     jlen = rle_voxels->jlen;
01186     klen = rle_voxels->klen;
01187     RLElen = rle_voxels->run_lengths;
01188     if (rle_voxels->scan_offsets_per_slice != jlen)
01189         return(VPERROR_BAD_VOLUME);
01190     if (rle_voxels->voxel_index == NULL) {
01191         Alloc(vpc, rle_voxels->voxel_index, VoxelLocation *,
01192               ilen * jlen * klen * sizeof(VoxelLocation), "voxel_index");
01193     }
01194     index = rle_voxels->voxel_index;
01195     voxel_size = vpc->rle_bytes_per_voxel;
01196     run_is_zero = 0;
01197     run_count = 0;
01198     for (k = 0; k < klen; k++) {
01199         for (j = 0; j < jlen; j++) {
01200             ASSERT(run_is_zero == 0);
01201             ASSERT(run_count == 0);
01202             len_offset = 0;
01203             data_offset = 0;
01204             for (i = 0; i < ilen; i++) {
01205                 /* record index for current voxel */
01206                 if (len_offset > 256) {
01207                     Dealloc(vpc, rle_voxels->voxel_index);
01208                     rle_voxels->voxel_index = NULL;
01209                     return(VPERROR_LIMIT_EXCEEDED);
01210                 }
01211                 index->run_count = run_count;
01212                 index->len_offset = len_offset;
01213                 if (run_is_zero)
01214                     index->data_offset = data_offset | INDEX_RUN_IS_ZERO;
01215                 else
01216                     index->data_offset = data_offset;
01217                 index++;
01218 
01219                 /* go on to next voxel */
01220                 while (run_count == 0) {
01221                     run_count = *RLElen++;
01222                     run_is_zero = !run_is_zero;
01223                     len_offset++;
01224                 }
01225                 run_count--;
01226                 if (!run_is_zero)
01227                     data_offset += voxel_size;
01228             }
01229             ASSERT(run_count == 0);
01230             if (run_is_zero) {
01231                 run_count = *RLElen++;
01232                 run_is_zero = !run_is_zero;
01233                 len_offset++;
01234             }
01235         }
01236     }
01237     return(VP_OK);
01238 }
01239 #endif /* INDEX_VOLUME */
01240 
01241 /*
01242  * VPComputeRLEScanOffsets
01243  *
01244  * Recompute the scan_offsets arrays for the classified volume data in
01245  * a context.  Return value is a result code.
01246  */
01247 
01248 vpResult
01249 VPComputeRLEScanOffsets(vpc)
01250 vpContext *vpc;
01251 {
01252     vpResult result;
01253 
01254     if (vpc->rle_x != NULL) {
01255         if ((result = ComputeScanOffsets(vpc, vpc->rle_x)) != VP_OK)
01256             return(result);
01257 #ifdef DEBUG
01258         VPCheckScanOffsets(vpc->rle_x, vpc->rle_bytes_per_voxel);
01259 #endif
01260     }
01261 
01262     if (vpc->rle_y != NULL) {
01263         if ((result = ComputeScanOffsets(vpc, vpc->rle_y)) != VP_OK)
01264             return(result);
01265 #ifdef DEBUG
01266         VPCheckScanOffsets(vpc->rle_y, vpc->rle_bytes_per_voxel);
01267 #endif
01268     }
01269 
01270     if (vpc->rle_z != NULL) {
01271         if ((result = ComputeScanOffsets(vpc, vpc->rle_z)) != VP_OK)
01272             return(result);
01273 #ifdef DEBUG
01274         VPCheckScanOffsets(vpc->rle_z, vpc->rle_bytes_per_voxel);
01275 #endif
01276     }
01277     return(VP_OK);
01278 }
01279 
01280 /*
01281  * ComputeScanOffsets
01282  *
01283  * Recompute the scan_offsets array for a classified volume.
01284  * Return value is a result code.
01285  */
01286 
01287 static vpResult
01288 ComputeScanOffsets(vpc, rle_voxels)
01289 vpContext *vpc;
01290 RLEVoxels *rle_voxels;
01291 {
01292     int ilen, jlen, klen;       /* size of volume */
01293     unsigned char *RLElen;      /* pointer to current run length */
01294     ScanOffset *scan_offset;    /* pointer to current scanline offset */
01295     int i, j, k;                /* current voxel coordinates */
01296     unsigned len_offset;        /* offset in bytes from beginning of
01297                                    run lengths to current run length */
01298     unsigned data_offset;       /* offset in bytes from beginning of
01299                                    voxel data to current voxel data */
01300     int voxel_size;             /* size of a voxel in bytes */
01301     int zerocount, nonzerocount;
01302 
01303     if (rle_voxels->mmapped)
01304         return(VPERROR_IO);
01305     ilen = rle_voxels->ilen;
01306     jlen = rle_voxels->jlen;
01307     klen = rle_voxels->klen;
01308     RLElen = rle_voxels->run_lengths;
01309     if (rle_voxels->scan_offsets_per_slice != jlen) {
01310         if (rle_voxels->scan_offsets != NULL)
01311             Dealloc(vpc, rle_voxels->scan_offsets);
01312         Alloc(vpc, rle_voxels->scan_offsets, ScanOffset *,
01313               klen * jlen * sizeof(ScanOffset), "scan_offsets");
01314         rle_voxels->scan_offsets_per_slice = jlen;
01315     }
01316     scan_offset = rle_voxels->scan_offsets;
01317     len_offset = 0;
01318     data_offset = 0;
01319     voxel_size = vpc->rle_bytes_per_voxel;
01320 
01321     for (k = 0; k < klen; k++) {
01322         for (j = 0; j < jlen; j++) {
01323             scan_offset->first_len = len_offset;
01324             scan_offset->first_data = data_offset;
01325             scan_offset++;
01326             for (i = 0; i < ilen; ) {
01327                 zerocount = *RLElen++;   /* get length of run of zeros */
01328                 nonzerocount = *RLElen++;/* get length of run of non-zeros */
01329                 len_offset += 2;
01330                 data_offset += nonzerocount * voxel_size;
01331                 i += zerocount + nonzerocount;
01332             }
01333             ASSERT(i == ilen);
01334         }
01335     }
01336     return(VP_OK);
01337 }
01338 
01339 #ifdef DEBUG
01340 /*
01341  * VPCheckScanOffsets
01342  *
01343  * Check the scan_offsets field of an RLEVolume for internal consistency.
01344  */
01345 
01346 void
01347 VPCheckScanOffsets(rle_voxels, rle_bytes_per_voxel)
01348 RLEVoxels *rle_voxels;
01349 {
01350     int i, j, k;
01351     int ilen, jlen, klen;
01352     int run_length;
01353     int is_non_zero;
01354     unsigned char *run_length_ptr;
01355     int length_offset;
01356     int data_offset;
01357     int scan_offsets_per_slice;
01358     ScanOffset *scan_offset;
01359 
01360     scan_offsets_per_slice = rle_voxels->scan_offsets_per_slice;
01361     if (scan_offsets_per_slice == 0)
01362         return;
01363     ilen = rle_voxels->ilen;
01364     jlen = rle_voxels->jlen;
01365     klen = rle_voxels->klen;
01366     run_length_ptr = rle_voxels->run_lengths;
01367     run_length = 0;
01368     is_non_zero = 1;
01369     length_offset = 0;
01370     data_offset = 0;
01371     for (k = 0; k < klen; k++) {
01372         for (j = 0; j < jlen; j++) {
01373             if (j < scan_offsets_per_slice) {
01374                 scan_offset = &rle_voxels->scan_offsets[
01375                                k*scan_offsets_per_slice + j];
01376                 if (scan_offset->first_len != length_offset) {
01377                     printf("Bad length offset on slice %d, scanline %d: ",k,j);
01378                     printf("%d should be %d\n", scan_offset->first_len,
01379                            length_offset);
01380                 }
01381                 if (scan_offset->first_data != data_offset) {
01382                     printf("Bad data offset on slice %d, scanline %d: ",k,j);
01383                     printf("%d should be %d\n", scan_offset->first_data,
01384                            data_offset);
01385                 }
01386             }
01387             for (i = 0; i < ilen; i++) {
01388                 while (run_length == 0) {
01389                     run_length = *run_length_ptr++;
01390                     is_non_zero = !is_non_zero;
01391                     length_offset++;
01392                 }
01393                 run_length--;
01394                 if (is_non_zero)
01395                     data_offset += rle_bytes_per_voxel;
01396             }
01397             if (run_length != 0) {
01398                 printf("Run did not terminate at end of scanline ");
01399                 printf("on slice %d, scanline %d\n", k, j);
01400             }
01401             if (!is_non_zero) {
01402                 if (*run_length_ptr++ != 0) {
01403                     printf("Missing zero run at end of scanline ");
01404                     printf("on slice %d, scanline %d\n", k, j);
01405                 }
01406                 is_non_zero = !is_non_zero;
01407                 length_offset++;
01408             }
01409         }
01410     }
01411 }
01412 
01413 /*
01414  * VPValidateClassifiedVolume
01415  *
01416  * Compare the classified volume to the unclassified volume.
01417  */
01418 
01419 void
01420 VPValidateClassifiedVolume(vpc)
01421 vpContext *vpc;
01422 {
01423     if (vpc->raw_voxels == NULL)
01424         return;
01425     if (vpc->rle_z != NULL) {
01426         printf("Checking Z view....\n");
01427         ValidateRLEVoxels(vpc, vpc->rle_z, vpc->xstride, vpc->ystride,
01428                           vpc->zstride, VP_Z_AXIS);
01429     }
01430     if (vpc->rle_y != NULL) {
01431         printf("Checking Y view....\n");
01432         ValidateRLEVoxels(vpc, vpc->rle_y, vpc->zstride, vpc->xstride,
01433                           vpc->ystride, VP_Y_AXIS);
01434     }
01435     if (vpc->rle_x != NULL) {
01436         printf("Checking X view....\n");
01437         ValidateRLEVoxels(vpc, vpc->rle_x, vpc->ystride, vpc->zstride,
01438                           vpc->xstride, VP_X_AXIS);
01439     }
01440 }
01441 
01442 static void
01443 ValidateRLEVoxels(vpc, rle, istride, jstride, kstride, axis)
01444 vpContext *vpc;
01445 RLEVoxels *rle;
01446 int istride, jstride, kstride;
01447 int axis;
01448 {
01449     char *rawvoxel;
01450     char *rlevoxel;
01451     unsigned char *lengths;
01452     int i, j, k;
01453     int count;
01454     int is_non_zero;
01455     int num_runs;
01456     float opacity;
01457     int ilen, jlen, klen;
01458     int founderror;
01459     int raw_opc_int;
01460     int rle_opc_int;
01461     int rle_bytes_per_voxel;
01462 
01463     rawvoxel = (char *)vpc->raw_voxels;
01464     rlevoxel = (char *)rle->data;
01465     lengths = rle->run_lengths;
01466     ilen = rle->ilen;
01467     jlen = rle->jlen;
01468     klen = rle->klen;
01469     rle_bytes_per_voxel = vpc->rle_bytes_per_voxel;
01470     founderror = 0;
01471     for (k = 0; k < klen; k++) {
01472         for (j = 0; j < jlen; j++) {
01473             count = 0;
01474             is_non_zero = 1;
01475             num_runs = 0;
01476             for (i = 0; i < ilen; i++) {
01477                 while (count == 0) {
01478                     count = *lengths++;
01479                     is_non_zero = !is_non_zero;
01480                     if (++num_runs > rle->ilen)
01481                         VPBug("runaway scan detected by ValidateRLEVoxels");
01482                 }
01483                 opacity = VPClassifyVoxel(vpc, rawvoxel);
01484                 if (is_non_zero) {
01485                     if (opacity <= vpc->min_opacity &&
01486                         fabs(opacity - vpc->min_opacity) > 0.001) {
01487                         printf("\n");
01488                         printf("**** zero rawvoxel in nonzero rlerun ****\n");
01489                         printf("voxel (i,j,k)=(%d,%d,%d), viewaxis %d\n",
01490                                i, j, k, axis);
01491                         printf("Actual opacity: %17.15f\n", opacity);
01492                         printf("Threshold:      %17.15f\n", vpc->min_opacity);
01493                         founderror = 1;
01494                     }
01495                     raw_opc_int = (int)rint(opacity*255.);
01496                     rle_opc_int = ByteField(rlevoxel, rle_bytes_per_voxel-1);
01497                     if (abs(raw_opc_int - rle_opc_int) > 1) {
01498                         printf("\n");
01499                         printf("**** rawvoxel and rlevoxel disagree ****\n");
01500                         printf("voxel (i,j,k)=(%d,%d,%d), viewaxis %d\n",
01501                                i, j, k, axis);
01502                         printf("Raw opacity: %3d\n", raw_opc_int);
01503                         printf("RLE opacity: %3d\n", rle_opc_int);
01504                         founderror = 1;
01505                     }
01506                     rlevoxel += rle_bytes_per_voxel;
01507                 } else {
01508                     if (opacity > vpc->min_opacity &&
01509                         fabs(opacity - vpc->min_opacity) > 0.001) {
01510                         printf("\n");
01511                         printf("**** nonzero rawvoxel in zero rlerun ****\n");
01512                         printf("voxel (i,j,k)=(%d,%d,%d), viewaxis %d\n",
01513                                i, j, k, axis);
01514                         printf("Actual opacity: %17.15f\n", opacity);
01515                         printf("Threshold:      %17.15f\n", vpc->min_opacity);
01516                         founderror = 1;
01517                     }
01518                 }
01519                 if (founderror) {
01520                     VPDumpClassifier(vpc);
01521                     VPBug("ValidateRLEVoxels found a problem");
01522                 }
01523                 rawvoxel += istride;
01524                 count--;
01525             }
01526             if (count != 0)
01527                 VPBug("Run did not terminate at end of scanline");
01528             if (!is_non_zero) {
01529                 if (*lengths++ != 0)
01530                     VPBug("Missing zero run at end of scanline");
01531                 is_non_zero = !is_non_zero;
01532             }
01533             rawvoxel += jstride - ilen*istride;
01534         }
01535         rawvoxel += kstride - jlen*jstride;
01536     }
01537 }
01538 #endif
01539 
01540 void
01541 VPDumpView(vpc)
01542 vpContext *vpc;
01543 {
01544     int c;
01545 
01546     printf("MODEL:\n");
01547     for (c = 0; c < 4; c++) {
01548         printf("    %12.6f %12.6f %12.6f %12.6f\n",
01549                vpc->transforms[VP_MODEL][c][0],
01550                vpc->transforms[VP_MODEL][c][1],
01551                vpc->transforms[VP_MODEL][c][2],
01552                vpc->transforms[VP_MODEL][c][3]);
01553     }
01554     printf("VIEW:\n");
01555     for (c = 0; c < 4; c++) {
01556         printf("    %12.6f %12.6f %12.6f %12.6f\n",
01557                vpc->transforms[VP_MODEL][c][0],
01558                vpc->transforms[VP_MODEL][c][1],
01559                vpc->transforms[VP_MODEL][c][2],
01560                vpc->transforms[VP_MODEL][c][3]);
01561     }
01562     printf("PROJECT:\n");
01563     for (c = 0; c < 4; c++) {
01564         printf("    %12.6f %12.6f %12.6f %12.6f\n",
01565                vpc->transforms[VP_MODEL][c][0],
01566                vpc->transforms[VP_MODEL][c][1],
01567                vpc->transforms[VP_MODEL][c][2],
01568                vpc->transforms[VP_MODEL][c][3]);
01569     }
01570 }
01571 
01572 void
01573 VPDumpClassifier(vpc)
01574 vpContext *vpc;
01575 {
01576     int c, d;
01577 
01578     for (d = 0; d < vpc->num_clsfy_params; d++) {
01579         printf("CLASSIFIER PARAM %d:\n   ", d);
01580         for (c = 0; c < vpc->field_max[vpc->param_field[d]]; c++) {
01581             printf(" %8.6f", vpc->clsfy_table[d][c]);
01582             if (c % 8 == 7)
01583                 printf("\n   ");
01584         }
01585         printf("\n");
01586     }
01587 }