00001 /*
00002  * vp_extract.c
00003  *
00004  * Routines to extract fields from a volume.
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:22 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 static int ExtractRawVolume ANSI_ARGS((vpContext *vpc, int x0, int y0, int z0,
00034     int x1, int y1, int z1, int field, void *dst, int dst_xstride,
00035     int dst_ystride, int dst_zstride));
00036 static int ClassifyRawVolume ANSI_ARGS((vpContext *vpc, int correct, 
00037     int x0, int y0, int z0, int x1, int y1, int z1, unsigned char *dst,
00038     int dst_xstride, int dst_ystride, int dst_zstride));
00039 static int ShadeRawVolume ANSI_ARGS((vpContext *vpc, int x0, int y0, int z0,
00040     int x1, int y1, int z1, unsigned char *dst, int dst_xstride,
00041     int dst_ystride, int dst_zstride));
00042 static float CorrectOpacity ANSI_ARGS((vpContext *vpc, int quant_opc,
00043     int x, int y, int z));
00044 static void ShadeVoxel ANSI_ARGS((vpContext *vpc, void *voxel, int x,
00045     int y, int z, float *dst));
00046 static int ExtractClassifiedVolume ANSI_ARGS((vpContext *vpc, int axis,
00047     int x0, int y0, int z0, int x1, int y1, int z1, int field, void *dst,
00048     int dst_xstride, int dst_ystride, int dst_zstride));
00049 
00050 /*
00051  * vpExtract
00052  *
00053  * Extract a field from a volume.
00054  */
00055 
00056 vpResult
00057 vpExtract(vpc, volume_type, x0, y0, z0, x1, y1, z1, field, dst, dst_size,
00058           dst_xstride, dst_ystride, dst_zstride)
00059 vpContext *vpc;         /* context */
00060 int volume_type;        /* which volume representation to extract from */
00061 int x0, y0, z0;         /* origin of extracted region */
00062 int x1, y1, z1;         /* opposite corner of extracted region */
00063 int field;              /* field to extract */
00064 void *dst;              /* buffer to store result into */
00065 int dst_size;           /* size of dst in bytes */
00066 int dst_xstride;        /* stride (in bytes) for destination array */
00067 int dst_ystride;
00068 int dst_zstride;
00069 {
00070     int field_size;
00071     int xrange, yrange, zrange;
00072     int retcode;
00073     int axis;
00074 
00075     /* check for errors */
00076     if (x0 < 0 || y0 < 0 || z0 < 0 || x1 >= vpc->xlen || y1 >= vpc->ylen ||
00077         z1 >= vpc->zlen || x0 > x1 || y0 > y1 || z0 > z1)
00078         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00079     if (field == VP_OPACITY_FIELD || field == VP_CORRECTED_OPAC_FIELD)
00080         field_size = 1;
00081     else if (field == VP_COLOR_FIELD)
00082         field_size = vpc->color_channels;
00083     else if (field < 0 || field >= vpc->num_voxel_fields)
00084         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00085     else if (volume_type != VP_RAW_VOLUME && field >= vpc->num_shade_fields)
00086         return(VPSetError(vpc, VPERROR_BAD_VALUE));
00087     else
00088         field_size = vpc->field_size[field];
00089     if (dst == NULL || dst_size != field_size*(x1-x0+1)*(y1-y0+1)*(z1-z0+1))
00090         return(VPSetError(vpc, VPERROR_BAD_SIZE));
00091 
00092     /* choose axis */
00093     switch (volume_type) {
00094     case VP_CLASSIFIED_VOLUME:
00095         xrange = x1 - x0;
00096         yrange = y1 - y0;
00097         zrange = z1 - z0;
00098         if (vpc->rle_z != NULL && zrange < xrange && zrange < yrange)
00099             axis = VP_Z_AXIS;
00100         else if (vpc->rle_x != NULL && xrange < yrange && xrange < zrange)
00101             axis = VP_X_AXIS;
00102         else if (vpc->rle_z != NULL && yrange < zrange && yrange < xrange)
00103             axis = VP_Y_AXIS;
00104         else if (vpc->rle_z != NULL && xrange >= yrange && xrange >= zrange)
00105             axis = VP_Z_AXIS;
00106         else if (vpc->rle_x != NULL && yrange >= zrange)
00107             axis = VP_X_AXIS;
00108         else if (vpc->rle_y != NULL)
00109             axis = VP_Y_AXIS;
00110         else if (vpc->rle_z != NULL)
00111             axis = VP_Z_AXIS;
00112         else if (vpc->rle_x != NULL)
00113             axis = VP_X_AXIS;
00114         else
00115             return(VPSetError(vpc, VPERROR_BAD_VOLUME));
00116         break;
00117     case VP_CLX_VOLUME:
00118         axis = VP_X_AXIS;
00119         break;
00120     case VP_CLY_VOLUME:
00121         axis = VP_Y_AXIS;
00122         break;
00123     case VP_CLZ_VOLUME:
00124         axis = VP_Z_AXIS;
00125         break;
00126     case VP_RAW_VOLUME:
00127         break;
00128     default:
00129         return(VPSetError(vpc, VPERROR_BAD_OPTION));
00130     }
00131 
00132     /* compute result */
00133     if (volume_type == VP_RAW_VOLUME) {
00134         if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
00135             return(retcode);
00136         if (field == VP_OPACITY_FIELD)
00137             return(ClassifyRawVolume(vpc, 0, x0, y0, z0, x1, y1, z1, dst,
00138                                      dst_xstride, dst_ystride, dst_zstride));
00139         else if (field == VP_CORRECTED_OPAC_FIELD)
00140             return(ClassifyRawVolume(vpc, 1, x0, y0, z0, x1, y1, z1, dst,
00141                                      dst_xstride, dst_ystride, dst_zstride));
00142         else if (field == VP_COLOR_FIELD)
00143             return(ShadeRawVolume(vpc, x0, y0, z0, x1, y1, z1, dst,
00144                                   dst_xstride, dst_ystride, dst_zstride));
00145         else
00146             return(ExtractRawVolume(vpc, x0, y0, z0, x1, y1, z1, field, dst,
00147                                     dst_xstride, dst_ystride, dst_zstride));
00148     } else {
00149         if ((retcode = VPCheckClassifiedVolume(vpc, axis)) != VP_OK)
00150             return(retcode);
00151         if (field == VP_COLOR_FIELD) {
00152             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00153         } else {
00154             return(ExtractClassifiedVolume(vpc, axis, x0, y0, z0, x1, y1, z1,
00155                         field, dst, dst_xstride, dst_ystride, dst_zstride));
00156         }
00157     }
00158 }
00159 
00160 /*
00161  * ExtractRawVolume
00162  *
00163  * Extract a field from a raw volume into an array.
00164  */
00165 
00166 static int
00167 ExtractRawVolume(vpc, x0, y0, z0, x1, y1, z1, field, dst,
00168                  dst_xstride, dst_ystride, dst_zstride)
00169 vpContext *vpc;         /* context */
00170 int x0, y0, z0;         /* origin of extracted region */
00171 int x1, y1, z1;         /* opposite corner of extracted region */
00172 int field;              /* field to extract */
00173 void *dst;              /* buffer to store result into */
00174 int dst_xstride;        /* stride (in bytes) for destination array */
00175 int dst_ystride;
00176 int dst_zstride;
00177 {
00178     int x, y, z;
00179     unsigned char *voxel, *dstptr;
00180     int field_size;
00181     int field_offset;
00182     int xstride, ystride, zstride;
00183     int retcode;
00184 
00185     field_size = vpc->field_size[field];
00186     field_offset = vpc->field_offset[field];
00187     xstride = vpc->xstride;
00188     ystride = vpc->ystride;
00189     zstride = vpc->zstride;
00190     voxel = vpc->raw_voxels;
00191     voxel += x0*xstride + y0*ystride + z0*zstride;
00192     dstptr = dst;
00193     for (z = z0; z <= z1; z++) {
00194         for (y = y0; y <= y1; y++) {
00195             for (x = x0; x <= x1; x++) {
00196                 if (field_size == 1)
00197                     ByteField(dstptr, 0) = ByteField(voxel, field_offset);
00198                 else if (field_size == 2)
00199                     ShortField(dstptr, 0) = ShortField(voxel, field_offset);
00200                 else
00201                     IntField(dstptr, 0) = IntField(voxel, field_offset);
00202                 dstptr += dst_xstride;
00203                 voxel += xstride;
00204             }
00205             dstptr += dst_ystride - (x1-x0+1)*dst_xstride;
00206             voxel += ystride - (x1-x0+1)*xstride;
00207         }
00208         dstptr += dst_zstride - (y1-y0+1)*dst_ystride;
00209         voxel += zstride - (y1-y0+1)*ystride;
00210     }
00211     return(VP_OK);
00212 }
00213 
00214 /*
00215  * ClassifyRawVolume
00216  *
00217  * Classify a portion of a raw volume, quantize the result, and store
00218  * as an array of 8-bit opacities.
00219  */
00220 
00221 static int
00222 ClassifyRawVolume(vpc, correct, x0, y0, z0, x1, y1, z1, dst,
00223                   dst_xstride, dst_ystride, dst_zstride)
00224 vpContext *vpc;         /* context */
00225 int correct;            /* if true then correct for view */
00226 int x0, y0, z0;         /* origin of extracted region */
00227 int x1, y1, z1;         /* opposite corner of extracted region */
00228 unsigned char *dst;     /* buffer to store result into */
00229 int dst_xstride;        /* stride (in bytes) for destination array */
00230 int dst_ystride;
00231 int dst_zstride;
00232 {
00233     float *opc;
00234     int num_voxels;
00235     int retcode;
00236 
00237     /* check for errors */
00238     if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
00239         return(retcode);
00240 
00241     /* compute opacities */
00242     num_voxels = (x1-x0+1)*(y1-y0+1)*(z1-z0+1);
00243     Alloc(vpc, opc, float *, num_voxels*sizeof(float), "opacity_block");
00244     VPClassifyBlock(vpc, correct, x0, y0, z0, x1, y1, z1, opc,
00245                     sizeof(float), (x1-x0+1)*sizeof(float),
00246                     (x1-x0+1)*(y1-y0+1)*sizeof(float));
00247 
00248     /* quantize opacities */
00249     VPQuantize(opc, x1-x0+1, y1-y0+1, z1-z0+1, 255., 255, dst,
00250                dst_xstride, dst_ystride, dst_zstride);
00251 
00252     Dealloc(vpc, opc);
00253     return(VP_OK);
00254 }
00255 
00256 /*
00257  * ShadeRawVolume
00258  *
00259  * Shade a portion of a raw volume, quantize the result, and store
00260  * as an array of 8-bit intensities.
00261  */
00262 
00263 static int
00264 ShadeRawVolume(vpc, x0, y0, z0, x1, y1, z1, dst,
00265                dst_xstride, dst_ystride, dst_zstride)
00266 vpContext *vpc;         /* context */
00267 int x0, y0, z0;         /* origin of extracted region */
00268 int x1, y1, z1;         /* opposite corner of extracted region */
00269 unsigned char *dst;     /* buffer to store result into */
00270 int dst_xstride;        /* stride (in bytes) for destination array */
00271 int dst_ystride;
00272 int dst_zstride;
00273 {
00274     float *shd;
00275     int num_colors;
00276     int retcode;
00277     int xstride, ystride, zstride;
00278 
00279     /* check for errors */
00280     if ((retcode = VPCheckShader(vpc)) != VP_OK)
00281         return(retcode);
00282 
00283     /* compute colors */
00284     num_colors = (x1-x0+1)*(y1-y0+1)*(z1-z0+1)*vpc->color_channels;
00285     Alloc(vpc, shd, float *, num_colors*sizeof(float), "color_block");
00286     xstride = vpc->color_channels * sizeof(float);
00287     ystride = xstride * (x1-x0+1);
00288     zstride = ystride * (y1-y0+1);
00289     VPShadeBlock(vpc, x0, y0, z0, x1, y1, z1, shd, xstride, ystride, zstride);
00290 
00291     /* quantize colors */
00292     VPQuantize(shd, x1-x0+1, y1-y0+1, z1-z0+1, 1., 255, dst,
00293                dst_xstride, dst_ystride, dst_zstride);
00294 
00295     Dealloc(vpc, shd);
00296     return(VP_OK);
00297 }
00298 
00299 /*
00300  * VPClassifyBlock
00301  *
00302  * Classify a block of the current raw volume.  The result is an
00303  * array of floating point opacities in the range 0.0-1.0.
00304  */
00305 
00306 vpResult
00307 VPClassifyBlock(vpc, correct, x0, y0, z0, x1, y1, z1, opc,
00308                 dst_xstride, dst_ystride, dst_zstride)
00309 vpContext *vpc;         /* context */
00310 int correct;            /* if true then correct for view */
00311 int x0, y0, z0;         /* origin of extracted region */
00312 int x1, y1, z1;         /* opposite corner of extracted region */
00313 float *opc;             /* buffer to store result into */
00314 int dst_xstride;        /* stride (in bytes) for destination array */
00315 int dst_ystride;
00316 int dst_zstride;
00317 {
00318     unsigned char *voxel;
00319     int xstride, ystride, zstride;
00320     int x, y, z;
00321     float opacity;
00322     int quant_opc;
00323     int retcode;
00324 
00325     if (correct) {
00326         if ((retcode = VPFactorView(vpc)) != VP_OK)
00327             return(retcode);
00328     }
00329     xstride = vpc->xstride;
00330     ystride = vpc->ystride;
00331     zstride = vpc->zstride;
00332     voxel = vpc->raw_voxels;
00333     voxel += x0*xstride + y0*ystride + z0*zstride;
00334     for (z = z0; z <= z1; z++) {
00335         for (y = y0; y <= y1; y++) {
00336             for (x = x0; x <= x1; x++) {
00337                 opacity = VPClassifyVoxel(vpc, voxel);
00338                 if (correct) {
00339                     quant_opc = opacity * 255.;
00340                     if (quant_opc > 255)
00341                         quant_opc = 255;
00342                     else if (quant_opc < 0)
00343                         quant_opc = 0;
00344                     opacity = CorrectOpacity(vpc, quant_opc, x, y, z);
00345                 }
00346                 *opc = opacity;
00347                 opc = (float *)((char *)opc + dst_xstride);
00348                 voxel += xstride;
00349             }
00350             opc = (float *)((char *)opc + dst_ystride - (x1-x0+1)*dst_xstride);
00351             voxel += ystride - (x1-x0+1)*xstride;
00352         }
00353         opc = (float *)((char *)opc + dst_zstride - (y1-y0+1)*dst_ystride);
00354         voxel += zstride - (y1-y0+1)*ystride;
00355     }
00356     return(VP_OK);
00357 }
00358 
00359 /*
00360  * VPClassifyVoxel
00361  *
00362  * Classify a single voxel.  Return value is an opacity.
00363  */
00364 
00365 float
00366 VPClassifyVoxel(vpc, voxel)
00367 vpContext *vpc;         /* context */
00368 void *voxel;            /* pointer to voxel */
00369 {
00370     int num_params;             /* number of parameters to classifier */
00371     int p;                      /* current parameter number */
00372     int field;                  /* field for the parameter */
00373     int field_size;             /* size of the field */
00374     int field_offset;           /* offset for the field */
00375     int index;                  /* index for table lookup */
00376     float opacity;              /* current value of the opacity */
00377 
00378     num_params = vpc->num_clsfy_params;
00379     opacity = 1;
00380     for (p = 0; p < num_params; p++) {
00381         /* get table index */
00382         field = vpc->param_field[p];
00383         field_offset = vpc->field_offset[field];
00384         field_size = vpc->field_size[field];
00385         index = VoxelField(voxel, field_offset, field_size);
00386 
00387         /* load table value */
00388         opacity *= vpc->clsfy_table[p][index];
00389     }
00390     return(opacity);
00391 }
00392 
00393 /*
00394  * CorrectOpacity
00395  *
00396  * Correct an opacity for the current view.
00397  * Return value is the corrected opacity.
00398  */
00399 
00400 static float
00401 CorrectOpacity(vpc, quant_opc, x, y, z)
00402 vpContext *vpc;         /* context */
00403 int quant_opc;          /* input opacity (0-255) */
00404 int x, y, z;            /* voxel coordinates in object space */
00405 {
00406     float opacity;
00407 
00408     if (vpc->affine_view) {
00409         opacity = vpc->affine_opac_correct[quant_opc];
00410     } else {
00411         /* XXX perspective rendering not available yet */
00412         opacity = (float)quant_opc / (float)255.;
00413     }
00414     return(opacity);
00415 }
00416 
00417 /*
00418  * VPShadeBlock
00419  *
00420  * Shade a block of the current raw volume.  The result is an
00421  * array of floating point colors in the range 0.0-255.0.
00422  */
00423 
00424 vpResult
00425 VPShadeBlock(vpc, x0, y0, z0, x1, y1, z1, shd,
00426              dst_xstride, dst_ystride, dst_zstride)
00427 vpContext *vpc;         /* context */
00428 int x0, y0, z0;         /* origin of extracted region */
00429 int x1, y1, z1;         /* opposite corner of extracted region */
00430 float *shd;             /* buffer to store result into */
00431 int dst_xstride;        /* stride (in bytes) for destination array */
00432 int dst_ystride;
00433 int dst_zstride;
00434 {
00435     unsigned char *voxel;
00436     int xstride, ystride, zstride;
00437     int x, y, z;
00438     int color_channels;
00439 
00440     color_channels = vpc->color_channels;
00441     xstride = vpc->xstride;
00442     ystride = vpc->ystride;
00443     zstride = vpc->zstride;
00444     voxel = vpc->raw_voxels;
00445     voxel += x0*xstride + y0*ystride + z0*zstride;
00446     for (z = z0; z <= z1; z++) {
00447         for (y = y0; y <= y1; y++) {
00448             for (x = x0; x <= x1; x++) {
00449                 ShadeVoxel(vpc, voxel, x, y, z, shd);
00450                 shd = (float *)((char *)shd + dst_xstride);
00451                 voxel += xstride;
00452             }
00453             shd = (float *)((char *)shd + dst_ystride - (x1-x0+1)*dst_xstride);
00454             voxel += ystride - (x1-x0+1)*xstride;
00455         }
00456         shd = (float *)((char *)shd + dst_zstride - (y1-y0+1)*dst_ystride);
00457         voxel += zstride - (y1-y0+1)*ystride;
00458     }
00459     return(VP_OK);
00460 }
00461 
00462 /*
00463  * ShadeVoxel
00464  *
00465  * Shade a voxel.
00466  */
00467 
00468 static void
00469 ShadeVoxel(vpc, voxel, x, y, z, dst)
00470 vpContext *vpc;         /* context */
00471 void *voxel;            /* voxel data */
00472 int x, y, z;            /* voxel coordinates */
00473 float *dst;             /* storage for result (1 or 3 intensities, 0-255) */
00474 {
00475     int num_materials;
00476     int color_channels;
00477     int color_index_size, color_index_offset, color_index, color_table_offset;
00478     int weight_index_size, weight_index_offset, weight_index;
00479     int weight_table_offset;
00480     int m;
00481     float r, g, b;
00482     float *color_table;
00483     float *weight_table;
00484 
00485     /* check shading mode */
00486     if (vpc->shading_mode == CALLBACK_SHADER) {
00487         if (vpc->color_channels == 1)
00488             vpc->shade_func(voxel, dst, vpc->client_data);
00489         else
00490             vpc->shade_func(voxel, dst, dst+1, dst+2, vpc->client_data);
00491         return;
00492     } else if (vpc->shading_mode != LOOKUP_SHADER) {
00493         VPBug("unknown shader type");
00494     }
00495 
00496     /* compute table indices */
00497     num_materials = vpc->num_materials;
00498     color_channels = vpc->color_channels;
00499     color_index_size = vpc->field_size[vpc->color_field];
00500     color_index_offset = vpc->field_offset[vpc->color_field];
00501     color_index = VoxelField(voxel, color_index_offset, color_index_size);
00502     color_table_offset = color_index * num_materials;
00503     weight_index_size = vpc->field_size[vpc->weight_field];
00504     weight_index_offset = vpc->field_offset[vpc->weight_field];
00505     weight_index = VoxelField(voxel, weight_index_offset, weight_index_size);
00506     weight_table_offset = weight_index * num_materials;
00507 
00508     /* look up values in tables */
00509     if (color_channels == 1) {
00510         color_table = vpc->shade_color_table + color_table_offset;
00511         weight_table = vpc->shade_weight_table + weight_table_offset;
00512         if (num_materials == 1) {
00513             r = *color_table;
00514         } else {
00515             r = 0;
00516             for (m = 0; m < num_materials; m++)
00517                 r += *color_table++ * *weight_table++;
00518         }
00519         *dst = r;
00520     } else {
00521         color_table = vpc->shade_color_table + 3*color_table_offset;
00522         weight_table = vpc->shade_weight_table + weight_table_offset;
00523         if (num_materials == 1) {
00524             r = *color_table++;
00525             g = *color_table++;
00526             b = *color_table;
00527         } else {
00528             r = 0;
00529             g = 0;
00530             b = 0;
00531             for (m = 0; m < num_materials; m++) {
00532                 r += *color_table++ * *weight_table;
00533                 g += *color_table++ * *weight_table;
00534                 b += *color_table++ * *weight_table;
00535             }
00536         }
00537         dst[0] = r;
00538         dst[1] = g;
00539         dst[2] = b;
00540     }
00541 }
00542 
00543 /*
00544  * VPQuantize
00545  *
00546  * Quantize a floating point array and store the result in a byte array.
00547  */
00548 
00549 void
00550 VPQuantize(src, xlen, ylen, zlen, scale, maxvalue, dst,
00551            dst_xstride, dst_ystride, dst_zstride)
00552 float *src;             /* floating point array */
00553 int xlen, ylen, zlen;   /* array dimensions */
00554 double scale;           /* scale to apply to each array element */
00555 int maxvalue;           /* clamp each array element to this value */
00556 unsigned char *dst;     /* store results here */
00557 int dst_xstride;        /* stride (in bytes) for destination array */
00558 int dst_ystride;
00559 int dst_zstride;
00560 {
00561     int value;
00562     int x, y, z;
00563 
00564     for (z = 0; z < zlen; z++) {
00565         for (y = 0; y < ylen; y++) {
00566             for (x = 0; x < xlen; x++) {
00567                 value = (int)rint(*src++ * scale);
00568                 if (value > maxvalue)
00569                     value = maxvalue;
00570                 else if (value < 0)
00571                     value = 0;
00572                 *dst = value;
00573                 dst += dst_xstride;
00574             }
00575             dst += dst_ystride - xlen*dst_xstride;
00576         }
00577         dst += dst_zstride - ylen*dst_ystride;
00578     }
00579 }
00580 
00581 /*
00582  * ExtractClassifiedVolume
00583  *
00584  * Extract a field from a classified volume into an array.
00585  */
00586 
00587 static int
00588 ExtractClassifiedVolume(vpc, axis, x0, y0, z0, x1, y1, z1, field, dst,
00589                         dst_xstride, dst_ystride, dst_zstride)
00590 vpContext *vpc;         /* context */
00591 int axis;               /* which axis to extract from */
00592 int x0, y0, z0;         /* origin of extracted region */
00593 int x1, y1, z1;         /* opposite corner of extracted region */
00594 int field;              /* field to extract */
00595 void *dst;              /* buffer to store result into */
00596 int dst_xstride;        /* stride (in bytes) for destination array */
00597 int dst_ystride;
00598 int dst_zstride;
00599 {
00600     int i, j, k;                /* voxel coordinates in rotated object space */
00601     int i0, j0, k0;             /* origin of extracted region */
00602     int i1, j1, k1;             /* opposite corner of extracted region */
00603     int dst_istride;            /* stride (in bytes) for destination array */
00604     int dst_jstride;
00605     int dst_kstride;
00606     int ilen, jlen, klen;       /* volume size */
00607     RLEVoxels *rle_voxels;      /* run-length encoded, classified volume */
00608     unsigned char *voxel;       /* pointer to current voxel in volume */
00609     unsigned char *dstptr;      /* pointer to destination */
00610     unsigned char *length;      /* pointer to current run length */
00611     int run_length;             /* length of current run */
00612     int is_non_zero;            /* true if current run is nonzero */
00613     int rle_bytes_per_voxel;    /* size of unclassified voxel */
00614     int value;                  /* value of parameter for current voxel */
00615     ScanOffset *slice_runs;     /* offsets to start of runs for a slice */
00616     int field_size;             /* size of field in bytes */
00617     int field_offset;           /* byte offset for voxel field */
00618 
00619     /* initialize */
00620     switch (axis) {
00621     case VP_X_AXIS:
00622         rle_voxels = vpc->rle_x;
00623         i0 = y0; j0 = z0; k0 = x0; i1 = y1; j1 = z1; k1 = x1;
00624         dst_istride = dst_ystride;
00625         dst_jstride = dst_zstride;
00626         dst_kstride = dst_xstride;
00627         break;
00628     case VP_Y_AXIS:
00629         rle_voxels = vpc->rle_y;
00630         i0 = z0; j0 = x0; k0 = y0; i1 = z1; j1 = x1; k1 = y1;
00631         dst_istride = dst_zstride;
00632         dst_jstride = dst_xstride;
00633         dst_kstride = dst_ystride;
00634         break;
00635     case VP_Z_AXIS:
00636         rle_voxels = vpc->rle_z;
00637         i0 = x0; j0 = y0; k0 = z0; i1 = x1; j1 = y1; k1 = z1;
00638         dst_istride = dst_xstride;
00639         dst_jstride = dst_ystride;
00640         dst_kstride = dst_zstride;
00641         break;
00642     default:
00643         return(VPSetError(vpc, VPERROR_BAD_OPTION));
00644     }
00645     if (rle_voxels == NULL)
00646         return(VPSetError(vpc, VPERROR_BAD_VOLUME));
00647     if (rle_voxels->scan_offsets_per_slice < 1)
00648         return(VPSetError(vpc, VPERROR_BAD_VOLUME));
00649     ilen = rle_voxels->ilen;
00650     jlen = rle_voxels->jlen;
00651     klen = rle_voxels->klen;
00652     rle_bytes_per_voxel = vpc->rle_bytes_per_voxel;
00653     if (field == VP_OPACITY_FIELD || field == VP_CORRECTED_OPAC_FIELD) {
00654         field_size = 1;
00655         field_offset = rle_bytes_per_voxel - 1;
00656     } else {
00657         field_size = vpc->field_size[field];
00658         field_offset = vpc->field_offset[field];
00659     }
00660 
00661     /* extract slice */
00662     dstptr = dst;
00663     for (k = k0; k <= k1; k++) {
00664         slice_runs = &rle_voxels->scan_offsets[k *
00665                         rle_voxels->scan_offsets_per_slice];
00666         voxel = (unsigned char *)rle_voxels->data + slice_runs->first_data;
00667         length = rle_voxels->run_lengths + slice_runs->first_len;
00668         run_length = 0;
00669         is_non_zero = 1;
00670         for (j = 0; j < jlen; j++) {
00671             for (i = 0; i < ilen; i++) {
00672                 while (run_length == 0) {
00673                     run_length = *length++;
00674                     is_non_zero = !is_non_zero;
00675                 }
00676                 run_length--;
00677                 if (i >= i0 && i <= i1 && j >= j0 && j <= j1) {
00678                     if (is_non_zero) {
00679                         if (field_size == 1)
00680                             ByteField(dstptr, 0) = ByteField(voxel,
00681                                                              field_offset);
00682                         else if (field_size == 2)
00683                             ShortField(dstptr, 0) = ShortField(voxel,
00684                                                                field_offset);
00685                         else
00686                             IntField(dstptr, 0) = IntField(voxel,field_offset);
00687                         voxel += rle_bytes_per_voxel;
00688                     } else {
00689                         if (field_size == 1)
00690                             ByteField(dstptr, 0) = 0;
00691                         else if (field_size == 2)
00692                             ShortField(dstptr, 0) = 0;
00693                         else
00694                             IntField(dstptr, 0) = 0;
00695                     }
00696                     dstptr += dst_istride;
00697                 } else {
00698                     if (is_non_zero)
00699                         voxel += rle_bytes_per_voxel;
00700                 }
00701             }
00702             if (j >= j0 && j <= j1)
00703                 dstptr += dst_jstride - (i1-i0+1)*dst_istride;
00704         }
00705         dstptr += dst_kstride - (j1-j0+1)*dst_jstride;
00706     }
00707     return(VP_OK);
00708 }