1 /*
2 * vp_extract.c
3 *
4 * Routines to extract fields from a volume.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.28 $
29 */
30
31 #include "vp_global.h"
32
33 static int ExtractRawVolume ANSI_ARGS((vpContext *vpc, int x0, int y0, int z0,
34 int x1, int y1, int z1, int field, void *dst, int dst_xstride,
35 int dst_ystride, int dst_zstride));
36 static int ClassifyRawVolume ANSI_ARGS((vpContext *vpc, int correct,
37 int x0, int y0, int z0, int x1, int y1, int z1, unsigned char *dst,
38 int dst_xstride, int dst_ystride, int dst_zstride));
39 static int ShadeRawVolume ANSI_ARGS((vpContext *vpc, int x0, int y0, int z0,
40 int x1, int y1, int z1, unsigned char *dst, int dst_xstride,
41 int dst_ystride, int dst_zstride));
42 static float CorrectOpacity ANSI_ARGS((vpContext *vpc, int quant_opc,
43 int x, int y, int z));
44 static void ShadeVoxel ANSI_ARGS((vpContext *vpc, void *voxel, int x,
45 int y, int z, float *dst));
46 static int ExtractClassifiedVolume ANSI_ARGS((vpContext *vpc, int axis,
47 int x0, int y0, int z0, int x1, int y1, int z1, int field, void *dst,
48 int dst_xstride, int dst_ystride, int dst_zstride));
49
50 /*
51 * vpExtract
52 *
53 * Extract a field from a volume.
54 */
55
56 vpResult
vpExtract(vpc,volume_type,x0,y0,z0,x1,y1,z1,field,dst,dst_size,dst_xstride,dst_ystride,dst_zstride)57 vpExtract(vpc, volume_type, x0, y0, z0, x1, y1, z1, field, dst, dst_size,
58 dst_xstride, dst_ystride, dst_zstride)
59 vpContext *vpc; /* context */
60 int volume_type; /* which volume representation to extract from */
61 int x0, y0, z0; /* origin of extracted region */
62 int x1, y1, z1; /* opposite corner of extracted region */
63 int field; /* field to extract */
64 void *dst; /* buffer to store result into */
65 int dst_size; /* size of dst in bytes */
66 int dst_xstride; /* stride (in bytes) for destination array */
67 int dst_ystride;
68 int dst_zstride;
69 {
70 int field_size;
71 int xrange, yrange, zrange;
72 int retcode;
73 int axis;
74
75 /* check for errors */
76 if (x0 < 0 || y0 < 0 || z0 < 0 || x1 >= vpc->xlen || y1 >= vpc->ylen ||
77 z1 >= vpc->zlen || x0 > x1 || y0 > y1 || z0 > z1)
78 return(VPSetError(vpc, VPERROR_BAD_VALUE));
79 if (field == VP_OPACITY_FIELD || field == VP_CORRECTED_OPAC_FIELD)
80 field_size = 1;
81 else if (field == VP_COLOR_FIELD)
82 field_size = vpc->color_channels;
83 else if (field < 0 || field >= vpc->num_voxel_fields)
84 return(VPSetError(vpc, VPERROR_BAD_VALUE));
85 else if (volume_type != VP_RAW_VOLUME && field >= vpc->num_shade_fields)
86 return(VPSetError(vpc, VPERROR_BAD_VALUE));
87 else
88 field_size = vpc->field_size[field];
89 if (dst == NULL || dst_size != field_size*(x1-x0+1)*(y1-y0+1)*(z1-z0+1))
90 return(VPSetError(vpc, VPERROR_BAD_SIZE));
91
92 /* choose axis */
93 switch (volume_type) {
94 case VP_CLASSIFIED_VOLUME:
95 xrange = x1 - x0;
96 yrange = y1 - y0;
97 zrange = z1 - z0;
98 if (vpc->rle_z != NULL && zrange < xrange && zrange < yrange)
99 axis = VP_Z_AXIS;
100 else if (vpc->rle_x != NULL && xrange < yrange && xrange < zrange)
101 axis = VP_X_AXIS;
102 else if (vpc->rle_z != NULL && yrange < zrange && yrange < xrange)
103 axis = VP_Y_AXIS;
104 else if (vpc->rle_z != NULL && xrange >= yrange && xrange >= zrange)
105 axis = VP_Z_AXIS;
106 else if (vpc->rle_x != NULL && yrange >= zrange)
107 axis = VP_X_AXIS;
108 else if (vpc->rle_y != NULL)
109 axis = VP_Y_AXIS;
110 else if (vpc->rle_z != NULL)
111 axis = VP_Z_AXIS;
112 else if (vpc->rle_x != NULL)
113 axis = VP_X_AXIS;
114 else
115 return(VPSetError(vpc, VPERROR_BAD_VOLUME));
116 break;
117 case VP_CLX_VOLUME:
118 axis = VP_X_AXIS;
119 break;
120 case VP_CLY_VOLUME:
121 axis = VP_Y_AXIS;
122 break;
123 case VP_CLZ_VOLUME:
124 axis = VP_Z_AXIS;
125 break;
126 case VP_RAW_VOLUME:
127 break;
128 default:
129 return(VPSetError(vpc, VPERROR_BAD_OPTION));
130 }
131
132 /* compute result */
133 if (volume_type == VP_RAW_VOLUME) {
134 if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
135 return(retcode);
136 if (field == VP_OPACITY_FIELD)
137 return(ClassifyRawVolume(vpc, 0, x0, y0, z0, x1, y1, z1, dst,
138 dst_xstride, dst_ystride, dst_zstride));
139 else if (field == VP_CORRECTED_OPAC_FIELD)
140 return(ClassifyRawVolume(vpc, 1, x0, y0, z0, x1, y1, z1, dst,
141 dst_xstride, dst_ystride, dst_zstride));
142 else if (field == VP_COLOR_FIELD)
143 return(ShadeRawVolume(vpc, x0, y0, z0, x1, y1, z1, dst,
144 dst_xstride, dst_ystride, dst_zstride));
145 else
146 return(ExtractRawVolume(vpc, x0, y0, z0, x1, y1, z1, field, dst,
147 dst_xstride, dst_ystride, dst_zstride));
148 } else {
149 if ((retcode = VPCheckClassifiedVolume(vpc, axis)) != VP_OK)
150 return(retcode);
151 if (field == VP_COLOR_FIELD) {
152 return(VPSetError(vpc, VPERROR_BAD_VALUE));
153 } else {
154 return(ExtractClassifiedVolume(vpc, axis, x0, y0, z0, x1, y1, z1,
155 field, dst, dst_xstride, dst_ystride, dst_zstride));
156 }
157 }
158 }
159
160 /*
161 * ExtractRawVolume
162 *
163 * Extract a field from a raw volume into an array.
164 */
165
166 static int
ExtractRawVolume(vpc,x0,y0,z0,x1,y1,z1,field,dst,dst_xstride,dst_ystride,dst_zstride)167 ExtractRawVolume(vpc, x0, y0, z0, x1, y1, z1, field, dst,
168 dst_xstride, dst_ystride, dst_zstride)
169 vpContext *vpc; /* context */
170 int x0, y0, z0; /* origin of extracted region */
171 int x1, y1, z1; /* opposite corner of extracted region */
172 int field; /* field to extract */
173 void *dst; /* buffer to store result into */
174 int dst_xstride; /* stride (in bytes) for destination array */
175 int dst_ystride;
176 int dst_zstride;
177 {
178 int x, y, z;
179 unsigned char *voxel, *dstptr;
180 int field_size;
181 int field_offset;
182 int xstride, ystride, zstride;
183 int retcode;
184
185 field_size = vpc->field_size[field];
186 field_offset = vpc->field_offset[field];
187 xstride = vpc->xstride;
188 ystride = vpc->ystride;
189 zstride = vpc->zstride;
190 voxel = vpc->raw_voxels;
191 voxel += x0*xstride + y0*ystride + z0*zstride;
192 dstptr = dst;
193 for (z = z0; z <= z1; z++) {
194 for (y = y0; y <= y1; y++) {
195 for (x = x0; x <= x1; x++) {
196 if (field_size == 1)
197 ByteField(dstptr, 0) = ByteField(voxel, field_offset);
198 else if (field_size == 2)
199 ShortField(dstptr, 0) = ShortField(voxel, field_offset);
200 else
201 IntField(dstptr, 0) = IntField(voxel, field_offset);
202 dstptr += dst_xstride;
203 voxel += xstride;
204 }
205 dstptr += dst_ystride - (x1-x0+1)*dst_xstride;
206 voxel += ystride - (x1-x0+1)*xstride;
207 }
208 dstptr += dst_zstride - (y1-y0+1)*dst_ystride;
209 voxel += zstride - (y1-y0+1)*ystride;
210 }
211 return(VP_OK);
212 }
213
214 /*
215 * ClassifyRawVolume
216 *
217 * Classify a portion of a raw volume, quantize the result, and store
218 * as an array of 8-bit opacities.
219 */
220
221 static int
ClassifyRawVolume(vpc,correct,x0,y0,z0,x1,y1,z1,dst,dst_xstride,dst_ystride,dst_zstride)222 ClassifyRawVolume(vpc, correct, x0, y0, z0, x1, y1, z1, dst,
223 dst_xstride, dst_ystride, dst_zstride)
224 vpContext *vpc; /* context */
225 int correct; /* if true then correct for view */
226 int x0, y0, z0; /* origin of extracted region */
227 int x1, y1, z1; /* opposite corner of extracted region */
228 unsigned char *dst; /* buffer to store result into */
229 int dst_xstride; /* stride (in bytes) for destination array */
230 int dst_ystride;
231 int dst_zstride;
232 {
233 float *opc;
234 int num_voxels;
235 int retcode;
236
237 /* check for errors */
238 if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
239 return(retcode);
240
241 /* compute opacities */
242 num_voxels = (x1-x0+1)*(y1-y0+1)*(z1-z0+1);
243 Alloc(vpc, opc, float *, num_voxels*sizeof(float), "opacity_block");
244 VPClassifyBlock(vpc, correct, x0, y0, z0, x1, y1, z1, opc,
245 sizeof(float), (x1-x0+1)*sizeof(float),
246 (x1-x0+1)*(y1-y0+1)*sizeof(float));
247
248 /* quantize opacities */
249 VPQuantize(opc, x1-x0+1, y1-y0+1, z1-z0+1, 255., 255, dst,
250 dst_xstride, dst_ystride, dst_zstride);
251
252 Dealloc(vpc, opc);
253 return(VP_OK);
254 }
255
256 /*
257 * ShadeRawVolume
258 *
259 * Shade a portion of a raw volume, quantize the result, and store
260 * as an array of 8-bit intensities.
261 */
262
263 static int
ShadeRawVolume(vpc,x0,y0,z0,x1,y1,z1,dst,dst_xstride,dst_ystride,dst_zstride)264 ShadeRawVolume(vpc, x0, y0, z0, x1, y1, z1, dst,
265 dst_xstride, dst_ystride, dst_zstride)
266 vpContext *vpc; /* context */
267 int x0, y0, z0; /* origin of extracted region */
268 int x1, y1, z1; /* opposite corner of extracted region */
269 unsigned char *dst; /* buffer to store result into */
270 int dst_xstride; /* stride (in bytes) for destination array */
271 int dst_ystride;
272 int dst_zstride;
273 {
274 float *shd;
275 int num_colors;
276 int retcode;
277 int xstride, ystride, zstride;
278
279 /* check for errors */
280 if ((retcode = VPCheckShader(vpc)) != VP_OK)
281 return(retcode);
282
283 /* compute colors */
284 num_colors = (x1-x0+1)*(y1-y0+1)*(z1-z0+1)*vpc->color_channels;
285 Alloc(vpc, shd, float *, num_colors*sizeof(float), "color_block");
286 xstride = vpc->color_channels * sizeof(float);
287 ystride = xstride * (x1-x0+1);
288 zstride = ystride * (y1-y0+1);
289 VPShadeBlock(vpc, x0, y0, z0, x1, y1, z1, shd, xstride, ystride, zstride);
290
291 /* quantize colors */
292 VPQuantize(shd, x1-x0+1, y1-y0+1, z1-z0+1, 1., 255, dst,
293 dst_xstride, dst_ystride, dst_zstride);
294
295 Dealloc(vpc, shd);
296 return(VP_OK);
297 }
298
299 /*
300 * VPClassifyBlock
301 *
302 * Classify a block of the current raw volume. The result is an
303 * array of floating point opacities in the range 0.0-1.0.
304 */
305
306 vpResult
VPClassifyBlock(vpc,correct,x0,y0,z0,x1,y1,z1,opc,dst_xstride,dst_ystride,dst_zstride)307 VPClassifyBlock(vpc, correct, x0, y0, z0, x1, y1, z1, opc,
308 dst_xstride, dst_ystride, dst_zstride)
309 vpContext *vpc; /* context */
310 int correct; /* if true then correct for view */
311 int x0, y0, z0; /* origin of extracted region */
312 int x1, y1, z1; /* opposite corner of extracted region */
313 float *opc; /* buffer to store result into */
314 int dst_xstride; /* stride (in bytes) for destination array */
315 int dst_ystride;
316 int dst_zstride;
317 {
318 unsigned char *voxel;
319 int xstride, ystride, zstride;
320 int x, y, z;
321 float opacity;
322 int quant_opc;
323 int retcode;
324
325 if (correct) {
326 if ((retcode = VPFactorView(vpc)) != VP_OK)
327 return(retcode);
328 }
329 xstride = vpc->xstride;
330 ystride = vpc->ystride;
331 zstride = vpc->zstride;
332 voxel = vpc->raw_voxels;
333 voxel += x0*xstride + y0*ystride + z0*zstride;
334 for (z = z0; z <= z1; z++) {
335 for (y = y0; y <= y1; y++) {
336 for (x = x0; x <= x1; x++) {
337 opacity = VPClassifyVoxel(vpc, voxel);
338 if (correct) {
339 quant_opc = opacity * 255.;
340 if (quant_opc > 255)
341 quant_opc = 255;
342 else if (quant_opc < 0)
343 quant_opc = 0;
344 opacity = CorrectOpacity(vpc, quant_opc, x, y, z);
345 }
346 *opc = opacity;
347 opc = (float *)((char *)opc + dst_xstride);
348 voxel += xstride;
349 }
350 opc = (float *)((char *)opc + dst_ystride - (x1-x0+1)*dst_xstride);
351 voxel += ystride - (x1-x0+1)*xstride;
352 }
353 opc = (float *)((char *)opc + dst_zstride - (y1-y0+1)*dst_ystride);
354 voxel += zstride - (y1-y0+1)*ystride;
355 }
356 return(VP_OK);
357 }
358
359 /*
360 * VPClassifyVoxel
361 *
362 * Classify a single voxel. Return value is an opacity.
363 */
364
365 float
VPClassifyVoxel(vpc,voxel)366 VPClassifyVoxel(vpc, voxel)
367 vpContext *vpc; /* context */
368 void *voxel; /* pointer to voxel */
369 {
370 int num_params; /* number of parameters to classifier */
371 int p; /* current parameter number */
372 int field; /* field for the parameter */
373 int field_size; /* size of the field */
374 int field_offset; /* offset for the field */
375 int index; /* index for table lookup */
376 float opacity; /* current value of the opacity */
377
378 num_params = vpc->num_clsfy_params;
379 opacity = 1;
380 for (p = 0; p < num_params; p++) {
381 /* get table index */
382 field = vpc->param_field[p];
383 field_offset = vpc->field_offset[field];
384 field_size = vpc->field_size[field];
385 index = VoxelField(voxel, field_offset, field_size);
386
387 /* load table value */
388 opacity *= vpc->clsfy_table[p][index];
389 }
390 return(opacity);
391 }
392
393 /*
394 * CorrectOpacity
395 *
396 * Correct an opacity for the current view.
397 * Return value is the corrected opacity.
398 */
399
400 static float
CorrectOpacity(vpc,quant_opc,x,y,z)401 CorrectOpacity(vpc, quant_opc, x, y, z)
402 vpContext *vpc; /* context */
403 int quant_opc; /* input opacity (0-255) */
404 int x, y, z; /* voxel coordinates in object space */
405 {
406 float opacity;
407
408 if (vpc->affine_view) {
409 opacity = vpc->affine_opac_correct[quant_opc];
410 } else {
411 /* XXX perspective rendering not available yet */
412 opacity = (float)quant_opc / (float)255.;
413 }
414 return(opacity);
415 }
416
417 /*
418 * VPShadeBlock
419 *
420 * Shade a block of the current raw volume. The result is an
421 * array of floating point colors in the range 0.0-255.0.
422 */
423
424 vpResult
VPShadeBlock(vpc,x0,y0,z0,x1,y1,z1,shd,dst_xstride,dst_ystride,dst_zstride)425 VPShadeBlock(vpc, x0, y0, z0, x1, y1, z1, shd,
426 dst_xstride, dst_ystride, dst_zstride)
427 vpContext *vpc; /* context */
428 int x0, y0, z0; /* origin of extracted region */
429 int x1, y1, z1; /* opposite corner of extracted region */
430 float *shd; /* buffer to store result into */
431 int dst_xstride; /* stride (in bytes) for destination array */
432 int dst_ystride;
433 int dst_zstride;
434 {
435 unsigned char *voxel;
436 int xstride, ystride, zstride;
437 int x, y, z;
438 int color_channels;
439
440 color_channels = vpc->color_channels;
441 xstride = vpc->xstride;
442 ystride = vpc->ystride;
443 zstride = vpc->zstride;
444 voxel = vpc->raw_voxels;
445 voxel += x0*xstride + y0*ystride + z0*zstride;
446 for (z = z0; z <= z1; z++) {
447 for (y = y0; y <= y1; y++) {
448 for (x = x0; x <= x1; x++) {
449 ShadeVoxel(vpc, voxel, x, y, z, shd);
450 shd = (float *)((char *)shd + dst_xstride);
451 voxel += xstride;
452 }
453 shd = (float *)((char *)shd + dst_ystride - (x1-x0+1)*dst_xstride);
454 voxel += ystride - (x1-x0+1)*xstride;
455 }
456 shd = (float *)((char *)shd + dst_zstride - (y1-y0+1)*dst_ystride);
457 voxel += zstride - (y1-y0+1)*ystride;
458 }
459 return(VP_OK);
460 }
461
462 /*
463 * ShadeVoxel
464 *
465 * Shade a voxel.
466 */
467
468 static void
ShadeVoxel(vpc,voxel,x,y,z,dst)469 ShadeVoxel(vpc, voxel, x, y, z, dst)
470 vpContext *vpc; /* context */
471 void *voxel; /* voxel data */
472 int x, y, z; /* voxel coordinates */
473 float *dst; /* storage for result (1 or 3 intensities, 0-255) */
474 {
475 int num_materials;
476 int color_channels;
477 int color_index_size, color_index_offset, color_index, color_table_offset;
478 int weight_index_size, weight_index_offset, weight_index;
479 int weight_table_offset;
480 int m;
481 float r, g, b;
482 float *color_table;
483 float *weight_table;
484
485 /* check shading mode */
486 if (vpc->shading_mode == CALLBACK_SHADER) {
487 if (vpc->color_channels == 1)
488 vpc->shade_func(voxel, dst, vpc->client_data);
489 else
490 vpc->shade_func(voxel, dst, dst+1, dst+2, vpc->client_data);
491 return;
492 } else if (vpc->shading_mode != LOOKUP_SHADER) {
493 VPBug("unknown shader type");
494 }
495
496 /* compute table indices */
497 num_materials = vpc->num_materials;
498 color_channels = vpc->color_channels;
499 color_index_size = vpc->field_size[vpc->color_field];
500 color_index_offset = vpc->field_offset[vpc->color_field];
501 color_index = VoxelField(voxel, color_index_offset, color_index_size);
502 color_table_offset = color_index * num_materials;
503 weight_index_size = vpc->field_size[vpc->weight_field];
504 weight_index_offset = vpc->field_offset[vpc->weight_field];
505 weight_index = VoxelField(voxel, weight_index_offset, weight_index_size);
506 weight_table_offset = weight_index * num_materials;
507
508 /* look up values in tables */
509 if (color_channels == 1) {
510 color_table = vpc->shade_color_table + color_table_offset;
511 weight_table = vpc->shade_weight_table + weight_table_offset;
512 if (num_materials == 1) {
513 r = *color_table;
514 } else {
515 r = 0;
516 for (m = 0; m < num_materials; m++)
517 r += *color_table++ * *weight_table++;
518 }
519 *dst = r;
520 } else {
521 color_table = vpc->shade_color_table + 3*color_table_offset;
522 weight_table = vpc->shade_weight_table + weight_table_offset;
523 if (num_materials == 1) {
524 r = *color_table++;
525 g = *color_table++;
526 b = *color_table;
527 } else {
528 r = 0;
529 g = 0;
530 b = 0;
531 for (m = 0; m < num_materials; m++) {
532 r += *color_table++ * *weight_table;
533 g += *color_table++ * *weight_table;
534 b += *color_table++ * *weight_table;
535 }
536 }
537 dst[0] = r;
538 dst[1] = g;
539 dst[2] = b;
540 }
541 }
542
543 /*
544 * VPQuantize
545 *
546 * Quantize a floating point array and store the result in a byte array.
547 */
548
549 void
VPQuantize(src,xlen,ylen,zlen,scale,maxvalue,dst,dst_xstride,dst_ystride,dst_zstride)550 VPQuantize(src, xlen, ylen, zlen, scale, maxvalue, dst,
551 dst_xstride, dst_ystride, dst_zstride)
552 float *src; /* floating point array */
553 int xlen, ylen, zlen; /* array dimensions */
554 double scale; /* scale to apply to each array element */
555 int maxvalue; /* clamp each array element to this value */
556 unsigned char *dst; /* store results here */
557 int dst_xstride; /* stride (in bytes) for destination array */
558 int dst_ystride;
559 int dst_zstride;
560 {
561 int value;
562 int x, y, z;
563
564 for (z = 0; z < zlen; z++) {
565 for (y = 0; y < ylen; y++) {
566 for (x = 0; x < xlen; x++) {
567 value = (int)rint(*src++ * scale);
568 if (value > maxvalue)
569 value = maxvalue;
570 else if (value < 0)
571 value = 0;
572 *dst = value;
573 dst += dst_xstride;
574 }
575 dst += dst_ystride - xlen*dst_xstride;
576 }
577 dst += dst_zstride - ylen*dst_ystride;
578 }
579 }
580
581 /*
582 * ExtractClassifiedVolume
583 *
584 * Extract a field from a classified volume into an array.
585 */
586
587 static int
ExtractClassifiedVolume(vpc,axis,x0,y0,z0,x1,y1,z1,field,dst,dst_xstride,dst_ystride,dst_zstride)588 ExtractClassifiedVolume(vpc, axis, x0, y0, z0, x1, y1, z1, field, dst,
589 dst_xstride, dst_ystride, dst_zstride)
590 vpContext *vpc; /* context */
591 int axis; /* which axis to extract from */
592 int x0, y0, z0; /* origin of extracted region */
593 int x1, y1, z1; /* opposite corner of extracted region */
594 int field; /* field to extract */
595 void *dst; /* buffer to store result into */
596 int dst_xstride; /* stride (in bytes) for destination array */
597 int dst_ystride;
598 int dst_zstride;
599 {
600 int i, j, k; /* voxel coordinates in rotated object space */
601 int i0, j0, k0; /* origin of extracted region */
602 int i1, j1, k1; /* opposite corner of extracted region */
603 int dst_istride; /* stride (in bytes) for destination array */
604 int dst_jstride;
605 int dst_kstride;
606 int ilen, jlen, klen; /* volume size */
607 RLEVoxels *rle_voxels; /* run-length encoded, classified volume */
608 unsigned char *voxel; /* pointer to current voxel in volume */
609 unsigned char *dstptr; /* pointer to destination */
610 unsigned char *length; /* pointer to current run length */
611 int run_length; /* length of current run */
612 int is_non_zero; /* true if current run is nonzero */
613 int rle_bytes_per_voxel; /* size of unclassified voxel */
614 int value; /* value of parameter for current voxel */
615 ScanOffset *slice_runs; /* offsets to start of runs for a slice */
616 int field_size; /* size of field in bytes */
617 int field_offset; /* byte offset for voxel field */
618
619 /* initialize */
620 switch (axis) {
621 case VP_X_AXIS:
622 rle_voxels = vpc->rle_x;
623 i0 = y0; j0 = z0; k0 = x0; i1 = y1; j1 = z1; k1 = x1;
624 dst_istride = dst_ystride;
625 dst_jstride = dst_zstride;
626 dst_kstride = dst_xstride;
627 break;
628 case VP_Y_AXIS:
629 rle_voxels = vpc->rle_y;
630 i0 = z0; j0 = x0; k0 = y0; i1 = z1; j1 = x1; k1 = y1;
631 dst_istride = dst_zstride;
632 dst_jstride = dst_xstride;
633 dst_kstride = dst_ystride;
634 break;
635 case VP_Z_AXIS:
636 rle_voxels = vpc->rle_z;
637 i0 = x0; j0 = y0; k0 = z0; i1 = x1; j1 = y1; k1 = z1;
638 dst_istride = dst_xstride;
639 dst_jstride = dst_ystride;
640 dst_kstride = dst_zstride;
641 break;
642 default:
643 return(VPSetError(vpc, VPERROR_BAD_OPTION));
644 }
645 if (rle_voxels == NULL)
646 return(VPSetError(vpc, VPERROR_BAD_VOLUME));
647 if (rle_voxels->scan_offsets_per_slice < 1)
648 return(VPSetError(vpc, VPERROR_BAD_VOLUME));
649 ilen = rle_voxels->ilen;
650 jlen = rle_voxels->jlen;
651 klen = rle_voxels->klen;
652 rle_bytes_per_voxel = vpc->rle_bytes_per_voxel;
653 if (field == VP_OPACITY_FIELD || field == VP_CORRECTED_OPAC_FIELD) {
654 field_size = 1;
655 field_offset = rle_bytes_per_voxel - 1;
656 } else {
657 field_size = vpc->field_size[field];
658 field_offset = vpc->field_offset[field];
659 }
660
661 /* extract slice */
662 dstptr = dst;
663 for (k = k0; k <= k1; k++) {
664 slice_runs = &rle_voxels->scan_offsets[k *
665 rle_voxels->scan_offsets_per_slice];
666 voxel = (unsigned char *)rle_voxels->data + slice_runs->first_data;
667 length = rle_voxels->run_lengths + slice_runs->first_len;
668 run_length = 0;
669 is_non_zero = 1;
670 for (j = 0; j < jlen; j++) {
671 for (i = 0; i < ilen; i++) {
672 while (run_length == 0) {
673 run_length = *length++;
674 is_non_zero = !is_non_zero;
675 }
676 run_length--;
677 if (i >= i0 && i <= i1 && j >= j0 && j <= j1) {
678 if (is_non_zero) {
679 if (field_size == 1)
680 ByteField(dstptr, 0) = ByteField(voxel,
681 field_offset);
682 else if (field_size == 2)
683 ShortField(dstptr, 0) = ShortField(voxel,
684 field_offset);
685 else
686 IntField(dstptr, 0) = IntField(voxel,field_offset);
687 voxel += rle_bytes_per_voxel;
688 } else {
689 if (field_size == 1)
690 ByteField(dstptr, 0) = 0;
691 else if (field_size == 2)
692 ShortField(dstptr, 0) = 0;
693 else
694 IntField(dstptr, 0) = 0;
695 }
696 dstptr += dst_istride;
697 } else {
698 if (is_non_zero)
699 voxel += rle_bytes_per_voxel;
700 }
701 }
702 if (j >= j0 && j <= j1)
703 dstptr += dst_jstride - (i1-i0+1)*dst_istride;
704 }
705 dstptr += dst_kstride - (j1-j0+1)*dst_jstride;
706 }
707 return(VP_OK);
708 }
709