1 /***********************************************************************/
2 /* Open Visualization Data Explorer                                    */
3 /* (C) Copyright IBM Corp. 1989,1999                                   */
4 /* ALL RIGHTS RESERVED                                                 */
5 /* This code licensed under the                                        */
6 /*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
7 /***********************************************************************/
8 
9 #include <dxconfig.h>
10 
11 
12 /*
13  * Header: /usr/people/gresh/code/svs/src/libdx/RCS/volume.c,v 5.0 92/11/12 09:09:08 svs Exp Locker: gresh
14  */
15 
16 
17 #include <string.h>
18 #include <stdlib.h>
19 #include <dx/dx.h>
20 #include "render.h"
21 #include <stdio.h>
22 
23 /*
24  * Make a sort routine called facesort()
25  */
26 
27 #define TYPE struct sortindex
28 #define LT(a,b) ((a)->z < (b)->z)
29 #define GT(a,b) ((a)->z > (b)->z)
30 #define QUICKSORT facesort
31 #include "qsort.c"
32 
33 
34 /*
35  * Irregular.  At this point, _dxf_Gather (called from Tile) has gathered
36  * into gather->sort all the volume and translucent faces.  We must
37  * now sort them and call _dxf_Triangle to render them.
38  *
39  * XXX - a better sort algorithm would be desirable here
40  *
41  * The size of the elements that go into the sort array depends on such
42  * factors as the element type (triangle, quad etc.); for faces derived
43  * from connection-dependent volume elements, there is an additional word
44  * naming the element that the face came from to identify the color.  This
45  * size appears in SIZE macro in shade.c, the INFO macro in gather.c,
46  * and the ADVANCE macro in volume.c, and these must be kept in sync.
47  */
48 
49 #define GET(s,n) (((int *)((long)(s) + sizeof(struct sort)))[n])
50 
51 #define ADVANCE(n) \
52     (s = (struct sort *)((long)s + sizeof(struct sort) + n * sizeof(int)))
53 
54 #define IMMEDIATE(s,type) ((type *)((long)(s) + sizeof(struct sort)))
55 
56 Error
_dxf_VolumeIrregular(struct buffer * b,struct gather * gather,int clip)57 _dxf_VolumeIrregular(struct buffer *b, struct gather *gather, int clip)
58 {
59     struct sortindex *sortindex = NULL, *si;
60     struct sort *s;
61     struct xfield *xf;
62     Quadrilateral quad;
63     Triangle tri;
64     Line line;
65     Point *p;
66     float z;
67     int i, j, g, n;
68 
69     /* anything to do? */
70     if (!gather->nthings)
71 	return OK;
72 
73     /* allocate buffer for sort index array */
74     sortindex = (struct sortindex *)
75 	DXAllocateLocal(gather->nthings * sizeof(struct sortindex));
76     if (!sortindex) {
77 	DXResetError();
78 	sortindex = (struct sortindex *)
79 	    DXAllocate(gather->nthings * sizeof(struct sortindex));
80 	if (!sortindex)
81 	    goto error;
82 	DXWarning("sort index array (%d bytes) "
83 		"is too big for local memory",
84 		gather->nthings * sizeof(struct sortindex));
85     }
86 
87     /* compute z values, sort the faces */
88     for (i=0, s=gather->sort; s<gather->current; i++) {
89 	xf = &gather->fields[s->field];
90 	p = xf->positions;
91 	sortindex[i].sort = s;
92 	switch (s->type) {
93 	case SORT_POINT:
94 	    z = p[GET(s,0)].z;
95 	    ADVANCE(1);
96 	    break;
97 	case SORT_LINE:
98 	    line = xf->c.lines[GET(s,0)];
99 	    z = (((double)p[line.p].z) + ((double)p[line.q].z)) / 2.0;
100 	    ADVANCE(1);
101 	    break;
102 	case SORT_POLYLINE:
103 	    line.p = xf->edges[GET(s,1)];
104 	    line.q = xf->edges[GET(s,1)+1];
105 	    z = (((double)p[line.p].z) + ((double)p[line.q].z)) / 2.0;
106 	    ADVANCE(2);
107 	    break;
108 	case SORT_TRI:
109 	    tri = xf->c.triangles[GET(s,0)];
110 	    z = (((double)p[tri.p].z) + ((double)p[tri.q].z) +
111 		 ((double)p[tri.r].z)) / 3.0;
112 	    ADVANCE(1);
113 	    break;
114 	case SORT_QUAD:
115 	    quad = xf->c.quads[GET(s,0)];
116 	    z = (((double)p[quad.p].z) + ((double)p[quad.q].z) +
117 		 ((double)p[quad.r].z) + ((double)p[quad.s].z)) / 4.0;
118 	    ADVANCE(1);
119 	    break;
120 	case SORT_INV_TRI_PTS:
121 	case SORT_TRI_PTS:
122 	    z = (((double)p[GET(s,0)].z) + ((double)p[GET(s,1)].z) +
123 	         ((double)p[GET(s,2)].z)) / 3.0;
124 	    ADVANCE(3);
125 	    break;
126 	case SORT_INV_TRI_FLAT:
127 	case SORT_TRI_FLAT:
128 	    z = (((double)p[GET(s,0)].z) + ((double)p[GET(s,1)].z) +
129 	   	 ((double)p[GET(s,2)].z)) / 3.0;
130 	    ADVANCE(4);
131 	    break;
132 	case SORT_INV_QUAD_PTS:
133 	case SORT_QUAD_PTS:
134 	    z = (((double)p[GET(s,0)].z) + ((double)p[GET(s,1)].z) +
135 		 ((double)p[GET(s,2)].z) + ((double)p[GET(s,3)].z)) / 4.0;
136 	    ADVANCE(4);
137 	    break;
138 	case SORT_INV_QUAD_FLAT:
139 	case SORT_QUAD_FLAT:
140 	    z = (((double)p[GET(s,0)].z) + ((double)p[GET(s,1)].z) +
141 		 ((double)p[GET(s,2)].z) + ((double)p[GET(s,3)].z)) / 4.0;
142 	    ADVANCE(5);
143 	    break;
144 	default:
145 	    DXErrorReturn(ERROR_INTERNAL, "unknown type in sort array");
146 	}
147 	sortindex[i].z = z;
148     }
149     DXMarkTimeLocal("start sort");
150     facesort(sortindex, gather->nthings);
151     DXMarkTimeLocal("end sort");
152 
153     /* render them in sorted order */
154     for (i=0, si=sortindex, n=gather->nthings; i<n; i++, si++) {
155 
156 	s = si->sort;
157 	xf = &gather->fields[s->field];
158 
159 	ASSERT(xf->opacities || xf->volume);
160 
161 	/*
162 	 * Merge identical surface faces from different fields into one
163 	 * non-surface face.  This avoids artifacts that occur when surface
164 	 * faces, e.g. in partitions, are mis-sorted.  (Identical is taken
165 	 * to mean identical z sort values, face type, x sum).
166 	 */
167 
168 
169 #define SUM(s,xf,x) (							\
170     p = xf->positions,							\
171     type==SORT_QUAD_PTS?						\
172 	p[GET(s,0)].x + p[GET(s,1)].x + p[GET(s,2)].x + p[GET(s,3)].x :	\
173 	p[GET(s,0)].x + p[GET(s,1)].x + p[GET(s,2)].x			\
174 )
175 
176 #define IGNORE (2)
177 
178 	if (s->surface==IGNORE)
179 	    continue;
180 	if (s->surface) {
181 	    int type = s->type;
182 	    if (type==SORT_QUAD_PTS || type==SORT_TRI_PTS) {
183 		float x = SUM(s,xf,x);
184 		float y = SUM(s,xf,y);
185 		float z = si->z;
186 		struct sortindex *si1;
187 		Point *p;
188 		for (j=i+1, si1=si+1;  j<n && si1->z==z;  j++, si1++) {
189 		    struct sort *s1 = si1->sort;
190 		    struct xfield *xf1 = &gather->fields[s1->field];
191 		    if (s1->surface && xf1!=xf && s1->type==type) {
192 			float x1 = SUM(s1,xf1,x);
193 			float y1 = SUM(s1,xf1,y);
194 			if (x==x1 && y==y1) {
195 			    s->surface = 0;
196 			    s1->surface = IGNORE;
197 			}
198 		    }
199 		}
200 	    }
201 	}
202 
203 #define FCOLORP(i) (Pointer)(xf->cmap? \
204     (Pointer)&(((unsigned char *)xf->fcolors)[xf->fcst ? 0 : i]) : \
205     (Pointer)&(((RGBColor *)xf->fcolors)[xf->fcst ? 0 : i]))
206 
207 #define BCOLORP(i) (Pointer)(xf->cmap? \
208     (Pointer)&(((unsigned char *)xf->bcolors)[xf->fcst ? 0 : i]) : \
209     (Pointer)&(((RGBColor *)xf->bcolors)[xf->fcst ? 0 : i]))
210 
211 #define OPACITYP(i) (Pointer)(xf->omap? \
212     (Pointer)&(((unsigned char *)xf->opacities)[xf->ocst ? 0 : i]) : \
213     xf->opacities? (Pointer)&(((float *)xf->opacities)[xf->ocst ? 0 : i]) : NULL)
214 
215 
216 	switch (s->type) {
217 	case SORT_POINT:
218 	    if (!_dxf_Point(b, xf, GET(s,0)))
219 		goto error;
220 	    break;
221 	case SORT_LINE:
222 	    g = GET(s,0);
223 	    if (xf->colors_dep == dep_positions) {
224 		if (!_dxf_Line(b, xf, 1, xf->c.lines, &g, clip, INV_VALID))
225 		    goto error;
226 	    } else if (xf->colors_dep == dep_connections) {
227 		if (!_dxf_LineFlat(b, xf, 1, xf->c.lines, &g,
228 				xf->fcolors, xf->opacities,
229 				clip, INV_VALID))
230 		    goto error;
231 	    }
232 	    break;
233 	case SORT_POLYLINE:
234             line.p = xf->edges[GET(s,1)];
235 	    line.q = xf->edges[GET(s,1)+1];
236 	    if (xf->colors_dep == dep_positions) {
237 		if (!_dxf_Line(b, xf, 1, &line, NULL, clip, INV_VALID))
238 		    goto error;
239 	    } else if (xf->colors_dep == dep_polylines) {
240 		if (!_dxf_LineFlat(b, xf, 1, &line, NULL,
241 			       xf->fcolors, xf->opacities,
242 			       clip, INV_VALID))
243 		    goto error;
244 	    }
245 	    break;
246 	case SORT_TRI:
247 	    g = GET(s,0);
248 	    if (xf->colors_dep == dep_positions || xf->lights) {
249 		if (!_dxf_Triangle(b, xf, 1, xf->c.triangles+g, &g,
250 					s->surface, clip, INV_VALID))
251 		    goto error;
252 	    } else if (xf->colors_dep == dep_connections) {
253 		if (!_dxf_TriangleFlat(b, xf, 1, xf->c.triangles+g, NULL,
254 				   FCOLORP(g), BCOLORP(g),
255 				   OPACITYP(g), s->surface, clip, INV_VALID))
256 		    goto error;
257 	    }
258 	    break;
259 	case SORT_QUAD:
260 	    g = GET(s,0);
261 	    if (xf->colors_dep == dep_positions || xf->lights) {
262 		if (!_dxf_Quad(b, xf, 1, xf->c.quads+g, &g,
263 						s->surface, clip, INV_VALID))
264 		    goto error;
265 	    } else if (xf->colors_dep == dep_connections) {
266 		if (!_dxf_QuadFlat(b, xf, 1, xf->c.quads+g, NULL,
267 			       FCOLORP(g), BCOLORP(g),
268 			       OPACITYP(g), s->surface, clip, INV_VALID))
269 		    goto error;
270 	    }
271 	    break;
272 	case SORT_TRI_PTS:
273 	    if (!_dxf_Triangle(b, xf, 1, IMMEDIATE(s,Triangle), NULL,
274 						s->surface, clip, INV_VALID))
275 		goto error;
276 	    break;
277 	case SORT_INV_TRI_PTS:
278 	    if (!_dxf_Triangle(b, xf, 1, IMMEDIATE(s,Triangle), NULL,
279 						s->surface, clip, INV_INVALID))
280 		goto error;
281 	    break;
282 	case SORT_TRI_FLAT:
283 	    g = GET(s,3);
284 	    if (!_dxf_TriangleFlat(b, xf, 1, IMMEDIATE(s,Triangle), NULL,
285 			       FCOLORP(g), BCOLORP(g),
286 			       OPACITYP(g), s->surface, clip, INV_VALID))
287 		goto error;
288 	    break;
289 	case SORT_INV_TRI_FLAT:
290 	    g = GET(s,3);
291 	    if (!_dxf_TriangleFlat(b, xf, 1, IMMEDIATE(s,Triangle), NULL,
292 			       FCOLORP(g), BCOLORP(g),
293 			       OPACITYP(g), s->surface, clip, INV_INVALID))
294 		goto error;
295 	    break;
296 	case SORT_QUAD_PTS:
297 	    if (!_dxf_Quad(b, xf, 1, IMMEDIATE(s,Quadrilateral), NULL,
298 			    s->surface, clip, INV_VALID))
299 		goto error;
300 	    break;
301 	case SORT_INV_QUAD_PTS:
302 	    if (!_dxf_Quad(b, xf, 1, IMMEDIATE(s,Quadrilateral), NULL,
303 			    s->surface, clip, INV_INVALID))
304 		goto error;
305 	    break;
306 	case SORT_QUAD_FLAT:
307 	    g = GET(s,4);
308 	    if (!_dxf_QuadFlat(b, xf, 1, IMMEDIATE(s,Quadrilateral), NULL,
309 			   FCOLORP(g), BCOLORP(g),
310 			   OPACITYP(g), s->surface, clip, INV_VALID))
311 		goto error;
312 	    break;
313 	case SORT_INV_QUAD_FLAT:
314 	    g = GET(s,4);
315 	    if (!_dxf_QuadFlat(b, xf, 1, IMMEDIATE(s,Quadrilateral), NULL,
316 			   FCOLORP(g), BCOLORP(g),
317 			   OPACITYP(g), s->surface, clip, INV_INVALID))
318 		goto error;
319 	    break;
320 	default:
321 	    DXErrorReturn(ERROR_INTERNAL, "unknown type in sort array");
322 	}
323     }
324 
325     DXFree((Pointer)sortindex);
326     return OK;
327 
328 error:
329     DXFree((Pointer)sortindex);
330     return ERROR;
331 }
332 
333 
334 /*
335  * Some macros to help us access counts (Ki), strides (Si), and
336  * delta vectors (Di), and also the same for the indices into
337  * connection-dependent colors.
338  */
339 
340 #define D0 xf->d[0]
341 #define D1 xf->d[1]
342 #define D2 xf->d[2]
343 
344 #define K0 xf->k[0]
345 #define K1 xf->k[1]
346 #define K2 xf->k[2]
347 
348 #define S0 xf->s[0]
349 #define S1 xf->s[1]
350 #define S2 xf->s[2]
351 
352 
353 #define CK0 xf->ck[0]
354 #define CK1 xf->ck[1]
355 #define CK2 xf->ck[2]
356 
357 #define CS0 xf->cs[0]
358 #define CS1 xf->cs[1]
359 #define CS2 xf->cs[2]
360 
361 
362 /*
363  *
364  */
365 
366 #define FLOOR(a) (((int)((a)+100000.0)-100000))
367 #define CEIL(a) (((int)((a)-100000.0)+100000))
368 #define ROUND(a) FLOOR((a)+0.5)
369 #define ABS(a) ((a)<0? -(a) : (a))
370 #define DOTXY(u,v) (u.x*v.x+u.y*v.y)
371 
372 static
373 void
reduce(Vector u,Vector v,Vector * up,Vector * vp)374 reduce(Vector u, Vector v, Vector *up, Vector *vp)
375 {
376     int i, j;
377     float d;
378     do {
379 	d = DOTXY(u,u);
380 	if (!d) break;
381 	d = DOTXY(u,v) / d;
382 	i = d>0? (int)(d+0.4999) : (int)(d-0.4999);
383 	v = DXSub(v, DXMul(u,i));
384 	d = DOTXY(v,v);
385 	if (!d) break;
386 	d = DOTXY(v,u) / d;
387 	j = d>0? (int)(d+0.4999) : (int)(d-0.4999);
388 	u = DXSub(u, DXMul(v,j));
389     } while (i!=0 || j!=0);
390     *up = u, *vp = v;
391 }
392 
393 static
394 float
side(Vector u,Vector v)395 side(Vector u, Vector v)
396 {
397     float a, b, minx, maxx, miny, maxy;
398 #   define mm(MAX,x) (a = MAX(0,u.x), b = MAX(v.x,u.x+v.x), MAX(a,b))
399     maxx = mm(MAX,x);
400     minx = mm(MIN,x);
401     maxy = mm(MAX,y);
402     miny = mm(MIN,y);
403     return MAX(maxx-minx, maxy-miny);
404 }
405 
406 static
407 float
area(Vector u,Vector v)408 area(Vector u, Vector v)
409 {
410     float a = u.x*v.y - u.y*v.x;
411     return ABS(a);
412 }
413 
414 
415 /*
416  * A number of the algorithms depend on enumerating planes of voxels
417  * or voxel faces that are the most parallel to the image plane.  This
418  * is judged by comparing the delta z's between the planes, as measured
419  * in the direction of the camera view (i.e. perpendicular to the
420  * image plane).  The orient() routine accomplishes this by reordering
421  * the axes so that the number 2 axis has the smallest delta z between planes,
422  * and returns that delta z.  After calling orient(), the rendering algorithm
423  * can then just be written to iterate through axis 2 in the outermost loop.
424  */
425 
426 static void
swap(struct xfield * xf,float * den,int a,int b)427 swap(struct xfield *xf, float *den, int a, int b) {
428     if (ABS(den[a]) > ABS(den[b])) {
429 	{float t; t=den[a]; den[a]=den[b]; den[b]=t;}
430 	{Vector t; t=xf->d[a]; xf->d[a]=xf->d[b]; xf->d[b]=t;}
431 	{int t; t=xf->k[a]; xf->k[a]=xf->k[b]; xf->k[b]=t;}
432 	{int t; t=xf->s[a]; xf->s[a]=xf->s[b]; xf->s[b]=t;}
433 	{int t; t=xf->ck[a]; xf->ck[a]=xf->ck[b]; xf->ck[b]=t;}
434 	{int t; t=xf->cs[a]; xf->cs[a]=xf->cs[b]; xf->cs[b]=t;}
435     }
436 }
437 
438 static
439 float
orient(struct xfield * xf,int print)440 orient(struct xfield *xf, int print)
441 {
442     float num, den[3];
443 
444     num = D0.x*D1.y*D2.z + D1.x*D2.y*D0.z + D2.x*D0.y*D1.z
445 	- D2.x*D1.y*D0.z - D0.x*D2.y*D1.z - D1.x*D0.y*D2.z;
446     den[0] = D1.x*D2.y - D2.x*D1.y;
447     den[1] = D2.x*D0.y - D0.x*D2.y;
448     den[2] = D0.x*D1.y - D1.x*D0.y;
449 
450     if (print)
451 	DXDebug("R", "%.3f / %.3f,%.3f,%.3f", num, den[0], den[1], den[2]);
452 
453     /* bubble largest den[i] (smallest dz) to number 2 position */
454     swap(xf, den, 0, 1);
455     swap(xf, den, 1, 2);
456 
457     if (print)
458 	DXDebug("R", "%.3f / %.3f,%.3f,%.3f", num, den[0], den[1], den[2]);
459 
460     return num / den[2];
461 }
462 
463 
464 /*
465  * Regular volume rendering.  In this case, _dxf_Gather only gathered the
466  * fields in gather->fields, because Tile set gather->sort to NULL.
467  * For each field, we determine the backmost corner and deltas and
468  * strides relative to that corner, and then sort the fields wrt the
469  * z value of the backmost corner.  Then we enumerate the fields in
470  * order of the z value of the backmost corner.  This is guaranteed
471  * to be the correct order for partitioned regular fields.  For each
472  * field we then paint the field according to one of a number of
473  * algorithms:
474  *
475  * 0 - automatically choose one
476  * 1 - all faces using irregular element algorithm
477  * 2 - one plane direction (plus edges) using irregular element algorithm
478  * 3 - composite planes, using smooth shaded quads
479  * 4 - composite planes, using flat shaded quads w/ average of four corners
480  * 5 - composite planes, using flat shaded quads w/ corner value
481  * 6 - patch
482  * 7 - pixel average
483  */
484 
485 static
486 int
487 #ifdef OS2
488 _Optlink
489 #endif
compare(const void * p1,const void * p2)490 compare(const void *p1, const void * p2)
491 {
492     const struct xfield **a = (const struct xfield **) p1;
493     const struct xfield **b = (const struct xfield **) p2;
494     if ((*a)->o.z < (*b)->o.z)
495 	return -1;
496     else if ((*a)->o.z > (*b)->o.z)
497 	return 1;
498     else return 0;
499 }
500 
501 static
502 void
step(struct xfield * xf,int i)503 step(struct xfield *xf, int i)		    /* step along axis i: */
504 {
505     Vector d, m;
506     d = xf->d[i];
507     if (d.z < 0) {			    /* if it takes us towards back */
508 	m = DXMul(d, xf->k[i]-1);		    /* delta to move to new origin */
509 	xf->o = DXAdd(xf->o, m);		    /* move origin */
510 	xf->d[i] = DXNeg(d);		    /* adjust delta vector */
511 	xf->n += (xf->k[i]-1) * xf->s[i];   /* adjust origin index */
512 	xf->s[i] = -xf->s[i];		    /* adjust stride */
513 	xf->cn += (xf->ck[i]-1) * xf->cs[i];/* same thing for */
514 	xf->cs[i] = -xf->cs[i];		    /* connection-dep indices */
515     }
516 }
517 
518 
519 /*
520  *
521  */
522 
523 /* --CHECKME-- the following was prototyped but never defined
524 static Error VolumeRegularQuad(struct buffer *, struct xfield *, int);
525  */
526 
527 static Error VolumeRegularFace(struct buffer *, struct xfield *, int);
528 static Error VolumeRegularFace3(struct buffer *, struct xfield *, int);
529 static Error VolumeRegularPlane(struct buffer *, struct xfield *, int);
530 Error
_dxf_VolumeRegular(struct buffer * b,struct gather * gather,int clip)531 _dxf_VolumeRegular(struct buffer *b, struct gather *gather, int clip)
532 {
533     int i /*, patch_size = 0*/;
534     struct xfield **xfields, *xf = NULL;
535     Vector u, v;
536     float a, s;
537 
538     /* anything to do? */
539     if (!gather->nfields)
540 	return OK;
541 
542     /* allocate an array of pointers to xfields */
543     xfields = (struct xfield **)
544 	DXAllocateLocal(gather->nfields * sizeof(struct xfield *));
545     if (!xfields)
546 	goto error;
547 
548     /*
549      * Collect counts (->k[i]) for each axis; compute strides (->s[i])
550      * from counts; get origin (->o), get delta vectors (->d[i]) for each
551      * axis; step along edges of cube to get to back corner (new origin
552      * ->o and point index of this origin ->n), adjusting strides and
553      * deltas as necessary.
554      */
555     for (i=0; i<gather->nfields; i++) {
556 
557 	int n;
558 
559 	/* xf */
560 	xf = &gather->fields[i];
561 	xfields[i] = xf;
562 
563 	/* perspective? */
564 	if (xf->tile.perspective)
565 	    DXErrorReturn(ERROR_BAD_PARAMETER,
566 			"perspective volume rendering is not supported");
567 
568 	/* counts and strides */
569 	DXQueryGridConnections(xf->connections_array, &n, xf->k);
570 	ASSERT(n==3);
571 	S2 = 1;				    /* stride for axis 2 */
572 	S1 = K2;			    /* stride for axis 1 */
573 	S0 = K2 * K1;			    /* stride for axis 0 */
574 
575 
576 	/* for conn-dep colors */
577         CK0 = K0 - 1;			    /* one fewer in ea dim */
578         CK1 = K1 - 1;			    /* for conn-dep colors */
579         CK2 = K2 - 1;
580 	CS2 = 1;			    /* stride for axis 2 */
581 	CS1 = CK2;			    /* stride for axis 1 */
582 	CS0 = CK2 * CK1;		    /* stride for axis 0 */
583 
584 	/* origin and delta vectors */
585 	DXQueryGridPositions(xf->positions_array, &n, NULL,
586 			   (float *)(&xf->o), (float *)(xf->d));
587 	ASSERT(n==3);
588 
589 	/* step along axes to get to back corner */
590 	step(xf, 0);
591 	step(xf, 1);
592 	step(xf, 2);
593     }
594 
595     /* sort the fields by z value of back corner */
596     qsort((char *)xfields, gather->nfields, sizeof(*xfields), compare);
597 
598     /* compute parameters */
599     if (xf->volAlg == 4)
600     {
601 	xf = xfields[0];
602 	orient(xf, 1);
603 	a = area(D0, D1);
604 	s = side(D0, D1);
605 	DXDebug("R", "(%.2f,%.2f) x (%.2f,%.2f): area %.2f, side %.2f",
606 		    D0.x, D0.y, D1.x, D1.y, a, s);
607 	reduce(D0, D1, &u, &v);
608 	a = area(u, v);
609 	s = side(u, v);
610 	DXDebug("R", "(%.2f,%.2f) x (%.2f,%.2f): area %.2f, side %.2f",
611 	      u.x, u.y, v.x, v.y, a, s);
612 	/*patch_size = CEIL(s);*/
613 	/* XXX - warp doesn't do z yet, or color multipliers */
614 	if (s<=1)
615 	    b->vol = /*"warp"*/ "plane";
616 	else if (a<50)
617 	    b->vol = "plane";
618 	else if (a < 75)
619 	    b->vol = "face";
620 	else
621 	    b->vol = "face3";
622     }
623     else if (xf->volAlg == 1)
624 	b->vol = "face";
625     else if (xf->volAlg == 2)
626 	b->vol = "face3";
627     else
628 	b->vol = "plane";
629 
630     /*
631      * If there are invalid connections, cannot use plane algorithm
632      */
633     if (xf->iElts && !strcmp(b->vol, "plane"))
634        b->vol = "face";
635 
636     /* for each field, in sorted order, call appropriate rendering routine */
637     for (i=0; i<gather->nfields; i++) {
638 	struct xfield *xf = xfields[i];
639 	Error rc;
640 	if (strcmp(b->vol,"face")==0)
641 	    rc = VolumeRegularFace(b, xf, clip);
642 	else if (strcmp(b->vol,"face3")==0)
643 	    rc = VolumeRegularFace3(b, xf, clip);
644 	else if (strcmp(b->vol,"plane")==0)
645 	    rc = VolumeRegularPlane(b, xf, clip);
646 	else
647 	    DXErrorGoto(ERROR_BAD_PARAMETER, "unknown volume rendering method");
648 	if (xf->positions_local) {
649 	    DXFreeArrayDataLocal(xf->positions_array, (Pointer)xf->positions);
650     	    xf->positions = NULL;
651 	    xf->positions_local = 0;
652 	}
653 	if (!rc) goto error;
654     }
655 
656     DXFree((Pointer)xfields);
657     return OK;
658 
659 error:
660     DXFree((Pointer)xfields);
661     return ERROR;
662 }
663 
664 
665 /*
666  * Looping macros
667  */
668 
669 #define FOR2(k) for (i2=0,n2=xf->n; i2<k;  i2++,n2+=S2)
670 #define FOR1(k) for (i1=0,n1=n2;    i1<k;  i1++,n1+=S1)
671 #define FOR0(k) for (i0=0,n =n1;    i0<k;  i0++,n +=S0)
672 
673 
674 /*
675  * DXRender all faces using the same face-based algorithm as for
676  * irregular grids
677  */
678 
679 /* XXX - for now */
680 #define EXPAND_POSITIONS						    	\
681     xf->positions = (Point *) DXGetArrayDataLocal(xf->positions_array);	    	\
682     xf->positions_local = 1;						        \
683     if (!xf->positions)								\
684 	return ERROR;
685 
686 #define Q(o, i, j, surface) {							\
687 	Quadrilateral quad;							\
688         quad.p = o, quad.q = o+i, quad.r = o+j, quad.s = o+i+j;			\
689         if (!_dxf_Quad(b, xf, 1, &quad, NULL, surface, clip, INV_VALID))	\
690 	    return ERROR;							\
691     }
692 
693 #define Q_INVALID(o, i, j, surface) {						\
694 	Quadrilateral quad;							\
695         quad.p = o, quad.q = o+i, quad.r = o+j, quad.s = o+i+j;			\
696         if (!_dxf_Quad(b, xf, 1, &quad, NULL, surface, clip, INV_INVALID))	\
697 	    return ERROR;							\
698     }
699 
700 #define QF_INVALID(o, i, j, col, surface) {					\
701 	Quadrilateral quad;							\
702         quad.p = o, quad.q = o+i, quad.r = o+j, quad.s = o+i+j;			\
703         if (!_dxf_QuadFlat(b, xf, 1, &quad, NULL, 				\
704 		       FCOLORP(col), BCOLORP(col), OPACITYP(col),		\
705 		       surface, clip, INV_INVALID))				\
706 	    return ERROR;							\
707     }
708 
709 #define QF(o, i, j, col, surface) {						\
710 	Quadrilateral quad;							\
711         quad.p = o, quad.q = o+i, quad.r = o+j, quad.s = o+i+j;			\
712         if (!_dxf_QuadFlat(b, xf, 1, &quad, NULL,				\
713 		       FCOLORP(col), BCOLORP(col), OPACITYP(col),		\
714 		       surface, clip, INV_VALID))				\
715 	    return ERROR;							\
716     }
717 
718 #define FORWARDS  1
719 #define BACKWARDS 2
720 
721 
722 static Error
VolumeRegularFace3(struct buffer * b,struct xfield * xf,int clip)723 VolumeRegularFace3(struct buffer *b, struct xfield *xf, int clip)
724 {
725     int i0, i1, i2, n0, n1, n2;
726     int col0, col1, col2;
727     InvalidComponentHandle ich = xf->iElts;
728 
729     EXPAND_POSITIONS;
730     orient(xf, 0);
731 
732     if (xf->colors_dep==dep_positions) {
733 
734 	if (ich)
735 	{
736 	    n1 = xf->n, col1 = xf->cn;
737 	    for (i1 = 0; i1 < CK1; i1++)
738 	    {
739 	        col0 = col1;
740 	        n0   = n1;
741 
742 	        for (i0 = 0; i0 < CK2; i0++)
743 	        {
744 		    Q(n0, S1, S2, 1);
745 
746 		    col0 += CS2;
747 		    n0   +=  S2;
748 		}
749 
750 		col1 += CS1;
751 		n1   +=  S1;
752 	    }
753 
754 	    n1 = xf->n, col1 = xf->cn;
755 	    for (i1 = 0; i1 < CK0; i1++)
756 	    {
757 	        col0 = col1;
758 	        n0   = n1;
759 
760 	        for (i0 = 0; i0 < CK2; i0++)
761 	        {
762 		    Q(n0, S0, S2, 1);
763 
764 		    col0 += CS2;
765 		    n0   +=  S2;
766 		}
767 
768 		col1 += CS0;
769 		n1   +=  S0;
770 	    }
771 
772 	    /*
773 	     * Now the front faces of all the cubes
774 	     */
775 	    n2 = xf->n, col2 = xf->cn;
776 	    for (i2 = 0; i2 <= CK2; i2++)
777 	    {
778 	        int surface = i2 == 0 || i2 == CK2;
779 
780 	        col1 = col2;
781 	        n1   = n2;
782 
783 	        for (i1 = 0; i1 < CK1; i1++)
784 	        {
785 		    col0 = col1;
786 		    n0   = n1;
787 
788 		    for (i0 = 0; i0 < CK0; i0++)
789 		    {
790 			if (i2 != 0 && DXIsElementInvalid(ich, col0-CS2))
791 			{
792 			    Q_INVALID(n0, S0, S1, surface);
793 			    Q_INVALID(n0-S2+S1, S0, S2, i1==(CK1-1));
794 			    Q_INVALID(n0+S0-S2, S2, S1, i0==(CK0-1));
795 			}
796 			else
797 			{
798 			    Q(n0, S0, S1, surface);
799 			    if (i2 != 0)
800 			    {
801 				Q(n0-S2+S1, S0, S2, i1==(CK1-1));
802 				Q(n0+S0-S2, S2, S1, i0==(CK0-1));
803 			    }
804 			}
805 
806 			col0 += CS0;
807 			n0   +=  S0;
808 		    }
809 
810 		    col1 += CS1;
811 		    n1   +=  S1;
812 		 }
813 
814 		 col2 += CS2, n2 += S2;
815 	    }
816 	}
817 	else
818 	{
819 	    n1 = xf->n, col1 = xf->cn;
820 	    for (i1 = 0; i1 < CK1; i1++)
821 	    {
822 	        col0 = col1;
823 	        n0   = n1;
824 
825 	        for (i0 = 0; i0 < CK2; i0++)
826 	        {
827 		    Q(n0, S1, S2, 1);
828 
829 		    col0 += CS2;
830 		    n0   +=  S2;
831 		}
832 
833 		col1 += CS1;
834 		n1   +=  S1;
835 	    }
836 
837 	    n1 = xf->n, col1 = xf->cn;
838 	    for (i1 = 0; i1 < CK0; i1++)
839 	    {
840 	        col0 = col1;
841 	        n0   = n1;
842 
843 	        for (i0 = 0; i0 < CK2; i0++)
844 	        {
845 		    Q(n0, S0, S2, 1);
846 
847 		    col0 += CS2;
848 		    n0   +=  S2;
849 		}
850 
851 		col1 += CS0;
852 		n1   +=  S0;
853 	    }
854 
855 	    /*
856 	     * Now the front faces of all the cubes
857 	     */
858 	    n2 = xf->n, col2 = xf->cn;
859 	    for (i2 = 0; i2 <= CK2; i2++)
860 	    {
861 	        int surface = i2 == 0 || i2 == CK2;
862 
863 	        col1 = col2;
864 	        n1   = n2;
865 
866 	        for (i1 = 0; i1 < CK1; i1++)
867 	        {
868 		    col0 = col1;
869 		    n0   = n1;
870 
871 		    for (i0 = 0; i0 < CK0; i0++)
872 		    {
873 			Q(n0, S0, S1, surface);
874 			if (i2 != 0)
875 			{
876 			    Q(n0-S2+S1, S0, S2, i1==(CK1-1));
877 			    Q(n0+S0-S2, S2, S1, i0==(CK0-1));
878 			}
879 
880 			col0 += CS0;
881 			n0   +=  S0;
882 		    }
883 
884 		    col1 += CS1;
885 		    n1   +=  S1;
886 		 }
887 
888 		 col2 += CS2, n2 += S2;
889 	    }
890 	}
891 
892     } else if (xf->colors_dep==dep_connections) {
893 
894 	if (ich)
895 	{
896 	    n1 = xf->n, col1 = xf->cn;
897 	    for (i1 = 0; i1 < CK1; i1++)
898 	    {
899 	        col0 = col1;
900 	        n0   = n1;
901 
902 	        for (i0 = 0; i0 < CK2; i0++)
903 	        {
904 		    QF(n0, S1, S2, col0, 1);
905 
906 		    col0 += CS2;
907 		    n0   +=  S2;
908 		}
909 
910 		col1 += CS1;
911 		n1   +=  S1;
912 	    }
913 
914 	    n1 = xf->n, col1 = xf->cn;
915 	    for (i1 = 0; i1 < CK0; i1++)
916 	    {
917 	        col0 = col1;
918 	        n0   = n1;
919 
920 	        for (i0 = 0; i0 < CK2; i0++)
921 	        {
922 		    QF(n0, S0, S2, col0, 1);
923 
924 		    col0 += CS2;
925 		    n0   +=  S2;
926 		}
927 
928 		col1 += CS0;
929 		n1   +=  S0;
930 	    }
931 
932 	    /*
933 	     * Now the front faces of all the cubes
934 	     */
935 	    n2 = xf->n, col2 = xf->cn;
936 	    for (i2 = 0; i2 <= CK2; i2++)
937 	    {
938 	        int surface = i2 == 0 || i2 == CK2;
939 
940 	        col1 = col2;
941 	        n1   = n2;
942 
943 	        for (i1 = 0; i1 < CK1; i1++)
944 	        {
945 		    col0 = col1;
946 		    n0   = n1;
947 
948 		    for (i0 = 0; i0 < CK0; i0++)
949 		    {
950 			if (i2 != 0 && DXIsElementInvalid(ich, col0-CS2))
951 			{
952 			    QF_INVALID(n0, S0, S1, col0-CS2, surface);
953 			    QF_INVALID(n0-S2+S1, S0, S2, col0-CS2, i1==(CK1-1));
954 			    QF_INVALID(n0+S0-S2, S2, S1, col0-CS2, i0==(CK0-1));
955 			}
956 			else if (i2 == 0)
957 			{
958 			    QF(n0, S0, S1, col0, surface)
959 			}
960 			else
961 			{
962 			    QF(n0, S0, S1, col0-CS2, surface);
963 			    QF(n0-S2+S1, S0, S2, col0-CS2, i1==(CK1-1));
964 			    QF(n0+S0-S2, S2, S1, col0-CS2, i0==(CK0-1));
965 			}
966 
967 			col0 += CS0;
968 			n0   +=  S0;
969 		    }
970 
971 		    col1 += CS1;
972 		    n1   +=  S1;
973 		 }
974 
975 		 col2 += CS2, n2 += S2;
976 	    }
977 
978 	}
979 	else
980 	{
981 	    n1 = xf->n, col1 = xf->cn;
982 	    for (i1 = 0; i1 < CK1; i1++)
983 	    {
984 	        col0 = col1;
985 	        n0   = n1;
986 
987 	        for (i0 = 0; i0 < CK2; i0++)
988 	        {
989 		    QF(n0, S1, S2, col0, 1);
990 
991 		    col0 += CS2;
992 		    n0   +=  S2;
993 		}
994 
995 		col1 += CS1;
996 		n1   +=  S1;
997 	    }
998 
999 	    n1 = xf->n, col1 = xf->cn;
1000 	    for (i1 = 0; i1 < CK0; i1++)
1001 	    {
1002 	        col0 = col1;
1003 	        n0   = n1;
1004 
1005 	        for (i0 = 0; i0 < CK2; i0++)
1006 	        {
1007 		    QF(n0, S0, S2, col0, 1);
1008 
1009 		    col0 += CS2;
1010 		    n0   +=  S2;
1011 		}
1012 
1013 		col1 += CS0;
1014 		n1   +=  S0;
1015 	    }
1016 
1017 	    /*
1018 	     * Now the front faces of all the cubes
1019 	     */
1020 	    n2 = xf->n, col2 = xf->cn;
1021 	    for (i2 = 0; i2 <= CK2; i2++)
1022 	    {
1023 	        int surface = i2 == 0 || i2 == CK2;
1024 
1025 	        col1 = col2;
1026 	        n1   = n2;
1027 
1028 	        for (i1 = 0; i1 < CK1; i1++)
1029 	        {
1030 		    col0 = col1;
1031 		    n0   = n1;
1032 
1033 		    for (i0 = 0; i0 < CK0; i0++)
1034 		    {
1035 			if (i2 == 0)
1036 			{
1037 			    QF(n0, S0, S1, col0, surface)
1038 			}
1039 			else
1040 			{
1041 			    QF(n0, S0, S1, col0-CS2, surface);
1042 			    QF(n0-S2+S1, S0, S2, col0-CS2, i1==(CK1-1));
1043 			    QF(n0+S0-S2, S2, S1, col0-CS2, i0==(CK0-1));
1044 			}
1045 
1046 			col0 += CS0;
1047 			n0   +=  S0;
1048 		    }
1049 
1050 		    col1 += CS1;
1051 		    n1   +=  S1;
1052 		 }
1053 
1054 		 col2 += CS2, n2 += S2;
1055 	    }
1056 
1057 	}
1058 
1059     }
1060 
1061     return OK;
1062 }
1063 
1064 
1065 /*
1066  * Find out which plane has smallest delta z spacing along view
1067  * axis between planes, and render only faces in that plane (plus
1068  * the edge faces that form the boundary).
1069  */
1070 
1071 static Error
VolumeRegularFace(struct buffer * b,struct xfield * xf,int clip)1072 VolumeRegularFace(struct buffer *b, struct xfield *xf, int clip)
1073 {
1074     int i0, i1, i2, n, n0, n1, n2, col0, col1, col2;
1075     InvalidComponentHandle ich = xf->iElts;
1076 
1077     EXPAND_POSITIONS;
1078     orient(xf, 0);
1079 
1080     if (xf->colors_dep==dep_positions) {
1081 
1082 	if (ich)
1083 	{
1084 	    n1 = xf->n, col1 = xf->cn;
1085 	    for (i1 = 0; i1 < CK1; i1++)
1086 	    {
1087 	        col0 = col1;
1088 	        n0   = n1;
1089 
1090 	        for (i0 = 0; i0 < CK2; i0++)
1091 	        {
1092 		    Q(n0, S1, S2, 1);
1093 
1094 		    col0 += CS2;
1095 		    n0   +=  S2;
1096 		}
1097 
1098 		col1 += CS1;
1099 		n1   +=  S1;
1100 	    }
1101 
1102 	    n1 = xf->n, col1 = xf->cn;
1103 	    for (i1 = 0; i1 < CK0; i1++)
1104 	    {
1105 	        col0 = col1;
1106 	        n0   = n1;
1107 
1108 	        for (i0 = 0; i0 < CK2; i0++)
1109 	        {
1110 		    Q(n0, S0, S2, 1);
1111 
1112 		    col0 += CS2;
1113 		    n0   +=  S2;
1114 		}
1115 
1116 		col1 += CS0;
1117 		n1   +=  S0;
1118 	    }
1119 
1120 	    /*
1121 	     * Now all the faces perpendicular to axis 2 in back-to-front
1122 	     * order
1123 	     */
1124 	    n2 = xf->n, col2 = xf->cn;
1125 	    for (i2 = 0; i2 <= CK2; i2++)
1126 	    {
1127 	        int surface = i2 == 0 || i2 == CK2;
1128 
1129 	        col1 = col2;
1130 	        n1   = n2;
1131 
1132 	        for (i1 = 0; i1 < CK1; i1++)
1133 	        {
1134 		    col0 = col1;
1135 		    n0   = n1;
1136 
1137 		    for (i0 = 0; i0 < CK0; i0++)
1138 		    {
1139 			if (i2 != 0 && DXIsElementInvalid(ich, col0-CS2))
1140 			{
1141 			    if (i0 != 0)
1142 				Q(n0-S2, S2, S1, 0);
1143 
1144 			    if (i1 != 0)
1145 				Q(n0-S2, S0, S2, 0);
1146 
1147 			    Q_INVALID(n0, S0, S1, surface);
1148 
1149 			    if (i1 != (CK1-1))
1150 				Q_INVALID(n0-S2+S1, S0, S2, 0);
1151 
1152 			    if (i0 != (CK0-1))
1153 				Q_INVALID(n0+S0-S2, S2, S1, 0);
1154 			}
1155 			else
1156 			    Q(n0, S0, S1, surface);
1157 
1158 			col0 += CS0;
1159 			n0   +=  S0;
1160 		    }
1161 
1162 		    col1 += CS1;
1163 		    n1   +=  S1;
1164 		 }
1165 
1166 		 col2 += CS2, n2 += S2;
1167 	    }
1168 
1169 	    /*
1170 	     * front face perpendicular to axis 0
1171 	     */
1172 	    n1 = xf->n+CK0*S0, col1 = xf->cn+(CK0-1)*CS0;
1173 	    for (i1 = 0; i1 < CK1; i1++)
1174 	    {
1175 	        col0 = col1;
1176 	        n0   = n1;
1177 
1178 	        for (i0 = 0; i0 < CK2; i0++)
1179 	        {
1180 		    if (DXIsElementInvalid(ich, col0))
1181 			Q_INVALID(n0, S1, S2, 1)
1182 		    else
1183 			Q(n0, S1, S2, 1);
1184 
1185 		    col0 += CS2;
1186 		    n0   +=  S2;
1187 		}
1188 
1189 		col1 += CS1;
1190 		n1   +=  S1;
1191 	    }
1192 
1193 	    /*
1194 	     * front face perpendicular to axis 1
1195 	     */
1196 	    n1 = xf->n+CK1*S1, col1 = xf->cn+(CK1-1)*CS1;
1197 	    for (i1 = 0; i1 < CK0; i1++)
1198 	    {
1199 	        col0 = col1;
1200 	        n0   = n1;
1201 
1202 	        for (i0 = 0; i0 < CK2; i0++)
1203 	        {
1204 		    if (DXIsElementInvalid(ich, col0))
1205 			Q_INVALID(n0, S0, S2, 1)
1206 		    else
1207 			Q(n0, S0, S2, 1);
1208 
1209 		    col0 += CS2;
1210 		    n0   +=  S2;
1211 		}
1212 
1213 		col1 += CS0;
1214 		n1   +=  S0;
1215 	    }
1216 	}
1217 	else
1218 	{
1219 	    FOR2 (K2) {
1220 		FOR1 (K1-1)
1221 		    FOR0 (K0-1)
1222 			Q(n, S0, S1, i2==0 || i2==K2-1)
1223 		if (i2 < K2-1) {
1224 		    FOR1 (K1) {
1225 			if (i1==0 || i1==K1-1)
1226 			    FOR0 (K0-1)
1227 				Q(n, S2, S0, 1)
1228 			if (i1 < K1-1) {
1229 			    Q(n1, S1, S2, 1)
1230 			    Q(n1+(K0-1)*S0, S1, S2, 1)
1231 			}
1232 		    }
1233 		}
1234 	    }
1235 	}
1236 
1237     } else if (xf->colors_dep==dep_connections) {
1238 
1239 	if (ich)
1240 	{
1241 	    n1 = xf->n, col1 = xf->cn;
1242 	    for (i1 = 0; i1 < CK1; i1++)
1243 	    {
1244 	        col0 = col1;
1245 	        n0   = n1;
1246 
1247 	        for (i0 = 0; i0 < CK2; i0++)
1248 	        {
1249 		    QF(n0, S1, S2, col0, 1);
1250 
1251 		    col0 += CS2;
1252 		    n0   +=  S2;
1253 		}
1254 
1255 		col1 += CS1;
1256 		n1   +=  S1;
1257 	    }
1258 
1259 	    n1 = xf->n, col1 = xf->cn;
1260 	    for (i1 = 0; i1 < CK0; i1++)
1261 	    {
1262 	        col0 = col1;
1263 	        n0   = n1;
1264 
1265 	        for (i0 = 0; i0 < CK2; i0++)
1266 	        {
1267 		    QF(n0, S0, S2, col0, 1);
1268 
1269 		    col0 += CS2;
1270 		    n0   +=  S2;
1271 		}
1272 
1273 		col1 += CS0;
1274 		n1   +=  S0;
1275 	    }
1276 
1277 	    /*
1278 	     * Now all the faces perpendicular to axis 2 in back-to-front
1279 	     * order
1280 	     */
1281 	    n2 = xf->n, col2 = xf->cn;
1282 	    for (i2 = 0; i2 <= CK2; i2++)
1283 	    {
1284 	        int surface = i2 == 0 || i2 == CK2;
1285 
1286 	        col1 = col2;
1287 	        n1   = n2;
1288 
1289 	        for (i1 = 0; i1 < CK1; i1++)
1290 	        {
1291 		    col0 = col1;
1292 		    n0   = n1;
1293 
1294 		    for (i0 = 0; i0 < CK0; i0++)
1295 		    {
1296 			if (i2 != 0 && DXIsElementInvalid(ich, col0-CS2))
1297 			{
1298 			    if (i0 != 0)
1299 				QF(n0-S2, S2, S1, col0-CS0, 0);
1300 
1301 			    if (i1 != 0)
1302 				QF(n0-S2, S0, S2, col0-CS1, 0);
1303 
1304 			    QF_INVALID(n0, S0, S1, col0-CS2, surface);
1305 
1306 			    if (i1 != (CK1-1))
1307 				QF_INVALID(n0-S2+S1, S0, S2, col0-CS2, 0);
1308 
1309 			    if (i0 != (CK0-1))
1310 				QF_INVALID(n0+S0-S2, S2, S1, col0-CS2, 0);
1311 			}
1312 			else if (i2 == 0)
1313 			    QF(n0, S0, S1, col0, surface)
1314 			else
1315 			    QF(n0, S0, S1, col0-CS2, surface);
1316 
1317 			col0 += CS0;
1318 			n0   +=  S0;
1319 		    }
1320 
1321 		    col1 += CS1;
1322 		    n1   +=  S1;
1323 		 }
1324 
1325 		 col2 += CS2, n2 += S2;
1326 	    }
1327 
1328 	    /*
1329 	     * front face perpendicular to axis 0
1330 	     */
1331 	    n1 = xf->n+CK0*S0, col1 = xf->cn+(CK0-1)*CS0;
1332 	    for (i1 = 0; i1 < CK1; i1++)
1333 	    {
1334 	        col0 = col1;
1335 	        n0   = n1;
1336 
1337 	        for (i0 = 0; i0 < CK2; i0++)
1338 	        {
1339 		    if (DXIsElementInvalid(ich, col0))
1340 			QF_INVALID(n0, S1, S2, col0, 1)
1341 		    else
1342 			QF(n0, S1, S2, col0, 1);
1343 
1344 		    col0 += CS2;
1345 		    n0   +=  S2;
1346 		}
1347 
1348 		col1 += CS1;
1349 		n1   +=  S1;
1350 	    }
1351 
1352 	    /*
1353 	     * front face perpendicular to axis 1
1354 	     */
1355 	    n1 = xf->n+CK1*S1, col1 = xf->cn+(CK1-1)*CS1;
1356 	    for (i1 = 0; i1 < CK0; i1++)
1357 	    {
1358 	        col0 = col1;
1359 	        n0   = n1;
1360 
1361 	        for (i0 = 0; i0 < CK2; i0++)
1362 	        {
1363 		    if (DXIsElementInvalid(ich, col0))
1364 			QF_INVALID(n0, S0, S2, col0, 1)
1365 		    else
1366 			QF(n0, S0, S2, col0, 1);
1367 
1368 		    col0 += CS2;
1369 		    n0   +=  S2;
1370 		}
1371 
1372 		col1 += CS0;
1373 		n1   +=  S0;
1374 	    }
1375 	}
1376 	else
1377 	{
1378 	    col2 = xf->cn;
1379 	    FOR2 (K2) {
1380 		col1 = col2;
1381 		FOR1 (K1-1) {
1382 		    col0 = col1;
1383 		    FOR0 (K0-1) {
1384 			QF(n, S0, S1, col0, i2==0 || i2==K2-1)
1385 			col0 += CS0;
1386 		    }
1387 		    col1 += CS1;
1388 		}
1389 		if (i2 < K2-1) {
1390 		    col1 = col2;
1391 		    FOR1 (K1) {
1392 			if (i1==0 || i1==K1-1) {
1393 			    col0 = col1;
1394 			    FOR0 (K0-1) {
1395 				QF(n, S2, S0, col0, 1)
1396 				col0 += CS0;
1397 			    }
1398 			}
1399 			if (i1 < K1-1) {
1400 			    QF(n1, S1, S2, col1, 1)
1401 			    QF(n1+(K0-1)*S0, S1, S2, col1+(CK0-1)*CS0, 1)
1402 			}
1403 			if (i1 < K1-2)
1404 			    col1 += CS1;
1405 		    }
1406 		}
1407 		if (i2 < K2-2)
1408 		    col2 += CS2;
1409 	    }
1410 	}
1411 
1412     }
1413 
1414     return OK;
1415 }
1416 
1417 
1418 static Error
VolumeRegularPlane(struct buffer * buf,struct xfield * xf,int clip)1419 VolumeRegularPlane(struct buffer *buf, struct xfield *xf, int clip)
1420 {
1421     int i2, n2, i;
1422     float dz, omul, cmul;
1423     Pointer op = xf->opacities;
1424 
1425 
1426 #if 0
1427     if (xf->colors_dep==dep_connections)
1428 	DXErrorReturn(ERROR_BAD_PARAMETER,
1429 		    "hi-res volumes must have position-dependent colors");
1430 #endif
1431 
1432     dz = orient(xf, 0);
1433     cmul = dz * xf->tile.color_multiplier;
1434     omul = dz * xf->tile.opacity_multiplier;
1435 
1436     /*
1437      * If connections-dep, copy conn-dep indices ->cn, ->ck, and ->cs
1438      * into ->n, ->k, and ->s and then do the same loop.
1439      * Couldn't do it earlier because we didn't know if we
1440      * would end up here or in say VolumeRegularFace which needs
1441      * both sets of numbers.  Also expand grid to fill the space.
1442      */
1443     if (xf->colors_dep==dep_connections) {
1444 	D0 = DXMul(D0, (double)K0/(double)CK0);
1445 	D1 = DXMul(D1, (double)K1/(double)CK1);
1446 	D2 = DXMul(D2, (double)K2/(double)CK2);
1447 	xf->n = xf->cn;
1448 	for (i=0; i<3; i++) {
1449 	    xf->k[i] = xf->ck[i];
1450 	    xf->s[i] = xf->cs[i];
1451 	}
1452     }
1453 
1454     FOR2 (K2) {
1455 	if (i2==0) continue;
1456 	_dxf_CompositePlane(buf,
1457 			buf->ox + xf->o.x + i2*D2.x,
1458 			buf->oy + xf->o.y + i2*D2.y,
1459 			xf->o.z + i2*D2.z,
1460 			clip,
1461 			D0.x, D0.y, D0.z, D1.x, D1.y, D1.z,
1462 			xf->cmap ?
1463 			    xf->fcst ?
1464 				(Pointer)((unsigned char *)(xf->fcolors) + 0) :
1465 				(Pointer)((unsigned char *)(xf->fcolors) + n2) :
1466 			    xf->fcst ?
1467 				(Pointer)((RGBColor *)(xf->fcolors) + 0) :
1468 				(Pointer)((RGBColor *)(xf->fcolors) + n2),
1469 			xf->fcst, xf->cmap,
1470 			xf->omap ?
1471 			    xf->ocst ?
1472 				(Pointer)(op? ((unsigned char *)op)+0 : NULL) :
1473 				(Pointer)(op? ((unsigned char *)op)+n2 : NULL) :
1474 			    xf->ocst ?
1475 				(Pointer)(op? ((float *)op)+0 : NULL) :
1476 				(Pointer)(op? ((float *)op)+n2 : NULL),
1477 			xf->ocst, xf->omap,
1478 			cmul, omul,
1479 			S0, S1, K0, K1);
1480     }
1481 
1482     return OK;
1483 }
1484