home *** CD-ROM | disk | FTP | other *** search
/ Collection of Hack-Phreak Scene Programs / cleanhpvac.zip / cleanhpvac / FAQSYS18.ZIP / FAQS.DAT / QUATERN.TXT < prev    next >
Text File  |  1996-08-18  |  16KB  |  437 lines
  1.  
  2. A Quaternion is a 'Julia-Set' in 4D. It has the same formula, Z(n) =
  3. Z(n-1)^2 + C, but Z and C are quaternion numbers, not complex. A quaternion
  4. number has 1 real value, and 3 imaginary, often showed as C = (a, bi, cj,
  5. dk)
  6. NOTE: Quaternion calculus is not commutativ ! (That is, bi*cj is not equal
  7. to cj*bi)
  8.  
  9. As I showed in my Mandelbrot-tutorial,
  10.  
  11. (a, bi)^2.real      = a*a - b*b
  12. (a, bi)^2.imaginary = 2*a*b
  13.  
  14. With Quaternions, you get the same thing:
  15.  
  16. (a, bi, cj, dk)^2.real         = a*a - b*b - c*c - d*d
  17. (a, bi, cj, dk)^2.imaginary(i) = 2*a*b
  18. (a, bi, cj, dk)^2.imaginary(j) = 2*a*c
  19. (a, bi, cj, dk)^2.imaginary(k) = 2*a*d
  20.  
  21. The Quaternion-Set is found the same way as Mandelbrot- and Julia-Sets, by
  22. iterating the function Z(n) = Z(n-1)^2 + C , and see if Z(n)->infinity. The
  23. formula would be something like this : (In 'C')
  24.  
  25. double calc_l(double x, double y, double z, double w)
  26.   double length;
  27.   double temp;
  28.   int m=0;
  29.   do {
  30.     temp=x+x;
  31.     x=x*x-y*y-z*z-w*w+cr;
  32.     y=temp*y+ci;
  33.     z=temp*z+cj;
  34.     w=temp*w+ck;
  35.  
  36.     m++;
  37.     length=x*x+y*y+z*z+w*w;
  38.   } while (!((m>iter) && (length>4)));
  39.  return length;
  40. }
  41.  
  42. This function is called with the four quaternion values, and uses the
  43. global constants cr, ci, cj and ck in the calculation. It returns the
  44. length of the vector represented by (x,y,z,w), the check to see if this is
  45. 'infinite' is done where we call the function. (Of course you can make this
  46. function as you want. This is just an example..)
  47.  
  48. Everyone used to making 2D-fractals (like Mandelbrot and Julia) might see
  49. the problem with Quaternions: We are using 4-dimensional values. What do we
  50. do if a 4D point is inside the quaternion-set ? (With 'standard' 2D
  51. fractals we only gave the correct 2D-pixel on screen a specific color). It
  52. is the 4D-thing that makes Quaternions so exciting, but also difficult to
  53. make. We are actually making a film of a 3D object! (But, usually we keep
  54. the 4th dimension constant, so we 'only' make a still-frame of a 3D object,
  55. like the one below)
  56.  
  57. [Image]QUATERN1.GIF
  58.  
  59. What I usually do when calculating Quaternions, is scanning a 3D room
  60. (having the 4th dimension constant), building up a Z-Buffer, and then
  61. 'ray-tracing' it.
  62.  
  63. The 3D-room has positive X-values to the right, positive Y-values at the
  64. bottom, and positive Z-values inside the screen. (See illustration)
  65.  
  66. [Image]
  67.  
  68. Since this room might be rotated, I often calculate the unit-vectors ex, ey
  69. and ez. (They could be the vectors on the illustration, but usually they
  70. are a lot smaller. I guess you want images with better resolution than 2*2
  71. or something :-)
  72. For the X-vector, this would be calculating the upper left corner and the
  73. upper right corner, and divide the delta x, y and z by the horisontal
  74. screen-size.
  75.  
  76. dx.x=(rightx-leftx)/screenx;dx.y=(righty-lefty)/screenx;dx.z=(rightz-leftz)/screenx;
  77. The other two would be:
  78. dy.x=(rightx-leftx)/screeny;dy.y=(righty-lefty)/screeny;dy.z=(rightz-leftz)/screeny;
  79. dz.x=(rightx-leftx)/screenz;dz.y=(righty-lefty)/screenz;dz.z=(rightz-leftz)/screenz;
  80.  
  81. Now, a given point (x,y,z) would in our rotated room be
  82. ((x*dx.x+y*dy.x+z*dz.x),(x*dx.y+y*dy.y+z*dz.y),(x*dx.z+y*dy.z+z*dz.z)).
  83. This enables us to always trace trough a standard space (x,y,z), even if
  84. the images is to be rotated!
  85.  
  86. For Mandelbrot-Sets, we had:
  87.  
  88. FOR EVERY y:
  89.   FOR EVERY x:
  90.     COMPUTE Z(n)=Z(n-1)^2 + C for (x,y);
  91.  
  92. With Quaternions, we get:
  93.  
  94. FOR EVERY y:
  95.   FOR EVERY x:
  96.     FOR EVERY z UNTIL INSIDE SET:
  97.       COMPUTE Z(n)=Z(n-1)^2 + C for ROTATED(x,y,z);
  98.  
  99. [Image]Difference between 2D and 3D
  100.  
  101. When we, for a given (x,y), have found a z that result in (x,y,z) being
  102. inside the set, we put that z-value into our Z-Buffer. (If we reach zmax
  103. without getting into the QuaternionSet, we put zmax into the Z-Buffer).
  104. When we have enough values in the Z-Buffer, we can 'ray-trace' it, and we
  105. get our final, increadible-looking, image.
  106. To 'ray-trace' a given point, we need to know 4 points surrounding it.
  107. Knowing 4 points, we can calculate 2 vectors in 3D space, and calculate the
  108. vector-normal to those two vectors.
  109.  
  110. [Image]
  111.  
  112. In the left part of the illustration above, we are 'ray-tracing' the point
  113. represented by the thick red line (just under the intersection of all the
  114. black lines). If we imagine this point being (0,0), we also know point
  115. (-1,0), (1,0), (0,-1) and (0,1). You can see that I have connected (-1,0)
  116. and (1,0), and (0,-1) and (0,1). This black lines are the two I was talking
  117. about above. You can also see a black arrow where the two lines intersect.
  118. That is supposed to be the vector that is normal to both of the lines. (Not
  119. easy to explain, this one..).
  120. OK. That black arrow should also be normal to the Quaternion set in point
  121. (0,0). (Since fractals don't have a continous surface, it don't have a
  122. normal, but our arrow is a good aproximation)
  123.  
  124. The reason why I have been talking about 'ray-tracing', is that we have a
  125. light-source in a given point. When we know the normal of a given point, we
  126. can also calculate the angle between that normal and the line from our
  127. point to the light-source. (As I try to show you in the right part of the
  128. illustration above..)
  129. The angle will be a number between -pi and pi, but we only need the
  130. absolute value of it, so we have a number between 0 and pi. If the angle
  131. between them is 0, they are parallell, and the point is facing away from
  132. the lightsource. If the angle is pi/2, the vector is normal to the ray from
  133. the lightsource, and if the angle is pi, it is pointing directly at the
  134. lightsource. Values from 0 - pi/2 should be colored black, and values from
  135. pi/2 - pi would result in brighter and brighter color. That is the whole
  136. magic behing Quaternions. Now it's just to code it :-)
  137.  
  138. [Image]QUATERN2.GIF
  139.  
  140. [Image]QUATERN3.GIF
  141.  
  142. Here is an explanation of some of the things I do in my code
  143.  
  144.   int zant=250;
  145.   int zant1=25;
  146.   int pixsize=2;
  147.   int vissize=1;
  148.  
  149. zant is the number of steps I trace in the Z-direction. Bigger zant ->
  150. better resolution, but also longer computing-time.
  151. zant1 is the number of steps I divide one z-step into. The reason why I do
  152. this, is to get a good resolution without having a -very- big zant. I trace
  153. into to Z-axis with big steps until I hit the Quaternion, and then I trace
  154. out again till I'm out of the set, this time using small steps.
  155. pixsize is used in preview-calculations. If pixsize is 4, I only compute
  156. every 4th pixel, and will quick get a rough idea about how the image will
  157. be.
  158. vissize tells me how big each pixel is. If I use pixsize=4 and vissize=1, I
  159. will get a detailed image covering only 1/16th of the screen. If I use
  160. pixsize=4 and vissize=4, I will get a low-detail image covering the entire
  161. screen.
  162.  
  163.   double xmin=-1.66, xmax=1;   //this values define
  164.   double ymin=-1, ymax=1;       //the 3D space
  165.   double zmin=-1.5, zmax=1.5; //to search for the QuaternionSet
  166.  
  167.   int iter=30;     //Number of iterations. With Quaternions this can be a small number
  168.  
  169. The lines above tell the program where to look for the Quaternion, and how
  170. many iterations to use. One of the cool things about Quaternions, is that
  171. you can get nice images even with few iterations. 2 of the images I show
  172. you on this page is computed with 10 iterations !!
  173.  
  174.   double lightx=-1, lighty=1, lightz=-3;  //Location of lightsource
  175.  
  176.   double vx=0, vy=0, vz=0;     //Rotationangle (in degrees) around x-, y- and z-axis
  177.  
  178. Here we define where the lightsource is positioned, and how the image is to
  179. be rotated.
  180.  
  181.   double cr=-0.46;  //constant real value
  182.   double ci=0.20;  //constant imaginary(1) value
  183.   double cj=0.6;  //constant imaginary(2) value
  184.   double ck=0.25;   //constant imaginary(3) value
  185.  
  186.   double wk=0.022;   //4th dimension
  187.  
  188. To compute a Julia-Set, you use 2 constant values. With Quaternions we have
  189. to use 4. Since the Quaternion is 4D, we also have to keep one of the
  190. dimensions constant. (I have chosen the 4th)
  191.  
  192.   int background = 0;  //0 -> no background.
  193.  
  194. This line simply tells the program if it shall raytrace the background, or
  195. just ignore it. I prefere no background..
  196.  
  197. [.. Some lines of code deleted]
  198.  
  199. void rotate3D(double &x,double &y,double &z)
  200. {
  201.   x-=origx;y-=origy;z-=origz;
  202.   double xy=y*cosx-z*sinx;
  203.   double xz=y*sinx+z*cosx;
  204.   double xx=x;
  205.   x=xx;
  206.   y=xy;
  207.   z=xz;
  208.   double yx=x*cosy+z*siny;
  209.   double yz=-x*siny+z*cosy;
  210.   x=yx;
  211.   z=yz;
  212.   double zx=x*cosz-y*sinz;
  213.   double zy=x*sinz+y*cosz;
  214.   x=zx;
  215.   y=zy;
  216.   x+=origx;y+=origy;z+=origz;
  217. }
  218.  
  219. This routine simply rotates a point in 3D. If you don't understand it,
  220. never mind..
  221.  
  222. void rotatevalues()
  223. {
  224.   rminx=xmin;rminy=ymin;rminz=zmin;
  225.   rotate3D(rminx, rminy, rminz);
  226.   tempx=xmax;tempy=ymin;tempz=zmin;
  227.   rotate3D(tempx, tempy, tempz);
  228.   dxx=(tempx-rminx)/sx;dxy=(tempy-rminy)/sx;dxz=(tempz-rminz)/sx;
  229.   tempx=xmin;tempy=ymax;tempz=zmin;
  230.   rotate3D(tempx, tempy, tempz);
  231.   dyx=(tempx-rminx)/sy;dyy=(tempy-rminy)/sy;dyz=(tempz-rminz)/sy;
  232.   tempx=xmin;tempy=ymin;tempz=zmax;
  233.   rotate3D(tempx, tempy, tempz);
  234.   dzx=(tempx-rminx)/zant;dzy=(tempy-rminy)/zant;dzz=(tempz-rminz)/zant;
  235.   dzx1=dzx/zant1;dzy1=dzy/zant1;dzz1=dzz/zant1;
  236. }
  237.  
  238. This routine is called during setup, and rotates the 3D-room and calculates
  239. the ex-, ey- and ez-vectors. This is only done for speed purposes, you
  240. could rotate every point when you calculate. (I like it 'Quick'n'dirty')
  241.  
  242. double calc_l(double x, double y, double z)
  243.   double lengde;
  244.   double temp;
  245.   double w=wk;
  246.   int m=0;
  247.   do {
  248.         temp=x+x;
  249.         x=x*x-y*y-z*z-w*w+cr;
  250.         y=temp*y+ci;
  251.         z=temp*z+cj;
  252.         w=temp*w+ck;
  253.  
  254.         m++;
  255.         lengde=x*x+y*y+z*z+w*w;
  256.   } while (!((m>iter) && (lengde>2)));
  257.  return lengde;
  258. }
  259.  
  260. Now we are getting somewhere!! This routine calculates a given 3D-point,
  261. and returns the length of the Quaternion vector. I had to invert the
  262. while-condition, because Netscape thought it was a hypertext-command..
  263. Please note that the number '2' in the while-condition is the
  264. bailout-value. This can be changed.
  265.  
  266. double calc_angle(double w,double e,double n,double s,double cx,double cy,double cz)
  267. {
  268.   double lightdx=cx-lightx;
  269.   double lightdy=cy-lighty;
  270.   double lightdz=cz-lightz;
  271.  
  272.   double lightlength=sqrt(lightdx*lightdx+lightdy*lightdy+lightdz*lightdz);
  273.  
  274.   double fx=/*(0)*(s-n)*/-(e-w)*(dy+dy);
  275.   double fy=/*(e-w)*(0)*/-(dx+dx)*(s-n);
  276.   double fz=(dx+dx)*(dy+dy)/*-(0)*(0)*/;
  277.  
  278.   double flength=sqrt(fx*fx+fy*fy+fz*fz);
  279.   double c=(fx*lightdx+fy*lightdy+fz*lightdz)/(flength*lightlength);
  280.   return c;
  281. }
  282.  
  283. Here I calculate the light-to-point-vector, the length of it, the direction
  284. of the vector normal to the Quaternion-Set in a given point, its length and
  285. finally the angle between the lightbeam and the normal-vector. This should
  286. be standard calculus. If you don't understand it, read a book on
  287. vector-calculus.(and give it to me when you are finished)
  288.  
  289. void show_buffer(int y)
  290. {
  291.   double a;
  292.  
  293.   for (int t=1; tzmax) && (background==0)) {
  294.                          setfillstyle(1,0);
  295.                   } else if (a<0) {
  296.                          setfillstyle(1,1);
  297.                   } else {
  298.                          setfillstyle(1,1+15*a);
  299.                   }
  300.                   bar(t*vissize,(y+i)*vissize,t*vissize+vissize-1,(y+i)*vissize+vissize-1);
  301.                 }
  302.          }
  303.   }
  304.  
  305.   for (t=0; t<640; t++) {
  306.          z_buffer[t][0]=z_buffer[t][8];
  307.          z_buffer[t][1]=z_buffer[t][9];
  308.   }
  309.   buffer_y=2;
  310. }
  311.  
  312. Now THIS is a routine I hate ! Here I 'ray-trace' the Z-Buffer. Since our
  313. 'beloved' IBM-compatibles won't accept arrays bigger than 64Kb, this
  314. routine is quite messy. Sorry.
  315. (What I do, is for a every 10th line, trace line 1-8, copy line 8 to line
  316. 0, copy line 9 to line 1, and continue computing the fractal)
  317.  
  318. void main()
  319. {
  320.   int pz, pz1;
  321.   double l;
  322.  
  323.         int gdriver = VGA, gmode = VGAHI, errorcode;
  324.         errorcode = registerbgidriver(EGAVGA_driver);
  325.         if (errorcode < 0) {
  326.                 printf("Graphics error: %s\n", grapherrormsg(errorcode));
  327.                 exit(1);
  328.         }
  329.  
  330.         initgraph(&gdriver, &gmode, "");
  331.  
  332.         errorcode = graphresult();
  333.         if (errorcode != grOk) {
  334.                 printf("Graphics error: %s\n", grapherrormsg(errorcode));
  335.                 exit(1);
  336.         }
  337.  
  338.   for (int i=1; i<16; i++) {
  339.          setrgbpalette(i, 0, i*4, 0);
  340.          setpalette(i, i);
  341.   }
  342.   setrgbpalette(0,0,0,63);
  343.   setpalette(0, 0);
  344.  
  345. Here I open a 640*480, 16-color screen (remember to link the
  346. egavga.obj-file!!). What, you don't like 640*480, 16-color ?!? Then change
  347. it !!
  348. I also initialize those 'pretty' green-colors on blue background..
  349.  
  350.   sx=getmaxx()/pixsize;sy=getmaxy()/pixsize;
  351.   dx=(xmax-xmin)/sx;
  352.   dy=(ymax-ymin)/sy;
  353.   dz=(zmax-zmin)/zant;
  354.  
  355.   origx=(xmin+xmax)/2;
  356.   origy=(ymin+ymax)/2;
  357.   origz=(zmin+zmax)/2;
  358.  
  359.   vx=vx/180*3.14159265;
  360.   vy=vy/180*3.14159265;
  361.   vz=vz/180*3.14159265;
  362.  
  363.   cosx=cos(vx);cosy=cos(vy);cosz=cos(vz);
  364.   sinx=sin(vx);siny=sin(vy);sinz=sin(vz);
  365.  
  366.   rotatevalues();
  367.  
  368. Calculates the screen-width, and the steps in our non-rotated 3D-space. I
  369. also calculate the origin, so I can rotate around other points than
  370. (0,0,0).
  371. I also convert the angles from degrees to radians, calculating the sinus
  372. and cosinus (precalced for speed purposes) and rotates the now-so-famous
  373. 3D-room
  374.  
  375.   buffer_y=0;
  376.   for (int py=0; py<=sy; py++) {
  377.          for (int px=0; px<=sx; px++) {
  378.                 pz=0;
  379.                 tempx=rminx+px*dxx+py*dyx/*+pz*dzx*/;
  380.                 tempy=rminy+px*dxy+py*dyy/*+pz*dzy*/;
  381.                 tempz=rminz+px*dxz+py*dyz/*+pz*dzz*/;
  382.                 do {
  383.                   tempx+=dzx;
  384.                   tempy+=dzy;
  385.                   tempz+=dzz;
  386.                   l=calc_l(tempx,tempy,tempz);
  387.                   pz++;
  388.                 } while ((l>2) && (!(pz>=zant)));
  389.                 pz1=0;
  390.                 do {
  391.                          pz1++;
  392.                          tempx-=dzx1;
  393.                          tempy-=dzy1;
  394.                          tempz-=dzz1;
  395.                          l=calc_l(tempx,tempy,tempz);
  396.                   } while (!(l>2));
  397.                 if (pz < zant)
  398.                   z_buffer[px][buffer_y]=zmin+pz*dz-pz1*dz/zant;
  399.                 else
  400.                   z_buffer[px][buffer_y]=zmax+dz;
  401.                 setfillstyle(1,15-pz/10);
  402.                 bar(px*vissize,py*vissize,px*vissize+vissize-1,py*vissize+vissize-1);
  403.                 if (kbhit()) break;
  404.          }
  405.          buffer_y++;
  406.          if (buffer_y==10) show_buffer(py-9);
  407.          if (kbhit()) break;
  408.   }
  409.   if (!kbhit()) {
  410.          show_buffer(py-buffer_y);
  411.          cout  '\7';
  412.   }
  413.   getch();
  414.   closegraph();
  415. }
  416.  
  417. The main-routine. Here I trace throug every Y, every X, and the necesarry
  418. number of Z-values, put the result into the Z-Buffer, traces this when it
  419. is full, exits if you press a key, beeps, and close the graphics-screen
  420. when finished. Phew !
  421.  
  422. If you are still reading, that means you must be -very- interested in
  423. Quaternions. Perhaps so interested that you want the C++ Source-code ?
  424. (Never mind if you don't have a math-prossessor. SX's are too slow..)
  425.  
  426. If there are any (if ??) confusing things here, please let me know. I will
  427. try to make it more understandable, but I don't have much time. (I join the
  428. army the 11th July..)
  429.  
  430. [Image]QUATERN4.GIF
  431.  
  432. [Image]QUATERN5.GIF
  433.  
  434.                                   [Image]
  435.                      visitors sice 18th November 1995.
  436. Back to homepage
  437.