+/* Clipping planes */
+#define PLANES 5
+static double plane_orig[][2][3] = {
+ /* X goes into screen, Y goes right, Z goes down(up?) */
+ /* {Normal}, {Point} */
+ { {1.0, 0, 0}, {0.01, 0, 0} },
+ { {1.0, 1.0, 0.0}, {0, 0, 0} },
+ { {1.0,-1.0, 0.0}, {0, 0, 0} },
+ { {1.0, 0.0, 1.0}, {0, 0, 0} },
+ { {1.0, 0.0,-1.0}, {0, 0, 0} }
+};
+static double plane[PLANES][2][3];
+static double plane_d[PLANES];
+
+#define BOX_P 32
+#define BOX_L 36
+#define MIN_BOX (3)
+#define MAX_BOX (MIN_BOX + BOX_L)
+/* Points that make up the box (normalized coordinates) */
+static double box_orig[][3] = {
+ {1,1,1}, /* 0 */
+ {1,1,-1}, /* 1 */
+ {1,-1,-1}, /* 2 */
+ {1,-1,1}, /* 3 */
+ {-1,1,1}, /* 4 */
+ {-1,1,-1}, /* 5 */
+ {-1,-1,-1},/* 6 */
+ {-1,-1,1}, /* 7 */
+ {1, .8, .8},
+ {1, .8,-.8},
+ {1,-.8,-.8},
+ {1,-.8, .8},
+ { .8,1, .8},
+ { .8,1,-.8},
+ {-.8,1,-.8},
+ {-.8,1, .8},
+ { .8, .8,1},
+ { .8,-.8,1},
+ {-.8,-.8,1},
+ {-.8, .8,1},
+ {-1, .8, .8},
+ {-1, .8,-.8},
+ {-1,-.8,-.8},
+ {-1,-.8, .8},
+ { .8,-1, .8},
+ { .8,-1,-.8},
+ {-.8,-1,-.8},
+ {-.8,-1, .8},
+ { .8, .8,-1},
+ { .8,-.8,-1},
+ {-.8,-.8,-1},
+ {-.8, .8,-1}
+};
+
+/* Container for scaled box points */
+static double box[BOX_P][3];
+
+/* Lines connecting the box dots */
+static double lines[][2] = {
+ {0,1}, {1,2}, {2,3}, {3,0}, /* box */
+ {4,5}, {5,6}, {6,7}, {7,4},
+ {0,4}, {1,5}, {2,6}, {3,7},
+ {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
+ {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
+ {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
+ {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
+ {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
+ {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
+};
+
+/* Boundaries of bees */
+double xmin, xmax;
+double ymin, ymax;
+double zmin, zmax;
+
+void init_clip(flowstruct *sp)
+{
+ int i;
+
+ /* Scale the planes to the screen. I had to invert the projection
+ * algorithms so that when projected they would be right at the edge of the
+ * screen. */
+
+ /* #### jwz: I'm not really sure what it means when sp->view.depth is 0
+ in here -- what's the right thing to do? */
+
+ double width = (sp->view.depth
+ ? sp->size/sp->view.depth/2
+ : 1);
+ double height = (sp->view.depth
+ ? (sp->size/sp->view.depth/2*
+ sp->view.height/sp->view.height)
+ : 1);
+ for (i = 0; i < PLANES; i++) {
+ /* Copy orig planes into planes, expanding <-> clippings */
+ plane[i][0][0] = plane_orig[i][0][0];
+ plane[i][0][1] = plane_orig[i][0][1] / width;
+ plane[i][0][2] = plane_orig[i][0][2] / height;
+ plane[i][1][0] = plane_orig[i][1][0];
+ plane[i][1][1] = plane_orig[i][1][1];
+ plane[i][1][2] = plane_orig[i][1][2];
+
+ /* Calculate the 'd' part of 'ax + by + cz = d' */
+ plane_d[i] = - plane[i][0][0] * plane[i][1][0];
+ plane_d[i] -= plane[i][0][1] * plane[i][1][1];
+ plane_d[i] -= plane[i][0][2] * plane[i][1][2];
+ }
+ xmin = X(0, i); xmax = X(0,i);
+ ymin = Y(0, i); ymax = Y(0,i);
+ zmin = Z(0, i); zmax = Z(0,i);
+}
+
+/* Scale the box defined above to fit around all points */
+void create_box(flowstruct *sp)
+{
+ int i = MAX_BOX;
+ double xmid, ymid, zmid;
+ double xsize, ysize, zsize;
+ double size;
+
+ /* Count every 5th point for speed.. */
+ for (; i < sp->beecount; i += 5) {
+ if ( X(0,i) < xmin ) xmin = X(0, i);
+ else if ( X(0,i) > xmax ) xmax = X(0, i);
+ if ( Y(0,i) < ymin ) ymin = Y(0, i);
+ else if ( Y(0,i) > ymax ) ymax = Y(0, i);
+ if ( Z(0,i) < zmin ) zmin = Z(0, i);
+ else if ( Z(0,i) > zmax ) zmax = Z(0, i);
+ }
+ xmid = (xmax+xmin)/2;
+ ymid = (ymax+ymin)/2;
+ zmid = (zmax+zmin)/2;
+ xsize = xmax - xmin;
+ ysize = ymax - ymin;
+ zsize = zmax - zmin;
+ size = xsize;
+ if (ysize> size) size = ysize;
+ if (zsize > size) size = zsize;
+ size /= 2;
+
+ /* Scale box */
+ for (i = 0; i < BOX_P; i++) {
+ box[i][0] = box_orig[i][0] * size + xmid;
+ box[i][1] = box_orig[i][1] * size + ymid;
+ box[i][2] = box_orig[i][2] * size + zmid;
+ }
+
+}
+
+/* Returns true if point is infront of the plane (rather than behind) */
+int infront_of(double x, double y, double z, int i)
+{
+ double sum = plane[i][0][0]*x + plane[i][0][1]*y + plane[i][0][2]*z + plane_d[i];
+ return sum >= 0.0;
+}
+
+/* Returns true if line was behind a clip plane, or clips the line */
+int clip(double *x1, double *y1, double *z1, double *x2, double *y2, double *z2)
+{
+ int i;
+ for (i = 0; i < PLANES; i++) {
+ double t;
+ double x, y, z; /* Intersection point */
+ double dx, dy, dz; /* line delta */
+ int front1, front2;
+ front1 = infront_of(*x1, *y1, *z1, i);
+ front2 = infront_of(*x2, *y2, *z2, i);
+ if (!front1 && !front2) return 1;
+ if (front1 && front2) continue;
+
+ dx = *x2 - *x1;
+ dy = *y2 - *y1;
+ dz = *z2 - *z1;
+
+ /* Find t in line equation */
+ t = ( plane_d[i] -
+ plane[i][0][0]*(*x1) - plane[i][0][1]*(*y1) - plane[i][0][2]*(*z1) )
+ /
+ ( plane[i][0][0]*dx + plane[i][0][1]*dy + plane[i][0][2]*dz );
+
+ x = *x1 + dx * t;
+ y = *y1 + dy * t;
+ z = *z1 + dz * t;
+ /* Make point that was behind to be the intersect */
+ if (front2) {
+ *x1 = x;
+ *y1 = y;
+ *z1 = z;
+ } else {
+ *x2 = x;
+ *y2 = y;
+ *z2 = z;
+ }
+ }
+ return 0;
+}
+
+