+
+ /* Fill the segment lists. */
+
+ if(sp->view.depth) { /* perspective view has special points */
+ if(b==0){ /* point of view */
+ sp->centre.x = X(0, b) * pp + sp->circle[0].x * pc;
+ sp->centre.y = Y(0, b) * pp + sp->circle[0].y * pc;
+ sp->centre.z = Z(0, b) * pp + sp->circle[0].z * pc;
+ /*printf("center: (%3.3f,%3.3f,%3.3f)\n",sp->centre.x, sp->centre.y, sp->centre.z);*/
+ }else if(b==1){ /* neighbour: used to compute local axes */
+ double x[3], p[3], x2=0, xp=0, C[3][3];
+ int j;
+
+ /* forward */
+ x[0] = X(0, 0) - X(1, 0);
+ x[1] = Y(0, 0) - Y(1, 0);
+ x[2] = Z(0, 0) - Z(1, 0);
+
+ /* neighbour */
+ p[0] = X(0, 1) - X(1, 0);
+ p[1] = Y(0, 1) - Y(1, 0);
+ p[2] = Z(0, 1) - Z(1, 0);
+
+ for(i=0; i<3; i++){
+ x2+= x[i]*x[i]; /* X . X */
+ xp+= x[i]*p[i]; /* X . P */
+ M[0][i] = x[i]; /* X */
+ }
+
+ for(i=0; i<3; i++) /* (X x P) x X */
+ M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
+
+ M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
+ M[2][1] = -x[0]*p[2] + x[2]*p[0];
+ M[2][2] = x[0]*p[1] - x[1]*p[0];
+
+ /* normalise axes */
+ for(j=0; j<3; j++){
+ double A=0;
+ for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
+ A=sqrt(A);
+ for(i=0; i<3; i++) M[j][i]/=A;
+ }
+
+ X(0, 1)=X(0, 0)+M[1][0]; /* adjust neighbour */
+ Y(0, 1)=Y(0, 0)+M[1][1];
+ Z(0, 1)=Z(0, 0)+M[1][2];
+
+ /* Look at trained bee into C matrix */
+ /* forward */
+ x[0] = 0 - sp->circle[0].x;
+ x[1] = 0 - sp->circle[0].y;
+ x[2] = 0 - sp->circle[0].z;
+
+ /* neighbour */
+ p[0] = sp->circle[0].x - sp->circle[1].x;
+ p[1] = sp->circle[0].y - sp->circle[1].y;
+ p[2] = sp->circle[0].z - sp->circle[1].z;
+
+ for(i=0; i<3; i++){
+ x2+= x[i]*x[i]; /* X . X */
+ xp+= x[i]*p[i]; /* X . P */
+ C[0][i] = x[i]; /* X */
+ }
+
+ for(i=0; i<3; i++) /* (X x P) x X */
+ C[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
+
+ C[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
+ C[2][1] = -x[0]*p[2] + x[2]*p[0];
+ C[2][2] = x[0]*p[1] - x[1]*p[0];
+
+ /* normalise axes */
+ for(j=0; j<3; j++){
+ double A=0;
+ for(i=0; i<3; i++) A+=C[j][i]*C[j][i]; /* sum squares */
+ A=sqrt(A);
+ for(i=0; i<3; i++) C[j][i]/=A;
+ }
+
+ /* Interpolate between Center and Trained Bee matrices */
+ /* This isn't very accurate and leads to weird transformations
+ * (shearing, etc), but it works. Besides, sometimes they look
+ * cool :) */
+ pp = pp * pp; /* Don't follow bee's direction until very close */
+ pc = 1 - pp;
+ for (i = 0; i < 3; i++)
+ for (j = 0; j < 3; j++)
+ M[i][j] = M[i][j] * pp + C[i][j] * pc;
+
+
+#if 0 /* display local axes for testing */
+ X(1, b)=X(0, 0);
+ Y(1, b)=Y(0, 0);
+ Z(1, b)=Z(0, 0);
+ }else if(b==2){
+ X(0, b)=X(0, 0)+0.5*M[0][0];
+ Y(0, b)=Y(0, 0)+0.5*M[0][1];
+ Z(0, b)=Z(0, 0)+0.5*M[0][2];
+ X(1, b)=X(0, 0);
+ Y(1, b)=Y(0, 0);
+ Z(1, b)=Z(0, 0);
+ }else if(b==3){
+ X(0, b)=X(0, 0)+1.5*M[2][0];
+ Y(0, b)=Y(0, 0)+1.5*M[2][1];
+ Z(0, b)=Z(0, 0)+1.5*M[2][2];
+ X(1, b)=X(0, 0);
+ Y(1, b)=Y(0, 0);
+ Z(1, b)=Z(0, 0);
+#endif
+ /* Draw a box... */
+ }
+ }
+ if (b >= MIN_BOX && b < MAX_BOX) {
+ int p1 = lines[b-MIN_BOX][0];
+ int p2 = lines[b-MIN_BOX][1];
+ X(0, b) = box[p1][0]; Y(0, b) = box[p1][1]; Z(0, b) = box[p1][2];
+ X(1, b) = box[p2][0]; Y(1, b) = box[p2][1]; Z(1, b) = box[p2][2];
+ }
+
+
+#if 0 /* Original code */
+ for(i=0; i<2; i++){
+ double x=X(i,b)-sp->centre.x;
+ double y=Y(i,b)-sp->centre.y;
+ double z=Z(i,b)-sp->centre.z;
+ double X=M[0][0]*x + M[0][1]*y + M[0][2]*z;
+ double Y=M[1][0]*x + M[1][1]*y + M[1][2]*z;
+ double Z=M[2][0]*x + M[2][1]*y + M[2][2]*z+sp->view.height;
+ double absx, absy;
+ if(sp->view.depth){
+ if(X <= 0) break;
+ absx=SCALE_X(sp->view.depth*Y/X);
+ absy=SCALE_Y(sp->view.depth*Z/X);
+ if(absx < -sp->width || absx > 2*sp->width ||
+ absy < -sp->height || absy > 2*sp->height)
+ break;
+ }else{
+ absx=SCALE_X(X);
+ absy=SCALE_Y(Y);
+ }
+ if(i){
+ sp->csegs[ix].x1 = (short) absx;
+ sp->csegs[ix].y1 = (short) absy;
+ }else{
+ sp->csegs[ix].x2 = (short) absx;
+ sp->csegs[ix].y2 = (short) absy;
+ }
+ }
+ if(i == 2) /* both assigned */
+ sp->cnsegs[col]++;
+#else
+ /* Chalky's code w/ clipping */
+ if (b < ((sp->mode & FLOW_BOX) ? 2 : MAX_BOX))
+ continue;
+ do {
+ double x1=X(0,b)-sp->centre.x;
+ double y1=Y(0,b)-sp->centre.y;
+ double z1=Z(0,b)-sp->centre.z;
+ double X1=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
+ double Y1=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
+ double Z1=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1+sp->view.height;
+ double absx1, absy1;
+ double x2=X(1,b)-sp->centre.x;
+ double y2=Y(1,b)-sp->centre.y;
+ double z2=Z(1,b)-sp->centre.z;
+ double X2=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
+ double Y2=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
+ double Z2=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2+sp->view.height;
+ double absx2, absy2;
+ if(sp->view.depth){
+ /* Need clipping if: is part of box, or close to viewer */
+ if ( (b >= MIN_BOX && b < MAX_BOX) || X1 <= 0.1 || X2 < 0.1)
+ if (clip(&X1, &Y1, &Z1, &X2, &Y2, &Z2))
+ break;
+ if (X1 <= 0 || X2 <= 0) break;
+ absx1=SCALE_X(sp->view.depth*Y1/X1);
+ absy1=SCALE_Y(sp->view.depth*Z1/X1);
+ if(absx1 < -sp->width || absx1 > 2*sp->width ||
+ absy1 < -sp->height || absy1 > 2*sp->height)
+ break;
+ absx2=SCALE_X(sp->view.depth*Y2/X2);
+ absy2=SCALE_Y(sp->view.depth*Z2/X2);
+ if(absx2 < -sp->width || absx2 > 2*sp->width ||
+ absy2 < -sp->height || absy2 > 2*sp->height)
+ break;
+ }else{
+ absx1=SCALE_X(X1);
+ absy1=SCALE_Y(Y1);
+ absx2=SCALE_X(X2);
+ absy2=SCALE_Y(Y2);
+ }
+
+ sp->csegs[ix].x1 = (short) absx1;
+ sp->csegs[ix].y1 = (short) absy1;
+ sp->csegs[ix].x2 = (short) absx2;
+ sp->csegs[ix].y2 = (short) absy2;
+
+ sp->cnsegs[col]++;
+ } while (0);
+#endif